Hartree-Self Consistant Field (SCF)

Application for He ground state

Coulomb  potential

kounter=1;
pots=Table[i,{i,1,10}];

vc[r_]=-Z/r /. Z->2;
veff[r_]=vc[r]+l(l+1)/2/r^2 /. l->0

Print["kounter=",kounter]

x0=10^(-12);
xend=4.0;
pots[[kounter]]=Plot[veff[r],{r,x0,xend},PlotRange->{-10,0}]

If[kounter==1,kounter2=kounter,kounter2=kounter-1];
Show[pots[[kounter]],pots[[kounter2]]]

eigcoulomb=-Z^2/n^2/2 /. {Z->2,n->1}

eig=eigcoulomb;
eig=-2.0
Clear[y];
Print["Eigenguess=",eig]
wrule=NDSolve[{-1/2 y''[x]+ veff[x]y[x]- eig y[x]==0,y[x0]==0.0,y'[x0]==0.1},y,{x,x0,xend},WorkingPrecision->6,MaxSteps->6000,Compiled->True][[1,1]]
"keep on guessing for the correct eigenvalue";

fun[x_]:=y[x] /. wrule

Plot[fun[x],{x,x0,xend},PlotRange->{-0.5,0.5}]

If we are satisfied with our numerical radial wave function we normalize it

norm1=NIntegrate[fun[x]^2,{x,x0,xend}];
nfun[x_]=fun[x]/Sqrt[norm1];

Lets construct the Hartree Self-Consistant Field for the 1s orbital

We  now  iterate

kounter = kounter + 1;

\(\(\(va[r_] := 
      If[r > x0, 1/r\ NIntegrate[\ nfun[z]^2, {z, x0, r}], 0] + 
	
        If[r < xend, NIntegrate[\ nfun[z]^2/z, {z, r, xend}], 0]\)\(
\)
    \)\)

veff[r_] := 
      vc[r] +  va[r] +  l (l + 1)/2/r^2  /.  l -> 0

"We now have a new effective potential";

Print["kounter=", kounter]

Lets manually loop back to the beginning and calculate a new wavefunction & eigenvalue 
and iterate again... and so on

Part II

The eigenvalue obtained above can be interpreted as a laboratory measurable quantity:
It's negative value is equal to the IONIZATION POTENTIAL of the atom. It is the
energy required to remove an electron from the atom (for example in hydrogen that is
just equal to -13.6 eV).   This fact is a consequence of  what is called  Koopman's Theorem.

In the ground state both electrons see the same Self consistant
potential that we have evaluated above

If we now use this potential as a mean-background field
we have the following Hamiltonian

H=((Subscript[p, 1]^2)/(2  m))+((Subscript[p, 2]^2)/(2  m))+U((Subscript[r, 1]))+U((Subscript[r, 2]))+H'

where U is the self consistant potential gotten above

H'=((e^2)/(Subscript[r, 12])) - U((Subscript[r, 1])) - U((Subscript[r, 2])) - ((Ze^2)/(Subscript[r, 1])) - ((Ze^2)/(Subscript[r, 2]))

ΔE=<φ|H'|φ>

H.W. Show that this energy shift is equal to the following integral

We find

ΔE=-∫d(Subscript[r, 1])∫d(Subscript[r, 1]) |φ((Subscript[r, 1]))\( | \^2\) |φ((Subscript[r, 2]))\( | \^2\) ((e^2)/(Subscript[r, 12]))= -I(1s,1s)

Note: This integral is just the direct integral we considered in the previous Mathematica  notebook,
except that the Coulomb radial wavefunctions are replaced by the wavefunctions obtained above by
the SCF procedure, thus

I1s1s = NIntegrate[nfun[x]^2 va[x], {x, x0, xend}]

The total energy of the He(1  (s^2))  ground state is equal to

E = ϵ(1s)+ϵ(1s)-I(1s,1s)

Energy = 2 eig  -  I1s1s