(* F2C=O F2-formal, with "orthogonal" σ,π LPs, 19.01.2013 *)

Clear[ZC,ZN,ZO,sig1,sig2,sig3,sig4,sig5,sig6,sig7,sig8,sig9,k1,k2,k3,k4k5,k6,k7,k8,k9,
k10,k11,nc,w1,w2,w3,pFL,pch,vee,vne,vnn,xn,yn,zn,xc,yc,zc,c1,c2,c3,cs1,ss1,cs2,ss2,cs3,
ss3,pi,i,j,d,t];
Clear[Cr,Fr,Oa,CF,CO,OL,FL,Fpi,Cpi,Opi,Osg,w1,w2];
ZC=6.0; ZF=9.0; ZO=8.0; sig1=0.3; sig2=0.3; sig3=0.3; sig4=0.3; sig5=0.3; sig6=0.3;
sig7=0.3; sig8=0.3; sig9=0.3;
nc=22; nn=4; pi=0.0;
k1=0.98; k2=0.98; k3=0.98; k4=1.3; k5=1.65; k6=1.25; k7=1.24; k8=2.1; k9=1.1; k10=1.1; 
k11=0.75; k12=0.92;  
bohr=0.529177; rad=57.29578; pFL=0.0;

(* 13 Variables: Cr,Fr,Oa,CF,CO,OL,FL,Fpi,Cpi,Opi,Osg,w1,w2 *)

(* Cr=0.26; Fr=0.22; Oa=0.186; FL=1.2; CH=1.25; CO=0.928; CF=1.04; OL=1.0;
Fpi=1.1; Cpi=1.4; Osg=1.05; w1=1.9; w2=1.0; w3=1.0; pch=0.58;pFL=0.5; *)

(* C He shell *)
Ekin = 2.25*k1/Cr^2;
vee=3.0*sig1/Cr;
vne=-3.0*ZC/Cr;

(* 2* F He shell *)
Ekin = Ekin+2*2.25*k2/Fr^2;
vee=vee+2*3.0*sig2/Fr;
vne=vne-2*3.0*ZN/Fr;

(* O He shell *)
Ekin = Ekin+2.25*k3/Oa^2;
vee=vee+3.0*sig3/Oa;
vne=vne-3.0*ZO/Oa;

(* bonding and lone pairs *)
Ekin = Ekin + 2.25*(4*k4/FL^2+2*k5/CF^2+k6/CO^2+2*k8/Fpi^2+k9/OL^2+k10/Osg^2+k11/Cpi^2+k12/Opi^2);

vee=vee+3.0*(2*sig4/FL+sig5/CF+sig6/CO+sig7/CH+sig9/OL);

(*   1  2  3  4  5  6  7  8  9 10 11 12  13  14  15  16  17  18  19  20  21  22*)
(*   C  F  F  O  FL FL FL FL  B  B  B  L2 L1 L1 L1 L1  Bp   Bp   Bp   Bp L1 L1  *)
as={" ",C,F,F,O};    (* atom string *)  
rr={Cr,Fr,Fr,Oa,CF,CF,FL,FL,FL,FL,CO,OL,Fpi,Fpi,Fpi,Fpi, Cpi, Cpi, Opi, Opi,Osg,Osg};
cc={-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2, -1, -1, -1, -1,-1/2,-1/2,-1/2,-1/2, -1, -1};
cn={ZC,ZF,ZF,ZO};
csw1=Cos[w1];    ssw1=Sin[w1];
csw2=Cos[w2];    ssw2=Sin[w2];

dcx=0; dcy=0;
dcfx=(Cr+Fr+2*CF)*csw1;      dcfy=(Fr+CF)*ssw1;
dccfx=(Cr+CF)*csw1;          dccfy=(Cr+CF)*ssw1;
dcox=(Cr+Oa+2*CO);           dcoy=0;
dFLx=dCFx+(Fr+FL)*csw1;
dchx=(Cr+CH)*csw3;           dcCFx=Cr+CF; 
pFLx=dCFx+(Fr+FL+pFL)*csw1;  pchx=(Cr+CH+pch)*csw3;
      dFLy=(Fr+FL)*ssw1;       dchy=(Cr+CH)*ssw3;  
pFLy=(Fr+FL+pFL)*ssw1;       pchy=(Cr+CH+pch)*ssw3;
dccox=(Cr+CO)*csw2;          dccoy=(Cr+CO)*ssw2;      
dolx=(Cr+2*CO+2*Oa+OL)*csw2; doly=(Cr+2*CO+2*Oa+OL)*ssw2;
dOsg1x=dcox-(Oa+Osg)*ssw2;     dOsg2x=dcox+(Oa+Osg)*ssw2;
dOsg1y=dcoy+(Oa+Osg)*csw2;     dOsg2y=dcoy-(Oa+Osg)*csw2;

(* 4 nuclei:  *)
xn={dcx,-dcfx,-dcfx,dcox};
yn={dcy, dcfy,-dcfy,dcoy};
zn={0,0,0,0,0,0};

(* 22 clouds: cOpird *)        
xc={dcx,-dcfx,-dcfx,dcox,-dflx,-dflx,-dflx,-dflx,dcox,dolx,-dcfx,-dcfx,-dcfx,dcfx,0,0,dcox,dcox,dcox,dcox};
yc={0,dcfy,-dcfy,dcoy,dFLy,-dFLy,-dchy,0,dccoy,doly,0,0,0,0,dcoy,dcoy,dOsg1y,dOsg2y};
zc={0,0,0,0,0,0,0,0,0,Fpi,-Fpi,Cpi,-Cpi,Opi,-Opi,0,0};

(* epot terms *)
vne=vne-(3-(pch/CH)^2)/CH;
vnn = 0.0;

(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + cc[[i]]*cc[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])^2]]]

(*nn*)
For[i = 1, i < nn, i++,
  For[j = i+1, j < nn+1, j++,
    vnn = vnn + cn[[i]]*cn[[j]]/Sqrt[(xn[[i]]-xn[[j]])^2+(yn[[i]]-yn[[j]])^2+(zn[[i]]-zn[[j]])^2]]]

(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nn+1, j++,
    If[i ≠ j,
      vne = vne + cc[[i]]*cn[[j]]/Sqrt[(xc[[i]]-xn[[j]])^2+(yc[[i]]-yn[[j]])^2+(zc[[i]]-zn[[j]])^2]]]]

Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{Cr,0.26},{Fr,0.17},{Oa,0.19},{FL,0.96},{CH,1.22},{CO,0.91},{CF,0.92},{OL,0.90},{Fpi,0.70},{Cpi,1.42},{Opi,1.24},{Osg,0.84},{w1,2.41},{w2,1.05},{w3,1.03},{pch,0.57},{Method -> "Newton"}, {MaxIterations -> 800}]

{-212.3745773652353`,{Cr->0.25065062786145126`,Fr->0.16491191734154403`,Oa->0.1850880649193057`,FL->0.8984627963791239`,CH->1.2303713758008552`,CO->0.8774159467098952`,CF->1.0263884327126138`,OL->0.9240774122334838`,Fpi->0.8723369907371266`,Cpi->1.3151236824317798`,Opi->0.9807741650885555`,Osg->0.8619330161750159`,w1->2.394624355134884`,w2->1.0611598406083182`,w3->1.0280583469886226`,pch->0.5650391618101004`}}

Ekin /. t[[2]]
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
w1*rad /. t[[2]]
w2*rad /. t[[2]]
w3*rad /. t[[2]]
(Cr+Oa+2*CO)*bohr /. t[[2]]
(Cr+Fr+2*CF)*bohr /. t[[2]]
(Cr+CH+pch)*bohr /. t[[2]]
(Fr+FL+pFL)*bohr /. t[[2]]

plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},Cr],      Circle[{xc[[2]],yc[[2]]},Fr],Circle[{xc[[3]],yc[[3]]},Oa],Circle[{xc[[4]],yc[[4]]},FL],
Circle[{xc[[5]],yc[[5]]},FL],Circle[{xc[[6]],yc[[6]]},CH],Circle[{xc[[7]],yc[[7]]},CF],Circle[{xc[[8]],yc[[8]]},CO],Circle[{xc[[9]],yc[[9]]},OL],Circle[{xc[[10]],yc[[10]]},Fpi],Circle[{xc[[11]],yc[[11]]},Fpi],Circle[{xc[[12]],yc[[12]]},Cpi],Circle[{xc[[13]],yc[[13]]},Cpi],Circle[{xc[[14]],yc[[14]]},Opi],Circle[{xc[[15]],yc[[15]]},Opi],Circle[{xc[[16]],yc[[16]]},Osg],Circle[{xc[[17]],yc[[17]]},Osg],Disk[{xn[[1]],yn[[1]]},0.08], Disk[{xn[[2]],yn[[2]]},0.08], Disk[{xn[[3]],yn[[3]]},0.08],(*Disk[{xn[[4]],yn[[4]]},0.08],Disk[{xn[[5]],yn[[5]]},0.08],*) Disk[{xn[[6]],yn[[6]]},0.08],
{Thickness[0.008],
Line[{{xn[[2]],yn[[2]]},{xn[[1]],yn[[1]]}}],
Line[{{xn[[3]]+0.02,yn[[3]]-0.07},{xn[[1]]+0.03,yn[[1]]-0.07}}],
Line[{{xn[[3]]-0.02,yn[[3]]+0.06},{xn[[1]]-0.02,yn[[1]]+0.06}}],
(*Line[{{xn[[2]],yn[[2]]},{xn[[4]],yn[[4]]}}],
Line[{{xn[[2]],yn[[2]]},{xn[[5]],yn[[5]]}}],*)
Line[{{xn[[6]],yn[[6]]},{xn[[1]],yn[[1]]}}]}} ] /. t[[2]];

(* xy-plane *)
Show[plot1,{AspectRatio -> AutOsgatic,Axes -> True,GridLines -> AutOsgatic, PlotRange -> AutOsgatic, Frame -> True}]

plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},Cr],      Circle[{xc[[2]],zc[[2]]},Fr],Circle[{xc[[3]],zc[[3]]},Oa],Circle[{xc[[4]],zc[[4]]},FL],
Circle[{xc[[5]],zc[[5]]},FL],Circle[{xc[[6]],zc[[6]]},CH],Circle[{xc[[7]],zc[[7]]},CF],Circle[{xc[[8]],zc[[8]]},CO],Circle[{xc[[9]],zc[[9]]},OL],Circle[{xc[[10]],zc[[10]]},Fpi],Circle[{xc[[11]],zc[[11]]},Fpi],Circle[{xc[[12]],zc[[12]]},Cpi],Circle[{xc[[13]],zc[[13]]},Cpi],Circle[{xc[[14]],zc[[14]]},Opi],Circle[{xc[[15]],zc[[15]]},Opi],Circle[{xc[[16]],zc[[16]]},Osg],Circle[{xc[[17]],zc[[17]]},Osg],Disk[{xn[[1]],zn[[1]]},0.08], Disk[{xn[[2]],zn[[2]]},0.08], Disk[{xn[[3]],zn[[3]]},0.08],(*Disk[{xn[[4]],zn[[4]]},0.08], Disk[{xn[[5]],zn[[5]]},0.08],*) Disk[{xn[[6]],zn[[6]]},0.08]} ] /. t[[2]];

(* xz-plane *)
Show[plot2,{AspectRatio -> AutOsgatic,Axes -> True,GridLines -> AutOsgatic, PlotRange -> AutOsgatic, Frame -> True}]

plot3=Graphics3D[{
Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],      
Sphere[{xc[[2]],yc[[2]],zc[[2]]},Fr],
Sphere[{xc[[3]],yc[[3]],zc[[3]]},Oa],
{Opacity[0.6],
Sphere[{xc[[4]],yc[[4]],zc[[4]]},FL],
Sphere[{xc[[5]],yc[[5]],zc[[5]]},FL],
Sphere[{xc[[6]],yc[[6]],zc[[6]]},CH],
Sphere[{xc[[7]],yc[[7]],zc[[7]]},CF],
Sphere[{xc[[8]],yc[[8]],zc[[8]]},CO],
{Lighter[Blue],
Sphere[{xc[[9]],yc[[9]],zc[[9]]},OL]},
{Lighter[Blue,0.6],
Sphere[{xc[[16]],yc[[16]],zc[[16]]},Osg],
Sphere[{xc[[17]],yc[[17]],zc[[17]]},Osg]}},
{Opacity[0.4],{Lighter[Green],
Sphere[{xc[[10]],yc[[10]],zc[[10]]},Fpi],
Sphere[{xc[[11]],yc[[11]],zc[[11]]},Fpi],
Sphere[{xc[[12]],yc[[12]],zc[[12]]},Cpi],
Sphere[{xc[[13]],yc[[13]],zc[[13]]},Cpi],
Sphere[{xc[[14]],yc[[14]],zc[[14]]},Opi],
Sphere[{xc[[15]],yc[[15]],zc[[15]]},Opi]}},
{Darker[Red,1],
Sphere[{xn[[1]],yn[[1]],zn[[1]]},0.08],
Sphere[{xn[[2]],yn[[2]],zn[[2]]},0.08],
Sphere[{xn[[3]],yn[[3]],zn[[3]]},0.08],
(*Sphere[{xn[[4]],yn[[4]],zn[[4]]},0.08],
Sphere[{xn[[5]],yn[[5]],zn[[5]]},0.08],*)
Sphere[{xn[[6]],yn[[6]],zn[[6]]},0.08],
Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],
Cylinder[{{xn[[3]],yn[[3]]-0.14,zn[[3]]},{xn[[1]],yn[[1]]-0.14,zn[[1]]}},0.05],
Cylinder[{{xn[[3]],yn[[3]]+0.14,zn[[3]]},{xn[[1]],yn[[1]]+0.14,zn[[1]]}},0.05],
(*Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[4]],yn[[4]],zn[[4]]}},0.05],
Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[5]],yn[[5]],zn[[5]]}},0.05],*)
Cylinder[{{xn[[6]],yn[[6]],zn[[6]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05]}}] /. t[[2]];

Show[plot3,{AspectRatio -> AutOsgatic,Boxed->False}]

as;
x=bohr*xn /. t[[2]];
y=bohr*yn /. t[[2]];
z=bohr*zn /. t[[2]];
s=Table[{as[[l+1]],x[[l]],y[[l]],z[[l]]},{l,6}];
(* Cartesian COpirdinates of stationary structure *)
u=MatrixForm[s]

(* Distance matrix Å *)
Dij=Table[Sqrt[(x[[i]]-x[[j]])^2+(y[[i]]-y[[j]])^2],{i,6},{j,6}];
MatrixForm[Dij]

plot4=Graphics3D[{
Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],      
Sphere[{xc[[2]],yc[[2]],zc[[2]]},Fr],
Sphere[{xc[[3]],yc[[3]],zc[[3]]},Oa],
{Lighter[Blue],
Sphere[{xc[[9]],yc[[9]],zc[[9]]},OL]},
{Lighter[Blue,0.6],
Sphere[{xc[[16]],yc[[16]],zc[[16]]},Osg],
Sphere[{xc[[17]],yc[[17]],zc[[17]]},Osg]},
{Darker[Red,1],
Sphere[{xn[[1]],yn[[1]],zn[[1]]},0.08],
Sphere[{xn[[2]],yn[[2]],zn[[2]]},0.08],
Sphere[{xn[[3]],yn[[3]],zn[[3]]},0.08],
(*Sphere[{xn[[4]],yn[[4]],zn[[4]]},0.08],
Sphere[{xn[[5]],yn[[5]],zn[[5]]},0.08],*)
Sphere[{xn[[6]],yn[[6]],zn[[6]]},0.08],
Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],
Cylinder[{{xn[[3]],yn[[3]]-0.14,zn[[3]]},{xn[[1]],yn[[1]]-0.14,zn[[1]]}},0.05],
Cylinder[{{xn[[3]],yn[[3]]+0.14,zn[[3]]},{xn[[1]],yn[[1]]+0.14,zn[[1]]}},0.05],
(*Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[4]],yn[[4]],zn[[4]]}},0.05],
Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[5]],yn[[5]],zn[[5]]}},0.05],*)
Cylinder[{{xn[[6]],yn[[6]],zn[[6]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05]}}] /. t[[2]];

(* Two lone pairs of O atOsg; one of (Subscript[π, x]), the other of σ symmetry *)
Show[plot4,{AspectRatio -> AutOsgatic,Boxed->False}]

plot5=Graphics3D[{Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],Sphere[{xc[[2]],yc[[2]],zc[[2]]},Fr],Sphere[{xc[[3]],yc[[3]],zc[[3]]},Oa],{Opacity[0.4],{Lighter[Green,0.3],Sphere[{xc[[10]],yc[[10]],zc[[10]]},Fpi],Sphere[{xc[[11]],yc[[11]],zc[[11]]},Fpi],Sphere[{xc[[12]],yc[[12]],zc[[12]]},Cpi],Sphere[{xc[[13]],yc[[13]],zc[[13]]},Cpi],Sphere[{xc[[14]],yc[[14]],zc[[14]]},Opi],Sphere[{xc[[15]],yc[[15]],zc[[15]]},Opi]}},{Darker[Red,1],Sphere[{xn[[1]],yn[[1]],zn[[1]]},0.08],Sphere[{xn[[2]],yn[[2]],zn[[2]]},0.08],Sphere[{xn[[3]],yn[[3]],zn[[3]]},0.08],(*Sphere[{xn[[4]],yn[[4]],zn[[4]]},0.08],Sphere[{xn[[5]],yn[[5]],zn[[5]]},0.08],*)Sphere[{xn[[6]],yn[[6]],zn[[6]]},0.08],Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],Cylinder[{{xn[[3]],yn[[3]]-0.14,zn[[3]]},{xn[[1]],yn[[1]]-0.14,zn[[1]]}},0.05],Cylinder[{{xn[[3]],yn[[3]]+0.14,zn[[3]]},{xn[[1]],yn[[1]]+0.14,zn[[1]]}},0.05],(*Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[4]],yn[[4]],zn[[4]]}},0.05],Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[5]],yn[[5]],zn[[5]]}},0.05],*)Cylinder[{{xn[[6]],yn[[6]],zn[[6]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05]}}]/.t[[2]];

(* Overlapping π system, forcing the N atOsg and FL bonds into plane of all 6 atOsgs *)
Show[plot5,{AspectRatio->AutOsgatic,Boxed->False}]

plot6=Graphics3D[{Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],Sphere[{xc[[2]],yc[[2]],zc[[2]]},Fr],Sphere[{xc[[3]],yc[[3]],zc[[3]]},Oa],{Opacity[0.6],Sphere[{xc[[4]],yc[[4]],zc[[4]]},FL],Sphere[{xc[[5]],yc[[5]],zc[[5]]},FL],Sphere[{xc[[6]],yc[[6]],zc[[6]]},CH],Sphere[{xc[[7]],yc[[7]],zc[[7]]},CF],Sphere[{xc[[8]],yc[[8]],zc[[8]]},CO]},{Darker[Red,1],Sphere[{xn[[1]],yn[[1]],zn[[1]]},0.08],Sphere[{xn[[2]],yn[[2]],zn[[2]]},0.08],Sphere[{xn[[3]],yn[[3]],zn[[3]]},0.08],(*Sphere[{xn[[4]],yn[[4]],zn[[4]]},0.08],Sphere[{xn[[5]],yn[[5]],zn[[5]]},0.08],Sphere[{xn[[6]],yn[[6]],zn[[6]]},0.08],*)Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],Cylinder[{{xn[[3]],yn[[3]]-0.14,zn[[3]]},{xn[[1]],yn[[1]]-0.14,zn[[1]]}},0.05],Cylinder[{{xn[[3]],yn[[3]]+0.14,zn[[3]]},{xn[[1]],yn[[1]]+0.14,zn[[1]]}},0.05],(*Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[4]],yn[[4]],zn[[4]]}},0.05],Cylinder[{{xn[[2]],yn[[2]],zn[[2]]},{xn[[5]],yn[[5]],zn[[5]]}},0.05],*)Cylinder[{{xn[[6]],yn[[6]],zn[[6]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05]}}]/.t[[2]];

(* the sigma skeleton of the molecule *)
Show[plot6,{AspectRatio->AutOsgatic,Boxed->False}]