8.9. Density Functional Theory (DFT)¶
8.9.1. Many Body Schrödinger Equation¶
We use (Hartree) atomic units in this whole section about DFT. We use the BornOppenheimer approximation, which says that the nuclei of the treated atoms are seen as fixed. A stationary electronic state (for electrons) is then described by a wave function fulfilling the manybody Schrödinger equation
where
is the kinetic term,
is the electronelectron interaction term and
is the interaction term between electrons and nuclei, where are positions of nuclei and the number of nucleons in each nucleus (we are using atomic units). So for one atomic calculation with the atom nucleus in the origin, we have just .
gives the probability density of measuring the first electron at the position , the second at , dots and the Nth electron at the position . The normalization is such that . The is antisymmetric, i.e. etc.
Integrating over the first electrons is the probability density that the th electron is at the position . Thus the probability density that any of the N electrons (i.e the first, or the second, or the third, dots, or the th) is at the position is called the particle (or number) density and is therefore given by:
(8.9.1.1)¶
Thus gives the number of particles in the region of integration . Obviously .
Note that the number density and potential in the Schroedinger equation is related to the electron charge density and electrostatic potential energy by:
where is the particle elementary charge, which for electrons is in atomic units. The amount of electronic charge in the region is given by:
The energy of the system is given by
(8.9.1.2)¶
where
(8.9.1.3)¶
It needs to be stressed, that generally is not a functional of alone, only the is. In the next section we show however, that if the is a ground state (of any system), then becomes a functional of .
8.9.2. The HohenbergKohn Theorem¶
The Schrödinger equation gives the map
where is the ground state. C is bijective (onetoone correspondence), because to every we can compute the corresponding from Schrödinger equation and two different and (differing by more than a constant) give two different , because if and gave the same , then by substracting
from
we would get , which is a contradiction with the assumption that and differ by more than a constant.
Similarly, from the ground state wavefunction we can compute the charge density giving rise to the map
which is also bijective, because to every we can compute from (8.9.1.1) and two different and give two different and , because different and give
adding these two inequalities together gives
which for gives , which is nonsense, so .
So we have proved that for a given ground state density (generated by a potential ) it is possible to calculate the corresponding ground state wavefunction , in other words, is a unique functional of :
so the ground state energy is also a functional of
We define an energy functional
(8.9.2.1)¶
where is any ground state wavefunction (generated by an arbitrary potential), that is, is a ground state density belonging to an arbitrary system. which is generated by the potential can then be expressed as
and for we have (from the Ritz principle)
and one has to minimize the functional :
(8.9.2.2)¶
The term
in (8.9.2.1) is universal in the sense that it doesn’t depend on . It can be proven [DFT], that is a functional of for degenerated ground states too, so (8.9.2.2) stays true as well.
The ground state densities in (8.9.2.1) and (8.9.2.2) are called purestate vrepresentable because they are the densities of (possible degenerate) ground state of the Hamiltonian with some local potential . One may ask a question if all possible functions are vrepresentable (this is called the vrepresentability problem). The question is relevant, because we need to know which functions to take into account in the minimization process (8.9.2.2). Even though not every function is vrepresentable [DFT], every density defined on a grid (finite of infinite) which is strictly positive, normalized and consistent with the Pauli principle is ensemble vrepresentable. Ensemble vrepresentation is just a simple generalization of the above, for details see [DFT].
The functional in (8.9.2.2) depends on the particle number , so in order to get , we need to solve the variational formulation
so
(8.9.2.3)¶
Let the be the solution of (8.9.2.3) with a particle number and the energy :
The Lagrangian multiplier is the exact chemical potential of the system
becuase
so
8.9.3. The KohnSham Equations¶
Consider an auxiliary system of noninteracting electrons (noninteracting gas):
the Schrödinger equation then becomes:
and the total energy is:
where
So:
The total energy is the sum of eigenvalues (energies of the individual independent particles) as expected. From the last equation it follows:
In other words, the kinetic energy of the noninteracting particles is equal to the sum of eigenvalues minus the potential energy coming from the total effective potential used to construct the single particle orbitals .
From (8.9.2.3) we get
(8.9.3.1)¶
Solution to this equation gives the density .
Now we want to express the energy in (8.9.1.2) using and for convenience, where is the classical electrostatic interaction energy of the charge distribution , defined using following relations  we start with a Poisson equation in atomic units
and substitute , and we use the fact that in atomic units:
or equivalently by expressing using the Green function:
(8.9.3.2)¶
and finally is related to using:
so we get:
Using the rules for functional differentiation, we can check that:
Using the above relations, we can see that
So from (8.9.2.1) we get
(8.9.3.3)¶
The rest of the energy is denoted by and it is called is the exchange and correlation energy functional. From (8.9.2.3)
From (8.9.3.2) we have
from (8.9.1.3) we get
we define
(8.9.3.4)¶
so we arrive at
(8.9.3.5)¶
Solution to this equation gives the density . Comparing (8.9.3.5) to (8.9.3.1) we see that if we choose
(8.9.3.6)¶
then . So we solve the KohnSham equations of this auxiliary noninteracting system
(8.9.3.7)¶
which yield the orbitals that reproduce the density of the original interacting system
(8.9.3.8)¶
The sum is taken over the lowest energies. Some of the can be degenerated, but it doesn’t matter  the index counts every eigenfunction including all the degenerated. In plain words, the trick is in realizing, that the ground state energy can be found by minimizing the energy functional (8.9.2.1) and in rewriting this functional into the form (8.9.3.3), which shows that the interacting system can be treated as a noninteracting one with a special potential.
8.9.4. The XC Term¶
The exchange and correlation functional
can always be written in the form
where the is called the XC energy density. The XC potential is defined as:
8.9.5. Total Energy¶
We already derived all the necessary things above, so we just summarize it here. The total energy is given by:
where
This is the correct, quadratically convergent expression for the total energy. We use the whole input potential and its associated eigenvalues to calculate the kinetic energy , this follows from the derivation of the expression for . Then we use the calculated charge density to express , and .
If one is not careful about the potential associated with the eigenvalues, i.e., confusing with , one gets a slowly converging formula for the total energy. By expanding using (8.9.3.6):
So is equal to:
And then the slowly converging form of total energy is:
The reason it is slowly converging is because the new formula for kinetic energy is mixing with , so it is not as precise (see above) and converges much slower with SCF iterations. Once selfconsistency has been achieved (i.e. ), the two expressions for total energy are equivalent.
8.9.6. XC Approximations¶
All the expressions above are exact (no approximation has been made so far). Unfortunately, no one knows exactly (yet). As such, various approximations for it exist.
LDA¶
The most simple approximation is the local density approximation (LDA), for which the xc energy density at is taken as that of a homogeneous electron gas (the nuclei are replaced by a uniform positively charged background, density ) with the same local density:
The xc potential defined by (8.9.3.4) is then
which in the LDA becomes
(8.9.6.1)¶
The xc energy density of the homogeneous gas can be computed exactly:
where the is the electron gas exchange term given by
the rest of is hidden in for which there doesn’t exist an analytic formula, but the correlation energies are known exactly from quantum Monte Carlo (QMC) calculations by Ceperley and Alder [pickett]. The energies were fitted by Vosko, Wilkes and Nussair (VWN) with and they got accurate results with errors less than in , which means that is virtually known exactly. VWN result:
where , , , , , (note that the value of is wrong in [pickett]), and is the electron gas parameter, which gives the mean distance between electrons (in atomic units):
The xc potential is then computed from (8.9.6.1):
Some people also use Perdew and Zunger formulas, but they give essentially the same results. The LDA, although very simple, is surprisingly successful. More sophisticated approximations exist, for example the generalized gradient approximation (GGA), which sometimes gives better results than the LDA, but is not perfect either. Other options include orbitaldependent (implicit) density functionals or a linear response type functionals, but this topic is still evolving. The conclusion is, that the LDA is a good approximation to start with, and only when we are not satisfied, we will have to try some more accurate and modern approximation.
RLDA¶
Relativistic corrections to the energydensity functional (RLDA) were proposed by MacDonald and Vosko:
where
We now calculate :
(8.9.6.2)¶
where the derivative can be evaluated as follows:
And in exactly the same manner:
So we can write
where
where we used the derivative of , which after a tedious, but straightforward differentiation is:
Plugging this back in, we get:
For we get , and as expected, because
Code:
>>> from sympy import limit, var, sqrt, log
>>> var("beta")
beta
>>> limit((beta*sqrt(1+beta**2)  log(beta+sqrt(1+beta**2)))/beta**2, beta, 0)
0
8.9.7. Radial DFT Problem¶
KohnSham Equations¶
For spherically symmetric potentials, we write all eigenfunctions as:
and we need to solve the following KohnSham equations:
With normalization:
For Schroedinger equation, the charge density is calculated by adding all “(n, l, m)” states together, counting each one twice (for spin up and spin down):
where we have introduced the occupation numbers by
Normalization of the charge density is:
So we can see, that it must hold:
where is the atomic number (number of electrons), and as such, are indeed the occupation numbers (integers). The factor is explicitly factored out, as it cancels with the spherical harmonics: assuming all states are occupied, this can be simplified to:
We can also use this machinery to prescribe “chemical occupation numbers”, that don’t necessarily correspond to the DFT ground state. For example for atom we get:
By summing all these , we get 92 as expected:
But this isn’t the DFT ground state, because some KS energies are skipped, for example there is only one state for , , but there are nine more states with the same energy — instead two more states are occupied in , , but those have higher energy. So this corresponds to excited DFT state, strictly speaking not physically valid in the DFT formalism, but in practice this approach is often used. One can also prescribe fractional occupation numbers (in the Dirac case).
Total Energy¶
The total energy is given by:
where
doing the integrals a bit we get:
We can also express everything using the charge density :
8.9.8. DFT As a Nonlinear Problem¶
The task is to find such a charge density , so that all the equations below hold (e.g. are selfconsistent):
This is a standard nonlinear problem, except that the Jacobian is dense, as shown below.
Reformulation¶
Let’s write everything in terms of explicitly:
Now we can write everything as just one (nonlinear) equation:
FE Discretization¶
The correspondig discrete problem has the form
where
Here is the vector of unknown coefficients for the th wavefunction . Our equation can then be written in the compact form
where with
Jacobian¶
The Jacobi matrix has the elements:
The only possible dense term is:
Now we can see that we have in there the following term:
which is dense in , as can be easily seen be fixing and writing
so for each there is some contribution from the integral for such where is nonzero, thus making the Jacobian dense.
8.9.9. ThomasFermiDirac Theory¶
There are two ways to derive equations for ThomasFermiDirac theory. One way is to start from grand potential and derive all equations from it. The other way is to start with low level equations and build our way up. Will start with the former.
Top Down Approach¶
We start with a grand potential for fermions:
The potential is the total potential that the electrons experience (it contains nuclear, Hartree, and XC terms) and is the Hartree energy:
For simplicity, we assume here that only contains the exchange of the homogeneous electron gas. For a general XC functional, the relation is nonlinear and one must simply numerically calculate the XC energy density and calculate the XC energy using:
In our case here, we have , which is only true for the exchange in homogeneous electron gas. Otherwise the relation is nonlinear. In the general case, the correction that must be applied is:
The density is a functional derivative with respect to :
By defining the function :
we can express the grand potential using as follows:
Now we can calculate the free energy:
where we used the fact that , i.e. the left hand side is a constant, thus the sum of the terms on the right hand side is also constant (even though the individual terms are not).
We can calculate the entropy as follows:
The total energy is then equal to:
From which we can see that the kinetic energy is equal to:
The relation between the total energy and free energy can be also written as:
But it gives the same result as we obtained above.
To determine the kinetic part of the free energy, we set all potentials equal to zero () and obtain:
If the potentials are zero, then the pressure can be calculated from:
If the potentials are not zero, then one can calculate the pressure using:
Summary:
where:
and is calculated as follows:
So can also be expressed as:
Bottom Up Approach¶
The other way to derive these equations is to use the following considerations. The number of states in a box of side is given by:
We use atomic units, so . The electronic particle density is:
(8.9.9.1)¶
where we used the relation for Fermi energy . The potential is the total potential that the electrons experience (it contains Hartree, nuclear and XC terms). At finite temperature we need to use the Fermi distribution and this generalizes to:
Now we use the relation and substitutions , to rewrite this using the FermiDirac Integral:
At low temperature () we have , and we obtain:
Identical with (8.9.9.1). We can see that the chemical potential becomes the Fermi energy in the limit . In the finitetemperature case, is determined from the normalization condition for the number of electrons :
The kinetic energy is
From the last formula it can be shown that the kinetic energy is equal to
The potential energy is equal to:
The internal energy is equal to:
The entropy is equal to:
where is the number of states at energy . We used the following calculation expressing one of the sums in terms of the kinetic energy:
where we used .
The free energy is equal to:
The grand potential is equal to:
We can now express the free energy functional as a function of the density:
8.9.10. Orbital Free Density Functional Theory¶
The orbitalfree electronic free energy is given by:
where the kinetic energy can be written in a few different equivalent ways as
where is a special function of one variable, composed of a FermiDirac Integral of order and its inverse of order :
the electronnuclei term has the form
The electronelectron (Hartree) term takes the form:
and the exchange and correlation functional is given by the PerdewZunger LDA:
is the (positive) electron density, is the (positive) nuclei density.
We minimize this free energy under the condition of particle conservation. The constrained functional is (we use from now on):
The variational solution is:
Or:
(8.9.10.1)¶
Finally we get:
(8.9.10.2)¶
The individual terms are:
and
and
and
All together the Hamiltonian is:
We can also introduce an artificial orbital as follows:
and minimize with respect to :
We will use tilde to denote functions in terms of . So . Using the relation
we obtain:
So the equation (8.9.10.1) gets multiplied by :
as well as the equation (8.9.10.2):
So the Hamiltonian expressed using and the Hamiltonian expressed using are related by .
Free Energy Minimization¶
For clarity, we will be using from equation (8.9.10.2) as our main quantity, but we will also write the final relations using for completeness.
We start with some initial guess for (it must be normalized as ). Let’s calculate :
We calculate the steepestdescent (SD) vector :
The conjugategradient (CG) vector is calculated as:
To satisfy the normalization constraint of , the CG vector is further orthogonalized to and normalized to (this step is one particular, but not the only way to impose the normalization constraint):
That is, now and . The new CG vector is then updated as usual in CG by , but then it must be normalized. As such, equivalently, it is updated by a linear combination of and :
such that it remains normalized:
So , are any real numbers satisfying the equation , whose parametric solution is , with :
where is determined by minimizing the free energy as a function of .
8.9.11. References¶
[DFT]  (1, 2, 3)

[pickett]  (1, 2)
