The interacting Hamiltonian for many body Schrödinger equation is (see the general QFT notes for derivation):

where are spin orbitals (thus the integration over below) and:

We would like to minimize the energy using the following basis for electrons:

We express the energy in this basis:

We minimize it with the constrain :

We obtain:

(1)

in the -representation:

And writing the individual terms explicitly (in this section, all orbitals are
*spin* orbitals):

we get the Hartree-Fock equations:

(2)

Let’s introduce the number density , Hartree potential and nonlocal exchange potential with its kernel :

then we can write the HF equations as:

The Hartree potential can be calculated by solving the Poisson equation:

where:

The application of the exchange potential on any function can be calculated by:

Starting from (1) and integrating over spins we get (here , are spatial orbitals, not spin orbitals):

(3)

We introduce basis functions by (below the greek letters are basis functions, latin letters are spatial orbitals):

substitute into (3) and multiply by from the left:

(4)

Now we expand the functions :

(5)

we introduce the density matrix:

and get:

(6)

introducing:

the equation (6) is:

(7)

These are the Roothaan equations. It is a generalized eigenvalue problem.

Total energy is given by (the in the first equation are spin orbitals, in the other equations are spatial orbitals):

The same thing can be derived in -representation starting from (2) and introducing spatial orbitals:

(8)

We introduce basis functions :

substitute into (8) and also multiply the whole equation by and integrate over :

(9)

This can be written as:

where:

Introducing the density matrix and density:

Expanding the functions and using the density matrix we get for :

or

In physical and chemistry notation this is written as:

Note that this notation implicitly assumes the factor, so for example actually means and one has to understand this from the context.

The two particle matrix element is:

(10)

The is called the physicists’ notation because the and kets are:

The is called the chemists’ notation. From (10) there are two types of symmetries — interchanging of the dummy variables:

and taking complex conjugate:

If the matrix elements are real, then:

In general those are the only symmetries (4 total).

If however, the functions are real, then there are additional symmetries: an exchange and . The symmetries of are exchange of with or with or the and pairs (8 total):

So if we view as two boxes then we can permute the labels in the given box “”, as well as exchange the boxes (the only thing we cannot do is to take one particle from one box and put it into the other). As such the box “” is a pair of two electrons (in any order) and the two electron integral assigns a unique number to a pair of such boxes (in any order). The symmetries of the symbol are:

Example I: the Slater integral has all 8 symmetries (Slater integral uses physical notation).

Example II: In spherical symmetry, the 2-particle integrals can be written as (see below for derivation):

They only have 4 symmetries, because spherical harmonics are complex. In particular:

and

We used the symmetries of the Slater integrals as well as the coefficients that change a sign, but thanks to the , the overall sign does not change. The other two symmetries are missing, i.e. and .

Spherical symmetry is this particular choice of a basis:

It can be shown that the solutions are of the form:

We can now write:

where

Also we get:

From which it follows:

The index runs over all combinations of . In particular, here is an example of one possible way to index the basis of 12 radial functions for each :

So the radial index always starts from 1 for each .

The G matrix is given by:

Now we use:

and

where

Where the sum over all occupied orbitals can be written as:

Where the sum over is the outer sum, both and depend on . We get:

The first part is the direct term, the second part the exchange term. Let’s treat the direct term first:

For the exchange term we get:

For this can be written as:

All together we get:

where

Note: performing the sum over and we get:

The Roothaan equations are:

The eigenvalues will be degenerate with respect to and so the radial Roothaan equations are to be solved for each :

The total energy is:

We use the following functions for :

And the multipole expansion:

And we get:

In the last step we used the fact that the symbols are zero unless and , from which it follows that and so one of the symbols is zero unless , which is expressed by . Given this condition, the sum over must be such that one is equal to , which means that otherwise the symbols will be zero. Finally, must also satisfy the conditions and . The sign factor is equal to both and so we just used the former.

We can write this using the symbols as:

We can also couple the angular momenta as follows:

and we get for the matrix elements:

Where we used the symbol:

Where we have renamed to .

In this section we express the matrix elements in the STO basis. It turns out that all integrals that we need can be expressed in terms of the following simple integral (where ):

(11)

The STO basis function for the radial Schrödinger equation for is:

(12)

Where the normalization constant is such that the STO orbital is normalized as the radial wavefunction :

from which we get:

Note that for we get the following STO basis function:

(13)

One uses either (12) or (13) depending on whether one solves the radial Schrödinger equation for or for .

In this section we also need the following integral:

where

is the incomplete gamma function.

The Slater integral is

where is the integral over the lower triangle (assuming is the -axis and is the -axis), that is :

where:

Much more tedious method is the following:

The Slater integral is given by:

where

and we get:

Putting everything together we get:

where:

In this section we express the matrix elements in the GTO basis. It turns out that all integrals that we need can be expressed in terms of the following simple integral (where ):

(14)

The GTO basis function for the radial Schrödinger equation for is:

(15)

However, unlike STO, the GTO functions must satisfy the condition where (this condition is later used in (14) to determine whether is even or odd). The normalization constant is such that the GTO orbital is normalized as the radial wavefunction :

from which we get:

Note that for we get the following GTO basis function:

(16)

One uses either (12) or (13) depending on whether one solves the radial Schrödinger equation for or for .

We will need the following integral:

Just like for STO, we get:

where is the integral over the lower triangle (assuming is the -axis and is the -axis), that is :

where and:

Let’s calculate the exchange integral

for the particular choice of the functions :

We use multipole expansion:

And we get:

We have a sum over electron states like this:

where are some functions that depend on the state numbers (for example squares of the wavefunctions). Then there are two options — either there is a way to sum over the and degrees of freedom, so that the sum can be written exactly as:

where (that don’t depend on and ) will in general be different to , but the sum will be the same. Or we have to approximate the sum (for example by averaging over the angles, or in some other way) as:

In either case, the occupation numbers are simply the number of times the functions appear in the sum for the given and :

So for closed shells atoms, it is always:

because there are two spins, and possibilities for , for open shell atoms, is anything between and .

As an example, let’s say that after some calculation for closed shell systems we get exactly:

Then because there are exactly states in the shell, we write the above as:

Then we do similar calculation for the open shell system, and we have to use some approximations to get the following formula, where the happen to be exactly the same as for the closed shell system:

Then we denote by the number of electrons in the shell (at least one of them will be open, for which we have ), and we can write the above as:

The usual chemical occupation numbers for the Uranium atom are:

So the , and , shells are open, all others are closed. By summing all these , we get 92 states as expected:

Code:

```
def f_nl(n, l):
if n < 5 or (n == 5 and l <= 2):
return 2*(2*l+1)
else:
d = {
(5, 3): 3,
(6, 0): 2,
(6, 1): 6,
(6, 2): 1,
(7, 0): 2,
}
if (n, l) in d:
return d[n, l]
else:
return 0
print "Sum f_nl =", sum([f_nl(n, l) for n in range(8) for l in range(n)])
```

prints:

`Sum f_nl = 92`

Hartree screening function is defined as:

and it occurs in many formulas in the Hartree Fock theory, so this section shows how to calculate it. It depends on and a function .

We first do the integral:

where:

Now we differentiate :

Also , so we get the following set of first order differential equations with boundary conditions:

One way to calculate the Hartree screening function is to integrate the second equation from the left using the boundary condition and then integrate the first equation from the right, using the boundary condition .

Another way is to obtain one second order equation. Expressing from the first equation:

and substituting into the second equation we get:

With boundary condition on the left:

and on the right:

which for becomes:

but in practise, it’s better to use the former Newton (Robin) boundary condition. We have obtained one second order equation for

with boundary conditions:

The weak formulation is:

The boundary term can be simplified using the boundary conditions as:

so we get

where the test functions have the constrain on the left boundary and no constrain on the right.

For both open and closed shell atoms we get exactly:

For closed shell atoms we use the fact, that

and the second term disappears, and for open shell atoms we have to use the central field approximation: we average the integral for over the angles:

and using the fact, that

the second term disappears as well. We got the same expression for both open shell (with central field approximation) and closed shell (no approximation) atoms. The radial charge density is:

So we got:

The Hartree screening function is given by the equation:

So satisfies the radial Poisson equation:

Similarly, we calculate:

Functions with different spins don’t contribute to the sum, so there is no
multiplication by 2. We assumed closed shells atoms (we summed over all in
the above). We used the result of the integral in
*Example VI* and also:

(17)

Using the above integrals, the HF equations become:

with:

Using the Hartree screening functions, the HF equations are:

with:

The total energy is:

where:

For Helium atom, the only nonzero occupation numbers are:

and the sum over simplifies to:

so we only need to solve for the state and we get:

with:

We can combine the equations:

and we obtain:

The weak formulation is ():

for closed shell atoms:

or (here we use the index to label all functions for the given )

Introducing radial basis (where , labels all basis functions for the given ):

we get (here is again restricted for the subset corresponding to the given ):

This can be written as

where

The indices go over all occupied orbitals . Introducing density:

We introduce the density matrix (where as before , run over basis functions for the given only):

where the coefficients are the same for all corresponding to the given . The index runs over all occupied states for the given . We can write as

Finally we get:

and

So we get:

The density matrix is zero if there are no occupied orbitals for the given .

The total energy is:

where we used:

and

The 4-index transformation is a way to convert the two particle integrals over basis functions into two particle integrals over atomic (or molecular) orbitals :

The self energy up to a second order is given by:

The , are occupied orbitals, , are virtual orbitals. The Dyson equation is:

(18)

where

The and are HF energies and eigenvectors, and are interacting energies and eigenvectors. The matrix is a diagonal matrix of the eigenvalues, is a diagonal matrix of the eigenvalues. Any Green’s function (interacting or not) can be written using the spectral density function as follows:

where

From (18) we get:

The poles of the Green’s function are then given by:

or equivalently:

and from the theory of matrices:

one obtains that

and

The number of particles can be calculated as follows ( are occupied orbitals, are all orbitals):

The countour only encloses poles corresponding to occupied orbitals . Similarly for the total energy :

For doubly filled orbitals we multiply the expressions by 2.