(*Fermi holes in H2- molecule, corresponding to the examples in the book *) 
R=2  (*internuclear distance*)
Clear[x,y]
a1s[x_,y_]:=(Sqrt[(1/π)]) Exp[-Sqrt[x^2+y^2]]  (*the 1s orbital on nucleus  a*)
Plot3D[a1s[x,y],{x,-1,3},{y,-2,2},PlotPoints->50]  (*the  1s orbital on nucleus  a*)
b1s[x_,y_]:=(Sqrt[(1/π)]) Exp[-Sqrt[(((x-R))^2)+(y^2)]]  (*the 1s orbital on nucleus b*)
Plot3D[b1s[x,y],{x,-1,3},{y,-2,2},PlotPoints->50](*the 1s orbital on nucleus b*)
S=(E^(-R)) (((R^2)/3)+R+1) (*overlap integral*)
fi1=((a1s[x,y]+b1s[x,y])/(Sqrt[2 (1+S)])) (*bonding orbital, doubly occupied*)
Plot3D[fi1,{x,-1,3},{y,-2,2},PlotPoints->50]  (*bonding orbital, doubly occupied*)
fi2=((a1s[x,y]-b1s[x,y])/(Sqrt[2 (1-S)]))(*antibonding orbital, singly occupied*)
Plot3D[fi2,{x,-1,3},{y,-2,2},PlotPoints->50]  (*antibonding orbital*)

fi10=((a1s[0,0]+b1s[0,0])/(Sqrt[2 (1+S)]))//N  (*the value of the bonding orbital at nucleus a*)
fi1R=((a1s[R,0]+b1s[R,0])/(Sqrt[2 (1+S)])) //N(*the value of the bonding orbital at nucleus b*)
fi20=((a1s[0,0]-b1s[0,0])/(Sqrt[2 (1-S)]))//N (*the value of the antibonding orbital at nucleus a*)
fi2R=((a1s[R,0]-b1s[R,0])/(Sqrt[2 (1-S)]))//N(*the value of the antibonding orbital at nucleus b*)
(*Example 1, p.598, We will compute the conditional probability of finding  electron 3 (its spin equal 1/2), when electron 1 of spin 1/2 is at the nucleus a,  
electron 2 is at b with the spin -1/2*)
Psi1=((fi10*fi1R*fi2-fi1R*fi20*fi1)/(Sqrt[6]))   (*Example 1, p.598*)
Psikw1=(Psi1^2)
Plot3D[Psikw1,{x,-1,3},{y,-2,2},PlotPoints->121,PlotRange->{0.0,.0032}]
(*Example 2, p.600: the samew but electron 1 is at nucleus  b*)
Psi2=((fi1R*fi1R*fi2-fi1R*fi2R*fi1)/(Sqrt[6]))  (*Example 2, p.600*)
Psikw2=(Psi2^2)
Plot3D[Psikw2,{x,-1,3},{y,-2,2},PlotPoints->121,PlotRange->{0.0,.0032}]
(*Example 3, p.601, a dilemma for electron 3: electron 1 at the center of the bond*)
fi1R5=((a1s[R/2,0]+b1s[R/2,0])/(Sqrt[2 (1+S)]))
Psi3=((fi1R5*fi1R*fi2-fi1R*fi1*0)/(Sqrt[6]))  (*Example 3, p.601*)
Plot3D[Psi3,{x,-1,3},{y,-2,2},PlotPoints->121] 
Plot3D[Psi3^2,{x,-1,3},{y,-2,2},PlotPoints->121,PlotRange->{0.0,.0005}]
fi1R6=((a1s[R/3,0]+b1s[R/3,0])/(Sqrt[2 (1+S)]))(*Example with electron 3 at 1/3R*)
fi2R6=((a1s[R/3,0]-b1s[R/3,0])/(Sqrt[2 (1-S)]))
Psi4=((fi1R6*fi1R*fi2-fi1R*fi2R6*fi1)/(Sqrt[6]))  (*Przykład 3 s.536*)
Plot3D[Psi4,{x,-1,3},{y,-2,2},PlotPoints->121]
Plot3D[Psi4^2,{x,-1,3},{y,-2,2},PlotPoints->121,PlotRange->{0.0,.0008}]
Show[
Plot3D[120*Psi4^2,{x,-1.,3.},{y,-2.,2.},PlotPoints->24,
PlotStyle->Directive[Opacity[0.8]],PlotRange->{0.0,.02}],
Graphics3D[{Sphere[{0,0,0},0.3]}],
Graphics3D[{Sphere[{2,0,0},0.3]}]
]

(13/(3 (E^2)))