(* O2 with Linnet's clouds  13.05.2013 *)
Clear[a,z,sig1,sig2,sig3,sig4,k1,k2,k3,k4,k5,nc,vee,vne,vnn,R1,R2,R3,R4,R5,
xc,yc,zc,xn,yn,zn,ow,cn,rr,d1,d2,d3,d4,d5,i,j,d,t,sq,Ekin,Epot,pi,p];
(* distance of bananas from C-C axis, pi =0 is regular tetraedric *)

nc=13;  (* number of clouds *) z=9.0; pi=0.0;
sig1=0.3; sig2=0.3; sig3=sig2; (* screening const. from e-e interaction in doubly occ. clouds *)
c={k1->0.8828, k2->1.4, k3->1.35, k4->1.6, k5->1.6};(* parameters for kinetic energy of clouds; *)(* k=1.0 Kimball's lowest value *)
bohr=0.529177; 
R1=0.1486;
R2=0.98; R3=1.166;R4=1.18; R5=1.18; 
(* C He-shells *)
ow={-2,-2,-2,-2,-2,-1,-1,-1,-1,-1,-1,-1,-1};
cn={ z,z};
Ekin = 2*(2.25*k1/R1^2) /. c;
vee=2*(3.0*sig1/R1);
vne=-2*(3.0*z/R1);
(* bonding pairs *)
Ekin = Ekin + 2.25*(2*k2/R2^2+k3/R3^2+2*k4/R4^2+2*k5/R5^2)/. c;
vee=vee+3.0*(2*sig2/R2+sig3/R3) /. c;
d1=R1+R3;   p=d1; q=d1;
d2=d1+R1+R2;
xn={-d1,d1,0,0,0,0,0,0,0,0,0,0,0};
yn={0,0,0,0,0,0,0,0,0,0,0,0,0};
zn={0,0,0,0,0,0,0,0,0,0,0,0,0};
(* cloud coordinates in terms of radii *)        
xc={d1,-d1,d2,-d2,0,-p,     -p,    p,     p,   q,q,-q,-q};
yc={0,0,0,0,                0,pi+R4,-pi-R4, 0,0, pi+R5,    -pi-R5,  0,      0};
zc={0,0,0,0,                0,0,      0,   pi+R4, -pi-R4, 0,0,   pi+R5,-pi-R5};
(* cc: sum of cloud-cloud potential energies *)
For[i = 1, i < nc, i++,
  For[j = i+1, j ≤ nc, j++,
vee = vee + ow[[i]]*ow[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])^2]]]
(*vee=vee-2*0.25/(2*R5);*)
(* nn: sum of nuclei-nuclei potential energies *)
vnn =z^2/(2*d1);
(* cn: sum of cloud-nuclei potential energies *)
For[i = 2, i ≤ nc, i++,
  sq = Sqrt[(xc[[i]]-xn[[1]])^2+yc[[i]]^2+zc[[i]]^2];
    If[ sq ≠ 0,vne = vne + 2*z*ow[[i]]/sq]]
Epot=vne+vee+vnn  /. c;
func=Ekin+Epot /. c
(*t = FindMinimum[func,{R1,0.19},{R2,0.97},{R3,0.95},{R4,1.15},{R5,1.15},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,R2,R3,R4,R4,R4,R4,R5,R5,R5,R5} (*(* /. t[[2]] *)*);

(* projection on xy-plane of molecule *)

plot0=Table[Graphics3D[{Darker[Green,0.3],Sphere[{xc[[i]],yc[[i]],zc[[i]]},rr[[i]]]}],{i,1,2}] (* /. t[[2]] *);
plot1=Table[Graphics3D[{Sphere[{xc[[i]],yc[[i]],zc[[i]]},rr[[i]]]}],{i,3,5}]  (* /. t[[2]] *);
plot2=Table[Graphics3D[{Opacity[0.5],Lighter[Green,0.4],Sphere[{xc[[i]],yc[[i]],zc[[i]]},rr[[i]]]}],{i,6,9}] (* /. t[[2]] *);
plot2a=Table[Graphics3D[{Opacity[0.5],Lighter[Cyan,0.3],Sphere[{xc[[i]],yc[[i]],zc[[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[{xc[[i]],yc[[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[{xc[[i]],zc[[i]]},rr[[i]]]}],{i,1,nc}] (* /. t[[2]] *);
Show[plot4,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]