(*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.5; sig5=0.3;
nc=11; nn=6; pi=0.0;
k1=1.02; k4=1.058; k5=1.12; k11=1.12;
bohr=0.529177; rad=57.29578;
(* 6 Variables: Cr,CH,CC,Cpi,pch,w1 *)
(* 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;
dcpi=0.5*Cr+Cpi;
dpix=dccx-0.15;
(* 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,dpix,dpix,-dpix,-dpix};
yc={0,0,dchy,-dchy,dchy,-dchy,0,0,0,0,0};
zc={0,0,0,0,0,0,0,dcpi,-dcpi,dcpi,-dcpi};
(* 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.99},{Cpi,1.28},{w1,2.04},{pch,0.55},{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}]
