For questions and suggestions write to:

reibaretti2004@yahoo.com

References:

1. Density-Functional Theory of Atoms and Molecules (International Series of Monographs on Chemistry) by Robert G. Parr and Yang Weitao

2. http://en.wikipedia.org/wiki/Density_functional_theory

3.Introduction to Chemical Physics by J. C. Slater

4. http://en.wikipedia.org/wiki/Functional_derivative

The main goal of density functional theory (DFT)  is to find the ground state properties of atoms and molecules working with the electron density  ρ(r)   instead of the wave function Ψ (r₁ , r₂ ,   ..r_n ) . Just as the simplification is immense so are the various schemes and approximations proposed through the years.

The first connection between the electron density and the total kinetic energy was obtained in the free electron model of metals. This would later be incorporated in the first versions of DFT . A relation is established between the maximum kinetic energy and the electron density.

Consider electrons in a metal as free particles in a box of side L³  . The standing wave are constrained to vanish at the boundaries.

Thus         k_x L = n_x π ,    k_y L = n_y π ,       k_z L = n_z π .  Letting    p_x = h' k_x ,   p_y = h' k_y ,   p_x = h' k_x , ( h' =h/(2π)  )  ,

 one has a minimal volume  in momentum space equal to

                                    ∆p_x  ∆p_y   ∆p_z   = (π /h') ³(1/L³ ) (∆n_x  ∆n_y   ∆n_z )                        .  (1)

But the minimal change in the quantum numbers is , ∆n_x = 1, ∆n_y =1  ,  ∆n_z =1  . Hence the minimal volume in p space is

                                                               dτ_p   =   π ³/ L³  =  π ³/ V    ,  ( adopting h'=1)                 (2)

 where V=L³ is the box volume.  

Consider now  (1/8) - th  of the volume of  the spherical shell     (4π p² dp ) / 8. Dividing  by equation (2) one would have the number of states in the range p , p+dp

         number of states between , p and p+dp    ~  { (π p² dp ) / 2 } (V/ π ³ )           .                  (3)

 Each momentum state can be occupied by two electrons (in the ground state at T=0 absolute). Multiplying (3) by the factor 2  one gets

                                                                 dN (p)  =  (V/ π ² )   p² dp                        .    (4)

Introduce now  the kinetic energy expression ,      p² = ( 2 m ε ) = 2 ε   , dp = dε /(2ε )^1/2  , to obtain

                                                             dN(ε ) = V (2^1/2 /π ) ε ^1/2 dε    .                        (5)

Integrating (5) establishes a relation between the total number , N electrons , and the maximum value of ε .

 The number of electrons is

                                                         N = ( V 2^3/2 /3π² ) ε_max ^3/2                 .              (6)

Introducing the density  ρ = N/V  ( perhaps maximum density ) one gets

                                                    ε_max =  ( 3^2/3 π^4/3 / 2 ) ρ ^2/3                     .           (7)              

 The kinetic energy of the free electron gas is given by the integration

                                  T = ∫ ε  dN( ε ) = V (2^1/2 /π ) ∫  ε ^3/2 dε     ,    0 <  ε   <  ε_max  .    (8)

ε_max is called the Fermi energy level.  It corresponds to the electron with the shortest wavelength , maximum k value.

It follows from (7) that

                                                        T = { V 2^3/2 /(5π² ) } ε_max ^5/2                         .   (9)

Inserting (7)  in (9)  yields  the "kinetic energy density "  

                                  t_ρ   =   T/V  =  ( 3^5/3 π^4/3 /10 ) ρ ^5/3                         

                                                     = C_k ρ ^5/3                                            ,              (10)

where   C_k =  ( 3^5/3 π^4/3 /10 ) = 2.871234000188192 .

The first approximation in DFT to the kinetic energy in an atom was

                                                            KE =  ∫ C_k ρ ^5/3  d³ r                             .     (11)

In the Thomas Fermi model the total energy of the atom is

E_ρ  =  ∫ C_k ρ(r) ^5/3 d³ r  -Z ∫  (ρ(r )/ r ) d³ r  + (1/2) ∫ ∫ ρ(r )ρ(r' )d³ r  d³ r^' / abs(r-r')  .   (12)

A small check of formula (11) may be done using the well known  helium atom ground state. A first approximation is to take the  wave function as

                                                         Ψ_1s (r) = ( z³ /π )^1/2 exp(-zr)           ,              (13)

where z=2 ( the atomic number) . The electron  density is

                                                      ρ(r) = 2 { Ψ_1s (r) }²                           .          (14) 

Using (11) the kinetic energy estimate turns out to be   1.28 .

In the Hartree-Fock scheme ,the total energy is  -2.86 au  and by the virial theorem  KE =(1/2) /E/ = 2.86/2 =1.43 au .

Hence our functional density estimate is in error by nearly 10% .                           

MAXIMA code

1.278783881235112

By taking the functional derivative of (12)  ,see ref. 4 , an integral equation is found  for the density

                      δ E / δ ρ =  (5/3) C_k ρ(r) ^2/3 -  Z/r    +   ∫ { ρ(r') / abs(r-r') } d³ r'                  .   (15)

Since we are looking for the ground state , the minimum energy , it makes sense to equate δ E/ δ ρ =0 ,  just as in ordinary Calculus.

A non linear integral equation results

                           (5/3) C_k ρ(r) ^2/3 =   Z/r    -   ∫ { ρ(r') / abs(r-r') } d³ r'                          .   (16)  

A few steps have to be taken to get from (16) a manageable differential equation  known as the Thomas Fermi equation.

The right hand side of (16) is the electrostatic potential   Φ(r)   , and one writes

                               (5/3) C_k ρ(r) ^2/3 =  Φ(r)  =  Z/r    -   ∫ { ρ(r') / abs(r-r') } d³ r'            (17)

from which   one obtains  

                                                   ρ(r) =  { 3 Φ(r)/(5C_k ) }^3/2                              .         (18)

Recall that the laplacian operator is such that when applied to 1/r  ,

                     ∆² (1/ r )  = - 4πδ( r)   and   ∆² (1/ abs (r-r') )  = - 4πδ( r -r ').

Applying the laplacian operator  to  Φ(r) in (17) yields 

                                                  ∆²  Φ(r) = 4π  ρ(r) - 4π Z δ(r)                         ,      (19)

and using (18)

                                         ∆²  Φ(r) = 4π  { 3 Φ(r)/(5C_k ) }^3/2   - 4π Z δ(r)          ,(20)

or                             ∆²  Φ(r) = 4π  { Φ(r)/C  }^3/2   - 4π Z δ(r)             ,            (21)  

where  C= 5C_k/3 =  4.785390000313654.

To get rid of the nasty delta function introduce the function  screening  funtion  χ (r)  , through the definition

                                                              Φ(r) = (Z/r) χ (r)                             (22)

with boundary conditions χ (0) = 1    ,   χ ( ∞) = 0 .

From (21) ,

                           ∆²  Φ = ∂²Φ/∂ r²   + (2/r) ∂Φ /∂r    = (Z/r) χ ''   = (4π /C^3/2) (Z/r)^3/2   χ ^3/2    ,  (23)

and ths Thomas Fermi equation is reached,

Fig. 1.   The initial slope is about  -1.55  , in order to have    χ ( u_final = 5) = 0 .

Fig. 2 .   The Chi curve can be approximated by   χ (u) = 0.7772 exp( -.5019 u)  in order to simplify future calculations.

The energy functional (12) scales like Z^7/3 . That is    E_total = -constant Z^7/3 . Consider for example the attractive potential

  V_N = -Z ∫  (ρ(r )/ r ) d³ r .  If we go to the second  a-dimensional scale ,eq (25) , the Z dependence of ρ(r) is,

 ρ(r) =  { 3 Φ(r)/(5C_k ) }^3/2  ~ Z^3/2 L^-3/2 ~ Z^3/2 ( Z^-1/3 )^-3/2 ~ Z² .            

Also in the integrand

(1/r) d³r ~  L² ~ (Z^-1/3 )²   ~   Z^-2/3  . Hence

V_N ~  -Z (Z² ) Z^-2/3  ~ - Z^7/3  .

The same dependence on Z holds for the kinetic energy C_k ρ(r) ^5/3 d³ r  and the repulsive energy 

V_ee = (1/2) ∫ ∫ ρ(r )ρ(r' )d³ r  d³ r^' / abs(r-r') .

Calculation of E_ρ   

We start with T_ρ  =  ∫ C_k ρ(r) ^5/3 d³ r  . Using the relations between the density and the potential and  χ this reduces to

                                                T_ρ = 0.720253 Z^5/2 ∫  r^-1/2 χ(r) dr             .                    ( 27 ) 

Change now to the dimensionless variable x using the relation ( see eq .25-b)

     r = 0.885341 Z^-1/3 x   ,   dr = 0.885341 Z^-1/3 dx . whereby (27)  leads to

                              T_ρ = 0.677704 Z^7/3 ∫  x^-1/2 χ(x)^5/2 dx     ,              0 < x ∞ =      .    (28)

The integral   ∫  x^-1/2 χ(x)^5/2 dx  =  0.84261330803541  where we have used  the numerical fit  χ (u) = 0.7772 exp( -.5019 u) .

MAXIMA code

chi(u):= 0.7772*exp( -.5019*u);

 integrate(u^(-1/2)*chi(u)^(5/2),u,0,inf),numer;

0.84261330803541

Hence

                                                    T_ρ =  0.571042  Z^7/3                                             .  (29)

Going to the attractive potential term one gets

    V_N = -Z ∫  (ρ(r )/ r ) d³ r = -4π (3/5C_k )^3/2 Z^5/2  ∫  r^-1/2 χ(r) ^3/2 dr

                                        =  -1.200422 Z^5/2 ∫  r^-1/2 χ(r) ^3/2 dr                                      . (30)

Introducing the variable x , we get

 V_N = -1.200422 Z^5/2 ∫ (0.885341 Z^-1/3 x ) ^-1/2  χ(x) ^3/2  ( 0.885341 Z^-1/3 dx )

                                  =  -  1.12951   Z^7/3  ∫ x^-1/2 χ(x) ^3/2 dx      ,   0 < x < ∞   ,

                                 =   -  1.12951 *  Z^7/3 *  1.39965,

                                V_N   =    - 1.5809   Z^7/3                                                 .  (31) 

Multipltying(17) by ρ(r) and integrating gives

               (5/3)  C_k  ∫ ρ^5/3  d³ r  =   ∫  (Z/r) ρ(r)  d³ r  -  ∫ ρ(r )ρ(r' )d³ r  d³ r^' / abs(r-r')   , or what is equivalent ,

                                            (5/3) KE  = - V_Ne - 2 V_ee                          .        (32)

The repulsive energy is

                                          V_ee =  -(5/6)KE - (1/2) V_Ne                                . (33)

In terms of Z^7/3

                                           V_ee   =  0.31458 Z^7/3                              .       (34)

Addition of  (29, (31) and (34) gives the Thomas-Fermi formula for the ground state energy,

                                             E  =  -0.695278 Z^7/3                              .  (35)

In the literature (ref. 1 , page 113) one finds ,

                                           E =  - 0.7687 Z^7/3                                        . (36)

We are reminded that we  employed a smaller initial slope χ' (0) = 1.55 and after solving for χ  used an exponential fit assumed to be valid from 0 to infinity.   

Fig 3.  The approximation    E  =  -0.695278 Z^7/3   on the left ,  E (Hartree - Fock ) on the right.

integrate(u^(-1/2)*chi(u)^(3/2),u,0,inf),numer;

 1.399651445065135             

 (.570816-2.1412+.5949),numer;
(%o2) - 0.975484

FORTRAN code

  dimension rho(0:5000) ,phi(0:5000) ,chi(0:20000)
data nstep ,ck/15000,2.87123/
data niter/1/
etf(z)=-.7687*z**(7./3.)
c ehf has a variable factor multilying Z**(7./3.)
c ehf(z)= -.5678*z**(7./3.)
pi=2.*asin(1.)
c=5.*ck/3.
uzero=0.
ufinal=5.
du=(ufinal-uzero)/float(nstep)
delta=.40/float(niter)
d1chi=-1.555
do 50 it=1,niter
chi(0)=1.
chi(1)=chi(0)+d1chi*du
do 10 i=2,nstep
u=uzero+du*float(i)
chi(i)=2.*chi(i-1)-chi(i-2)+du**2*(1./(u-du)**.5)*chi(i-1)**1.5
if(chi(i).lt.0.)then
chifin=chi(i)
goto 20
endif
10 continue
chifin=chi(nstep)
20 print 100,u,d1chi,chifin
d1chi=d1chi+delta
50 continue
do 60 i=0,nstep,int(float(nstep)/60.)
u=uzero+du*float(i)
print 110,u,chi(i)
60 continue
100 format('u,d1chi,chi(nstep)=',3(3x,e10.3))
110 format('u,chi(i)=',2(3x,e10.3))
stop
end