(* Kimball model and energy derived from Facio-PDB source, ES February 2013 *) 
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

t1=1.889726*{{0.779000,6.964000,-0.134000},{-0.635000,6.488000,0.167000},{-1.620000,7.528000,-0.364000},{-1.381000,8.740000,0.350000},{-0.104000,9.287000,0.052000},{0.977000,8.324000,0.518000},{0.942000,7.032000,-1.231000},{1.737000,6.076000,0.457000},{-0.775000,6.374000,1.261000},{-0.907000,5.270000,-0.529000},{-3.084000,7.130000,-0.166000},{-1.453000,7.696000,-1.450000},{-0.083000,10.493000,0.819000},{-0.011000,9.521000,-1.026000},{2.253000,8.828000,0.104000},{0.948000,8.223000,1.624000},{2.889000,8.151000,0.304000},{1.596000,5.219000,0.067000},{-0.762000,4.087000,0.249000},{-3.740000,8.018000,-0.271000},{-3.294000,6.541000,1.126000},{-3.397000,6.416000,-0.958000},{-3.253000,5.592000,1.019000},{-0.225000,10.283000,1.731000},{-1.875000,3.130000,-0.139000},{-1.614000,1.754000,0.451000},{-0.232000,1.297000,0.007000},{0.790000,2.330000,0.477000},{2.234000,1.944000,0.150000},{-3.098000,3.633000,0.418000},{-2.624000,0.881000,-0.074000},{0.110000,0.065000,0.644000},{0.485000,3.560000,-0.179000},{2.338000,1.398000,-1.174000},{-0.766000,4.286000,1.338000},{-1.954000,3.066000,-1.244000},{-1.675000,1.789000,1.561000},{-0.192000,1.214000,-1.099000},{0.716000,2.468000,1.578000},{2.895000,2.831000,0.229000},{2.611000,1.206000,0.890000},{-3.758000,2.969000,0.252000},{-2.436000,0.008000,0.259000},{-0.079000,-1.099000,-0.152000},{2.321000,0.447000,-1.094000},{1.064000,-2.056000,0.142000},{0.770000,-3.425000,-0.450000},{-0.574000,-3.895000,0.088000},{-1.633000,-2.866000,-0.302000},{-3.047000,-3.262000,0.123000},{2.242000,-1.534000,-0.492000},{1.815000,-4.304000,-0.013000},{-0.943000,-5.128000,-0.533000},{-1.293000,-1.641000,0.347000},{-3.056000,-3.779000,1.463000},{-0.148000,-0.876000,-1.233000},{1.220000,-2.135000,1.238000},{0.751000,-3.373000,-1.559000},{-0.534000,-3.984000,1.192000},{-1.635000,-2.715000,-1.404000},{-3.725000,-2.385000,0.069000},{-3.461000,-4.020000,-0.576000},{2.919000,-2.189000,-0.375000},{1.581000,-5.178000,-0.311000},{-0.713000,-6.288000,0.260000},{-3.061000,-4.731000,1.401000},{-1.885000,-7.232000,0.065000},{-1.563000,-8.576000,0.689000},{-0.270000,-9.093000,0.081000},{0.840000,-8.072000,0.320000},{2.186000,-8.499000,-0.259000},{-3.016000,-6.681000,0.755000},{-2.634000,-9.474000,0.365000},{0.056000,-10.316000,0.757000},{0.454000,-6.852000,-0.319000},{2.012000,-8.989000,-1.595000},{-0.559000,-6.055000,1.331000},{-2.111000,-7.348000,-1.016000},{-1.463000,-8.478000,1.791000},{-0.413000,-9.272000,-1.006000},{0.970000,-7.895000,1.409000},{2.896000,-7.645000,-0.257000},{2.637000,-9.298000,0.365000},{-3.701000,-7.338000,0.710000},{-2.359000,-10.336000,0.651000},{0.219000,-10.106000,1.667000},{1.857000,-9.923000,-1.555000}};

n=Length[t1]

u= {"C","C","C","O","C","C","H","O","H","O","C","H","O","H","O","H","H","H","C","H","O","H","H","H","C","C","C","C","C","O","O","O","O","O","H","H","H","H","H","H","H","H","H","C","H","C","C","C","C","C","O","O","O","O","O","H","H","H","H","H","H","H","H","H","C","H","C","C","C","C","C","O","O","O","O","O","H","H","H","H","H","H","H","H","H","H","H"};

un=Table[If[u[[i]]=="H","Z",u[[i]]],{i,1,n}]

num=Table[i,{i,1,n}];

um=Table[{un[[i]],num[[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.1,If[us[[i,1]]=="N",0.24,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="Cl",0.5]]]]],{i,1,n}]
u2=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.6,If[us[[i,1]]=="N",0.22522,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="Cl",0.5]]]]],{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]]=="Cl",17.0]]]]],{i,1,n}]

bc=0;
For[i=1, i<n+1,i++,If[u[[i]]≠"H", bc=bc+1]];
bc
bh=n-bc

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

t=Table[t1[[num[[i]]]],{i,1,n}];

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

dc=Table[Sqrt[(t[[i,1]]-t[[j,1]])^2+(t[[i,2]]-t[[j,2]])^2+(t[[i,3]]-t[[j,3]])^2],{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 maTetrx Bohr *)
Dij=Table[If[i≥j,If[(dc[[i,j]]<3.0 ||(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+1,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 *)

For[i=2,i<n+2,i++,
For[j=2,j≤i,j++,
If[(Pij[[i,j]]>0 && ( u[[i-1]]=="O" ||  u[[j-1]]=="O" )) ,Print[u[[i-1]],i-1," ",u[[j-1]],j-1,"  ",Extract[Pij,{i,j}]]]]]

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

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[Orange,0.6],{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}];

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}];
Show[Plot1,Plot7,Plot8,Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]
Show[Plot1,Plot2,Plot9,{AspectRatio->Automatic,Boxed->False, 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 *)
twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}]; (* H-bonding clouds *)
tw=Join[tc,twc,twh];
bt=Length[tw]
MatrixForm[tw]
(* coordinates of clouds *)

(* cloud-cloud distances *)
dw=Table[Sqrt[(tw[[i,1]]-tw[[j,1]])^2+(tw[[i,2]]-tw[[j,2]])^2+(tw[[i,3]]-tw[[j,3]])^2],{i,1,bt},{j,1,bt}];
(*MatrixForm[dw]*)

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

tv1=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3= Table[{x[[i]],y[[i]],z[[i]]},{i,bc+1,n}];
tv  = Join[tv1,tv2,tv3];
MatrixForm[tv]

(* nucleus-cloud distances *)
dcw=Table[sq=Sqrt[(tw[[i,1]]-tv[[j,1]])^2+(tw[[i,2]]-tv[[j,2]])^2+(tw[[i,3]]-tv[[j,3]])^2];If[sq> 0.1,-2/sq,0],{i,1,bt},{j,1,bt}];
dcw=Table[If[i== j,dcw[[i,j]]=0.0,dcw[[i,j]]],{i,1,bt},{j,1,bt}];
(*MatrixForm[dcw]*)

cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}]]

Vee=0.0; Vnn=0.0;

(* all clouds repulsion *)
For[i=1,i<bt,i++,
For[j=i+1,j<bt+1,j++,
Vee=Vee+4/dw[[i,j]]]]
Vee

sc=0.24;   (* core clouds *)
For[i=1 ,i<bc+1,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.24;  (* 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=Apply[Plus,cnv.dcw]

(* 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+1,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]
Vnn

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