(* N2 with Linnet type clouds:  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.21;
R2=0.98; R3=0.86; R4=1.18; R5=1.18;

nc=9;  (* number of clouds *) z=7.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; *)
(* cloud charges *)
ow={-2,-2,-2,-2,-2,-1,-1,-1,-1};
(* nucleus charges *)
cn={ z,z,0,0,0,0,0,0,0};
(* core contributions *)
Ekin = 2*(2.25*k1/R1^2) /. c;
vee=2*(3.0*sig1/R1);
vne=-2*(3.0*z/R1);
(* bonding clouds, Ekin and internal vee *)
Ekin = Ekin + 2.25*(2*k2/R2^2+k3/R3^2+2*k4/R4^2)/. c;
vee=vee+3.0*(2*sig2/R2+sig3/R3) /. c;
(* nuclear coordinates in terms of radii *)
d1=R1+R3;   p=d1; q=d1;  w4=pi+R4;
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 *)        
xw={d1,-d1,d2,-d2,0,-d1,     -d1,    d1,     d1,  };
yw={0,0,0,0,                0,w4,-w4, 0,0};
zw={0,0,0,0,                0,0,      0,   w4, -w4};
(* 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 && i≠j),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 && i≠j),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≠j),1/sq,0],{i,1,nc},{j,1,nc}];
(*MatrixForm[dcc]*)
(* sum of cloud-nucleus potential energies *)
vne=vne+cn.(ow.dcw);
(* sum of cloud-cloud potential energies, symm. matrix, hence 0.5 *)
vee=vee+0.5*ow.(ow.dw);
(* sum of nucleus-nucleus potential energies, symm matrix hence 0.5 *)
vnn =0.5*cn.(cn.dcc);

Epot=vne+vee+vnn  /. c;
func=Ekin+Epot /. c;
t = FindMinimum[func,{R1,0.21},{R2,0.99},{R3,0.86},{R4,1.17},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} (*  /. t[[2]]  *);

(* projection on xy-plane of molecule *)

plot0=Table[Graphics3D[{Darker[Blue,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,5}]    /. t[[2]]  ;
plot2=Table[Graphics3D[{Opacity[0.5],Lighter[Green,0.4],Sphere[{xw[[i]],yw[[i]],zw[[i]]},rr[[i]]]}],{i,6,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[{xw[[i]],zw[[i]]},rr[[i]]]}],{i,1,nc}]   /. t[[2]]  ;
Show[plot4,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]