(* Kimball model and energy derived from Facio-PDB source, ES 08 May 2013 *) 
(* Graphite stack of 4 layers 50 atoms each, dangling bonds on the outside *)
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

t1=1.889726*{{4.141,3.305,-3.641},{3.254,3.762,-2.621},{1.962,0.786,-4.295},{2.087,0.286,-2.965},{-0.217,-1.734,-4.949},{-1.104,-1.278,-3.929},{-2.397,-4.254,-5.603},{-2.272,-4.753,-4.273},{2.265,4.743,-2.925},{1.378,5.199,-1.906},{0.085,2.223,-3.579},{0.210,1.724,-2.249},{-2.094,-0.297,-4.233},{-2.981,0.160,-3.214},{-4.273,-2.817,-4.887},{-4.148,-3.316,-3.557},{0.388,6.180,-2.210},{-0.499,6.636,-1.190},{-1.792,3.660,-2.864},{-1.666,3.161,-1.533},{-3.971,1.141,-3.518},{-4.858,1.597,-2.498},{-6.150,-1.379,-4.172},{-6.025,-1.878,-2.841},{-1.489,7.617,-1.494},{-2.376,8.074,-0.474},{-3.668,5.098,-2.148},{-3.543,4.598,-0.818},{-5.848,2.578,-2.802},{-6.734,3.034,-1.782},{-8.027,0.058,-3.456},{-7.902,-0.441,-2.126},{-3.366,9.055,-0.778},{-4.252,9.511,0.242},{-5.545,6.535,-1.432},{-5.420,6.036,-0.102},{-7.724,4.015,-2.086},{-8.611,4.472,-1.066},{-9.904,1.496,-2.740},{-9.779,0.996,-1.410},{5.300,1.876,-1.994},{4.413,2.332,-0.974},{3.121,-0.644,-2.648},{3.246,-1.143,-1.317},{0.941,-3.164,-3.302},{0.055,-2.707,-2.282},{-1.238,-5.684,-3.956},{-1.113,-6.183,-2.625},{3.423,3.313,-1.278},{2.536,3.769,-0.258},{1.244,0.793,-1.932},{1.369,0.294,-0.601},{-0.935,-1.727,-2.586},{-1.822,-1.270,-1.566},{-3.115,-4.246,-3.240},{-2.990,-4.746,-1.910},{1.546,4.750,-0.562},{0.660,5.207,0.458},{-0.633,2.231,-1.216},{-0.508,1.731,0.114},{-2.812,-0.289,-1.870},{-3.699,0.167,-0.850},{-4.991,-2.809,-2.524},{-4.866,-3.308,-1.194},{-0.330,6.188,0.154},{-1.217,6.644,1.173},{-2.510,3.668,-0.500},{-2.385,3.169,0.830},{-4.689,1.148,-1.154},{-5.576,1.605,-0.135},{-6.868,-1.371,-1.808},{-6.743,-1.871,-0.478},{-2.207,7.625,0.869},{-3.094,8.082,1.889},{-4.386,5.105,0.215},{-4.261,4.606,1.546},{-6.566,2.586,-0.439},{-7.452,3.042,0.581},{-8.745,0.066,-1.093},{-8.620,-0.433,0.238},{6.459,0.446,-0.346},{5.572,0.902,0.674},{4.279,-2.074,-1.000},{4.404,-2.573,0.330},{2.100,-4.594,-1.654},{1.213,-4.137,-0.634},{-0.079,-7.113,-2.308},{0.046,-7.613,-0.978},{4.582,1.883,0.370},{3.695,2.340,1.389},{2.403,-0.637,-0.284},{2.528,-1.136,1.046},{0.223,-3.156,-0.938},{-0.664,-2.700,0.081},{-1.956,-5.676,-1.592},{-1.831,-6.175,-0.262},{2.705,3.321,1.085},{1.818,3.777,2.105},{0.526,0.801,0.431},{0.651,0.302,1.762},{-1.653,-1.719,-0.223},{-2.540,-1.262,0.797},{-3.833,-4.239,-0.877},{-3.708,-4.738,0.454},{0.828,4.758,1.801},{-0.058,5.214,2.821},{-1.351,2.238,1.147},{-1.226,1.739,2.477},{-3.530,-0.281,0.493},{-4.417,0.175,1.513},{-5.710,-2.801,-0.161},{-5.584,-3.300,1.169},{-1.048,6.195,2.517},{-1.935,6.652,3.537},{-3.228,3.676,1.863},{-3.103,3.176,3.193},{-5.407,1.156,1.209},{-6.294,1.612,2.229},{-7.586,-1.364,0.555},{-7.461,-1.863,1.885},{7.617,-0.984,1.301},{6.730,-0.528,2.321},{5.438,-3.504,0.647},{5.563,-4.003,1.978},{3.259,-6.023,-0.007},{2.372,-5.567,1.013},{1.079,-8.543,-0.661},{1.204,-9.042,0.670},{5.740,0.453,2.017},{4.854,0.910,3.037},{3.561,-2.066,1.363},{3.686,-2.566,2.694},{1.382,-4.586,0.709},{0.495,-4.130,1.729},{-0.797,-7.106,0.055},{-0.672,-7.605,1.385},{3.864,1.891,2.733},{2.977,2.347,3.753},{1.684,-0.629,2.079},{1.809,-1.128,3.409},{-0.495,-3.149,1.425},{-1.382,-2.692,2.445},{-2.674,-5.668,0.771},{-2.549,-6.168,2.101},{1.987,3.328,3.449},{1.100,3.785,4.468},{-0.192,0.809,2.795},{-0.067,0.309,4.125},{-2.372,-1.711,2.141},{-3.258,-1.255,3.160},{-4.551,-4.231,1.487},{-4.426,-4.730,2.817},{0.110,4.766,4.164},{-0.777,5.222,5.184},{-2.069,2.246,3.510},{-1.944,1.747,4.841},{-4.248,-0.274,2.856},{-5.135,0.183,3.876},{-6.428,-2.794,2.202},{-6.303,-3.293,3.533},{8.776,-2.414,2.949},{7.889,-1.957,3.969},{6.597,-4.933,2.295},{6.722,-5.433,3.625},{4.417,-7.453,1.641},{3.531,-6.997,2.661},{2.238,-9.973,0.987},{2.363,-10.472,2.317},{6.899,-0.976,3.665},{6.012,-0.520,4.684},{4.720,-3.496,3.011},{4.845,-3.995,4.341},{2.541,-6.016,2.357},{1.654,-5.559,3.376},{0.361,-8.535,1.703},{0.486,-9.035,3.033},{5.022,0.461,4.380},{4.136,0.918,5.400},{2.843,-2.059,3.726},{2.968,-2.558,5.057},{0.664,-4.578,3.072},{-0.223,-4.122,4.092},{-1.515,-7.098,2.418},{-1.390,-7.597,3.749},{3.146,1.898,5.096},{2.259,2.355,6.116},{0.966,-0.621,4.442},{1.091,-1.120,5.772},{-1.213,-3.141,3.788},{-2.100,-2.685,4.808},{-3.392,-5.661,3.134},{-3.267,-6.160,4.464},{1.269,3.336,5.812},{0.382,3.792,6.832},{-0.910,0.816,5.158},{-0.785,0.317,6.488},{-3.090,-1.704,4.504},{-3.976,-1.247,5.524},{-5.269,-4.223,3.850},{-5.144,-4.723,5.180}};

n=Length[t1]

u=Table["C",{i,1,n}];

u2=Table[0.25514,{i,1,n}]; u1=u2;

cn=Table[6.0,{i,1,n}]

bc=n;
bh=0;

DeleteFile[ "C:\\Kimb_Mathem\\tmp3"];
For[i=1,i<n+1,i++,
{u[[i]],Extract[t1,{i}]}>>> "C:\\Kimb_Mathem\\tmp3"]

t=t1;

(*x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];
m=Table[If[u[[i]]=="C",12.0107,If[u[[i]]=="H",1.00794]],{i,1,n}];
M=Sum[m[[i]],{i,n}]*)

(*coords=(m.t)/M*)
(*no=Table[1,{i,1,n}];
coords=(no.t)/n
x=Table[t[[i,1]]-coords[[1]],{i,1,n}]; 
y=Table[t[[i,2]]-coords[[2]],{i,1,n}];
z=Table[t[[i,3]]-coords[[3]],{i,1,n}];
t=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];*)

(*w3=Sqrt[3.0]/2; ws=0.5;
x=Table[t[[i,1]]*w3-ws*t[[i,2]],{i,1,n}]; 
y=Table[t[[i,1]]*ws+w3*t[[i,2]]-2.8727498,{i,1,n}];
t=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];
MatrixForm[t]*)

(*coords=(m.t)/M*)

dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];

(* Distance Matrix Bohr *)
Dij=Table[(*If[i≥j,*)If[(dc[[i,j]]<2.71 ||(u[[i]]=="Cl" && j==10 && i== 17)),dc[[i,j]],0](*,0]*),{i,1,n},{j,1,n}];
MatrixForm[Dij]

For[i=1 ,i≤n,i++,
For[j=1, j<i,j++,
If[(Dij[[i,j]]≠0 && u[[i]]== "H"),If[u[[j]]=="C",u2[[i]]=0.5693,If[u[[j]]=="F",u2[[i]]=0.5,If[u[[j]]=="N",u2[[i]]=0.5541,If[u[[j]]=="O",u2[[i]]=0.5241]]]]]]]
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]*)

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: + 8.7616 *)

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," "]];

(*MatrixForm[Pij]*)

(* sigma bond cloud radii *)

con=Table[{i},{i,1,n}]; conj=Table[0,{bc}];
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 *)
con1=con2=con3=con4=Table[0,{n}];

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==2(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con1[[j]]=i]];
con1=Take[con1,{1,j}]
n1=Length[con1]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==3(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con2[[j]]=i]];
con2=Take[con2,{1,j}]
n2=Length[con2]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==4(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con3[[j]]=i]];
con3=Take[con3,{1,j}]
n3=Length[con3]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==5(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con4[[j]]=i]];
con4=Take[con4,{1,j}]
n4=Length[con4]

MatrixForm[conn, TableAlignments->Left]

Clear[c,v,p,Q,P,v1,v2,nk,nl,i]; (* from formamide3d.nb *)
rpi=rpo=1.4;
P=3; c={0,0,0,0,0}; 
(*--------------------------------------------------------------------------*)(*next:π-clouds for C=C and C=O,half occ each above and below*)
connc={{10,50,58},{18,58,66},{26,66,74}};
cj=Length[connc]
tpo=Table[{0,0,0},{cj}];
tpu=tpo;
t1=Table[Join[t[[i]],{1}],{i,1,n}];

<< "C:\\Kimb_Mathem\\PIflattrans1.m"

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}];

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

data1=Table[a[[i,3]],{i,1,nb}];
Histogram[data1,{1.092,1.094,0.00005}]

rcc=1.093;  L=u2[[1]]+rcc;
db1=Table[{0,0,0},{i,1,n1}];
db2=Table[{0,0,0},{i,1,n2}];
For[i=1,i≤n1,i++,db1[[i]]=t[[con1[[i]]]]];
For[i=1,i≤n2,i++,db2[[i]]=t[[con2[[i]]]]];

Plot1=Table[Graphics3D[{If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="O",Lighter[Red],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="B",Lighter[Orange],Lighter[Brown]]]]],Sphere[{t[[i,1]],t[[i,2]],t[[i,3]]},u1[[i]]]}],{i,1,n}];
Plot2=Table[Graphics3D[Text[i,{t[[i,1]]-0.6,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,1,n}];

wxc=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}];
wyc=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}];   
wzc=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[{Darker[Gray,0.8],{Opacity[0.4],If[wxc[[i]]≠0,Sphere[{wxc[[i]],wyc[[i]],wzc[[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}];*)

rpo=1.4;
px=Table[x[[i]],{i,1,bc}];
py=Table[y[[i]],{i,1,bc}];
pz=Table[z[[i]]+rpo,{i,1,bc}];
Plot8a=Table[Graphics3D[{Opacity[0.3],Lighter[Cyan],Sphere[{px[[i]],py[[i]],pz[[i]]},rpo]}],{i,1,bc}];
rpi=1.4;
qx=Table[x[[i]],{i,1,bc}];
qy=Table[y[[i]],{i,1,bc}];
qz=Table[z[[i]]-rpi,{i,1,bc}];
Plot8b=Table[Graphics3D[{Opacity[0.3],Lighter[Yellow],Sphere[{qx[[i]],qy[[i]],qz[[i]]},rpi]}],{i,1,bc}];

Plot8c=Table[Graphics3D[{Darker[Brown,0.5],{Opacity[0.3],Sphere[db1[[i]],rcc]}}],{i,1,n1}];
Plot8d=Table[Graphics3D[{Darker[Orange,0.5],{Opacity[0.3],Sphere[db2[[i]],rcc]}}],{i,1,n2}];

Plot9= Table[Graphics3D[{(*Darker[Red,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.15]}],{i,1,nb}];

(*Plot10=Graphics3D[{Darker[Red,0.5],Sphere[{coords[[1]],coords[[2]],coords[[3]]},0.2]}];*)

Show[Plot1,Plot7,(*Plot8,*)Plot8a,Plot8b,Plot8c,Plot8d,Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]
Show[Plot1,Plot2,Plot8c,Plot8d,(*Plot8e,*)Plot9,(*Plot10,*){AspectRatio->Automatic,Axes->Automatic,Boxed->True, AxesLabel -> {xl, yl, zl}, SphericalRegion->True}]

tc=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];   (* core clouds *)
twc=Table[{wxc[[i]],wyc[[i]],wzc[[i]]},{i,1,cw}];  (* Bond cöouds between cores *)
tpo=Table[{px[[i]],py[[i]],pz[[i]]},{i,1,bc}]; (* pi clouds above *)
tpu=Table[{qx[[i]],qy[[i]],qz[[i]]},{i,1,bc}]; (* pi clouds below *) 
tw=Join[tc,twc,tpo,tpu,db1,db2];
bt=Length[tw]
(*MatrixForm[tw]*)
(* coordinates of clouds *)

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

(* nucleus-nucleus distances *)
Length[dc]

(*MatrixForm[dc] *)

tv1=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3 =Table[{0,0,0},{i,1,bc}]; 
tv4 =Table[{0,0,0},{i,1,bc}];
tv5=Table[{0,0,0},{i,1,n1}]; 
tv6=Table[{0,0,0},{i,1,n2}]; 
tv  = Join[tv1,tv2,tv3,tv4,tv5,tv6];
(*MatrixForm[tv]*)
(*For[i=1,i≤40,i++,cn[[onebond[[i]]]]=4.0];
For[i=1,i≤32,i++,cn[[twobond[[i]]]]=5.0];*)

ow=Join[Table[-2.0,{i,1,bc+nb}],Table[-0.5,{i,1,2*bc}],Table[-2.0,{i,1,n1}],Table[-1.0,{i,1,n2}]]
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*bc}],Table[0.0,{i,1,n1+n2}]]

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

(*cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}]]*)
Elcharge=Sum[ow[[i]],{i,1,Length[ow]}]
Length[ow]
Nuccharge=Sum[cnv[[i]],{i,1,Length[cnv]}]
Length[cnv]

Vee=0.0; Vnn=0.0;

(* all clouds repulsion *)
For[i=1,i<bt,i++,
For[j=i+1,j≤bt,j++,
If[dw[[i,j]]≠ 0,Vee=Vee+ow[[i]]*ow[[j]]/dw[[i,j]],Vee]]]
Vee
(* half cloud with other common halfcloud: subtract WW *)
(*For [i=nb+bc+1,i<nb+2*bc+1,i++,Vee=Vee-ow[[i]]*ow[[i+bc]]/dw[[i,i+bc]]];*)
Vee=Vee-bc*0.5*0.25/rpo;
Vee=Vee+0.5*n1/rcc;
Vee

sc=0.385;   (* core clouds *)
For[i=1 ,i<bc+1,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.385;  (* bonding clouds *)
For[i=1 ,i<nb+1,i++,Vee=Vee+3*sw/a[[i,3]]]
Vee

(* all nuclei with all clouds, except i = j *)
Vne=cnv.(ow.dcw)

Vne=Vne-(2*n1+n2)*6.0/L

(* core nuclei in their cloud centers *)

For[i=1,i<bc+1,i++,Vne=Vne-3*cn[[i]]/u2[[i]]]
Vne

(* protons in their bonding cloud, distance u2[[a[[i,1]]]] to cloud center *)
For[i=cw+1,i<nb+1,i++,Vne=Vne-(3-(u2[[a[[i,1]]]]/a[[i,3]])^2)/a[[i,3]]]
Vne

(* nucleus-nucleus *)
For[i=1,i<n,i++,
For[j=i+1,j≤n,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]
Vnn

Vee
Vnn
Vne
Tpi =2.15*bc*0.5*1.125*(1/rpo^2+1/rpi^2) (* pi kin e added *)
Tdb = (2*n1+n2)*1.125/rcc^2
Vtot=Vne+Vee+Vnn
Ttot = Ttot+Tpi+Tdb
Vir   =-Vtot/Ttot
Etot=Ttot+Vtot