(* LiH molecule with Li[1s2 2s], overlap, parametrized ES 31.12.2012 *)
(* Variables: P,Q,R,S *)
(* with exchange *)

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

Z=3.0;

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

Vne=-3.0*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);

Vee=3.0*s1/P+3.0*s2/Q+(4/Q)*(1.5-0.5*(R/Q)^2-0.3*(P/Q)^2);
Vex  = -1.2*(P/Q)^3/Q;
Vee = Vee+Vex;

Vnn=Z/(R+S);

(* Parameters in list c *)

(* c={k1 -> 1.0, k2-> 1.0, s1 -> 0.3, s2 -> 0.3}; *) 
c={k1 -> 1.01613950, k2 -> 0.65858068,s1-> 0.25742249,s2-> 0.22958899};

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

t=FindMinimum[func,{P,0.557},{Q,2.99},{R,1.49}, {S, 1.49}]

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[{-R,0},P],Disk[{-R,0},0.1],
			Circle[{0,0},Q],Disk[{S,0},0.1],
			{Thickness[0.012],Line[{{-R,0},{S,0}}]}}]/. u;

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

plot2=Graphics3D[{{Opacity[0.6],Sphere[{-R,0,0},P]},{Darker[Red,1],Sphere[{-R,0,0},0.1]},
	   {Opacity[0.4],Sphere[{0,0,0},Q]},
{Darker[Red],Sphere[{0,0,0},0.04],Sphere[{S,0,0},0.1]}}]/. u;
  Show[plot2, {AspectRatio->Automatic, Boxed -> False}]