(* O2 with:  12.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];
z=8.0;
(*pi=0.0;*) (* distance of bananas from C-C axis, pi =0 is regular tetraedric *)
nc=13;  (* number of clouds *)
sig1=0.3; sig2=0.3; sig3=sig2; (* screening const. from e-e interaction in doubly occ. clouds *)
k1=1.0; k2=1.0; k3=1.0; k4=1.0; k5=1.0;(* parameters for kinetic energy of clouds; k=1.0 Kimball's lowest value *)
bohr=0.529177; 
(*R1=0.18; R2=1.0; R3=0.96; R4=1.0; R5=1.0;*)

(* C He-shells *)
ow={-2,-2,-2,-2,-2,-1,-1,-1,-1,-0.5,-0.5,-0.5,-0.5};
cn={ z,z};
Ekin = 2*(2.25*k1/R1^2);
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+k5/R5^2);
vee=vee+3.0*(2*sig2/R2+sig3/R3);
d1=R1+R3;  
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,d1,  d1,-d1,-d1,d1,   d1,-d1,-d1};
yc={0,0,0,0,                0,R4,-R4,  0,      0,   R5,-R5,   0,     0};
zc={0,0,0,0,                0,0,      0,     R4,-R4,0,      0,      R5,-R5};
(* cc: sum of cloud-cloud potential energies *)
For[i = 1, i < nc, i++,
  For[j = i+1, j ≤ nc, j++,sq=Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])^2];
If[sq ≠ 0,vee = vee + ow[[i]]*ow[[j]]/sq]]]
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 = 1, i ≤ nc, i++,
  For[j = 1, j ≤ 2, j++,sq = Sqrt[(xc[[i]]-xn[[j]])^2+yc[[i]]^2+zc[[i]]^2];
    If[(i ≠ j && sq ≠ 0), vne = vne + ow[[i]]*cn[[j]]/sq]]]
Epot=vne+vee+vnn ;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.17},{R2,0.8},{R3,0.95},{R4,0.9},{R5,0.9},Method ->"Newton", MaxIterations -> 500]

(* func *)
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 *)

plot1=Table[Graphics[{Circle[{xc[[i]],yc[[i]]},rr[[i]]]}],{i,1,13}] /. t[[2]];
Show[plot1,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]

(* projection on xz-plane, perpendicular to molecular plane *)

plot2=plot1=Table[Graphics[{Circle[{xc[[i]],zc[[i]]},rr[[i]]]}],{i,1,13}] /. t[[2]];(*Graphics[{Circle[{xc[[1]],zc[[1]]},rr[[1]]],Circle[{xc[[2]],zc[[2]]},rr[[1]]],Circle[{xc[[3]],zc[[3]]},rr[[2]]],Circle[{xc[[4]],zc[[4]]},rr[[2]]],Circle[{xc[[5]],zc[[5]]},rr[[3]]],Circle[{xc[[6]],zc[[6]]},rr[[4]]],Circle[{xc[[7]],zc[[7]]},rr[[4]]],Circle[{xc[[8]],zc[[8]]},rr[[4]]],Circle[{xc[[9]],zc[[9]]},rr[[4]]],Circle[{xc[[10]],zc[[10]]},rr[[5]]],Circle[{xc[[11]],zc[[11]]},rr[[5]]],Circle[{xc[[12]],zc[[12]]},rr[[5]]],Circle[{xc[[13]],zc[[13]]},rr[[5]]]} ]  /. t[[2]];*)

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