(* 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*{{5.814000,7.901000,-30.790000},{3.149000,9.127000,-32.553000},{7.561000,7.561000,-33.767000},{8.734000,7.364000,-29.079000},{4.284000,8.853000,-27.804000},{5.414000,8.329000,-32.810000},{4.176000,8.724000,-33.336000},{4.234000,8.885000,-34.812000},{5.508000,8.110000,-35.193000},{6.231000,7.974000,-33.857000},{5.069000,6.778000,-35.737000},{4.414000,10.364000,-34.998000},{5.765000,10.849000,-34.445000},{5.916000,12.319000,-34.970000},{5.142000,12.710000,-35.755000},{6.849000,13.100000,-34.365000},{7.823000,7.602000,-31.299000},{8.325000,7.402000,-32.599000},{9.709000,6.944000,-32.504000},{10.063000,6.898000,-31.198000},{8.840000,7.323000,-30.421000},{10.619000,6.520000,-33.613000},{11.404000,6.476000,-30.694000},{12.230000,5.860000,-31.424000},{11.980000,6.812000,-29.317000},{6.414000,8.128000,-28.791000},{7.622000,7.715000,-28.293000},{7.607000,7.550000,-26.801000},{6.354000,8.377000,-26.431000},{5.580000,8.411000,-27.700000},{7.295000,6.054000,-26.600000},{6.117000,9.115000,-25.348000},{6.780000,8.846000,-23.951000},{4.151000,8.893000,-30.252000},{3.559000,9.101000,-29.003000},{2.191000,9.503000,-29.146000},{1.987000,9.548000,-30.540000},{3.201000,9.181000,-31.159000},{1.206000,9.693000,-28.061000},{1.006000,9.719000,-31.567000},{-0.232000,9.845000,-31.472000},{1.731000,9.571000,-32.936000},{0.954000,8.467000,-33.637000},{0.490000,8.533000,-34.647000},{0.870000,7.311000,-32.986000},{0.000000,6.330000,-33.596000},{7.224000,14.416000,-34.971000},{8.030000,14.064000,-36.222000},{9.229000,14.456000,-36.548000},{10.000000,15.424000,-35.643000},{9.845000,13.847000,-37.852000},{10.210000,14.811000,-39.006000},{10.445000,13.908000,-40.267000},{11.728000,14.401000,-40.937000},{11.656000,15.933000,-41.054000},{11.893000,13.769000,-42.326000},{13.295000,14.105000,-42.858000},{13.478000,13.468000,-44.244000},{14.922000,13.678000,-44.734000},{15.274000,15.178000,-44.701000},{15.069000,13.129000,-46.168000},{16.549000,13.150000,-46.600000},{16.688000,12.603000,-48.036000},{18.176000,12.549000,-48.442000},{18.306000,12.130000,-49.919000},{18.932000,11.550000,-47.546000},{8.063000,7.337000,-34.708000},{9.632000,7.088000,-28.527000},{3.754000,9.033000,-26.869000},{3.346000,8.469000,-35.288000},{6.100000,8.660000,-35.924000},{5.945000,6.194000,-36.020000},{4.506000,6.241000,-34.974000},{4.437000,6.933000,-36.612000},{3.637000,10.934000,-34.488000},{4.379000,10.647000,-36.050000},{6.503000,10.162000,-34.858000},{5.668000,10.782000,-33.361000},{10.549000,7.234000,-34.434000},{11.646000,6.486000,-33.248000},{10.326000,5.531000,-33.966000},{12.343000,7.840000,-29.315000},{11.203000,6.700000,-28.561000},{12.805000,6.136000,-29.093000},{8.518000,7.936000,-26.344000},{7.262000,5.829000,-25.534000},{6.330000,5.820000,-27.050000},{8.072000,5.454000,-27.073000},{5.046000,9.035000,-25.159000},{6.425000,10.129000,-25.603000},{6.426000,9.584000,-23.232000},{6.512000,7.847000,-23.607000},{7.864000,8.920000,-24.042000},{1.096000,10.757000,-27.852000},{0.244000,9.282000,-28.367000},{1.551000,9.181000,-27.163000},{1.753000,10.514000,-33.483000},{-0.003000,5.422000,-32.993000},{-1.013000,6.729000,-33.656000},{0.359000,6.099000,-34.599000},{7.811000,14.925000,-34.206000},{6.277000,14.909000,-35.190000},{7.425000,14.425000,-37.054000},{8.131000,12.979000,-36.191000},{10.832000,15.857000,-36.199000},{9.333000,16.219000,-35.309000},{10.384000,14.885000,-34.777000},{9.116000,13.146000,-38.259000},{10.770000,13.344000,-37.570000},{11.113000,15.333000,-38.688000},{9.357000,15.480000,-39.121000},{9.577000,14.032000,-40.914000},{10.542000,12.882000,-39.911000},{12.575000,14.104000,-40.318000},{12.563000,16.305000,-41.530000},{10.791000,16.212000,-41.655000},{11.563000,16.370000,-40.060000},{11.131000,14.178000,-42.990000},{11.774000,12.689000,-42.236000},{14.045000,13.711000,-42.172000},{13.401000,15.187000,-42.935000},{12.793000,13.925000,-44.958000},{13.274000,12.398000,-44.193000},{15.615000,13.136000,-44.090000},{16.297000,15.319000,-45.049000},{14.591000,15.726000,-45.350000},{15.183000,15.551000,-43.681000},{14.494000,13.737000,-46.867000},{14.706000,12.102000,-46.219000},{17.143000,12.531000,-45.927000},{16.932000,14.170000,-46.572000},{16.156000,13.247000,-48.736000},{16.271000,11.598000,-48.097000},{18.635000,13.529000,-48.312000},{19.360000,12.095000,-50.196000},{17.861000,11.145000,-50.058000},{17.789000,12.854000,-50.549000},{19.981000,11.519000,-47.840000},{18.856000,11.865000,-46.505000},{18.494000,10.558000,-47.657000}};

n=Length[t1]

u={"Mg","C","C","C","C","N","C","C","C","C","C","C","C","C","O","O","N","C","C","C","C","C","C","O","C","N","C","C","C","C","C","C","C","N","C","C","C","C","C","C","O","C","C","O","O","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","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]]=="Mg",0.35]]]]],{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]]=="Mg",0.35]]]]],{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]]=="Mg",12.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]]=="Mg",Lighter[Orange],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