(* BeH2 molecule with Be[1s2 2s2], overlap, parametrized ES 02.11.2012/5.3.2013 *)
(* Variables: P,Q,R,S *)
(* with exchange *)

Clear[k1,k2,s1,s2, Z, P,Q,R,S];

Z=4.0;

T=2.25*k1/P^2+18.0*k2/Q^2;

Vne=-1.5*Z/P-(2.0*Z/Q)*(1.5-0.5*(R/Q)^2)-
(2/Q)*(1.5-0.5*(S/Q)^2)-2.0/(R+S);
Vne=2*Vne-4/(2*R+S);

Vee=1.5*s1/P+3.0*s2/Q+(2/Q)*(3.0-(R/Q)^2-0.6*(P/Q)^2);
Vex  = -1.2*(P/Q)^3/Q;
Vee = 2*(Vee+Vex)+(4/Q)*(1.2-0.5*(2*R/Q)^2+
(3/16)*(2*R/Q)^3-(1/160)*(2*R/Q)^5);

Vnn=2*Z/(R+S)+1/(2*(R+S));

(* Parameters in list c *)

(* c={k1 -> 1.0, k2-> 1.0, s1 -> 0.3, s2 -> 0.3}; *) 
c={k1 -> 1.02098865, k2 -> 0.64883648,s1-> 0.26743072,s2-> 0.21850553};

func=T+Vne+Vee+Vnn /. c;

t=FindMinimum[func,{P,0.4116},{Q,2.3},{R,1.29},{S, 1.27},Method->"Newton"]

u=t[[2]];

Vne /.u

Vee /.u  /. c

Vex /. u

Vnn /. u

-(Vne+Vee+Vnn)/T /. u /.c

(R+S)*0.529177 /. u

plot1=Graphics[{Circle[{0,0},P],Disk[{0,0},0.1],
			Circle[{-R,0},Q],Circle[{R,0},Q],Disk[{-(R+S),0},0.1],Disk[{R+S,0},0.1],
			{Thickness[0.012],Line[{{-(R+S),0},{R+S,0}}]}}]/. u;

Show[plot1, {AspectRatio->Automatic, PlotRange->{{-4,4},{-3,3}},Axes->True,GridLines->Automatic,Frame->True}]

plot2=Graphics3D[{{Darker[Cyan],Sphere[{0,0,0},P],Sphere[{0,0,0},0.12]},
			{Opacity[0.4],Sphere[{-R,0,0},Q],Sphere[{R,0,0},Q],Sphere[{-(R+S),0,0},0.1],Sphere[{R+S,0,0},0.08]},
			Tube[{{-(R+S),0,0},{R+S,0,0}},0.1]}]/. u;

Show[plot2, {AspectRatio->Automatic, PlotRange->Automatic,Boxed->False,SphericalRegion->True}]