(* H2 molecule stable equilibrium:  15.10.2012 ES d,r measured from origin *)
Clear[R,d,i,T,Vne,Vnn,Vee,r,r1,r2,func,Force];
sig=0.4165;
d = 0; (* H2 in stable equilibrium *)
T = (9/4)/R^2;
Vne = 2*(-(1.5-0.5*((r-d)/R)^2)/R-(1.5-0.5*((d+r)/R)^2)/R);
If[d+r > R, Vne=2*(-(1.5-0.5*((r-d)/R)^2)/R-1/(r+d))];
Vnn=1/(2*r);
Vee=(1.2-0.5*(2*d/R)^2+(3/16)*(2*d/R)^3-(1/160)*(2*d/R)^5)/R;
If[d > R, Vee=1/(2*d)];
func = T + Vne + Vnn + Vee; 
t=FindMinimum[func, {R, 1.1.36363}, {r,0.681799}]

N[t[[1]]]
u=t[[2]];
Vir=-(Vne+Vnn+Vee)/T /. u
Force = -1/(2*r)^2+2*(r/R)^3/(r)^2;
d-r /. u
Force /. u

(* Q=Table[t[[{i},2,1]],{i,1,81}];
P=Table[t[[{i},2,2]],{i,1,81}]; *)

Show[Graphics[{Circle[{-d /. u,0},R /. u],Circle[{d /. u,0},R /. u],Circle[{-d /. u,0},0.08],Circle[{d /. u,0},0.08], Disk[{-r /. u,0},0.08],Disk[{r /. u,0},0.08]}],{AspectRatio -> Automatic,PlotRange -> {{-2.5,2.5},{-2,2}}, Frame -> True}]