(* Kimball model and energy derived from experimental molecule data source, ES end january 2013 *)
ChemicalData["Fullerene","MoleculePlot"]

(* t=0.01889726*ChemicalData["Fullerene","AtomPositions"];*) (* convert to bohr *)
t=1.8897*ReadList["C:\\Kimb_Mathem\\C60.cc2"]

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

n=Length[u]

60

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,If[us[[i,1]]=="B",5.0]]]]]],{i,1,n}]

u1={0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555};

(*For[i=1,i<15,i++,If[u[[i]]="C",u1[[i]]=0.26,u1[[i]]=0.18]]
u1*)

ChemicalData["Fullerene","MolecularFormulaString"]

M=ChemicalData["Fullerene","MolecularWeight"]

m=Table[If[u[[i]]=="C",12.0107,If[u[[i]]=="H",1.00794]],{i,1,n}];
coords=(m.t)/M[[1]]

(*lp=ListPlot[v, Axes->False];*)

(*Plot1=Table[Graphics[{Circle[{v[[i,1]]/0.5292,v[[i,2]]/0.5292},If[i<5,0.18,0.26]]}],{i,1,13}];*)
Plot1=Table[Graphics3D[{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.5,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,1,n}];
(*Show[Plot2,{AspectRatio->Automatic,Axes -> False,Frame -> Automatic,PlotRange -> {{0,14},{-6,6}}}];*)

Show[Plot1,Plot2,{AspectRatio->Automatic}]

ChemicalData["Fullerene","CHStructureDiagram"]

(* Distance matrix Bohr *)

x=Table[t[[i,1]]-coords[[1]],{i,1,n}]; 
y=Table[t[[i,2]]-coords[[2]],{i,1,n}];
z=Table[t[[i,3]]-coords[[3]],{i,1,n}];
t=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];
coords=(m.t)/M[[1]]

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]," "],{i,1,n},{j,1,n}];
(*Dij=Table[If[sq=Sqrt[(x[[i]]-x[[j]])^2+(y[[i]]-y[[j]])^2+(z[[i]]-z[[j]])^2]; sq<3.0,sq," "],{i,1,n},{j,1,n}];*)
(*MatrixForm[Dij]*)

(*Eij=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,0.529177*sq,0]," "],{i,1,n},{j,1,n}];
MatrixForm[Eij];*)

u2=u1;

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[j<61,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];*)

(* Total kinetic energy in Eh, = |Etot| *)
Tsigma=Sum[If[i≤j,Nij[[i,j]],0],{i,1,n},{j,1,n}]+Sum[2.25/u2[[i]]^2,{i,1,n}]

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]

wx=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,nb}];
wy=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,nb}];   wz=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,nb}];
Plot7=Table[Graphics3D[{Opacity[0.5],Sphere[{wx[[i]],wy[[i]],wz[[i]]},a[[i,3]]]}],{i,1,nb}];

rpo=1.7;
px=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];x[[i]]*(d+rpo)/d,{i,1,n}];
py=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];y[[i]]*(d+rpo)/d,{i,1,n}];
pz=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];z[[i]]*(d+rpo)/d,{i,1,n}];
Plot8=Table[Graphics3D[{Opacity[0.4],Lighter[Cyan],Sphere[{px[[i]],py[[i]],pz[[i]]},rpo]}],{i,1,n}];
rpi=1.2;
qx=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];x[[i]]*(d-rpi)/d,{i,1,n}];
qy=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];y[[i]]*(d-rpi)/d,{i,1,n}];
qz=Table[d=Sqrt[x[[i]]^2+y[[i]]^2+z[[i]]^2];z[[i]]*(d-rpi)/d,{i,1,n}];
Plot8a=Table[Graphics3D[{Opacity[0.6],Lighter[Yellow],Sphere[{qx[[i]],qy[[i]],qz[[i]]},rpi]}],{i,1,n}];

Plot9= Table[Graphics3D[{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.1]}],{i,1,nb}];

Show[Plot1,(*Plot7,*)Plot8,Plot8a,Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]

Show[Plot1,Plot7,(*Plot8,Plot8a,*)Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]

Plot8=Table[Graphics3D[{Opacity[0.2],Lighter[Cyan],Sphere[{px[[i]],py[[i]],pz[[i]]},rpo]}],{i,1,n}];
Plot8a=Table[Graphics3D[{Opacity[0.2],Lighter[Yellow],Sphere[{qx[[i]],qy[[i]],qz[[i]]},rpi]}],{i,1,n}];
Show[Plot1,Plot7,Plot8,Plot8a,Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]

Show[Plot1,Plot2,Plot9,{AspectRatio->Automatic,SphericalRegion->True}]

tc=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];   (* core clouds *)
tpo=Table[{px[[i]],py[[i]],pz[[i]]},{i,1,n}]; (* pi clouds above *)
tpu=Table[{qx[[i]],qy[[i]],qz[[i]]},{i,1,n}]; (* pi clouds below *) 
twc=Table[{wx[[i]],wy[[i]],wz[[i]]},{i,1,nb}];  (* Bond clouds between cores *)
(*twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}];*) (* H-bonding clouds *)
tw=Join[tc,twc,tpo,tpu];
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 *)
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]*)

tv1=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];
tv2=Table[{0,0,0},{i,1,nb}];
(*tv3= Table[{x[[i]],y[[i]],z[[i]]},{i,bc+1,n}];*)
tv4 =Table[{0,0,0},{i,1,n}]; 
tv5 =Table[{0,0,0},{i,1,n}]; 
tv  = Join[tv1,tv2,(*tv3,*)tv4,tv5];
(*MatrixForm[tv]*)

ow=Join[Table[-2.0,{i,1,n+nb}],Table[-0.5,{i,1,2*n}]]

cnv=Join[Table[cn[[i]],{i,1,n}],(*Table[0.0,{i,1,cw}],*)Table[0.0,{i,1,nb}],Table[0.0,{i,1,2*n}]]

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,1/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 *)

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+ow[[i]]*ow[[j]]/dw[[i,j]]]];
(* half cloud with other common halfcloud: subtract WW *)
For [i=nb+1,i<nb+n+1,i++,Vee=Vee-ow[[i]]*ow[[i+n]]/dw[[i,i+n]]];
Vee

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

(* core nuclei in their cloud centers *)

For[i=1,i<n+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
Tpi = 1.64*n*0.5*1.125*(1/rpo+1/rpi) (* pi kin e added *)
Vtot=Vne+Vee+Vnn
Ttot=Tsigma+Tpi
-Vtot/Ttot
Etot=Ttot+Vtot