ChemicalData["Aspirin","MoleculePlot"]

t=0.01889726*ChemicalData["Aspirin","AtomPositions"];

ChemicalData["Aspirin","BondTally"]

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

u1={0.1887,0.1887,0.1887,0.1887,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04};

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

ChemicalData["Aspirin","MolecularFormulaString"]

v=0.01889726*ChemicalData["Aspirin","VertexCoordinates"];

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[Graphics[{Circle[{v[[i,1]],v[[i,2]]},u1[[i]]]}],{i,1,21}];

Plot2=Table[Graphics[Text[i,{v[[i,1]]-0.5,v[[i,2]]} ,{1,0}]],{i,1,21}];
(*Show[Plot2,{AspectRatio->Automatic,Axes -> False,Frame -> Automatic,PlotRange -> {{0,14},{-6,6}}}];*)

Show[Plot1,lp,Plot2,{AspectRatio->Automatic,Axes -> False,Frame -> Automatic,PlotRange -> {{0,14},{-6,6}}}]

ChemicalData["Aspirin","CHStructureDiagram"]

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

(* 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,21},{j,1,21}];
(*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,21},{j,1,21}];*)
(*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,21},{j,1,21}];
MatrixForm[Eij];*)

u2={0.1887,0.1887,0.1887,0.1887,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.2555,0.563,0.563,0.563,0.563,0.563,0.563,0.563,0.50};

Fij=Dij-u2;
(*MatrixForm[Fij];*)
Gij=Table[If[i≥j,If[Fij[[i,j]]>0,Fij[[i,j]],0]," "],{i,1,21},{j,1,21}];
(*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,21},{j,1,21}];
(*MatrixForm[Lij];*)

Lij=Table[If[i≤j,If[Lij[[i,j]]>0,If[j<14,0.5*Lij[[i,j]],Lij[[i,j]]],0]," "],{i,1,21},{j,1,21}]; Qij=Lij;
(*MatrixForm[Lij];*)

Nij=Table[If[i≤j,If[Lij[[i,j]]>0,2.25/Lij[[i,j]]^2,0]," "],{i,1,21},{j,1,21}];
(*MatrixForm[Nij];*)

Tsigma=Sum[If[i≤j,Nij[[i,j]],0],{i,1,21},{j,1,21}]+Sum[2.25/u2[[i]]^2,{i,1,13}]

Mij=Prepend[Qij,u];
Oij=Transpose[Mij];
Pij=Prepend[Oij,Prepend[u," "]];

(*MatrixForm[Pij];*)

(* sigma bond cloud radii *)

For[i=2,i<23,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<23,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"]

a1=a[[11,1]];
a2=a[[11,2]];
a3=a[[11,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot4=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[14,1]];
a2=a[[14,2]];
a3=a[[14,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot5=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[3,1]];
a2=a[[3,2]];
a3=a[[3,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot6=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[5,1]];
a2=a[[5,2]];
a3=a[[5,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot7=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[1,1]];
a2=a[[1,2]];
a3=a[[1,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot8=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[2,1]];
a2=a[[2,2]];
a3=a[[2,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot9=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[4,1]];
a2=a[[4,2]];
a3=a[[4,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot10=Graphics[{Circle[{wx,wy},0.7*a3]}];

a1=a[[6,1]];
a2=a[[6,2]];
a3=a[[6,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot11=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[7,1]];
a2=a[[7,2]];
a3=a[[7,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot12=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[8,1]];
a2=a[[8,2]];
a3=a[[8,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot13=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[9,1]];
a2=a[[9,2]];
a3=a[[9,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot14=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[10,1]];
a2=a[[10,2]];
a3=a[[10,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot15=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[12,1]];
a2=a[[12,2]];
a3=a[[12,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot16=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[13,1]];
a2=a[[13,2]];
a3=a[[13,3]];
q=(a3+u2[[a1]])/(u2[[a1]]+2*a3+u2[[a2]]);
wx=q*(v[[a2,1]]-v[[a1,1]])+v[[a1,1]];
wy=q*(v[[a2,2]]-v[[a1,2]])+v[[a1,2]];
Plot17=Graphics[{Circle[{wx,wy},0.7*a3]}];a1=a[[4,1]];

a1=a[[15,1]];
a2=a[[15,2]];
a3=a[[15,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot18=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[16,1]];
a2=a[[16,2]];
a3=a[[16,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot19=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[17,1]];
a2=a[[17,2]];
a3=a[[17,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot20=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[18,1]];
a2=a[[18,2]];
a3=a[[18,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot21=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[19,1]];
a2=a[[19,2]];
a3=a[[19,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot22=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[20,1]];
a2=a[[20,2]];
a3=a[[20,3]];
a3=a3-0.2555;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot23=Graphics[{Circle[{wx,wy},0.56*a3]}];

a1=a[[21,1]];
a2=a[[21,2]];
a3=a[[21,3]];
a3=a3-0.1887;
q=(a3+u2[[a2]])/(u2[[a1]]+a3+u2[[a2]]);
wx=q*(v[[a1,1]]-v[[a2,1]])+v[[a2,1]];
wy=q*(v[[a1,2]]-v[[a2,2]])+v[[a2,2]];
Plot24=Graphics[{Circle[{wx,wy},0.56*a3]}];

Show[Plot1,lp,Plot2,Plot4,Plot5,Plot6,Plot7,Plot8,Plot9,Plot10,Plot11,Plot12,Plot13,Plot14,Plot15,Plot16,Plot17,Plot18,Plot19,Plot20,Plot21,Plot22,Plot23,Plot24,{AspectRatio->Automatic,Axes -> False,Frame -> Automatic,PlotRange -> {{0,14},{-6,6}}}]