(* 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=7.0;
(*pi=0.0;*) (* distance of bananas from C-C axis, pi =0 is regular tetraedric *)
nc=9;  (* number of clouds *)
sig1=0.3; sig2=0.3; sig3=sig2; (* screening const. from e-e interaction in doubly occ. clouds *)
k1=0.8; k2=0.9; k3=0.9; 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 *)
(* N,N He shell*)
T = 2*2.25*k1/R1^2;
vee=2*3.0*sig1/R1;
vne=-2*3.0*z/R1;

a=Table[{0,0,0},{nc}];
b=Table[{0,0,0},{nc}];
c=Table[{0,0,0},{nc}];
(* bonding and lone pairs *)
T = T + 2.25*(2*k2/R2^2+k3/R3^2)+4*1.125*k4/R4^2);
vee=vee+3.0*(2*sig2/R2+sig3/R3);
c[[1,1]]=-2;      c[[1,2]]=z;  c[[1,3]]=R1;
c[[2,1]]=-2;      c[[2,2]]=z;  c[[2,3]]=R1;
c[[3,1]]=-2;      c[[3,2]]=0;  c[[3,3]]=R2;
c[[4,1]]=-2;      c[[4,2]]=0;  c[[4,3]]=R2;
c[[5,1]]=-2;      c[[5,2]]=0;  c[[5,3]]=R3;
c[[6,1]]=-1;      c[[6,2]]=0;  c[[6,3]]=R4;
c[[7,1]]=-1;      c[[7,2]]=0;  c[[7,3]]=R4;
c[[8,1]]=-1;      c[[8,2]]=0;  c[[8,3]]=R4;
c[[9,1]]=-1;      c[[9,2]]=0;  c[[9,3]]=R4;
d1=R1+R3;              b[[1,1]]=d1; 
d2=R1+R3+R1+R2;        b[[2,1]]=-d1;
                       a[[1,1]]=d1;
                       a[[2,1]]=-d1;
                       a[[3,1]]=d2;
		       a[[4,1]]=-d2;
                       a[[5,1]]=0.0;             
                       a[[6,1]]=d1;     a[[6,2]]=R4;
                       a[[7,1]]=d1;                 a[[7,3]]=-R4;
                       a[[8,1]]=-d1;    a[[8,2]]=R4;
		       a[[9,1]]=-d1;                 a[[9,3]]=-R4
For[ i=1,i<nc,i++,
  For[ j=i+1,j<=nc,j++,
    d=Sqrt[(a[[i,1]]-a[[j,1]])^2+(a[[i,2]]-a[[j,2]])^2+(a[[i,3]]-a[[j,3]])^2];
    If[ d<>0, vee=vee+c[[i,1]]*c[[j,1]]/d;
vnn=z*z/(2*d1);
For[ i=1,i<=nc,i++,
  For[ j=i+1,j<=nc,j++,
    If[ i<>j,
      d=Sqrt[(a[[i,1]]-b[[j,1]])^2+(a[[i,2]]-b[[j,2]]^2)+(a[[i,3]]-b[[j,3]])^2];
      If[ d<>0, vne=vne+c[[i,1]]*c[[j,2]]/d;
Ekin=T;
Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.21},{R2,0.8},{R3,0.9},{R4,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} /. t[[2]];

(* projection on xy-plane of molecule *)

plot1=Table[Graphics[{Circle[{xc[[i]],yc[[i]]},rr[[i]]]}],{i,1,9}] /. 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,9}] /. 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}]