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

t=0.01889726*ChemicalData["VitaminB12","AtomPositions"] (* convert to bohr *)

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

n=Length[u]
(*u1=StringJoin[u];
u2=StringReplace[u1,{"C"->"0.2556,","H"->"0.1,","N"->"0.213,","O"->"0.178,","P"->"0.35,","Co"->"0.5,"}]*)

u1={0.5,0.35,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1};

u2={0.5,0.35,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.178,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.213,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.2556,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,0.564,
0.564,0.564,0.564};

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

ChemicalData["VitaminB12","MolecularFormulaString"]

(*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[{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["VitaminB12","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]," "],{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];*)

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}]; 
(*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,2}]

Qij=Lij;
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,101}];     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,101}];      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,101}];
Plot7=Table[Graphics3D[{Darker[Gray],{Opacity[0.3],Sphere[{wx[[i]],wy[[i]],wz[[i]]},a[[i,3]]]}}],{i,1,101}];

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

(*PLot9= Table[Graphics3D[{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.1]}],{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}]