Relating the Energy Minimization to the Differential Equation

We begin by minimizing the energy given in Rydbergs (h/2π =1, 2m=1, e²=2, E_Hyd=-1) by

(1.1)

the functional derivative of the energy is

(1.2)

Note that the part coming from the variation of the denominator includes the energy to yield

(1.3)

Integration by parts utilizing δψ=0 on the boundaries yields

(1.4)
Start with

(1.5)

Express the Laplacian in spherical coordinates to yield

(1.6)

Let (1.7)

Then for V=V(r), independent of Θ and , the radial Schroedinger equation with u(r) subject to boundary conditions is

(1.8)

The Hartree Fock equations can be derived in the same way by assuming that

(1.9)

HFm.htm also requires an assumption that the Φ’s are orthogonal.

Leading to

(1.10)

The second potential term is non-local and thus much harder than the first.

Examining the equation at small r

a) Angular momentum > 0

For l > 0, the dominant term as r 0 is the l(l+1)/r² term which must be exactly cancelled by the second derivative of u.

Try the form so that the first two terms of (1.8) become

(1.11)

This implies so that for small r

(1.12)

b) Angular momentum = 0

In this case the dominant term in (1.8) is the 1/r singularity. Try the form yielding

The term in r 0 for small r, so that the first and third terms of (1.8) become

This implies that α=Z.

Examining the equation at large r

If we assume that for sufficiently large r, the potential and the angular momentum terms (L(L+1)/r²) are negligible, then equation (1.8) reduces to

(1.13)

For a solution of the form u=exp(-αr), so that this becomes

(1.14)

(1.15)

The - solution is the one that we want, but the + solution is also an exact solution to the differential equation. Suppose that at some step r₀, the region where the 1/r² term and the V(r) can be neglected has been reached. Then the numerical solution consist of exp(-αr) + εexp(+αr).
(1.16)

is solved at r₀ by

(1.17)

And also at all larger values of r. For sufficiently large values of r, the decaying exponential is negligible leaving only the exploding εexp(+αr) in the solution. The analytic result for u(r>r₀) and z(r>r₀) can thus be written as
(1.18)

In each step beyond r₀ the second term gets larger and larger. (1.19)

With ε=10^-7 and α=1 as for hydrogens’s ground state becomes 1 for r-r₀ = 16 Bohr radii. Which is well beyond the <r> of 1. For the first excited state α=.25 and the distance is 64 versus a <r> of 4, so things are not bad. But do not expect to integrate out 100 times the region of interest and have anything left over. Note that for scattering states E>0, the asymptotic solutions are of the form exp(i+-E^1/2) and thus ε times the unwanted solution stays of order ε and the problem above does not exist.

The shooting method (Press 16.1 old edition 17.1 new)

The Schroedinger equation for a bound state is a boundary value problem in which the object is to find a value for E such that u(0)=u(∞)=0. The E’s for which this is possible are eigenvalues. In a gross oversimplification their variation as a function of nuclear positions become the potentials in the Boltzman probabilities for finding the nuclei in these positions and hence give rise to the properties of the systems of nuclei and electrons. Shooting methods are best, except in the case of subtle boundary conditions such as “trying to maintain a decaying exponential in the presence of growing exponentials.” (Press 1st ed p. 581) Note that this is a problem for the Schoedinger equation. >>>>>> Beware of running into your authors in dark alleys, but “you might wish to follow the advice of your authors as notorious computer gunslingers: We always shoot first, and only then relax.”

(Press 1st ed p.581)<<<<<<<<<

The following code was used to produce the graphs below.

cpp\SHOOT.CPP

#include <stdio.h>

double shoot(double E)

// du/dr=x

// dz/dr=(veff(r)-E)*u This is f0 at beg of region f1 at end

{int n=1000,i,L=0,loop,iwrite=1;

double rmax=16,u0=0,z0=1,up,uc,zp,zc,f0,f1,h,r,veff,v1;

FILE *out;

if(iwrite == 1)out=fopen("u.out","w+");

h=rmax/n;

f0=0;

for (i=1;i <= n;++i){

r=h*i;

veff=-2/r+L*(L+1)/(r*r);

v1=veff-E;

zp=z0+h*f0;

up=u0;

loop=0;

predcorr:

f1=v1*up;

zc=z0+.5*h*(f0+f1);

uc=u0+.5*h*(z0+zp);

if(zp != zc || up != uc){

loop=loop+1;

up=uc;

zp=zc;

if(loop < 100)goto predcorr;}

if(iwrite == 1)fprintf(out," %lg %lg /n",r,up);

f0=f1;

z0=zp;

u0=up;

}

fclose(out);

printf(" from shoot E=%lg U(%lg) = %lg/n",E,rmax,u0);

return up;}

void main(void)

{double temp,Eig=7;

repeat:

printf(" enter a value for E -11 to exit/n");

scanf("%lf",&Eig);

if(Eig < -10)return;

temp=shoot(Eig);

printf("u(16) = %f /n",temp);

goto repeat;}

Note that the

shooting accuracy was achieved by human variation of the code