(* Kimball model and energy derived from Wolfram and MP2/6-31G(d) data source, ES 20 June 2013 *)
Clear[a,u,u1,u2,t];
Plot0=ChemicalData["Glycine","MoleculePlot"]

t1=0.01889726*ChemicalData["Glycine","AtomPositions"];
MatrixForm[t1]

u=ChemicalData["Glycine","VertexTypes"]

n=Length[u]

num1=Table[i,{i,1,n}];
dc1=Table[Norm[t1[[i]]-t1[[j]]],{i,1,n},{j,1,n}];
MatrixForm[Prepend[dc1,num1]]

cn1=Table[If[u[[i]]=="C",6.0,If[u[[i]]=="H",1.0,If[u[[i]]=="N",7.0,If[u[[i]]=="O",8.0,If[u[[i]]=="S",6.0,If[u[[i]]=="P",5.0]]]]]],{i,1,n}]
Vn=0.0;
For[i=1,i<n,i++,
For[j=i+1,j≤n,j++,Vn=Vn+cn1[[i]]*cn1[[j]]/dc1[[i,j]]]]
Vn

un=Table[If[u[[i]]=="H","Z",u[[i]]],{i,1,n}]
um=Table[{un[[i]],num1[[i]]},{i,1,n}];
us=Sort[um];
num=Table[us[[i,2]],{i,1,n}];

For[i=1,i<n+1,i++,If[us[[i,1]]=="Z",us[[i,1]]="H"]];
u1=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.12,If[us[[i,1]]=="N",0.24,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="S",0.5,If[us[[i,1]]=="P",0.5]]]]]],{i,1,n}]
u2=Table[If[us[[i,1]]=="C",(*0.25618*)0.2674,If[us[[i,1]]=="H",0.6,If[us[[i,1]]=="N",0.2126,If[us[[i,1]]=="O",0.1786,If[us[[i,1]]=="S",0.5244,If[us[[i,1]]=="P",0.5785]]]]]],{i,1,n}]
u=Table[us[[i,1]],{i,1,n}]

cn=Table[If[u[[i]]=="C",6.0,If[u[[i]]=="H",1.0,If[u[[i]]=="N",7.0,If[u[[i]]=="O",8.0,If[u[[i]]=="S",6.0,If[u[[i]]=="P",5.0]]]]]],{i,1,n}]

bc=0;
For[i=1, i<n+1,i++,If[u[[i]]≠"H", bc=bc+1]];
bc
bh=n-bc
t=Table[t1[[num[[i]]]],{i,1,n}];
DeleteFile[ "C:\\Kimb_Mathem\\tmp3"];
For[i=1,i<n+1,i++,
{u[[i]],0.529177*Extract[t,{i}]}>>> "C:\\Kimb_Mathem\\tmp3"]

Plot1=Table[Graphics3D[{Sphere[{t[[i]]},u1[[i]]]}],{i,1,n}];
Plot1o=Table[Graphics3D[{Opacity[0.8],If[u[[i]]=="N",Lighter[Blue,0.6],If[u[[i]]=="H",Lighter[White],If[u[[i]]=="O",Lighter[Red,0.6],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="B",Lighter[Orange],Darker[Brown]]]]]],Sphere[{t[[i]]},u1[[i]]]}],{i,1,n}];

Plot2=Table[Graphics3D[Text[i,{t[[i,1]]-0.5,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,1,n}];
Show[Plot1o,Plot2,{AspectRatio->Automatic,SphericalRegion->True}]

(*Distance Matrix Bohr*)
dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];

x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];

(* Distance matrix Bohr *)
Dij=Table[If[(dc[[i,j]]<3.0||(u[[i]]=="Cl" && dc[[i,j]]<3.24)),dc[[i,j]],0],{i,1,n},{j,1,n}];
MatrixForm[Prepend[Dij,num1]]

For[i=1 ,i<n+1,i++,
For[j=1, j<bc+1,j++,
If[(Dij[[i,j]]≠0  && u[[i]]== "H"),If[u[[j]]=="C",u2[[i]]=0.59,If[u[[j]]=="F",u2[[i]]=0.5,If[u[[j]]=="N",u2[[i]]=0.5541,If[u[[j]]=="O",u2[[i]]=0.5]]]]]]]
u2

Fij=Dij-u2;
(*MatrixForm[Fij];*)
Gij=Table[If[i≥j,If[Fij[[i,j]]>0,Fij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Gij];*)

Hij=Transpose[Gij];
(*MatrixForm[Hij];*)

Kij=Hij-u2;
(*MatrixForm[Kij];*)
Lij=Table[If[i≤j,If[Kij[[i,j]]>0,Kij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Lij];*)

Lij=Table[If[i≤j,If[Lij[[i,j]]>0,If[u[[j]]≠"H",0.5*Lij[[i,j]],Lij[[i,j]]],0]," "],{i,1,n},{j,1,n}];  
Qij=Lij;
MatrixForm[Lij]

Nij=Table[If[i≤j,If[Lij[[i,j]]>0,2.25/Lij[[i,j]]^2,0]," "],{i,1,n},{j,1,n}];
MatrixForm[Nij]

(* Total kinetic energy in Eh, = |Etot| *)
Tval=Sum[If[i≤j,Nij[[i,j]],0],{i,1,n},{j,1,n}]
Tcor=Sum[If[u[[i]]≠"H",2.25/u2[[i]]^2,0],{i,1,n}]
Ttot=Tval+Tcor  (* plus π and lone pairs *)
TvalH=Sum[If[(i≤j && u[[j]]== "H"),Nij[[i,j]],0],{i,1,n},{j,1,n}]
TvalC=Tval-TvalH

Mij=Prepend[Qij,u];
Oij=Transpose[Mij];
Pij=Prepend[Oij,Prepend[u," "]];

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,1,bc-1},{j,1,bc-1}];

nc=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="C",nc=nc+1,nc]]
cnt=Table[0,{nc}];
nc
no=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="O",no=no+1,no]]
no

For[j=1,j≤nc,j++,
For[i=bc+1,i≤n,i++,If[ Dij[[i,j]]≠0 ,cnt[[j]]=cnt[[j]]+1,0]]];
cnt

con=Table[{i},{i,1,n}]; conj=Table[0,{nc}];
For[i=1,i≤n,i++,
For[j=1,j≤n,j++,If[ Dij[[i,j]]≠0 ,con[[i]]={con[[i]],j}]]];
conn=Table[Flatten[con[[i]]],{i,1,n}];

(* identify C atoms at double bonds *)
j=0;
For[i=1,i≤nc,i++,If[(Length[conn[[i]]]==4 && u[[conn[[i,4]]]]≠"O"),j=j+1;conj[[j]]=i]];
conj=Take[conj,{1,j}]
nj=Length[conj]
MatrixForm[conn, TableAlignments->Left]

j=0; h=0;
coco=Table[{0,0,0,0},{no}];
coch=coco;
For[i=nc+1, i≤nc+no, i++,If[Length[conn[[i]]]==2,j=j+1;coco[[j]]=Join[{conn[[i,1]]},Take[conn[[conn[[i,2]]]],3]],h=h+1;coch[[h]]=conn[[i]]]];
coco=Take[coco,{1,j}];
coch=Take[coch,{1,h}];
oj=Length[coco]
hj=Length[coch]

Clear[c,v,p,Q,P,v1,v2,nk,nl,i]; 
rpo={1.2961,1.0875}; (* from formamide3d.nb *)
rpi=rpo;

(*--------------------------------------------------------------------------*)
(*next: π-clouds for C=C and C=O,half occ each above and below*)

k=1;
connc={{2,1,4},{5,2,4}};
cj=Length[connc]
tpo=Table[{0,0,0},{cj}];
tpu=tpo;

P=3; c={0,0,0,0,0};

Do[(c = connc[[k]]; (*c = Take[c, P];*) 
v = Table[Join[t[[c[[i]]]], {1}], {i, 1, P}]; 
x2 = t[[c[[2]],1]]; 
y2 = t[[c[[2]],2]]; 
z2 = t[[c[[2]],3]];
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {-x2, -y2, -z2, 1}}; 
p = v.Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}];
T1 = ArcTan[p[[1,3]]/p[[1,2]]]; 
U = p; 
T = RotationMatrix[T1, {1, 0, 0}]; 
p = U.T; 
T2 = ArcTan[p[[1,2]]/p[[1,1]]]; 
U = p; 
T = RotationMatrix[T2, {0, 0, 1}]; 
p = U.T; 
T3 = ArcTan[p[[3,3]]/p[[3,2]]]; 
U = p; 
T = RotationMatrix[T3, {1, 0, 0}]; 
p = U.T; 
mk = p; 
ml = p; 
mk[[1,3]] = rpo[[k]]; 
ml[[1,3]] = -rpi[[k]]; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T3, {1, 0, 0}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T2, {0, 0, 1}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T1, {1, 0, 0}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
Clear[v, v1, v2]; 
v = Table[Join[p[[i]], {1}], {i, 1, P}];
v1 = Table[Join[mk[[i]], {1}], {i, 1, P}]; 
v2 = Table[Join[ml[[i]], {1}], {i, 1, P}]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {x2, y2, z2, 1}}; 
p = v.Q; 
mk = v1.Q; 
ml = v2.Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
mk = Table[Take[mk[[i]], 3], {i, 1, P}]; 
ml = Table[Take[ml[[i]], 3], {i, 1, P}];
tpo[[k]] = mk[[1]]; 
tpu[[k]] = ml[[1]]), {k, cj}];

MatrixForm[Table[Join[{i},{u[[i]]},t[[i]]],{i,1,n}]]
x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];

(*-------------------------------------------------------------------------*)
(* next O-atoms: 2 doubly occupied LP's C=O *)

k=1;
coco={{1,5,4}};
oj=Length[coco];
rlo = 1.125;
tplo=Table[{0,0,0},{oj}];
tplu=Table[{0,0,0},{oj}];
Okk={{-1,1,1}};
Okl={{1,1,1}};
<< "C:\\Kimb_mathem\\OD1trans.m"

(*-------------------------------------------------------------------------*)
(*next:O atoms in C-O-H,C-O-C,O-O-H:two doubly occ LP*)
(*tpso=Table[{0,0,0},{hj}];
tpsu=tpso;
P=3.0;
rlos={1.119};
coch
Ok={{-1,-1,1}};*)
Clear[coch];
coch = {{4,2,8}};
hj=Length[coch]
tpso=Table[{0,0,0},{hj}];
tpsu=tpso;
rlos={1.125,1.125};
Ok={{1,-1,1}};

<< "C:\\Kimb_mathem\\OBtrans.m"

(*-------------------------------------------------------------------------*)
(* next: N3 in C-NH2; one doubly occ LP as pyramidal ligand: Ok *)

k = 1;
Clear[P,tpln];
P = 4;
tpln = {0.0, 0.0, 0.0};
cnch=Table[{0,0,0,0},{1}];
cnch[[1]]= conn[[3]];
cnch[[1]]=RotateLeft[cnch[[1]],1]
rln = 1.246;

Do[(c = cnch[[k]];
   a1 = dc[[c[[2]], c[[1]]]]; a2 = dc[[c[[2]], c[[3]]]]; d1 = dc[[c[[1]], c[[3]]]];
   α = ArcCos[(a1^2 + a2^2 - d1^2)/(2*a1*a2)];
   L = (u2[[c[[2]]]] + rln);
   β = ArcSin[rln/L];
   γ = π/2 - α/2; (*γ=3*π/2+α/2;*)
   v = Table[Join[t[[c[[i]]]], {1}], {i, 1, P}];
   x2 = t[[c[[2]], 1]];
   y2 = t[[c[[2]], 2]]; 
   z2 = t[[c[[2]], 3]];
   Q = {{1, 0, 0, 0},
     {0, 1, 0, 0},
     {0, 0, 1, 0},
     {-x2, -y2, -z2, 1}};
   p = v.Q;
   p = Table[Take[p[[i]], 3], {i, 1, P}];
   T1 = ArcTan[p[[1, 3]]/p[[1, 2]]];
   U = p;
   T = RotationMatrix[T1, {1, 0, 0}];
   p = U.T;
   T2 = ArcTan[p[[1, 2]]/p[[1, 1]]];
   U = p;
   T = RotationMatrix[T2, {0, 0, 1}];
   p = U.T;
   T3 = ArcTan[p[[3, 3]]/p[[3, 2]]];
   U = p;
   T = RotationMatrix[T3, {1, 0, 0}];
   p = Chop[U.T];
   nk = p;
   (*nk[[4,2]]=L*Cos[γ]; nk[[4,1]]=-L*Sin[γ];*)
   nk[[4,3]] = nk[[4,3]] - u2[[c[[4]]]] - rln; 
   U = p;
   U1 = nk;
   T = RotationMatrix[-T3, {1, 0, 0}];
   p = U.T;
   nk = U1.T;
   U = p;
   U1 = nk;
   T = RotationMatrix[-T2, {0, 0, 1}];
   p = U.T;
   nk = U1.T;
   U = p;
   U1 = nk;
   T = RotationMatrix[-T1, {1, 0, 0}];
   p = U.T;
   nk = U1.T;
   Clear[v, v1, v2];
   v  = Table[Join[p[[i]], {1}], {i, 1, P}];
   v1 = Table[Join[nk[[i]], {1}], {i, 1, P}];
   Q = {{1, 0, 0, 0},
     {0, 1, 0, 0},
     {0, 0, 1, 0},
     {x2, y2, z2, 1}};
   p = v.Q;
   nk = v1.Q;
   p = Table[Take[p[[i]], 3], {i, 1, P}];
    nk = Table[Take[nk[[i]], 3], {i, 1, P}];
   tpln[[k]] = nk[[4]]), {k, 1}];

(* sigma bond cloud radii *)

DeleteFile[ "C:\\Kimb_Mathem\\tmp4"];
For[i=2,i<n+2,i++,
For[j=2,j≤i,j++,
If[Pij[[i,j]]>0,{i-1,j-1,Extract[Pij,{i,j}]}>>> "C:\\Kimb_Mathem\\tmp4"]]]

a=ReadList["C:\\Kimb_Mathem\\tmp4"]
nb=Length[a]
cw=nb-bh

wx=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,2]]]]-x[[a[[i,1]]]])+x[[a[[i,1]]]],{i,1,cw}];      wy=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,2]]]]-y[[a[[i,1]]]])+y[[a[[i,1]]]],{i,1,cw}];  wz=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,2]]]]-z[[a[[i,1]]]])+z[[a[[i,1]]]],{i,1,cw}];
Plot7=Table[Graphics3D[{Opacity[0.5],Sphere[{wx[[i]],wy[[i]],wz[[i]]},a[[i,3]]]}],{i,1,cw}];

wxh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,1]]]]-x[[a[[i,2]]]])+x[[a[[i,2]]]],{i,cw+1,nb}];    wyh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,1]]]]-y[[a[[i,2]]]])+y[[a[[i,2]]]],{i,cw+1,nb}];    wzh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,1]]]]-z[[a[[i,2]]]])+z[[a[[i,2]]]],{i,cw+1,nb}];

Plot8=Table[Graphics3D[{Opacity[0.5],Sphere[{wxh[[i]],wyh[[i]],wzh[[i]]},a[[i+cw,3]]]}],{i,1,bh}];
Plot8a=Table[Graphics3D[{Opacity[0.3],Lighter[Cyan],Sphere[tpo[[i]],rpo[[i]]]}],{i,1,cj}];
Plot8b=Table[Graphics3D[{Opacity[0.3],Lighter[Yellow],Sphere[tpu[[i]],rpi[[i]]]}],{i,1,cj}];
Plot8d=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tplo[[i]],rlo]}],{i,1,oj}];
Plot8e=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tplu[[i]],rlo]}],{i,1,oj}];
Plot8f=Graphics3D[{Opacity[0.4],Lighter[Blue],Sphere[tpln[[1]],rln]}];
Plot8ib=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tpso[[i]],rlos[[i]]]}],{i,1,hj}];
Plot8jb=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tpsu[[i]],rlos[[i]]]}],{i,1,hj}];
Plot9= Table[Graphics3D[{Lighter[Gray,0.8],Tube[{{x[[a[[i,1]]]],y[[a[[i,1]]]],z[[a[[i,1]]]]},{x[[a[[i,2]]]],y[[a[[i,2]]]],z[[a[[i,2]]]]}},0.07]}],{i,1,nb}];

Show[Plot1o,Plot2,Plot7,Plot8,Plot8a,Plot8b,Plot8d,Plot8e,Plot8f,Plot8ib,Plot8jb,Plot9,{AspectRatio->Automatic,Boxed->False,SphericalRegion->True}]

tc=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];   (* core clouds *)
twc=Table[{wx[[i]],wy[[i]],wz[[i]]},{i,1,cw}];  (* Bond cöouds between cores *)
twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}]; (* H-bonding clouds *) 
tw=Join[tc,twc,twh,tpo,tpu,tplo,tplu,tpso,tpsu,{tpln[[1]]}];
bt=Length[tw]

(* cloud-cloud distances *)
dw=Table[Norm[tw[[i]]-tw[[j]]],{i,1,bt},{j,1,bt}];
MatrixForm[tw]

(* nucleus-nucleus distances *)
tv1=Table[t[[i]],{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3= Table[t[[i]],{i,bc+1,n}];
tv4 =Table[{0,0,0},{i,1,2*cj}]; 
tv5=Table[{0,0,0},{i,1,2*oj}];
tv6 =Table[{0,0,0},{i,1,2*hj}];
tv  = Join[tv1,tv2,tv3,tv4,tv5,tv6,{{0,0,0}}];
MatrixForm[tv]
Length[tv]

ow=Join[Table[-2.0,{i,1,bc+nb}],Table[-0.5,{i,1,2*cj}],Table[-2.0,{i,1,2*oj}],Table[-2.0,{i,1,2*hj}],{-2.0}];
 ow
cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}],Table[0.0,{i,1,2*cj}],Table[0.0,{i,1,2*oj}],Table[0.0,{i,1,2*hj}],{0.0}]

(* nucleus-cloud distances *)
dcw=Table[sq=Norm[tw[[i]]-tv[[j]]];If[(sq> 0.1 && i≠j),1/sq,0],{i,1,bt},{j,1,bt}];
Length[dcw]

Elcharge=Sum[ow[[i]],{i,1,Length[ow]}]
Length[ow]
Nuccharge=Sum[cnv[[i]],{i,1,Length[cnv]}]
Length[cnv]

Vee=0.0;
For[i=1,i<bt,i++,
For[j=i+1,j≤bt,j++,
Vee=Vee+ow[[i]]*ow[[j]]/dw[[i,j]]]]
Vee
Vee=Vee-0.5*0.25*Sum[1/rpo[[i]],{i,1,cj}];

sc=0.418;   (* core clouds *)
For[i=1 ,i≤bc,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.418;  (* bonding clouds *)
For[i=1 ,i≤nb,i++,Vee=Vee+3*sw/a[[i,3]]]
For[i=1 ,i≤hj,i++,Vee=Vee+2*3*sw/rlos[[i]]]
For[i=1 ,i≤oj,i++,Vee=Vee+2*3*sw/rlo]
Vee=Vee+3*sw/rln;
Vee

(* all nuclei with all clouds, except i = j *)
Vne=cnv.(ow.dcw)
For[i=1,i≤bc,i++,Vne=Vne-3.06122*cn[[i]]/u2[[i]]]
Vne

For[i=cw+1,i≤nb,i++,Vne=Vne-(3-(u2[[a[[i,1]]]]/a[[i,3]])^2)/a[[i,3]]]
Vne

(* nucleus-nucleus *)

(*Vnn=0.0;
For[i=1,i<n,i++,
For[j=i+1,j<n+1,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]*)
Vnn=Vn;

ratio=-(Vne+2*Vnn)/(Vee-Vnn)

ratiomp2=-(-1024.6230504966+2*179.3207697874)/(278.9647594823-179.3207697874)

Vee
Vnn
Vne
Tpi = 2.2*1.125*(0.5*Sum[(1/rpo[[i]]^2+1/rpi[[i]]^2),{i,1,cj}])
Tlp =1.533*2.25*(2*oj/rlo^2+Sum[2/rlos[[i]]^2,{i,1,hj}]+1/rln)

Vtot=Vne+Vee+Vnn
Ttot =Ttot+Tlp
Vir   =-Vtot/Ttot
Etot=Ttot+Vtot