(* Kimball model and energy derived from experimental molecule data source, ES FePenary 2013 *) 
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

ChemicalData["Eicosane","MoleculePlot"]

(*t1=0.01889726*ChemicalData["Eicosane","AtomPositions"]*)
(*MatrixForm[t1]*)
t1={{-0.062566938134`,2.235545858`,-4.047793091999999`},{1.51223433424`,3.7346654937999997`,-2.1839563382`},{-2.8767298897999996`,2.7520079737999996`,-3.8034515202`},{4.2224037744`,2.7565433162`,-1.9666378481999998`},{-4.5083193182`,0.8268496113`,-5.233596156999999`},{5.6493358769999995`,3.8932135052`,0.30891350922`},{-7.2967990037999995`,0.9503621026599999`,-4.4760050036`},{8.1150503618`,2.4921706487999997`,0.9020607061`},{-8.813682064`,-1.24980808462`,-5.5518260154000005`},{9.2987747282`,3.3514290609999997`,3.3839323481999997`},{-11.464400724199999`,-1.4563551364199998`,-4.3692354846`},{11.474038326799999`,1.63733419544`,4.211254391`},{-12.8754591284`,-3.8491828893999998`,-5.205628212200001`},{12.458396600199999`,2.2909148298`,6.8600833252`},{-15.373109982599999`,-4.2127661718`,-3.8113883694`},{14.4656635574`,0.38983157654`,7.778112216`},{-16.691004895`,-6.6799924374`,-4.5340195918`},{15.3222763532`,0.9638169517799999`,10.467192313999998`},{-19.112688763999998`,-7.109338184599999`,-2.9935149566`},{17.184223381`,-1.00822551278`,11.4577866832`},{0.2310190035`,0.27484174944`,-3.6505726868`},{0.52179114312`,2.7196936592`,-5.9220233388`},{1.60445296304`,5.6387534114`,-2.8580216024`},{0.6379148058199999`,3.6246834405999997`,-0.3641502002`},{-3.2319983777999997`,4.5691684954`,-4.6158447276`},{-3.375995499`,2.6915367418`,-1.8448828020199999`},{4.129807200399999`,0.7553990712399999`,-1.69226853026`},{5.2199901298`,3.2333211859999995`,-3.6596433715999996`},{-3.8282069308`,-1.01901584824`,-4.766266917199999`},{-4.358841991599999`,1.19710362648`,-7.215729758399999`},{6.1308380618`,5.7990021762`,-0.16485969624`},{4.4489819218`,3.8047743283999997`,1.9335676431999997`},{-8.06913002`,2.662623934`,-5.2247144448`},{-7.4262452348`,0.9062181033`,-2.4585335259999996`},{7.6790905736`,0.5252493417`,1.076198957`},{9.4304886304`,2.8304316028`,-0.5959061968399999`},{-7.81401701`,-2.9543976283999998`,-5.123425131199999`},{-9.011914321399999`,-0.9873062459599999`,-7.5468097536`},{10.0463503338`,5.209029718999999`,3.102930092`},{7.8769448857999995`,3.338200979`,4.8218248616`},{-12.539465845599999`,0.138439437034`,-4.9932230098`},{-11.275050178999999`,-1.4864017798199998`,-2.3562993494`},{10.774650734199998`,-0.25894915378`,4.2666233628`},{12.99186625`,1.8320893569999999`,2.889391054`},{-11.6979708578`,-5.4352299211999995`,-4.7730699308`},{-13.2541602188`,-3.7174689871999997`,-7.1873838683999995`},{13.3218124096`,4.1171460362`,6.7724000388`},{10.905608746`,2.2773088026`,8.155301525599999`},{-16.595195786799998`,-2.6873793446`,-4.3291732934`},{-15.009715672799999`,-4.2293957606`,-1.82230057632`},{13.641931993999998`,-1.45671418436`,7.7586480382`},{16.0668283972`,0.48080298618`,6.5465777818`},{-15.422053886`,-8.2044344016`,-4.141334529`},{-17.169105573`,-6.6142299726`,-6.497633878399999`},{16.2607142848`,2.7548425628`,10.453019369`},{13.701080417799998`,0.99486514996`,11.6752941458`},{-20.39959217`,-5.585652110799999`,-3.323083171`},{-18.654808154199998`,-7.190974347799999`,-1.0257432728`},{-19.970624368`,-8.8509096662`,-3.5585430305999997`},{18.820348151799998`,-1.0468137177`,10.2702828648`},{16.304555928`,-2.8287308494`,11.452117505199999`},{17.731865975799998`,-0.5300681429999999`,13.344678094199999`}};

MatrixForm[t1]

(*ChemicalData["Eicosane","BondTally"]*)

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

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

ChemicalData["Eicosane","MolecularFormulaString"]

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

Plot1=Table[Graphics3D[{If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="O",Lighter[Red],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.3,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,1,n}];

ChemicalData["Eicosane","CHStructureDiagram"]

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[i≥j,If[sq=Sqrt[(x[[i]]-x[[j]])^2+(y[[i]]-y[[j]])^2+(z[[i]]-z[[j]])^2]; sq<3.0,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]]]]]]]
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[If[u[[a[[i,1]]]]≠"H",(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]]]],0],{i,1,nb}];

wyc=Table[If[u[[a[[i,1]]]]≠"H",(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]]]],0],{i,1,nb}];

wzc=Table[If[u[[a[[i,1]]]]≠"H",(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]]]],0],{i,1,nb}];

Plot7=Table[Graphics3D[{Darker[Gray],{Opacity[0.5],If[wxc[[i]]≠0,Sphere[{wxc[[i]],wyc[[i]],wzc[[i]]},a[[i,3]]]]}}],{i,1,nb}];

wxh=Table[If[u[[a[[i,1]]]]=="H",(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]]]],0],{i,1,nb}];

wyh=Table[If[u[[a[[i,1]]]]=="H",(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]]]],0],{i,1,nb}];

wzh=Table[If[u[[a[[i,1]]]]=="H",(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]]]],0],{i,1,nb}];

Plot8=Table[Graphics3D[{Opacity[0.5],If[wxh[[i]]≠0,Sphere[{wxh[[i]],wyh[[i]],wzh[[i]]},a[[i,3]]]]}],{i,1,nb}];

Plot9= Table[Graphics3D[{Darker[Red,0.8],Cylinder[{{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.05]}],{i,1,nb}];
Show[Plot1,Plot7,Plot8,Plot9,{AspectRatio->Automatic,Boxed->False}]
Show[Plot1,Plot2,Plot9,{AspectRatio->Automatic}]

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+cw]],wyh[[i+cw]],wzh[[i+cw]]},{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

(*For[i=1,i<bt+1,i++,
For[j=1,j<bt+1,j++,(*sq=Sqrt[(tw[[i,1]]-t[[j,1]])^2+(tw[[i,2]]-t[[j,2]])^2+(tw[[i,3]]-t[[j,3]])^2];*)If[ i≠j,Vne=Vne+cnv[[j]]*dcw[[i,j]],0]]]
Vne*)
(*For[i=1,i<bt+1,i++,
For[j=1,j<bt+1,j++,Vne=Vne-If[(dcw[[i,j]]≠0 && i≠j),
2*cnv[[j]]/dcw[[i,j]],0]]]
Vne*)

(* protons with bonding clouds*)
(*For[i=bc+1,i<bt+1,i++,
For[j=bc+1,j<n+1,j++,If[i≠j,
Vne=Vne-2/Sqrt[(tw[[i,1]]-t[[j,1]])^2+(tw[[i,2]]-t[[j,2]])^2+(tw[[i,3]]-t[[j,3]])^2],0]]]
Vne*)

(* protons with core clouds*)
(*For[i=1,i<bc+1,i++,
For[j=bc+1,j<n+1,j++,
Vne=Vne-2/Sqrt[(t[[i,1]]-t[[j,1]])^2+(t[[i,2]]-t[[j,2]])^2+(t[[i,3]]-t[[j,3]])^2]]]
Vne*)

(* core nuclei with core clouds*)
(*For[i=1,i<bc+1,i++,
For[j=i+1,j<bc+1,j++,
Vne=Vne-2*cnv[[j]]/Sqrt[(t[[i,1]]-t[[j,1]])^2+(t[[i,2]]-t[[j,2]])^2+(t[[i,3]]-t[[j,3]])^2]]]
Vne*)

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