(* Kimball model and energy derived from MP2/6-31g(d) opt, ES 19 April 2013 *) 
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

Plot0=ChemicalData["Ethane","MoleculePlot"]
t1=1.889727*{{0.6922245314,0.1305130785,-0.2933526108},{-0.6922148293,-0.1305340513,0.2933592995},{1.4667473928,-0.3924513130,0.2739295845},{0.7524315730,-0.2101397044,-1.3304280679},{0.9326338768,1.1969101538,-0.2792275020},{-1.4667423894,0.3924572298,-0.2739009182},{-0.7524467881,0.2101714838,1.3304148740},{-0.9326333672,-1.1969238771,0.2792053408}};

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

n=Length[t1]

MatrixForm[Table[Join[{u[[i]]},t1[[i]]],{i,1,n}]]

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

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

um=Table[{un[[i]],num1[[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.12,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.25618*)0.2562,If[us[[i,1]]=="H",0.6,If[us[[i,1]]=="N",0.22522,If[us[[i,1]]=="O",0.1786,If[us[[i,1]]=="Cl",0.48698]]]]],{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",7.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["Ethane","MolecularFormulaString"]

t=Table[t1[[num[[i]]]],{i,1,n}];
Plot1o=Table[Graphics3D[{Opacity[0.8],If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="H",Lighter[White],If[u[[i]]=="O",Lighter[Red],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="B",Lighter[Orange],Darker[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}];

ChemicalData["Ethane","CHStructureDiagram"]

(* Distance Matrix Bohr *)

dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];
(*MatrixForm[Prepend[dc,num1]]*)

x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];

(* connected pairs: distance criteria < 3.0 bohr *)
Dij=Table[If[(dc[[i,j]]<3.0||(u[[i]]=="Cl" && dc[[i,j]]<3.24)),dc[[i,j]],0],{i,1,n},{j,1,n}];
MatrixForm[Prepend[Dij,num1]]

For[i=1 ,i<n+1,i++,
For[j=1, j<bc+1,j++,
If[(Dij[[i,j]]≠0  && u[[i]]== "H"),If[u[[j]]=="C",u2[[i]]=0.59,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

(* subtract core radii from pair distances *)
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];*)

(* determine bond cloud radii *)
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];*)

(* determine kinetic energy *)
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," "]];

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,1,bc+1},{j,1,bc+1}];
(*CD=Select[Flatten[Eij],(#≥2.8 &&#≤2.9) &]
Mean[CD] ± StandardDeviation[CD]*)
CC=Select[Flatten[Dij],(#>2.8 &&#≤3.0) &]
Mean[CC] ± StandardDeviation[CC]

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,bc+1,n},{j,1,bc}];
CH=Select[Flatten[Eij],(#>2.0 &&#≤2.11) &]
Mean[CH] ± StandardDeviation[CH]

(* count C atoms *)
nc=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="C",nc=nc+1,nc]]
cnt=Table[0,{nc}];
pri=cnt; sec=cnt;ter=cnt; qua=cnt;
nc

(* count H atoms bonded to each C atom *)
For[j=1,j≤nc,j++,
For[i=bc+1,i≤n,i++,If[ Dij[[i,j]]≠0 ,cnt[[j]]=cnt[[j]]+1,0]]];
cnt

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,2,nb}];
Histogram[data1,{1.219735,1.219750,0.0000002}]

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[Brown,0.4],Opacity[0.5],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[{Lighter[Orange,0.5],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[Plot1o,Plot2,Plot7,Plot8,(*Plot8a,Plot8b,*)(*Plot8d,Plot8e,*)(*Plot8f,Plot8g,*)Plot9,{AspectRatio->Automatic,Boxed->False,(* Axes->Automatic,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 *)
twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}]; (* H-bonding clouds *) 
tw=Join[tc,twc,twh(*tpo,tpu,tplo,tplu,tpso,tpsu*)];
bt=Length[tw]
(*MatrixForm[tw]*)
(* coordinates of clouds *)

(* cloud-cloud reciprocal distances *)
dwr=Table[sq=Norm[tw[[i]]-tw[[j]]];If[ i≠j,1/sq,0],{i,1,bt},{j,1,bt}];
(*MatrixForm[dw]*)
Length[dwr]

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

tv1=Table[t[[i]],{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3= Table[t[[i]],{i,bc+1,n}];
(*tv4 =Table[{0,0,0},{i,1,cj}]; 
tv5 =Table[{0,0,0},{i,1,cj}]; 
tv6 =Table[{0,0,0},{i,1,8}]; *)
tv  = Join[tv1,tv2,tv3(*,tv4,tv5*)(*,tv6*)];
(*MatrixForm[tv]*)
Length[tv]

ow=Join[Table[-2.0,{i,1,bc+nb}](*,Table[-0.5,{i,1,2*cj}]*)];
ow
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*cj}]*)]

(* nucleus-cloud reciprocal distances *)
dcwr=Table[sq=Norm[tw[[i]]-tv[[j]]];If[ i≠j,1/sq,0],{i,1,bt},{j,1,bt}];
Length[dcwr]
(*MatrixForm[dcw] *)

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 *)
Vee=0.5*ow.(ow.dwr)

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

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

(* 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 *)
dcr=Table[If[i≠j,1/dc[[i,j]],0],{i,1,n},{j,1,n}];
Vnn=0.5*cn.(cn.dcr)

(* now add the energy terms of Cl[Ne]7+ from Ne_s24b2_atoms.nb *)
Vee
Vnn
Vne
Vtot=Vne+Vee+Vnn
Ttot
Vir   =-Vtot/Ttot
Etot=Ttot+Vtot
Mean[CC]*0.529177
Mean[CH]*0.529177