initial setup of thesis repository

[cluster_expansion_thesis.git] / General.Theory / DFT / DFT.Introduction.tex

\section{Quantum mechanics\cite{Adamson}}

It can be said that most of the properties of a chemical compound is 
determined by its electronic structure, which can be described
by the wavefunction, $\Psi$, of the 
system. Therefore access to the detailed electronic structure should 
enable prediction or evaluation of the behaviour of the system under 
investigation.

Energy is the fundamental quantity of a system. Despite this fact, 
molecular systems are usually described with other features, for 
example the 
molecular structure or the dipole moment. These 
properties are determined by external parameters, which may become 
variables on which a potential energy surface can be mapped. 
Therefore, analytical derivatives of the energy with respect to these 
variables yield the molecular properties. 
For example, the derivatives of the energy with respect to the nuclear 
positions give the nuclear forces, which allows minimisation of the 
energy with respect to nuclear coordinates, providing the molecular 
structure. The energy is thus the basis for the calculation of all 
molecular properties, and a method for evaluating this property is of 
fundamental interest.

Computational methods lead in different ways to the energy and 
electronic structure of the investigated compound. This introduction 
will describe Hartree-Fock methods and then move to Density Functional 
Theory (DFT), the method of choice in all calculations of this project.

\subsection{The Schrödinger equation; Born interpretation}

The Schrödinger equation\cite{Schroedinger} is able to describe all non-relativistic chemical 
phenomena. The time-dependent Schrödinger equation is
\[
\hat{H}\Psi=-i\hbar \frac{\partial \Psi}{\partial t},
\label{schroedinger-time}
\]
where $\Psi$ is the wavefunction of the molecular system and $\hat{H}$ is 
the Hamiltonian operator, which describes the 
energy of the system with the wavefunction. $\Psi$ is a 
function of the nuclear and electron positions, the electron spins and 
time. Positions $(\mathbf{r})$\footnote{Bold identifiers are 
vectors.} and spins $(s)$ are often combined in 
space-spin coordinates $(\mathbf{x})$. A system can be fully described with the wave function.

There are some different interpretations of $\Psi$, 
but the most common is due to Born\cite{Born}, who stated that:

\begin{quote}
   The probability of finding a system with the 
    wavefunction $\Psi$ in a given state is proportional to $\Psi^* 
    \cdot \Psi$. %\footnote{Wavefunctions with a superscript asterix are 
    %complex conjugate.}
\end{quote}


\subsection{Wavefunction properties}

The wavefunction is subject to some constraints.

As a result of the Born interpretation the integral of the system's exact wavefunction
over all space-spin 
coordinates must equal one, because the 
electron must be located somewhere. This is:
\[
\int \Psi^* \cdot \Psi \, \mathrm{d}\tau = 1. 
\label{normalized}
\]
Equation (\ref{normalized}) is called the {\em normalisation} criteria.

Wavefunctions are considered to be orthogonal, the {\em 
orthogonality} criteria. That means:
\[
\int \Psi^*_{m} \cdot \Psi_{n}\, \mathrm{d}\tau = 0 \Leftrightarrow m 
\not= n.
\label{orthogonal}
\]

Both properties are often summed up under {\em orthonormality} and 
expressed with the Kronecker delta:
\[
\int \Psi^*_{m} \cdot \Psi_{n} \, \mathrm{d}\tau = \delta_{mn} = 
\left\{ \begin{array}{lcr}
        0 & \Leftrightarrow & m \not= n \\
        1 & \Leftrightarrow & m = n
        \end{array} \right.
\label{Kronnecker}
\]

Because electrons are fermions, the wavefunction must be 
{\em antisymmetric} with respect to exchange of two electrons, i. e.
the wavefunction has to change sign if two electrons are 
exchanged. This is the Pauli Exclusion Principle.\cite{Pauli}

For most problems, a separation of the wave function into a 
time-dependent part and a time-independent part is possible. The 
wavefunction is then of the form:
\[
\Psi(t)=\Psi \cdot \mathrm{e}^{\left[\frac{-iEt}{\hbar}\right]}.
\]
Considering only the time independent part, the Schrödinger equation (\ref{schroedinger-time}) 
reduces to
\[
\hat{H}\Psi = E \Psi,
\label{schroedinger-time-independent}
\]
an eigenvalue equation, where $E$, the energy, is the eigenvalue for 
the operator $\hat{H}$ and 
$\Psi$ is the eigenfunction.

\subsection{The expectation value}

In general, considering $\hat{A}$ as a time-independent operator for a 
physical quantity and $A$ the resulting eigenvalue, the equation 
(\ref{schroedinger-time-independent}) changes to
\[
\hat{A}\Psi = A \Psi.
\label{equation1}
\]
Multiplying equation (\ref{equation1}) with the complex conjugate of 
$\Psi$, which is $\Psi^*$, and integrating over all coordinates leads 
to
\[
\int \Psi^* \hat{A} \Psi \, \mathrm{d}\tau = \int \Psi^* A \Psi \, \mathrm{d}\tau.
\]
$A$ is a constant scalar and so can be taken outside the integral
\[
\int \Psi^* \hat{A} \Psi \,  \mathrm{d}\tau =  A\int \Psi^* \Psi \,  \mathrm{d}\tau.
\] 
Dividing by $\int \Psi^* \Psi \,  \mathrm{d}\tau$ gives the {\em 
expectation 
value} of the physical quantity (which can be considered to be the 
average value of the physical quantity connected with $\hat{A}$ that 
can be observed):
\[
A=\frac{\int \Psi^* \hat{A} \Psi \, \mathrm{d}\tau}{\int \Psi^* \Psi 
\, \mathrm{d}\tau}.
\label{expecting-value}
\]
For normalised wavefunctions, $\int \Psi^* \Psi \, \mathrm{d}\tau=1$, 
 and so equation (\ref{expecting-value}) simplifies to $A=\int \Psi^* \hat{A} \Psi \, \mathrm{d}\tau$ or 
written in {\em Bra-Ket} notation;
\[
A= \langle \Psi | \hat{A} | \Psi \rangle =\int \Psi^* \hat{A} \Psi \, \mathrm{d}\tau 
\]
One can say $A$ is a {\em functional} of $\Psi$, written as 
$A[\Psi]$, where $\Psi$ is itself a function of space-spin coordinates.


\subsection{The Hamiltonian operator}

If we consider a system of $M$ nuclei of charge $Z_{A}$ and $N$ 
electrons without an external field, the Hamiltonian, in atomic 
units, is given by
\[
\hat{H}=\hat{T}+\hat{V}
\]
where,
\[
\hat{T}= \underbrace{-\frac{1}{2} \sum\limits_{A}^{M} 
\nabla^2_{A}}_{\mbox{\tiny{for the nuclei}}} \underbrace{ - 
\frac{1}{2} \sum\limits_{i}^{N} \nabla^2_{i}}_{\mbox{\tiny{for the 
electrons}}}
\]
for the kinetic energy of nuclei and electrons, and
\[
\begin{array}{lll}
\hat{V}& = & \underbrace{-\sum\limits_{A}^{M} \sum\limits_{i}^{N} 
\frac{Z_{A}}{\left|\mathbf{R}_{A}-\mathbf{r}_{i}\right|}}_{\mbox{\tiny{nuclear-electron (ne) attraction}}}\\
&&\\
& & \underbrace{+\sum\limits_{A}^{M} \sum\limits_{B>A}^{M} 
\frac{Z_{A}Z_{B}}{\left|\mathbf{R}_{A}-\mathbf{R}_{B}\right|}}_{\mbox{\tiny{nuclear-nuclear (nn) repulsion}}}\\
&&\\
& & \underbrace{-\sum\limits_{i}^{N} \sum\limits_{j>i}^{N} 
\frac{1}{\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|}}_{\mbox{\tiny{electron-electron (ee) repulsion}}}
\end{array}
\]
for the potential energy of attractive and repulsive interactions of 
nuclei and electrons.

Thus the total energy of the system, $E_{\mbox{\tiny{tot}}}$, is 
 the sum of all partial energies.

\subsection{Born-Oppenheimer approximation}
It is common to treat quantum mechanical problems from the point of 
view of the Born-Oppenheimer approximation.\cite{Born-Oppenheimer}

The masses of the nuclei are much greater than the masses of the 
electrons, hence the electrons can respond almost instantaneously to 
any change in the nuclear positions. Thus, to a good approximation, we 
can think of the electrons moving in a field of fixed nuclei. As a 
consequence, the terms introducing the nuclear kinetic energy can neglected and 
the nuclear-nuclear repulsion can then be considered as constant. Removing these terms leads to the electronic Hamiltonian,
the electronic wavefunction and the electronic Schrödinger equation. 
Only these will be considered in the following treatment.

The resulting energy will be the electronic energy. In order to regain the 
total energy, the nuclear-nuclear repulsion terms must be added. 

This reduces the Hamilatonian to:
\[
\hat{H}=-\frac{1}{2} \sum\limits_{i}^{N} \nabla^2_{i}-\sum\limits_{A}^{M} \sum\limits_{i}^{N} 
\frac{Z_{A}}{\left|\mathbf{R}_{A}-\mathbf{r}_{i}\right|} + -\sum\limits_{i}^{N} \sum\limits_{j>i}^{N} 
\frac{1}{\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|}
\]
Introducing the external potential $v(\mathbf{r}_{i})$
\[
v(\mathbf{r}_{i})= - \sum\limits_{A}^{N}\frac{Z_{A}}{\left|\mathbf{R}_{A}-\mathbf{r}_{i}\right|}
\]
and the shorthand form
\[
r_{ij}= \sum\limits_{i}^{N}\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|
\]
follows
\[
\hat{H}={-\frac{1}{2} \sum\limits_{A}^{M}\nabla^2_{A}} + 
\sum\limits_{i}^{N} v(\mathbf{r}_{i}) + 
\sum\limits_{j>i}^{N}\frac{1}{r_{ij}}.
\label{Hamiltonian-simplyfied}
\]


\subsection{The Variation Method}

As stated above, the expectation value of an operator is given 
by equation (\ref{expecting-value}) and the operator for the molecular 
energy is the Hamiltonian, $\hat{H}$. Consider an approximate
wavefunction $\Psi_{\mbox{\tiny{app}}}$ and the associated expectation 
value of the energy 
$E_{\mbox{\tiny{app}}}$, which can be obtained with
\[
E_{\mbox{\tiny{app}}}=\frac{\int \Psi^*_{\mbox{\tiny{app}}} \hat{H} 
\Psi^{}_{\mbox{\tiny{app}}}\, \mathrm{d}\tau}
{\int \Psi^*_{\mbox{\tiny{app}}} \Psi^{}_{\mbox{\tiny{app}}}\, \mathrm{d}\tau}.
\]
The varation theorem states that the energy calculated from a trial $\Psi_{\mbox{\tiny{app}}}$ is an upper bond to the 
true ground state energy. Minimisation 
with respect to all allowed $N$-electron wavefunctions will give the 
true ground state energy and the exact wavefunction, which is minimization of the functional 
$E_{\mbox{\tiny{app}}}[\Psi_{\mbox{\tiny{app}}}]$, known as 
the {\em Rayleigh-Ritz} minimal principle. 
Therefore
\[
E_{\mbox{\tiny{app}}} = E_{0} \Leftrightarrow 
\Psi_{\mbox{\tiny{app}}} = \Psi_{0} \qquad .
\]
This minimisation with respect to all allowed $\Psi$s is not practical. What is usually
done is to develop the trial  
wavefunction as a linear combination of well-defined wavefunctions of 
component parts of the system.

If $\Psi_{\mbox{\tiny{app}}}$ is developed 
as a linear combination of wavefunctions,
\[
\Psi_{\mbox{\tiny{app, i}}}=\sum\limits_{i} c_{i}\Psi_{i},
\]
$E_{i}$ is the energy for the $i$th eigenstate of $\hat{H}$ resulting in the combination 
of the $i$ wavefunctions. As stated above, every trial wavefunction yields an energy 
greater than or equal the true ground state energy. As a resuslt, the process of energy minimisation 
is reduced to a process of finding an optimum set of coefficients with this ansatz. This 
task can be solved via matrix diagonalisation.
%
%$E_{0}\ge E_{1} \ge E_{2} \ge \cdots E[\Psi]$, is always greater than 
%or equal to $E_{0}$. Thus the 
%problem is  reduced, one has to find the optimum sets of coefficients, $c_{i}$, which 
%can be achieved via matrix diagonalisation.



% this last paragraph has to be reviewed once again due to some 
% unclearness in its meaning and purpose


\subsection{Molecular wavefunction, atomic orbitals and basis sets}

The wavefunction $\Psi$, which is used here, is the molecular 
wavefuntion, describing the complete electronic structure of a given 
molecule. For a molecular wavefunction $\Psi_{\mbox{\tiny{app}}}$, the 
most obvious choice of component wavefunctions,  
$\Psi_{\mbox{\tiny{i}}}$, are those of the atoms that make up the 
molecule. 

As a consequence, the molecular wavefunction, $\Psi$, can be conceived as a combination 
of one-electron wavefunctions, $\chi_{i}(\mathbf{x})$, which may be 
regarded as atomic orbitals. This is called the molecular orbital 
approximation.

The $\chi_{i}$ consists of a spatial 
orbital $\psi_{i}(\mathbf{r})$ and one of the two spin functions 
$\alpha(s)$, spin-up, or $\beta(s)$, spin-down, thus we get
$\chi_{i}(\mathbf{x})=\psi_{i}(\mathbf{r})\alpha(s_{i})$ or 
$\chi_{i}(\mathbf{x})=\psi_{i}(\mathbf{r})\beta(s_{i})$.

%The orbital approximation states that you can consider all electrons 
%as occupying 1-electron orbitals similar to those of 

%in this paragraph we have to introduce the concepts of:
%- one electron molecular wave funtions , so calles orbitals,
%- separation of these functions into on spatial and one spin 
%dependent part
%- introduce the concept of MO theory..
%.. which in my opinion has nothing to do with atomic wave functions, 
%depending of the atoms building up the molecule!!


\subsubsection{Basis sets\cite{Computational-Chemistry}}

Because the exact mathematical description of the spatial orbitals 
are even almost unknown for more complex molecules than hydrogen, they are usually developed in a set of 
known functions, the {\em basis set}.
\[
\psi_{i}=\sum\limits_{\mu}c_{\mu i}\phi_{\mu}
\]
This is not an approximation as long as the basis set is complete, i. e.
it contains an infinite number of functions. This is not possible in 
practice, so the 
introduction of the basis set necessarily introduces some error.

The most convenient way to define a basis set for any nuclear 
configuration is to define a particular set of functions for each
nucleus, depending only its charge.

In principle, any kind of function, or combination of functions can 
be used. Their behaviour should agree with the physics of the 
problem, e. g.
\begin{itemize}
    \item{For bound atomic and molecular systems the function should 
    tend towards zero as the distance between nuclei and electron 
    becomes large.}
    \item{The energy should converge reasonably rapidly as more basis 
    functions are added.}
\end{itemize}

The best compromise between the flexibility of the basis set, which 
increases with the number of basis functions, and the computational 
cost has to be made. The basis set should be kept as small as possible 
to minimize the computational cost.

There are two main types in use today. The Slater-Type Atomic Orbitals (STOs) introduced in 
1930 by Slater\cite{Slater-Orbitals} and the Gaussian-Type Atomic 
Orbitals (GTOs), introduced in 1950 by Boys\cite{Gaussian-Orbitals}.

\paragraph{STOs} are of the form
\[
\phi_{a}(\mathbf{r})=(x-A_{x})^{a_{x}}(y-A_{y})^{a_{y}}(z-A_{z})^{a_{z}}\cdot \mathrm{e}^{-\alpha |\mathbf{r}-\mathbf{A}|}
\]
with a center $\mathbf{A}(A_{x},A_{y},A_{z})$, angular momentum 
$\mathbf{a}(a_{x},a_{y},a_{z})$, a nuclei dependent exponent $\alpha$ 
and exponential radial parts, which takes care of the decay for large 
distances. Like the real wavefunctions they have cusps at the nuclei. 
Their major disadvantage is that integrals over STOs are expensive to 
compute.

To overcome this difficulty, Boys suggested, in 1950, the construction 
of basis functions with Gaussian-Type Atomic 
Orbitals, GTOs.

\paragraph{GTOs} are of the form
\[
\phi_{a}(\mathbf{r})=(x-A_{x})^{a_{x}}(y-A_{y})^{a_{y}}(z-A_{z})^{a_{z}}\cdot \mathrm{e}^{-\alpha |\mathbf{r}-\mathbf{A}|^2}.
\]
GTOs decay too fast and have incorrect nuclear cusps, but the speed 
with which integrals over GTOs can be computed compensates for this. 
As a result, if the same accuracy as with STOs is required, they can be 
approximated by a sum of Gaussian-Type orbitals. This leads to the
STO-$n$G basis sets, where one STO is represented by so-called 
{\em contraction} of $n$ Gaussian Type Orbitals. A contraction is a linear 
combination of Gaussian Orbitals where the coefficients are not 
changed during the calculation.

Having specified the type of functions to be used and the locations, 
it must be decided what number of functions are to be used. The 
smallest number of functions possible is a {\em minimal basis set}, 
where only one function for each electron of the neutral atoms are 
used. Improvements in accuracy and flexibility can be achieved through 
using double the amount of basis functions needed, a 
{\em Double Zeta} basis set, indicated by double-$\zeta$. The next 
step is the use of three times more functions than the number of 
electrons, triple-$\zeta$. A variation of the basis set which only increases 
the number of valence orbitals produces a {\em split valence basis set}.

This is reasonable, because core electrons are not significantly
influenced by the bonding environment, especially for atoms with a 
large number of electrons. Functions with a higher angular momentum 
are denoted {\em polarisation functions}. 
%Usually the next orbital function of the last occupied valence orbital 
%is used to improve the 
%description of the electronic distribution in a certain direction, e.g. p-functions to polarize s-functions, 
%f-functions to polarize d-functions, etc.

If basis sets of the form {\em k-nlm}G{\em w} are referred to, then some 
special kinds of split-valence basis sets, designed by {\em Pople et. 
al}\cite{Pople-basis-sets1, Pople-basis-sets2, Pople-basis-sets3, 
Pople-basis-sets4} have been used. The {\em k} indicates how many 
primitive GTOs have been contracted to represent the core orbitals, 
the number of values after the dash indicates the grade of splitting 
of the valence orbitals, {\em nlm} is thus a triple split valence 
basis set. Thus the valence basis is represented with three functions, 
each of which is constructed from the number of Gaussian Type Orbitals given.
Values after the G indicates polarisation functions.



\subsection{The Hartree approximation, Hartree-Fock Theory, Slater 
Determinant}

If you want to construct a molecular wavefunction built up from 
one-electron wavefuntions a very simple approach is to construct your 
molecular wavefuntion as a combination of {\em independent} 
one-electron wavefunctions,

The simplest way to construct a molecular wavefunction of one-electron 
wavefunction is the independent-particle, or Hartree-approximation, 
where the total wavefunction is built up as a product of orthonormal 
one-electron wavefunctions\cite{Hartree1, Hartree2, Hartree3}. This 
approximation assumes that each 
electron moves independently in its own orbital, only influenced by 
the {\em average} field generated by all the other electrons and explicit electron 
interaction is neglected.

The 
Hartree wavefunction for a system of $N$ electrons is
\[
\Psi = \chi_{1}(\mathbf{x}_{1})\chi_{2}(\mathbf{x}_{2}) \ldots 
\chi_{N-1}(\mathbf{x}_{N-1})\chi_{N}(\mathbf{x}_{N}).
\]

This wavefunction is an invalid one, because it is not antisymmetric with 
respect to electron interchange. 

This was first mentioned by Fock in 
1930\cite{Fock} and one solution he suggested 
was to add and substract all possible permutations of the Hartree 
product wavefunction, forming the Hartree-Fock (HF) wavefunction. 
Slater found later that the result is the determinant of a matrix, 
called the Slater determinant\cite{Slater1, Slater2}
\[
\Psi_{\mbox{\tiny{SD}}} = \frac{1}{\sqrt{N!}} \left|
\begin{array}{cccc}
  \chi_{1}(\mathbf{x}_{1}) & \chi_{2}(\mathbf{x}_{1}) & \cdots & 
  \chi_{N}(\mathbf{x}_{1}) \\
  \chi_{1}(\mathbf{x}_{2}) &  \chi_{2}(\mathbf{x}_{2}) & \cdots & 
  \chi_{N}(\mathbf{x}_{2})\\
  \vdots & \vdots & & \vdots \\
  \chi_{1}(\mathbf{x}_{2})& \chi_{1}(\mathbf{x}_{N})& \cdots &  
  \chi_{N}(\mathbf{x}_{N})
\end{array} \right|.
\]
The prefactor $1/\sqrt{N!}$ normalizes the HF wavefunction.

Usually the trial wavefunction is considered as 
a {\em single} Slater determinant, which implies that 
electron-interaction is neglected (the electron-electron repulsion is 
only included as an average effect). If additional approximations are 
made, this leads to {\em semi-empirical} methods. If additional 
determinants are included, this leads to more exact {\em post-HF} methods.

\subsubsection{Hartree-Fock equations}

Hartree-Fock calculations start with an initial guess of the 
wavefunction and then successively improve this wavefunction with 
respect to the energy following the variational principle and 
optimizing the set of coefficients.


Considering the Born-Oppenheimer approximation the Hamiltonian is 
given as
\[
\hat{H}=\hat{T}_{\mbox{\tiny{e}}} + \hat{V}_{\mbox{\tiny{ne}}} + 
\hat{V}_{\mbox{\tiny{ee}}}\ ,
\]
grouping according to the number of electron indices (e) leads to
\begin{eqnarray}
\hat{h}_{i}&=& \hat{T}_{\mbox{\tiny{e}}}+\hat{V}_{\mbox{\tiny{ne}}} \nonumber \\
\hat{g}_{ij}&=& \hat{V}_{\mbox{\tiny{ee}}}\ .
\end{eqnarray}

The one-electron operator $\hat{h}_{i}$ describes the motion of 
electron $i$ in the field of all nuclei and the electron-electron 
repulsion is described with the two-electron operator $\hat{g}_{ij}$.

To evaluate the energy of a single Slater determinant, we apply these 
operators to the determinant. For each of the operators we only get 
 a non-zero value if there is no contribution of overlap integrals, i. e. 
for the one-electron operator, if both one electron wavefunctions of 
the determinant are of the same electron. That yields for the 
expectation value
\begin{eqnarray}
 & {} & \bra \chi_{i}(\mathbf{x}_{i}) \chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) |\hat{h}_{i}| \chi_{i}(\mathbf{x}_{i})
\chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) \ket  \nonumber \\
 & = & \bra \chi_{i}(\mathbf{x}_{i})|\hat{h}_{i}|\chi_{i}(\mathbf{x}_{i}) \ket 
\bra \chi_{i+1}(\mathbf{x}_{i+1})|\chi_{i+1}(\mathbf{x}_{i+1}) \ket \ldots
\bra \chi_{N}(\mathbf{x}_{N})|\chi_{N}(\mathbf{x}_{N}) \ket \\ 
 & = & \sum\limits_{i}^{N} h_{i} \ . \nonumber
\end{eqnarray}
For the two electron operator, the value is non-zero only if 
the same orbital combinations contribute and every possible permutation 
in which only two orbitals contribute. For 
the first case, if identical orbitals contribute, it leads to
\begin{eqnarray}
 & {} & \bra \chi_{i}(\mathbf{x}_{i}) \bra \chi_{j}(\mathbf{x}_{j}) \chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) |\hat{g}_{ij}| \chi_{i}(\mathbf{x}_{i}) 
\chi_{j}(\mathbf{x}_{j})
\chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) \ket \nonumber \\
 & = & \bra \chi_{i}(\mathbf{x}_{i}) \bra 
 \chi_{j}(\mathbf{x}_{j})|\hat{g}_{ij}|\chi_{i}(\mathbf{x}_{i}) \bra \chi_{j}(\mathbf{x}_{j})\ket 
\bra \chi_{i+1}(\mathbf{x}_{i+1})|\chi_{i+1}(\mathbf{x}_{i+1}) \ket \ldots
\bra \chi_{N}(\mathbf{x}_{N})|\chi_{N}(\mathbf{x}_{N}) \ket \nonumber \\ 
 & = & \sum\limits_{i}^{N} \sum\limits_{j> i}^{N} J_{ij}\ . 
\end{eqnarray}
This is called {\em Coulomb} integral. It represents a classical 
repulsion between two charge distributions described by the $i$th and 
$j$th wavefunction.

For the second case the integral is
\begin{eqnarray}
 & {} & \bra \chi_{i}(\mathbf{x}_{i}) \bra \chi_{j}(\mathbf{x}_{j}) \chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) |\hat{g}_{ij}| \chi_{j}(\mathbf{x}_{j})  
\chi_{i}(\mathbf{x}_{i})
\chi_{i+1}(\mathbf{x}_{i+1}) \ldots 
\chi_{N}(\mathbf{x}_{N}) \ket \nonumber \\
 & = & \bra \chi_{i}(\mathbf{x}_{i}) \bra 
 \chi_{j}(\mathbf{x}_{j})|\hat{g}_{ij}|\chi_{j}(\mathbf{x}_{j}) \bra 
 \chi_{i}(\mathbf{x}_{i})\ket 
\bra \chi_{i+1}(\mathbf{x}_{i+1})|\chi_{i+1}(\mathbf{x}_{i+1}) \ket \ldots
\bra \chi_{N}(\mathbf{x}_{N})|\chi_{N}(\mathbf{x}_{N}) \ket \nonumber \\ 
 & = & \sum\limits_{i}^{N} \sum\limits_{j> i}^{N} K_{ij} \ . 
\end{eqnarray}
It is termed {\em exchange} integral and has no classical analogue.

The energy of a single Slater determinant can then be written as
\begin{eqnarray}
E_{\mbox{\tiny{SD}}} &=&  \sum\limits_{i}^{N} h_{i} + 
\sum\limits_{i}^{N} \sum\limits_{j> i}^{N} J_{ij} -
\sum\limits_{i}^{N} \sum\limits_{j> i}^{N} K_{ij}  \nonumber \\
&=& \sum\limits_{i}^{N} h_{i} + \sum\limits_{i}^{N} \sum\limits_{j\ge 
i}^{N} \left( J_{ij} -  K_{ij} \right) \\
&=& \sum\limits_{i}^{N} h_{i} + \frac{1}{2}\sum\limits_{i}^{N} 
\sum\limits_{j}^{N} \left( J_{ij} -  K_{ij} \right) \ .  \nonumber
\end{eqnarray}

One remarkable point is that for $i=j$, the Coulomb operator 
generates a spurious self-interaction of the electron, which is exactly 
compensated by the coresponding exchange operator, because
$J_{ii}=K_{ii}$. 

For Hartree-Fock calculation, the  energy is expressed in terms of operators, the identifiers are the 
same as for the integrals. Minimizing the energy with respect to the 
orbitals is then carried out following the variation principle. This 
must be carried out in a way that the orbitals remain orthonormal. This 
constraint optimization is usually done self-consistently with the method of {\em 
Lagrange} multipliers, as small changes in the orbitals should not 
change the Lagrange function. The {\em Fock} operator is then 
defined, which is an effective one-electron operator and associated 
with the variation of the energy. The result is a set of orbitals satisfying the {\em 
Hartree-Fock equations}. A normal Hartree-Fock 
calculation starts with a guess of the set of orbitals, and then 
iteratively minimises the energy with respect to the orbitals. 
Obviously a good starting guess is necessary to gain a fast convergence.

The principle of {\em self-consistent} calculations is that based on 
the results of the first set, a second set is generated, then a third 
one based on the results of the second one and so on. This is carried 
on until the energy variation is beneath a predefined value and then 
is called {\em self-consistent} or stationary. Subsequent iterations would give a 
variation of the energy under this threshold. The {\em Hartree-Fock} 
limit is the lowest energy can be reached 
with HF calculations in the limit of the chosen basis set.

Explicit 
electron-electron interaction is neglected. Electrons of the same 
spin are partially but not completely, correlated. No attemp is made 
to include electron-correlation of opposite-spin electrons.

 Because of this neglect, the energy of HF calculations, 
$E_{\mbox{\tiny{HF}}}$ is  always above the exact non-relativistic 
value, $E_{\mbox{\tiny{exact}}}$, the difference is defined as the 
{\em correlation energy}, 
$E_{\mbox{\tiny{corr}}}=E_{\mbox{\tiny{exact}}}-E_{\mbox{\tiny{HF}}}$.\cite{Loewdin}
This is the main deficiency of HF theory.

\subsubsection{Post-HF methods}

The major aim of post-HF methods is to include electron correlation. 
There are different approaches to describe the interactions of the 
electrons correctly,
\begin{itemize}
    \item{via Configuration Interaction (CI), in which a molecular 
    wavefunction is constructed as a linear combination of wavefunctions 
    of the HF determinant and determinants corresponiding to 
    excitation of electron. 
    If the expansion of the molecular wavefunction as a linear 
    combination is complete, every possible excited state is 
    included, the exact correlation energy within the 
    basis set approximation is achieved.}
    
%    \item{via Coupled Cluster (CC) Theory, in which the coefficients 
%   of the additional 
%    wavefunctions are constructed of by the exponentiated operator of 
%    the excitations. In 
%    this way, higher excitations are only partially included and their 
%    coefficients are determined by lower order excitations. This 
%    reveals typically $95\%$ of the correlation energy.}
    \item{via M{\o}ller-Plesset Pertubation Theory, in which the 
    exact Hamiltonian is treated as a {\em small} perturbation of the 
    HF-Hamiltonian.}
\end{itemize}

All these approaches are based on wavefunction-expansions of the many 
electron wavefunction $\Psi(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})$.  This is 
combined with an increase in the number of integrals to solve 
according to the $3N$-coordinate dependency of the wavefunction. If 
spin is included, this becomes $4N$.  


\section{Density Functional Theory\cite{Parr-Yang}}

\subsection{The electron density}

For a given state with wave function $\Psi$ the correlated electron density 
is given by
\[
\rho(\mathbf{r})=N \int \Psi^{*} \cdot \Psi \, \mathrm{d}s_{1} 
\mathrm{d}\mathbf{x}_{2}\mathrm{d}\mathbf{x}_{3}\ldots\mathrm{d}\mathbf{x}_{N}
\label{electron-density}
\]
which is the square sum over all occupied orbitals
\[
\rho(\mathbf{r})= \sum\limits_{i}^{\mbox{\tiny{occ.}}} 
\psi_{i}(\mathbf{r})\cdot \psi_{i}^{*}(\mathbf{r})
\]
and it is
\[
\int \rho(\mathbf{r}) \, \mathrm{d}\mathbf{r} = N \ .
\]

If the exact electron density is known, then the cusps in 
$\rho(\mathbf{r})$ will provide the position of the nuclei. The slope 
of $\rho(\mathbf{r})$ at the nucleus $A$ must obey
\[
\left. \frac{\partial}{\partial 
\mathbf{r}_{A}}\bar{\rho}(\mathbf{r}_{A}) \right|_{\mathbf{r}_{A}} = 
-2Z_{A}\cdot \bar{\rho}(0)
\]
giving the charge at the nucleus, $Z_{A}$ ($\bar{\rho}$ denotes the 
spherical average of the density). As the Hamiltonian is defined 
completely by the nuclear charges and positions, it is fully known. 
Therefore the wavefunction and energy can be found and the system can be 
completely described by the electron density. A Hamiltonian like
equation (\ref{Hamiltonian-simplyfied}) is completely determined by 
the external potential $v(\mathrm{r})$.


\subsection{The first Hohenberg-Kohn theorem\cite{HK1}}
The theoream says:

For non-degenerate ground-states
\begin{quote}
    the external potential $v(\mathrm{r})$  is
    determined by the electron density, $\rho(\mathbf{r})$.
\end{quote}

More precisely, there is a one-to-one correspondence between the electron 
density of a system and its energy and vice versa. This theorem can be proven with the 
variational principle.\cite{Kohn-Rev.Mod.Phys.-1999}

%Consider a system with the potential $v_{a}(\mathbf{r})$ and the corresponding 
%electron density $\rho$. Consider another system with the same 
%electron density $\rho$ and a different potential $v_{b}(\mathbf{r})$. 
%A Hamiltonian $\hat{H}_{a}$ and a wavefunction 
%$\Psi_{a}$ are connected 
%with potential $a$. Because potential $b$ is different from potential $a$, there 
%must be a different wavefunction $\Psi_{b}$ with the Hamiltonian 
%$\hat{H}_{b}$. Using the variational principle,
%\begin{eqnarray}
%E_{a}^0 < \bra \Psi_{b}|\hat{H}_{a}|\Psi_{b} \ket & = & \bra 
%\Psi_{b}|\hat{H}_{b}|\Psi_{b} \ket + \bra 
%\Psi_{b}|\hat{H}_{a}-\hat{H}_{b}|\Psi_{b} \ket \nonumber \\
%&=& E_{b}^0 + \int \rho(\mathbf{r}) 
%[v_{a}(\mathbf{r})-v_{b}(\mathbf{r})] \, \mathrm{d}\mathbf{r}.
%\end{eqnarray}
%Similarly
%\[
%E_{b}^0= E_{a}^0 + \int \rho(\mathbf{r}) 
%[v_{b}(\mathbf{r})-v_{a}(\mathbf{r})] \, \mathrm{d}\mathbf{r},
%\]
%which leads to the contradiction
%\[
%E_{a}^0 + E_{b}^0 < E_{b}^0 + E_{a}^0 \ .
%\]

Hence the external potential is determined by the density and we can 
represent the energy as a functional of the density. In comparison 
with the wave mechanics approach, the energy may be divided into 
several parts: kinetic energy, $T[\rho]$ and attraction between the 
nuclei and electrons, $V_{\mbox{\tiny{ne}}}[\rho]$.

The last may be written in 
terms of the external potential and the electron density as
\[
V_{\mbox{\tiny{ne}}}[\rho] = \int \rho(\mathbf{r})v(\mathbf{r})\ \mathrm{d}\mathbf{r}.
\]
Further, the electron-electron repulsion -- divided into Coloumb part, $J[\rho]$, 
and an exchange part, $K[\rho]$:
\[
E[\rho]= T[\rho] + V_{\mbox{\tiny{ne}}}[\rho] + J[\rho] + K[\rho]\ .
\]

The theorem has since been extended to non-degenerate ground 
states.\cite{Dreizler}

\subsection{The second Hohenberg-Kohn theorem\cite{HK1}}

The second Hohenberg-Kohn theorem adopts the variational principle 
into density functional theory.
\begin{quote}
    For a trial density $\tilde{\rho}(\mathbf{r})$, such that 
    $\tilde{\rho}(\mathbf{r})\ge 0$ and $ \int 
    \tilde{\rho}(\mathbf{r}) \, \mathrm{d}\mathbf{r} = N$, an energy  
    $E[\tilde{\rho}] \ge E^0$ arises.
\end{quote}
This implies, an incorrect density will yield a energy which is higher than 
the exact ground state energy for every trial density supplied by a 
trial wavefunction.


%Assume a trial wavefunction $\Psi_{\rho_{0}}$, that integrates to the 
%ground state densisty $\rho_{0}(\mathbf{r})$, and the real ground state 
%wavefunction $\Psi_{0}$. The varational principle gives
%\[
%\bra \Psi_{\rho_{0}} | \hat{H} | \Psi_{\rho_{0}} \ket \ge \bra 
%\Psi_{0} | \hat{H} | \Psi_{0} \ket = E^0 \ . 
%\]
%Splitting the Hamiltonian into the usual parts and taking into account 
%that the potential energy due to the external potential is a 
%functional of the density leads to
%\begin{eqnarray}
%\bra \Psi_{\rho_{0}} | \hat{T} + \hat{J} + \hat{K} | \Psi_{\rho_{0}} 
%\ket + \int \rho_{0}(\mathbf{r})v(\mathbf{r})\, 
%\mathrm{d}\mathbf{r}  &\ge&  \bra \Psi_{\rho_{0}} | \hat{T} + \hat{J} + \hat{K} | \Psi_{\rho_{0}} 
%\ket + \int \rho_{0}(\mathbf{r})v(\mathbf{r})\, 
%\mathrm{d}\mathbf{r} \nonumber \\
%\bra \Psi_{\rho_{0}} | \hat{T} + \hat{J} + \hat{K} | \Psi_{\rho_{0}} 
%\ket & \ge & \bra \Psi_{\rho_{0}} | \hat{T} + \hat{J} + \hat{K} | \Psi_{\rho_{0}} 
%\ket \ .
%\end{eqnarray}
%Thus the $\Psi_{0}$ is the wavefunction that integrates to $\rho_{0}$ 
%and minimizes the expectation value of $\hat{T} + \hat{J} + \hat{K}$. 
%Defining a universal functional as
%\[
%F[\rho]=\min_{\scriptscriptstyle{\Psi \to \rho}} \bra \Psi | 
%\hat{T} + \hat{J} + \hat{K} | \Psi \ket \ ,
%\]
%which searches all $\Psi$'s that yields the input density $\rho$, 
%allows the energy to be expressed as
%\begin{eqnarray}
%E^0 & = & \min_{\scriptscriptstyle{\rho}} \left[ F[\rho] +  \int \rho(\mathbf{r})v(\mathbf{r})\, 
%\mathrm{d}\mathbf{r} \right] \nonumber \\
% &=& \min_{\scriptscriptstyle{\rho}} \ E[\rho] \ ,
%\end{eqnarray}
%which is a search over all $N$-representable densities. $N$-representable 
%means that the density must be obtainable from some antisymmetric 
%wavefunction.

The task is now to find an accurate functional $E[\rho]$ and minimize it with respect to 
the electron densitiy and the corresponding energy.


\subsection{The functionals\cite{Computational-Chemistry}}

The energy functional can be split according to the wave mechanics 
approach into the usual parts
\[
E[\rho]= T[\rho] + V_{\mbox{\tiny{ne}}}[\rho] + J[\rho] + K[\rho]\ .
\]
The functionals for the electron-nucleus attraction and for the 
Coulomb electron interaction can be expressed in their classical form
\begin{eqnarray}
 V_{\mbox{\tiny{ne}}}[\rho]& = & \int \rho(\mathbf{r})v(\mathbf{r})\, 
\mathrm{d}\mathbf{r}   \\
J[\rho] & = & \frac{1}{2} \int \!\!\! \int 
\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\, 
\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}'
\end{eqnarray}

There are no known ways to express the functionals for the 
kinetic energy and the exchange energy classically.



\subsubsection{Thomas-Fermi-Dirac functionals}

An early approach to obtain approximate expressions for the unknown 
functionals are by {\em Thomas}\cite{Thomas}, {\em Fermi}\cite{Fermi1, Fermi2, 
Fermi3}, {\em Block}\cite{Block} and {\em  Dirac}\cite{Dirac} considering 
a non-interacting electron gas. For this model the functionals are 
given as
\begin{eqnarray}
    T[\rho] &=& \frac{3}{10}(3\pi^2)^{2/3} \int 
    \rho^{5/3}(\mathbf{r}) \, \mathrm{d}\mathbf{r} \\
    K[\rho] &=& \frac{3}{4}\left( \frac{3}{\pi} \right)^{1/3} \int 
    \rho^{4/3}(\mathbf{r}) \, \mathrm{d}\mathbf{r} 
\end{eqnarray}

The assumption of a non-interacting electron gas is not good for 
atomic and molecular systems. It may be a usable model for solid state 
calculations, and may give reasonable results there.

For molecular systems, for example, this model does not predict any 
bonding, which is obviously not reasonable and can therefore not 
compete with the results of wave mechanic approaches.

The main problem in Thomas-Fermi models is that the kinetic energy 
part is poorly represented.


\subsubsection{Kohn-Sham Theory}

The foundation for the use of DFT in computational chemistry was the 
introduction of orbitals by {\em Kohn} and {\em 
Sham}.\cite{Kohn-Sham} With this introduction, the treatment of the main problem in 
Thomas-Fermi Theory can be significantly improved. The basic idea is 
to minimize the error in the kinetic energy term by splitting the 
kinetic energy into two parts, one of which can be calculated exactly, 
and a small correction term.

Consider a system of {\em non-interacting} electrons. For this system, 
the kinetic energy is described exactly with the functional of the
Hartree-Fock-Slater theory, that is
\begin{eqnarray}
T_{\mbox{\tiny{S}}} & = & \sum \limits_{i}^N \bra \chi_{i} | \hat{T} | 
\chi_{i} \ket \nonumber \\
 & = & \sum \limits_{i}^N \bra \chi_{i} | -\frac{1}{2} \nabla^2 | 
\chi_{i} \ket \ .
\end{eqnarray}
The electron density for this system is given exactly by
\[
\rho(\mathbf{r})=\sum\limits_{i}^{N} \left| \chi_{i}(\mathbf{r}) 
\right|^{2}
\]
and the total energy is given by
\[
E[\rho]=T_{\mbox{\tiny{S}}}[\rho] + V_{\mbox{\tiny{ne}}}[\rho] + 
J[\rho] \ .
\]
In this way the major contribution of the kinetic enery is described exactly.

Kohn and Sham reformulated the interaction problem so that the 
difference between the exact kinetic energy and the one calculated 
for non-interacting orbitals is absorbed into an {\em 
exchange-correlation} term, including also the exchange part $K[\rho]$.

In general, a DFT energy expression can be written as
\[
E_{\mbox{\tiny{DFT}}}= T_{\mbox{\tiny{S}}}[\rho] + 
V_{\mbox{\tiny{ne}}}[\rho] + J[\rho] + E_{\mbox{\tiny{xc}}}[\rho] \ .
\label{DFT-energy}
\]
The exchange term,$E_{\mbox{\tiny{xc}}}[\rho]$ is defined by equation (\ref{DFT-energy}).

This exchange correlation energy can not completely be compared with 
the terms of Hartree-Fock theory, because the definitions are different. 
The major problem in this is the proper handling of long and short 
range terms of the different functionals. But despite that, the problem 
remains in principle the same. One has to iteratively find the set of 
orbitals, yielding the density, which minimizes the energy to a 
certain threshold.

Again the method of choice is the use of Lagrange 
multipliers. The equations are arranged with a one-electron operator, 
similar to the Fock-operator and the set of coefficients making the 
eigenvalue equations stationary with respect to the energy must be 
found. This yields the Kohn-Sham orbitals of the system, which will 
describe the electron density correctly and, as a result, yield  
the correct energy. The eigenvalue equations are known as {\em 
Kohn-Sham equations}.

Because the KS-orbitals only have to decribe the correct density as 
opposed to the HF-orbitals, which have to describe the complete 
electronic properties, the 
requirements for their accuracy are not as high as in Hartree-Fock, 
which allows the use of lower order basis sets and therefore in this 
treatment some computational cost can be saved.

The problem which still remains is to find the correct exchange functional. 
There is no systematic way to find such a functional. If the exact 
functional was known, DFT would provide the exact energy.
The major advance of DFT is, that the more exact description of the 
electron behaviour is included in the Hamiltonian as a 
correction term, as  opposed to post-HF methods, where the 
wavefunction is extended. That leads to a much simpler treatment, as the functional 
can be defined as most practical and accurate for ones purposes. 

As already stated, there is no systematic way for defining 
correct functionals, but there are a number  of well known 
approximations 
which yield results in the range of high order wave mechanics.

\subsubsection{Local Density Approximation (LDA)}

A famous and
very simple approximation is the {\em Local Density 
Approximation}. It is assumed that the density locally can be 
treated as that of an uniform electron gas, which is equivalent to 
the view 
that the density is a slowly varying function. Then the exchange 
functional is defined by
\[
E_{\mbox{\tiny{xc}}}^{\mbox{\tiny{LDA}}}[\rho]=\int 
\epsilon_{\mbox{\tiny{xc}}}(\rho(\mathbf{r})) \, \rho(\mathbf{r}) \, 
\mathrm{d}\mathbf{r} \ . 
\label{LDA}
\]
Here $\epsilon_{\mbox{\tiny{xc}}}(\rho(\mathbf{r}))$ is the 
exchange-correlation energy of a uniform interacting electron gas of 
density $\rho(\mathbf{r})$. This quantity is known to a very high 
accuracy ($\approx 0.1\%$).\cite{Kohn-JPhysChem-1996}.

It is\cite{Kohn-Rev.Mod.Phys.-1999}
\[
\epsilon_{\mbox{\tiny{x}}}(\rho(\mathbf{r}))=-\frac{0.458}{r_{\mbox{\tiny{s}}}} \ ,
\]
the {\em exchange} part, given in atomic units, and
\[
\epsilon_{\mbox{\tiny{c}}}(\rho(\mathbf{r}))=-\frac{0.44}{r_{\mbox{\tiny{s}}} +7.8} \ , 
\]
the {\em correlation} part. $r_{\mbox{\tiny{s}}}$ is the sphere 
containing one electron and given by 
$r_{\mbox{\tiny{s}}}=\sqrt[3]{3/4\pi}$.

LDA leads to extremly useful results for most applications. Experience 
has
shown\cite{Kohn-Rev.Mod.Phys.-1999}
that LDA gives ionisation energies of atoms and dissociation energies 
of molecules with a fair accuracy, typically of $10$--$20\%$. 
Bond lengths and thus the geometry of molecules have typically 
accuracies of $\approx 1\%$ for main-group molecules. 
The results are not so good for transition metals.

The solution of the Kohn-Sham equations in the LDA
is only minimally 
more difficult than the the solution of the Hartree equation and very much 
easier than the solution of the Hartree-Fock equations. The accuracy 
for the exchange part is about $10\%$, while the much smaller correlation 
part is generally overestimated with by factor of two. The two errors 
typically cancel partially.

Based on the LDA more sophisticated exchange functionals have been 
developed, for which accuracy is much higher. The principle is the 
so-called {\em generalised gradient correction} (GGA), where a more 
accurate behaviour of the electron density and its variation with 
respect to the position is taken into account. Some recent and well used 
functionals are {\em hybrid} functionals, where some degree of {\em 
``exact''} exchange energy is incorporated using traditional wave 
mechanics. 

\subsection{The method of choice}

DFT becomes more favourable with a greaer number of electrons in the system. For this reason DFT calculations are often 
performed for systems including metal centres, or very large systems. For 
the evaluation of molecular geometries, even the very simple LDA 
obtains useful results. The more simple the calculations are, the faster 
they are.

For all of these reasons, the calculations presented in this report 
have been performed with DFT methods in the LDA.