Simplified and Biased Introduction

to Density Functional Approaches in Chemistry

THIS WAS WRITTEN IN 1996 - IT IS OBSOLETE!!!

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!!!.

Jan K. Labanowski (jkl@ccl.net)

Ohio Supercomputer Center, 1224 Kinnear Rd, Columbus, OH 43221-1153
(Unreviewed, uncorrected, and unfinished draft currently in preparation)

There are many approaches of computational chemistry which are popular in molecular modeling in biological sciences:

• 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.

• 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.

• 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:
  • 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.

  • Nonempirical methods - do not require empirical parameters and can be used for any molecular system.
    • Traditional ab initio - use Hartree-Fock method as a starting point, i.e., wave function is used to describe electronic structure.
    • Density Functional Methods - electron density as a primary way of describing the system.

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 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 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 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 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 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.

Wave Functions
Since inception of quantum mechanics by Heisenberg, Born, and Jordan in 1925, and Schrödinger 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 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ödinger equation for the system:

$$\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)$$ (1)

(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:

1.

they should be continuous functions,

2.

they should be at least doubly differentiable,

3.

its square should be integrable,

4.

they should vanish at infinity (for finite systems),

When the Schrödinger equation is solved exactly (e.g., for the hydrogen atom), the resulting eigenfunctions $E_k$ form a 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:

$$\int \Psi_k^\ast \Psi_l d^N{\bf r} = 0 {\rm if} k\ne l$$ (2)

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ödinger 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:

$$E = \frac{\int \int \cdots \int \Psi^\ast({\bf r}_1, {\bf r}... ... \cdots {\bf r}_N) d {\bf r_1} d {\bf r_2} \cdots d {\bf r_N}}$$ (3)

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 bra ( $\left<\vert\right.$) and ket ( $\left.\vert\right>$) notation to save space (and to confuse the innocent):

$$E=\frac{\left<\Psi\vert H\vert\Psi\right>}{\left<\Psi\vert\Psi\right>}$$ (4)

And if $\left<\Psi\vert\Psi\right> = 1$, i.e., the wave function is normalized, the equation looks even simpler:

$$E=\left<\Psi\vert H\vert\Psi\right>$$ (5)

Figure 1: Volume element for a particle

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:

$$\vert\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)\vert^2 d{\bf r}_1 d{\bf r}_2 \cdots d{\bf r}_N$$ (6)

or

$$\Psi^\ast({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N)\Psi({\bf r... ...}_2, \cdots {\bf r}_N) d{\bf r}_1 d{\bf r}_2 \cdots d{\bf r}_N$$ (7)

or

$$\left.\vert\Psi\right>\left<\Psi\vert\right. d V$$ (8)

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 $\vert\Psi({\bf r})\vert^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:

$$\Psi_{normalized} = \overbrace{\frac{1}{\sqrt{\left<\Psi_{un... ...lized}\right>}}}^{Normalization constant} \Psi_{unnormalized}$$ (9)

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:

$$\rho({\bf r}) = N\left<\Psi\vert\delta({\bf r} - {\bf r}_i)\vert\Psi\right>$$ (10)

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:

$$\rho({\bf r}_1) = N \underbrace{\int \int \cdots \int}_{N-1}... ...t^2 \underbrace{d{\bf r}_2 d{\bf r}_3 \cdots d{\bf r}_N}_{N-1}$$ (11)

It is interesting to note, that for the wave function which describes the system containing only a single electron (but only then !!!):

$$\rho({\bf r}) = \Psi^\ast({\bf r})\Psi({\bf r}) = \vert\Psi({\bf r})\vert^2 = \left.\vert\Psi\right>\left<\Psi\right.\vert$$ (12)

i.e., logically, the electron density and the probability density of finding the single electron are the same thing. Function, Operator, Functional
Before we move any further, let us introduce a few definitions.

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). 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{... ...al}{\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} + \... ...rtial^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:

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$).

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$$

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}$$

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ödinger equation (1) 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).

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 (3) 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:

$$\delta F[f] = F[f + \delta f] - F[f] = \int \frac{\delta F}{\delta f(x)} \delta f(x) dx$$ (13)

The functional derivatives have properties similar to traditional function derivatives, e.g.:

$$\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)}$$ (14)

$$\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$$ (15)

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 variational principle or variational theorem which leads to a method (variational method or 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.

$$E = \frac{\left<\Psi\vert H\vert\Psi\right>}{\left<\Psi\vert... ...si_0\vert H\vert\Psi_0\right>}{\left<\Psi_0\vert\Psi_0\right>}$$ (16)

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ödinger 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:

$$\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)$$ (17)

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: $+\leavevmode \kern.1em \raise .5ex \hbox{\the\scriptfont0 1}\kern-.1em$/ $\kern-.15em \lower .25ex \hbox{\the\scriptfont0 2}$ (spin up, or $\alpha$ or $\uparrow$) or $-\leavevmode \kern.1em \raise .5ex \hbox{\the\scriptfont0 1}\kern-.1em$/ $\kern-.15em \lower .25ex \hbox{\the\scriptfont0 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...

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.

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:

$$\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_i, {\bf r}_{i+1}, ... ... {\bf r}_2, \cdots {\bf r}_{i+1}, {\bf r}_i, \cdots {\bf r}_N)$$ (18)

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)$.

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:

$$\hat{H}_{tot} = \hat{T}_{nucl} + \hat{T}_{e} + \hat{U}_{nucl} + \hat{V}_{ext} + \hat{U}_{ee}$$ (19)

where

• $\hat{T}_{nucl}$ is the operator of kinetic energy of nuclei,
• $\hat{T}_{e}$ represents the kinetic energy of electrons,
• $\hat{U}_{nucl}$ is the interaction energy of nuclei (Coulombic repulsion),
• $\hat{V}_{ext}$ is the external potential (in this case, the electrostatic potential coming form nuclei, with which electrons interact),
• $\hat{U}_{ee}$ denotes electrostatic repulsion between electrons.

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.

$$\hat{H}_{el} = \hat{T}_{e} + \hat{V}_{ext} + \hat{U}_{ee}$$ (20)

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Å; 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).

$$\hat{T}_{e} = \sum_{i=1}^{N}-\frac{1}{2}\nabla^2_i$$ (21)

$$\hat{V}_{ext} = \sum_{i=1}^{N}\hat{v}_i = \sum_{i=1}^{N}\un... ...-Z_\alpha}{\vert{\bf r}_i - {\bf R}_\alpha\vert}\right)}_{v_i}$$ (22)

$$\hat{U}_{ee} = \sum_{i=1}^{N-1}\sum_{j=i+1}^{N} \frac{1}{\vert{\bf r}_i - {\bf r}_j\vert}$$ (23)

Here, the $Z_\alpha$ is the charge of an $\alpha$-th nucleus (atomic number), $\vert{\bf r}_i - {\bf R}_\alpha\vert$ 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), $\vert{\bf r}_i - {\bf r}_j\vert$ is the distance between electron $i$ and $j$. Note, that $\vert{\bf r}_i - {\bf R}_\alpha\vert$ and $\vert{\bf r}_i - {\bf r}_j\vert$ can be expressed in cartesian coordinates:

$$\vert{\bf r}_i - {\bf R}_\alpha\vert = \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}$),

$$\vert{\bf r}_i - {\bf r}_j\vert = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2}$$

We can rewrite again the electronic hamiltionian:

$$\hat{H}_{el} = \sum_{i=1}^{N} \hat{h}_i + \hat{U}_{ee}$$ (24)

where operator

$$\hat{h}_i = -\frac{1}{2}\nabla^2_i + \hat{v}_i$$ (25)

and $\hat{h}_i$ only depends on coordinates of ${\bf r}_i$ of a given $i$-th electron. In the product function from equation (17), $\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ödinger 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 (17)) we can calculate densities corresponding to each electron using equation (12):

$$\rho_i({\bf r}) = \vert\phi_i({\bf r})\vert^2$$ (26)

And the total density of the electrons will be a sum of individual electron densities:

$$\rho_{tot}({\bf r}) = \sum_{i=1}^{N}\rho_i({\bf r}) = \sum_{i=1}^{N}\vert\phi_i({\bf r})\vert^2$$ (27)

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}$:

$$\rho^{(k)}({\bf r}) = \rho_{tot}({\bf r}) - \rho_k({\bf r}) ... ...t^2 = \sum_{i=1 \atop i \ne k}^{N} \vert\phi_i({\bf r})\vert^2$$ (28)

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:

$$\hat{g}_k({\bf r}) = \int \rho^{(k)}({\bf r^\prime})\frac{1}{\vert{\bf r} - {\bf r}^{\prime}\vert} d{\bf r}^\prime$$ (29)

Now, if we make such an approximation, we can write:

$$\hat{U}_{ee} \approx \sum_{i=1}^N \hat{g}_i({\bf r})$$ (30)

(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:

$$\hat{H}_{el} \approx \sum_{i=1}^{N} (-\frac{1}{2}\nabla_i^2 + \hat{v}_i + \hat{g}_i)$$ (31)

and the many electron Schrödinger equation can be solved as $N$ independent one electron equations:

$$(-\frac{1}{2}\nabla_i^2 + \hat{v}_i + \hat{g}_i)\phi_i({\bf r}) = \epsilon_i \phi_i({\bf r})$$ (32)

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 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 (20). 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:

$$E = \sum_{i=1}^{N} \epsilon_i + \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} J_{ij}$$ (33)

where $J_{ij}$'s represent Coulomb interaction of electron $i$ and $j$. They are called Coulomb integrals and are defined as:

$$J_{ij} = \int \int \frac{\rho_i({\bf r}_1) \rho_j({\bf r}_2)... ..._2\vert} \vert\phi_j({\bf r}_2)\vert^2 d {\bf r}_1 d {\bf r}_2$$ (34)

$$J _{ij} = \int \int \phi_i^\ast({\bf r}_1) \phi_j^\ast({\bf ... ...} \phi_i({\bf r}_1) \phi_j({\bf r}_2) d {\bf r}_1 d {\bf r}_2$$ (35)

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):

$$\Psi({\bf r}_1, {\bf r}_2, \cdots {\bf r}_N) = \frac{1}{\sqr... ... r}_N) & \cdots & \phi_N({\bf r}_N)\\ \end{array} \right\vert$$ (36)

Let us examine some properties of this function taking a 2-electron case:

$$\Psi({\bf r}_1, {\bf r}_2) = \frac{1}{\sqrt{2}}\left\vert \b... ...\phi_2({\bf r}_2) - \phi_2({\bf r}_1) \phi_1({\bf r}_2)\right]$$ (37)

If we change electron labels $1 \rightarrow 2$ and $2 \rightarrow 1$ we get

$$\Psi({\bf r}_2, {\bf r}_1) = \frac{1}{\sqrt{2}}\left\vert \b... ...\phi_2({\bf r}_1) - \phi_2({\bf r}_2) \phi_1({\bf r}_1)\right]$$ (38)

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:

$$\Psi({\bf r}_2, {\bf r}_1) = \frac{1}{\sqrt{2}}\left\vert \b... ...f r}_1) - \phi({\bf r}_2) \phi({\bf r}_1)\right] \equiv {\bf0}$$ (39)

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, 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:

$$E = \left<\Psi\vert H\vert\Psi\right> = \sum_{i=1}^{N} H_i + \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N (J_{ij} - K_{ij})$$ (40)

where:

$$H_i = \int \phi_i^\ast({\bf r})[-\frac{1}{2}\nabla_i^2 + \hat{v}_i]\phi_i({\bf r})d ({\bf r})$$ (41)

is an element of one-electron operator $\hat{h}_i$ which was defined by equation (25), the $J_{ij}$ is the Coulomb integral defined by equation (34), and $K_{ij}$ is the new term, called exchange integral, which is defined by

$$K_{ij} = \int \int \phi_i^\ast({\bf r}_1) \phi_j({\bf r}_1) ... ...\phi_i({\bf r}_2) \phi_j^\ast({\bf r}_2) d{\bf r}_1 d{\bf r}_2$$ (42)

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$:

$$\phi_i({\bf r}) = \sum_{k=1}^n C_{ki} \chi_k({\bf r})$$ (43)

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 (36) which corresponds to the lowest possible value of energy. It is done by the variational calculus. The condition of minimum of energy is:

$$\delta E_0 \equiv \delta \left<\Psi_0\vert H\vert\Psi_0\right> = 0$$ (44)

Another condition is that orbitals $\phi$ have to be orthonormal, i.e.:

$$\left<\phi_i\vert\phi_j\right> = \delta_{ij}$$ (45)

where $\delta_{ij}$ is called a 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:

$$\delta \left[ \left<\Psi_0\vert H\vert\Psi_0\right> - E\left(\left<\Psi_0\vert\Psi_0\right> - 1\right)\right] = 0$$ (46)

As a result, one gets Hartree-Fock eigenequations.

$$\hat{f}({\bf r}_1)\phi_i({\bf r}_1) = \epsilon_i \phi_i({\bf r}_1)$$ (47)

where the Fock operator is defined as:

$$\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]$$ (48)

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 (25). The Coulomb operator,

$$\hat{j}_a({\bf r}_1)\phi_b({\bf r}_1) = \phi_b({\bf r}_1) \i... ...\bf r}_2) \frac{1}{\vert{\bf r}_1 - {\bf r}_2\vert}d {\bf r}_2$$ (49)

and the exchange operator

$$\hat{k}_a({\bf r}_1)\phi_b({\bf r}_1) = \phi_a({\bf r}_1) \... ...{\bf r}_2)\frac{1}{\vert{\bf r}_1 - {\bf r}_2\vert}d {\bf r}_2$$ (50)

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 (43), the eigenproblems become a set of algebraic equations. Substituting equation (43) into (47) one gets:

$$\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)$$ (51)

and multiplying on the left by $\chi_l({\bf r}_1)$ and integrating over ${\bf r}_1$ one gets:

$$\sum_{k=1}^n C_{ki} \int \chi_l^\ast({\bf r}_1) \hat{f}({\bf... ..._{ki} \int \chi_l^\ast({\bf r}_1)\chi_k({\bf r}_1) d {\bf r}_1$$ (52)

Integrals are further denoted as:

$$({\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$$ (53)

which is an element of a $n\times n$ Fock matrix ${\bf F}$, and

$$({\bf S})_{lk} = S_{lk} = \int \chi_l^\ast({\bf r}_1)\chi_k({\bf r}_1) d {\bf r}_1$$ (54)

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:

$${\bf FC = SC\epsilon}$$ (55)

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).

$${\bf X}^\dagger{\bf S}{\bf X} = {\bf 1}$$ (56)

Applying this matrix to equation (55) and rearranging it we get

$${\bf C}^{\prime\dagger}{\bf F}^{\prime}{\bf C}^\prime = {\bf\epsilon}$$ (57)

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:

$$\left<ij\vert kl\right> = \int \int \chi_i^\ast({\bf r}_1) \... ...t} \chi_k({\bf r}_1) \chi_l({\bf r}_2) d {\bf r}_1 d {\bf r}_2$$ (58)

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<ij\vert kl\right>$ 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 ab initio methods see: Bartlett & Stanton, 1994). 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)¹ 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:

$$E_{X\alpha}[\rho_\uparrow, \rho_\downarrow] = -\frac{9}{4}\... ...}{3}({\bf r}) + \rho_\downarrow^\frac{4}{3}({\bf r})]d {\bf r}$$ (59)

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 $\leavevmode \kern.1em \raise .5ex \hbox{\the\scriptfont0 2}\kern-.1em$/ $\kern-.15em \lower .25ex \hbox{\the\scriptfont0 3}$ (Gáspár, 1954). 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:

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.

II.

The ground state density can be calculated, in principle exactly, using the variational method involving only density,

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 (20). 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ödinger 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$:

$$N = \int \rho({\bf r}) d {\bf r}$$ (60)

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:

1.

Assume that we have an exact ground state density $\rho({\bf r})$,

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),

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\vert H\vert\Psi\right>$ and $E_0^\prime = \left<\Psi^\prime\vert H^\prime\vert\Psi^\prime\right>$, respectively.

4.

Now, let us calculate the expectation value of energy for the $\Psi^\prime$ with the hamiltonian $\hat{H}$ and use variational theorem:

$$E_0 < \left<\Psi^\prime\vert H\vert\Psi^\prime\right> = \ove... ...t \rho({\bf r})[\hat{V}_{ext} - \hat{V}_{ext}^\prime]d {\bf r}$$ (61)

5.

Now let us calculate the expectation value of energy for the $\Psi$ with the hamiltonian $\hat{H}^\prime$ and use varational theorem:

$$E_0^\prime < \left<\Psi\vert H^\prime\vert\Psi\right> = \und... ...t \rho({\bf r})[\hat{V}_{ext} - \hat{V}_{ext}^\prime]d {\bf r}$$ (62)

6.

By adding equations (61) and (62) by sides we obtain:

$$E_0 + E_0^\prime < E_0^\prime + E_0$$ (63)

and it leads to a contradiction.

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²:

$$E[\rho] = T_e[\rho] + V_{ext}[\rho] + U_{ee}[\rho]$$ (64)

Additionally, HK grouped together all functionals which are secondary (i.e., which are responses) to the $V_{ext}[\rho]$:

$$E[\rho] = V_{ext}[\rho] + F_{HK}[\rho] = \int \rho({\bf r})\hat{V}_{ext}({\bf r})d{\bf r} + F_{HK}[\rho]$$ (65)

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})$³. For a trial density $\tilde{\rho}({\bf r})$ such that $\tilde{\rho}({\bf r}) \ge {\bf0}$ and for which $\int \tilde{\rho}({\bf r}) d {\bf r} = N$,

$$E_0 \le E[\tilde{\rho}]$$ (66)

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$-representability, 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 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}$:

$$\tilde{\rho} \rightarrow \hat{\tilde{H}}_{el} \rightarrow \t... ...{\Psi}\left.\right> = E[\tilde{\rho}] \ge E[\rho_0] \equiv E_0$$ (67)

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⁴. 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:

$$constraint = \int \rho({\bf r}) d {\bf r} - N = 0$$ (68)

These constraints are then multiplied by an undetermined constants and added to a minimized function or functional.

$$E[\rho({\bf r})] - \mu \left[\int \rho({\bf r}) d {\bf r} - N\right]$$ (69)

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).

$$\delta \left\{ E[\rho({\bf r})] - \mu \left[\int \rho({\bf r}) d {\bf r} - N\right]\right\} = 0$$ (70)

Solving this differential equation will provide us with a prescription of finding a minimum which satisfies the constraint. In our case it leads to:

$$\delta E[\rho({\bf r})] - \mu \delta \left\{ \int \rho({\bf r}) d {\bf r}\right\} = 0$$ (71)

since $\mu$ and $N$ are constants. Using the definition of the differential of the functional (see, e.g., Parr & Yang, 1989):

$$F[f + \delta f] - F[f] = \delta F= \int \frac{\delta F}{\delta f(x)}\delta f(x) dx$$ (72)

and the fact that differential and integral signs may be interchanged, we obtain

$$\int \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})}\de... ...f r}) d {\bf r} - \mu \int \delta \rho({\bf r}) d {\bf r} = 0$$ (73)

Since integration runs over the same variable and has the same limits, we can write both expressions under the same integral:

$$\int \left\{ \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} - \mu \right\} \delta \rho({\bf r}) d {\bf r} = 0$$ (74)

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 (65):

$$\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} =... ...r}) + \frac{\delta F_{HK}(\rho({\bf r})}{\delta \rho({\bf r})}$$ (75)

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. 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:

$$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]$$ (76)

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.

$$\hat{U}_{cl}({\bf r})=\int \frac{\rho({\bf r}^\prime)}{\vert{\bf r}^\prime - {\bf r}\vert}d{\bf r}^\prime$$ (77)

is a pure Coulomb (``classical'') interaction between electrons. It includes electron self-interaction explicitly, since the corresponding energy is

$$E_{cl}[\rho] = \int \int \frac{\rho({\bf r}^\prime) \rho({\b... ...}{\vert{\bf r}^\prime - {\bf r}\vert}d {\bf r} d{\bf r}^\prime$$ (78)

and it represents interaction of $\rho$ with itself. $\hat{V}_{ext}({\bf r})$ is the external potential, i.e., potential coming from nuclei:

$$\hat{V}_{ext} = \sum_{\alpha} \frac{-Z_\alpha}{\vert{\bf R}_\alpha - {\bf r}\vert}$$ (79)

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.:

• electron exchange
• 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öwdin for ab initio methods.
• 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]$.
• correction for self-interaction introduced by the classical coulomb potential.

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 (75) we may apply the variational principle and obtain:

$$\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} =... ...}) + \frac{\delta E_{xc}[\rho({\bf r})]}{\delta \rho({\bf r})}$$ (80)

There is another trick which can be done with this equation, namely:

$$\mu = \frac{\delta E[\rho({\bf r})]}{\delta \rho({\bf r})} =... ...\rho({\bf r})]}{\delta \rho({\bf r})} + \hat{V}_{eff}({\bf r})$$ (81)

where we lumped together all terms, except our noninteracting electron kinetic energy, into an effective potential depending upon ${\bf r}$:

$$\hat{V}_{eff}({\bf r}) = \hat{V}_{ext}({\bf r}) + \hat{U}_{cl}({\bf r}) + \hat{V}_{xc}({\bf r})$$ (82)

where the exchange correlation potential is defined as a functional derivative of exchange correlation energy: