(*CH2CH2 Ethene, with "orthogonal" σ,π LPs, 22.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,pnh,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,Nr,Cpi,NH,CH,CO,CN,CC,NM,OC,OO,OM,w1,w2,w3,pch,pnh];
ZC=6.0; sig1=0.3; sig4=0.3; sig5=0.3;
nc=11; nn=6; pi=0.0;
k1=1.034; k4=1.22; k5=1.1; k11=1.5;  
bohr=0.529177; rad=57.29578;

(* 6 Variables: Cr,CH,CC,Cpi,pch,w1 *)

(* Cr=0.26; Nr=0.22; Oa=0.186; NH=1.2; CH=1.25; CO=0.928; CN=1.04; OL=1.0;
NM=1.1; OC=1.4; OM=1.05; w1=1.9; w2=1.0; w3=1.0; pch=0.58;pnh=0.5; *)

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

(* bonding and lone pairs *)
Ekin = Ekin + 2.25*(4*k4/CH^2+k5/CC^2+k11/Cpi^2);

vee=vee+3.0*(4*sig4/CH+sig5/CC);

(*  1  2  3  4  5  6  7  8  9  10 11 *)
(*  C  C  Hc Hc Hc Hc CC L1 L1 L1 L1 *)
as={" ",C,C,H,H,H,H};    (* atom string *)  
rr={Cr,Cr,CH,CH,CH,CH,CC, Cpi ,Cpi, Cpi, Cpi};
oc={-2,-2,-2,-2,-2,-2,-2,-1/2,-1/2,-1/2,-1/2};
ch={ZC,ZC,1, 1, 1, 1};
csw1=Cos[w1/2];  ssw1=Sin[w1/2];

dccx=Cr+CC;             
dchx=dccx+(Cr+CH)*csw1;
dchy=(Cr+CH)*ssw1;  
pchx=dccx+(Cr+CH+pch)*csw1;
pchy=(Cr+CH+pch)*ssw1;

(* 6 nuclei: coord *)
xn={dccx,-dccx,pchx,pchx,-pchx,-pchx};
yn={0,0,pchy,-pchy,pchy,-pchy};
zn={0,0,0,0,0,0};

(* 11 clouds: coord *)        
xc={dccx,-dccx,dchx,dchx,-dchx,-dchx,0,dccx-0.145,dccx-0.145,-dccx+0.145,-dccx+0.145};
yc={0,0,dchy,-dchy,dchy,-dchy,0,0,0,0,0};
zc={0,0,0,0,0,0,0,Cpi,-Cpi,Cpi,-Cpi};

(* epot terms *)
vne=vne-4*(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 + oc[[i]]*oc[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])^2]]]
    
vee=vee-0.25/Cpi;

(*nn*)
For[i = 1, i < nn, i++,
  For[j = i+1, j < nn+1, j++,
    vnn = vnn + ch[[i]]*ch[[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 + oc[[i]]*ch[[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},{CH,1.25},{CC,0.98},{Cpi,1.31},{w1,2.02},{pch,0.54},{Method -> "Newton"}, {MaxIterations -> 800}]

Ekin /. t[[2]]
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
w1*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]]+0.06},{xn[[1]],yn[[1]]+0.06}}],
Line[{{xn[[2]],yn[[2]]-0.06},{xn[[1]],yn[[1]]-0.06}}],
Line[{{xn[[3]],yn[[3]]},{xn[[1]],yn[[1]]}}],
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]],yc[[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[[3]],zn[[3]]},{xn[[1]],zn[[1]]}}],
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}]