(*NH2COH Formamide, with "orthogonal" σ,π LPs, 13.01.2013, not final optimization: slight overlaps *)

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,Oa,NH,CH,CO,CN,OL,NM,OC,OO,OM,w1,w2,w3,pch,pnh];
ZC=6.0; ZN=7.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=17; nn=6; pi=0.0;
k1=1.02; k2=1.02; k3=1.02; k4=1.2; k5=1.25; k6=1.3; k7=1.26; k8=1.85; k9=1.1; k10=1.1; k11=0.78; k12=0.99;  
bohr=0.529177; rad=57.29578;

(* 17 Variables: Cr,Nr,Or,NH,CH,CN,CO,OL,NM,OC,OO,OM,pnh,pch,w1,w2,w3 *)

(* 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.25*k1/Cr^2;
vee=3.0*sig1/Cr;
vne=-3.0*ZC/Cr;

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

(* 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*(2*k4/NH^2+k5/CN^2+k6/CO^2+k7/CH^2+k8/NM^2+k9/OL^2+k10/OM^2+k11/OC^2+k12/OO^2);

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

(*  1  2  3  4  5  6  7  8  9  10 11   12   13   14   15 16 17  *)
(*  C  N  O  Hn Hn Hc B  B  L2 L1 L1   Bp   Bp   Bp   Bp L1 L1  *)
as={C,N,O,H,H,H};    (* atom string *)  
rr={Cr,Nr,Oa,NH,NH,CH,CN,CO,OL,NM,NM,  OC,  OC,  OO,  OO,OM,OM};
oc={-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,-1/2,-1/2,-1/2,-1/2,-1,-1};
ch={ZC,ZN,ZO, 1, 1, 1, 0, 0, 0, 0, 0,   0,   0,   0,   0, 0, 0};
csw1=Cos[w1/2];  ssw1=Sin[w1/2];
csw2=Cos[w2];    ssw2=Sin[w2];
csw3=Cos[w3];    ssw3=Sin[w3];

dcnx=Cr+Nr+2*CN;             dcox=(Cr+Oa+2*CO)*csw2;  dnhx=dcnx+(Nr+NH)*csw1; dchx=(Cr+CH)*csw3;           dccnx=Cr+CN; 
pnhx=dcnx+(Nr+NH+pnh)*csw1;  pchx=(Cr+CH+pch)*csw3;
dcoy=(Cr+Oa+2*CO)*ssw2;      dnhy=(Nr+NH)*ssw1;       dchy=(Cr+CH)*ssw3;  
pnhy=(Nr+NH+pnh)*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;
dom1x=dcox-(Oa+OM)*ssw2;     dom2x=dcox+(Oa+OM)*ssw2;
dom1y=dcoy+(Oa+OM)*csw2;     dom2y=dcoy-(Oa+OM)*csw2;

(* 6 nuclei: coord *)
xn={0,-dcnx,dcox,-pnhx,-pnhx,pchx};
yn={0,0,dcoy,pnhy,-pnhy,-pchy};
zn={0,0,0,0,0,0};

(* 17 clouds: coord *)        
xc={0,-dcnx,dcox,-dnhx,-dnhx,dchx,-dccnx,dccox,dolx,-dcnx,-dcnx,0,0,dcox,dcox,dom1x,dom2x};
yc={0,0,dcoy,dnhy,-dnhy,-dchy,0,dccoy,doly,0,0,0,0,dcoy,dcoy,dom1y,dom2y};
zc={0,0,0,0,0,0,0,0,0,NM,-NM,OC,-OC,OO,-OO,0,0};

(* epot terms *)
vne=vne-2*(3-(pnh/NH)^2)/NH;
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 + oc[[i]]*oc[[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 + 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},{Nr,0.22},{Oa,0.186},{NH,1.11},{CH,1.12},{CO,0.92},{CN,1.03},{OL,0.905},{NM,0.765},{OC,1.36},{OO,1.3},{OM,0.86},{w1,2.15},{w2,1.0},{w3,1.12},{pch,0.51},{pnh,0.5},{Method -> "Newton"}, {MaxIterations -> 800}]

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+Nr+2*CN)*bohr /. t[[2]]
(Cr+CH+pch)*bohr /. t[[2]]
(Nr+NH+pnh)*bohr /. t[[2]]

plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},Cr],      Circle[{xc[[2]],yc[[2]]},Nr],Circle[{xc[[3]],yc[[3]]},Oa],Circle[{xc[[4]],yc[[4]]},NH],
Circle[{xc[[5]],yc[[5]]},NH],Circle[{xc[[6]],yc[[6]]},CH],Circle[{xc[[7]],yc[[7]]},CN],Circle[{xc[[8]],yc[[8]]},CO],Circle[{xc[[9]],yc[[9]]},OL],Circle[{xc[[10]],yc[[10]]},NM],Circle[{xc[[11]],yc[[11]]},NM],Circle[{xc[[12]],yc[[12]]},OC],Circle[{xc[[13]],yc[[13]]},OC],Circle[{xc[[14]],yc[[14]]},OO],Circle[{xc[[15]],yc[[15]]},OO],Circle[{xc[[16]],yc[[16]]},OM],Circle[{xc[[17]],yc[[17]]},OM],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]} ] /. t[[2]];

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]]},Nr],Circle[{xc[[3]],zc[[3]]},Oa],Circle[{xc[[4]],zc[[4]]},NH],
Circle[{xc[[5]],zc[[5]]},NH],Circle[{xc[[6]],zc[[6]]},CH],Circle[{xc[[7]],zc[[7]]},CN],Circle[{xc[[8]],zc[[8]]},CO],Circle[{xc[[9]],zc[[9]]},OL],Circle[{xc[[10]],zc[[10]]},NM],Circle[{xc[[11]],zc[[11]]},NM],Circle[{xc[[12]],zc[[12]]},OC],Circle[{xc[[13]],zc[[13]]},OC],Circle[{xc[[14]],zc[[14]]},OO],Circle[{xc[[15]],zc[[15]]},OO],Circle[{xc[[16]],zc[[16]]},OM],Circle[{xc[[17]],zc[[17]]},OM],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]];

Show[plot2,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> Automatic, Frame -> True}]

plot3=Graphics3D[{
Sphere[{xc[[1]],yc[[1]],zc[[1]]},Cr],      
Sphere[{xc[[2]],yc[[2]],zc[[2]]},Nr],
Sphere[{xc[[3]],yc[[3]],zc[[3]]},Oa],
{Opacity[0.6],
Sphere[{xc[[4]],yc[[4]],zc[[4]]},NH],
Sphere[{xc[[5]],yc[[5]],zc[[5]]},NH],
Sphere[{xc[[6]],yc[[6]],zc[[6]]},CH],
Sphere[{xc[[7]],yc[[7]],zc[[7]]},CN],
Sphere[{xc[[8]],yc[[8]],zc[[8]]},CO],
{Lighter[Blue],
Sphere[{xc[[9]],yc[[9]],zc[[9]]},OL],
Sphere[{xc[[16]],yc[[16]],zc[[16]]},OM],
Sphere[{xc[[17]],yc[[17]],zc[[17]]},OM]}},
{Opacity[0.4],{Lighter[Green],
Sphere[{xc[[10]],yc[[10]],zc[[10]]},NM],
Sphere[{xc[[11]],yc[[11]],zc[[11]]},NM],
Sphere[{xc[[12]],yc[[12]],zc[[12]]},OC],
Sphere[{xc[[13]],yc[[13]],zc[[13]]},OC],
Sphere[{xc[[14]],yc[[14]],zc[[14]]},OO],
Sphere[{xc[[15]],yc[[15]],zc[[15]]},OO]}},
{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 -> Automatic,Boxed->False}]

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

s[[3,2]]