(* F2 identical to F2_F2.PAS. ES 23.05.2013 after pas *)
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,x1,x2,x3,x4,Ekin,Epot,func,nc,z];
sq3=Sqrt[3];
nc=9; z=9.0;
a=Table[{0,0,0},{nc}];
b=a; c=a; 
g=Table[1.21283096,{nc}];  (* kin param *)
x1=1.1913455; x2=0.1593004; x3=0.8003616; x4=1.374017;  (* radii and angle *)
g[[1]]=2.35018989;g[[2]]=0.95875665;g[[6]]=g[[2]];
sig=Table[0.30820426,{nc}];
sig[[1]]=0.31403448; sig[[2]]:=0.30821367;sig[[6]]=sig[[2]];

ss=Sin[x4];cs=Cos[x4];
c[[1,1]]=-2.;c[[1,2]]=0;c[[1,3]]=x1;
(*F right*)
c[[2,1]]=-2.;c[[2,2]]=z;c[[2,3]]=x2;

c[[3,1]]=-2.;c[[3,2]]=0;c[[3,3]]=x3;
c[[4,1]]=-2.;c[[4,2]]=0;c[[4,3]]=x3;
c[[5,1]]=-2.;c[[5,2]]=0;c[[5,3]]=x3;
(*F left*)
c[[6,1]]=-2.;c[[6,2]]=z;c[[6,3]]=x2;

c[[7,1]]=-2.;c[[7,2]]=0;c[[7,3]]=x3;
c[[8,1]]=-2.;c[[8,2]]=0;c[[8,3]]=x3;
c[[9,1]]=-2.;c[[9,2]]=0;c[[9,3]]=x3;

(*coord with zero value not explicitely written!*)
(*nuclei*)
b[[2,1]]=x1+x2;
b[[6,1]]=-b[[2,1]];
(*clouds*)

a[[2,1]]=b[[2,1]];

a[[3,1]]=b[[2,1]]+(x2+x3)*cs;a[[3,2]]=(x2+x3)*ss;
a[[4,1]]=a[[3,1]];a[[4,2]]=-0.5*a[[3,2]];a[[4,3]]=0.5*sq3*a[[3,2]];
a[[5,1]]=a[[3,1]];a[[5,2]]=a[[4,2]];a[[5,3]]=-a[[4,3]];

a[[6,1]]=-b[[2,1]];
a[[7,1]]=-a[[3,1]];a[[7,2]]=-a[[3,2]];
a[[8,1]]=-a[[3,1]];a[[8,2]]=-a[[4,2]];a[[8,3]]=a[[4,3]];
a[[9,1]]=-a[[3,1]];a[[9,2]]=-a[[4,2]];a[[9,3]]=a[[5,3]];

vnn=z*z/(2*(x1+x2));   vee=0.0;vne=0.0;ekin=0.0;
(*WI/WJ EXCL.I=J*)
For[ i=1,i< nc ,i++,
For[ j=i+1,j<nc+1,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]]]

(*KJ/WI I=J*)
(*vne=1.5*c[[2,1]]*c[[2,2]]/c[[2,3]];
vne=vne+1.5*c[[6,1]]*c[[6,2]]/c[[6,3]];*)
vne=-2*3.0*z/x2;
(*KJ/WI I<>J inc;*)
For[ i=1,i<nc+1,i++,
For[j=1,j<nc+1,j++,
d=Sqrt[(a[[i,1]]-b[[j,1]])^2+a[[i,2]]^2+a[[i,3]]^2];
If[ d>0,vne=vne+c[[i,1]]*c[[j,2]]/d]]];
(*d=Sqrt[(a[[i,1]]-b[[6,1]])^2+a[[i,2]]^2+a[[i,3]]^2];
If[ d>0 ,vne=vne+c[[i,1]]*z/d]];*)
For[ i=1,i< nc+1,i++,
      ekin=ekin-1.125*g[[i]]*c[[i,1]]/c[[i,3]]^2];
For[ i=1,i< nc+1,i++,
vee=vee+3.0*sig[[i]]/c[[i,3]]];
epot=vne+vee+vnn;
func=ekin+epot

(*t=FindMinimum[func,{x1,1.14,1.167},{x2,0.14,0.156},{x3,0.9,1.0},{=,1.21,1.3},Method->Automatic,MaxIterations->500];*)

ekin
vne
vee
vnn
vir =-epot/ekin
2*(x1+x2)*0.529177

rr={x1,x2,x3,x3,x3,x2,x3,x3,x3}  (* /. t[[2]] *) ;

plot0=Table[Graphics3D[Sphere[{a[[i,1]],a[[i,2]],a[[i,3]]},rr[[i]]]],{i,1,1}]    (* /. t[[2]] *)  ;
plot1=Table[Graphics3D[{Darker[Green,0.3],Sphere[{a[[i,1]],a[[i,2]],a[[i,3]]},rr[[i]]]}],{i,2,2}]     (* /. t[[2]] *)  ;
plot2=Table[Graphics3D[{Opacity[0.5],Lighter[Green,0.4],Sphere[{a[[i,1]],a[[i,2]],a[[i,3]]},rr[[i]]]}],{i,3,5}]    (* /. t[[2]] *)  ;
plot2a=Table[Graphics3D[{Darker[Green,0.3],Sphere[{a[[i,1]],a[[i,2]],a[[i,3]]},rr[[i]]]}],{i,6,6}]     (* /. t[[2]] *)  ;
plot2b=Table[Graphics3D[{Opacity[0.5],Lighter[Green,0.4],Sphere[{a[[i,1]],a[[i,2]],a[[i,3]]},rr[[i]]]}],{i,7,nc}]    (* /. t[[2]] *)  ;

Show[plot0,plot1,plot2,plot2a,plot2b,{AspectRatio -> Automatic,Axes -> True,PlotRange -> All, AxesLabel -> {xl, yl, zl},SphericalRegion->True}]

plot3=Table[Graphics[{Circle[{a[[i,1]],a[[i,2]]},rr[[i]]]}],{i,1,nc}]    (* /. t[[2]] *)  ;Show[plot3,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]

plot4=Table[Graphics[{Circle[{a[[i,2]],a[[i,3]]},rr[[i]]]}],{i,1,5}]    (* /. t[[2]] *)  ;
Show[plot4,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-4,4},{-3,3}}, Frame -> True}]