Computational Chem 6

Lecture Notes
in
Computational Chemistry
Electronic Structure Theory
Jürg Hutter
Physical Chemistry Institute
University of Zurich
Winterthurerstrasse 190
8057 Zurich, Switzerland
hutter@pci.unizh.ch
August 11, 2005

ii
Contents
I Lectures 1
1 Basic Quantum Mechanics 3
1.1 Operators in quantum mechanics . . . . . . . . . . . . . . . . 3
1.1.1 Eigenfunctions and eigenvalues . . . . . . . . . . . . . 3
1.1.2 Commutators . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.3 Construction of operators . . . . . . . . . . . . . . . . 5
1.2 Postulates of quantum mechanics . . . . . . . . . . . . . . . . 5
1.2.1 Postulates . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Separation of the Schrödinger equation . . . . . . . . . 9
1.3 Hydrogen Atom . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Angular momentum operators . . . . . . . . . . . . . . 10
1.3.2 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.3 Radial function of hydrogen atom . . . . . . . . . . . . 17
2 Basic Mathematical Review 21

2.1 Linear algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.1 Vector spaces . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.2 Operators . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.1.3 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.1.4 Determinants . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.5 Dirac notation . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.6 Change of basis . . . . . . . . . . . . . . . . . . . . . . 27
2.1.7 Eigenvalue problem . . . . . . . . . . . . . . . . . . . . 29
2.1.8 Functions of matrices . . . . . . . . . . . . . . . . . . . 32
2.2 Perturbation theory . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.2 Energy corrections . . . . . . . . . . . . . . . . . . . . 35
(1) (2)
2.2.3 Corrections Ψn and En . . . . . . . . . . . . . . . . 36
2.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3 The variation method . . . . . . . . . . . . . . . . . . . . . . . 37
2.3.1 The variation principle . . . . . . . . . . . . . . . . . . 38
iii
2.3.2 Linear variations . . . . . . . . . . . . . . . . . . . . . 38
3 Molecular Hamiltonian 41
3.1 Time-independent non-relativistic Schrödinger equation . . . . 41
3.2 Born–Oppenheimer Approximation . . . . . . . . . . . . . . . 42
3.2.1 Part I . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2.2 Part II . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2.3 Part III . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2.4 Order of magnitude of effects for H2 . . . . . . . . . . 45
3.3 Stationary points . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4 Internal coordinates . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4.1 Z-matrix . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4 Two-Electron Systems and Spin 53

4.1 Helium atom: Ground state . . . . . . . . . . . . . . . . . . . 53
4.1.1 Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1.2 Orbital model . . . . . . . . . . . . . . . . . . . . . . . 54
4.1.3 Formal separation . . . . . . . . . . . . . . . . . . . . . 54
4.1.4 Perturbation theory . . . . . . . . . . . . . . . . . . . . 55
4.1.5 Variational calculations . . . . . . . . . . . . . . . . . . 55
4.1.6 Overview of ground state energies . . . . . . . . . . . . 56
4.2 Electron spin . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2.1 Electrons are indistinguishable . . . . . . . . . . . . . . 59
4.3 Pauli principle . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Hartree–Fock Approximation 63
5.1 Slater-determinant . . . . . . . . . . . . . . . . . . . . . . . . 63
5.1.1 Energy expectation value for Slater determinants . . . 64
5.1.2 Unitary transformation of spin orbitals . . . . . . . . . 66
5.2 Hartree–Fock method . . . . . . . . . . . . . . . . . . . . . . . 66
5.2.1 Solution of the Hartree-Fock equations . . . . . . . . . 70
5.2.2 Orbital energies and Koopmans’ theorem . . . . . . . . 70
5.2.3 Brillouin’s theorem . . . . . . . . . . . . . . . . . . . . 73
6 Molecular Orbital Theory 75

6.1 LCAO-Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.2 Roothaan-Hall equations . . . . . . . . . . . . . . . . . . . . . 76
6.3 Unrestricted open-shell Hartree–Fock . . . . . . . . . . . . . . 79
6.3.1 Spin contamination . . . . . . . . . . . . . . . . . . . . 80
6.3.2 Pople–Nesbet equations . . . . . . . . . . . . . . . . . 81
iv
7 Correlation Energy 85
7.1 Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.2 Matrixelements . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.3 Configuration interaction method (CI) . . . . . . . . . . . . . 88
7.4 Dynamical and non-dynamical electron
correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8 Coupled Cluster Approaches 93

8.1 Size Consistency . . . . . . . . . . . . . . . . . . . . . . . . . 93
8.2 Coupled Cluster Methods . . . . . . . . . . . . . . . . . . . . 95
8.3 Approximate Coupled Cluster methods . . . . . . . . . . . . . 96
8.4 Accuracy in Coupled Cluster calculations . . . . . . . . . . . . 98
9 Møller–Plesset Perturbation Theory 99

9.1 Electron correlation with perturbation
theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
9.2 Convergence of MP series . . . . . . . . . . . . . . . . . . . . 102
9.3 The dissociation problem . . . . . . . . . . . . . . . . . . . . . 103
10 Density Functional Theory: Part I 105

10.1 Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
10.1.1 Hohenberg–Kohn Theorems . . . . . . . . . . . . . . . 105
10.1.2 Kohn–Sham Theory . . . . . . . . . . . . . . . . . . . 108
11 Density Functional Theory: Part II 111

11.1 Exchange and correlation energy . . . . . . . . . . . . . . . . . 111
11.1.1 Density operators . . . . . . . . . . . . . . . . . . . . . 111
11.1.2 Electron-electron interaction . . . . . . . . . . . . . . . 113
11.1.3 One-determinant wave functions . . . . . . . . . . . . . 114
11.1.4 Exchange and correlation energy in the Kohn–Sham
method . . . . . . . . . . . . . . . . . . . . . . . . . . 116
12 Density Functional Theory: Part III 119

12.1 Local density approximation . . . . . . . . . . . . . . . . . . 119
12.2 Gradient corrections . . . . . . . . . . . . . . . . . . . . . . . 120
12.2.1 Exchange functional of Becke (1988) . . . . . . . . . . 121
12.2.2 Correlation functionals of Perdew et al. . . . . . . . . . 121
12.2.3 Lee-Yang-Parr (LYP) correlation functional (1988) . . 121
12.3 Hybrid functionals . . . . . . . . . . . . . . . . . . . . . . . . 122
12.4 Basis sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
12.5 Model calculations . . . . . . . . . . . . . . . . . . . . . . . . 124
v
12.6 Coulomb energy . . . . . . . . . . . . . . . . . . . . . . . . . . 124
12.6.1 Resolution of identity method . . . . . . . . . . . . . . 125
12.6.2 Local expansion of density . . . . . . . . . . . . . . . . 126
12.7 Numerical integration . . . . . . . . . . . . . . . . . . . . . . . 127
12.7.1 Voronoi polyhedra . . . . . . . . . . . . . . . . . . . . 127
12.7.2 Radial integration . . . . . . . . . . . . . . . . . . . . . 129
12.7.3 Angular integration . . . . . . . . . . . . . . . . . . . . 130
13 Molecular Properties 133

13.1 Wavefunction analysis . . . . . . . . . . . . . . . . . . . . . . 133
13.1.1 Basisfunction based population analysis . . . . . . . . . 133
13.1.2 Electrostatic potential derived charges . . . . . . . . . 135
13.1.3 Density based charges . . . . . . . . . . . . . . . . . . 136
13.2 Response properties . . . . . . . . . . . . . . . . . . . . . . . . 136
13.2.1 Derivative techniques . . . . . . . . . . . . . . . . . . . 138
14 NMR Chemical Shielding 141

14.1 General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
14.2 Chemical shifts . . . . . . . . . . . . . . . . . . . . . . . . . . 143
14.3 Linear response for Hartree–Fock . . . . . . . . . . . . . . . . 144
14.3.1 Calculation of response orbitals . . . . . . . . . . . . . 145
vi
Part I
Lectures
1
Lecture 1
Basic Quantum Mechanics
1.1 Operators in quantum mechanics

An observable is a dynamical, measurable variable of a system. In classical
physics observables are described using functions, in quantum mechanics
operators are used. An operator is a symbol for a recipe how an operation
on a function has to be performed.
1.1.1 Eigenfunctions and eigenvalues

In general an operator creates a new function when it is applied to a function.
In special cases the new function is the same as the original one multiplied
by a constant. A function f is an eigenfunction of an operator, if
Ω̂f = ωf ,
where ω is a constant. Such an equation is called an eigenvalue equation. ω
is the eigenvalue of the operator Ω̂.
The set of all eigenfunctions of an operator is a complete basis, this means an
arbitrary function can be written as a linear combination of eigenfunctions.
Ω̂fn = ωn fn .
A function g is expanded in eigenfunctions
X
g= cn f n ,
n
with cn the expansion coefficients and the sum runs over all functions fn . We
can easily calculate the effect of the operator Ω̂ on the function g
X X X
Ω̂g = Ω̂ cn f n = cn Ω̂fn = c n ωn f n
n n n
3
Quantum mechanics can be developed without explicitely defining operators.
However, it is instructive and intuitive to use special forms of operators. The
most common uses the following form for position and momentum operator
(position representation)
h ∂
x −→ x· , px −→ .
2πi ∂x
The special form of the momentum operator follows from the postulates
introduced in the next section.
1.1.2 Commutators
The order in which operators are applied to a function is important. Gener-
ally, we have
ABf 6= BAf ,
or simply AB 6= BA. The quantity AB − BA is called the commutator of A
and B.
[A, B] = AB − BA .
Operators with the same set of eigenfunctions commute
ABg = BAg
X X
A ωnB cn fn
= B ωnA cn fn
n n
X X
ωnA ωnB cn fn = ωnB ωnB cn fn
n n
The opposite is also true; operators that commute can be diagonalized with
the same set of functions.
The commutator for position and momentum operator in the special repre-
sentation introduced above is
h ∂f h ∂
[x, px ]f = (xpx − px x)f = x · − · xf
2πi ∂x 2πi ∂x
h ∂f hf h ∂f
=x· − −x·
2πi ∂x 2πi 2πi ∂x
ih
= ·f
2π
This is true for any function f , therefore we write
ih
[x, px ] =
2π
4
1.1.3 Construction of operators
Many operators we will need can be constructed from the operators for
position and momentum. For example the operator of the kinetic energy
T = p2 /2m is in one dimension
2
p2x 1 ~ ∂ ~2 ∂ 2
T̂ = = =−
2m 2m i ∂x 2m ∂x2
and in three dimensions
( 2 2 2 )
p2x + p2y + p2z 1 ~ ∂ ~ ∂ ~ ∂
T̂ = = + +
2m 2m i ∂x i ∂y i ∂z

~2 ∂2 ∂2 ∂2
= − + +
2m ∂x2 ∂y 2 ∂z 2
~2 2
= − ∇
2m
where ~ = h/2π and the operator ∇2 is the Laplacian.
The potential energy of a particle V (x) is a multiplicative function in the
position representation. This is also the case for the operator of the potential
energy, for example the operator of the Coulomb energy of an electron in the
field of a nuclei with charge Z
Ze2
V̂ = − ·
4π0 r
The operator of the total energy of a system is the Hamilton operator
Ĥ = T̂ + V̂
The general recipe for the construction of quantum mechanical operators is:
1. Write down the classical expression for the observable using the position
and momentum variables of the system.
2. Replace the different functions with their quantum mechanical opera-
tors.
1.2 Postulates of quantum mechanics

1.2.1 Postulates
The following postulates are not complete and in their most strict form, but
they are sufficient for our purposes.
5
Table 1.1: Transition from classical expressions to quantum mechanical op-
erators.
classical quantum mechanical

Position
e.g. x, y x· , y·
Functions of position multiplicative operators
e.g. x2 y x2 y·
or general f(x,y,z) f (x, y, z)·
Momentum
~ ∂
e.g. px i ∂x
Functions of momentum differential operators
∂2
e.g. (px )2 −~2 ∂x 2
~ ∂ ~ ∂ ~ ∂
or general g(px ,py ,pz ) g( i ∂x , i ∂y , i ∂z )
Functions of position and momentum

e.g. ~l = ~r × p~ (x, y, z) × ~
i
∂
, ∂, ∂
∂x ∂y ∂z
Time t t·
∂
Energy E i~ ∂t
Postulate 1
The state of a system is completely characterized by its wave function
Ψ(r1 , r2 , . . . , t).
The positions r1 , r2 , . . . are associated with particles 1, 2, . . . . In general

additional internal coordinates (spin) are needed. From the wave function
all information on the properties of a system can be derived, that are also
accessible by experiments.
Postulate 2
Observables are represented by operators which have the following prop-
erties
[q, pq0 ] = i~δq,q0 , [q, q 0 ] = 0 , [pq , pq0 ] = 0 ,
where q and q 0 are one of the coordinates x, y, z and pq and pq0 are the
corresponding momenta.
6
With this postulate all operators that are functions of position and momen-
tum coordinates are defined. Internal (spin) coordinates have to be defined
separately.
Postulate 3
For a system characterized by its wave function Ψ the average value of
a series of measurements of the observable Ω is given by the expectation
value of the corresponding operator.
The expectation value of an operator Ω̂ for a general state Ψ is defined by

R ∗
Ψ Ω̂Ψdτ
hΩi = R ∗ .
Ψ Ψdτ
If the wave function is normalized

Z
Ψ∗ Ψdτ = 1 ,
we get for the expectation value

Z
hΩi = Ψ∗ Ω̂Ψdτ .
We will assume that state functions are normalized. If Ψ is an eigenstate of

Ω with eigenvalue ω, we have
Z Z Z
hΩi = Ψ Ω̂Ψdτ = Ψ ωΨdτ = ω Ψ∗ Ψdτ = ω .
∗ ∗
A series of measurements on identical systems will give the averaged value

ω. Assuming that Ψ is an eigenstate of the system, but not of Ω̂, we can
expand Ψ in the eigenfunctions of Ω̂
X
Ψ= cn Ψ n , where ΩΨn = ωn Ψn gilt.
n
The expectation value is

Z !∗ ! Z
X X X
hΩi = cm Ψ m Ω cn Ψn dτ = c∗m cn Ψ∗m Ω̂Ψn dτ
m n nm
X Z
= c∗m cn ωn Ψ∗m Ψn dτ .
nm
7
Because the eigenfunctions are an orthonormal basis, the last integral van-
ishes for m 6= n and we get
X Z X X
hΩi = cm cn ωn Ψ∗m Ψn dτ =
∗
c∗n cn ωn = |cn |2 ωn .
nm n n
The expectation value is therefore a weighted sum of the eigenvalues of Ω.

The contribution of a certain eigenvalue is given by the square of the absolute
value of the corresponding expansion coefficient.
Postulate 3’
If Ψ is an eigenfunction of Ω̂ all measurements will give the value ω. Is
Ψ not an eigenfunction of Ω̂ one will get with each measurement a value
ωn , an eigenvalue of Ω̂. The probability to get a certain eigenvalue is
proportional to |cn |2 , where cn is the overlap of Ψ with the eigenfunction
Ψn of Ω̂.
Postulate 4
The probability to find a particle in the volume element dτ at position r
is proportional to |Ψ(r)|2 dτ .
This probability interpretation of the wave function is due to Max Born.

From this it follows that |Ψ(r)|2 is a probability density. The wave function
is a probability amplitude. From this interpretation it also follows that the
wave functions have to be square integrable
Z
|Ψ|2 dτ < ∞ ,
Because the probability to find a particle at any place in space has to be

finite.
Postulate 5
The time evolution of a wave function Ψ(r1 , r2 , . . . , t) is governed by the
Schrödinger equation
∂Ψ
i~ = ĤΨ .
∂t
This is a partial differential equation. Ĥ is the Hamilton operator of the

system. For example the Schrödinger equation for a particle in one dimension
is
∂Ψ ~2 ∂ 2 Ψ
i~ =− + V (x)Ψ .
∂t 2m ∂x2
8
1.2.2 Separation of the Schrödinger equation
Often the Schrödinger equation can be separated into an equation for the
positions and an equation for the time evolution. If the Hamilton operator
does not depend on time, we get in one dimension
~2 ∂ 2 Ψ ∂Ψ
ĤΨ = − 2
+ V (x)Ψ = EΨ = i~
2m ∂x ∂t
For stationary states the energy is a constant and we can solve the ordinary
differential equation using

E
Ψ(x, t) = ψ(x) exp t
i~
Inserting this into the Schrödinger equation we get
Ĥψ = Eψ ,
the time-independent Schrödinger equation which determines the function

ψ(x).
1.3 Hydrogen Atom

The hydrogen atom consists of a proton of mass mp and an electron of mass
me . Reducing the two particle problem to a single particle equation can be
achieved by separation of the center of mass movement:
• 3 center of mass coordinates (translation)
• 3 inner coordinates (relative movement)
This separation leads to the reduced mass for the relative movement:
me mp
µ= = 0.999456me ≈ me
me + m p
~2
Ĥ = − ∇2 + V (r)
2µ
The relative movement can be cast into spherical coordinates (r, θ, φ). The
Coulomb potential for this special case is a central field problem
Ze2
V (r) = − .
r
The nuclear charge is Z = 1 for the hydrogen atom.
9
1.3.1 Angular momentum operators
1.3.2 Definitions
Classical angular momentum:
~l = ~r × p~ = (lx , ly , lz )
lx = ypz − zpy
ly = zpx − xpz
lz = xpy − ypx
Operators in quantum mechanics:

ˆlx = ~ ∂ ∂
y −z
i ∂z ∂y

ˆly = ~ z ∂ − x ∂
i ∂x ∂z

ˆlz = ~ x ∂ − y ∂
i ∂y ∂x
Total angular momentum :
classical:
| ~l |2 = ~l · ~l = lx2 + ly2 + lz2
in quantum mechanics:
ˆl2 = ˆl2 + ˆl2 + ˆl2 = ˆlx ˆlx + ˆly ˆly + ˆlz ˆlz
x y z
Commutators
For two components:

ˆlx ˆly f = −~2 y ∂ − z ∂ z
∂
−x
∂
f
∂z ∂y ∂x ∂z

2 ∂2f ∂f ∂2f 2 ∂ f
2
∂2f
= −~ yz +y − yx 2 − z + zx
∂z∂x ∂x ∂z ∂y∂x ∂y∂z
2 2 2 2

ˆly ˆlx f = −~ zy
2 ∂ f ∂f ∂ f 2 ∂ f ∂ f
+x − xy 2 − z + xz
∂x∂z ∂y ∂z ∂x∂y ∂z∂y
Differences

ˆ ˆ ˆ ˆ 2 ∂f ∂f
(lx ly − ly lx )f = −~ y −x
∂x ∂y

~ ∂f ∂f
= i~ x −y
i ∂y ∂x
= i~ˆlz f
10
Commutator-equations:
[ˆlz , ˆlx ] = i~ˆly
For a single component and l2 :
[ˆlx , ˆl2 ] = lˆx ˆl2 − ˆl2 ˆlx

= ˆlx (ˆlx ˆlx + ˆly ˆly + ˆlz ˆlz ) − (ˆlx ˆlx + ˆly ˆly + ˆlz ˆlz )ˆlx
= (ˆlx lˆy )ˆly + (ˆlx ˆlz )ˆlz − ˆly (ˆly ˆlx ) − ˆlz (ˆlz ˆlx )
= (i~ˆlz + ˆly ˆlx )ˆly + (−i~ˆly + ˆlz ˆlx )ˆlz − ˆly (−i~ˆlz + ˆlx ˆly ) − ˆlz (i~ˆly + ˆlx ˆlz )
= 0
[ˆlx , ˆl2 ] = [ˆly , ˆl2 ] = [ˆlz , ˆl2 ] = 0
From the commutators we learn that it is possible to measure the absolute

value of the angular momentum and the component along one special di-
rection at the same time. However, it is not possible to measure all three
components at the same time.
Central field problems: The potential V (r) is spherical symmetric and
only depends on the distance from the origin. In this case we have
[ˆlx , T̂ ] = [ˆly , T̂ ] = [ˆlz , T̂ ] = 0

[ˆlx , V̂ (r)] = [ˆly , V̂ (r)] = [ˆlz , V̂ (r)] = 0
[ˆlx , Ĥ] = [ˆly , Ĥ] = [ˆlz , Ĥ] = 0
[ˆl2 , Ĥ] = 0
For central field problems Ĥ, ˆl2 and ˆlz have a common set of eigenfunctions.
Every stationary state is characterized by unique expectation values for Ĥ, ˆl2
and ˆlz .
Angular momentum and Hamilton operators in spherical coordi-

nates.
Central field problems are best treated in spherical coordinates (R, θ, φ). The
operators in these coordinates are:
(x, y, z) −→ (R, θ, φ) .
11
We have to consider the rules of calculus, especially the chain rule.
Results (without deduction):
ˆlz = ~ ∂
i ∂φ
2

ˆl = −~
2 2 1 ∂ ∂ 1 ∂
sin θ +
sin θ ∂θ ∂θ sin2 θ ∂φ2
ˆl2
~2 ∂ 2 2 ∂
Ĥ = − + + + V (r)
2m ∂r 2 r ∂r 2mr 2
Therefore  
ˆlz   φ
ˆl2 is a function of θ, φ
 
Ĥ r, θ, φ
Common set of eigenfunctions Ψ(r, θ, φ).
 
φ   ˆlz
θ independent part from eigenvalue problem for ˆl2
 
r Ĥ
Ψ(r, θ, φ) = R(r) Θ(θ) Φ(φ)

| {z }
ˆ
lz
| {z }
l̂2
| {z }
Ĥ
Eigenvalue problem for the operator ˆlz
ˆlz Φ(φ) = ~ ∂ Φ(φ) = lz Φ(φ)

i ∂φ
Eigenfunctions for Φ(φ) can be found easily:
Φ(φ) = Ceimφ constants C, m

ˆlz Φ(φ) = ~ imCeimφ = m~Φ(φ) eigenvalue m~
i
Initial conditions:
Φ(φ) = Φ(φ + 2π) periodicity

Ceimφ = Ceim(φ+2π) = Ceimφ eim2π
eim2π = 1 m integer
12
Normalization:
Z 2π Z 2π Z 2π
∗ 2 −imφ imφ 2
Φ Φdφ = C e e dφ = C dφ = 2πC 2 = 1
0 0 0
−1/2
C = (2π)
Φ(φ) = (2π)−1/2 eimφ
The z-component of the angular momentum can take the values m~ (m in-
teger).
Eigenvalue problem for the operator ˆl2

Eigenfunctions Y (θ, φ) and eigenvalues λ~2

ˆl2 Y (θ, φ) = −~2 1 ∂ ∂ 1 ∂2
sin θ + Y (θ, φ)
= λ~2 Y (θ, φ)
Y (θ, φ) = Θ(θ)Φ(φ)
In the product function the eigenfunctions Φ(φ) of ˆlz are known. We get:
∂2
Φ(φ) = −m2 Φ(φ)
∂φ2

2 1 ∂ ∂ 1 ∂2
−~ sin θ + Θ(θ)Φ(φ) = λ~2 Θ(θ)Φ(φ)
The operator is no longer dependent on φ and therefore we can divide by
Φ(φ). We get a differential equation for Θ. This equation can be solved
using similar methods as used for the harmonic oscillator.
• Substitution x = cos θ
• Power series ansatz
• Recursion formula
• Investigation of convergence, i.e. are there divergencies
• Divergencies can only be avoided at special values of λ where we get a

finite power series.
We get λ = l(l + 1) with integer l ≥ 0 and l ≥ |m|. The square of the
absolute value of the angular momentum can only take values l(l + 1)~2 with
(l = 0, 1, 2, 3, . . . and l ≥ |m|).
13
Angular momentum eigenfunctions Ylm (θ, φ)
Spherical harmonics Ylm (θ, φ) in product form:
Ylm (θ, φ) = Θ(θ)Φ(φ)
The θ dependent part of Θ(θ) is proportional to the Legendre polynomials

|m|
Pl for m = 0 and the associated Legendre functions Pl for m 6= 0.
We get:
x = cos θ
1 dl 2
Pl (x) = (x − 1)l (l = 0, 1, . . .)
2l l! dxl
|m| d|m|
Pl (x) = (1 − x2 )|m|/2 |m| Pl (x) (|m| = 0, 1, 2, . . . , l)
dx
Eigenvalue equation:
ˆl2 Y m (θ, φ) = l(l + 1)~2 Y m (θ, φ) (l = 0, 1, . . .)
l l
ˆlz Y (θ, φ) = m~Y (θ, φ)
m m
(−l ≤ m ≤ l)
l l
Spherical harmonics
1/2
m 2l + 1 (l − |m|)! |m|
Yl (θ, φ) = Pl (cos θ)eimφ
4π (l + |m|)!
Some special spherical harmonics:
1/2
1
Y00 = s
4π
1/2
3
Y10 = cos θ p
4π
1/2
3
Y1±1 = sin θe±iφ p
8π
1/2
5
Y20 = (3 cos2 θ − 1) d
16π
1/2
15
Y2±1 = sin θ cos θe±iφ d
8π
1/2
15
Y2±2 = sin2 θe±2iφ d
32π
14
Linear combinations of complex spherical harmonics Ylm give the real spher-
(m)
ical harmonics Sl .
(
√1 Y −m (θ, φ) + Y m (θ, φ) m>0
(m) 2 l l
Sl (θ, φ) = 1 −m m

√
i 2
Yl (θ, φ) − Yl (θ, φ) m<0
The pair Ylm , Yl−m differs only in the factor eimφ , e−imφ . Using the Euler
formulas:
1
cos mφ = √ e−imφ + eimφ m>0
2
1
sin mφ = √ e−imφ − eimφ m<0
i 2
we get using the definition of Ylm the real functions

(m) |m| cos mφ m>0
Sl (θ, φ) = Nlm Pl (cos θ)
sin mφ m<0
General ansatz:
Ĥ(r, θ, φ)Ψ(r, θ, φ) = EΨ(r, θ, φ)
Ψ(r, θ, φ) = R(r)Ylm (θ, φ)
The angular part of the eigenfunctions is given by the spherical harmonics
Ylm (θ, φ) . With this ansatz and using the eigenvalue equation for ˆl2 we get
the radial Schrödinger equation for the hydrogen atom.

~2 ∂ 2 2 ∂ l(l + 1)~2 Ze2
− + R(r) + R(r) − R(r) = ER(r)
2µ ∂r 2 r ∂r 2µr 2 r
Solution of the radial Schrödinger equation
• Using a new radial function u(r) = rR(r) results in the new equation

~2 d2 u(r) l(l + 1)~2 Ze2
− + − u(r) = Eu(r)
2µ dr 2 2µr 2 r
• Investigation of the asymptotic behavior for r → ∞ and r → 0 leads

to the following form of the solution
u(r) = (κr)l+1 e−κr w(κr)
r
−2µE
κ=
~2
κ is real for E < 0.
15
• Inserting this definition into the radial Schrödinger equation leads to
the Laguerre differential equation for w(κr).
• Solutions of this differential equation can be found with a power series

(recursion formula).
• The power series is not diverging only for special values of E: Quanti-
zation of energy.
Ze2 κ
=k+l+1=n
2|E|
k, l, n integer; k ≥ 0, l ≥ 0, n > l
• Solution for E
Z 2 µe4
En = −
2~2 n2
n ≥ 1 integer, n > l
Energies
Z 2 e2 1
En = − (n = 1, 2, 3, . . .)
2a n2
2
~
a= 2 Bohr Radius = 0.529 Å
µe
There is an infinite number of bound states with En < 0.

E = 0 corresponds to proton and electron infinitely separated.
E > 0 non-bound states with additional kinetic energy.
Ground state of hydrogen atom (Z = 1, n = 1):
e2
E1 = − ≈ −13.6 eV
2a
Excited states of hydrogen atom (Z = 1, En ∝ − n12 ):
E2 ≈ −3.4eV, E3 ≈ −1.5eV, E4 ≈ −0.85eV
Each state Ψnlm is characterized by three quantum numbers, however, the

energy depends only on n.
Energies are n2 times degenerate.
16
n l m Notation States Degeneracy
1 0 0 1s 1 1
2 0 0 2s 1
1 -1,0,1 2p 3 4
3 0 0 3s 1
1 -1,0,1 3p 3
2 -2,-1,0,1,2 3d 5 9
1.3.3 Radial function of hydrogen atom

" 3 #1/2 l
(n − l − 1)! 2Z 2Zr − Zr 2Zr
Rnl (r) = − e na L2l+1
n+l
2n[(n + l)!]3 na na na
Laguerre polynomials
dp
Lp (x) = ex p
xp e−x
dx
Generalized Laguerre polynomials
dq
Lqp (x) = Lp (x)
dxq
Explicite formulas
Zr
ρ=
na
3/2
Z
R10 =2 e−ρ 1s
a
3/2
−1/2 Z
R20 =2 (1 − ρ)e−ρ 2s
a
3/2
−1/2 Z
R21 =6 ρe−ρ 2p
a
3/2
2 Z
R30 = √ (3 − 6ρ + 2ρ2 )e−ρ 3s
9 3 a
3/2
4 Z
R31 = √ (2ρ − ρ2 )e−ρ 3p
9 6 a
3/2
4 Z
R32 = √ ρ2 e−ρ 3d
9 30 a
17
Electron density distribution:
Probability to find an electron within the volume dτ :
Ψ∗nlm Ψnlm dτ = Ψ∗nlm (r, θ, φ)Ψnlm (r, θ, φ)r 2 sin θdrdθdφ
Probability to find an electron within r and r + dr, independent of the angles

(θ, φ):
Z Z
∗ 2
Pnl (r)dr = Rnl Rnl r dr Ylm∗ (θ, φ)Ylm (θ, φ) sin θdθdφ
The integral is one (Normalization). For real Rnl one gets for the radial
probability density
2 2
Pnl (r) = Rnl r
Probability to find an electron within the space angle dΩ (within θ and θ+dθ,
φ and φ + dφ), independent of the distance (r):
Z
Wnl (θ, φ)dΩ = Yl (θ, φ)Yl (θ, φ)dΩ Rnl (r)r 2 dr
m∗ m
The integral is one (Normalization). With the well known formulas for the
spherical harmonics Ylm we get:
2 |m|
Wnl (θ, φ) = Nlm [Pl (cos θ)]2 = Wlm (θ)
The angular distribution of the probability density is for the eigenfunctions

of the hydrogen atom independent of φ and therefore rotational symmetric
w.r.t. the z axis.
18
Unit values in a.u. values in SI-units
Charge elemental charge e=1 e= 1.602189 · 10−19 C
Mass electron mass me = 1 me =9.109534 · 10−31 kg
Angular momentum Planck constant ~ ~=1 ~ = 1.054589 · 10−34 Js

~2
Length Bohr radius a0 = me e 2
=1 a0 = 0.529177 · 10−10 m
me e 4
Energy 2|E1 (H)| with µ = me Hartree= ~2
=1 Hartree = 4.359814 · 10−18 J
19
20
Lecture 2
Basic Mathematical Review
2.1 Linear algebra

2.1.1 Vector spaces
A n-dimensional vector can be represented by its components w.r.t. a set of
orthogonal and normalised basis vectors {~ei }
n
X
~a = ~e1 a1 + ~e2 a2 + · · · + ~en an = ~ei ai
i=1
For a given basis a vector is specified by its components. We represent a

vector  
a1
 a2 
 
a =  .. 
.
an
as a column matrix. The scalar product of two vectors ~a and ~b is defined as

n
X
~a · ~b = a∗i bi
i=1
This is due to ~ei · ~ej = δij where

(
1 if i = j
δij =
0 otherwise
21
2.1.2 Operators
An operator acts on a vector and transforms it into another vector
O~a = ~b
An operator is linear if

O x~a + y~b = xO~a + yO~b
A linear operator is completely determined if its effect on every possible

vector is known. Because any vector can be expressed in terms of the basis,
it is sufficient to know what O does to the basis vectors. O~ei is a vector and
can be expanded in the basis
n
X
O~ei = ~ej Oji i = 1, 2, . . . n
j=1
The numbers Oji form a n × n matrix

 
O11 O12 · · · O1n
 O21 O22 · · · O2n 
 
O =  .. .. .. 
 . . . 
On1 On2 · · · Onn
O is the matrix representation of the operator O in the basis {~ei }.

If A and B are the matrix representations of the operators A and B, the
matrix representation of the operator C = AB can be found as follows:
n
X
C~ej = ~ei Cij
i=1
= AB~ej
X n
=A ~ek Bkj
k=1
n
X
= ~ei Aik Bkj
i,k=1
so that n
X
Cij = Aik Bkj
k=1
22
which is the definition of matrix multiplication and hence
C = AB
The order in which two operators or two matrices are multiplied is crucial.
In general AB 6= BA and AB 6= BA. This means two operators do not
necessarily commute. The commutator of two operators is defined by
[A, B] = AB − BA
and can easily also be written in the matrix representation
[A, B] = AB − BA
2.1.3 Matrices
A general rectangular n × m Matrix A has n rows and m columns
 
A11 A12 · · · A1m
 A21 A22 · · · A2m 
 
A =  .. .. .. 
 . . . 
An1 An2 · · · Anm
If n = m the matrix is square. Two matrices A and B can be multiplied

C = AB if the number of columns of A is equal to the number of rows of B.
M
X
Cij = Aik Bkj
k=1
Matrix-vector multiplication is a special case of matrix-matrix multiplication.

Some important definitions and properties of matrices.
The adjoint of an N × M matrix A, denoted by A† , is a M × N matrix with
elements
(A† )ij = A?ji
i.e. we take the complex conjugate of each of the matrix elements of A and
interchange rows and columns. If the elements of A are real, then the adjoint
of A is called the transpose of A.
The adjoint of the product of two matrices, is equal to the product of the
adjoint matrices in reversed order.
(AB)† = B† A†
23
Proof
(AB)†ij = (AB)?ji
X X
= (Ajk Bki )? = A?jk Bki
?
k k
X
= A†kj Bik
†
= (B A† )ij
†
The adjoint of a column matrix (vector) is a row matrix containing the

complex conjugates of the elements.
a† = (a?1 , a?2 , · · · , a?m )
If a and b are two column matrices with m elements, then

 
b1
 b2  X m
† ? ? ? 
a?i bi

a b = (a1 , a2 , · · · , am )  ..  =
 .  i=1
bm
This is the scalar product of two vectors.

Special properties of square matrices
• A matrix A is diagonal if all its off-diagonal elements are zero
Aij = Aii δij .
• The trace of A is the sum of its diagonal elements

X
trA = Aii .
i
• The unit matrix is defined by (1)ij = δij and
1A = A1 = A
• The inverse of A, denoted A−1 is a matrix such that
A−1 A = AA−1 = 1
• A unitary matrix A is one whose inverse is its adjoint, a real unitary

matrix is called orthogonal
A−1 = A†
24
• A Hermitian matrix is self-adjoint, a real Hermitian matrix is called
symmetric
A† = A
Further properties
• trAB = trBA
• (AB)−1 = B−1 A−1
• If U is unitary and B = U† AU, then A = UBU†
• If A, B, C Hermitian matrices, and C = AB, then A and B commute.
2.1.4 Determinants
The determinant of an n × n matrix A is a number obtained from

A11 · · · A1n
X n!
.. .
.. = (−1)pi Pi A11 A22 · · · Ann

det(A) = |A| = .

An1 · · · Ann i=1
where Pi is a permutation operator that permutes the column indices 1, 2, 3, . . . , n

and the sum runs over all n! permutations of the indices; pi is the number of
transpositions required to restore a given permutation to natural order.
Important properties of determinants
• If each element in a row or in a column is zero the value of the deter-
minant is zero.
Q
• If (A)ij = Aii δij , then |A| = i Aii = A11 A22 . . . Ann .
• A single interchange of any two rows (or columns) of a determinant
changes its sign.
• |A| = (|A† |)?
• |AB| = |A||B|
• If any two rows (or columns) of a determinant are equal, the value of
the determinant is zero.
• |A−1 | = (|A|)−1
• If AA† = 1, then |A|(|A|)? = 1.
• If U† OU = Ω and U† U = UU† = 1, then |O| = |Ω|.
25
2.1.5 Dirac notation
Dirac introduced a notation for vector spaces that allows to express results
in a concise and simple manner. We consider n basis vectors denoted by the
symbol |ii, i = 1, 2, . . . , n, which are called ket vectors or simply kets. We
assume this basis is complete so that any ket vector |ai can be written as
n
X
|ai = |ii ai .
i=1
The column matrix a is the matrix representation of the abstract vector |ai
in the basis {|ii}. We introduce an abstract bra vector ha| whose matrix
representation is a† . The scaler product between a bra ha| and a ket |bi is
defined as
 
b1
 b2  X n
† ? ? ?  
ha||bi ≡ ha | bi = a b = (a1 , a2 , . . . , an )  ..  = a ? bi
 .  i=1 i
bn
The square of the length of a vector

n
X n
X
ha | ai = a?i ai = |ai |2
i=1 i=1
is always real and positive. In analogy to the ket basis we can introduce
the bra basis {hi|} that is complete in the sense that any bra vector can be
written as a linear combination of the bra basis vectors.
X
ha| = a?i hi| .
i
The scalar product now becomes

X
ha | bi = a?i hi | jibj .
ij
For this to be identical with the previous definition we must have that
hi | ji = δij
which is a statement of the orthonormality of the basis.
26
For a given ket |ai or bra ha| the components of its matrix representation
with respect to a basis {|ii or {hi|} can be determined by multiplication from
the left or right with a general basis vector |ji or hj|
X X
hj | ai = hj | ii ai = δji ai = aj
i i
X X
ha | ji = a?i hi | ji = a?i δij = a?j
i i
Using these results we can write

X X
|ai = |ii ai = |iihi | ai
i i
and X X
ha| = a?i hi| = ha | iihi|
i i
which suggests we write X
1= |iihi|
i
which is a statement of the completeness of the basis.
An operator O is defined as an entity which when acting on a ket |ai converts
it into a ket |bi
O | ai = |bi .
The operator is completely determined by its action on all basis vectors
X X
O | ii = Oji | ji = (O)ji | ji
j j
so that O is the matrix representation of the operator O in the basis {|ii}.

The matrix representation is found by multiplying on the left by hk|
X X
hk | O | ii = hk | ji(O)ji = δkj (O)ji = (O)ki
j j
which provides a useful expression for the matrix elements of O.
2.1.6 Change of basis

We specify two different basis sets, one, {|ii} where we use latin letters
i, j, k, . . . for the vectors and another {|αi} where we use Greek letters α, β, γ, . . ..
Thus we have X
hi | ji = δij , |iihi| = 1
i
27
and X
hα | βi = δαβ , |αihα| = 1 .
α
Since the basis {|ii} is complete we can express any ket in the basis {|αi} as
a linear combination of kets in the basis {|ii} and vice versa. That is,
X X X
|αi = 1 | αi = |iihi | αi = |iiUiα = |ii(U)iα
i i i
where we have defined the elements of a transformation matrix U as
hi | αi = Uiα = (U)iα .
Transforming in the opposite direction, we have

X X X
?
|ii = 1 | ii = |αihα | ii = |αiUiα = |αi(U† )αi
α α α
where we have used
hα | ii = hi | αi? = Uiα
?
= (U† )αi .
As a consequence of the orthonormality of the bases the transformation ma-

trix U is unitary
X
δij = hi | ji = hi | αihα | ji
α
X
= Uiα (U† )αj = (UU† )ij
α
which is in matrix notation

1 = UU† .
Starting from δαβ = hα | βi one can show that
1 = U† U
and hence U is unitary. Thus we arrive at the important result that two
orthonormal bases are related by a unitary matrix. The elements of the
transformation matrix U are scalar products between the two bases.
We now investigate the matrix representations of an operator for the two
different bases. Suppose O is the matrix representation of O in the basis
{|ii}, while Ω is its matrix representation in the basis {|αi}
Oij = hi | O | ji ; Ωαβ = hα | O | βi .
28
We find the relation between O and Ω from
Ωαβ = hα | O | βi = hα | 1 O 1 | βi
X
= hα | iihi | O | jihj | βi
ij
X
= (U† )αi (O)ij (U)jβ .
ij
Thus
Ω = U† OU
or, multiplying on the left by U and on the right by U†
O = UΩU† .
These equations show that the matrices O and Ω are related by a unitary
transformation. The importance lies in the fact that for any operator with a
Hermitian matrix representation which is not diagonal in a given basis {|ii},
it is always possible to find a basis {|αi} in which the matrix representation
of the operator is diagonal, i.e.
Ωαβ = ωα δαβ .
2.1.7 Eigenvalue problem

If the result of the operation O | αi is simply a constant times |αi, i.e.
O | αi = ωα |αi
then we say that |αi is an eigenvector of the operator O with an eigenvalue ωα .

Without loss of generality we can choose the eigenvectors to be normalised
hα | αi = 1
• The eigenvalues of a Hermitian operator are real.

This follows from
hα | O | αi = hα | O † | αi = hα | O | αi?
Multiplying the eigenvalue relation by hα| and using the above relations
we have
ωα = ωα?
which is the required result.
29
• The eigenvectors of a Hermitian operator are orthogonal.
We consider
O | βi = ωβ |βi
and its adjoint
hβ | O † = hβ|ωβ?
Since O is Hermitian and ωβ is real, we have
hβ | O = hβ|ωβ
Multiplying O | αi = ωα |αi by hβ| and hβ | O = hβ|ωβ by |αi and

subtracting the results we obtain
(ωβ − ωα )hβ | αi = 0
so that hβ | αi = 0 if ωα 6= ωβ . Thus orthogonality follows imme-

diately if the two eigenvalues are not the same (nondegenerate). It
can be shown, that degenerate eigenvectors can always be chosen to be
orthogonal (we omit this here).
The eigenvalue problem can be posed as follows. Given the matrix represen-
tation, O, of a Hermitian operator O in the orthonormal basis {|ii}, we wish
to find the orthonormal basis {|αi} in which the matrix representation, Ω of
O is diagonal. In other words we wish to diagonalise the matrix O. The two
representations of the operator O are related by an unitary transformation
Ω = U† OU
Thus the problem of diagonalising the Hermitian matrix O is equivalent to

the problem of finding the unitary matrix u that converts O into a diagonal
matrix  
ω1
 ω2 0 
U† OU = ω = 
 
..  .
 0 . 
ωn
The problem can be reformulated as follows. Given the n × n Hermitian
matrix O, find all n distinct column vectors c (the eigenvectors) and corre-
sponding n numbers ω (the eigenvalues) such that
Oc = ωc .
This equation can be rewritten as
(O − ω1) c = 0 .
30
Nontrivial solutions (c 6= 0) only exist when
|O − ω1| = 0 .
This determinant is a polynomial of degree n in the unknown ω. It has n

roots ωα which in this case will be the eigenvalues of the matrix O. Once
we have found the eigenvalues, we can find the corresponding eigenvectors
by substituting each ωα into the above equation and solving the resulting
equation for cα . In this way cα can be found within a multiplicative constant,
which is finally determined by requiring the eigenvectors to be normalised
X
(cαi )? cαi = 1 .
i
In this way we find n solutions to
Ocα = ωα cα α = 1, 2, . . . , n
Since O is Hermitian, the eigenvalues are real and the eigenvectors are or-
thogonal. We construct a matrix U defined as Uiα = cαi , i.e.
 
c11 c21 · · · cn1
 c1 c2 · · · cn 
 2 2 2 1 2 n
U =  .. .. ..  = (c , c , · · · , c )
. . .
c1n c2n · · · cnn
Thus the αth colum of U is just the column matrix cα . Then from the
orthogonality relation we find
X X
?
Uiα U iβ = (U† )αi (U)iβ = δαβ
i i
which in matrix notation is

U† U = 1
Multiplying OU = Uω by U† we have
U† OU = ω
which shows that U is the unitary transformation between the original and
new basis for the diagonal representation of O.
31
2.1.8 Functions of matrices
Given a Hermitian matrix A, we can define a function of A, i.e. f (A), in
much the same way we define a function f (x) of a simple variable x. For
example, the square root of a matrix A, which we denote by A1/2 , is simply
that matrix which when multiplied by itself gives A
A1/2 A1/2 = A
The sine or the exponential of a matrix are defined by the Taylor series of
the function
1 1 1
exp(A) = 1 + A + A2 + A3 + · · ·
1 2 3
or in general
∞
X
f (A) = c k Ak
k=0
We are still faced with the problem of calculating A1/2 or exp(A). If A is a

diagonal matrix
(A)ij = ai δij
everything is simple, since
 
ak1
 ak1 0 
Ak = 
 
.. 
 0 . 
ak1
so that
P 
k
k ck a 1 P
∞ k
k ck a 2 0
X  
k  
f (A) = ck A =  .. 
k=0
 0 . 
P k
k ck a n
 
f (a1 )

 f (a2 ) 0 

= .. 
 0 . 
f (an )
Similarly, the square root of a diagonal matrix is
 1/2 
a1
1/2
 a2 0 
A1/2 = 
 
. ..

 0 
1/2
an
32
If A is Hermitian we can always find a unitary transformation that diago-
nalises it
U† AU = a
The reverse transformation is
A = UaU†
Now notice that

A2 = UaU† UaU† = Ua2 U†
or in general
Ak = Uak U†
so that
∞
!
X X
f (A) = c k Ak = U c k ak U† = U f (a) U†
k=0 k
 
f (a1 )

 f (a2 ) 0 
 †
= U .. U
 0 . 
f (an )
Thus to calculate any function of a Hermitian matrix A, we first diagonalise

A to obtain a, the diagonal matrix containing all the eigenvalues of A. We
then calculate f (a), which is easy because a is diagonal. Finally we transform
f (a) back using the matrix U formed by the eigenvectors of A. For example,
we can find the square root of a matrix A as
A1/2 = Ua1/2 U†
since
A1/2 A1/2 = Ua1/2 U† Ua1/2 U† = Ua1/2 a1/2 U† = UaU† = A
2.2 Perturbation theory

For an un-perturbed model system with the Hamiltonian Ĥ (0) the exact wave
(0) (0)
functions Ψn and exact energies En are known.
Ĥ (0) Ψ(0) (0) (0)

n = E n Ψn
33
For the system of interest, a perturbed system with the Hamiltonian Ĥ, the
wave functions Ψn and energies En cannot be calculated exactly.
ĤΨn = En Ψn
Perturbation theory describes how the exact solutions of the model system
change under the influence of a perturbation λĤ (1) (perturbation parameter
λ, perturbation operator Ĥ (1) ). In a systematic procedure corrections to the
energy and the wave functions in first, second , third, . . . order are defined:
En(1) , En(2) , En(3) , . . . and Ψ(1) (2) (3)
n , Ψn , Ψn , . . .
The perturbation λĤ (1) has to be small and we have

Ĥ = Ĥ (0) + λĤ (1)
)
(0)
En −→ En
(0) for λĤ (1) −→ 0
Ψn −→ Ψn
2.2.1 Derivation
Ansatz:
E = E (0) + λE (1) + λ2 E (2) + · · ·
Ψ = Ψ(0) + λΨ(1) + λ2 Ψ(2) + · · ·
0 = hΨ(0) | Ψ(p) i with p = 1, 2, . . .
Formal power series in the perturbation parameter λ. Orders of magnitude:
|E (0) | |λE (1) | |λ2 E (2) | . . .
Insert in the Schrödinger equation:

(0) (1)
Ĥ + λĤ Ψ(0) + λΨ(1) + λ2 Ψ(2) + · · · =

E (0) + λE (1) + λ2 E (2) + · · · Ψ(0) + λΨ(1) + λ2 Ψ(2) + · · ·
Combine w.r.t. the orders of powers of λ:
h i
Ĥ (0) Ψ(0) − E (0) Ψ(0) +
h i
(0) (1) (1) (0) (0) (1) (1) (0)
λ Ĥ Ψ + Ĥ Ψ − E Ψ − E Ψ
h
+ λ2 Ĥ (0) Ψ(2) + Ĥ (1) Ψ(1) − E (0) Ψ(2)

−E (1) Ψ(1) − E (2) Ψ(0) + · · · = 0
34
For an arbitrary perturbation parameter λ each order has to be zero.
Ĥ (0) Ψ(0) =E (0) Ψ(0)
Ĥ (0) Ψ(1) + Ĥ (1) Ψ(0) =E (0) Ψ(1) + E (1) Ψ(0)
Ĥ (0) Ψ(2) + Ĥ (1) Ψ(1) =E (0) Ψ(2) + E (1) Ψ(1) + E (2) Ψ(0)
..
.
Ĥ (0) Ψ(p) + Ĥ (1) Ψ(p−1) =E (0) Ψ(p) + E (1) Ψ(p−1) + · · · + E (p) Ψ(0)
The solution for the un-perturbed problem is known. For the calculation of
perturbative corrections the equations above have to be solved.
∗
Solutions: Equations are multiplied from left with Ψ(0) and then integrated
over all coordinates.
2.2.2 Energy corrections

First order:
hΨ(0) | Ĥ (0) Ψ(1) i +hΨ(0) | Ĥ (1) Ψ(0) i = E (0) hΨ(0) | Ψ(1) i +E (1) hΨ(0) | Ψ(0) i
| {z } | {z } | {z }
=0 =0 =1
The right hand side can be simplified using that the unperturbed wave func-
tion Ψ(0) is normalised and orthogonal to the corrections Ψ(p) . The first term
on the left hand side is zero because Ĥ (0) is Hermitian and the orthogonality
condition.
hΨ(0) | Ĥ (0) Ψ(1) i = hĤ (0) Ψ(0) | Ψ(1) i = E (0) hΨ(0) | Ψ(1) i = 0
Therefore we get
E (1) = hΨ(0) | Ĥ (1) Ψ(0) i
The energy correction in first order is equal to the expectation value of the
perturbation operator w.r.t. the un-perturbed wave functions.
Higher orders:
From using the same derivation we get:
E (2) = hΨ(0) | Ĥ (1) Ψ(1) i
..
.
E (p) = hΨ(0) | Ĥ (1) Ψ(p−1) i
The wavefunction in (p − 1)-th order determines the energy in p-th order
within perturbation theory.
Note: A more thorough analysis shows that the wave function in p-th order
is sufficient to calculate the energy up to (2p + 1)-th order.
35
(1) (2)
2.2.3 Corrections Ψn and En
(1)
Ψn is expanded in a complete set of solutions to the un-perturbed problem
(0)
Ψj .
X (0)
Ψ(1)
n = aj Ψj
j6=n
(0) (1)
With this orthogonality of Ψn w.r.t. Ψn is guaranteed.
X (0)
hΨ(0) (1)
n | Ψn i = aj hΨ(0)
n | Ψj i = 0
j6=n
From the general equation in first order we get:
Ĥ (1) Ψ(0) (1) (0) (0) (1)

n − En Ψn = En Ψn − Ĥ
(0) (1)
Ψ
Xn
(0)
= En(0) − Ĥ (0) aj Ψj
j6=n
X
(0) (0)
= aj En(0) − Ej Ψj
j6=n
(0) ∗
Now we use the usual procedure: Multiplication from left with Ψi (i 6= n)
and integration.

(0) (0) (1)
hΨi | Ĥ (1) Ψ(0)
n i = a i E (0)
n − E i = Hin
(1)
This defines the matrix element Hin of the perturbation operator. With this
definition we get:
(1)
Hin
ai = (0) (0)
En − E i
(1)
X Hin (0)
Ψ(1)
n = (0) (0)
Ψi
i6=n En − Ei
Energy correction in second order:

(1)
(0)
X Hin (0)
En(2) = hΨi | Ĥ (1) Ψ(1)
n i = (0) (0)
hΨi | Ĥ (1) Ψ(0)
n i
i6=n En − Ei
(1)
X |Hin |2
= (0) (0)
i6=n En − E i
36
2.2.4 Summary
(0)
The equations of perturbation theory contain the energies En and wave
(0) (1)
functions Ψn of the un-perturbed problem and the matrix elements Hin of
the perturbation operator:
(1) (0)
Hin = hΨi | Ĥ (1) Ψ(0)
n i
Perturbation series of the energy:

(1)
X |Hin |2
En = En(0) + (1)
Hnn + (0) (0)
+···
i6=n En − E i
Perturbation series of the wave function:

(1)
X Hin (0)
Ψn = Ψ(0)
n + (0) (0)
Ψi + · · ·
i6=n En − Ei
Higher order corrections are not used very often.
2.3 The variation method

The expectation value of the energy EΨ calculated with an arbitrary (valid)
wave function Ψ is an upper bound for the exact energy E0 of the ground
state of the system.
EΨ = hΨ | ĤΨi ≥ E0 Ψ normalised
Proof: Expand Ψ in the complete set of the exact eigenfunctions Ψn of the

system.
ĤΨn = En Ψn
X∞
Ψ= an Ψn
n=0
Normalisation:
∞ X
X ∞
hΨ | Ψi = a∗n am hΨn | Ψm i
n=0 m=0
∞
X
hΨ | Ψi = |an |2 = 1
n=0
37
Expectation value of the energy:
∞ X
X ∞
EΨ = a∗n am hΨn | ĤΨm i
n=0 m=0
∞ X
X ∞
= a∗n am Em hΨn | Ψm i
n=0 m=0
∞
X
= |an |2 En
n=0
In the last equation we replace all eigenvalues En with the lowest eigenvalue
E0 (ground state), with E0 < E1 < E2 . . .. Then we have
∞
X
EΨ ≥ E 0 |an |2 = E0
n=0
2.3.1 The variation principle

Within all valid wave functions Ψ for the ground state we are looking for the
function with the lowest expectation value for the energy EΨ .
EΨ −→ Minimum
Practical procedure: Choose an Ansatz for Ψ with free parameters αi and

determine the parameters using the variational principle. This gives the best
possible solution within this Ansatz.
This is equivalent to a extremal value problem. The variation δEΨ has to be
zero for any arbitrary variation of the parameters αi .
∂EΨ ∂EΨ
δEΨ (α1 , α2 , . . .) = dα1 + dα2 + · · · = 0
∂α1 ∂α2
Because the variations dαi are arbitrary, we have:
∂EΨ
=0 for all i
∂αi
2.3.2 Linear variations

Variational function Ψ is chosen as a linear combination of n fixed basis
functions Φi :
Xn
Ψ= ci Φi
i=1
38
The coefficients ci are varied until the minimal value of the energy expectation
value EΨ is found. For not normalised functions Ψ we have:
hΨ | ĤΨi A
EΨ = =
hΨ | Ψi B
XX
∗
A= ck cj hΦk | ĤΦj i
k j
XX
B= c∗k cj hΦk | Φj i
k j
Notation:
∗
Hkj = hΦk | ĤΦj i =Hjk (Matrix element of Ĥ)
∗
Skj = hΦk | Φj i =Sjk (Overlap integral)
Minimisation of EΨ w.r.t. the variational parameters ci :
∂EΨ
=0 for i = 1, 2, . . . , n
∂ci

∂EΨ 1 ∂A ∂B 1 ∂A A ∂B
= 2 B −A = − =0
∂ci B ∂ci ∂ci B ∂ci B ∂ci
Because B 6= 0 (Normalisation integral) and with EΨ = A/B we get:
∂A ∂B
− EΨ =0
∂ci ∂ci
Calculation of the derivatives:
XX
A= c∗k cj Hkj
k j
∂A X
=2 cj Hij
∂ci j
XX
B= c∗k cj Skj
k j
∂B X
=2 cj Sij
∂ci j
Inserting the derivatives we get the generalised eigenvalue equation:

X
(Hij − EΨ Sij ) cj = 0 for i = 1, 2, . . . , n
j
39
Known: Integrals Hij and Sij
Unknown: Coefficients cj and corresponding energies EΨ .
Linear algebra: Non-trivial solution of a homogeneous system of linear equa-
tions, for the case of vanishing determinant:
|Hij − EΨ Sij | = 0
The solution of this equation gives n energies Ek , where the lowest value is
the variational energy of the ground state. For each energy Ek we get after
inserting in the system of linear equations the coefficients cjk and therefore
the corresponding wave function Ψk .
Additional index to distinguish the individual solutions:
n
X
Ψk = cjk Φj
j=1
X X
Hij cij = Sij cjk Ek for i, k = 1, 2, . . . , n
j j
Using matrices:
H = {Hij }, S = {Sij }, C = {Cij }, E = {Ei δij }

HC = SCE
Using the n × n -matrices H, C, S and E we can write the equation in matrix

form. H and S are real for real basis functions Φj and symmetric (hermitian
otherwise), E is a real diagonal matrix.
Important special case: Orthonormal basis
Sij = δij
HC = CE
H is the matrix representation of the Hamilton operator in the orthonormal

basis. Through diagonalisation of H we get the energy eigenvalues (diagonal
elements of E) and eigenvectors (columns of C). We have
E = C −1 HC = C T HC
n
X
T
C C=I or |cjk |2 = 1 for k = 1, 2, . . . , n
j=1
The inverse and transpose of the matrix C are identical. Eigenvectors are
orthonormal.
40
Lecture 3
Molecular Hamiltonian
3.1 Time-independent non-relativistic Schrö-

dinger equation
For a system with N electrons and K nuclei the time-independent Schrödinger
equation is
~ n (~r, ~s, R,
Ĥ(~r, R)Ψ ~ ~σ ) = En Ψn (~r, ~s, R,
~ ~σ ) .
Notation:
~ri = (xi , yi , zi ) Coordinates for electron i
~si Spin - coordinate for electron i
~ A = (xA , yA , zA )
R Coordinates for nuclei A
~σA Spin - coordinate for nuclei A
~r = (~r1 , . . . , ~rN ) Coordinates for all electrons
~s = (~s1 , . . . , ~sN ) Spin - coordinates for all electrons
~ = (R
R ~ A, . . . , R ~ K) Coordinates for all nuclei
~σ = (~σA , . . . , ~σK ) Spin - coordinates for all nuclei
Conventions:
Electrons : Number N ; Indices i, j, . . .; Mass m
Nuclei: Number K ; Indices A, B, . . .; Mass MA
Non-relativistic Hamilton operator
~ = T̂N (R)
Ĥ(~r, R) ~ + T̂el (~r) + V̂Ne (~r, R)
~ + V̂ee (~r) + V̂NN (R)
~
41
Kinetic energy:
X ~2 ∂ 2 ∂2 ∂2
X ~2
~ =−
T̂N (R) + 2 + 2 =− ∇2A
2
A
2M A ∂x A ∂yA ∂zA A
2M A

~2 X ∂ 2 ∂2 ∂2 ~ 2 X
T̂el (~r) = − + 2+ 2 =− ∇2i
2m i ∂x2i ∂yi ∂zi 2m i
Potential energy:
~ =−
XX e2
V̂Ne (~r, R) ZA
A i
rAi
X e2
V̂ee (~r) = +
r
i<j ij
~ =+
X e2
V̂NN (R) ZA ZB
A<B
RAB
Where ZA is the charge of nuclei A and rij = |ri − rj |.

Ĥ includes the non-relativistic kinetic energy and the potential energy of the
Coulomb interaction.
Not included are:

• Relativistic correction to the kinetic energy
• Interaction of magnetic moments

(orbit/orbit, Spin/orbit, Spin/Spin)
• Interaction with external electric and magnetic fields
3.2 Born–Oppenheimer Approximation

3.2.1 Part I
Goal: Separation of nuclear motion and electron motion.
Idea: Electrons are much lighter than nuclei and move therefore much faster.
They adjust to each new position of the nuclei almost instantaneously.
~ is calculated for fixed posi-
Ansatz: The electronic wave function Ψel (~r; R)
~ of the nuclei. It directly depends on ~r and parametrically on
tion (R)
~
R.
42
Electronic Schrödinger equation
~ el (~r; R)
Ĥel (~r; R)Ψ ~ = Eel (R)Ψ
~ el (~r; R)
~
with the electronic Hamilton operator

~ = T̂el (~r) + V̂Ne (~r; R)
Ĥel (~r; R) ~ + V̂ee (~r) .
For the derivation we have used

~ = 0 non-moving nuclei
• T̂N (R)
~ does not act on ~r (multiplicative operator) and can be neglected
• V̂NN (R)
for the solution
Interpretation of the equations
• The calculation of Ψel takes into account the electronic interaction of
nuclei and electrons, however not their dynamic coupling. There is no
exchange of kinetic energy.
• Macroscopic analogy: No heat exchange, adiabaticity.
• The solution of the electronic Schrödinger equation is difficult but pos-

sible (Quantum chemistry).
~ are a complete orthonormal set of func-
• The eigenfunctions Ψnel (~r; R)
tions.
3.2.2 Part II
~ is expanded in the electronic wave functions
The total wave function Ψ(~r, R)
n ~
Ψel (~r; R) X
~ =
Ψ(~r, R) ~ nel (~r; R)
χn (R)Ψ ~ .
n
Insertion in the Schrödinger equation

X X
~
Ĥ(~r, R) ~ n (~r; R)
χn (R)Ψ ~ =E ~ n (~r; R)
χn (R)Ψ ~
el el
n n
Multiplication from left with Ψm∗ ~ and integration over all electronic
el (~
r; R)
coordinates ~r
h i X
~ + Eelm (R)
T̂N (R) ~ + V̂NN (R)
~ − E χm (R)~ =− ~ n (R)
Ĉmn (R)χ ~ ,
n
43
and
X ~2
~ =−
Ĉmn (R) hΨm 2 n m n
el | ∇A | Ψel i + 2hΨel | ∇A | Ψel i · ∇A .
A
2MA
~ are due to the action of
Note for the derivation: The coupling terms Ĉmn (R)
~ n ~
T̂N (R) on Ψel (~r; R).
Comparison of relevant matrix ~
h elements shows that Ĉmni(R) is normally van-
~ + E m (R)
ishingly small compared to T̂N (R) ~ + V̂NN (R)
~ .
el
Possible approximations
Ĉmn = 0 for m 6= n : Adiabatic approximations
Ĉmn = 0 for all m, n : Born–Oppenheimer approximation
Both approximations lead to non-coupled systems of equations for the nuclear
~
functions χn (R).
3.2.3 Part III

After the decoupling of the nuclear Schrödinger equation we get for each
electronic state an effective Schrödinger equation for the nuclear motion:
~ m (R)
ĤN (R)χ ~ = Eχm (R)
~ ,
with
~ = T̂N (R)
ĤN (R) ~ + V̂eff (R)
~ ,
where we have been using the adiabatic approximation
~ = E m (R)
V̂eff (R) ~ + V̂NN (R)
~ + Ĉmm (R)
~
el
and the Born–Oppenheimer approximation, respectively.

~ = E m (R)
V̂eff (R) ~ + V̂NN (R)
~
el
~ is small and difficult to calculate. Most

The diagonal correction term Ĉmm (R)
calculations are done using the Born–Oppenheimer potential energy surface:
~ = E m (R)
U (R) ~ + VNN (R)
~
el
~ and the interaction energy of the nuclei are

The electronic energy Eelm (R)
independent of the nuclear masses, therefore we see that
~
Isotopes have the same potential energy surface U (R).
The Born–Oppenheimer approximation is often valid to high accuracy. Ex-

ceptions occur if different electronic states are significantly coupled through
nuclear motion. Special care is necessary for
44
• close lying potential energy surfaces (crossings and avoided crossings
of electronically excited states)
• rapid movements of nuclei (impact of fast particles).
3.2.4 Order of magnitude of effects for H2

a) Relativistic effects in H2 (State X 1 Σ+
g ) at the equilibrium geometry Re
−1
Change of D0 : −0.48 cm
D1 : −0.34 cm−1
E(ν = 1) - E(ν = 0): −0.14 cm−1
b) Non-adiabatic effects in H2 (State X 1 Σ+

g ) at the equilibrium geometry
Re
Approximately calculated change in energy
H2 : 0.65 cm−1 , HD : 0.42 cm−1 , D2 : 0.23 cm−1
All effects are certainly smaller than 1 cm−1 ,

D0 = 36118.3 ± 0.5 cm−1 (experiment)
D0 = 36118.0 cm−1 (theory; adiabatic + relativistic)
c) Born–Oppenheimer effects, more precisely: energy change due to the di-

agonal adiabatic correction Cnn
H2 X 1 Σ +
g 120 cm−1
a 3 Σ−
g 100 cm−1
3 +
B Σu 200 cm−1
E, F 1 Σ+
g 430 cm−1 (avoided crossing)
This could be important for the exact calculation of electron spectra
(80 cm−1 ≈ 0.01eV).
Not important for vibrational spectroscopy, Born–Oppenheimer effects
on calculated transition energies (∆ν = 1) are smaller than 1 cm −1 for
H2 (State X 1 Σ+
g ) for ν = 1 − 10.
3.3 Stationary points

N Atoms, 3N degrees of freedom (coordinates)
For non-linear molecules:
• 3 degrees of freedom for translation
• 3 degrees of freedom for rotation
45
• 3N − 6 degrees of freedom for internal coordinates
Coordinates (e.g. Cartesian coordinates):
~ = (R1 , R2 , . . . , R3N ) = {Ri }
R
Definition
~
∂U (R)
gi = components of gradient
∂Ri
2 ~
∂ U (R)
fij = harmonic force constants
∂Ri ∂Rj
Taylor expansion within the vicinity of a reference point with coordinates
~ 0 = {Ri,0 }:
R
X 1 XX
U = U (0) + g i xi + fij xi xj + · · ·
i
2 i j
with
xi = Ri − Ri,0
or in matrix form
1
U = U (0) + g T x + xT F x + · · ·
2
with
   
x1 g1  
F1,1 · · · F1,3N
 x2   g2 
F =  ... .. .. 
    
x= .. , g= .. , . . 
 .   . 
F3N,1 · · · F3N,3N
x3N g3N
Definition:
Points with the coordinates R ~ e = {Ri,e } are called stationary if all compo-
nents of the gradient are zero.
~ e) = 0
g i (R for i = 1, 2, . . . , 3N
At a stationary point we have the Taylor expansion

1 XX
U = U (0) + fij xi xj + · · ·
2 i j
46
with
xi = Ri − Ri,e .
Characterisation of stationary points:

The force constant matrix F (R~ e ) has 3N real eigenvalues, whereof
• 3 for translation (equal zero)
• 3 for rotation (equal zero)
• 3N − 6 for internal motion (different from zero)
The number of negative eigenvalues is used to characterise the stationary

point
• 0 for a minimum
• 1 for a transition state
• 2, 3, . . . for a saddle point of higher order
3.4 Internal coordinates

Changes in bond length and in the angles between chemical bonds provide
the most significant and physically meaningful set of coordinates for the
description of the potential energy. These coordinates are called ”internal
coordinates” since they describe (being unaffected by translation and rota-
tions of the molecule as a whole) just the internal motions of the molecule,
i.e. the molecular vibrations.
The types of internal coordinates which are generally used are the following:
• Bond stretching coordinate : variation of the length of a chemical bond.
r
A B
• In-plane bending coordinate : variation of the angle between two chem-

ical bonds having one atom in common.
47
A B
α
• Out-of-plane bending coordinate : variation of the angle between the

plane defined by two bonds with one atom in common and a third bond
connected to the common

atom.

Α

α

Β

D

C
• Torsion coordinate : variation in the dihedral angle between the planes

determined by three consecutive bonds connecting four atoms.
48
Α
C
B
D
Α
D
α
B C
Since the variations in bond length and angles in a molecule can be written
in terms of the variation in the Cartesian displacement coordinates of the
atoms, the problem we have to face is that of writing the explicit form of
the transformation from Cartesian coordinates to the new set of internal
coordinates.
Let us call sk a generic internal coordinate. The most general relation be-
tween sk and the Cartesian displacement coordinates can be written in the
form X 1X k
sk = Bik xi + B xi xj + higher terms
i
2 ij ij
where the coefficients Bik , Bijk , etc., are determined by the molecular geome-
try. A drastic simplification can be achieved if we restrict our treatment to
the case of infinitesimal amplitudes of vibration where we can drop all terms
not linear in x.
If we call s and x the vectors whose components are the internal and the
Cartesian coordinates, the linear equation in matrix notation is
s = Bx
where B is a matrix whose elements are the coefficients Bik .

The matrix B is not, in general, a square matrix. Let us consider, for example
a bent triatomic molecule and choosing as internal coordinates the two bond
length and the bond angle. There are also 3n − 6 = 3 vibrational degrees
of freedom and with the 9 Cartesian coordinates the B-matrix will have 3
rows and 9 columns. The translations and rotations are not included in the
set of internal coordinates which describe, by definition a molecular motion
49
in which the relative positions of the atoms are changed. Since B is not a
square matrix, it cannot be inverted.
We can always include in the s-vector six additional coordinates to describe
the three translations and the three rotations. These six coordinates are
most conveniently chosen in the form
1 X√ X mI 1/2
R1 = √ mI qIx R4 = yI0 qIz − zI0 qIy
M I I
Ix
1 X√ X mI 1/2
R2 = √ mI qIy R5 = zI0 qIz − x0I qIy
M I I
Iy
1 X√ X mI 1/2
R3 = √ mI qIz R6 = x0I qIz − yI0 qIy
M I I
Iz
P
where M = I mI is the total mass of the molecule and Ix , Iy , and Iz are
the moments of inertia. x0I , yI0 , and zI0 are the coordinates of the equilibrium
position of the Ith atom relative to the centre of mass and
√ √ √
qIx = mI ∆xI qIy = mI ∆yI qIz = mI ∆zI
In this way B is a square and non-singular matrix (i.e. |B| 6= 0) and therefore
can be inverted. We can thus define the inverse transformation to get from
internal to Cartesian coordinates
x = B−1 s
where B−1 is the inverse matrix of B including translation and rotation.
3.4.1 Z-matrix
A special choice of internal coordinates often used as input to quantum chem-
ical programs is the Z-matrix. The Z-matrix of a molecule is build from a
list of the atoms. The coordinates of each new atom (A) are characterised
by
1. The distance to another, already defined atom B
2. The angle between the bond to atom B and the bond of atom B to an
already defined atom C.
3. The dihedral angle between A, B, C, and an already defined atom D.
50
The first three atoms are special cases with
• The first atom is at the origin of the Cartesian coordinate system.
• The second atom is defined by the distance to atom 1 and is oriented

along a Cartesian coordinate axis.
• The third atom is defined by the distance to atom 1 or 2 and the angle
(3-2-1) or the angle (3-1-2)
The Z-matrix of a molecule is not unique. There are many different possibil-
ities to order the atoms and to define distances, angles and dihedrals.
Example: Z-matrix for staggered conformation of ethane.
3
7
1
6 5
2
4
1 C
2 C 1.54 1
3 H 1.01 1 109.5 2
4 H 1.01 2 109.5 1 180.0 3
5 H 1.01 1 109.5 2 60.0 4
6 H 1.01 2 109.5 1 -60.0 5
7 H 1.01 1 109.5 2 180.0 6
8 H 1.01 2 109.5 1 60.0 7
In the first line of the Z-matrix atom 1, a carbon atom, is defined. Atom
number 2 is also a carbon that is a distance of 1.54 Å from atom 1 (columns
3 and 4). Atom 3 is a hydrogen atom that is bonded to atom 1 with a bond
length of 1.01 Å. The angle formed by atoms 2-1-3 is 109.5o , information that
is specified in columns 5 and 6. The fourth atom is a hydrogen, a distance
of 1.01 Å from atom 2, the angle 4-2-1 is 109.5o , and the torsion angle for
atoms 4-2-1-3 is 180o . Thus for all except the first three atoms, each atom has
51
three internal coordinates: the distance of the atom from one of the atoms
previously defined, the angle formed by the atom and two of the previous
atoms, and the torsion angle defined by the atom and three of the previous
atoms.
The Z-matrix can often be written in many different ways as there are many
combinations of internal coordinates. As a rule internal coordinates should
be chosen to have minimal coupling. However, it should be stressed that the
final results of a quantum mechanical calculations are not dependent on the
actual coordinate system chosen.
52
Lecture 4
Two-Electron Systems and Spin
4.1 Helium atom: Ground state

4.1.1 Basics
Electronic Hamilton operator within Born–Oppenheimer approximation:
~ = T̂e (~r) + VNe (~r; R)

Ĥel (~r; R) ~ + Vee (~r)
Coordinate system: nuclei at origin, spherical coordinates for electron 1

(r1 , θ1 , φ1 ) and electron 2 (r2 , θ2 , φ2 ).
Electronic Hamilton operator for He (Z=2):
~2 Ze2 Ze2 e2
Ĥel (~r1 , ~r2 ) = − (∇21 + ∇22 ) − − +
2me r1 r2 r12
Corresponding Schrödinger equation:
Ĥel (~r1 , ~r2 )Ψel (~r1 , ~r2 ) = Eel Ψel (~r1 , ~r2 )
For this differential equation with 6 variables there is no exact solution. We

are looking for approximations, first for the ground state.
Experiment:
The energy E0 of the ground state is equal to the sum of ionisation potentials
for He (24.59 eV) and He+ (54.42 eV).
E0 (exp) = −79.01 eV
This is our reference to assess different approximations.
53
4.1.2 Orbital model
Ansatz:
Ψel (~r1 , ~r2 ) = Ψ1 (~r1 )Ψ2 (~r2 )
One electron functions are called orbitals. In the orbital model the N-electron
function is approximated by a product of N orbitals.
Probability interpretation:
Ψ∗el Ψel dτ = [Ψ1 (~r1 )∗ Ψ1 (~r1 )dτ1 ] [Ψ2 (~r2 )∗ Ψ2 (~r2 )dτ2 ]
W12 = W1 W2
The probability W12 to find at the same time electron 1 at ~r1 and electron 2
at ~r2 is equal to the product of the single probabilities W1 and W2 .
In the orbital model the probability distribution of an electron is independent
from the other electrons (’independent particle model’).
Note:
More accurate methods that go beyond the orbital model have to include the
correlation of electrons.
4.1.3 Formal separation

The Hamilton operator consists of one- and two-electron parts
e2
Ĥel = Ĥ (0) + V̂ee = Ĥ (0) +
r12
2
X
Ĥ (0) = ĥi (~ri )
i=1
~2 2 Ze2
ĥi = − ∇ −
2me i ri
The one-electron operator ĥi is the same as the Hamilton operator of a

hydrogen-like atom with Z = 2 for He. Eigenvalues and eigenfunctions Ψi
are known.
ĥi Ψk (~ri ) = k Ψk (~ri )

Z 2 e2 1
k = − with ni = 1, 2, 3, . . .
2a n2k
Zr1
Ψk = π −1/2 (Z/a)3/2 e− a for nk = 1
Crude Approximation: With V̂ee = 0 the Hamilton operator separates like

Ĥ (0) . The two-electron wave function can in this case be written as a product
54
of orbitals (the solutions for the hydrogen atom with Z = 2), and we get for
the energy
E (0) = 1 (n1 ) + 2 (n2 )
Ground state of the He atom with n1 = n2 = 1:
(0) Z 2 e2
E0 =− = −108.84 eV
a
(0) Z3 Z
Ψ0 = 3 e− a (r1 +r2 )
πa
Because of the assumption V̂ee = 0 we get a huge discrepancy to the experi-
mental ground state energy (−79.01 eV).
4.1.4 Perturbation theory

We take the electron-electron interaction V̂ee as the perturbation operator
Ĥ (1) . This is not really a valid choice as V̂ee is not small compared to the other
terms in the Hamilton operator. Nevertheless we calculate the correction to
the energy to first order:
(1) (0) (0) (0) e2 (0)
E0 = hΨ0 | V̂ee | Ψ0 i = hΨ0 | | Ψ0 i
r12
(0)
Inserting the analytic solutions for Ψ0 leads to a difficult six-dimensional
integral for the variables (r1 , θ1 , φ1 , r2 , θ2 , φ2 ), that can be solved analytically.
(1) 5Z e2
E0 =
8 a
For He with Z = 2 we get:
(1)
E0 = 34.01 eV
(0) (1)
E0 + E0 = −74.83 eV
With the perturbation correction to first order the energy is off by 5.3% from
the experimental value (-79.01 eV).
Corrections to higher order are possible, but much more complicated than
simple variational calculations.
4.1.5 Variational calculations

Idea: In the solutions used up to know, both electrons were seeing the field
of the nuclei with Z = 2. This neglects the partial screening of the
nuclear charge due to the other electron.
55
Ansatz: We replace the nuclear charge Z = 2 in the wave function for the
ground state by an effective nuclear charge η < 2.
η 3 − η (r1 +r2 )
Ψ(r1 , r2 ) = e a
πa3
Variational calculation for the determination of the optimal value of η:
dEΨ d
= hΨ | Ĥel | Ψi = 0
dη dη
The calculation of the expectation value of Ĥel w.r.t. to Ψ affords the solution
of several integrals (difficult but analytic solutions exist).
Solution:

2 5 e2
hΨ | Ĥel | Ψi = η − 2ηZ + η
8 a
2
dEΨ 5 e
= 2η − 2Z + =0
dη 8 a
5
ηopt = Z −
16
2
2 2
2 e 5 e
Eopt = −ηopt = − Z −
a 16 a
For He with Z = 2 we get:
ηopt (He) = 1.6875

Eopt (He) = −77.49 eV
The nuclear charge is screened by 0.3e. The error in the ground state energy
is reduced from 5.3% with η = 2 and a perturbative calculation to 1.9% with
η = ηopt .
4.1.6 Overview of ground state energies

Goal: Assess accuracy of approximative solutions by comparison to experi-
ment.
56
E0 (a.u.) error %
Experiment -2.903784(3)
Orbital model
(0)
Zeroth order (without V̂ee ), E0 -4.00 -37.8
(0) (1)
Perturbation theory 1. order, E0 + E0 -2.75 5.3
Simple variational calculation, Eopt -2.8477 1.9
Best variational calculation -2.8617 1.4
Including electron correlation
Simple variational calculation, 40 Parameters -2.90276 0.04
Best variational calculation, 1078 Parameters -2.9037462 0.001
With additional relativistic corrections -2.9037848
4.2 Electron spin

Experimental fact:
Electrons have a proper angular momentum (Spin), the components of which
can have only two distinct values in any arbitrary direction:
~
sz = ±
2
Quantum mechanical description with spin operators:
ŝ2 = ŝ2x + ŝ2y + ŝ2z = ŝ+ ŝ− − ŝz + ŝ2z = ŝ− ŝ+ + ŝz + ŝ2z
ŝ+ = ŝx + iŝy
ŝ− = ŝx − iŝy
The linear combinations ŝ+ and ŝ− are useful for derivations. Spin operators
have the same commutation relations as the angular momentum operators:
[ŝx , ŝy ] = iŝz ; [ŝy , ŝz ] = iŝx ; [ŝz , ŝx ] = iŝy

[ŝ2 , ŝx ] = [ŝ2 , ŝy ] = [ŝ2 , ŝz ] = 0
Eigenfunctions:
For a one electron system there are two orthonormal eigenfunctions |αi and
|βi. Application of operators:
57
|αi |βi
3 3
ŝ2 4
|αi 4
|βi
1
ŝz 2
|αi − 12 |βi
1 1
ŝx 2
|βi 2
|αi
i
ŝy 2
|βi − 2i |αi
ŝ+ 0 |αi
ŝ− |βi 0
General eigenvalue equation for spin:
ŝz |σi = ms ~|σi

ŝ2 |σi = s(s + 1)~2 |σi
Spin eigenfunction |σi, quantum numbers s and ms for the total spin and its
z-component with −s ≤ ms ≤ s.
For a single electron we have:
1 1
ms = ± ; s=
2 2
Spin orbitals: For the complete description of an electron we need three
spatial variables and one spin variable. The spin orbitals Ψ(~r)|σi are products
of a spatial orbital Ψ(~r) and a spin eigenfunction |σi, |αi oder |βi.
Two electron systems:

Operators and eigenvalue equation:
Ŝx = ŝx1 + ŝx2 ; Ŝy = ŝy1 + ŝy2 ; Ŝz = ŝz1 + ŝz2

Ŝ 2 = Ŝx2 + Ŝy2 + Ŝz2
Because we have the general commutator relations for angular moment op-
erators, the eigenvalue equations have the form:
Ŝz |σ1 σ2 i = MS ~|σ1 σ2 i (−S ≤ MS ≤ S)

2 2
Ŝ |σ1 σ2 i = S(S + 1)~ |σ1 σ2 i
58
Spin eigenfunction |σ1 σ2 i, quantum numbers S and MS for total spin and its
z-component.
Relation to the Hamilton operator:

We use a non-relativistic Hamilton operator without any spin components,
therefore Hamilton operator and spin operators commute.
[Ĥ, Ŝ 2 ] = [Ĥ, Ŝz ] = 0
Conclusions:
• There are exact expectation values for Ŝ 2 and Ŝz . The system can be
characterised by the spin quantum numbers S and MS .
• For |σ1 σ2 i we can use a product function.
4.2.1 Electrons are indistinguishable

In microscopic systems similar particles (e.g. electrons) can experimentally
not be distinguished. The probability density Ψ∗ Ψ is in such many particle
systems not allowed to be dependent on the arbitrary enumeration of the
particles. A change of the numeration, exchange of particles, is not allowed
to change Ψ∗ Ψ.
Conclusion: With an exchange of two particles of the same type, the many-
body wavefunction Ψ either stays the same or changes sign. Wave functions
are either symmetric or anti-symmetric w.r.t. the exchange of two particles
of the same kind.
Spin eigenfunctions for two-electron systems:

A naive product function leads to the following combination


 α(1)α(2)

α(1)β(2)
|σ1 σ2 i =

 β(1)α(2)

β(1)β(2)
In the second and third product we can distinguish the two electrons, this is
not possible for the products one and four. In the first case symmetric wave
functions can be constructed through linear combinations:
2−1/2 [α(1)β(2) + β(1)α(2)]

2−1/2 [α(1)β(2) − β(1)α(2)]
59
An acceptable set of product functions is
α(1)α(2);
2−1/2 [α(1)β(2) + β(1)α(2)];
2−1/2 [α(1)β(2) − β(1)α(2)];
β(1)β(2)
Electron spin for two-electron systems: Singlets and triplets
Spin eigenfunction Exchange MS S Notation
α(1)α(2) symmetric 1 1 triplet
2−1/2 [α(1)β(2) + β(1)α(2)] symmetric 0 1 triplet
β(1)β(2) symmetric -1 1 triplet
2−1/2 [α(1)β(2) − β(1)α(2)] antisymmetric 0 0 singlet
Simple example:
Ŝz α(1)α(2) = (ŝz1 + ŝz2 )α(1)α(2)

= ŝz1 α(1)α(2) + ŝz2 α(1)α(2)
1 1
= ~α(1)α(2) + ~α(1)α(2)
2 2
= ~α(1)α(2) → MS = 1
4.3 Pauli principle

All in nature realisable states of many-electron systems are described by
space-spin-wave functions that are anti-symmetric w.r.t. electron exchange.
Special case of the Pauli principle: Every spatial orbital can be occupied by
at most two electrons with opposite spin.
Generalisation to other particles:
Fermions = Particles with half-integer spin

Bosons = Particles with integer spin
60
All systems with several equal fermions (bosons) are described by wave func-
tions, that are anti-symmetric (symmetric) w.r.t. exchange of two equal
fermions (bosons).
Example for fermions: electron, proton, neutron

Example for bosons: H-Atom, He-nuclei
61
62
Lecture 5
Hartree–Fock Approximation
5.1 Slater-determinant
We are looking for a n-electron wave function Ψ within the orbital model that
satisfies the Pauli principle. This wave function is build from orthonormal
spin orbitals η.
η(~r, σ) = Ψ(~r)|σi
hηi |ηj i = δij
Short hand for electron 1:

Ψα(1) Ψ(1)
η(1) = Ψ(1)|σ(1)i = =
Ψβ(1) Ψ̄(1)
Ansatz

η1 (1) η1 (2) · · · η1 (N )

1 η2 (1) η2 (2) · · · η2 (N )

Ψ(1, 2, . . . , N ) = √ .. .. ..
N ! . . .

ηN (1) ηN (2) · · · ηN (N )
Exchange of two electrons corresponds to an exchange of two columns. This

leads to a change of sign for the determinant. The determinantal wave func-
tions is therefore anti-symmetric w.r.t. an exchange of electrons (Pauli prin-
ciple).
A determinant vanishes if two rows are the same. In non-trivial wave func-
tions Ψ 6= 0 all spin orbitals have to be different. A given spatial orbital can
only appear twice in the determinant, combined with either α or β.
Pauli principle in the orbital model: A spatial orbital can at most be occupied
by two electrons with opposite spin.
63
Writing out the determinant leads to N ! products of the form η1 η2 · · · ηN . All
possible permutations of the N electrons in the N spin orbitals appear. In
this approach the electrons are evenly spread out over all spin orbitals and
therefore indistinguishable.
Short notation: Often we omit the normalisation constant (N !)−1/2 and we
only write the diagonal elements.
Ψ(1, 2, . . . , N ) = |η1 (1)η2 (2) · · · ηN (N )|
Slater determinants are always eigenfunctions of Ŝz , however, they are not
necessarily eigenfunctions of Ŝ 2 . For the general case there are always linear
combination of determinants that are eigenfunctions of Ŝz and Ŝ 2 at the same
time. Such spin-adapted linear combination of determinants (configurations)
are needed to describe open-shell systems.
Important special case: In closed-shell systems the spins of all electrons are
paired.

Ψ(1, 2, . . . , N ) = Ψ1 (1)Ψ̄1 (2), Ψ2 (3)Ψ̄2 (4), . . . ΨN/2 (N − 1)Ψ̄N/2 (N )
This Slater determinant is an eigenfunction of Ŝz and Ŝ 2 with S = MS = 0

(Singlet).
5.1.1 Energy expectation value for Slater determinants
The electronic energy Eel is the expectation value of the non-relativistic elec-
tronic Hamilton operators Ĥel w.r.t. the N -electron wave function, written
as a Slater determinant of orthonormal spin orbitals.
Eel = hΨ|Ĥel |Ψi

N N
X X e2
Ĥel = ĥ(i) +
i=1
r
i<j ij
~2 2 X Z A e 2
ĥ(i) = − ∇ −
2me i A
riA
64
Using the Slater determinant function we get
N
X N
X
(η) (η) (η)
Eel = hii + Jij − Kij
i=1 i<j
Z
(η)
hii = ηi∗ (1)ĥ(1)ηi (1)dτ1
Z Z
(η) e2
Jij = ηi∗ (1)ηi (1) ηj∗ (2)ηj (2)dτ1 dτ2
r12
Z Z
(η) e2
Kij = ηi∗ (1)ηj (1) ηj∗ (2)ηi (2)dτ1 dτ2
r12
(η) (η) (η)

hii is the one-electron integral, Jij Coulomb integral, and Kij an exchange
integral. In these equations all integrals are over spin orbitals η. By inte-
gration over the spin coordinates we get to a description in terms of spatial
orbitals.
Example for the transition to spatial orbitals:
ηi = Ψi |αi
(η) (Ψ) (Ψ)
hii = hii hα|αi = hii
Z
(Ψ)
hii = Ψ∗ (~r1 )ĥ(1)Ψi (~r1 )d~r1
A similar separation of spatial and spin parts occurs in all integrals. This is
due to the fact that Ĥel does not contain any spin coordinates. A systematic
investigation of all integrals leads to:
(η) (Ψ)
hii = hii
(η) (Ψ)
Jij = Jij
(Ψ)
(η) Kij equal
Kij = if spin functions in ηi , ηj
0 opposite
Interaction of two electrons in spatial orbitals Ψi , Ψj :
Jij with opposite spin
Jij − Kij with equal spin (smaller, because Kij > 0)
Important special case: closed shell systems

N/2 N/2 N/2
(Ψ)
X (Ψ)
X X (Ψ) (Ψ)
Eel = 2 hii + 2Jij − Kij
i=1 i=1 j=1
65
5.1.2 Unitary transformation of spin orbitals
In the definition of the Slater determinant we regard the spin orbitals ηi (j)
formally as elements Aji of a matrix A.
Ψ = (N !)−1/2 det(A)
A unitary transformation U of the spin orbitals leaves them orthonormal:

0
X
ηk = ηi Uik
i
−1 †
U =U , (U −1 )ik = (U † )ik = (U ∗ )ki
In matrix form:
A0 = AU
Ψ0 = (N !)−1/2 det(A0 )
= (N !)−1/2 det(AU )
= (N !)−1/2 det(A) det(U )
Ψ0 = Ψ det(U )
Calculation of det(U ) for a unitary transformation:
U †U = I
det(U † U ) = det(U † ) det(U ) = [det(U )]∗ det(U ) = |det(U )|2 = 1
det(U ) = eiΦ
Ψ0 differs from Ψ at most by a phase factor eiΦ with ±1 for the real case.
For all expectation values we get
hΨ|Ô|Ψi = hΨ0 |Ô|Ψ0 i
In this sense determinantal wave functions are invariant w.r.t. unitary trans-
formation of spin orbitals.
5.2 Hartree–Fock method

Goal: Calculation of optimal orbitals
Approach: Variational method, minimisation of the energy expectation

value for the Slater determinant under the condition that orbitals stay
orthonormal.
66
Derivation for spin and spatial orbitals is similar. We study the more general
case of spin orbitals. For orthonormal spin orbitals we have the energy in
simplified notation:
N N N
X 1 XX
Eel = hii + (Jij − Kij )
i=1
2 i=1 j=1
Sij = hηi |ηj i = δij , Sij = Sij∗
All integrals are for spin orbitals. Spin orbitals may be complex. Optimal
spin orbitals minimise the energy expression while maintaining the constraint
condition. This can be achieved by the method of Lagrange multipliers:
N X
X N
L[{ηi }] = Eel [{ηi }] − ji (Sij − δij )
i=1 j=1
ji = ?ij
δL = 0
Lagrange multipliers ji have to be Hermitian in order for the Lagrange-

Function L to be real. Variation of the spin orbitals must lead to a vanishing
variation of L:
N X
X N
δL = δEel − ji δSij = 0
i=1 j=1
We investigate the variation of the integrals that is caused by the variation

of the spin orbitals (c.c. = conjugate complex)
Z Z
δSij = δηi∗ (1)ηj (1)dτ1 + ηi∗ (1)δηj (1)dτ1
N N
X 1 XX
δEel = δhii + (δJij − δKij )
i
2 i=1 j=1
Z Z
δhii = δηi (1)ĥ(1)ηj (1)dτ1 + ηi∗ (1)ĥ(1)δηj (1)dτ1
∗
Z
δhii = δηi∗ (1)ĥ(1)ηj (1)dτ1 + c.c.
67
Because of Hermiticity of ĥ, the second term in the equation for δhii is the
complex conjugate of the first term.
Z Z
e2 ∗
δJij = dτ1 dτ2 δηi∗ (1)ηi (1) η (2)ηj (2) A
r12 j
Z Z
∗ e2 ∗
+ dτ1 dτ2 ηi (1)δηi (1) ηj (2)ηj (2) B=A∗
r12
Z Z
∗ e2 ∗
+ dτ1 dτ2 ηi (1)ηi (1) δηj (2)ηj (2) C
r12
Z Z
e2
+ dτ1 dτ2 ηi∗ (1)ηi (1) ηj∗ (2)δηj (2) D=C∗
r12
N X
N N X
N Z Z
X X e2
δJij = 2 dτ1 dτ2 δηi∗ (1)ηi (1) ηj∗ (2)ηj (2) + c.c.
i=1 j=1 i=1 j=1
r12
In the double sum in the last equation the terms A and C appear in equivalent
pairs and can be combined (same is true for δKij ):
N X
N N X N Z Z
X X e2
δKij = 2 dτ1 dτ2 δηi∗ (1)ηj (1) ηj∗ (2)ηi (2) + c.c.
i=1 j=1 i=1 j=1
r12
For the remaining terms we get:
N X
X N N X
X N Z
ji δSij = ji δηi∗ (1)ηj (1)dτ1 + c.c.
i=1 j=1 i=1 j=1
A further simplification can be achieved with the introduction of the Coulomb

operator Jˆj and exchange operator K̂j :
Z
e2
Jˆj (1)ηi (1) = dτ2 ηj∗ (2)
ηj (2)ηi (1)
r12
Z
e2
K̂j (1)ηi (1) = dτ2 ηj∗ (2) ηi (2)ηj (1)
r12
Z
Jij = dτ1 ηi∗ (1)Jˆj (1)ηi (1)
Z
Kij = dτ1 ηi∗ (1)K̂j (1)ηi (1)
68
Taking all terms together we get:
XZ h
dτ1 δηi∗ (1) ĥ(1)ηi (1)
i
#
X X
+ Jˆj (1) − K̂j (1) ηi (1) − ji ηj (1) + c.c. = 0
j j
Variations δηi∗are arbitrary and therefore every single integral to an index i

has to disappear and the integrand [. . .] has to be zero.
" #
X X
ĥ(1) + Jˆj (1) − K̂j (1) ηi (1) = ji ηj (1) for i = 1, . . . , N
j j
The same steps performed on the conjugate-complex part lead to the same
result.
Introduction of the Fock-operator:

N
X
F̂ (1) = ĥ(1) + Jˆj (1) − K̂j (1)
j=1
X
F̂ (1)ηi (1) = ji ηj (1)
j
The current derivation did not lead to the desired eigenvalue equation for spin
orbitals. We now investigate unitary transformations of the spin orbitals.
We have already shown that determinantal wave functions and their expec-
tation values are invariant w.r.t. such unitary transformations. The same is
true for the Fock operator. For the Lagrange multipliers we get ji :
X
hηk (1)|F̂ (1)|ηi (1)i = ji hηk (1)|ηj (1)i = ki
j
The multipliers are equivalent to the matrix elements of the Fock operator.
A transformation of the spin orbitals leads to (in matrix form):
0 = U † U
is an Hermitian matrix. Therefore there is a special unitary matrix U ,
which transforms to diagonal form. For this special set of transformed spin
orbitals we get
0 0 0
F̂ (1)ηi (1) = i ηi (1)
Normally one omits the special mark for these canonic spin orbitals and one
writes:
F̂ ηi = i ηi for i = 1, . . . , N
69
5.2.1 Solution of the Hartree-Fock equations
Problem: The Fock operator consists of Coulomb Jˆj and exchange operators
K̂j . Therefore, it depends on the spin orbitals ηi , we want to calculate
with the Hartree-Fock equations.
Solution: Iterative method starting from an arbitrary but reasonable initial

(0)
set of orbitals ηi .
Hartree-Fock equations have to be solved iteratively.

Scheme
(0) (1) (2)
{ηi } −→ F̂ (1) −→ {ηi } −→ F̂ (2) −→ {ηi } −→ etc.
This sequence is stopped, if the solutions in two consecutive iterations don’t

change (within a given tolerance). This method is called SCF (Self-Consistent-
Field).
For a known (self-consistent) set of occupied spin orbitals the Fock operator F̂
is a well-defined Hermitian operator with an infinite number of eigenfunctions
ηi (N occupied spin orbitals with lowest orbital energies and un-occupied
(virtual) orbitals with higher orbital energies).
5.2.2 Orbital energies and Koopmans’ theorem

For an N -electron system, minimisation of the energy of the determinant
|Ψi = |η1 η2 · · · ηN i leads to an eigenvalue equation F |ηa i = a |ηa i for the n
occupied spin orbitals {ηA }. The Fock operator has a functional dependence
on these occupied spin orbitals, but once the occupied spin orbitals are known
the Fock operator becomes a well-defined Hermitian operator, which will have
an infinite number of eigenfunctions
F |ηj i = j |ηj i j = 1, 2, . . . , ∞
Each of the solutions has a spin orbital energy j . The N spin orbitals with
the lowest orbital energies are just the spin orbitals occupied in |Ψi for which
we use the indices a, b, . . . . the remaining infinite number of spin orbitals
with higher energies are the virtual or un-occupied spin orbitals, which we
label with the indices r, s, . . . .
The matrix elements of the Fock operator in the basis of the spin orbital
eigenfunctions is diagonal with diagonal elements equal to the orbital energies
hηi | F | ηj i = j hηi | ηj i = j δij
70
The orbital energies can be expressed
X
i = hηi | F | ηi i = hηi | h | ηi i + hηi | (Jb − Kb ) | ηi i
b
X X
= hηi | h | ηi i + hηi | Jb | ηi i − hηi | Kb | ηi i
b b
X
= hi | h | ii + ((ii, bb) − (ib, bi)) ,
b
where we have introduced the short form (ij, kl) for the two-electron integrals
Z Z
e2
(ij, kl) = dτ1 dτ2 δηi∗ (1)ηj (1) ηk∗ (2)ηl (2) .
r12
In particular we get for occupied orbitals
X
a = ha | h | ai + ((aa, bb) − (ab, ba))
b
X
= ha | h | ai + ((aa, bb) − (ab, ba))
b6=a
and for virtual orbitals

X
r = hr | h | ri + ((rr, bb) − (rb, br))
b
The orbital energy a represents the energy of an electron in the spin orbital
|ηa i. This energy is the kinetic energy and the attraction to the nuclei (ha |
h | ai) plus a Coulomb and exchange interaction with each of the remaining
N − 1 electrons. The result for r has a different character. It includes the
kinetic energy and the nuclear attraction of an electron in |ηr i, but includes
Coulomb and exchange interactions with n electrons of the Hartree–Fock
ground state. It is as if an electron has been added to |Ψ0 i to produce a
(N + 1)-electron state.
If we simply add up the orbital energies of the occupied states, we get
N
X N
X N X
X N
a = ha | h | ai + ((aa, bb) − (ab, ba))
a a a b
If we compare this with the total Hartree–Fock energy

N N N
X 1 XX
EHF = ha | h | ai + ((aa, bb) − (ab, ba))
a
2 a b
71
we see that
N
X
EHF 6= a
a
and the total energy of the state |Ψ0 i is not just the sum of the orbital
energies. The reason is as follows. The energy a includes Coulomb and
exchange interactions between an electron in orbital a and electrons in all
other occupied spin orbitals, in particular orbital b. But b includes Coulomb
and exchange interactions between an electron in orbital b and electrons in
all other occupied spin orbitals, in particular orbital a. Thus if we add a
and b we include the electron-electron interactions between an electron in a
and one in b twice. The sum of orbital energies counts the electron-electron
interactions twice. This is the reason for the factor 21 in the correct energy
expression for the total energy relative to the sum of orbital energies.
If the total energy is not the sum of the orbital energies what physical sig-
nificance can we attach to orbital energies?
Koopmans’ theorem
Given an N -electron Hartree–Fock single determinant |N Ψ0 i with occupied
and virtual spin orbital energies a and r , then the ionisation potential to
produce an (N − 1)-electron single determinant |N −1 Ψa i with identical spin
orbitals, obtained by removing an electron from spin orbital ηa , and the
electron affinity to produce an (N + 1)-electron single determinant |N +1 Ψr i
with identical spin orbitals, obtained by adding an electron to spin orbital
ηr , are just −a and −r , respectively.
Proof: The ionisation potential is
IP =N −1 Ea −N E0
where N −1 Ea and N E0 are the expectation values of the energy of the two
relevant single determinants
N
E 0 = h N Ψ 0 | H |N Ψ 0 i
N −1
Ea = hN −1 Ψa | H |N −1 Ψa i
72
The ionisation potential for orbital c is therefore
IP =N −1 Ec −N E0
X 1 XX
= ha | h | ai + ((aa, bb) − (ab, ba))
a6=c
2 a6=c b6=c
X 1 XX
− ha | h | ai + ((aa, bb) − (ab, ba))
a
2 a b
1X 1X
= −hc | h | ci − ((aa, cc) − (ac, ca)) − ((cc, bb) − (bc, cb))
2 a 2 b
X
= −hc | h | ci − ((cc, bb) − (bc, cb))
b
= −c
Occupied orbital energies are generally negative and ionisation potentials are
positive.
With a very similar calculation we can also calculate
EA =N E0 −N +1 E r = −r
This result is consistent with the previous observation that r included inter-
actions with all N other electrons of the ground state and thus describes an
N + 1th electron.
Koopmans’ ionisation potentials are reasonable first approximations to ex-
perimental ionisation potentials (missing are orbital relaxation and corre-
lation effects). Koopmans’ electron affinities are unfortunately often bad.
Many neutral molecules will add an electron to form a stable negative ion.
Hartree–Fock calculations on neutral molecules, however, almost always give
positive orbital energies for all the virtual orbitals.
5.2.3 Brillouin’s theorem

The Hartree–Fock equation produces a set of spin orbitals. The single de-
terminant |Ψ0 i, formed from the N lowest spin orbitals is the Hartree–Fock
approximation to the ground state. There are many other determinants with
N functions that can be formed from the full set of spin orbitals. Brillouin’s
theorem makes a statement of a subset of these determinants. This subset
is the set of singly excited determinants |Ψra i obtained from |Ψ0 i by a single
replacement of ηa with ηr .
Brillouin’s theorem
Singly excited determinants |Ψra i will not interact directly with a reference
73
Hartree–Fock determinant |Ψ0 i, i.e., hΨ0 | H | Ψra i = 0 .
Calculating the matrix element

X
hΨ0 | H | Ψra i = ha | h | ri + ((ar, bb) − (ab, br))
b
The right hand side of this equation is just the matrix element of the Fock
operator in the spin orbital basis. Therefore,
hΨ0 | H | Ψra i = hηa | F | ηr i
The matrix element that mixes singly excited determinants with |Ψ0 i is thus
equal to an off-diagonal element of the Fock matrix. Now, by definition, solv-
ing the Hartree–Fock eigenvalue problem requires the off-diagonal elements
to satisfy hηi | F | ηj i = 0 (i 6= j). One can then say that solving the
Hartree–Fock eigenvalue equation is equivalent to ensuring that |Ψ0 i will not
mix with any singly excited determinants.
74
Lecture 6
Molecular Orbital Theory
6.1 LCAO-Ansatz
LCAO: Linear Combination of Atomic Orbitals
In the LCAO approach the molecular orbitals (MOs) Ψi are written as linear
combinations of atomic orbitals (AOs) Φµ .
X
Ψi = cµi Φµ
µ
In the following the equations can also be applied to other types of basis sets
that are not necessarily of atomic type.
Conventions:
i, j, k, l, . . . for MOs
µ, ν, λ, σ, . . . for AOs
We restrict the following derivations to closed-shell systems, singlet states
with S = MS = 0, with all orbitals doubly occupied.
Energy formula:
N/2 N/2 N/2
X X X
Eel = 2 hii + (2Jij − Kij )
i=1 i=1 j=1
Z
hii = Ψ∗i (1)ĥ(1)Ψi (1)dτ1 = hΨi |ĥΨi i
Z Z
e2
Jij = Ψ∗i (1)Ψi (1) Ψ∗j (2)Ψj (2)dτ1 dτ2 = (ii, jj)
r12
Kij = (ij, ij)
We insert the LCAO equation and express all MO integrals with AO integrals.
We assume that MOs and AOs are real functions.
75
Connection of MO and AO integrals:
* +
X X XX
hii = cµi Φµ | ĥ cνi Φν = cµi cνi hµν
µ ν µ ν
hµν = hΦµ |ĥΦν i

XXXX
Jij = cµi cνi cλj cσj (µν, λσ)
µ ν λ σ
XXXX
Kij = cµi cνj cλi cσj (µν, λσ)
µ ν λ σ
Z Z
e2
(µν, λσ) = Φµ (1)Φν (1) Φλ (2)Φσ (2)dτ1 dτ2
r12
Energy expression in AO basis:
XXX
Eel = 2 cµi cνi hµν
i µ ν
XXXXXX
+ cµi cνi cλj cσj [2(µν, λσ) − (µλ, νσ)]
i j µ ν λ σ
Introduction of the density matrix Pµν :

N/2
X
Pµν = 2 cµi cνi
i=1
Simplification of the energy expression in the AO basis:

XX 1 XXXX
Eel = Pµν hµν + Pµν Pλσ [2(µν, λσ) − (µλ, νσ)]
µ ν
4 µ ν λ σ
For a given basis {Φµ } all the AO integrals hµν and (µν, λσ) are (in principle)
known.
The LCAO coefficients cµi and the density matrix elements Pµν are not
known. They can be calculated from a linear variational calculation (Root-
haan-Hall method).
6.2 Roothaan-Hall equations

Instead of a complete derivation it is easier to insert the LCAO approach
in the general Hartree–Fock equations and multiply from left with Φ∗µ and
76
integrate:
! !
X X
F̂ cνi Φν = i cνi Φν
ν ν
X X
cνi hΦµ |F̂ Φν i = i cνi hΦµ |Φν i
ν ν
Introduction of the Fock matrix F and the overlap matrix S in AO basis:
Fµν = hΦµ |F̂ Φν i

Sµν = hΦµ |Φν i
Eigenvalue equation:
X X
Fµν cνi = i Sµν cνi i = 1, 2, . . . , K
ν ν
F C = SC Matrix-Notation
F = {Fµν }, C = {cνi }, S = {Sµν }, = {i δij }
For K basis functions Φµ there are K equations. The matrices have the
dimension K × K; F and S are symmetric, is diagonal.
Fock matrix elements in AO basis:
N/2
X
Fµν = hΦµ |ĥ + (2Jˆj − K̂j )|Φν i
i=1
N/2
X
Fµν = hµν + [2(µν, jj) − (µj, jν)]
j=1
N/2
XXX
Fµν = hµν + cλj cσj [2(µν, λσ) − (µλ, σν)]
j=1 λ σ
XX 1
Fµν = hµν + Pλσ [(µν, λσ) − (µλ, σν)]
λ σ
2
Solution of the Roothaan equations

As for the general Hartree–Fock equations, the Roothaan equations can only
be solved iteratively: The Fock matrix elements Fµν depend on the density
matrix elements Pλσ . The density matrix has to be calculated from the
unknown LCAO coefficients cνi .
A LCAO-MO-SCF calculation:
1. Choose a basis of atomic orbitals {Φµ }
77
2. Calculate the AO integrals hµν , (µν, λσ) and Sµν .
(0)
3. Choose a reasonable initial density matrix {Pµν }.
4. Calculate the corresponding Fock matrix.
5. Solve the eigenvalue equation (Roothaan–Hall).

6. Calculate a new density matrix.
7. Iterate (4)-(5)-(6) until self consistency is achieved.
Solving the eigenvalue equations:
We examine a linear transformation of the non-orthogonal basis {Φµ } into
0
an orthogonal basis {Φµ }.
0
X
Φµ = Xνµ Φν
ν
0 0
X X
δµν = hΦµ |Φν i = h Xλµ Φλ | Xσν Φσ i
λ σ
XX
∗
= Xλµ Sλσ Xσν
λ σ
†
I = X SX
The overlap matrix S in the basis {Φµ } is Hermitian and can be diagonalised
using a unitary transformation U .
U † SU = s
The diagonal matrix s contains the eigenvalues of S (sii > 0). For the
transformation matrix X the following approaches are used:
X = S −1/2 = U s−1/2 U † symmetric orthogonalisation

−1/2
X = Us canonical orthogonalisation
The symmetric orthogonalisation is used most often. The canonical orthog-

onalisation has the advantage that linear dependencies in the basis (sii ≈ 0)
can easily be eliminated.
Connection between the different sets of LCAO coefficients:
" #
X 0 0 X 0 X X X 0
X
Ψi = cµi Φµ = cµi Xνµ Φν = Xνµ cµi Φν = cνi Φν
µ µ ν ν µ ν
0 0 −1
C = XC , C =X C
78
With these definitions we can rewrite the eigenvalue equations:
0 0
F XC = SXC
0 0 0
(X † F X)C = (X † SX)C = C
0 0
Using the Fock matrix F in the basis {Φµ } we get a standard eigenvalue
equation:
0
F = X †F X
0 0 0
F C =C
Modified LCAO-MO-SCF procedure:
1. Choose a basis of atomic orbitals {Φµ }
2. Calculate AO integrals hµν , (µν, λσ) and Sµν . In addition calculate the
matrix X.
(0)
3. Choose a reasonable initial density matrix {Pµν }.
4. Calculate the Fock matrix.
5. Solve the Roothaan–Hall equations.
(a) Transform the Fock matrix

(b) Solve eigenvalue equation
(c) Back transform coefficients
6. Calculate a new density matrix
7. Iterate (4)-(5)-(6) until self-consistency is reached
6.3 Unrestricted open-shell Hartree–Fock

Up to now we have first derived the Hartree–Fock equation for general spin
orbitals η
F η i = i ηi
and then have had a look at the special case of closed shell orbitals. The
special case lead to the restricted Hartree–Fock (RHF) equations. This spe-
cial case has the advantage that the wavefunction is an eigenfunction not
only of the Hartree–Fock Hamiltonian, but also of the spin operators Sz and
79
S 2 . For the case of unpaired spins there are two possible methods avail-
able. Either the electrons are partitioned in paired and non-paired sets and
the paired electrons treated with equal spatial orbitals (the restricted open-
shell Hartree–Fock approach, ROHF) or all electrons are treated on the same
footing using general spin orbitals in an unrestricted Hartree–Fock approach
(UHF). The first method has the advantage that the wavefunction is still
an eigenfunction of the spin operators but the disadvantage that it needs a
rather involved theory and that the solution is not variationally stable w.r.t.
an unrestricted wavefunction. The UHF method is conceptually simpler and
leads to the variational ground state, but might lead to wavefunctions that
are no longer eigenfunctions of S 2 .
6.3.1 Spin contamination

We start from a Slater determinant
|ΨUHF i = |η1 η2 · · · ηN i
where the spin orbitals are
(
φα (r)α(ω)
η(x) =
φβ (r)β(ω)
either occupied with an α or β electron. There are Nα α and Nβ β electrons
in the system.
N = N α + Nβ
If we assume that Nα > Nβ then we can easily find that the expectation
value of Sz is
1
hSz i = (Nα − Nβ )
2
The exact expectation value of a corresponding wavefunction for S 2 would
be
2 1
hSexact i = (Nα − Nβ ) (Nα − Nβ + 2)
4
However, if we calculate the expectation value for an unrestricted Slater
determinant, we get
Nβ
Nα X
X
2
hSUHF i = 2
hSexact i + Nβ − |Sijαβ |2
i j
where
Sijαβ = hφαi | φβj i
The value for hSUHF i should always be monitored when working with unre-
stricted wavefunctions.
80
6.3.2 Pople–Nesbet equations
The Hartree–Fock equations for a spin orbital is
F φσi σ(ω1 ) = σi φσi σ(ω1 )
where σ stands for a general spin function. If we multiply from the left with
σ(ω1 ) and integrate over the spin variable ω1 we get
F α φαi = αi φαi

F β φβi = βi φβi
where the spatial Fock operators are defined by

Z
F (r1 ) = dω1 α? (ω1 )F (r1 , ω1 )α(ω1 )
α
Z
F (r1 ) = dω1 β ? (ω1 )F (r1 , ω1 )β(ω1 )
β
The operator F α (1) is the kinetic energy, nuclear attraction, and effective
potential of an electron of α spin. The effective interactions of an electron of
α spin include a Coulomb and exchange interaction with all other electrons
of α spin plus only a Coulomb interaction with electrons of β spin. Thus
Nα Nβ
X X
α
F (1) = h(1) + [Jaα (1) − Kaα (1)] + Jaβ (1)
a a
The sum over the Nα orbitals formally includes the interaction of an α elec-
tron with itself. However, since
[Jaα (1) − Kaα (1)] φαa (1) = 0
this self-interaction is eliminated. The corresponding Fock operator for elec-

trons of β spin is
Nβ Nα
X X
β
F (1) = h(1) + Jaβ (1) − Kaβ (1) + Jaα (1)
a a
The kinetic energy and nuclear attraction of an electron in one of the unre-
stricted orbitals φαi or φβi is the expectation value
hαii = hφαi | h | φαi i or hβii = hφβi | h | φβi i
81
The Coulomb interaction of an electron in φαi with one of opposite spin in
φβj is
Jijαβ = Jjiβα = hφαi | Jjβ | φαi i = hφβj | Jiα | φβj i = (φαi φαi , φβj φβj )
The corresponding Coulomb interactions between electrons of the same spin

are
Jijαα = hφαi | Jjα | φαi i = hφαj | Jiα | φαj i = (φαi φαi , φαj φαj )
and
Jijββ = hφβi | Jjβ | φβi i = hφβj | Jiβ | φβj i = (φβi φβi , φβj φβj )
The exchange interactions between electrons of parallel spin are
Kijαα = hφαi | Kjα | φαi i = hφαj | Kiα | φαj i = (φαi φαj , φαj φαi )
and
Kijββ = hφβi | Kjβ | φβi i = hφβj | Kiβ | φβj i = (φβi φβj , φβj φβi )
There is no exchange interaction between electrons of opposite spin.
The total energy for a unrestricted determinantal wavefunction is
Nα Nβ Nα XαN
X X 1X
E0 = hαaa + hβaa + αα
(Jab αα
− Kab )
a a
2 a b
Nβ Nβ α X N Nβ
1 X X ββ ββ
X
αβ
+ Jab − Kab + Jab
2 a b a b
The factor of 1/2 in the third and fourth terms eliminates the double counting
αα αα
in the free summation. The self-interaction disappears since Jaa − Kaa =
ββ ββ
Jaa − Kaa = 0.
Basis set expansion of the orbitals (LCAO)
K
X
ψiα = α
Cµi χµ i = 1, 2, . . . , K
µ=1
K
X
ψiβ = β
Cµi χµ i = 1, 2, . . . , K
µ=1
Notice that both sets of orbitals are expanded in only one set of basis function.
The two eigenvalue equations for the orbital sets guarantee that the eigen-
functions for α and β sets individually form orthogonal sets. There is no
reason, however, that a member of one set need be orthogonal to a member
82
of the other set. Even though the two sets of spatial orbitals overlap with
each other, the set of 2K spin orbitals will form an orthonormal set, either
from spatial orthogonality (αα and ββ case) or spin orthogonality (αβ case).
Substituting the expansions of the orbitals into the Hartree–Fock equations
gives X X
σ
Cνj F σ (1)χν (1) = σj σ
Cνj χν (1)
ν ν
If we multiply this equation by χ?µ (1)

and integrate we get
X X
σ σ
Fµν Cνj = σj σ
Sµν Cνj
ν ν
where S is the overlap matrix of the basis functions and Fσ is the matrix
representation of F σ
Z
Fµν = dr1 χ?µ (1)F σ (1)χν (1)
σ
Therefore we end up with two coupled matrix equations, the Pople–Nesbet

equations
Fα Cα = SCα α
Fβ Cβ = SCβ β
Unrestricted density matrices

The charge density of the system is
Nα Nβ
X X
T
ρ (r) = |ψaα (r)|2 + |ψaβ (r)|2 = ρα (r) + ρβ (r)
a=1 a=1
The corresponding spin density is defined by
ρS (r) = ρα (r) − ρβ (r)
By substituting the basis set expansion for the orbitals we get

Nα
X XX
α
ρ (r) = |ψaα (r)|2 = α
Pµν χµ (r)χν (r)
a=1 µ ν
and
Nβ
X XX
β
ρ (r) = |ψaβ (r)|2 = β
Pµν χµ (r)χν (r)
a=1 µ ν
83
where the α and β electron density matrices are used
Nα
X
α α α ?
Pµν = Cµa (Cνa )
a
Nβ
X
β β β ?
Pµν = Cµa (Cνa )
a
In addition to these two density matrices we can also define a total and spin
density matrix
PT = P α + P β
PS = P α − P β
With these definitions we can write the Fock matrices in compact form
XX
α
Fµν = hcore
µν + T
Pλσ α
(µν, σλ) − Pλσ (µλ, σν)
λ σ
XX β
β
Fµν = hcore
µν + T
Pλσ (µν, σλ) − Pλσ (µλ, σν)
λ σ
and the total energy is

1 X X T core α α β β

E0 = Pµν hµν + Pµν Fµν + Pµν Fµν
2 µ ν
84
Lecture 7
Correlation Energy
7.1 Basics
For many systems the Hartree–Fock method accounts for about 99 % of the
total energy. However, often the remaining 1 % are important for physical
and chemical properies of the system.
Correlation:
Reduced probability to find an electron with different spin close to another
electron: Coulomb hole.
Reduced probability to find an electron with the same spin close to another
electron: Fermi hole.
Hartree–Fock: Variational optimal one determinant wave function.
Improvements are only possible with multi determinant wave functions.
X
Ψ = a0 ΦHF + a i Φi
i
Excited Slater determinants:

Hartree–Fock determinant
ΦHF = |Φ1 Φ̄1 · · · ΦN/2 Φ̄N/2 |
Not used MOs (virtual, unoccupied MOs)
ΦN/2+1 , · · · , ΦK
Produce new determinants by exchange of MOs in ΦHF with one, two, three,
. . . virtual MOs. Singly, doubly, triply, . . . excited determinants.
Φ0 , Φra , Φrs rst

ab , Φabc , . . .
85
This procedure generates determinants with all possible eigenvalues of Ŝz ,
linear combination of these determinants that are also eigenfunctions of Ŝ 2
are called CSF (Configuration State Functions).
The complete optimization problem to find an exact wave function with this
method is two dimensional. Enlargement of the one-electron basis to get
a converged Hartree–Fock wave function also increases the number of vir-
tual MOs. This means an increased number of excited determinants. The
complete solution for a given basis set (with K functions) is called Full Con-
figuration Interaction (FCI).
Basis exact
complete exact HF solution
minimal
HF FCI
86
7.2 Matrixelements
Examples of determinants with two orbitals
1
Φ0 : √ [a(1)b(2) − b(1)a(2)]
2
1
Φra : √ [r(1)b(2) − b(1)r(2)]
2
1
Φrs
ab : √ [r(1)s(2) − s(1)r(2)]
2
1
hΦ0 | Ô | Φ0 i = ha(1)b(2) − b(1)a(2) | Ô | a(1)b(2) − b(1)a(2)i
2
1 1
= ha(1)b(2) | Ô | a(1)b(2)i − hb(1)a(2) | Ô | b(1)a(2)i
2 2
1 1
− ha(1)b(2) | Ô | b(1)a(2)i + hb(1)a(2) | Ô | a(1)b(2)i
2 2
Operators:
Ô0 = 1
N
X
Ô1 = ĥ(i)
i=1
N XN
X e2
Ô2 =
i=1 j>i
rij
Case 1:
hΦ|Φi = 1
X
hΦ|Ô1 |Φi = hi|ĥ|ii
i
1 XX
hΦ|Ô2 |Φi = (ii|jj) − (ij|ji)
2 i j
Case 2:
hΦ|Φra i = 0
hΦ|Ô1 |Φra i = ha|ĥ|ri
X
hΦ|Ô2 |Φra i = (ar|jj) − (aj|jr)
j
87
Case 3:
hΦ|Ô1 |Φrs
ab i = 0
hΦ|Ô2 |Φrs
ab i = (ar|bs) − (as|br)
Case 4:
hΦ|Ô2 |Φrst
abc i = 0
Energies of determinants:
X
hΦ | Ĥ | Φi = Hii + EJ + EK
i
hi i : contribution from each electron

EJ : contribution from each unique electron pair
EK : contribution from each unique electron pair with same spin
Examples
Orbital
2
a b c d
a) h11 + h22 + J12 − K12

b) h11 + h22 + J12
c) 2h11 + J11
d) h11 + 2h22 + 2J12 + J22 − K12
7.3 Configuration interaction method (CI)
Ansatz X X
ΨCI = a0 ΦHF + a S Φs + a D ΦD + · · ·
S D
Variational calculation with normalisation constraint
hΨCI |ΨCI i = 1
88
Lagrange function
L = hΨCI |Ĥ|ΨCI i − λ [hΨCI |ΨCI i − 1]
XX
hΨCI |Ĥ|ΨCI i = ai aj hΦi |Ĥ|Φj i
i=0 j=0
X XX
= a2i Ei + ai aj hΦi |Ĥ|Φj i
i=0 i=0 j6=j
XX
hΨCI |ΨCI i = ai aj hΦi |Φj i
i=0 j=0
X
= a2i hΦi |Φi i
i=0
X
= a2i
i=0
Variation:
∂L X
=2 aj hΦi |Ĥ|Φj i − 2λai = 0
∂ai j
X
ai hΦi |Ĥ|Φj i − λ + aj hΦi |Ĥ|Φj i = 0
j6=i
X
ai (Ei − λ) + aj hΦi |Ĥ|Φj i = 0
j6=i
Multiplication of the last equation with ai and summation over all equations
leads to X X XX
a2i Ei − a2i λ + ai aj hΦi |Ĥ|Φj i = 0
i i i j6=i
and therefore
λ = hΨCI |Ĥ|ΨCI i .
The Lagrange multiplier is the CI energy (λ = ECI ). We have a CI eigenvalue
equation in matrix form
Ha = Ea
Size of the CI matrix:

N electorns (indistinguishable) in K MOs.

K K!
Number of Determinants = =
N N !(K − N )!
89
Many determinants have a wrong (not desired) MS and can be neglected. In
addition the CI matrix is reduced in the CSF basis
K!(K + 1)!
Number of CSF(Singlet) =
(N/2)!(N/2 + 1)!(K − N/2)!(K − N/2 + 1)!
Examplel:
H2 O 10 electrons in 38 AOs (= 38 MOs)
Excitation CSF total CSF

0 1 1
1 71 72
2 2485 2557
3 40040 42597
4 348530 391127
5 1723540 2114667
6 5033210 7147877
7 8688080 15836557
8 8653645 24490202
9 4554550 29044752
10 1002001 30046752
Sizes:
CI-vector ≈ 30 · 106 elements (≈ 240 Mbytes)
CI-matrix ≈ (30 · 106 )2 = 900 · 1012 elements (≈ 7200 Tbytes)
but most Hij = 0, number of non-zero matrix elements is proportional to the
matrix size.
Solution of the CI eigenvalue equation:

Only lowest eigenvalues are of interest.
Iterative diagonalization for some of the lowest eigenvalues
• Lanczos method
• Davidson method
• Residuum methode
The CI matrix is sparely occupied. The number of matrix elements Hij 6= 0

is proportional to the dimension of the matrix. H is not explicitely calculated
but only Ha.
90
Reduced CI methods:
Idea: Not all determinants are equally important.
Ansatz: Only allow excitations from a subset of MOs into a subset of virtual
MOs (active space).
Allow only a maximal number of excitations.
Ne-Atom: Weight of CSF
Excitation Weight
0 0.9644945573
1 0.0009864929
2 0.0336865893
3 0.0003662339
4 0.0004517826
5 0.0000185090
6 0.0000017447
7 0.0000001393
8 0.0000000011
Active Space: restriction of orbitals for excitations
active
space
91
7.4 Dynamical and non-dynamical electron
correlation
For the majority of molecules, for example all moleculaes that can be assigned
a single Lewis structure, the main error from in the Hartree–Fock approxi-
mation comes from ignoring the correlated motion of each electron with all
the other electrons. This is called dynamical correlation because it refers to
the dynamical character of the electron-electron interactions. Even though
the correlation energy may be large, the Hartree–Fock determinant domi-
nates the configuration interaction expansion in these cases. The correction
is made up from many small contributions from excited state determinants.
In certain cases, however, it happens that a set of nearly degenerate orbitals
occure in a molecule that can only be partally occupied. In a restricted
Hartree–Fock calculation one or several of these orbitals cannot be included
into the wavefunction. An equally good wavefunction could be generated by
including these orbitals and exculding the others. In a subsequent config-
uration interaction calculation we expect that these determinants will have
similar weights. The correlation energy gained for these cases is called non-
dynamical correlation.
92
Lecture 8
Coupled Cluster Approaches
8.1 Size Consistency
F1 F2
Method X is size consistent if
EX (R = ∞) = 2EX (F )
Example: Hartree-Fock
Ψr=∞ = ΨF1 · ΨF2
therefore E(R = ∞) = E(F1 ) + E(F2 ) = 2E(F )

Example: SDCI
Minimal basis set (H2 )2
93
H2 H2
r s
R
a b
Determinants for the full calculation:

HF aābb̄
singly ar̄bb̄
singly aābr̄
CT as̄bb̄
For R = ∞ the CT determinants are zero and we are left
CT aābr̄
doubly rr̄bb̄
doubly aāss̄
CT as̄br̄
with 5 determinants.
For the calculation of the separated systems we have three determinants per
system.
HF aā
singly ar̄
doubly rr̄
The total wavefunction will be a linear combination of all the
HF bb̄
singly bs̄
doubly ss̄
products. We get 9 determinants, including triply and quadruply excited
ones. The energy of these two calculations will certainly differ. Truncated
CI methods are not size consistent.
2E(F ) < E(R = ∞)
The bigger the system the smaller the percentage of correlation energy is
included in a truncated CI calculation.
ECorr (SDCI)/ECorr(F CI) −→ 0 for Nel → ∞
A related concept is size extensivity. A method is size extensive if
ECorr (FN ) ∝ N · ECorr (F )
94
8.2 Coupled Cluster Methods
Perturbation theory based methods add all corrections (singly, doubly, triply
. . . excitations) up to a given order to the reference function. In Coupled
Cluster (CC) methods for a given type of correction all order will be included.
Ansatz for the CC wave function
ΨCC = eT̂ Φ0
∞
X 1 k
eT̂ = T̂
k=0
k!
The operator T̂ is defined as
T̂ = T̂1 + T̂2 + T̂3 + · · · + T̂N ,
where the operator T̂i applied to the Hartree–Fock wave function generates
all i-th excited Slater determinants.
occ X
X vir
T̂1 Φ0 = tra Φra
a r
occ X
X vir
T̂2 Φ0 = trs rs
ab Φab
a<b r<s
The expansion coefficients t are called amplitudes. They are our new varia-
tional parameters.
1 1
eT̂ = 1 + T̂1 + (T̂2 + T̂12 ) + (T̂3 + T̂2 T̂1 + T̂13 )
2 6
1 1 1
+ (T̂4 + T̂3 T̂1 + T̂22 + T̂2 T̂12 + T̂44 ) + · · ·
2 2 24
Insertion into the Schrödinger equation:
ĤeT̂ Φ0 = EeT̂ Φ0
multiplication from left with Φ∗0 and integration
hΦ0 |ĤeT̂ Φ0 i = ECC hΦ0 |eT̂ Φ0 i

= ECC hΦ0 |(1 + T̂1 + T̂2 + · · · )Φ0 i
= ECC
95
We make use of the fact that the Hamilton operator only contains one- and
two-electron terms:
1
ECC = hΦ0 |Ĥ(1 + T̂1 + T̂2 + T̂12 )Φ0 i
2
= hΦ0 |ĤΦ0 i + hΦ0 |Ĥ T̂1 Φ0 i
1
+ hΦ0 |Ĥ T̂2 Φ0 i + hΦ0 |Ĥ T̂12 Φ0 i
2
Xocc X
vir
= E0 + tra hΦ0 |ĤΦra i
a r
occ X
X vir
+ (trs r s s r rs
ab + ta tb − ta tb ) hΦ0 |ĤΦab i
a<b r<s
For Φ0 = ΦHF we have hΦ0 |ĤΦra i = 0 (Brillouins theorem) and the second
integral can easily be written with MOs:
occ X
X vir
ECC = E0 + (trs r s s r
ab + ta tb − ta tb ) [(ar|bs) − (as|br)]
a<b r<s
The CC energy is completely determined by the Hartree–Fock energy and

the one and two electron amplitudes.
Coupled (nonlinear) equations for the amplitude can be generated from the
Schrödinger equation by multiplication from left with Φr∗ rs∗
a , Φab , . . . and sub-
sequent integration.
hΦra | −→ Nocc · Nvirt equations
rs 2 2
hΦab | −→ Nocc · Nvirt equations
rst
hΦabc | −→ · · ·
8.3 Approximate Coupled Cluster methods

For the exact CC method all operators T̂i (i ≤ N ) have to be used.
T̂ = T̂1 + T̂2 + T̂3 + · · · T̂N
The method scales like M 2N +2 for M basis functions and N electrons. The
only practical method is to use approximations of the form T̂ ≈ T̂1 + T̂2
(CCSD).
T̂ = T̂1 + T̂2 CCSD M6
T̂ = T̂1 + T̂2 + T̂3 CCSDT M8
T̂ = T̂1 + T̂2 + T̂3 + T̂4 CCSDTQ M 10
96
The case T̂ = T̂1 does not improve on the Hartree–Fock energy. Multiplica-
tion of the Schrödinger equation from left with Φ∗0 , Φr∗ rs∗
a , Φab and integration
leads for CCSD to a system of coupled, non-linear equations for the ampli-
tudes tra und trs
ab . These equations have to be solved iteratively. The most
common improvement over the CCSD method is achieved by a perturba-
tive correction for T̂3 . Inclusion of T̂3 with perturbation theory leads to the
CCSD(T) method with M 7 scaling.

T̂1 +T̂2 1 2
e = 1 + T̂1 + T̂2 + T̂1
2

1 3 1 2 1 2 1 4
+ T̂2 T̂1 + T̂1 + T̂ + T̂2 T̂1 + T1 + · · ·
6 2 2 2 24
CCSD equations

1 2
ECCSD = E0 + hΦ0 | Ĥ T̂2 + T̂1 Φ0 i
2

1 1
hΦra | Ĥ 1 + T̂1 + T̂2 + T̂12 + T̂2 T̂1 + T̂13 Φ0 i = ECCSD hΦra | T̂1 Φ0 i
2 6

rs rs 1 2
hΦab | Ĥ (· · · ) Φ0 i = ECCSD hΦab | T̂2 + T̂1 Φ0 i
2
CCSD is non-variational. CCSD is size-consistent; higher excitations are

included through products (trs tu
ab tcd ).
97
Table 8.1: H2 O at the basis set limit
Exp. CCSD(T) HF
rOH [pm] 95.72 95.64 94.0
α(HOH) 104.52 104.2 106.2
µe [D] 1.8473 1.853 1.936
Atomization energy [kJ/mol] 975.3 975.5 652.3
Table 8.2: Basis set error: Statistics over many small molecules
Basis HF CCSD - HF CCSD(T) - CCSD

cc-pVDZ 5.8 23.8 1.3
cc-pVTZ 3.9 6.3 0.5
cc-pVQZ 1.5 1.5 0.2
cc-pV5Z 0.2 0.5 0.0
cc-pV6Z 0.0 0.4 0.0
8.4 Accuracy in Coupled Cluster calculations
property
experiment
intrinsic error basis set
limit
basis
98
Lecture 9
Møller–Plesset Perturbation
Theory
9.1 Electron correlation with perturbation

theory
Problem: What is Ĥ (0) and what is Ĥ (1) ?
Perturbation theory
Ĥ = Ĥ (0) + Ĥ (1)
(0)
Basis Ψn
Ĥ (0) Ψ(0) (0) (0)
n = E n Ψn
(0)
First order energy correction to ground state Ψg
Eg(1) = hΨ(0)
g | Ĥ
(1)
| Ψ(0)
g i
Second order energy correction
Eg(2) = hΨ(0)
g | Ĥ
(1)
| Ψ(1)
g i
and first order correction to wavefunction

(1)
X Hig (0)
Ψ(1)
n = (0) (0)
Ψi
i6=g Eg − Ei
and we have
(1)
X |Hig |2
Eg(2) = (0) (0)
i6=g Eg − E i
99
with
(1) (0) (0)
Hig = hΨi | Ĥ (1) | Ψj i
Definition of Hamiltonian for the unperturbed system

X
Ĥ (0) = F̂ (i)
i
The unperturbed Hamiltonian is chosen to be the sum of all Fock-operators

of all electrons of the system. This is a sum of one-electron operators with-
out coupling. Therefore the solution will be the sum of all energies of the
independent systems and the wavefunction will be the product of the one-
electron wavefunctions (orbitals). Because of the Pauli-principle we have to
anti-symmetrize this product and get a determinantal wavefunction.
Ψ(0)
g = |Φ1 Φ2 · · · ΦN | = ΦHF
The complete basis of all zeroth-order solution are all possible determinants
generated from the orbitals.
(0)
Ψi = |Φ1 Φ2 · · · Φr · · · ΦN | = Φra
and all other excited determinants (Φrs rst

ab , Φabc , . . .).
Energy
X
Eg(0) = i 6= EHF
i
The energy of the other determinants (basis functions) is also just the sum
of the orbital energies.
(0)
Ei = 1 + 2 + 3 + · · ·
Perturbation Hamiltonian
X
Ĥ (0) = F̂ (i)
i
X X
Ĥ (1) = Ĥ − Ĥ (0) = ĥ(i) + V̂ee − F̂ (i)
i i
X
= V̂ee − (Jˆj − K̂j )
j
100
Møller–Plesset perturbation theory
X
MP0 = E(MP0) = hΦHF |Ĥ (0) |ΦHF i = i
i
X
MP1 = E(MP0) + E(MP1) = i + hΦHF |Ĥ (1) |ΦHF i
i
X 1X e2 X
= i + hΦHF | |ΦHF i − ((ii|jj) − (ij|ij))
i
2 ij rij ij
X 1X X
= i + (Jij − Kij ) − (Jij − Kij )
i
2 ij ij
X 1X
= i − (Jij − Kij )
i
2 ij
= E(HF )
MP2 = E(HF) + E(MP2)
From general perturbation theory:
X hΦ0 |Ĥ (1) |Φi ihΦi |Ĥ (1) |Φ0 i
E(MP2) =
i6=0
E0 − E i
With
hΦ0 |Ĥ (1) |Φra i = 0 (Brillouin Theorem)
hΦ0 |Ĥ |Φrs
(1)
ab i = (ab|rs) − (ar|bs)
hΦ0 |Ĥ (1) |Φrst
abc i =0
and the energy denominator
E0 − E(Φrs
ab ) = a + b − r − s
we get
occ X
vir
X [(ab|rs) − (ar|bs)]2
E(MP2) =
a<b r<s
a + b − r − s
MP perturbation theory is size-consistent and size-extensive.
Order Scaling % correlation
1 M4 0
2 M5 ≈ 90
3 M6 ≈ 95
4 M7 ≈ 98
5 M8 -
6 M9 -
101
Table 9.1: Water, DZ basis set, RHF calculation
Method re 1.5re 1.5re

MP2 92.9 86.9 80.6
MP3 95.5 88.0 78.5
MP4 99.2 96.8 94.0
MP5 99.6 97.5 95.4
MP6 99.9 99.3 99.6
MP7 99.95 99.5 99.9
MP8 99.98 99.9 101.3
FCI -0.130085 -0.191996 -0.299864
MP2 efficient (≈ HF) implementation possible

recovers ≈ 90 % of correlation energy
MP3 often not more accurate than MP2

→ not a standard method
MP4 improves on MP2

similar cost, but inferior to CCSD(T)
MP5· · · only for benchmarking

Only MP2 is a commonly used method in quantum chemistry.
9.2 Convergence of MP series

No smooth convergence, or no convergence at all is possible.
Detection of difficult cases: HF determinant is not a good starting point.

Weight of Φ0 in perturbation wavefunction is less than 90 %. Example: H2 O
at a 2.0re geometry. The HF determinant has only a 58% weight in MP2.
The wavefunction has multi-configuration character.
MP based on unrestricted Hartree–Fock
• Slow convergence due to spin contamination problems
• To be used with care! Better use spin-projection methods.
102
MP based on multi-configurational wavefunctions (e.g. CAS-PT2)
• expensive
• often very accurate
• needs expert knowledge
9.3 The dissociation problem
Restricted vs. unrestricted solutions.
Example: H2 −→ H + H, homolytic bond breaking

Minimal basis Φ1 , Φ2 with overlap S12 = hΦ1 | Φ2 i
Ψ1 = [2(1 + S12 )]−1/2 (Φ1 + Φ2 )

Ψ2 = [2(1 − S12 )]−1/2 (Φ1 − Φ2 )
Symmetry therefore fully defines the restricted solution Ψ0 = |Ψ1 Ψ̄1 |.

Unrestricted case
Ψα1 = cos Θ Ψ1 + sin Θ Ψ2

Ψβ1 = cos Θ Ψ1 − sin Θ Ψ2
if Θ = 0 then RHF=UHF.
103
α β
Ψ = Ψ
small r
α β
Ψ Ψ
large r
Energy
RHF
UHF
Instability point
Bond breaking (homolytic)

• restricted wavefunctions fail
• unrestricted wavefunctions may have spin-contamination
• MP2 is highly affected
• CCSD(T) might solve problems partially
Ultimate solution: Multi-configurational SCF (MC-SCF)
|Ψ0 i = a0 Φ0 + a1 Φ1 + · · ·
The wavefunction is approximated by a small CI expansion. Variationally
optimize orbitals and CI-coefficients at the same time.
104
Lecture 10
Density Functional Theory:

Part I
10.1 Basics
In the following we investigate N -electron systems within the Born–Oppenheimer
approximation. The system is described by a Hamilton operator Ĥ
Ĥ = T̂ + V̂ext + V̂ee ,
where T̂ is the operator of the kinetic energy, V̂ee the electron-electron inter-
action, and V̂ext the external potential.
10.1.1 Hohenberg–Kohn Theorems

Theorem 1
The external potential V̂ext is, except for a trivial constant, uniquely defined
by the electron density ρ(~r).
Proof:
0
Assumption: There are two potentials V̂ and V̂ compatible with the electron
density ρ(~r).
1) ρ −→ V̂ −→ Ĥ −→ Ψ
0 0 0
2) ρ −→ V̂ −→ Ĥ −→ Ψ
Using the variation principle for wavefunctions we get
0 0 0 0 0 0 0 0
E0 < hΨ |Ĥ|Ψ i = hΨ |Ĥ |Ψ i + hΨ |Ĥ − Ĥ |Ψ i
Z
0 0
= E0 + ρ(r) V̂ (r) − V̂ (r) dr
105
and
0 0 0
E0 < hΨ|Ĥ |Ψi = hΨ|Ĥ|Ψi + hΨ|Ĥ − Ĥ|Ψi
Z
0
= E0 − ρ(r) V̂ (r) − V̂ (r) dr .
0 0
We have used |Ψi and |Ψ i as test functions with Ĥ and Ĥ and have made
use of the variation principle. After adding up the two equations we get the
following contradiction
0 0
E0 + E 0 < E 0 + E 0 ,
and therefore our initial assumption must be wrong. Since we have used the
variation principle this proof holds only for the ground state.
Total energy
Z
E0 = EV (ρ0 ) = T [ρ0 ] + Vext (r)ρ0 (r)dr + J[ρ0 ] + ENC [ρ0 ]
Classical Coulomb-energy:
Z Z
1 ρ0 (r1 )ρ0 (r2 )
J[ρ0 ] = dr1 dr2
2 |r1 − r2 |
Kinetic energy: T [ρ0 ]

Non-classical electron-electron interaction energy: ENC [ρ0 ]
Only Vext directly depends on the system, therefore we can write

Z
E0 = Vext (r)ρ0 (r)dr + FHK [ρ0 ]
where FHK is a universal functional of the electron density:
FHK [ρ] = T [ρ] + J[ρ] + ENC [ρ] .
Theorem 2: R
For a trial density ρ̃(r), with the properties ρ̃(r) ≥ 0 and ρ̃(r)dr = N
(number of electrons) we always have E0 ≤ EV (ρ̃).
Proof:
0 0 0
ρ̃ −−−→ V → Ĥ → Ψ
HK1
106
0
We use |Ψ i as a trial function with the Hamilton operator Ĥ. We have
(variation principal):
0 0
E0 ≤ hΨ |Ĥ|Ψ i = EV [ρ̃]
and therefore
E0 ≤ EV [ρ̃] .
0
This proof is based
R on the assumption that ρ̃ → V is valid. But not all
ρ̃(r) ≥ 0 with ρ̃(r)dr = N have a corresponding potential.
The set of all electron densities with a corresponding potential are called
v-representable.
The variation principle
E0 = Minρ̃ {EV [ρ̃]}
is only valid for v-representable electron densities.
The proof can be expanded to N -representable densities (Levi constraint
search method) A N -representable density is defined by:
• calculated from an anti-symmetric wavefunction
• fulfils Z Z
ρ(r) ≥ 0; ρ(r)dr = N ; |∇ρ1/2 |2 dr < ∞
Plausibility arguments by Wilson
1. Starting from ρ(r)
2. Calculate the number of electrons N in the system

Z
N = ρ(r)dr
3. Find the position and charge of the nuclei (cusp-condition)

∂
ρ̄(rA ) = −2ZA ρ̄(0) .
∂rA rA =0
where ρ̄(rA ) is the over all space directions averaged density.

P ZA e2
This defines the external potential Vext (r) = A |r−R A|
4. N and Vext define the Hamilton operator.
5. With Ĥ we know Ψ and therfore all properties of the system can be

calculated.
107
The electron density ρ(r) defines all properties of an N -electron system.
The Hohenberg-Kohn theorems give an existence proof for the functionals
T [ρ] and ENC [ρ], as well as a variation principle for the case of the known
functionals. What is missing is the concrete form of the functionals. Many
approximate functionals have been proposed. The accuracy, especially for
the kinetic energy is for applications in chemistry not sufficient. In real
application we will use another method.
10.1.2 Kohn–Sham Theory

We discuss a system with N non-interacting electrons in a local potential V̂s :
Ĥs = T̂ + V̂s .
From the 1. Hohenberg–Kohn theorem we get

Z
Es [ρ] = Ts [ρ] + Vs (r)ρ(r)dr
and with the 2. Hohenberg–Kohn theorem it follows that
δEs [ρs ] = 0 ,
This allows us to calculate the ground state density ρs .

The solution of the Schrödinger equation
Ĥs Ψs = Es Ψs
is a one-determinant wave function. This is due to the non-interacting elec-

trons that allow for a simple product wavefunction for the solution. At the
same time we also want the Pauli-principle to hold
Ψs = |Φ1 (r) · · · ΦN (r)|
The one-particle density matrix for this system is

N
X
ρs1 (r, r 0 ) = Φi (r)Φ∗i (r 0 ) ,
i=1
the density
N
X
ρs (r) = |Φi (r)|2 ,
i=1
108
and the kinetic energy
Z
1 2
Ts [ρ] = Φ∗i (r) − ∇ Φi (r)dr .
2
The functions Φi are orthonormal. With this is the functional Ts [ρ] for all
densities of the form
N
X
ρ(r) = |Φi (r)|2
i=1
defined.
Let’s now turn to a system of N interacting electrons and an external poten-

tial V̂0 . The corresponding energy functional is
Z
EV0 [ρ] = T [ρ] + V0 (r)ρ(r)dr + J[ρ] + ENC [ρ] .
Adding and subtracting Ts [ρ] we get

Z
EV0 [ρ] = Ts [ρ] + V0 (r)ρ(r)dr + J[ρ] + T [ρ] − Ts [ρ] + ENC [ρ]
Z
= Ts [ρ] + V0 (r)ρ(r)dr + J[ρ] + Exc [ρ]
This defines the exchange–correlation functional Exc [ρ]
Exc [ρ] = T [ρ] − Ts [ρ] + ENC [ρ] .
The corresponding potential is defined by the functional derivative

δExc [ρ]
Vxc [ρ](r) = .
δρ(r)
We assume that the optimal density ρ0 defined by
δEV0 [ρ]|ρ=ρ0 = 0
is non-interacting v-representable. Experience shows, that this is a reason-

able assumption for physical systems. With that we have a unique potential
Vs with
δEVs [ρ]|ρ=ρ0 = 0 .
Comparing the results we get
Z
ρ0 (r 0 ) 0
Vs (r) = V0 (r) + dr + Vxc [ρ0 ](r) ,
|r − r 0 |
109
the definition of the Kohn–Sham potential.
From the variation principle for the non-interacting system we get the eigen-
value equation
T̂s + V̂s Φi = i Φi ,
where the Lagrange multipliers i have been introduced for the orthogonality
constraint of the Kohn–Sham orbitals Φi . In analogy to the Hartree–Fock
equation the potential Vs depends on the solution of the Kohn–Sham equa-
tions and the Kohn–Sham orbitals have to be calculated using an iterative
procedure.
Exactly like the Hartree–Fock orbitals the Kohn–Sham orbitals can be ex-
panded in basis functions. We get a matrix representation of the Kohn–Sham
operator X
hµ|ĤKS |νi = hµν + Pλσ (µν, λσ) + hµ|V̂xc |νi ,
λσ
where P is the one-particle density matrix of the Kohn–Sham orbitals. For-

mally we can get the Kohn–Sham equations from the Hartree–Fock equations
by replacing the matrix representation of the exchange operator by the ma-
trix representation of the exchange–correlation potentials.
1X
− Pλσ (µλ, νσ) −→ hµ|V̂xc |νi
2 λσ
110
Lecture 11

Part II
11.1 Exchange and correlation energy

11.1.1 Density operators
For a system with the wave function ΨN (x1 , x2 , . . . , xN ) (xi are combined
spatial and spin coordinates) we define a density operator γN as
0 0 0 0 0 0
γN (x1 , x2 , . . . , xN , x1 , x2 , . . . , xN ) = ΨN (x1 , x2 , . . . , xN )Ψ?N (x1 , x2 , . . . , xN )
or in Dirac notation as a projection operator in position representation
γN = |ΨN ihΨN | .
The trace of γN is one
Z Z
Tr(γN ) = ··· ΨN Ψ?N dx1 · · · dxN = 1
and expectation values of operators Â can be calculated from

hAi = Tr(γN Â) = Tr(ÂγN ) .
For fermions we define reduced density operators of the form
Z Z
0 0 0 N
γp (x1 , x2 , . . . , xp , x1 , x2 , . . . , xp ) = · · · γN dxp+1 . . . dxN .
p
Specially we will need the two-particle density operator
Z Z
0 0 N (N − 1)
γ2 (x1 , x2 , x1 , x2 ) = · · · γN dx3 . . . dxN
2
111
and the one-particle density operator
Z Z
0
γ1 (x1 , x1 ) = N · · · γN dx2 . . . dxN .
Properties of γ1 and γ2 are

N (N −1)
• Tr(γ2 ) = 2
• Tr(γ1 ) = N
• γ1 and γ2 are positive semi-definite and Hermitian and can be diago-

nalised.
Z
0 0
γ1 (x1 , x1 )φi (x1 ) dx1 = ni φi (x1 )
Z Z
0 0 0 0
γ2 (x1 , x2 , x1 , x2 )θi (x1 , x2 ) dx1 dx2 = gi θi (x1 , x2 )
φi are the natural spin orbitals and θi the natural spin geminals of the system.
The eigenvalues ni and gi are called occupation numbers of the orbitals and
geminals, respectively. We have
0 ≤ n i , gi ≤ 1 .
and
X
ni = N
i
X N (N − 1)
gi =
i
2
Remember:
Orbital one-electron function
Geminal two-electron function
In the basis of the eigenfunctions we can write the density operators as

X
γ̂1 = ni |φi ihφi |
i
X
γ̂2 = gi |θi ihθi | .
i
112
Integration over spin variables:
Z
0 0
ρ1 (r1 , r1 ) = γ1 (x1 , x1 ) 0
ds1
s1 =s1
Z Z
0 0 0 0
ρ2 (r1 , r2 , r1 , r2 ) = γ2 (x1 , x2 , x1 , x2 )s1 =s0 ds1 ds2
1
0
s2 =s2
ρ1 and ρ2 are spin-free density operators. Information about spin has been
lost in ρ:
0 0
ρσ1 (r1 , r1 ) = γ1 (r1 σ, r1 σ)
ρ1 (r1 , r1 ) = ρα1 (r1 , r1 ) + ρβ1 (r1 , r1 )
0 0 0
Expectation values of operators that are not dependent on spin variables can
be calculated by
Z
0 0
hO1 i = Tr(O1 ρ1 ) = O1 (r1 )ρ1 (r1 , r1 ) 0
dr1
r1 =r1
Z Z
0 0
hO2 i = Tr(O2 ρ2 ) = O2 (r1 , r2 )ρ2 (r1 , r2 , r1 , r2 )r1 =r0 dr1 dr2
1
0
r2 =r2
The electronic Hamilton operator is a two-particle operator, therefore we

have
E = hHi = Tr(Hρ2 ) = E(ρ2 )
i.e. the total energy is a functional of the two-particle density operator.

Electron density:
ρ(r) = ρ1 (r, r)
Z
2 0 0
ρ(r) = ρ2 (r, r )dr
N −1
0 0 0
ρ2 (r, r ) = ρ2 (r, r , r, r )
11.1.2 Electron-electron interaction

We write the diagonal part of the two-particle density operator in the form
1
ρ2 (r1 , r2 ) = ρ(r1 )ρ(r2 )[1 + h(r1 , r2 )] ,
2
113
where h(r1 , r2 ) is the pair correlation function. The electron-electron inter-
action energy is
Z Z
1 1
Eee = ρ(r1 )ρ(r2 )[1 + h(r1 , r2 )] dr1 dr2
|r1 − r2 | 2
Z Z
1 1
= J[ρ] + ρ(r1 )ρ(r2 )h(r1 , r2 ) dr1 dr2
2 |r1 − r2 |
Z Z
1 1
= J[ρ] + ρ(r1 )ρxc (r1 , r2 ) dr1 dr2 .
2 |r1 − r2 |
where ρxc (r1 , r2 ) = ρ(r2 )h(r1 , r2 ) is the exchange–correlation hole of an elec-

tron at position r1 .
Sum rule:
Z
2
ρ(r1 ) = ρ2 (r1 , r2 )dr2
N −1
Z
2 1
= ρ(r1 )ρ(r2 )[1 + h(r1 , r2 )]dr2
N −1 2
Z Z
1
= ρ(r1 ) ρ(r2 )dr2 + ρ(r2 )h(r1 , r2 )dr2
N −1
Z
(N − 1)ρ(r1 ) = ρ(r1 ) N + ρxc (r1 , r2 )dr2
From this we get Z

ρxc (r1 , r2 )dr2 = −1
11.1.3 One-determinant wave functions

ΨS = |χ1 , χ2 , . . . , χN i
Reduced density operators
Z
0 0
γ1 (x1 , x1 ) = ΨS (x1 , . . .)ΨS? (x1 , . . .)dx2 . . . dxN
N
X
0 0
γ1 (x1 , x1 ) = χi (x1 )χ?i (x1 )
i=1
0 0 1h 0 0 0 0
i
γ2 (x1 , x2 , x1 , x2 ) = γ1 (x1 , x1 )γ1 (x2 , x2 ) − γ1 (x1 , x2 )γ1 (x2 , x1 )
2
γ2 is a simple function of γ1 .
114
Electron density:
N
X
ρ(r) = |χi (r)|2
i=1
Pair correlation function:
1
ρ2 (r1 , r2 , r1 , r2 ) = [ρ(r1 )ρ(r2 ) − ρ1 (r1 , r2 )ρ1 (r2 , r1 )]
2
Be careful with spin!
ρ1 (r1 , r2 )ρ1 (r2 , r1 ) = ρα1 (r1 , r2 )ρα1 (r2 , r1 ) + ρβ1 (r1 , r2 )ρβ1 (r2 , r1 )
ρ1 (r1 , r2 )ρ1 (r2 , r1 )
h(r1 , r2 ) = −
ρ(r1 )ρ(r2 )
Sum rule:
Z Z
ρ1 (r1 , r2 )ρ1 (r2 , r1 )
ρ(r2 )h(r1 , r2 )dr2 = − ρ(r2 ) dr2
ρ(r1 )ρ(r2 )
Z
1
=− ρ1 (r1 , r2 )ρ1 (r2 , r1 )dr2
ρ(r1 )
Z X
1
=− χi (r1 )χ?i (r2 )χj (r2 )χ?j (r1 )dr2
ρ(r1 ) ij
1 X
=− χi (r1 )χ?i (r1 )
ρ(r1 ) i
1
=− ρ(r1 )
ρ(r1 )
= −1
The exchange-correlation hole of a Slater determinant fulfils the sum rule

exactly.
Invariance to unitary transformations:
X
χ̄σj = χσi Uijσ
i
U σ is a unitary matrix (U U ? = 1). Spin symmetry can not be broken by the

transformation.
115
Example:
0
X 0
X
ρσ1 (r1 , r1 ) = χ̄σi (r1 ) χ̄σ?
i (r1 )
i i
X 0
= χσj (r1 )Ujiσ {χσk (r1 )Uki
σ ?
}
ijk
X 0
= χσj (r1 )Ujiσ Uki
σ? σ?
χk (r1 )
ijk
X 0
= χσj (r1 )χσ?
j (r1 )
j
With ρ1 is also ρ2 invariant and therefore the expectation values of the Hamil-
ton operator and the total energy.
Natural orbitals: X
γ1 = ni |φi ihφi |
i
and X
γ1 = |χi ihχi |
i
and therefore ni = 1 and φi = χi .

Idem-potency
( )( )
X X
γ1 γ1 = ni |φi ihφi | ni |φi ihφi |
i i
X
= n2i |φi ihφi |
i
X
= |φi ihφi |
i
= γ1
11.1.4 Exchange and correlation energy in the Kohn–

Sham method
Kohn–Sham exchange and correlation functionals:
EXC [ρ] = T [ρ] − Ts [ρ] + Vee [ρ] − J[ρ]
We want to make a connection of EXC [ρ] with the pair correlation function
from the electron-electron interaction energy. We consider Hamilton opera-
tors where the operator V̂ee has been multiplied by a scaling factor λ. For
116
λ = 1 we get a system of fully interacting electrons and for λ = 0 a system
of non-interacting electrons (similar to the Kohn–Sham reference system).
The wave function to these operators will be denoted by Ψλ . From the
Hohenberg–Kohn theorems we get for non-interactive v-representable ρ
Fλ=1 [ρ] = F [ρ] = T [ρ] + Vee [ρ]

Fλ=0 [ρ] = Ts [ρ] .
And therefore we have

= Fλ=1 [ρ] − Fλ=0 [ρ] − J[ρ]
Z 1
∂Fλ [ρ]
= dλ − J[ρ]
0 ∂λ
From a non-trivial derivation one gets finally

∂Fλ [ρ]
= hΨλ | Vee | Ψλ i
∂λ
and
Z 1
EXC [ρ] = dλhΨλ | Vee | Ψλ i − J[ρ]
Z0 Z
1
= ρ̄2 (r1 , r2 ) dr1 dr2 − J[ρ]
r
Z Z 12
1 1
= ρ(r1 )ρ(r2 )h̄(r1 , r2 ) dr1 dr2
2 r12
Z Z
1 1
= ρ(r1 )ρ̄xc (r1 , r2 ) dr1 dr2
2 r12
where the averaged exchange–correlation hole ρ̄xc is defined by
ρ̄xc (r1 , r2 ) = ρ(r2 )h̄(r1 , r2 )
and the averaged pair correlation function h̄

Z 1
1
dλ ρλ2 (r1 , r2 ) = ρ(r1 )ρ(r2 )[1 + h̄(r1 , r2 )]
0 2
ρλ2 (r1 , r2 ) is the diagonal part of the two-particle density matrix calculated
from the wave function Ψλ .
117
The averaged pair correlation function has the property
h̄(r1 , r2 ) = h̄(r2 , r1 )
and the sum rule
Z Z
ρ̄xc (r1 , r2 )dr2 = ρ(r2 )h̄(r1 , r2 )dr2 = −1 .
The exchange-correlation energy depends only on a spherically averaged

ρ̄xc (r1 , r2 ) Z Z ∞
1
EXC [ρ] = drρ(r) 4πsdsρSA xc (r, s) ,
2 0
where the spherically averaged exchange–correlation hole is defined by
Z
SA 1 0 0
ρxc (r, s) = ρ̄xc (r, r )dr
4π Ω
0
with the integration domain Ω : |r − r | = s. The sum rule gets
Z
4π s2 dsρSA xc (r, s) = −1 .
Definition of exchange and correlation energy in quantum chemistry:

EXC [ρ] = EXHF [ρHF ] + ECHF [ρ]
ECHF : Exact energy - Hartree–Fock energy.
EXHF : Hartree–Fock exchange energy calculated with Hartree–Fock orbitals.
For the Kohn–Sham theory we will use a slightly different separation:
EXC [ρ] = EXKS [ρ] + ECKS [ρ]
EXKS : Hartree–Fock exchange energy calculated with Kohn–Sham orbitals.
The exchange hole for spin compensated systems can be calculated from
1 |ρ1 (r1 , r2 )|
ρx (r1 , r2 ) = −
2 ρ(r1 )
and is Z
ρx (r1 , r2 )dr2 = −1 .
The correlation hole is defined as the difference to ρx
ρ̄xc (r1 , r2 ) = ρx (r1 , r2 ) + ρc (r1 , r2 )
and with that we have Z
ρc (r1 , r2 )dr2 = 0 .
The exchange energy is identical to the Coulomb interaction of the electrons
with a charge distribution with total charge one. The correlation energy is a
Coulomb interaction of the electrons with a neutral charge distribution.
118
Lecture 12

Part III
12.1 Local density approximation

Z
EXC [ρ] = drρ(r)XC [ρ(r)]
where XC (r) = f (ρ(r)) is a function (not functional) of the density.

Potential:
δXC [ρ(r)]
VXC (r) = XC (r) + ρ(r)
δρ(r)
The approximation is to have at each point in space the energy density with
the value of the homogenous electron gas of the local density.
Separation of XC energy density
XC (r) = X (r) + C (r)
Hartree–Fock exchange for the homogenous electron gas (Dirac 1930)

1/3
3 3
X (ρ(r)) = − ρ
4 π
Correlation energy for the homogenous electron gas:

Quantum-Monte-Carlo calculations (Ceperley, Alder 1980, Ortiz, Ballone
1994)
Analytic function fitted to QMC values:
1/3
3
rs =
4πρ
119
γc
C = √ rs > 1
1 + β 1 rs + β 2 rs
C = A ln rs + B + Crs ln rs + Drs rs < 1
γc = −0.1423, β1 = 1.0529, β2 = 0.3334

A = 0.0311, B = −0.048, C = 0.0020, D = −0.0116
The formula for EXC in the local density approximation gives
ρLDA
xc (r1 , r2 ) = ρ(r1 )h̄0 (|r1 − r2 |; ρ(r1 )) ,
where h̄0 (|r1 − r2 |; ρ(r1 )) is the averaged pair correlation function for the
homogenous electron gas. Contrary to the exact formula we have here the
density ρ(r1 ) as a factor and not ρ(r2 ). The LDA exchange correlation density
is Z
1 1 LDA
xc (ρ(r1 )) = ρ (r1 , r2 )dr2 .
2 r12 xc
ρLDA
xc is spherically symmetric and fulfils the sum rule
Z
ρLDA
xc (r1 , r2 )dr2 = −1 ,
but h0 (|r1 − r2 |; ρ(r1 )) is not symmetric in r1 and r2 . The exact function

ρxc is not spherical symmetric, but only the symmetric part is used in an
exchange correlation energy calculation. Therefore, a compensation of errors
can be expected from the LDA approximation.
12.2 Gradient corrections

Simple Taylor expansion don’t work, but ’empirical’ functionals of the type
Z
EXC = FXC [ρ, ∇ρ]dr
give clearly improved results.

Potential:
3
δFXC X ∂ δFXC
VXC (r) = −
δρ(r) α=1 ∂rα δ(∇α ρ)
120
12.2.1 Exchange functional of Becke (1988)
Z
EX = EXLDA [ρ(r)] − FX [s]ρ4/3 dr
|∇ρ|
s=
ρ4/3
s2
FX [s] = β 2−1/3
1 + 6βs sinh−1 s
β optimised using Hartree–Fock exchange energies of noble gas atoms.
β = 0.0042
Typical error in the exchange energy EXLDA ≈ 10% and EX (Becke) ≈ 0.1%.
12.2.2 Correlation functionals of Perdew et al.

(Perdew 1986, Perdew-Wang 1991, Perdew-Bourke-Ernzerhof (PBE) 1996)
These functionals model the correlation hole based on exact properties with-
out using empirical parameters.
The PBE functional
Z
PBE
Ec = ρ(r) [c + H(rs , t)] dr
r 1/3 1/3
|∇ρ| 4kf 9π 3
t= , ks = , kf = rs rs =
2ρks π 4 4πρ
2

β 1 + At
H(rs , t) = γ log 1 + t2
γ 1 + At2 + A2 t4
with
β
A = [exp(−c /γ) − 1]−1
γ
γ = 0.031090690869, β = 0.066724550
12.2.3 Lee-Yang-Parr (LYP) correlation functional (1988)

This functional is based on an approximative formula for the correlation en-
ergy of the Helium atom by Colle and Salvetti (1975). The basic assumption
is
ρ2 (r1 , r2 ) = CρHF
2 (r1 , r2 ) .
121
Then they replaced ρHF
2 with ρ1 and after a partial integration one gets
Z
LYP −a −cx

Ec = ρ 1 + bCF e dr
1 + dx
Z −cx

2 5 e 1 7 dx
+ ab|∇ρ| x 1+ cx + dr
1 + dx 24 3 1 + dx
where x = ρ−1/3 and
3
CF = (3π 2 )2/3 a = 0.04918 b = 0.132 c = 0.2533 d = 0.349
10
For |∇ρ| = 0 (i.e. a homogenous system) the correct LDA value is not
recovered. We have approximately
EcLDA ≈ 2EcLYP (|∇ρ| = 0)
The LYP functional has the special form
Fc (ρ, ∇ρ) = |∇ρ|2 G(ρ)
12.3 Hybrid functionals

Let’s go back to the definition of the exchange–correlation energy in the
Kohn–Sham method
= Fλ=1 [ρ] − Fλ=0 [ρ] − J[ρ]
Z 1
∂Fλ [ρ]
= dλ − J[ρ]
0 ∂λ
λ
and define the function EN C as the non-classical Coulomb energy at the value
λ. With this we get Z 1
λ
EXC [ρ] = EN C dλ .
0
For a value of λ = 0 (a non-interacting system) this energy is exactly known.
It’s the Hartree-Fock exchange energy, because the wavefunction for this
system is a single Slater determinant. If we now approximate the end point
of integration λ = 1 with a common density functional, we can estimate the
total integral with a weighted sum of the two end points. How to choose
the weights is again an open question. There are many possibilities, some of
them purely empirical and others based on theoretical considerations.
The two most often used hybrid functionals are
122
Table 12.1: Dependence on basis set for some bond length [Å]. LDA values
Bond 6-31G(d,p) 6-311++G(d,p) Limit Experiment
H −H –/0.765 –/0.765 –/0.765 –/0.741
H3 C − CH3 1.513/1.105 1.510/1.101 1.508/1.100 1.526/1.088
H2 C − CH2 1.330/1.098 1.325/1.094 1.323/1.093 1.339/1.085
HC − CH 1.212/1.078 1.203/1.073 1.203/1.074 1.203/1.061
• B3LYP
B3LYP
EXC = (1 − a)ExLDA + aExHF + bExB88 + cEcLYP + (1 − c)EcLDA
where ExB88 is the Becke exchange functional. The values for the pa-
rameters are
a = 0.20, b = 0.72, c = 0.81
• PBE0
PBE0 PBE
EXC = EXC + 0.25(ExHF − ExPBE )
This is a parameter free functional. Neither the PBE functional nor

the weight 0.25 are fit to any experimental data.
12.4 Basis sets

In density functional theory the same kind of basis sets can be used as in
Hartree–Fock calculations. In principal the parameters would have to be
newly optimised (e.g. for a 6-31G* type basis). In practice even individual
parameters for each density functional would have to be generated. However,
it turned out that such optimised parameters only have a minor effect on
results and the expense needed for the optimisation is not justifiable. It is
therefore common to use the same basis sets in density functional theory that
are also used in wavefunction calculations.
In Kohn–Sham density functional theory only simple orbitals have to be
described. Therefore, the demand on the basis is comparable to other inde-
pendent particle models like the Hartree–Fock model. Basis sets of double-ζ
quality with one set of polarisation functions are often sufficient to get good
results. Triple-ζ basis sets with two sets of polarisation functions are for
most cases close to converged results.
123
Table 12.2: Averaged and maximal absolute error of the atomisation energy
for a test set (G2) of molecules [kcal/mol]. Basis: 6-311+G(3df,2p)
Functional RMS max error
LDA 36.4 84
PBE 8.6 26
BLYP 4.7 15
B3LYP 2.4 10
PBE0 3.5 10
12.5 Model calculations

Comparing different functional for all kind of molecular properties is a favourite
past time of all quantum chemists. There is a plethora of research papers on
this topic.
12.6 Coulomb energy

Density functional programs are build exactly in the same way as Hartree–
Fock programs. The only difference is that exchange integrals are replace
by integrals with the potential from the XC functionals. Many DFT pro-
grams differ in another point from Hartree–Fock programs. In Hartree–Fock
program the Fock matrix elements for Coulomb and exchange are calculated
together.
X
c
Fµν = Pλσ (µν, λσ)
λσ
x 1X
Fµν =− Pλσ (µλ, σν)
2 λσ
Computational time is dominated by the calculation of the two-electron in-

tegrals. Therefore one uses
c
Fµν ←− 4Pλσ (µν, λσ)
c
Fλσ ←− 4Pµν (µν, λσ)
x
Fνσ ←− −Pµλ (µν, λσ)
x
Fµσ ←− −Pνλ (µν, λσ)
x
Fνλ ←− −Pµσ (µν, λσ)
x
Fµλ ←− −Pνσ (µν, λσ)
124
to achieve a minimal number of integral calculations. In DFT calculation
only the Coulomb part is needed, except for hybrid functionals. The Coulomb
contribution has a special structure that can be used together with further
approximations to increase the efficiency of the program.
For the calculation of the Coulomb part of the Kohn–Sham matrix we need
all integrals (µν, λσ) over all basis functions. This method scales like M 4
where M is the number of basis functions. For large systems the number
of integrals can be reduced drastically if only integrals with a significant
contribution to the energy are calculated. All modern programs are using
such screening procedures. It can be shown that asymptotically only M 2
integrals have to be calculated. However, this procedure has a rather large
pre-factor, i.e. the asymptotic scaling is only reached for very large systems.
For this reason other approximate method for the calculation of the Coulomb
energy are used.
12.6.1 Resolution of identity method

The electron density is defined as
X
ρ(r) = Pµν φµ (r)φν (r) .
µν
The product of two atom-centred basis functions is a new function

χ(µν) (r) = φµ (r)φν (r) .
The density is therefore described exactly in a basis of the size M 2 .
X
ρ(r) = Pµν χ(µν) (r)
(µν)
The basic idea is now to replace this exact description of the density by
an approximate description in a smaller basis, where the new basis is only
slightly bigger than the basis used for the orbitals.
X
ρ(r) ≈ cα χα (r)
α
With this we get for the contribution to the Kohn–Sham matrix

X Z Z φµ (r)φν (r)χα (r 0 )
c
Fµν = cα 0|
drdr 0
α
|r − r
X
= cα (µν, α)
α
125
This method scales like M 3 and asymptotically again like M 2 , but it has a
much smaller pre-factor. Special basis sets for the electron density have to be
used that match the basis sets for the wavefunctions. Typically these basis
sets have three times more functions than the wave function basis.
We haven’t yet discussed how the expansion coefficients cα are calculated. We
can transform this problem into an optimisation. We are looking for those
expansion coefficients which minimise the difference to the exact electron
density with respect to a given metric.
X
ρ̃(r) = cα χα (r)
α
∆(r) = ρ(r) − ρ̃(r)
Z Z
Ω({c}) = ∆(r)O(r, r 0 )∆(r 0 )drdr 0
c ← MIN{c} Ω({c})
where O(r, r 0 ) defines a distance metric. The solution of this problem is

X
−1
cα = Sαβ Tβ
β
with
Z Z
Sαβ = χα (r)O(r, r 0 )χβ (r 0 )drdr 0
Z Z
Tβ = χβ (r)O(r, r 0 )ρ(r 0 )drdr 0
1
The operators O(r, r 0 ) = δ(r − r 0 ), O(r, r 0 ) = |r − r 0 | and O(r, r 0 ) = |r−r 0|
have been investigated. Most of the programs use the Coulomb operator.
These methods are called RI (resolution of identity) method, because they
can be derived using the identity operator (for orthogonal basis sets).
X
1= |αihα|
α
This method is available in many quantum chemistry codes z.B. Turbomole,

Gaussian 2003, DeMON.
12.6.2 Local expansion of density

The local expansion of density are based on similar ideas as the RI methods,
but instead of looking for a global approximation to the density separate
126
approximations to all pair densities are used.
X
ρAB (r) = Pαβ Φα (r)Φβ (r)
αβ
α∈A,β∈B
X X
ρAB (r) ≈ c µ χA
µ (r) + c ν χB
ν (r)
µ∈A ν∈B
The pair density of atoms A and B is expanded in a basis located at the

atoms A and B. The advantage of this method is that the calculation of the
expansion coefficients is done locally and is therefore independent of the size
of the molecule. The disadvantage lies in the limitation of the expansion to
the functions on the same atoms. To achieve good results large basis sets
(with high angular momenta) have to be used. This method is used in the
program ADF.
12.7 Numerical integration

With Gaussian functions as basis sets, all integrals needed in standard quan-
tum chemistry calculations can be derived analytically. Integrals for the ex-
change and correlation functionals needed in density functional theory cannot
be calculated analytically even for special types of basis functions. They have
to be calculated using numerical procedures.
12.7.1 Voronoi polyhedra

We investigate a general integral over all space
Z
I= F (r)dr .
In a first step we divide the integral into different parts using a weight func-
tion wA (r). The function wA (r) should be close to unity near nucleus A and
close to zero near all other nuclei. The function should be continuous and
differentiable and defined in all space. We require further that
X
wA (r) = 1 .
A
127
With this we get
Z
I= F (r)dr
Z !
X
= wA (r) F (r)dr
A
XZ
= wA (r)F (r)dr
A
X
= IA
A
To calculate the weight functions we introduce confocal elliptic coordinates

(λ, µ, φ). These coordinates are defined by
r1 + r 2
λ=
R12
r1 − r 2
µ= .
R12
Dabei ist r1 die Distanz zum Zentrum 1 und r2 die Distanz zum Zentrum 2,
R12 der Abstand der beiden Zentren und φ der Winkel um die Verbindungslinie
der Zentren. Der Definitionsbereich der Koordinaten ist
φ: [0, 2π]
λ: [1, ∞]
µ: [−1, 1]
Die wichtige Koordinate ist µ. Für zwei Zentren A und B haben wir
rA − r B
µAB =
RAB
und wir sehen, dass µAB den Wert −1 am Zentrum A und auf der Achse von
A weg besitzt. µAB hat den Wert 1 am Zentrum B und auf der Achse von
B weg. Definieren wir eine neue Funktion s(µAB )

1 for − 1 ≤ µAB ≤ 0
s(µAB ) = .
0 for 0 < µAB ≤ 1
With this we can define the weigth functions
Y
PA (r) = s(µAB )
B6=A
PA (r)
wA (r) = P .
A PA (r)
128
The whole space is devided into Voronoi polyhedra with respect to the nu-
clei. The discontinuity of the function s(µAB ) at µAB = 0 is not suited for
numerical calculations. Therefore we choose
ds
= Am (1 − µ2 )m ,
dµ
where Am is taken such that s(−1) = 1 and s(1) = 0. A typical value for m
is 10.
It is advantageous to use also the relative size of the atoms, i.e. we make a
variable transformation such that different atoms get different weights.
νAB = µAB + aAB (1 − µ2AB )
The parameter aAB is calcuated from the Bragg-Slater radii (RA , RB ) of the
atoms.
uAB
aAB = 2
uAB − 1
χ−1
uAB =
χ+1
RA
χ=
RB
12.7.2 Radial integration

We calculate the integrals IA using numerical integration, where the integra-
tion is split into a radial part and an integration over the unit sphere.
X X
IA ≈ wr (i) wΩ (j) wA (r(i), Ω(j)) FA (r(i), Ω(j))
i j
There are many different integration methods for radial grids. The best meth-
ods for this problem work in two steps. First the integral [0, ∞] is mapped
onto [−1, 1] using a variable transformation. Possible transformations are
a+x
r = η√ ,
1 − x2
a+x
r=η ,
1−x
η a+1
r= ln ,
ln 2 1−x
where the transformations are normalised such that the midpoint of the x-
interval is maped onto r = η. The value of η depends on the type of atom.
129
H 0.8 He 0.9 Li 1.8 Be 1.4 B 1.3 C 1.1 N 0.9 O 0.9
F 0.9 Ne 0.9
Na 1.4 Mg 1.3 Al 1.3 Si 1.2 P 1.1 S 1.0 CL 1.0 Ar 1.0
K 1.5 Ca 1.4 Sc 1.3 Ti 1.2 V 1.2 Cr 1.2 Mn 1.2 Fe 1.2
Co 1.2 Ni 1.1
Cu 1.1 Zn 1.1 Ga 1.1 Ge 1.0 As 0.9 Se 0.9 Br 0.9 Kr 0.9
Finally we use a Gauss quadrature for the integral on the interval [−1, 1].
Z 1 N
X vj
f (x)dx ≈ f (xj )
−1 j=1
W (xj )
The form of the function W (x) depends on the orthogonal polynome that
was used as a base for the integration method.
Gauss-Legendre
W (x) = 1, −1 < x < 1
(j + 1)Pj+1 = (2j + 1)xPj − jPj−1
Gauss-Chebyshev
1
W (x) = √ , −1 < x < 1
1 − x2
Tj+1 = 2xTj − Tj−1
12.7.3 Angular integration

Two-dimensional integrals over the unit sphere
Z π Z 2π
f (θ, φ)sinθdθdφ
0 0
can be efficiently performed using Lebedev grids. Lebedev grids with a cer-
tain number of points allow for the exact integration of all spherical har-
monics up to a certain value of Lmax . These grids are based on octahedral
symmetry. The position of the special points and their integration weights
have to be taken from a table.
130
Gitter Lmax Inred Itot
1 3 1 6
2 5 2 14
3 7 3 26
4 9 3 38
5 11 4 50
6 15 5 86
7 17 6 110
8 19 7 146
9 23 9 194
10 29 12 302
11 35 16 434
12 41 20 590
13 47 25 770
14 53 30 974
131
132
Lecture 13
Molecular Properties
13.1 Wavefunction analysis

Wavefunction Ψ −→ chemical concepts
atomic charges
oxidation states
reactivity
bonding
13.1.1 Basisfunction based population analysis

P
For Ψ = |Φ1 · · · ΦN | or γ(r, r 0 ) = i fi Φi (r)Φi (r 0 ) (single determinant wave-
functions) we have
X X
ρ(r) = ρi (r) = fi |Φi (r)|2
i i
Basis set expansion

AO
X
Φi (r) = cαi φα (r)
α
Inserting this into the density we get
X
ρi (r) = fi |Φi (r)|2 = fi cαi cβi φα (r)φβ (r)
αβ
133
Calulating the total number of electrons
MO Z
X MO
X X Z
N= ρi (r) dr = fi cαi cβi φα (r)φβ (r) dr
i i αβ
MO
X X
= fi cαi cβi Sαβ
i αβ
X
= Pαβ Sαβ
αβ
= Tr (P S)
where we have used the density matrix
MO
X
Pαβ = fi cαi cβi
i
We will now look into schemes to partition the N electrons using (PS).
Mulliken population analysis

Pαα Sαα : fraction of electron in AO α
Now we sum over all α
Pαβ Sαβ : put 12 electron in α and 12 electron in β
associated with atom A. This only works with basis sets where the functions
can be easily subdivided using atoms.
AO X
X AO
ρA = Pαβ Sαβ
α∈A β
The total charge on atom A is then

qA = Z A − ρ A
Löwdin population analysis

This population analysis is based on a slightly different calculation of the
total number of electrons based on the invariance of the trace.
1 1

N = Tr (P S) = Tr P S 2 S 2
1 1

= Tr S 2 P S 2 = Tr P̄
X
= P̄αα
α
134
The number of electrons on atom A is
AO
X
ρA = P̄αα
α∈A
The total charge on atom A is again

qA = Z A − ρ A
These types of charges are
• strongly depend on basis functions
• can only be used for trends, not absolute values
• use only with small basis sets (no diffuse functions)
13.1.2 Electrostatic potential derived charges

Die electrostatic potential of an electronic charge distribution ρ(r) and a
collection of atomic cores with charge ZA is
Z
X ZA ρ(r 0 )
VESP (r) = −− 0
dr 0
A
|r − RA | |r − r |
The potential of a collection of atomic charges qA at the same positions RA
as the atomic cores is
q
X qA
VESP (r) =
A
|r − RA |
We now use the electrostatic potential to define the optimal atomic charges
qA with
• outside the charge distribution ρ(r) we require
q
VESP (r) = VESP (r)
• charge conservation
X
qA = total charge of molecule
A
• fix higher multipoles (dipole etc.)

We achieve this using a least square fit with constraints
X q 2
Minq (VESP (ri ) − VESP (ri )) + constraints
i
These type of charges are often used in force fields and QM/MM calculations.
135
13.1.3 Density based charges
Hirschfeld method
Definition Z
qA = Z A − ρ(r) dr
Atom
This requires a definition of an atomic volume. We start with a weighting

function
ρ̄A (r)
wA (r) =
ρ̄(r)
P
where ρ̄A (r) is the electron density of an isolated atom and ρ̄(r) = A ρ̄A (r).
The weighting function has the property
X
wA (r) = 1
A
for any choice of r. Therefore we can write

Z Z X
N = ρ(r) dr = ρ(r) wA (r) dr
A
XZ X
= ρ(r)wA (r) dr = ρA
A A
with the definition of ρA

Z
ρA = ρ(r)wA (r) dr
Atoms in molecules (AIM)

Bader analysis: based on the topology of ρ(r) (maxima, minima, saddle
points).
Maxima of ρ(r) are at nuclei = attractors of electronic charge. The gradi-
ent of density at each point in space points to strongest attractor (basin of
attraction = atom).
13.2 Response properties

Response of a molecule to external (internal) perturbation
• external electric field (F)
136
Table 13.1: Example: Atomic charge of oxygen in water
Basis set Mulliken Löwdin ESP AIM

STO-3G -0.39 -0.27 -0.65 -0.89
3-21G -0.74 -0.46 -0.90 -0.93
6-31G** -0.67 -0.44 -0.81 -1.24
6-311G(2d2p) -0.52 0.00 -0.74 -1.24
6-311++G(2d2p) -0.47 -0.12 -0.76 -1.25
• external magnetic field (B)
• nuclear magnetic moment (I)
• nuclear geometry (R)
time-independent perturbation : static property

time-dependent perturbation : dynamic property
Total property: contribution from electrons (QM) and nuclei
Taylor series expansion of energy
∂E 1 ∂2E 2
E(λ) = E(0) + λ+ λ +···
∂λ 2 ∂λ2
Energy derivatives = ”Properties”
Example: External electric field

Z
Eint = ρ(r)Vext (r) dr
Assume electric field F = − ∂V∂rext is uniform over molecule (Multipole expan-

sion)
1
Eint = qV − µ · F − Q : F 0 − · · ·
2
q = net charge
µ = dipole moment
Q = quadrupole moment
F 0 = field gradient
137
Field −→ wavefunction −→ dipole
1 1
µ = µ0 + αF + βF 2 + γF 3 + · · ·
2 6
µ0 = permanent dipole
α = polarisability
β, γ = hyper-polarisability
The total energy is

E(F ) = hΨ | H0 − µF | Ψi
with the total Hamiltonian
1
H = H0 − µ0 F − αF 2 − βF 3 − · · ·
2
Series expansion
∂E 1 ∂2E 2 1 ∂3E 3
E(F ) = E(0) + F+ 2
F + 3
F +···
|∂F
{z } |2 ∂F
{z } |6 ∂F
{z }
−µ0 −α − 12 β
General:
∂ n+m E
Property ∝
∂F n ∂Rm
n m Property
0 0 Energy
1 0 Electric dipole moment
0 1 Energy gradient
2 0 Polarisability
0 2 Harmonic vibrational frequency
1 1 Infrared absorption
2 1 Raman intensity
13.2.1 Derivative techniques

Hamiltonian
H = H0 + λP1 + λ2 P2
Expectation value
E(λ) = hΨ | H0 + λP1 + λ2 P2 | Ψi
138
First derivative
∂E(λ) ∂Ψ ∂Ψ
=h | H | Ψi + hΨ | P1 + λP2 | Ψi + hΨ | H | i
∂λ ∂λ ∂λ
For real wavefunctions we have λ → 0 then Ψ → Ψ0 .

∂E(λ) ∂Ψ
= hΨ0 | P1 | Ψ0 i + 2h | H0 | Ψ 0 i
∂λ λ=0 ∂λ λ=0
The wavefunction Ψ depends on two different set of parameters, the wave-

function expansion coefficients (c) and the basis functions (ψ). The derivative
wrt λ is therefore
∂Ψ ∂Ψ ∂ψ ∂Ψ ∂c
= +
∂λ ∂ψ ∂λ ∂c ∂λ
For fully variational wavefunctions we have
∂Ψ(α1 , α2 , . . . , αn
=0
∂αi
and we get the Hellmann-Feynman theorem
∂ ∂H
hΨ | H | Ψi = hΨ | | Ψi
∂λ ∂λ
For variationally optimised wavefunctions (HF, DFT, MCSCF, but not MPx,
CI, CCSD) we have
∂Ψ ∂Ψ
= 0 and 6= 0
∂c ∂ψ
and therefore

∂E(λ) ∂Ψ ∂ψ
= hΨ0 | P1 | Ψ0 i +2 h | H0 | Ψ 0 i
∂λ λ=0 |
{z } ∂λ λ=0
∂λ
Hellmann-Feynman | {z }
Pulay term
139
140
Lecture 14
NMR Chemical Shielding
14.1 General
Maxwell equations
~ = 4πρ
divE
~ =0
divB
~
1 ∂B
~ =−
rotE
c ∂t
~
~ = 4π J~ + 1 ∂ E
rotB
c c ∂t
~ with the property
Partial solution by introducing a vector field A
~ = rotA
B ~
A~ determines B~ uniquely, but A(

~ B)
~ is not unique. We can add a gradient
~ without change of the physical properties.
field to A
~0 = A
A ~ + gradλ
~ 0 = rotA
rotA ~ + rot(gradλ) = rotA
~
λ: gauge function
Insert into Maxwell equation
~
1 ∂B
~ =−
rotE
c ∂t
~
~ + 1 ∂ A = −gradΦ
E
c ∂t
141
Φ : scalar potential with
Φ0 = Φ + const.
~ Coulomb gauge
Special form of A:
~=0
divA
Poisson equation
∆Φ = −4πρ
2~
~ + 1 ∂ A + 1 grad ∂Φ = 4π J~
−∆A
c2 ∂t2 c2 ∂t c
In the non-relativistic limit (c → ∞) these equations simplify to
∆Φ = −4πρ
~=0
∆A
However, there is still a change in the gauge possible

~0 = A
A ~ + gradλ
as we get for
~ = div(A
divA ~ + gradλ)
= divgradλ
= ∆λ
This means that a change of gauge is possible for λ that fulfill the Laplace
equation
∆λ = 0
~
Special case of a homogeous external magnetic field B(r) =B ~0
~0 = 1 B
A ~ 0 × (~r − R)
~
2
1~ 1~ ~
= B 0 ×~r− B 0 ×R
2 2
1~
= B 0 ×~r + gradλ
2
1 ~ ~ r
λ = − (B 0 × R) · ~
2
~ is called the gauge origin.
The position R
142
14.2 Chemical shifts
~ loc (R
Nuclear spins couple to the local magnetic field B ~ A)
• level split
• transitions between levels can be induced by radiofrequency fields

~ loc (R
• measure energy difference ∝ B ~ A)
Local magnetic field

~ loc (~r) = B
B ~ ext (~r) + B
~ ind (~r)
~ ext (~r) = B
B ~0
Biot-Savart law Z ~
~ ind (~s) = − 1
B
j(~r) × (~r − ~s)
d~r
c |~r − ~s|
where ~j(~r) is the current density. We introduce the magnetic shielding con-
stant
~ ind (~s) = −σ(~s)B
B ~0
and the chemical shift
δ = σstd − σ(~s)
~ 0 and the
The chemical shift is independent on the external magnetic filed B
magnetic moment I of the nuclei.
Calculation: For molecule and standard calculate σ(~s) at atomic positions,

this needs the calculation of the current density ~j(~r).
Hamiltonian:
ext
~
B
H0 −−→H
~ ext
p~ → p~ + A
~ we have
For the Coulomb gauge and gauge origin R
Ĥ = H0 + A ~ ext p~ + 1 A
~ 2ext
2
1
~ ext = B~ ext × (~r − R) ~
A
2
~ ~ 1 ~

~ 1 ~ ~
2
Ĥ(Bext , R) = Ĥ0 + (~r − R) × p~ Bext + Bext × (~r − R)
2 8
143
For one-determinant wavefunctions and closed shells we have
Ψ = |Φ1 · · · Φn |
and for the current density
N
X N
X
~j = i ~ k
(Φ∗k ∇Φ ~ ∗ )Φk ) − 2A
− (∇Φ ~ ext Φ∗k Φk
k
k=1 k=1
14.3 Linear response for Hartree–Fock

Fock operator
2
ĥ = ĥ0 + iĥ1 + O(Bext )
Hartree–Fock orbitals
2
Φk = Φ0k + iΦ1k + O(Bext )
The current density vanishes for Bext = 0
~j = ~j 1 + O(Bext
2
)
We get for the first order current density
N
X N
X

~j 1 = 2 Φ1k ∇Φ0k − Φ0k ∇Φ1k ~ ext
− 2A Φ0k Φ0k
| k=1 {z } | k=1
{z }
paramagnetic part diamagnetic part
~ The paramagnetic and diamagnetic

~j 1 is independent of the gauge origin R.
parts depend on the gauge origin. This causes numerical problems if the two
terms are not calculated with the same accuracy. This is called the gauge
invariance problem of chemical shifts. It is effectively related to the problem
of non-complete basis sets. The ground state and linear response orbitals
have different needs for basis sets.
Solution: distributed gauge methods

~ k for every localised orbital
IGLO individual gauge origins R
GIAO gauge including atomic orbitals
CSGT continous set of gauge transformations
LORG localised orbital/local origin
144
Table 14.1: Example: chemical shifts for C calculated with different methods
(CH4 )
Basis common gauge origin IGLO GIAO

C H
TZ 25.58 62.01 32.48 32.57
TZ+3p 30.97 41.22 31.29 31.35
uncontracted 30.98 32.07 31.20 31.34
extended 31.23 31.26 31.24 31.25
14.3.1 Calculation of response orbitals

We expand the linear response orbitals in the set of ground state orbitals.
N
X
Φ1p = Yqp Φ0q
q=1
with Yqp = Ypq . Stationarity conditions for Hartree-Fock

i occupied
Fia = 0
o unoccupied
Expansion of Fock matrix
Fia = Fia0 + iFia1 + O(Bext

2
)
with X
Fia1 = hΦ01 | ĥ1 | Φ0a + Hia,jb Yjb
jb
and
0
Hia,jb = δij Fab − δab Fij0 + (ab, ij) − (aj, bi)
The condition Fia1 = 0 leads to a set of linear equations for Y of the form
HY = h
H depends on Φ0 and h depends on Φ0 and ĥ1 .
1 1 ~

ĥ = −i (~r − R) × p~ Bext
2
1 ~ ~

= − (~r − R) × ∇ Bext
2
145
For canonical orbitals Fpq = p δpq and we get
hΦ0i | ĥ1 | Φ0a i

Yia = −
i − a
The induced current density is
X X
j1 = 2 Yia Φ0i ∇Φ0a − Φ0a ∇Φ0i − B~ ext × (~r − R)
~ Φ0i Φ0i
ia i
Calculations of chemical shifts have been performed with Hartree–Fock, MP2,

CCSD(T), and DFT electronic structure methods. For heavy elements rela-
tivistic corrections are needed.
146

Computational Chem 6

Uploaded by

Copyright:

Available Formats

Computational Chem 6

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Computational Chem 6

Uploaded by

Copyright:

Available Formats

Lecture Notes

Electronic Structure Theory

August 11, 2005

2 Basic Mathematical Review 21

4 Two-Electron Systems and Spin 53

6 Molecular Orbital Theory 75

8 Coupled Cluster Approaches 93

9 Møller–Plesset Perturbation Theory 99

10 Density Functional Theory: Part I 105

11 Density Functional Theory: Part II 111

12 Density Functional Theory: Part III 119

13 Molecular Properties 133

14 NMR Chemical Shielding 141

Basic Quantum Mechanics

1.1 Operators in quantum mechanics

1.1.1 Eigenfunctions and eigenvalues

1.2 Postulates of quantum mechanics

classical quantum mechanical

Functions of position and momentum  

The positions r1 , r2 , . . . are associated with particles 1, 2, . . . . In general

[q, pq0 ] = i~δq,q0 , [q, q 0 ] = 0 , [pq , pq0 ] = 0 ,

The expectation value of an operator Ω̂ for a general state Ψ is defined by

If the wave function is normalized

we get for the expectation value

We will assume that state functions are normalized. If Ψ is an eigenstate of

A series of measurements on identical systems will give the averaged value

The expectation value is

The expectation value is therefore a weighted sum of the eigenvalues of Ω.

This probability interpretation of the wave function is due to Max Born.

Because the probability to find a particle at any place in space has to be

This is a partial differential equation. Ĥ is the Hamilton operator of the

the time-independent Schrödinger equation which determines the function

1.3 Hydrogen Atom

[ˆlz , ˆlx ] = i~ˆly

For a single component and l2 :

[ˆlx , ˆl2 ] = lˆx ˆl2 − ˆl2 ˆlx

[ˆlx , ˆl2 ] = [ˆly , ˆl2 ] = [ˆlz , ˆl2 ] = 0

From the commutators we learn that it is possible to measure the absolute

[ˆlx , T̂ ] = [ˆly , T̂ ] = [ˆlz , T̂ ] = 0

Angular momentum and Hamilton operators in spherical coordi-

Ψ(r, θ, φ) = R(r) Θ(θ) Φ(φ)

Eigenvalue problem for the operator ˆlz

ˆlz Φ(φ) = ~ ∂ Φ(φ) = lz Φ(φ)

Φ(φ) = Ceimφ constants C, m

Φ(φ) = Φ(φ + 2π) periodicity

Eigenvalue problem for the operator ˆl2

• Power series ansatz

• Investigation of convergence, i.e. are there divergencies

• Divergencies can only be avoided at special values of λ where we get a

Ylm (θ, φ) = Θ(θ)Φ(φ)

The θ dependent part of Θ(θ) is proportional to the Legendre polynomials

• Investigation of the asymptotic behavior for r → ∞ and r → 0 leads

• Solutions of this differential equation can be found with a power series

There is an infinite number of bound states with En < 0.

Excited states of hydrogen atom (Z = 1, En ∝ − n12 ):

E2 ≈ −3.4eV, E3 ≈ −1.5eV, E4 ≈ −0.85eV

Each state Ψnlm is characterized by three quantum numbers, however, the

1.3.3 Radial function of hydrogen atom

Ψ∗nlm Ψnlm dτ = Ψ∗nlm (r, θ, φ)Ψnlm (r, θ, φ)r 2 sin θdrdθdφ

Probability to find an electron within r and r + dr, independent of the angles

Functions of position and momentum