\documentstyle[12pt,epsf]{article}
\def\fract#1#2{\leavevmode
\kern.1em \raise .5ex \hbox{\the\scriptfont0 #1}%
\kern-.1em $/$%
\kern-.15em \lower .25ex \hbox{\the\scriptfont0 #2}}%
\def\fractm#1#2{\leavevmode
\kern.1em \raise .5ex \hbox{\the\scriptfont0 #1}%
\kern-.1em /
\kern-.15em \lower .25ex \hbox{\the\scriptfont0 #2}}%
\def\Lra{$\Longrightarrow$}
\def\Ua{$\rule[-0.5em]{0.0in}{2em}{\Large \Uparrow}$}
\def\Da{$\rule[-0.5em]{0.0in}{2em}\Downarrow$}
\begin{document}
\centerline{\Large Simplified and Biased Introduction}
\centerline{\Large to Density Functional Approaches in Chemistry}
\centerline{\Large\bf THIS WAS WRITTEN IN 1996 -- IT IS OBSOLETE!!!}
{\bf There are many typos in this paper -- Igor Vilfan from
J. Stefan Institute (igor.vilfan@ijs.si) found a score of them, but I
am sure there are more: Thank you Igor!!!.}
\vskip 0.2in
\centerline{\large Jan K. Labanowski ({\tt jkl@ccl.net})}
\vskip 0.2in
\noindent
{\bf Ohio Supercomputer Center, 1224 Kinnear Rd, Columbus, OH 43221-1153}\\
\vskip 0.2in
\noindent
{\sl (Unreviewed, uncorrected, and unfinished draft currently in preparation)}
\vskip 0.2in
There are many approaches of computational chemistry which are popular
in molecular modeling in biological sciences:
\begin{itemize}
\item {\bf Simple Comparative Approaches} -- graphical inspection,
molecular superposition, overlapping/nonoverlaping volume, topological
indices, traditional QSAR, rigid conformational search, ComFA, shape
analysis, etc. Used as a first step in scanning biologically active
molecules and useful in detecting characteristic molecular
features needed for activity. These methods are not quantitative,
since they do not consider energetics of receptor-ligand
interactions.
\item {\bf Emiprical approaches} -- molecular mechanics,
molecular dynamics. Using simple models of harmonic potential,
electrostatic interaction, and dispersion forces, allow for
basic comparisons of energetics and geometry optimization.
Solvent effects can be included explicitly or via empirical models.
Very useful and fast, compared to rigorous quantum calculations.
Major drawbacks: experimental information needed to "standardize"
models and parameters. In principle, these approaches are not able to
model chemical reactions, bond forming/breaking -- electronic structure
does not enter these models.
\item{\bf Quantum Approaches} -- based on explicit consideration
of the electronic structure. These methods are substantially
more computationally demanding then comparative and empirical
approaches for the molecules of the same size. They can be roughly
divided into:
\begin{itemize}
\item {\em Semiempirical methods} -- approximate methods in which
some quantities are taken from experiment, some small quantities
are neglected, and some quantities estimated by fitting to experimental
data. May only be used for chemical species for which they were parametrized.
For distorted, uncommon bonding situations produce unreliable results.
\item {\em Nonempirical methods} -- do not require empirical parameters
and can be used for any molecular system.
\begin{itemize}
\item {\tt Traditional ab initio} -- use Hartree-Fock method as a starting
point, i.e., wave function is used to describe electronic structure.
\item {\tt Density Functional Methods} -- electron density as a
primary way of describing the system.
\end{itemize}
\end{itemize}
\end{itemize}
There is a justified interest in Density Functional Approaches among chemists.
While they have been a domain of physicists for a long time, these methods have now
found their way into mainstream chamistry. The traditional {\em ab initio}
approaches, which are the workhorse of Quantum Chemistry offer a prescription
for calculated truth. In principle, one can calculate chemical properties
with any desired accuracy. The methods are known and proven to work.
The only problem is that for systems larger than hydrogen or lighter H$_{n}$X molecules,
the calculations involved in obtaining accurate results are frequently impractical.
We know how to compute it, but we do not have computational power to do it (probably
for many years to come, even with the spectacular progress in computer technology).
For this reason, the contemporary research in traditional {\em ab initio} methods is
concerned mainly with better approximations to a Full CI method, or infinite order
perturbational expansions, which would give good quality reliable results
with reasonable computational effort. Still, the most promising Coupled Cluster (CC)
and Complete Active Space SCF (CASSCF) calculations scale more than 5th power
in molecular size, and are impractical for molecules containing tens of atoms.
Density Functional Theory does not provide a prescription how to calculate
truth. It only provides the existence proof for possibility of obtaining
accurate results, but no prescription for systematic improvement. The DFT
is accurate, if we knew how to derive necessary relations between density
and energy. Unfortunately,
energy functionals relating electronic density to energy are unknown,
and there is no general way to improve them beside trying new ones and
judging their quality by the results. As Prof. Perdew once said:
"We are like a blind person who wants to find out how the elephant looks like
by touching its legs" (sorry, the quatation is from memory, and may not
be accurate). But the Density Functional Theory provides a hope for a method
which scales with the size as $N^3$ in the worst case, and possibly linearly for
larger molecules (Zhou, 1995; Yang, 1991). So it is used in the
spirit of late prof. Slater saying (about his X$\alpha$ method): Do you want to
calculate it, or do you want it to be accurate? The DFT results are in many
cases surprisingly good if one takes into account the crude approximations
on which some of them are based. For example, Local Spin Density calculations
yield results on many molecular properties which are of quality comparable
with higher order {\em ab initio} correlated methods, yet, the LSD assumes that
electrons in molecules behave like electrons in an uniform electron gas.
Moreover, DFT methods frequently provide reliable answer for cases which are
especially difficult to address with conventional {\em ab initio} methods,
like, e.g., transition metals. On the other hand, they frequently fail miserably,
e.g., in charge-transfer complexes.
The fact that more, and more {\em ab initio} packages provide options for
DFT calculations, is a sign of changing times, and the indication that
even the most vigorous opponents of this method see that it is just
another way of doing things. On the other hand, DFT is in principle only
applicable to the ground state, and there is little hope that it will
be extended in a practical way to excited states in a straightforward manner
any time soon.
\vskip 0.2in
{\large\bf Wave Functions}\\
Since inception of quantum mechanics by Heisenberg, Born, and Jordan
in 1925, and Schr\"odinger in 1926, there were basically two competing
approaches to find the energy of a system of electrons. One was rooted in
statisticial mechanics and the fundamental variable was the total electron
density $\rho({\bf r})$, i.e., the number of electrons per unit volume
at a given point in space (e.g., in cartesian coordinates:
${\bf r} = (x, y, z)$). In this approach,
electrons were treated as particles forming a special gas, called {\em electron gas}.
The special case, the uniform electron gas, corresponds to
the $\rho({\bf r})=const$.
Another approach was to derive the many particle
wave function $\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N, t)$
(where the ${\bf r}_1$ denotes the coordinates of the 1st electron,
${\bf r}_2$ the 2nd electron, and so on, and $t$ is time) and solve the
stationary (time-independent) Schr\"odinger equation for the system:
\begin{equation}
\hat{H}\Psi_k({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) =
E_k \Psi_k({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)
\label{shreq}
\end{equation}
(where $\hat{H}$ is the hamiltonian, i.e.,
the operator of the total energy for the system),
and calculate the set of possible
wave functions (called eigenfunctions) $\Psi_k$ and corresponding
energies (eigenvalues) $E_k$. The eigenfunctions have to be
physically acceptable, and for finite systems:
\begin{itemize}
\item[1.] they should be continuous functions,
\item[2.] they should be at least doubly differentiable,
\item[3.] its square should be integrable,
\item[4.] they should vanish at infinity (for finite systems),
\end{itemize}
When the Schr\"odinger equation is solved {\bf exactly}
(e.g., for the hydrogen atom), the resulting eigenfunctions
$E_k$ form a {\em complete set} of functions, i.e., there is an
infinite number of them, they are orthogonal to each other (or can be
easily made orthogonal through a linear transformation), and
any function which is of ``physical quality'' can be expressed through
the combination of these eigenfunctions. Orthogonal means, that:
\begin{equation}
\int \Psi_k^\ast \Psi_l d^N{\bf r} = 0\ \ {\rm if}\ k\ne l
\end{equation}
The eigenfunction $\Psi_0$ corresponding to the lowest energy $E_0$,
describes the ground state of the system,
and the higher energy values correspond to
excited states. Physicists like to call $\Psi$'s the states (since they
contain all possible information about the state), while chemists
usually call them wave functions.
Once the function
$\Psi$ (or its approximation, i.e., in case
when the Schr\"odinger equation is
solved only approximately) is known,
the corresponding energy of the system can be calculated
as an expectation value of the hamiltonian $\hat{H}$, as:
\begin{equation}
E = \frac{\int \int \cdots \int \Psi^\ast({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)
\hat{H} \Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) d {\bf r_1}
d {\bf r_2} \cdots d {\bf r_N}}
{\int \int \cdots \int \Psi^\ast({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)
\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) d {\bf r_1}
d {\bf r_2} \cdots d {\bf r_N}}
\label{expeq}
\end{equation}
where $\Psi^\ast$ denotes the complex conjugate of $\Psi$ since in general
these functions may produce complex numbers. This is needed since
the operator in this case represents a physical observable, and
the result has to be a real number.
This equation is frequently written using Dirac {\em bra} ($\left<|\right.$) and
{\em ket} ($\left.|\right>$) notation to save space (and to confuse the innocent):
\begin{equation}
E=\frac{\left<\Psi|H|\Psi\right>}{\left<\Psi|\Psi\right>}
\end{equation}
And if $\left<\Psi|\Psi\right> = 1$, i.e., the wave function is normalized, the
equation looks even simpler:
\begin{equation}
E=\left<\Psi|H|\Psi\right>
\end{equation}
\begin{figure}
\epsfbox{coord50.eps}
\caption{Volume element for a particle}
\end{figure}
\noindent
Once we know the wave function $\Psi$ for a given state of our system, we can
calculate the expectation value for any quantity for which we can write down
the operator. The wave function itself, does not correspond to any physical quantity,
but its square represents the probablity density. In other words:
\begin{equation}
|\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)|^2 d{\bf r}_1 d{\bf r}_2 \cdots d{\bf r}_N
\end{equation}
or
\begin{equation}
\Psi^\ast({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)
d{\bf r}_1 d{\bf r}_2 \cdots d{\bf r}_N
\end{equation}
or
\begin{equation}
\left.|\Psi\right>\left<\Psi|\right. d V
\end{equation}
represents the probablity that electron 1 is in the volume element $d{\bf r}_1$ around
point ${\bf r}_1$, electron 2 is in the volume element of the size $d{\bf r}_2$
around point ${\bf r}_2$, and so on. If $\Psi$ describes the system containing only a single
electron, the $|\Psi({\bf r})|^2 d{\bf r}$ simply represents the probability of finding
an electron in the volume element of a size $d{\bf r}$ centered around point ${\bf r}$.
If you use cartesian coordinates, then $d {\bf r} = dx dy dz$ and the volume element would be
a brick (rectangular parallelipiped) with dimensions $dx\times dy \times dx$ whose
vertex closes to the origine of coordinate system is located at $(x, y, z)$.
Now, if we integrate the function $\Psi$ over all the space for all the variables (i.e., sum
up the probablilities in all the elements $d{\bf r}_i$), we should get a probability
of finding our electrons anywhere in the Universe, i.e., 1. This is why it is a good
idea to normalize fundtion $\Psi$. If it is not normalized, it can easily be done by
multiplying it by a normalization constant:
\begin{equation}
\Psi_{normalized} = \overbrace{\frac{1}{\sqrt{\left<\Psi_{unnormalized}|\Psi_{unnormalized}\right>}}}^{Normalization\ constant}
\Psi_{unnormalized}
\end{equation}
Since square of $\Psi$ represents the probablility density of
finding electrons, one may suspect, that it should be easy to
calculate the total electron density from it. And actually it is:
\begin{equation}
\rho({\bf r}) = N\left<\Psi|\delta({\bf r} - {\bf r}_i)|\Psi\right>
\end{equation}
where $N$ is the total number of electrons, and $\delta({\bf r} - {\bf r^{\prime}})$
is the famous Dirac delta function. In cartesians, it simply amounts to
integrating over all electron positions vectors ${\bf r}_i$ but one.
Which one, is not important, since
electrons are indistinguishable, and a proper wave function has to reflect this:
\begin{equation}
\rho({\bf r}_1) = N \underbrace{\int \int \cdots \int}_{N-1} |\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N, )|^2
\underbrace{d{\bf r}_2 d{\bf r}_3 \cdots d{\bf r}_N}_{N-1}
\end{equation}
It is interesting to note, that for the wave function which describes the system containing only
a single electron (but only then !!!):
\begin{equation}
\rho({\bf r}) = \Psi^\ast({\bf r})\Psi({\bf r}) = |\Psi({\bf r})|^2 = \left.|\Psi\right>\left<\Psi\right.|
\label{onero}
\end{equation}
i.e., logically, the electron density and the probability density of finding the single
electron are the same thing.
\vskip 0.3in
{\large\bf Function, Operator, Functional}\\
Before we move any further, let us introduce a few definitions.\\
{\bf FUNCTION}\\
Functions is a prescription which maps one or more
numbers to another number.
For example: take a number and multiply it by itself: $y = f(x) = x^2$, or
take two numbers and add them together: $z = g(x, y) = x + y$.
Sometimes function does not have a value for some numbers, and only
certain numbers can be used as an argument for a function. E.g., square root
is only defined for nonnegative numbers (if you want to have a real number as
a result).
\vskip 0.2in
{\bf OPERATOR}\\
Operator (usually written with a hat, e.g., $ \hat{F} $ or in
calligraphic style $\cal F$ ) is
a prescription which maps one function to another function.
For example: take a function and square its value:
$\hat{F} = {}^2$, e.g., $\hat{F}sin(x) = sin^2(x)$. Or
calculate second derivative of a function versus
$x$: $\hat{F} = \frac{\partial^2}{\partial x^2}$, and
$\hat{F} f(x) = \frac{\partial^2}{\partial x^2} f(x) = \frac{\partial^2 f(x)}{\partial x^2}$.
Nabla is the popular differenial operator in 3 dimensions. In cartesians it is:
\[ \nabla = \left(\vec{\bf i}\frac{\partial}{\partial x} + \vec{\bf j}\frac{\partial}{\partial y} +
\vec{\bf k}\frac{\partial}{\partial z}\right)\]
It is used to calculate forces (which are vectors),
i.e., gradients of potential energy: $Force = - grad\,\, V = -\nabla$.
Its square is called laplacian, and
represents the sum of second derivatives:
\[\Delta = \nabla^2 = \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}
+ \frac{\partial^2}{\partial z^2}\right)\]
The $\nabla^2$ or $\Delta$ appears in the kinetic energy operator.
The prescription of forming the quantum mechanical operators, are called Jordan rules.
For cartesian coordinate representation (coordinate space), they are obtained as follows:
\begin{itemize}
\item[1.] write a classical expression for the physical quantity and rearrange it
in such a way, that everything depends either on coordinates, or momenta (e.g., if
something depends on the component of velocity, $v_{x_i}$, change it to $p_{x_i}/m_i$).
\item[2.] replace coordinates with the operators of multuplying by the coordinate:
\[x_i \rightarrow \hat{x_i} = x_i \cdot\]
\[y_i \rightarrow \hat{y_i} = y_i \cdot\]
\[z_i \rightarrow \hat{z_i} = z_i \cdot\]
\item[3.] replace components of momenta with their operators:
\[p_{x_i} \rightarrow \hat{p}_{x_i} = -i\hbar\frac{\partial}{\partial x}\]
\[p_{y_i} \rightarrow \hat{p}_{y_i} = -i\hbar\frac{\partial}{\partial y}\]
\[p_{z_i} \rightarrow \hat{p}_{z_i} = -i\hbar\frac{\partial}{\partial z}\]
\end{itemize}
The operators can also be obtained in momentum space -- the physicists like them this
way very much. Chemists are interested more in where is the electron, rather than
how fast it moves, so they use coordinate space representation as a rule.
The Schr\"odinger equation (\ref{shreq}) is an example of an eigenproblem, i.e., the
equation in which an operator acts on a function and as a result it returns
the same function multiplied by a constant. For some operators, there are no
nontrivial solutions (trivial means: $\hat{F}0 = C\cdot 0$),
but for operators which correspond to some physical
quantity of some physical system, these equations
have solutions in principle, i.e., one can find a set of
functions, and corresponding constants, which satisfy them.
While these eigenproblems
have solutions in principle, these equations may not be
easy to solve. They frequently have the form of partial second order
differential equations, or integro-differential equations, and they
may not be analytically solved in general, some special cases
may have an analytic solution (e.g., one particle, or two particles
which interact in a special way).
\vskip 0.2in
{\bf FUNCTIONAL}\\
Functional takes a function and provides a number. It is usually
written with the function in square brackets as $F[f] = a$. For example:
take a function and integrate it from $-\infty$ to $+\infty$:
$F[f] = \int\limits_{-\infty}^{+\infty} f(x) dx$. Note that the formula for
the expectation value (\ref{expeq}) is the total energy functional
$E[\Psi]$, since it takes some function $\Psi$ and returns the value of
energy for this $\Psi$.
Functionals can also have derivatives, which behave similarly to
traditional derivatives for functions. The differential of the functional is
defined as:
\begin{equation}
\delta F[f] = F[f + \delta f] - F[f] = \int \frac{\delta F}{\delta f(x)}
\delta f(x) dx
\end{equation}
The functional derivatives have properties similar to traditional function
derivatives, e.g.:
\begin{equation}
\frac{\delta}{\delta f(x)}(C_1F_1 + C_2F_2) =
C_1\frac{\delta F_1}{\delta f(x)} + C_2 \frac{\delta F_2}{\delta f(x)}
\end{equation}
\begin{equation}
\frac{\delta}{\delta f(x)}(F_1F_2) \frac{\delta F_1}{\delta f(x)}F_2 +
\frac{\delta F_2}{\delta f(x)}F_1
\end{equation}
\vskip 0.2in
{\bf Hamiltonian, variational principle, Hartree and Hartree-Fock method}\\
Let us say more about quantum mechanical approches to calculating energy of
the system of electrons. One very important principle of quantum mechanics
is called {\em variational principle} or {\em variational theorem} which leads to
a method ({\em variational method} or {\em variation method}).
It was discovered very early in the history of quanum mechanics. The variational
principle states that if we take some wave function $\Psi$ (e.g., some
approximate one) for a system and calculate the expectation value
of energy $E$, this energy will be either higher, or equal to the
ground state energy $E_0$ for this system.
\begin{equation}
E = \frac{\left<\Psi|H|\Psi\right>}{\left<\Psi|\Psi\right>}\ \ \ \ge \ \ \ E_0 = \frac{\left<\Psi_0|H|\Psi_0\right>}{\left<\Psi_0|\Psi_0\right>}
\end{equation}
The $E = E_0$ occurs only if $\Psi$ is equivalent to $\Psi_0$, but otherwise,
$E > E_0$. This theorem can be easily proven (see for example: Levine, 1983; Slater, 1968, or
Szabo \& Ostlund, 1989).
What is more important however, that this theorem provides us with the prescription
how to find the good wave function: try to go as low with the expectation value of
energy, as you can (you cannot go lower than the true value). And the lower you go,
the better you are. The word of caution here... The principle is true only
if we are applying a true hamiltonian in the expression for expectation energy.
Of course, if your hamiltonian does not represent a physical system, you can
get whatever you wish, $-\infty$ included!!!
As said before, we can rarely solve the Schr\"odinger equation exactly.
And we have to use approximations.
The first successful attempt to derive approximate wave functions
for atoms was devised by Hartree (1928).
In this approach the many-electron wave function $\Psi$ is approximated
by the product of one-electron functions $\phi$ for each of the $N$ electrons:
\begin{equation}
\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) = \phi_1({\bf r}_1)
\phi_2({\bf r}_2) \cdots \phi_N({\bf r}_N)
\label{hart}
\end{equation}
In this equation, ${\bf r}_i$ are the positional coordinates and a spin coordinate
for the $i$-th electron. For example, in cartesians ${\bf r}_i = (x_i, y_i, z_i, m_{s_i})$
though to solve this equation for atoms, the spherical coordinates are more
useful ${\bf r}_i = (r_i, \theta_i, \varphi_i, m_{s_i})$,
where $m_{s_i}$ can adopt only 2 values: $+\fract{1}{2}$ (spin up,
or $\alpha$ or $\uparrow$) or $-\fract{1}{2}$ (spin down, $\beta$, or
$\downarrow$). The individual one-electron functions $\phi_i$ are called orbitals,
(or spinorbitals) and describe each electron in the atom (or molecule).
Is this a reasonable approximation? Not really...
\begin{itemize}
\item[1.] It assumes that electrons in the atom can be described independently,
i.e., their movements do not depend upon each other and their
interaction is not pairwise, but each electron interacts with some
averaged field of other electrons. This is a neat idea, the problem
is that it is not true. Electrons have to avoid each another (correlate
their movements), since they repel each other being of the same charge.
\item[2.] The function does not have a proper symmetry for interchanging particle
indices for fermions. It was discovered long time ago that the
many electron wave function has to be antisymmetric to the exchange of neighboring
indices, i.e., change sign:
\begin{equation}
\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_i, {\bf r}_{i+1}, \cdots {\bf r}_N) =
- \Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_{i+1}, {\bf r}_i, \cdots {\bf r}_N)
\end{equation}
e.g., $\Psi({\bf r}_1, {\bf r}_2, {\bf r}_3) = - \Psi({\bf r}_1, {\bf r}_3, {\bf r}_2)
= - \Psi({\bf r}_2, {\bf r}_1, {\bf r}_3)$.
\end{itemize}
Does the Hartree method work? Yes, it works quite well, at least for atoms...
Now, let us describe how it works.
For an atom or molecule, the hamiltonian operator can be written as:
\begin{equation}
\hat{H}_{tot} = \hat{T}_{nucl} + \hat{T}_{e} + \hat{U}_{nucl} + \hat{V}_{ext} + \hat{U}_{ee}
\end{equation}
where\\
\begin{list}{}{\setlength{\itemindent}{-0.3in}\setlength{\labelwidth}{0.5in}\setlength{\leftmargin}{0.3in}}
\item $\hat{T}_{nucl}$ is the operator of kinetic energy of nuclei,
\item $\hat{T}_{e}$ represents the kinetic energy of electrons,
\item $\hat{U}_{nucl}$ is the interaction energy of nuclei (Coulombic repulsion),
\item $\hat{V}_{ext}$ is the external potential (in this case, the
electrostatic potential coming form nuclei, with which electrons interact),
\item$\hat{U}_{ee}$ denotes electrostatic repulsion between electrons.
\end{list}
Since nuclei are much heavier than electrons (a proton is
1836 times heavier than an electron),
they have much larger intertia than electrons, and we can consider that electrons will adjust
their positions instantly whenever nuclei move. Prof. Rod Bartlett once used an analogy of
flies around the wedding cakes -- you move the cakes, and the flies move with them.
Without much error, we can therefore separate the movement of electrons and nuclei, and
assume that the movement of electrons depends on positions of nuclei in a parametric way.
This is the contents of the Born-Oppenheimer (1927) approximation. It allows us to use the electronic
hamiltonian $\hat{H}_{el}$ to study electrons, and add the components of energy coming from nuclei
(i.e., their kinetic energy, $T_{nucl}$, and internuclear repulsion energy, $V_{nucl}$) at the
end of calculations.
\begin{equation}
\hat{H}_{el} = \hat{T}_{e} + \hat{V}_{ext} + \hat{U}_{ee}
\label{hel}
\end{equation}
The form of these operators will be given below. Atomic units are used here,
i.e., mass of electron $m_e = 1$;
$\hbar = 1$; length is expressed in bohrs (1 bohr = 0.529177249\AA{}; it is the
radius of the first Bohr orbit in hydrogen atom);
energy is in hartrees (1 hartree is 2 rydbergs, and the ground state energy
of hydrogen atom is $-1$ rydberg, or 0.5 hartree, 1 hartree = 627.51 kcal/mol = 2625.5 kJ/mol);
1 unit of charge is a charge of a proton, the gaussian electrostatic system is used (the unit
of permittivity is $4\pi\epsilon_0$ and it does not appear in the Coulomb law for vacuum).
\begin{equation}
\hat{T}_{e} = \sum_{i=1}^{N}-\frac{1}{2}\nabla^2_i
\end{equation}
\begin{equation}
\hat{V}_{ext} = \sum_{i=1}^{N}\hat{v}_i =
\sum_{i=1}^{N}\underbrace{\left(\sum_{\alpha=1}^{N_{nucl}}\frac{-Z_\alpha}{|{\bf r}_i - {\bf R}_\alpha|}\right)}_{v_i}
\end{equation}
\begin{equation}
\hat{U}_{ee} = \sum_{i=1}^{N-1}\sum_{j=i+1}^{N} \frac{1}{|{\bf r}_i - {\bf r}_j|}
\end{equation}
Here, the $Z_\alpha$ is the charge of an $\alpha$-th nucleus
(atomic number), $|{\bf r}_i - {\bf R}_\alpha|$ is the distance between electron $i$ and the nucleus
$\alpha$, $N_{nucl}$ is the total number of nuclei in the molecule
(but it is obviously equal to 1 for an atom), $|{\bf r}_i - {\bf r}_j|$ is the distance between
electron $i$ and $j$. Note, that $|{\bf r}_i - {\bf R}_\alpha|$ and $|{\bf r}_i - {\bf r}_j|$ can be expressed in cartesian coordinates:
\[|{\bf r}_i - {\bf R}_\alpha| = \sqrt{(x_i - X_\alpha)^2 + (y_i - Y_\alpha)^2 + (z_i - Z_\alpha)^2}\]
(note that nuclear coordinates for each of the $N_{nucl}$ nuclei:
$X_\alpha$, $Y_\alpha$, $Z_\alpha$ are constants/parameters rather
than variables in the electronic hamiltonian $\hat{H}_{el}$),
\[|{\bf r}_i - {\bf r}_j| = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2}\]
We can rewrite again the electronic hamiltionian:
\begin{equation}
\hat{H}_{el} = \sum_{i=1}^{N} \hat{h}_i + \hat{U}_{ee}
\end{equation}
where operator
\begin{equation}
\hat{h}_i = -\frac{1}{2}\nabla^2_i + \hat{v}_i
\label{oneh}
\end{equation}
and $\hat{h}_i$ only depends on coordinates of ${\bf r}_i$ of
a given $i$-th electron. In the product function from equation (\ref{hart}), $\hat{h}_i$
would only affect the
function for the $i$-th electron. Unfortunately, there is still this $\hat{U}_{ee}$ which depends on
pairs of electrons, and we cannot separate variables in the the Schr\"odinger equation,
(a second order differential equation). So, Hartree found an approximation, in which
electron does not interact with other electrons one by one,
but with the averaged density of electrons. If we assume for the moment
that we know the one-electron functions $\phi_i({\bf r})$, for each electron
(they make the total wave function in equation (\ref{hart}))
we can calculate densities corresponding to each electron
using equation (\ref{onero}):
\begin{equation}
\rho_i({\bf r}) = |\phi_i({\bf r})|^2
\end{equation}
And the total density of the electrons will be a sum of individual electron densities:
\begin{equation}
\rho_{tot}({\bf r}) = \sum_{i=1}^{N}\rho_i({\bf r}) = \sum_{i=1}^{N}|\phi_i({\bf r})|^2
\end{equation}
However, the $k$-th electron does not interact with the whole density $\rho_{tot}$, since
it is itself a part of it. Electron cannot interact with itself!
It is not a human being!!! This would be a self-interaction. So, if we want to
find a correct density with which the $k$-th electron interacts (let us denote
it by $\rho^{(k)}({\bf r})$), we need to subtract its own
density from $\rho_{tot}$:
\begin{equation}
\rho^{(k)}({\bf r}) = \rho_{tot}({\bf r}) - \rho_k({\bf r})
= \left(\sum_{i=1}^{N}\rho_i({\bf r})\right) - |\phi_k({\bf r})|^2 = \sum_{i=1 \atop i \ne k}^{N} |\phi_i({\bf r})|^2
\end{equation}
Now, let us assume that we know some approximate functions for the orbitals.
And we want to calculate the energy of interation of point charge (our electron)
located at position ${\bf r}$, with other electrons represented as a smeared electron density.
Out point charge is $-e$, which in atomic units is $-1$, and the electron density is also negative,
so we get a positive contribution to energy:
\begin{equation}
\hat{g}_k({\bf r}) = \int \rho^{(k)}({\bf r^\prime})\frac{1}{|{\bf r} - {\bf r}^{\prime}|} d{\bf r}^\prime
\end{equation}
Now, if we make such an approximation, we can write:
\begin{equation}
\hat{U}_{ee} \approx \sum_{i=1}^N \hat{g}_i({\bf r})
\end{equation}
(it is not entirely true, as we will see later, since we doubly count
interaction, but let us leave it at that for a moment),
and our $\hat{H}_{el}$ consists of sums of one-electron operators:
\begin{equation}
\hat{H}_{el} \approx \sum_{i=1}^{N} (-\frac{1}{2}\nabla_i^2 + \hat{v}_i + \hat{g}_i)
\end{equation}
and the many electron Schr\"odinger equation can be solved as $N$ independent
one electron equations:
\begin{equation}
(-\frac{1}{2}\nabla_i^2 + \hat{v}_i + \hat{g}_i)\phi_i({\bf r}) = \epsilon_i \phi_i({\bf r})
\end{equation}
The $\epsilon_i$ is the energy of the $i$-th electron.
In practice, we start with some approximate orbitals $\phi_i$ (e.g., from
hydrogen atom). We need them, since the $\hat{g}_i$ depends on them.
We solve all $N$ equations and obtain the $N$ new $\phi^\prime_i$'s.
The new $\phi^\prime_i$'s are different from the old $\phi_i$'s, and we believe
that they are better. And now we have
better approximations for orbitals, and use the new $\phi_i^\prime$'s as a starting point
and repeat the process. At some point, $\phi_i^\prime$'s do not change from iteration
to iteration, and we obtained the {\em self-consistent field} orbitals.
From these orbitals we can form a many electron wave function $\Psi$,
and then calculate the total energy $E$ of the ground state.
Note, that the total energy is not equal to the sum of orbital energies $\epsilon_i$.
We calculate the expectation value of energy using the accurate hamiltonian
$\hat{H}_{el}$ from equation (\ref{hel}).
When we solved equation for $\phi_1$ and $\epsilon_1$ we included the
Coulomb interactions between electron: (1, 2), (1, 3), (1, 4), etc.
When we solved second equation, we included interactions:
(2, 1), (2, 3), (2, 4), etc. But the interaction (1, 2) is the same
as (2, 1), i.e., adding the energies we would count interactions twice.
For this reason, the correct total energy can be represented as:
\begin{equation}
E = \sum_{i=1}^{N} \epsilon_i + \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} J_{ij}
\end{equation}
where $J_{ij}$'s represent Coulomb interaction of electron $i$ and $j$. They are
called Coulomb integrals and are defined as:
\begin{equation}
J_{ij} = \int \int \frac{\rho_i({\bf r}_1) \rho_j({\bf r}_2)}{|{\bf r}_1 - {\bf r}_2|} d {\bf r}_1 d {\bf r}_2 =
\int \int |\phi_i({\bf r}_1)|^2 \frac{1}{|{\bf r}_1 - {\bf r}_2|} |\phi_j({\bf r}_2)|^2 d {\bf r}_1 d {\bf r}_2
\label{coulint}
\end{equation}
\begin{equation}
J _{ij} = \int \int \phi_i^\ast({\bf r}_1) \phi_j^\ast({\bf r}_2) \frac{1}{|{\bf r}_1 - {\bf r}_2|}
\phi_i({\bf r}_1) \phi_j({\bf r}_2) d {\bf r}_1 d {\bf r}_2
\end{equation}
The Hartree approximation works well for atoms. The one-electron functions are
quite good, and allow us to produce an approximate many-electron function for the whole atom.
But the form of the function adopted by Hartree was basically wrong. It was known
at that time that interchanging the electron labels in the wave functions for the system
of electrons has to change its sign. This is not something which can be derived. It is simply
a law of nature. And it is important, since without this property, the function will not
correctly describe the system. Electrons (fermions) have this propery, and we
should use it. In the system
of fermions, no two particles can be described by the same one particle function.
So, few years later, Fock (1930), and independently Slater(1930) proposed a fix to the
Hartree method. They also used one-electron functions, but the total wave function for the
system was not a simple product of orbitals, but an antisymmetrized sum of all the products
which can be obtained by interchanging electron labels. It is conveniently represented as
a determinant (called Slater determinant):
\begin{equation}
\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) = \frac{1}{\sqrt{N!}}\left| \begin{array}{llcl}
\phi_1({\bf r}_1) & \phi_2({\bf r}_1) & \cdots & \phi_N({\bf r}_1)\\
\phi_1({\bf r}_2) & \phi_2({\bf r}_2) & \cdots & \phi_N({\bf r}_2)\\
\phi_1({\bf r}_3) & \phi_2({\bf r}_3) & \cdots & \phi_N({\bf r}_3)\\
\multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} \\
\multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} \\
\multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} & \multicolumn{1}{c}{.} \\
\phi_1({\bf r}_N) & \phi_2({\bf r}_N) & \cdots & \phi_N({\bf r}_N)\\
\end{array} \right|
\label{sladet}
\end{equation}
Let us examine some properties of this function taking a 2-electron case:
\begin{equation}
\Psi({\bf r}_1, {\bf r}_2) = \frac{1}{\sqrt{2}}\left| \begin{array}{ll}
\phi_1({\bf r}_1) & \phi_2({\bf r}_1)\\
\phi_1({\bf r}_2) & \phi_2({\bf r}_2)\\
\end{array} \right| = \frac{1}{\sqrt{2}}\left[\phi_1({\bf r}_1) \phi_2({\bf r}_2) - \phi_2({\bf r}_1) \phi_1({\bf r}_2)\right]
\end{equation}
If we change electron labels $1 \rightarrow 2$ and $2 \rightarrow 1$ we get
\begin{equation}
\Psi({\bf r}_2, {\bf r}_1) = \frac{1}{\sqrt{2}}\left| \begin{array}{ll}
\phi_1({\bf r}_2) & \phi_2({\bf r}_2)\\
\phi_1({\bf r}_1) & \phi_2({\bf r}_1)\\
\end{array} \right| = \frac{1}{\sqrt{2}}\left[\phi_1({\bf r}_2) \phi_2({\bf r}_1) - \phi_2({\bf r}_2) \phi_1({\bf r}_1)\right]
\end{equation}
that is: $\Psi({\bf r}_1, {\bf r}_2) = - \Psi({\bf r}_2, {\bf r}_1)$.
Moreover, if we assume that two electrons are described by the same spinorbital:
$\phi_1 = \phi_2 = \phi$ we get:
\begin{equation}
\Psi({\bf r}_2, {\bf r}_1) = \frac{1}{\sqrt{2}}\left| \begin{array}{ll}
\phi({\bf r}_2) & \phi({\bf r}_2)\\
\phi({\bf r}_1) & \phi({\bf r}_1)\\
\end{array} \right| = \frac{1}{\sqrt{2}}\left[\phi({\bf r}_2) \phi({\bf r}_1) - \phi({\bf r}_2) \phi({\bf r}_1)\right]
\equiv {\bf 0}
\end{equation}
i.e., such wave function whould be zero everywhere, and hence, the probablility of finding
such electrons is zero. So, Slater determinant satisfies the Pauli exclusion
principle that each electron has to be described by a different wave function.
It was originally formulated by Pauli (1925) for electrons in atoms, that no two
electrons can have the same values of all four quantum numbers: $n$, $l$, $m_l$, and $m_s$.
The single determinant wave function is much better than the Hartree product
function as it naturally includes some basic fermion characteristics.
However, it complicates equations compared to Hartree method and introduces a new term,
{\em electron exchange}. The method of finding the best single determinant wave
function for the system is called Hartree-Fock method.
The expectation value of total energy for the single determinant wave function $\Psi$ is given by:
\begin{equation}
E = \left<\Psi|H|\Psi\right> = \sum_{i=1}^{N} H_i + \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N (J_{ij} - K_{ij})
\end{equation}
where:
\begin{equation}
H_i = \int \phi_i^\ast({\bf r})[-\frac{1}{2}\nabla_i^2 + \hat{v}_i]\phi_i({\bf r})d ({\bf r})
\end{equation}
is an element of one-electron operator $\hat{h}_i$ which was defined by equation (\ref{oneh}),
the $J_{ij}$ is the Coulomb integral defined by equation (\ref{coulint}), and $K_{ij}$ is the new
term, called {\em exchange integral}, which is defined by
\begin{equation}
K_{ij} = \int \int \phi_i^\ast({\bf r}_1) \phi_j({\bf r}_1) \frac{1}{|{\bf r}_1 - {\bf r}_2|}
\phi_i({\bf r}_2) \phi_j^\ast({\bf r}_2) d{\bf r}_1 d{\bf r}_2
\label{exch}
\end{equation}
Note, that $K_{ij}$ is similar in form to the $J_{ij}$ but the functions $\phi_i$ and $\phi_j$
were exchanged. Also, electrons $i$ and $j$ have to be of the same spin for $K_{ij}$ to be
nonzero, due to othogonality of their spin parts. It does not have a simple physical
interpretation like $J_{ij}$ (i.e., purely electrostatic interaction
of two charge densities for electron $i$ and $j$).
The exchange integral comes as a result of the determinantal form of $\Psi$ which
sums all possible products of permutations of electrons among orbitals.
$J_{ij}$ and $K_{ij}$ are equal only when $i=j$. This is very
important, since there is
no contribution to the total energy coming from electron self-interaction.
It can be shown that $J_{ij}$'s are larger (or equal to) $K_{ij}$'s and they are positive numbers.
There are essentially two approaches to HF method: purely numerical,
and the one in which orbitals $\phi$ are represented as combinations
of basis functions. The numerical HF was used at some
point to derive orbitals for many-electron atoms. There are also some
numerical HF programs for molecules. However, these methods are
quite intensive. It is much more popular to represents orbitals $\phi$
as a linear expansion into a set of some basis functions $\chi$:
\begin{equation}
\phi_i({\bf r}) = \sum_{k=1}^n C_{ki} \chi_k({\bf r})
\label{basis}
\end{equation}
and optimizing coefficients $C_{ki}$ rather than basis functions.
The derivation of specific equations for Hartree-Fock
method will not be given
here, but can be found in many books, (see, e.g.:
McWeeny \& Sutcliffe, 1969; Parr, 1963; Pilar, 1968;
Slater, 1968; Szabo \& Ostlung, 1989) and landmark papers
of Roothaan (1951, 1960).
In the SFC LCAO MO (Self Consistent Field -- Linear
Combination of Atomic Orbitals, Molecular Orbitals) method,
we are looking for the coefficients $C_{ki}$ which minimize the energy.
They are derived using varational principle, i.e.
the goal here is to find a single determinant function $\Psi$ as represented by
equation (\ref{sladet}) which corresponds to the lowest possible value of energy.
It is done by the variational calculus. The condition of minimum of energy is:
\begin{equation}
\delta E_0 \equiv \delta \left<\Psi_0|H|\Psi_0\right> = 0
\end{equation}
Another condition is that orbitals $\phi$ have to be
orthonormal, i.e.:
\begin{equation}
\left<\phi_i|\phi_j\right> = \delta_{ij}
\end{equation}
where $\delta_{ij}$ is called a {\em Kronecker delta}, and is equal to
1 only if $i=j$, but is zero otherwise.
The problem of minimization of a function with some subsidiary conditions
is done efficiently with the Lagrange's method of undetermined multipliers.
In this method, the constraints are converted into expressions whose
value is zero when the constraint is met. Such expressions are then
multiplied by an undetermined constant and added (or subtracted, according to
taste) to the minimized function. The resulting sum is then minimized:
\begin{equation}
\delta \left[ \left<\Psi_0|H|\Psi_0\right> - E\left(\left<\Psi_0|\Psi_0\right> - 1\right)\right] = 0
\end{equation}
As a result, one gets Hartree-Fock eigenequations.
\begin{equation}
\hat{f}({\bf r}_1)\phi_i({\bf r}_1) = \epsilon_i \phi_i({\bf r}_1)
\label{hfeq}
\end{equation}
where the Fock operator is defined as:
\begin{equation}
\hat{f}({\bf r}_1) = \hat{h}({\bf r}_1) + \sum_{a=1}^N\left[\hat{j}_a({\bf r}_1) - \hat{k}_a({\bf r}_1)\right]
\label{fock}
\end{equation}
In these equations, the fact that operators act on certain coordinates is stressed by writing explicitly
their dependence on coordinates of electron 1 (it could be any electron number, since it does not matter
how we number them). The $\hat{h}({\bf r}_1)$ was defined in equation (\ref{oneh}). The Coulomb
operator,
\begin{equation}
\hat{j}_a({\bf r}_1)\phi_b({\bf r}_1) = \phi_b({\bf r}_1) \int |\phi_a({\bf r}_2)|^2\frac{1}{|{\bf r}_1 - {\bf r}_2|}d {\bf r}_2 =
\phi_b({\bf r}_1) \int \rho_a({\bf r}_2) \frac{1}{|{\bf r}_1 - {\bf r}_2|}d {\bf r}_2
\end{equation}
and the exchange operator
\begin{equation}
\hat{k}_a({\bf r}_1)\phi_b({\bf r}_1) = \phi_a({\bf r}_1)
\int \phi_a^\ast({\bf r}_2) \phi_b({\bf r}_2)\frac{1}{|{\bf r}_1 - {\bf r}_2|}d {\bf r}_2
\end{equation}
are defined by the result of their action on a function.
If we now introduce basis, i.e., replace $\phi_i$'s with their expansions into basis functions
according to equation (\ref{basis}), the eigenproblems become a set of algebraic equations.
Substituting equation (\ref{basis}) into (\ref{hfeq}) one gets:
\begin{equation}
\hat{f}({\bf r}_1)\sum_{k=1}^n C_{ki} \chi_k({\bf r}_1) = \epsilon_i \sum_{k=1}^n C_{ki} \chi_k({\bf r}_1)
\end{equation}
and multiplying on the left by $\chi_l({\bf r}_1)$ and integrating over ${\bf r}_1$ one gets:
\begin{equation}
\sum_{k=1}^n C_{ki} \int \chi_l^\ast({\bf r}_1) \hat{f}({\bf r}_1) \chi_k({\bf r}_1) d {\bf r}_1 =
\epsilon_i \sum_{k=1}^n C_{ki} \int \chi_l^\ast({\bf r}_1)\chi_k({\bf r}_1) d {\bf r}_1
\end{equation}
Integrals are further denoted as:
\begin{equation}
({\bf F})_{lk} = F_{lk} = \int \chi_l^\ast({\bf r}_1) \hat{f}({\bf r}_1) \chi_k({\bf r}_1) d {\bf r}_1
\end{equation}
which is an element of a $n\times n$ Fock matrix ${\bf F}$, and
\begin{equation}
({\bf S})_{lk} = S_{lk} = \int \chi_l^\ast({\bf r}_1)\chi_k({\bf r}_1) d {\bf r}_1
\end{equation}
which is an element of a $n\times n$ overlap integrals matrix ${\bf S}$.
We can obtain $n$ (i.e., number of basis functions) of such equations, and it is convenient to
write them in a matrix form as:
\begin{equation}
{\bf FC = SC\epsilon}
\label{FC}
\end{equation}
The problem now is to find such marix ${\bf C}$ which diagonalizes ${\bf F}$, which is a standard
procedure in linear algebra. First step is to find a matrix ${\bf X}$
which diagonalizes overlap integrals ${\bf S}$, which is a step ortogonalizing the basis set
(i.e., the basis set is converted to linear combinations of original basis functions, which
are mutually orthogonal).
\begin{equation}
{\bf X}^\dagger{\bf S}{\bf X} = {\bf 1}
\end{equation}
Applying this matrix to equation (\ref{FC}) and rearranging it we get
\begin{equation}
{\bf C}^{\prime\dagger}{\bf F}^{\prime}{\bf C}^\prime = {\bf \epsilon}
\label{CFC}
\end{equation}
where ${\bf C} = {\bf X C}^\prime$ and ${\bf F}^\prime = {\bf X}^\dagger {\bf F} {\bf X}$. The efficient routines for
matrix diagonalization (i.e., obtaining ${\bf C}^\prime$ in this case) are widely available.
While they seem simple, the problem is that elements of ${\bf F}$ depend on orbitals $\phi$, i.e.,
one can solve them only through iterative process as in Hartree's method. Moreover, due to Coulomb
and exchange operators, the Fock matrix elements involve a massive number of two-electron integrals of a type:
\begin{equation}
\left = \int \int \chi_i^\ast({\bf r}_1) \chi_j^\ast({\bf r}_2) \frac{1}{|{\bf r}_1 - {\bf r}_2|}
\chi_k({\bf r}_1) \chi_l({\bf r}_2) d {\bf r}_1 d {\bf r}_2
\end{equation}
To escape these difficulties and complexities, people tried for many years to
decribe electron systems using density, rather than wave function.
In the HF approach the two-electron integrals $\left$ dominate the computational effort.
Moreover, the HF is an approximation, as it does not account for dynamic correlation
due to the rigid form of single determinant wave function. To solve the HF equations,
the assumption has to be made that electrons interact with the averaged potential coming
from other electrons, while in fact, the interactions between electrons are pairwise.
In reality, electrons correlate their movements trying to avoid each other, so there is least
amount of electrostatic repulsion.
To account for dynamic correlation, one has to go to correlated methods, which use multideterminant
wave functions, and these scale as fifth, or even greater powers with the size of a system.
While HF method is quite successful for geometries, it fails miserably to
describe bond breaking or forming (for review of correlated {\em ab initio} methods see: Bartlett \& Stanton, 1994).
\newpage
{\large\bf Density Functional Theory -- early developments}\\
For many years, the use of electron density as a fundamental
description of the system was based on intuition rather than
hard proof that this can be done. Electron density is more
attractive (depends only on $x,y,z$, and eventually, there may be
two densities for spin polarized systems, one for spin up electrons
$\rho_\uparrow({\bf r})$
and one for spin down electrons $\rho_\downarrow({\bf r})$,
as opposed to many particle
wave function which depends on all coordinates of all particles, i.e.,
for N electrons, it depends on $3N$ variables (or $4N$
if you count in spin).
The fact that the ground state properties are functionals of the
electron density $\rho({\bf r})$ was proved by Hohenberg and Kohn (1964)
and it provides the basic framework for modern Density Functional
methods.
More specifically, according to the theorem proved by them, the total
ground state energy of an electron system can be written as a functional
of the electronic density, and this energy is at minimum if the density is
an exact density for the groud state. The theorem of HK is an existence
proof of such a functional, but there is no prescription how to construct it.
If we knew the form of this functional accurately, and if it was
not complicated, quantum chemistry would be a done deal.
Unfortunately we do not know the exact form of energy functional.
It is necessary to use approximations regarding parts of
the functional dealing with kinetic energy and exchange and correlation
energies of the system of electrons.
The simplest approximation is the
local density approximation (LDA) which leads to a Thomas-Fermi
(Fermi, 1928; Thomas, 1927) term
for kinetic energy (for review see, e.g., Jones \& Gunnarsson, 1989;
Slater, 1968; March, 1957) and the Dirac (1930)\footnote{actually Bloch (1929) came with it first,
but Dirac is popularly blamed.} term for the exchange
energy. The corresponding functional is called
Thomas-Fermi-Dirac energy. As you see, these developments are not recent and were parallel to the
work done in the wave function approaches. These functionals can be further improved but the results
are not that encouraging for molecular systems. But, on the other hand, the
Thomas-Fermi-Dirac+improvments method is a true density functional method, since all components of
energy are expressed via density alone, without using many particle wave functions.
However, for the time being, it seems that there is no way to avoid wave functions in molecular
calculations and for accurate calculations they have to be used as a mapping step between
the energy and density. For example, the Thomas-Fermii theory does not predict chemical bonds.
While "pure" density functional theories are very usefull in studying solid phase (e.g., conductivity), they
fail to provide meaningful results for molecular systems.
The predecessor of the modern chemical approaches to the DFT was undoubtely the
Slater's X$\alpha$ method formulated in 1951 (Slater, 1951 \& 1974, for review see:
Johnson 1973 \& 1975). It was
developed as an approximate solution to the HF equations.
In this method, the HF exchange was approximated by:
\begin{equation}
E_{X\alpha}[\rho_\uparrow, \rho_\downarrow]
= -\frac{9}{4}\alpha\left(\frac{3}{4\pi}\right)^{\frac{1}{3}} \int [\rho_\uparrow^\frac{4}{3}({\bf r})
+ \rho_\downarrow^\frac{4}{3}({\bf r})]d {\bf r}
\label{LSD}
\end{equation}
The exchange energy $E_{X\alpha}$ is given here are a functional of densities for spin up ($\uparrow$) and
spin down ($\downarrow$) electrons and contains an adjustable parameter $\alpha$. This parameter, was empirically
optimized for each atom of the periodic table (see, e.g., Slater, 1974, Schwartz, 1972 \& 1974)
and its value was between 0.7--0.8 for most atoms. For a special case of homogenous electron gas, its value
is exactly $\fract{2}{3}$ (G\'asp\'ar, 1954).
\vskip 0.2in
{\bf Hohenberg and Kohn theorems}\\
The field of rigorous density functional theory was born in 1964 with
the publication of the Hohenberg and Kohn paper (1964).
They proved the following:
\begin{itemize}
\item[I.] Every observable of a stationary quantum mechanical system (including energy), can
be calculated, in principle exactly, from the ground-state density alone, i.e.,
every observable can be written as a functional of the ground-state density.
\item[II.] The ground state density can be calculated, in principle exactly,
using the variational method involving only density,
\end{itemize}
The original theorems refer to the time independent (stationary) ground state, but
are being extended to excited states and time-dependent potentials (for review see, e.g.,
Gross \& Kurth, 1994).
How these theorems were derived? By quite an original logic. Within a Born-Oppenheimer
approximation, the ground state of the system of electrons is a result of positions of
nuclei. This is obvious, if we look at the $\hat{H}_{el}$ hamiltonian in equation (\ref{hel}).
In this hamiltonian, the kinetic energy of electrons ($\hat{T}_{e}$) and the
electron-electron interaction ($\hat{U}_{ee}$) ``adjust'' themselves to the external
(i.e., coming from nuclei) potential $\hat{V}_{ext}$. Once the $\hat{V}_{ext}$ is in place,
everything else is, including electron density, which simply adjusts itself to give the
lowest possible total energy of the system. The external potential $\hat{V}_{ext}$ is the
only variable term in this equation, and everything else depends indirectly on it.
Hohenberg and Kohn posed a more interesting question, which is quite opposite to
the traditional logic: Is $\hat{V}_{ext}$ uniquely determined
from the knowledge of electron density $\rho({\bf r})$? Can we find out (in principle,
though it need not be easy) where and what the nuclei are, if we know the density $\rho({\bf r})$
in the ground state? Is there a precise mapping from $\rho({\bf r})$ to $\hat{V}_{ext}$? The answer to
this question is: Yes!. Actually the mapping is accurate within a constant, which would not
change anything, since Schr\"odinger equations with $\hat{H}_{el}$ and $\hat{H}_{el} + const$ yields
exactly the same eigenfunctions, i.e., states (it is easy to prove based on the linear property
of the hamiltonian), and the energies will be simply elevated by the value of this $const$. Note that
all energies are known only within some constant, which establishes the frame of reference (e.g.,
we do not include electron-Mars gravitational attraction within most calculations).
Why was this question so important? Because, if this is true, the knowledge
of density provides total information about the system, and formally if we know the density,
we know everything there is to known.
Since $\rho({\bf r})$ determines number of electrons, $N$:
\begin{equation}
N = \int \rho({\bf r}) d {\bf r}
\end{equation}
and $\rho$ determines the $\hat{V}_{ext}$, the knowledge of total density is as good, as
knowledge of $\Psi$, i.e., the wave function describing the state of the system.
They proved it through a contradiction:
\begin{itemize}
\item[1.] Assume that we have an exact ground state density $\rho({\bf r})$,
\item[2.] Assume that the ground state is nondegenerate (i.e., there is only one wave function
$\Psi$ for this ground state (though HK theorems can be easily
extended for degenerate ground
states, see, e.g., Dreizler \& Gross, 1990; Kohn, 1985),
\item[3.] Assume that for the density $\rho({\bf r})$ there are two possible external potentials:
$\hat{V}_{ext}$ and $\hat{V}_{ext}^\prime$, which obviously produce two different hamiltonians:
$\hat{H}_{el}$ and $\hat{H}_{el}^\prime$, respectively. They obviously produce two different
wave functions for the ground state: $\Psi$ and $\Psi^\prime$, respectively. They correspond to
energies:
$E_0 = \left<\Psi|H|\Psi\right>$ and $E_0^\prime = \left<\Psi^\prime|H^\prime|\Psi^\prime\right>$,
respectively.
\item[4.] Now, let us calculate the expectation value of energy
for the $\Psi^\prime$ with the
hamiltonian $\hat{H}$ and use variational theorem:
\begin{equation}
E_0 < \left<\Psi^\prime|H|\Psi^\prime\right> =
\overbrace{\left<\Psi^\prime|H^\prime|\Psi^\prime\right>}^{E_0^\prime} +
\left<\Psi^\prime|H - H^\prime|\Psi^\prime\right> =
E_0^\prime + \int \rho({\bf r})[\hat{V}_{ext} - \hat{V}_{ext}^\prime]d {\bf r}
\label{E0}
\end{equation}
\item[5.] Now let us calculate the expectation value of energy for the $\Psi$ with the
hamiltonian $\hat{H}^\prime$ and use varational theorem:
\begin{equation}
E_0^\prime < \left<\Psi|H^\prime|\Psi\right> = \underbrace{\left<\Psi|H|\Psi\right>}_{E_0} +
\left<\Psi|H^\prime - H|\Psi\right> =
E_0 - \int \rho({\bf r})[\hat{V}_{ext} - \hat{V}_{ext}^\prime]d {\bf r}
\label{E0P}
\end{equation}
\item[6.] By adding equations (\ref{E0}) and (\ref{E0P}) by sides we obtain:
\begin{equation}
E_0 + E_0^\prime < E_0^\prime + E_0
\end{equation}
and it leads to a contradiction.
\end{itemize}
Since now, we know that $\rho({\bf r})$ determines $N$ and $\hat{V}_{ext}$, it also
determines all properties of the ground state, including the kinetic
energy of electrons
$T_e$ and energy of interaction among electrons $U_{ee}$, i.e., the total ground state energy
is a functional of density with the following
components\footnote{Please, excuse my notation here. I use ``hats'' above operators/potentials (e.g., $\hat{V}_{ext}$
is the external potential)
and no ``hats'' above the corresponding energy components (e.g. $V_{ext}$ is the energy corresponding to
external potential). The exception is a hamiltonian for which $\hat{H}$ is the operator, while
$E$ is the value of energy. This may add confusion, since in much of the literature energy components
are denoted with capital letters or subscripted $E$'s, while the operators are in capital letters, and
potentials in lowercase letters. But this notation also leads to a conflict with notation for one-particle
operators, which are written often in lower case.}:
\begin{equation}
E[\rho] = T_e[\rho] + V_{ext}[\rho] + U_{ee}[\rho]
\end{equation}
Additionally, HK grouped together all functionals
which are secondary (i.e., which are responses)
to the $V_{ext}[\rho]$:
\begin{equation}
E[\rho] = V_{ext}[\rho] + F_{HK}[\rho] = \int \rho({\bf r})\hat{V}_{ext}({\bf r})d{\bf r} + F_{HK}[\rho]
\label{fdef}
\end{equation}
The $F_HK$ functional operates only on density and is universal, i.e., its form does not depend on the particular
system under consideration (note that N-representable densities integrate to N, and the information about
the number of electrons can be easily obtained from the density itself).
The second HK theorem provides variational principle in electron density representation $\rho({\bf r})$\footnote{actually,
in the original HK paper, and many papers with a physical slant, the density is denoted as $n({\bf r})$.}.
For a trial density $\tilde{\rho}({\bf r})$
such that $\tilde{\rho}({\bf r}) \ge {\bf 0}$ and for which
$\int \tilde{\rho}({\bf r}) d {\bf r} = N$,
\begin{equation}
E_0 \le E[\tilde{\rho}]
\end{equation}
where $E[\tilde{\rho}]$ is the energy functional. In other words, if some density represents the correct number
of electrons $N$, the total energy calculated from this density cannot be lower than the true energy of the
ground state.
As to the necessary conditions for this theorem, there is still some controversy concerning
the, so called, representability of density. The $N$-rep\-re\-sen\-tab\-ility, i.e., the fact that the trial $\tilde{\rho}$ has to
sum up to $N$ electrons is easy to achieve by simple rescaling. It is
automatically insured if $\rho({\bf r})$ can be mapped to some wave function
(for further discussion see: Parr \& Yang, 1989; Gilbert, 1975; Lieb, 1982;
and Harriman, 1980).
Assuring that the trial density has also $V_{ext}$-representability
(usually denoted in the literature as $v$-representability) is not that easy. Levy (1982) and Lieb (1983)
have shown, that there are some ``reasonable'' trial densities, which are not the ground state densities for any
possible $V_{ext}$ potential, i.e., they do not map to any external potential. Such densities do not correspond therefore
to any ground state, and their optimization will not lead to a ground state.
Moreover, during energy minimization, we may take a wrong turn, and get
stuck into some non $v$-representable density and never be able to converged
to a physically relevant ground state density.
For an interesting discussion, see
Hohenberg {\em et al.} (1990). Assuming that we restrict ourselved only to
trial densities which are both $N$ and $v$ representable, the variational principle for density is easly proven, since
each trial density $\tilde{\rho}$ defines a hamiltonian $\hat{\tilde{H}}_{el}$. From the hamiltonian
we can derive the corresponding wave function $\tilde{\Psi}$ for the ground state represented by this hamiltonian. And
according to the traditional variational principle, this wave function $\tilde{\Psi}$ will not be a ground state
for the hamiltonian of the real system $\hat{H}_{el}$:
\begin{equation}
\tilde{\rho} \rightarrow \hat{\tilde{H}}_{el} \rightarrow \tilde{\Psi};\ \ \left<\right.\tilde{\Psi}|H|\tilde{\Psi}\left.\right>
= E[\tilde{\rho}] \ge E[\rho_0] \equiv E_0
\end{equation}
where $\rho_0({\bf r})$ is the true ground state density of the real system.
The condition of minimum for the energy functional:
$\delta E[\rho({\bf r})] = 0$ needs to be constrained
by the $N$-representability of density which is optimized\footnote{
it also needs to be constrained by v-representability, but
we still do not know how to express v-representability in a closed
mathematical form. There exist, however, methods, e.g.,
{\em constrained search} (Levy, 1982) and {\em local-scaling
transformation} (Petkov {\em et al}, 1986) which assure v-representability
during density optimization, though their algorithmic implementation
needs to be done.}.
The Lagrange's method of undetermined multipliers is a very
convenient approach for the constrained minimization
problems. In this method
we represent constraints in such a way that their value is
exactly zero when they are satisfied. In our case, the
$N$ representability constraint can be represented as:
\begin{equation}
constraint = \int \rho({\bf r}) d {\bf r} - N = 0
\end{equation}
These constraints are then multiplied by an undetermined
constants and added to a minimized function or functional.
\begin{equation}
E[\rho({\bf r})] - \mu \left[\int \rho({\bf r}) d {\bf r} - N\right]
\end{equation}
where $\mu$ is yet undetermined Lagrange multiplier.
Now, we look for the minimum of this expression by requiring that its
differential is equal to zero (a necessary condition of minimum).
\begin{equation}
\delta \left\{ E[\rho({\bf r})] - \mu \left[\int \rho({\bf r}) d {\bf r} - N\right]\right\} = 0
\end{equation}
Solving this differential equation will provide us with a prescription of
finding a minimum which satisfies the constraint. In our case it leads to:
\begin{equation}
\delta E[\rho({\bf r})] - \mu \delta \left\{ \int \rho({\bf r}) d {\bf r}\right\} = 0
\end{equation}
since $\mu$ and $N$ are constants. Using the definition of the
differential of the functional (see, e.g., Parr \& Yang, 1989):
\begin{equation}
F[f + \delta f] - F[f] = \delta F= \int \frac{\delta F}{\delta f(x)}\delta f(x) dx
\end{equation}
and the fact that differential and integral signs may be interchanged, we obtain
\begin{equation}
\int \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})}\delta \rho({\bf r}) d {\bf r}
- \mu \int \delta \rho({\bf r}) d {\bf r} = 0
\end{equation}
Since integration runs over the same variable and has the same limits, we can
write both expressions under the same integral:
\begin{equation}
\int \left\{ \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})}
- \mu \right\} \delta \rho({\bf r}) d {\bf r} = 0
\end{equation}
which provides the condition for constrained minimisation and defines the value
of the Lagrange multiplier at minium. It is also expressed here via external potential from equation (\ref{fdef}):
\begin{equation}
\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})}
= \hat{V}_{ext}({\bf r}) + \frac{\delta F_{HK}(\rho({\bf r})}{\delta \rho({\bf r})}
\label{mudef}
\end{equation}
Density functional theory
gives a firm definition of the chemical potential $\mu$, and leads to
several important general conclusions. For review, please refer to Parr \& Yang (1989), chapters 4 and 5.
\vskip 0.2in
{\bf Kohn and Sham Method}\\
Above equations provide a prescription of minimizing energy by changing corresponding density.
Unfortunately, the expression relating kinetic energy to density is not known with satisfactory
accuracy. The current expressions, even those improved upon from the original
Thomas-Fermi theory, are quite crude and quite unsatisfactory for atoms and molecules in particular.
On the other hand, the kinetic energy is easily calculated from the wave function, provided that it is
known. For that reason, Kohn \& Sham (1965) proposed an ingenious method of marrying wave function and
density approach.
They repartitioned the total energy functional into following parts:
\begin{equation}
E[\rho] = T_0[\rho] + \int \left[\hat{V}_{ext}({\bf r}) + \hat{U}_{cl}({\bf r})\right]\rho({\bf r})d{\bf r} + E_{xc}[\rho]
\label{etot}
\end{equation}
where $T_0[\rho]$ is the kinetic energy of electrons in a system which has the same density $\rho$ as the real
system, but in which there is no electron-electron interactions. This is frequently called a system on noninteracting
electrons, but it may be misunderstood, since electrons still interact with nuclei.
\begin{equation}
\hat{U}_{cl}({\bf r})=\int \frac{\rho({\bf r}^\prime)}{|{\bf r}^\prime - {\bf r}|}d{\bf r}^\prime
\end{equation}
is a pure Coulomb (``classical'') interaction between electrons.
It includes electron self-interaction explicitly, since the corresponding
energy is
\begin{equation}
E_{cl}[\rho] = \int \int \frac{\rho({\bf r}^\prime) \rho({\bf r})}{|{\bf r}^\prime - {\bf r}|}d {\bf r} d{\bf r}^\prime
\end{equation}
and it represents interaction of $\rho$ with itself.
$\hat{V}_{ext}({\bf r})$ is the external potential, i.e., potential coming from nuclei:
\begin{equation}
\hat{V}_{ext} = \sum_{\alpha} \frac{-Z_\alpha}{|{\bf R}_\alpha - {\bf r}|}
\end{equation}
The last functional, $E_{xc}[\rho]$, called exchange-correlation energy, is defined by this equation.
$E_{xc}[\rho]$ includes all the energy contributions which were not accounted for by previous terms, i.e.:
\begin{itemize}
\item electron exchange
\item electron correlation, since non-interacting electrons do need to correlate their movements. Please note,
however, that this correlation component is not the same as defined by L\"owdin for {\em ab initio} methods.
\item a portion of the kinetic energy which is needed to correct $T_0[\rho]$ to obtain true kinetic energy of
a real system $T_e[\rho]$.
\item correction for self-interaction introduced by the classical coulomb potential.
\end{itemize}
In fact, all the difficult things were ``swept under the carpet'' in this functional. However, better and better
approximations for this functional are being published. To finish the derivation of Kohn-Sham equations, let us
assume for a while, that we know this energy functional reasonably well.
In a similar fashion, as was done for the equation defining chemical potential (\ref{mudef}) we may apply the
variational principle and obtain:
\begin{equation}
\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} = \frac{\delta T_0[\rho({\bf r})]}{\delta \rho({\bf r})} +
\hat{V}_{ext}({\bf r}) + \hat{U}_{cl}({\bf r}) + \frac{\delta E_{xc}[\rho({\bf r})]}{\delta \rho({\bf r})}
\end{equation}
There is another trick which can be done with this equation, namely:
\begin{equation}
\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} = \frac{\delta T_0[\rho({\bf r})]}{\delta \rho({\bf r})} +
\hat{V}_{eff}({\bf r})
\label{veff}
\end{equation}
where we lumped together all terms, except our noninteracting electron kinetic energy, into an effective potential
depending upon ${\bf r}$:
\begin{equation}
\hat{V}_{eff}({\bf r}) = \hat{V}_{ext}({\bf r}) + \hat{U}_{cl}({\bf r})
+ \hat{V}_{xc}({\bf r})
\end{equation}
where the exchange correlation potential is defined as a functional derivative of exchange correlation energy:
\begin{equation}
\hat{V}_{xc}({\bf r}) = \frac{\delta E_{xc}[\rho({\bf r})]}{\delta \rho({\bf r})}
\end{equation}
The form of equation (\ref{veff}) asks for a solution as a Schr\"odinger equation for noninteracting particles:
\begin{equation}
\left[ - \frac{1}{2}\nabla^2_i + \hat{V}_{eff}({\bf r})\right]\phi_i^{KS}({\bf r}) = \epsilon_i \phi_i({\bf r})^{KS}
\end{equation}
Yes!!! it is very similar to the eigenequation of the Hartree-Fock method, with one difference, it is much simpler.
The Fock operator in equation (\ref{fock}) contains the potential which is {\bf nonlocal}, i.e., different for each electron.
The Kohn-Sham operator depends only on ${\bf r}$, and not upon the index of the electron. It is the same for all
electrons. The Kohn-Sham orbitals, $\phi_i({\bf r})^{KS}$, which are quite easily
derived from this equation\footnote{if expansion into basis sets is used,
the matrix equation (\ref{CFC}) for expansion coefficients is identical as
in Hartree-Fock method.}, can be used immediately to compute the total density:
\begin{equation}
\rho({\bf r}) = \sum_{i=1}^N |\phi_i^{KS}({\bf r})|^2
\end{equation}
which can be used to calculate an improved potential $\hat{V}_{eff}({\bf r})$, i.e., lead to a new cycle of
{\em self-consistent field}\footnote{It has to be stressed that $\phi_i({\bf r})^{KS}$'s are not the real orbitals, and
they do not correspond to a real physical system. Their only role in the theory is to
provide a mapping between kinetic energy and density. For one, the total KS wave function is a single determinant
and is unable to model situations where more determinants are needed to describe the system (e.g., for
cases when molecules dissociate to atoms). Interesting discussion about the symmetry of this wave
function is given by Dunlap (1991, 1994). More studies are needed to asses physical relevance of these orbitals.}.
Density can also be used to calculate the total energy from equation (\ref{etot}), in which
the kinetic energy $T_0[\rho]$ is calculated from the corresponding orbitals, rather than density itself:
\begin{equation}
T_0[\rho] = \frac{1}{2} \sum_{i=1}^N \left<\phi_i^{KS}|\nabla_i^2|\phi_i^{KS}\right>
\end{equation}
and the rest of the total energy as:
\begin{equation}
V_{eff}[\rho] = \int \hat{V}_{eff}({\bf r})\rho({\bf r})d {\bf r}
\end{equation}
In practice, total energy is calculated more economically using orbital energies $\epsilon_i$ as:
\begin{equation}
E[\rho] = \sum_{i=1}^{N}\epsilon_i - \frac{1}{2} \int \int \frac{\rho({\bf r})\rho({\bf r}^\prime)}{|{\bf r}-{\bf r}^\prime|}
d {\bf r} d {\bf r}^\prime - \int \hat{V}_{xc}({\bf r})\rho({\bf r}) d {\bf r} + E_{xc}[\rho]
\end{equation}
It is a popular misconception to look at this method as describing noninteracting electrons moving
in a potential given by nuclei. In fact, they move in an effective potenitial $\hat{V}_{eff}({\bf r})$
which includes electron interaction, though in an artificial way. While this is more philosophical
than physical, when you look at the form of $\hat{V}_{eff}({\bf r})$, which in Kohn-Sham equations
takes a role of "external potential", the electron-electron interaction is replaced by interaction with some medium
which mimics the electron-electron interaction. This medium actually, in a way,
exaggerates the interaction between electrons, what is evidenced by the fact that the correction
which needs be added to $T_0$\footnote{$\Delta T = T_e - T_0$ is embedded in $E_{xc}$.}
is positive, i.e., the "noninteracting electrons" move slower
than the real, interacting ones.
\vskip 0.2in
{\bf Implementations of Kohn and Sham method}\\
The derivation of Kohn-Sham equations presented here is close to the original presentation, in which nonpolarized electron
density is used, and occupation numbers for Kohn-Sham orbitals are assumed one. However extensions exist both for
polarized spin densities (i.e., different orbitals for spin up, and spin down electrons), and for nonintegral occupation
numbers in the range $(0,1)$. For review see Parr \& Yang (1989) and Salahub {\em et al} (1994).
Several uneasy questions can be asked here. What are the Kohn-Sham orbitals and energies? Formally, they are artifacts with no
real physical significance. However, they are being used, and it seems, that they are quite close to the Hartee-Fock orbitals
(Zhao \& Parr, 1993).
Also, the Kohn-Sham formalism can be extended to fractional occupation numbers $0 \le n_i \le 1$
(Janak, 1978; Perdew \& Zunger, 1981;
Parr \& Yang, 1989). Orbital energies $\epsilon_i$ are in this case:
\begin{equation}
\epsilon_i = \frac{\partial E}{\partial n_i}
\end{equation}
and one application may be to integrate energy from $N-1$ to $N$ electrons, and calculate ionization potential.
The derivatives of energy versus occupation numbers provide also other response functions and can provide
rigurous definitions of chemical potential, electronegativity, softness and hardness,
(Parr, 1994; Parr \& Yang, 1989; Neshev \& Proynov, 1993)
First implementations of the Kohn-Sham method were using the local approximations to the exchange correlation energy.
The appropriate functionals were taken from data on homogenous electron gas.
There were two variants of the method: spin unpolarized (LDF/LDA -- Local Density Functional/Approximation)
and spin polarized (LSD -- Local Spin Density) where arguments require both $\alpha$ and $\beta$ electron
densities, rather than
a total density.
For historical reasons, the exchange correlation energy was partitioned into 2 parts:
\begin{equation}
E_{xc}[\rho] = E_x[\rho] + E_c[\rho]
\end{equation}
the exchange energy, and correlation energy. This partition is quite arbitrary, since exchange and correlation have
slightly different meaning than in {\em ab initio} approaches. The exchange energy in LDF/LSD was approximated with
the homogenous gas exchange result given by equation (\ref{LSD}) with $\alpha=\fract{2}{3}$. The correlation energy
was expressed as:
\begin{equation}
E_c[\rho] = \int \rho({\bf r}) \epsilon_c[\rho_\uparrow({\bf r})\rho_\downarrow({\bf r})]d {\bf r}
\end{equation}
where $\epsilon_c[\rho_\uparrow({\bf r})\rho_\downarrow({\bf r})]$ is the correlation energy per one electron in a gas
with spin densities $\rho_\uparrow({\bf r})$ and $\rho_\downarrow({\bf r})$. This function is not known analytically, but
is constantly improved on the basis of quantum Monte Carlo simmulations, and fitted to analytical expansion
(Vosko {\em et al}, 1980; von Barth \& Hedin, 1979).
The local functionals derived from electron gas data worked suprisingly well, taking into account that they
substantially underestimate the exchange energy (by as much as 15\%) and grossly overestimate the correlation energy,
sometimes by a 100\%. The error in exchange is however larger than the correlation error in absolute values.
LSD/LDF is known to overbind normal atomic bonds, on the other hand, it produces too weak hydrogen bonds.
Early attempts to improve functionals by GEA (Gradient Expansion Approximation), in which $E_{xc}[\rho]$ was
expanded in Taylor series versus $\rho$ and truncated at a linear term, did not improve results too much (see, e.g.
Langreth \& Vosko, 1990). Only GGA (Generalized Gradient Approximation) provided notable improvements by
expanding $E_{xc}[\rho]$. The expansion here is not a simple Taylor expansion, but tries
to find the right asymptotic behaviour and right scaling for the
usually nonlinear expansion. These enhanced functionals are frequently called nonlocal or gradient corrections,
since they depend not only upon the density, but also the magnitude of the gradient (i.e., first derivative)
of density at a given point.
While the term corrections has a historical meaning, where the corrections were added on the top of local
density functionals, it is probably no longer correct, since modern nonlocal functionals are quite complicated
functions in which the value of density and its gradient are integral parts of the formula. There are many different
corrections, and it will take some time when their respective merit, accuracy, and best domains of applications will be
established. For review of the existing functionals see e.g.: Clementi (1994), Johnson {\em et al} (1993),
Seminario \& Politzer (1995), Ziegler (1991). The most widely used local potentials are: Slater for exchange (Slater, 1974),
and VWN for correlation (Vosko {\em et al}, 1980).
Probably the most frequently in use today are:
\begin{itemize}
\item For exchange: B88 (Becke, 1988), PW86 (Perdew \& Wang, 1986)
\item For correlation: P86 (Perdew, 1986), LYP (Lee {\em at al}, 1988),
\end{itemize}
Beside different exchange and correlation functionals, there is recent interest in the
ACM (Adiabatic Connection Method) (see: Harris, 1984; Becke, 1993) to include at least
partially the exact exchange from Hartree-Fock type calculation, or exchange
Kohn-Sham orbitals. This method looks very promising
but its application is not without cost. In the worst case, ACM calculation of
exact exchange via exchange integrals (see equation \ref{exch}) scales as $N^4$.
Most molecular DFT codes use basis functions, with a notable exception of
a totally basis free NUMOL (Becke \& Dickinson, 1990) program. Quite impressive
variety of basis sets is used. Since most energy components have to be
calculated via numerical integration, the use of contracted
gaussians as a basis is not that essential. But most program use them,
as it allows a substantial code reuse from the well developed traditional
{\em ab initio} techniques. The well known {\em deMon} (Salahub {\em et al}, 1992), {\em, DGauss}
(Andzelm \& Wimmer, 1992), and DeFT (a program similar to deMon written aby Alain St-Amant {\em et al.}
and available at\\
{\tt ftp://www.ccl.net/pub/chemistry/software/SOURCES/FORTRAN/DeFT})\\
all use gaussian basis functions. Some traditional {\em ab initio} packages also
include options to perform DFT calculations: ACES2 (Bartlett \& Stanton, 1994),
Cadpac5 (Cambridge Analytic Derivative Package, e.g., Murray {\em et al}, 1993), Gaussian (Gaussian, Inc.),
Spartan (Wavefunction, Inc.), and I hear that many others.
The AMOL, recently renamed to ADF (Amsterdam Density Functional) (Boerrigter {\em et al} 1998)
uses Slater type basis.
The DVM (Ellis \& Painter, 1970) and DMol (Delley, 1995) use numerical basis sets, which
are given as spline functions for radial part (the angular part is taken as appropriate
spherical harmonics).
There are also programs for extended systems like: Corning (Teter {\em et al}, 1989),
Crystal94 (Univ. Turin \& SERC Daresbury Lab), and Wien95 (Blaha {\em et al}, 1995)
which use LAPW (Linearized Augmented Plane Wave's) as basis functions.
{\large\bf NOTE: This stuff was written in 1996 -- there are many new
codes and new versions out there}
\vskip 0.2in
{\bf Typical Program Organization for SCF-KS equations}\\
The single geometry SCF cycle or geometry optimization
involve following steps:
\begin{itemize}
\item[1.] Start with a density (for the 1st iteration, superposition of atomic
densities is typically used).
\item[2.] Establish grid for charge density and exchanger correlation potential
\item[3.] Compute KS matrix (equivalent to the ${\bf F}$ matrix in Hartree-Fock
method in equation (\ref{CFC})) elements and overlap integrals matrix.
\item[4.] Solve the equations for expansion coefficients to obtain KS orbitals.
\item[5.] Calculate new density $\rho = \sum_{i=occ} |\phi_i({\bf r})|^2$.
\item[6.] If density or energy changed substantially, go to step 1.
\item[7.] If SCF cycle converged and geometry optimization is not requested, go to step 10.
\item[8.] Calculate derivatives of energy vs. atom coordinates, and update atom
coordinates. This may require denser integration grids and recomputing of
Coulomb and exchange-correlation potential.
\item[9.] If gradients are still large, or positions of nuclei moved appreciably, go to step 1.
\item[10.] Calculate properties and print results.
\end{itemize}
Of course, there may be other variants of this method (e.g., when one computes
vibrational frequencies from the knowledge of gradients and energies only).
It is quite popular to limit expense of numerical integration during the SCF cycle.
It is frequently done by fitting auxiliary functions to charge density and exchange
correlation potential. This allows for much faster integral evaluation.
These auxiliary fitting functions are usually uncontracted gaussians (though quite different
from the atomic basis sets) for which the integrals required for KS matrix can be
calculated analytically. Different auxilliary sets are used for fitting charge density
and exchange-correlation potential (see e.g., Dunlap \& R\"osch, 1990). The need for
fitting is recently questioned (see e.g., Johnson, 1995) since it scales
as $N^3$ even for very large systems, however, it is still very popular in
DFT codes. The fitting procedures are in general non sparse, while for large
molecules many contributions coming from distant portions may be neglected
leading to less steep scaling with molecular size.
Early DFT codes were impaired by the lack of analytical gradients.
Currently, expressions for first and second derivatives exist
(see e.g.: Dunlap \& Andzelm, Komornicki \& Fitzgerald, 1993)
and are implemented in many programs, thus facilitating geometry
optimization and vibrational frequency calculations.
{\bf Performance of DFT}\\
This will be a very short list of DFT applications, as there are many excellent reviews on this topic
(see e.g.: Labanowski \& Andzelm, 1991; Parr \& Yang, 1989; Seminario \& Politzer 1995; Ziegler, 1991; Bartolotti and Flurchick, 1996; St-Amant. 1996).
The G1 database of Pople and coworkers is a benchmark for accuracy of the traditional {\em ab initio} methods.
The database containes 55 molecules for which experimental values of atomization energies are known within
1 kcal/mol. Curtiss {\em et al} (1991) achieved the 1.2 kcal/mol mean absolute error for these 55 atomization
energies using the G2 procedure, which is a quite involved prescription incorporating higher order correlated methods.
Becke (1992) was able to reproduce values in this database with a mean absolute error of 3.7 kcal/mol using his NUMOL
program with
gradient corrected functionals. This result was additionally improved by Becke (1993) to 2.4 kcal/mol when portion of
the electron exchange entering the echange-correlation energy, $E_{xc}$ was calculated exactly form Kohn-Sham orbitals
While the error in DFT is considered still too big, these results were obtained with a method
which is substantially less computationally demanding than
original correlated {\em ab initio} procedures used by Pople and coworkers. Obviusly, the errors refer to absolute
atomization energies, which in general are very difficult to calculate with good occuracy
(for review see, e.g., Hehre {\em et al}, 1986). The relative differences are usually reproduced much better
with DFT methods.
Even without gradient corrections DFT results for bond dissociation energies
are usually much better then the Hartree-Fock (which routinely underbinds) results,
though they have an overbinding tendency. The LDA results are approximately of MP2
quality. The inclusion of gradient corrections to DFT provides bond dissociaction
energies which pair in occuracy with MP4 and CC results.
Molecular geometries even with LSD are much better than corresponding HF results
and are of the MP2 quality. However, LSD fails to correctly treat hydrogen
bonding. This deficiency is removed when one uses gradient corrected DFT approaches.
Quite estonishingly, DFT provides excellent results for molecules which
are notoriously difficult for traditional {\em ab initio} approaches
like FOOF, FON, and metalorganic or inorganic moities. There seems to be a funny
regularity: "If something does not work with {\em ab initio}, try it with
DFT, and vice versa".
Transition states of organic molecules are frequently not reproduced
well with "pure" DFT. However, it seems that admixture of exact
exchange (see Becke, 1993) via ACM dramatically improves
the problem cases. DFT is however a metod of
choice for transition states in metalorganic reactions.
These systems are notoriously difficult to treat with even
high quality {\em ab initio} and have problems with convergence.
Vibrational frequencies are well reproduced even by LSD, though
gradient corrections improve agreement with experiment even
further.
Ionization potentials, electron affinities, and proton
affinities are reproduced fairly well within
gradient corrected DFT.
Recently, there is much interest in using DFT for
high spin species, since preliminary results
are very promissing. On the other hand good
performance of DFT in this field comes as a surprise,
since high multiplets are poorly described by
a single determinant wave function. For all the wrong reasons,
with KS wave function of broken symmetry, the multiplets
splitting are quite well reproduced by DFT.
The scope of applcations for DFT grows rapidly
with calculations of new molecular properties being
added to actively developed software. Recent extensions
include parameters for NMR and ESR spectroscopy,
diamagnetic properties, polarizabilities, relativistic
calculations, and others.
The DFT methodology is still far from mature. It is encouraging
that several groups which were primarily focused on the
development of traditional {\em ab initio} approaches
are now actively working in the DFT area. It will bring
new methodological developments, as well as, more
carefull assesment of applicability.
With all the hopes which DFT brings, one has to keep
in mind, that it is primarily a ground state oriented
method, and cannot compete with semiempirical and
correlated {\em ab initio} methods in calculations concerning
excited states.
\vskip 0.3in
{\bf Advantages of DFT in biological chemistry}\\
\begin{itemize}
\item Computational demands are much less severe than
for {\em ab initio} methods of similar quality -- hence
the method is applicable to much larger molecules.
\item Metals are frequently present in active centers
of enzymes. Traditional {\em ab initio} methods have
severe problems with transition metals. In fact, it can
be proved that Hartree-Fock equation cannot be solved for
the true metalic state. It is related to the fact that
there is a difficulty to converge HF when highest occupied
orbitals are very close in energy (the situation very
popular for transition metals).
\item The DFT, similar to {\em ab initio} methods,
is nonparametric, i.e., applicable to any molecule.
While some say that basis sets are {\em parameters}
for {\em ab initio} and DFT methods, this is an
exaggeration. Basis sets are easily derived from
atomic calculations, and beside, they were
derived long time ago for all elements of periodic
table.
\item The restriction of DFT being applicable to
the ground state only is not usually a problem, unless
you study interaction of radiation with biological
molecules (e.g., UV induced mutations).
\end{itemize}
\section*{References}
\begin{list}{}{\setlength{\itemindent}{-0.3in}\setlength{\labelwidth}{0.5in}}
\item Andzelm J. and Wimmer E. (1992), {\em J. Chem. Phys.}, {\bf 96}, 1280.
\item Barth, von U., Hedin L. (1979), {\em Phys. Rev. A}, {\bf 20}, 1693.
\item Bartlett R.J., Stanton J.F. (1994), In: {\em Reviews in Computational Chemistry, Volume V},
Lipkowitz, K.B. and Boyd D.B. Editors, pp. 65-169, VCH Publishers, New York.
\item Bartolotti L.J., Flurchick K. (1996), In: {\em Reviews in
Computational Chemistry, Volume VII},
Lipkowitz, K.B. and Boyd D.B. Editors, pp. 187-216, VCH Publishers, New York.
\item Becke A.D. (1988), {\em Phys. Rev. A}, {\bf 38}, 3098.
\item Becke A.D. (1992), {\em J. Chem. Phys.}, {\bf 96}, 2155.
\item Becke A.D. (1993), {\em J. Chem. Phys.}, {\bf 98}, 1372; {\bf 98}, 5648;
\item Becke A.D. and Dickinson R.M. (1990), {\em J. Chem. Phys.}, {\bf 92}, 3610.
\item Blaha P., Schwarz K., Dufek P., Augustyn R. (1995), WIEN95, Technical University of Vienna.
\item Bloch F. (1929), {\em Z. Physik}, {\bf 57}, 545.
\item Boerrigter P.M., te Velde G., Baerends E.J., {\em Int. J. Quant. Chem.},
{\bf 33}, 87.
\item Born M., Oppenheimer J.R. (1927), {\em Ann. Physik}, {\bf 84}, 457.
\item Clementi E., Chakravorty S.J., Corongiu G., Flores J.R., Sonnad V. (1991),
In: {\rm Modern Techniques in Computational Chemistry -- MOTECC-91}, Clementi E., Editor,
pp. 23--113, Escom, Leiden.
\item Clementi E., Corongiu G., Stradella O.G. (1991), {\bf ibid.}, pp. 295--379.
\item Curtis L.A., Raghavachari K., Trucks G.W., Pople J.A. (1991), {\em J. Chem. Phys.}, {\bf 94}, 7221.
\item Dirac P.A.M. (1930), {\em Proc. Cambridge Phil. Soc.}, {\bf 26}, 376.
\item Dreizler R.M. and Gross E.K.U (1990), {\em Density Functional Theory: An Approach to the
Quantum Many-Body Problem}, Springer-Verlag, Berlin.
\item Dunlap B.I. and R\"osch N. (1990), {\em Adv. Quant, Chem.}, {\bf 21}, 317-339.
\item Dunlap B.I. (1991), In: {\em Density Functional Methods in Chemistry}, J.K. Labanowski, J.W. Andzelm, Editors;
pp. 49--60, Springer-Verlag, New York.
\item Dunlap B.I. and Andzelm J., {\em Physical Review A}, {\bf 45}, 81.
\item Dunlap B.I. (1995), In: {\em Modern Density Functional Theory -- A Tool for Chemistry,
Theoretical and Computational Chemistry, Volume 2},
Seminario J. and Politzer P., Editors, pp. 151-167, Elsevier, Amsterdam, 1995.
\item Ellis D.E. \& Painter G.S. (1970), {\em Phys. Rev. B}, {\bf 2}, 2887/
\item G\'asp\'ar R. (1954), {\em Acta Phys. Acad. Sci. Hung.}, {\bf 3}, 263.
\item Hartree D.R. (1928), {\em Proc. Cambridge Phil. Soc.}, {\bf 24}, 89.
\item Fermi E. (1928), {\em Z. Physik}, {\bf 48}, 73.
\item Fock V. (1930), {\em Z. Physik}, {\bf 61}, 126 (1930); {\bf 62}, 795.
\item Gilbert T.L. (1975), {\em Phys. Rev. B}, {\bf 12}, 2111.
\item Gross E.K.U and Kurth S. (1994), In: {\em Relativistic and Electron
Correlation Effects in Molecules, NATO ASI Series, Ser. B, Volume 318},
Malli G.L., Editor, pp. 367-409, Plenum Press, New York.
\item Harriman J.E. (1980), {\em Phys. Rev. A}, {\bf 24}, 680.
\item Harris J. (1984), {\em Phys. Rev. A}, {\bf 29}, 1648.
\item Hehre W.J., Radom L., Schleyer P.v.R, Pople J.A. (1986), {\em Ab Initio
Molecular Orbital Theory}, John Willey \& Sons, New York.
\item Hohenberg P., Kohn W. (1964), {\em Phys. Rev. B}, {\bf 136}, 864.
\item Hohenberg P.C., Kohn W., Sham L.J. (1990), {\em Adv. Quant. Chem.}, {\bf 21}, 7.
\item Janak J.F. (1979), {\em Phys. Rev. B}, {\bf 18}, 7165.
\item Johnson B.G. (1995), In: {\em Modern Density Functional Theory -- A Tool for Chemistry,
Theoretical and Computational Chemistry, Volume 2},
Seminario J. and Politzer P., Editors, pp. 169 -- 219, Elsevier, Amsterdam, 1995.
\item Johnson B.G., Gill P.M.W, Pople J.A. (1993), {\bf J. Chem. Phys.}, {\bf 98}, 5612.
\item Johnson K.H. (1973) {\em Adv. Quant. Chem.}, {\bf 7}, 143.
\item Johnson K.H. (1975) {\em Ann. Rev. Phys. Chem.}, {\bf 26}, 39.
\item Jones R.O. and Gunnarsson O. (1989), {\em Rev Mod. Phys.}, {\bf 61}, 689-746.
\item Kohn W. (1985), In: {\em Highlights of Condensed Matter Theory}, Bassani F., Fumi F.,
Tosi M.P., Editors, p. 1, North-Holland, Amsterdam.
\item Kohn W. and Sham L.J. (1965), {\em Phys. Rev. A}, {\bf 140}, 1133.
\item Komornicki A. and Fitzgerald G. (1993), {\em J. Chem. Phys.}, {\bf 98}, 1398.
\item Labanowski J.K. and Andzelm J.W. (1991) {\em Density Functional Methods in Chemistry},
Springer-Verlag, New York.
\item Langreth D.C., Vosko S.H. (1990), {\em Adv. Quant. Chem.}, {\bf 21}, 175.
\item Lee C., Yang W., Parr R.G. (1988), {\em Phys. Rev. B}, {\bf 37}, 785.
\item Levine I.N. (1983), {\em Quantum Chemistry}, Allyn and Bacon, Boston.
\item Levy M. (1982), {\em Phys. Rev. A}, {\bf 26}, 1200.
\item Lieb E. (1983), {\em Int. J. Quant. Chem.}, {\bf 24}, 243.
\item March N.H. (1957), {\em Adv. Phys.}, {\bf 6}, 1.
\item McWeeny R. and Sutcliffe B.T. (1969), {\em Methods of Molecular Quantum
Mechanics}, Academic Press, London.
\item Murray C.W., Handy N.C., Laming G.C. (1993), {\em Mol. Phys.}, {\bf 78}, 997; {\bf 80}, 1121.
\item Neshev N.M. and Proynov E.I. (1993), {\em J. Mol. Catal.}, {\bf 54}, 484.
\item Parr R.G. (1963), {\em Quantum Theory of Molecular Electronic Structure},
Benjamin, New York.
\item Parr R.G. and Yang W. (1989), {\em Density-Functional Theory of Atoms
and Molecules}, Oxford University Press, New York.
\item Parr R.G. (1994), {\em Phil. Mag.}, {\bf 69}, 737.
\item Pauli W., Jr. (1925), {\em Z. Physik}, {\bf 31}, 765.
\item Perdew J.P. (1986), {\em Phys. Rev. B}, {\bf 33}, 8822; {\bf 34}, 7406E.
\item Perdew J.P. and Wang Y. {\em Phys. Rev. B}, {\bf 33}, 8800.
\item Perdew J.P. and Zunger A. (1981), {\em Phys. Rev. B.}, {\bf 23}, 5041.
\item Petkov I.Z, Stoitsov M., Kryachko E. (1986), {\em Int. J. Quant. Chem.}, {\em 29}, 146.
\item Pilar F.L. (1968), {\em Elementary Quantum Mechanics}, McGraw-Hill, New York.
\item Seminario J. and Politzer P., Editors, (1995), {\em Modern Density Functional Theory -- A Tool for Chemistry,
Theoretical and Computational Chemistry, Volume 2}, Elsevier, Amsterdam, 1995.
\item Roothaan C.C.J. (1951), {\em Rev. Mod. Phys.}, {\bf 23}, 69 (1951);
\item Roothaan C.C.J. (1960), {\em Rev. Mod. Phys.}, {\bf 32}, 179 (1960).
\item Salahub D.R., Fournier R., M\l ynarski P., Papai I., St-Amant A., Ushio J. (1991),
In: {\em Density Functional Methods in Chemistry}, Labanowski J. and Andzelm J., Eds.,
pp. 77-100, Springer-Verlag, New York.
\item Salahub D.R., Castro M., Proynov E.I. (1994), In: {\em Relativistic and Electron
Correlation Effects in Molecules, NATO ASI Series, Ser. B, Volume 318},
Malli G.L., Editor, pp. 411-445, Plenum Press, New York.
\item Schr\"odinger E. (1926), {\em Ann. Phys.}, {\bf 79}, 361, 486;
{\bf 80}, 437; {\bf 81}, 109.
\item Schwartz K. (1972), {\em Phys. Rev. B}, {\bf 5}, 2466.
\item Schwartz K. (1974), {\em Theor. Chim. Acta}, {\bf 34}, 225.
\item Slater J.C. (1930), {\em Phys. Rev.}, {\bf 35}, 210.
\item Slater J.C. (1951), {\em Phys. Rev.}, {\bf 81}, 385.
\item Slater J.C. (1968), {\em Quantum Theory of Matter}, McGraw-Hill, New York.
\item Slater J.C. (1974), {\em Quantum Theory of Molecules and Solids, Volume 4,
The self-consistent field for molecules and solids}, McGraw-Hill, New York.
\item St-Amant, A. (1996), In: {\em Reviews in
Computational Chemistry, Volume VII},
Lipkowitz, K.B. and Boyd D.B. Editors, pp. 217-259, VCH Publishers, New York.
\item Szabo A., and Ostlund N.S. (1989), {\em Modern Quantum Chemistry}, NcGraw-Hill,
New York.
\item Teter M.P., Payne M.C., Allan D.C. (1989), {\em Phys. Rev. B}, {\bf 40}, 12255.
\item Thomas L.H. (1927), {\em Proc. Cambridge Phil. Soc.}, {\bf 23}, 542.
\item Vosko S.J., Wilk L., Nusair M. (1980), {\em Can. J. Chem.}, {\bf 58}, 1200.
\item Yang W. (1991), {\em Phys. Rev. Lett.}, {\bf 66}, 1438.
\item Zhou Z. (1995), {\em Modern Density Functional Theory -- A Tool for Chemistry,
Theoretical and Computational Chemistry, Volume 2}, Seminario J. and Politzer P., Editors,
pp. 125 -- 150, Elsevier, Amsterdam, 1995.
\item Zhow Q., Parr R.G. (1933), {\em J. Chem. Phys.}, {\bf 98}, 1.
\item Ziegler T. (1991), {\em Chem. Rev.}, {\bf 91}, 651.
\end{list}
\end{document}