(* N2 with:  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];

z=7.0;
pi=0.35; (* distance of bananas from C-C axis, pi =0 is regular tetraedric *)
nc=9;  (* number of clouds *)
sig1=0.38; sig2=0.38; sig3=sig2; (* screening const. from e-e interaction in doubly occ. clouds *)
k1=0.9373; k2=1.14; k3=1.30; k4=1.25; (* parameters for kinetic energy of clouds; k=1.0 Kimball's lowest value *)
bohr=0.529177; 
(*R1=0.208; R2=1.1; R3=0.83; R4=1.2;*)
(* C He-shells *)
ow={-2,-2,-2,-2,-2,-1,-1,-1,-1};
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);
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};
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 *)        
xc={-d1,d1,d2,-d2,0,0,0,0,0};
yc={0,0,0,0,                0,pi+R4,-pi-R4,  0,      0};
zc={0,0,0,0,                0,0,      0,     pi+R4,-pi-R4};
(* 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]]]
(* 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*ow[[i]]*z/sq]]
Epot=vne+vee+vnn ;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.208},{R2,1.2},{R3,0.83},{R4,1.25},Method ->"Newton", MaxIterations -> 500]

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

plot0=Table[Graphics3D[{Darker[Blue,0.2],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[Orange,0.4],Sphere[{xc[[i]],yc[[i]],zc[[i]]},rr[[i]]]}],{i,6,nc}]   /. t[[2]]  ;

Show[plot0,plot1,plot2,{AspectRatio -> Automatic,Axes -> True,PlotRange -> All,SphericalRegion->True}]

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

plot3=Table[Graphics[{Circle[{xc[[i]],zc[[i]]},rr[[i]]]}],{i,1,nc}]   /. t[[2]]  ;

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