(* B2H6 with spherical bananas, but not TSM  28.06.2012 *)
Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,
a1,a2,a3,b1,b2,b3,c1,c2,c3,cs,ss,d1,d2,d3,d4,d5,pi,pb,i,j,d,t];
z=5.0; sig1=0.3; sig2=0.3; sig4=0.3; nc=8; pi=0.0;
k1=1.0; k2=1.31; k4=1.66; bohr=0.529177; rad=57.29578;

(* C He shell, i=j each *)
Ekin = 2*(2.25*k1/R1^2);
vee=2*(3.0*sig1/R1);
vne=-2*(3.0*z/R1);
R4=R2;
(* bonding pairs *)
Ekin = Ekin + 2.25*(4*k2/R2^2+2*k4/R4^2);
vee=vee+3.0*(4*sig2/R2+2*sig4/R4);
(* Cloud charges, nuclei charges, cloud radii *)
c1={-2,-2,-2,-2,-2,-2,-2,-2};
c2={z,z,1,1,1,1,1,1};
c3={R1,R1,R2,R2,R2,R2,R4,R4};
cs=Cos[w];  ss=Sin[w];
d1=R1+R3; d2=d1+(R1+R2+p)*cs; d3=(R1+R2+p)*ss; 
d4=d1+(R1+R2)*cs;  d5=(R1+R2)*ss; d6=R4+pb;
(* Nuclear positions; 1,2,3 = x,y,z *)
b1={-d1,d1,-d2,-d2,d2,d2,0,0};
b2={0,0,d3,-d3,-d3,d3,0,0};
b3={0,0,0,0,0,0,d6,-d6};        
(* Cloud positions *)
a1={-d1,d1,-d4,-d4,d4,d4,0,0};
a2={0,0,d5,-d5,-d5,d5,0,0};
a3={0,0,0,0,0,0,pi+R4,-pi-R4};
vne=vne-4*(3-(p/R2)^2)/R2;
vne=vne-2*(3-(pb/R4)^2)/R4;
vnn = 0.0;
(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + c1[[i]]*c1[[j]]/Sqrt[(a1[[i]]-a1[[j]])^2+(a2[[i]]-a2[[j]])^2+(a3[[i]]-a3[[j]])^2]]]
(*nn*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vnn = vnn + c2[[i]]*c2[[j]]/Sqrt[(b1[[i]]-b1[[j]])^2+(b2[[i]]-b2[[j]])^2+(b3[[i]]-b3[[j]])^2]]]
(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nc+1, j++,
    If[i != j,
      vne = vne + c1[[i]]*c2[[j]]/Sqrt[(a1[[i]]-b1[[j]])^2+(a2[[i]]-b2[[j]])^2+(a3[[i]]-b3[[j]])^2]]]]

Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.315},{R2,1.49},{R3,0.85},{p,0.6},{pb,0.8},{w,1.15},{Method -> Automatic}, {MaxIterations -> 500}]

vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
2*d1*bohr /. t[[2]]
(R1+R2+p)*bohr /. t[[2]]
Sqrt[d1^2+(R4+pb)^2]*bohr /. t[[2]]
2*w*rad /. t[[2]]
2*ArcTan[(R4+pb)/d1]*rad /. t[[2]]

plot1=Graphics[{Circle[{a1[[1]],a2[[1]]},R1],      Circle[{a1[[2]],a2[[2]]},R1],Circle[{a1[[3]],a2[[3]]},R2],Circle[{a1[[4]],a2[[4]]},R2],Circle[{a1[[5]],a2[[5]]},R2],Circle[{a1[[6]],a2[[6]]},R2],Circle[{a1[[7]],a2[[7]]},R4],Circle[{a1[[8]],a2[[8]]},R4],Disk[{b1[[1]],b2[[1]]},0.08], Disk[{b1[[2]],b2[[2]]},0.08], Disk[{b1[[3]],b2[[3]]},0.08], Disk[{b1[[4]],b2[[4]]},0.08],Disk[{b1[[5]],b2[[5]]},0.08],Disk[{b1[[6]],b2[[6]]},0.08]} ] /. t[[2]]

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

plot2=Graphics[{Circle[{a1[[1]],a2[[1]]},R1],      Circle[{a1[[2]],a3[[2]]},R1],Circle[{a1[[3]],a3[[3]]},R2],Circle[{a1[[4]],a3[[4]]},R2],Circle[{a1[[5]],a3[[5]]},R2],Circle[{a1[[6]],a3[[6]]},R2],Circle[{a1[[7]],a3[[7]]},R4],Circle[{a1[[8]],a3[[8]]},R4],Disk[{b1[[1]],b3[[1]]},0.08], Disk[{b1[[2]],b3[[2]]},0.08], Disk[{b1[[3]],b3[[3]]},0.08], Disk[{b1[[4]],b3[[4]]},0.08],Disk[{b1[[5]],b3[[5]]},0.08],Disk[{b1[[6]],b3[[6]]},0.08],Disk[{b1[[7]],b3[[7]]},0.08],Disk[{b1[[8]],b3[[8]]},0.08]} ] /. t[[2]]

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