(*the pseudoberyllium atom. This is to show that any single Slater determinant does not describe the dynamic correlation, while their linear combination does*)
dzeta=1. (* orbital exponent, an arbitrary value*)
a2s=Sqrt[dzeta^5/(3*Pi)]*Sqrt[x^2+y^2]*Exp[-dzeta*Sqrt[x^2+y^2]] (*section z=0 of the normalized  2s Slater orbital taken in the form N r exp(-dzeta r)*)
a2px=dzeta*x*a2s (*section z=0 of the normalized 2px Slater orbital*)
a2s1=Sqrt[dzeta^5/(3*Pi)]*Exp[-dzeta]  (*the value of  2s  at (-1,0,0) or (1,0,0)*)
a2px1=a2s1*dzeta  (*value of  2px at (1,0,0)*)
rho0=0.5*a2s1^2*a2s^2 (*conditional density of  electron 2 in configuration 2s^2*)
rho1=0.5*a2px1^2*a2px^2 (*conditional density of  electron 2 in configuration 2px^2*)
kappa=-0.7 (*mixing of the two configurations as an arbitrary  negative number*)
rhominus=rho0+kappa^2*rho1-kappa*a2s1*a2px1*a2s*a2px (*conditional density of  electron 2 in configuration 2s^2+kappa*2px2, el1 at -1,0,0*)
rhoplus=rho0+kappa^2*rho1+kappa*a2s1*a2px1*a2s*a2px (*conditional density of  electron 2 in configuration 2s^2, el1 at +1,0,0*)
Show[
Plot3D[rho0,{x,-4.,4.},{y,-4.,4.},PlotPoints->41,PlotRange->{0.,0.007},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{-1.,0,0}, 0.05]]
] 
Show[
Plot3D[rho0,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,PlotRange->{0.,0.007},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{1.,0,0}, 0.05]]
](*no correlation in the configuration  2s^2*)
Show[
ContourPlot[rho0,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{-1.,0.},0.1]}]
]
Show[
ContourPlot[rho0,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{1.,0.},0.1]}]
]
(*no correlation in the configuration 2s^2*)
Show[
Plot3D[rho1,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,PlotRange->{0.,0.001},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{-1.,0,0}, 0.05]]
]
Show[
Plot3D[rho1,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,PlotRange->{0.,0.001},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{1.,0,0}, 0.05]]
]
(*no correlation in the configuration 2px^2*)
Show[
ContourPlot[rho1,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{-1.,0.},0.1]}]
]
Show[
ContourPlot[rho1,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{1.,0.},0.1]}]
]
(*no correlation in the configuration 2px^2*)
Show[
Plot3D[rhominus,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,PlotRange->{0.,0.007},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{-1.,0,0}, 0.05]]
]
Show[
ContourPlot[rhominus,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{-1.,0.},0.1]}]
]
Show[
Plot3D[rhoplus,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,PlotRange->{0.,0.007},AxesLabel->{"x","y","rho0"}],
Graphics3D[Sphere[{1.,0,0}, 0.05]]
]
Show[
ContourPlot[rhoplus,{x,-3.,3.},{y,-3.,3.},PlotPoints->41,Epilog->Inset[Graphics3D[Sphere[],ImageSize->110,Boxed->False],{0.505,Center}]],
Graphics[{Red,Disk[{1.,0.},0.1]}]
]
(*we see the correlation for the wave function  2s^2+kappa*2px^2*)