(* B2H6 not TSM, overlaps exchange corrected ES 08.06.2013 *)
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,r2,ow,cn,c,w,n,dwr,dwnr,dnr,c2];

z=5.0;  nc=8; pi=-0.6;
bohr=0.529177; rad=57.29578;
k={1.0135,1.0135,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=pi+R4+pb;
r=pi+R4; r2=2*r;
(* 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,r,-r};
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); 
delvee=(4/R4)*(1.2-0.5*(r2/R4)^2+(3/16)*(r2/R4)^3-(1/160)*(r2/R4)^5);
vee=vee-4/r2+delvee-(2/R4)*(-0.5*(r2/R4)^2+(3/16)*(r2/R4)^3-(1/160)*(r2/R4)^5);
vnn = 0.5*cn . (cn . dnr); 
vne = cn . (ow . dwnr); 
ratio=-(vne+2*vnn)/(vee-vnn);
Epot = vne + vee + vnn; 
vir=-Epot/Ekin;
func = Ekin + Epot;

t = NMinimize[{func,2*d1-3.3 ==0 && R1+R2+p == 2.249 && d1/d6==0.89946 && wi==1.062},{{R1,0.31,0.332},{R2,1.42,1.45},{R3,1.3,1.35},{R4,1.7,1.90},{p,0.49,0.61},{pb,0.58,0.65},{wi,1.05,1.1}}]

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

ratio /. 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+d6^2]*bohr /. t[[2]]    (* B-H1  1.31 *)

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

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

plot1=Table[Graphics[{Circle[{a1[[i]],a2[[i]]},c[[i]]]}],{i,1,nc}] /.t[[2]];
plot1a=Table[Graphics[{Disk[{b1[[i]],b2[[i]]},0.08]}],{i,1,6}] /.t[[2]];
plot1b=Table[Graphics[{Circle[{a1[[i]],a2[[i]]},0.08]}],{i,1,nc}] /.t[[2]];
Show[plot1,plot1a,plot1b,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, 
PlotRange -> {{-4,4},{-4,4}}, Frame -> True}]
xy projection

plot2=Table[Graphics[{Circle[{a1[[i]],a3[[i]]},c[[i]]]}],{i,1,nc}] /.t[[2]];
plot2a=Table[Graphics[{Disk[{b1[[i]],b3[[i]]},0.08]}],{i,1,nc}] /.t[[2]];
plot2b=Table[Graphics[{Circle[{a1[[i]],a3[[i]]},0.08]}],{i,1,nc}] /.t[[2]];
Show[plot2,plot2a,plot2b,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, 
PlotRange -> {{-4,4},{-4,4}}, Frame -> True}]
xz projection

plot3=Table[Graphics3D[{Opacity[0.4],Lighter[Gray,0.4],Sphere[{a1[[i]],a2[[i]],a3[[i]]},c[[i]]]}],{i,1,nc-2}] /.t[[2]];
plot3a=Table[Graphics3D[{Darker[Red,1.0],Sphere[{b1[[i]],b2[[i]],b3[[i]]},0.08]}],{i,1,nc}] /.t[[2]];
plot3aa=Table[Graphics3D[{Opacity[0.8],Lighter[White,0.2],Sphere[{b1[[i]],b2[[i]],b3[[i]]},0.3]}],{i,3,nc}] /.t[[2]];
plot3b=Table[Graphics3D[{Opacity[0.6],Lighter[Blue,0.7],Sphere[{a1[[i]],a2[[i]],a3[[i]]},0.9*c[[i]]]}],{i,nc-1,nc}] /.t[[2]];
plot3c=Table[Graphics3D[{Opacity[0.6],Lighter[Green,0.2],Sphere[{a1[[i]],a2[[i]],a3[[i]]},c[[i]]]}],{i,1,2}] /.t[[2]];
Plot9= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{b1[[1]],b2[[1]],b3[[1]]},{b1[[i]],b2[[i]],b3[[i]]}},0.1]}],{i,3,4}] /.t[[2]];
Plot91= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{b1[[2]],b2[[2]],b3[[2]]},{b1[[i]],b2[[i]],b3[[i]]}},0.1]}],{i,5,6}] /.t[[2]];
Plot92= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{b1[[1]],b2[[1]],b3[[1]]},{b1[[i]],b2[[i]],b3[[i]]}},0.1]}],{i,7,8}] /.t[[2]];
Plot93= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{b1[[2]],b2[[2]],b3[[2]]},{b1[[i]],b2[[i]],b3[[i]]}},0.1]}],{i,7,8}] /.t[[2]];
Plot94= Graphics3D[{Lighter[Green,0.2],Tube[{{b1[[2]],b2[[2]],b3[[2]]},{b1[[1]],b2[[1]],b3[[1]]}},0.1]}] /.t[[2]];

Show[plot3,plot3aa,plot3b,plot3c,Plot9,Plot91,Plot92,Plot93,Plot94,plot3a,{AspectRatio -> Automatic,(*Axes -> True,*)Boxed->False,SphericalRegion->True}]