(* Ethene, 21.01.2013 *)

Clear[Cr,CC,CH,Cpi,nc,nn,xn,yn,zn,xc,yc,zc,cx,chx,chy,
cchx,cchy,Ekin,vne,vee,vnn,Epot,w,csw,ssw,
k1,k2,k3,k4,s1,s2,s3,s4,pch];
nc=11;  nn=6; Z=6.0;
bohr=0.529177;rad=57.29578;
k1=1.032; k2=1.12; k3=1.23; k4=1.6;
s1=0.3;s2=0.3; s3=0.3;

(* 6 Variable: Cr,CC,CH,Cpi,pch,w *)

(*Cr=0.26; CC=1.0; CH=1.23; Cpi=1.4; w=1.02; pch=0.56;*)

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

(*bonding clouds and pairs *)

Ekin=Ekin+ 2.25*(k2/CC^2+4*k3/CH^2+k4/Cpi^2);
vee=vee+3.0*(s2/CC+4*s3/CH);
vne=vne-4*(3.0-(pch/CH)^2)/CH;

csw=Cos[w];    ssw=Sin[w];
cx=Cr+CC;
chx=   cx+(Cr+CH+pch)*csw;     chy=  (Cr+CH+pch)*ssw;
cchx=cx+(Cr+CH)*csw;                  cchy=(Cr+CH)*ssw;
ccpi=Cpi;

cn={Z,        Z,     1,      1,         1,        1};
xn={cx,-cx,chx,   chx,-chx,-chx};
yn={0,       0, chy,-chy,   chy,-chy};
zn={0,0,0,0,0,0};

cc={-2,  -2,      -2,        -2,        -2,       -2,-2,-1/2,-1/2,-1/2,-1/2};
xc={cx,-cx, cchx,   cchx,-cchx,-cchx,   0,       cx,       cx,    -cx,    -cx};
yc={  0,      0, cchy,-cchy,   cchy,-cchy,  0,          0,          0,         0,          0};
zc={  0,      0,         0,          0,          0,           0,  0,   ccpi,-ccpi,ccpi,-ccpi};

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

vee=vee-(1/4)/ccpi;

vnn=0.0;

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

(*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]])^2]]]];

(*{Method->"Newton"},{MaxIterations->800}]*)

Epot=vne+vee+vnn;
func=Ekin+Epot;
t=FindMinimum[func,{Cr,0.26},{CC,1.0},{CH,1.25},{Cpi,1.41},{w,1.05},{pch,0.55},{Method->"Newton",MaxIterations->800}]

(*t={func,{Cr->0.26,CC->1.0,CH->1.23,Cpi->1.4,w->1.02,pch->0.56}};*)

Ekin/.t[[2]]
vne/.t[[2]]
vee/.t[[2]]
vnn/.t[[2]]
-Epot/Ekin/.t[[2]]
2*w*rad/.t[[2]]
(2*Cr+2*CC)*bohr/.t[[2]]
(Cr+CH+pch)*bohr/.t[[2]]

plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},Cr],
Circle[{xc[[2]],yc[[2]]},Cr],
Circle[{xc[[3]],yc[[3]]},CH],
Circle[{xc[[4]],yc[[4]]},CH],
              Circle[{xc[[5]],yc[[5]]},CH],
Circle[{xc[[6]],yc[[6]]},CH],
Circle[{xc[[7]],yc[[7]]},CC],
Circle[{xc[[8]],yc[[8]]},Cpi],
Circle[{xc[[9]],yc[[9]]},Cpi],
Circle[{xc[[10]],yc[[10]]},Cpi],
Circle[{xc[[11]],yc[[11]]},Cpi],
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[[1]],yn[[1]]},{xn[[3]],yn[[3]]}}],
Line[{{xn[[1]],yn[[1]]},{xn[[4]],yn[[4]]}}],
Line[{{xn[[2]],yn[[2]]},{xn[[5]],yn[[5]]}}],
Line[{{xn[[2]],yn[[2]]},{xn[[6]],yn[[6]]}}]}} ] /. t[[2]];

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

plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},Cr],
Circle[{xc[[2]],zc[[2]]},Cr],
Circle[{xc[[3]],zc[[3]]},CH],
Circle[{xc[[4]],zc[[4]]},CH],
              Circle[{xc[[5]],zc[[5]]},CH],
Circle[{xc[[6]],zc[[6]]},CH],
Circle[{xc[[7]],zc[[7]]},CC],
Circle[{xc[[8]],zc[[8]]},Cpi],
Circle[{xc[[9]],zc[[9]]},Cpi],
Circle[{xc[[10]],zc[[10]]},Cpi],
Circle[{xc[[11]],zc[[11]]},Cpi],
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],
{Thickness[0.008],
Line[{{xn[[2]],zn[[2]]},{xn[[1]],zn[[1]]}}],
Line[{{xn[[1]],zn[[1]]},{xn[[3]],zn[[3]]}}],
Line[{{xn[[1]],zn[[1]]},{xn[[4]],zn[[4]]}}],
Line[{{xn[[2]],zn[[2]]},{xn[[5]],zn[[5]]}}],Line[{{xn[[2]],zn[[2]]},{xn[[6]],zn[[6]]}}]}} ] /. t[[2]];

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

plot3=Graphics3D[{
Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],      
Sphere[{xc[[2]],yc[[2]],zc[[2]]},Cr],
{Opacity[0.6],
Sphere[{xc[[3]],yc[[3]],zc[[3]]},CH],
Sphere[{xc[[4]],yc[[4]],zc[[4]]},CH],
Sphere[{xc[[5]],yc[[5]],zc[[5]]},CH],
Sphere[{xc[[6]],yc[[6]],zc[[6]]},CH],
Sphere[{xc[[7]],yc[[7]],zc[[7]]},CC],
{Opacity[0.4],{Lighter[Green],
Sphere[{xc[[8]],yc[[8]],zc[[8]]},Cpi],
Sphere[{xc[[9]],yc[[9]],zc[[9]]},Cpi],
Sphere[{xc[[10]],yc[[10]],zc[[10]]},Cpi],
Sphere[{xc[[11]],yc[[11]],zc[[11]]},Cpi]}},
{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]]-0.1,zn[[2]]},{xn[[1]],yn[[1]]-0.1,zn[[1]]}},0.05],
Cylinder[{{xn[[2]],yn[[2]]+0.1,zn[[2]]},{xn[[1]],yn[[1]]+0.1,zn[[1]]}},0.05],
Cylinder[{{xn[[3]],yn[[3]],zn[[3]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],
Cylinder[{{xn[[1]],yn[[1]],zn[[1]]},{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[[2]],yn[[2]],zn[[2]]}},0.05]}}}] /. t[[2]];

Show[plot3,{AspectRatio -> Automatic,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 Coordinates 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]]},Cr],
{Opacity[0.6],
{Opacity[0.4],{Lighter[Green],
Sphere[{xc[[8]],yc[[8]],zc[[8]]},Cpi],
Sphere[{xc[[9]],yc[[9]],zc[[9]]},Cpi],
Sphere[{xc[[10]],yc[[10]],zc[[10]]},Cpi],
Sphere[{xc[[11]],yc[[11]],zc[[11]]},Cpi]}},
{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[[3]],yn[[3]],zn[[3]]},{xn[[1]],yn[[1]],zn[[1]]}},0.05],
Cylinder[{{xn[[1]],yn[[1]],zn[[1]]},{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[[2]],yn[[2]]-0.1,zn[[2]]},{xn[[1]],yn[[1]]-0.1,zn[[1]]}},0.05],
Cylinder[{{xn[[2]],yn[[2]]+0.1,zn[[2]]},{xn[[1]],yn[[1]]+0.1,zn[[1]]}},0.05],Cylinder[{{xn[[6]],yn[[6]],zn[[6]]},{xn[[2]],yn[[2]],zn[[2]]}},0.06]}}}] /. t[[2]];

(* overlapped π system *)
Show[plot4,{AspectRatio -> Automatic,Boxed->False}]

plot6=Graphics3D[{Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],
Sphere[{xc[[2]],yc[[2]],zc[[2]]},Cr],
{Opacity[0.6],
Sphere[{xc[[3]],yc[[3]],zc[[3]]},CH],
Sphere[{xc[[4]],yc[[4]],zc[[4]]},CH],
Sphere[{xc[[5]],yc[[5]],zc[[5]]},CH],
Sphere[{xc[[6]],yc[[6]],zc[[6]]},CH],
Sphere[{xc[[7]],yc[[7]],zc[[7]]},CC]},
{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]]-0.1,zn[[2]]},{xn[[1]],yn[[1]]-0.1,zn[[1]]}},0.05],Cylinder[{{xn[[2]],yn[[2]]+0.1,zn[[2]]},{xn[[1]],yn[[1]]+0.1,zn[[1]]}},0.05],Cylinder[{{xn[[1]],yn[[1]],zn[[1]]},{xn[[4]],yn[[4]],zn[[4]]}},0.05],Cylinder[{{xn[[1]],yn[[1]],zn[[1]]},{xn[[3]],yn[[3]],zn[[3]]}},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[[2]],yn[[2]],zn[[2]]}},0.05]}}]/.t[[2]];

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