(* B2H6 with spherical bananas, but not TSM  28.06.2012 *)
Clear[z,nc,R1,R2,R3,R4,wi,p,vee,vne,vnn,σ,a1,a2,a3,b1,b2,b3,cs,ss,d1,d2,d3,d4,d5,pi,pb,i,j,d,t,q,r,ow,cn,c,w,n,dwr,dwnr,dnr,c2];

z=5.0;  nc=8; pi=0.0;
bohr=0.529177; rad=57.29578;
k={1.0,1.0,1.3,1.3,1.3,1.3,1.6,1.6};
σ={0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3};

(* C He shell, i=j each *)

dwr=Table[0,{i,1,nc},{j,1,nc}];
dnr=dwr; dwnr=dwr;
c2=Table[0,{i,1,nc}];
(*R4=R2;*)

ow={-2,-2,-2,-2,-2,-2,-2,-2};
cn={z,z,1,1,1,1,1,1};
c={R1,R1,R2,R2,R2,R2,R4,R4};
ec={0,0,p,p,p,p,pb,pb};
cs=Cos[wi];  ss=Sin[wi];
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};        
n=Table[{b1[[i]],b2[[i]],b3[[i]]},{i,1,nc}];
(* 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};
w=Table[{a1[[i]],a2[[i]],a3[[i]]},{i,1,nc}];

For[i = 1, i ≤ nc, i++,
     c2[[i]] = k[[i]]/c[[i]]^2;
     For[j = 1, j ≤ nc, j++,
         d = Norm[w[[i]] - w[[j]]];
         dwr[[i,j]] = If[i ≠ j, 1/d,1.5*(σ[[i]]/c[[i]])];
        d= Norm[n[[i]] - n[[j]]];
        dnr[[i,j]] = If[i ≠ j, 1/d, 0];
        d= Norm[w[[i]] - n[[j]]];
        dwnr[[i,j]] = If[i ≠ j, 1/d,(1.5 - 0.5*(ec[[i]]/c[[i]])^2)/c[[i]]]]];

Ekin = -(9/8)*ow . c2; 
vee = 0.5*ow . (ow . dwr); 
vnn = 0.5*cn . (cn . dnr); 
vne = cn . (ow . dwnr); 
Epot = vne + vee + vnn; 
vir=-Epot/Ekin;
func = Ekin + Epot;

t = NMinimize[{func,2*d1-3.3 ==0},{{R1,0.31,0.332},{R2,1.49,1.56},{R3,0.7,0.9},{R4,0.75,1.80},{p,0.59,0.61},{pb,0.78,0.80},{wi,1.14,1.15}}]

vee /. t[[2]]
vne /. t[[2]]
vnn /. t[[2]]
vir /. t[[2]]

2*d1*bohr /. t[[2]]                                                     (* B-B 1.753 *)

(R1+R2+p)*bohr /. t[[2]]                                     (* B-H2 1.19 *)

Sqrt[d1^2+(b3[[7]])^2]*bohr /. t[[2]]    (* B-H1  1.31 *)

2*wi*rad /. t[[2]]                                                         (* H2BH2 121.7 *)

2*ArcTan[d1/b3[[7]]]*rad /. t[[2]]             (* BH1B 83.94 *)

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},{-6,6}}, Frame -> True}]