(* 2 H atoms  20.10.2012/19.03.2013 ES *)
Clear[R,d,i,T,Vne,Vnn,Vee,r,r1,r2,func,s,n,u];
sig=0.4165;  n=21;
d = (1.025-0.025*i)*R;
T = (9/4)/R^2;
Vne = -2*(1.5-0.5*((d-r)/R)^2)/R-2/(r+d);
If[d+r < R, Vne = -2*(1.5-0.5*((d-r)/R)^2)/R-2*(1.5-0.5*((d+r)/R)^2)/R];
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, Vner = -2/d]; *)
func = T + Vne + Vnn + Vee; 
t=Table[FindMinimum[func, {R, 0.996}, {r,0.5}, {Method->"Newton", MaxIterations -> 500}],{i,1,n}]

Q=Table[t[[{i},2,1]],{i,1,n}];
P=Table[t[[{i},2,2]],{i,1,n}];
B=Table[t[[i,1]],{i,1,n}]
Vir=Table[-(Vne+Vnn+Vee)/T /. t[[{i},2]],{i,1,n}]

Plot1=Table[Graphics[{Circle[{-d /. Q[[i]],0},R /. Q[[i]]],Circle[{d /. Q[[i]],0},R  /. Q[[i]]],Circle[{-d /. Q[[i]],0},0.06],Circle[{d /. Q[[i]],0},0.06], Disk[{-r /. P[[i]],0},0.06],Disk[{r /. P[[i]],0},0.06],Disk[{-d /. Q[[i]],10*(B[[i]]+0.96)-1.21},0.06],Disk[{d /. Q[[i]],10*(B[[i]]+0.96)-1.21},0.06]}],{i,1,n}];

Plot3=Table[Graphics[{Disk[{-d/.Q[[i]],10*(B[[i]]+0.96)-1.21},0.03],
Disk[{d/.Q[[i]],10*(B[[i]]+0.96)-1.21},0.03]}],{i,1,n}];


Show[Plot1,{AspectRatio -> Automatic,PlotRange -> {{-3,3},{-4,2}}, Frame -> True}]

sig=0.4165; u=21; s=0.525;
d = (s-0.025*i)*R;
T = (9/4)/R^2;
Vnn=1/(2*r);
Vee=1.0*(1.2-0.5*(2*d/R)^2+(3/16)*(2*d/R)^3-(1/160)*(2*d/R)^5)/R;
Vne = 2*(-(1.5-0.5*((d-r)/R)^2)/R-(1.5-0.5*((d+r)/R)^2)/R);

func = T + Vne + Vnn + Vee; 
t=Table[FindMinimum[func, {R,1.45}, {r,0.72}, {Method->"Newton", MaxIterations -> 500}], {i,1,u}]

Vir=Table[-(Vne+Vnn+Vee)/T /. t[[{i},2]],{i,1,u}]

Q=Table[t[[{i},2,1]],{i,1,u}];
P=Table[t[[{i},2,2]],{i,1,u}];
B=Table[t[[i,1]],{i,1,u}]

Plot2 = Table[Graphics[{Circle[{-(s-0.025*i)*R /. Q[[i]],0},R /. Q[[i]]],Circle[{(s-0.025*i)*R /. Q[[i]],0},R /. Q[[i]]],Circle[{-(s-0.025*i)*R /. Q[[i]],0},0.06],Circle[{(s-0.025*i)*R /. Q[[i]],0},0.06], Disk[{r /. P[[i]],0},0.06],Disk[{-r /. P[[i]],0},0.06],Disk[{(s-0.025*i)*R /. Q[[i]],10*(B[[i]]+0.96)-1.21},0.06],Disk[{-(s-0.025*i)*R /. Q[[i]],10*(B[[i]]+0.96)-1.21},0.06]}],{i,1,u}];

Show[Plot2,{AspectRatio -> Automatic,PlotRange -> {{-3,3},{-4,2}}, Frame -> True}]

Show[Plot1,Plot2,{AspectRatio->Automatic,PlotRange->{{-3,3},{-4,2}},Frame->True}]

Plot4= Table[Graphics[{Disk[{(s-0.025*i)*R /. Q[[i]],10*(B[[i]]+0.96)-1.21},0.03],Disk[{-(s-0.025*i)*R /. Q[[i]],10*(B[[i]]+
0.96)-1.21},0.03]}],{i,1,u}];
Plot5=Graphics[{Circle[{-(s-0.025*u)*R/.Q[[u]],0},R/.Q[[u]]],Circle[{(s-0.025*u)*R/.Q[[u]],0},0.06],Disk[{r/.P[[u]],0},0.06],Disk[{-r/.P[[u]],0},0.06],Disk[{(s-0.025*u)*R/.Q[[u]],10*(B[[u]]+0.96)-1.21},0.06],Disk[{-(s-0.025*u)*R/.Q[[u]],10*(B[[u]]+0.96)-1.21},0.06]}];


Show[Plot3,Plot4,Plot5,{AspectRatio -> Automatic,PlotRange -> {{-3,3},{-4,2}}, Frame -> True}]