(* F2 with tetrahedric clouds: Linnet cannot represent F2 !  18.05.2013 *)
Clear[a,z,sig1,sig2,sig3,sig4,k1,k2,k3,k4,k5,nc,vee,vne,vnn,R1,R2,R3,R4,R5,
xw,yw,zw,xn,yn,zn,ow,cn,rr,d1,d2,d3,d4,d5,i,j,d,t,sq,Ekin,Epot,pi,p,cc,dv,dw,dcw,dcc,tw,tv,w4,w5,q];
bohr=0.529177; 
R1=0.1486;α=1.23;
R3=0.96; R2=1.166;

nc=9;  (* number of clouds *) z=9.0; pi=0.0;
sig1=0.325; sig2=0.325; sig3=sig2; (* screening const. from e-e interaction in doubly occ. clouds *)
g=Table[1.2,{nc}];
g[[1]]=0.8828; g[[2]]=0.8828; g[[3]]=1.33;(* parameters for kinetic energy of clouds; *)
(* cloud charges *)
ow={-2,-2,-2,-2,-2,-2,-2,-2,-2};
(* nucleus charges *)
cn={ z,z,0,0,0,0,0,0,0};
(* nuclear coordinates in terms of radii *)
d1=R1+R2;   d2=R1+R3;  d3=d1+d2*Cos[α]; d4=d2*Sin[α]; d5=d4*Sqrt[3]/2;
(*p=d1+d3; q=d1;  w4=pi+R4; w5=pi+R5;*)
xn={d1,-d1,0,0,0,0,0,0,0};
yn={0,0,0,0,0,0,0,0,0};
zn={0,0,0,0,0,0,0,0,0};
(* cloud coordinates in terms of radii *)        
xw={d1,-d1,0,d3,     d3,    d3,   -d3,-d3,-d3};
yw={0,0,0,d4,  -d4/2,-d4/2,-d4, d4/2,d4/2};
zw={0,0,0,0,      d5,-d5,      0,   d5, -d5};
(* cloud distances and reciprocal *)
tw=Table[{xw[[i]],yw[[i]],zw[[i]]},{i,1,nc}]; 
(*MatrixForm[tw]*)
dw=Table[sq=Norm[tw[[i]]-tw[[j]]];If[(sq>0 && j≤i),1/sq,0],{i,1,nc},{j,1,nc}];
(*MatrixForm[dw]*)
(* nuclear distances and cloud-nucleus reciprocal distances *)
tv=Table[{xn[[i]],yn[[i]],zn[[i]]},{i,1,nc}];
(*MatrixForm[tv]*)
dcw=Table[sq=Norm[tw[[i]]-tv[[j]]];If[(sq> 0.0 && j≤2),1/sq,0],{i,1,nc},{j,1,nc}];
(*MatrixForm[dcw]*)

(* reciprocal nuclear distances *)
dcc=Table[sq=Norm[tv[[i]]-tv[[j]]];If[(sq> 0.0 && i≤2 && j≤2),1/sq,0],{i,1,nc},{j,1,nc}];
(*MatrixForm[dcc]*)
(* sum of cloud-nucleus potential energies *)
vne=-2*(3.0*z/R1);
vne=vne+cn.(ow.dcw);

(*For[ i=1,i<nc+1,i++,
For[j=1,j<nc+1,j++,
d=Sqrt[(xw[[i]]-xn[[j]])^2+(yw[[i]]-yn[[j]])^2+(zw[[i]]-zn[[j]])^2];
If[ d>0,vne=vne+ow[[i]]*cn[[j]]/d]]];*)
(* sum of cloud-cloud potential energies, dw symm. matrix, hence 0.5 *)
vee=3.0*(2.0*sig1/R1+sig2/R2+6*sig3/R3)+ow.(ow.dw);
(*For[ i=1,i< nc ,i++,
For[ j=i+1,j<nc+1,j++,
d=Sqrt[(xw[[i]]-xw[[j]])^2+(yw[[i]]-yw[[j]])^2+(zw[[i]]-zw[[j]])^2];
If[ d>0 ,vee=vee+ow[[i]]*ow[[j]]/d]]];*)
(* sum of nucleus-nucleus potential energies, dcc symm matrix hence 0.5 *)
vnn =0.5*cn.(cn.dcc);
(*vnn=z*z/(2*d1);*)
rr={R1,R1,R2,R3,R3,R3,R3,R3,R3} (* /. t[[2]] *);
ekin=0.0;
For[ i=1,i< nc+1,i++,
      ekin=ekin-1.125*g[[i]]*ow[[i]]/rr[[i]]^2];
(*Ekin = 2.25*(2*k1/R1^2+k2/R2^2+6*k3/R3^2)/. c;*)
epot=vne+vee+vnn ;
func=ekin+epot 
(*t = FindMinimum[func,{R1,0.148},{R2,1.166},{R3,1.1},{α,1.23},Method ->"Newton", MaxIterations -> 500]*)

(* func *)
ekin    (* /. t[[2]] *)  
vne     (* /. t[[2]] *)   
vee     (* /. t[[2]] *)   
vnn     (* /. t[[2]] *)  
vir = -epot/ekin     (* /. t[[2]] *)  
2*d1*bohr     (* /. t[[2]] *)

(*rr={R1,R1,R2,R3,R3,R3,R3,R3,R3} *)(* /. t[[2]] *)(*;*)

(* projection on xy-plane of molecule *)

plot0=Table[Graphics3D[{Darker[Green,0.3],Sphere[{xw[[i]],yw[[i]],zw[[i]]},rr[[i]]]}],{i,1,2}]    (* /. t[[2]] *)  ;
plot1=Table[Graphics3D[{Sphere[{xw[[i]],yw[[i]],zw[[i]]},rr[[i]]]}],{i,3,3}]     (* /. t[[2]] *)  ;
plot2=Table[Graphics3D[{Opacity[0.5],Lighter[Green,0.4],Sphere[{xw[[i]],yw[[i]],zw[[i]]},rr[[i]]]}],{i,4,9}]    (* /. t[[2]] *)  ;
(*plot2a=Table[Graphics3D[{Opacity[0.5],Lighter[Cyan,0.3],Sphere[{xw[[i]],yw[[i]],zw[[i]]},rr[[i]]]}],{i,10,nc}]    (* /. t[[2]] *)  ;*)

Show[plot0,plot1,plot2,(*plot2a,*){AspectRatio -> Automatic,Axes -> True,PlotRange -> All, AxesLabel -> {xl, yl, zl},SphericalRegion->True}]

(* projection on xy-plane *)

plot3=Table[Graphics[{Circle[{xw[[i]],yw[[i]]},rr[[i]]]}],{i,1,nc}]    (* /. t[[2]] *)  ;Show[plot3,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]

(* projection on xz-plane *)
plot4=Table[Graphics[{Circle[{yw[[i]],zw[[i]]},rr[[i]]]}],{i,1,6}]    (* /. t[[2]] *)  ;
Show[plot4,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]