Computational Physics

112 Density functional theory

the radial wave function before integrating Poisson’s equation! If you take for the
normalisation ∫

dru^2 (r)=

∫

drr^2 R^2 (r)=1, (5.76)

you have already included a factor 4πinto the density (arising from the angular
integrations) and the factor 4πin Poisson’s equation drops out:

U′′(r)=−

u^2 (r)

r

. (5.77)

We shall use the normalisation (5.76) throughout this section.
The solution of Eq. (5.77) contains two integration constants which have to be
fixed by the boundary conditions. We takeU( 0 )=0 as the first boundary condition.
Elementary electrostatics then leads to the second condition

VH′(rmax)=

qmax

rmax^2

, (5.78)

whereqmaxis the electron charge contained in a sphere of radiusrmax:

qmax=

∫max

0

dru^2 (r). (5.79)

For largermax,qmax is the total electron charge. In that case we see from the
asymptotic form of (5.78), using the fact thatU(rmax)is now constant as a function
ofrmax), thatU(rmax)=qmax. When carrying out the numerical integration, we
take for the first starting conditionU( 0 )=0. The second starting condition, for
U(h), is not known at the beginning – we takeU(h)=h. As the solutionU(r)=αr
solves the homogeneous differential equation,U′′(r)=0, we can add this solution
to the numerical solution found, withαtaken such as to satisfy the end condition
U(rmax)=qmax, without violating the starting conditionU( 0 )=0.

programming exercise
Add an extra integration to your program which solves Eq. (5.77).
It is useful to check for correctness by using the hydrogen atom as an example.
The normalised ground state density (in the sense of (5.76)), found atE=−0.5 a.u.,
is 4e−^2 rand we must solve

U′′(r)=− 4 re−^2 r, (5.80)

with the boundary conditionsU( 0 )=0,U(∞)=1, so

U(r)=−(r+ 1 )e−^2 r+1. (5.81)

CheckCheck whether your program produces these results