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

t1=1.889726*{{-1.920,2.730,2.009},{-3.015,2.354,1.285},{-3.258,3.144,0.079},{-1.959,4.110,0.067},{-1.271,3.687,1.284},{-3.068,2.386,-1.157},{-1.969,2.739,-1.885},{-1.286,3.672,-1.156},{-1.602,2.313,3.096},{-1.673,2.330,-2.979},{-3.972,1.358,1.747},{-0.049,4.340,1.736},{-4.068,1.441,-1.633},{-0.061,4.308,-1.622},{1.150,3.644,-1.159},{1.801,2.685,-1.880},{2.907,2.322,-1.165},{3.131,3.098,0.056},{-3.726,0.022,1.225},{-4.398,-0.465,0.021},{-3.012,-0.924,1.905},{-3.074,-2.078,1.178},{-3.894,-2.002,-0.030},{-3.821,0.070,-1.210},{-3.137,-0.837,-1.968},{-3.102,-2.002,-1.260},{-2.476,-0.782,2.975},{-2.674,-0.650,-3.066},{1.844,4.080,0.051},{1.162,3.672,1.278},{1.812,2.718,2.007},{2.899,2.329,1.278},{1.482,2.265,-2.965},{1.502,2.317,3.101},{3.898,1.374,-1.660},{3.690,0.007,-1.203},{3.067,-0.945,-1.958},{-2.523,-3.342,1.647},{3.070,-2.098,-1.227},{-2.465,-3.217,-1.751},{3.826,-2.037,0.021},{4.292,-0.488,0.035},{3.620,0.003,1.238},{2.901,-0.937,1.921},{2.973,-2.098,1.207},{3.862,1.344,1.751},{2.617,-0.799,-3.067},{2.353,-0.784,2.984},{2.453,-3.336,-1.677},{2.403,-3.357,1.675},{1.080,-3.490,-1.215},{0.001,-3.034,-1.917},{-1.117,-3.420,-1.235},{-0.863,-4.249,-0.056},{0.755,-4.281,-0.031},{1.038,-3.583,1.221},{-0.062,-3.276,1.970},{-1.154,-3.572,1.206},{0.028,-2.448,-2.970},{-0.068,-2.835,3.091},{-4.591,3.911,0.128},{-2.309,5.609,0.054},{4.474,3.849,0.064},{-5.930,-0.328,0.085},{-5.043,-3.024,-0.065},{2.213,5.574,0.035},{4.999,-3.032,0.064},{5.820,-0.310,0.068},{-1.504,-5.646,-0.150},{1.340,-5.705,-0.060},{-5.006,1.683,1.521},{-3.974,1.328,2.862},{-4.124,1.485,-2.745},{-5.084,1.760,-1.330},{-2.439,-3.204,-2.864},{-3.101,-4.098,-1.541},{-3.168,-4.190,1.349},{-2.571,-3.379,2.759},{2.441,-3.390,2.788},{3.046,-4.211,1.386},{3.861,1.322,2.864},{4.894,1.672,1.524},{-0.047,4.401,2.848},{-0.044,5.404,1.431},{-0.051,5.378,-1.339},{-0.062,4.351,-2.735},{4.922,1.712,-1.412},{3.906,1.398,-2.774},{3.069,-4.204,-1.373},{2.481,-3.391,-2.791},{-4.723,4.558,-0.766},{-4.668,4.535,1.046},{-5.485,3.255,0.149},{-1.424,6.277,0.041},{-2.881,5.902,0.962},{-2.897,5.881,-0.850},{4.577,4.485,0.970},{4.594,4.482,-0.843},{5.360,3.181,0.077},{-6.418,-0.773,-0.811},{-6.287,0.721,0.117},{-6.342,-0.809,1.000},{-4.706,-4.080,-0.120},{-5.686,-2.880,-0.961},{-5.672,-2.962,0.851},{2.800,5.857,0.937},{1.336,6.254,0.034},{2.793,5.839,-0.875},{5.594,-2.916,0.997},{5.669,-2.911,-0.815},{4.688,-4.096,0.043},{6.263,-0.757,0.986},{6.149,0.750,0.068},{6.302,-0.761,-0.826},{-1.228,-6.280,0.722},{-1.208,-6.166,-1.087},{-2.613,-5.631,-0.152},{1.002,-6.299,0.817},{2.448,-5.735,-0.017},{1.060,-6.242,-0.993}};

(*MatrixForm[t1]*)

u={"C","H","C","C","H","H","C","H","O","O","C","C","C","C","H","C","H","C","H","C","C","H","C","H","C","H","O","O","C","H","C","H","O","O","C","H","C","C","H","C","C","C","H","C","H","C","O","O","C","C","H","C","H","C","C","H","C","H","O","O","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"};

n=Length[u]

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]]=="P",0.5,If[us[[i,1]]=="S",0.5,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]]=="P",0.5,If[us[[i,1]]=="S",0.5,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]]=="P",15.0,If[u[[i]]=="S",16.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]]=="P",Lighter[Magenta],If[u[[i]]=="S",Lighter[Yellow],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}];

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[sq=Sqrt[(x[[i]]-x[[j]])^2+(y[[i]]-y[[j]])^2+(z[[i]]-z[[j]])^2]; sq<3.3(*||(u[[i]]=="S" && j==21 && i== 48))*),sq,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,If[u[[j]]=="P",u2[[i]]=0.5241,If[u[[j]]=="S",u2[[i]]=0.4]]]]]]]]]
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,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[Gray],{Opacity[0.2],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 *)

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]*)
(* cloud-cloud distances *)

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}];
(*MatrixForm[dc]*)  (* nucleus-nucleus distances *)

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

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]*)  
(* nucleus-cloud distances *)

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