Constructing Icosahedral Fullerene Graphs

Question

I am interested in constructing graphs of icosahedral fullerenes, following the construction explained in the first few pages of this article (alternate link). Mathematica has built-in information about the 60-site graph, $C_{60}$, in GraphData["TruncatedIcosahedralGraph"] and ChemicalData["Fullerene60"] and information about the 20-site dodecahedron, but there's a whole zoo of icosahedral fullerenes.

First, I construct an underlying hexagonal grid of ions:

Li = 10;
Lj = 10;
offset = {-3 - 1/2, Sqrt[3]/2};
a1 = {3/2, +Sqrt[3]/2};
a2 = {3/2, -Sqrt[3]/2};
ions = Flatten[Outer[offset + #1 a1 + #2 a2 + #3 (a1 + a2)/3 &,
    Range[-Li, Li], Range[-Lj, Lj], {0, 1}
    (* A/B shift comes last to ensure even/odd sites alternate *)], 2];

with a face centered on the origin. Then I compute the locations of the pentagonal defects:

Gi = 1;
Gj = 2;
b1 = {3/2, Sqrt[3]/2};
b2 = {0, Sqrt[3]};
upRight = {+1, +1} Gj b2 + {+1, -1}*Gi b1;
downRight = {+1, -1} Gi b2 + {+1, +1}*Gj b1 - upRight;
horizontal = upRight + downRight;
pentagonCenters = Flatten[{# horizontal, # horizontal + upRight, # horizontal + 
    2 upRight, # horizontal + downRight} & /@ Range[0, 4], 1] ~Join~
   {5 horizontal, 5 horizontal + upRight};

So that you can visually see what I've got,

Show[ListPlot[{ions, pentagonCenters}, AspectRatio -> 1, 
              PlotRange -> {{0, 20}, {-10, 10}}, 
              PlotMarkers -> ({#, Large} & /@ {\[FilledSmallCircle], ★})],
     Graphics[Line /@ (pentagonCenters[[#]] & /@ {
     {1, 3}, {4, 7}, {8, 11}, {12, 15}, {16, 19}, {20, 22}, {1, 4}, {2, 8}, {3, 12},
     {7, 16}, {11, 20}, {15, 21}, {19, 22}, {1, 21}, {2, 22}})]
 ]

which, for Gi = 1; Gj = 2 yields the $C_{60}$ graph where the blue circles are locations of the ions and the stars will be the centers of the pentagonal defects. The triangles with stars as their vertices are the faces of the icosahedron.

Running again with Gi=1;Gj=3 yields a 170-ion graph

My questions are all aimed at constructing the adjacency matrix for these icosahedral fullerenes.

1. How can I easily "cut out" the ions that aren't inside the black triangles. Also: the sites I generate are definitely too many...

   [EDIT]: This can be accomplished with by construction a region:

    region = Region[Polygon[pentagonCenters[[
        {1, 2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 22, 
         21, 20, 17, 16, 13, 12, 9, 8, 5, 4}]]]];
   

   and then using RegionMember. I still have to pick Li and Lj large enough to cover the area of interest. [Updated question]: how can I just generate the "needed" hexagonal grid to begin with, or just the needed ones plus a little extra, that can then be discarded via RegionMember?

2. How can I "stitch up" the sides of the black triangles? On a hexagonal grid I can detect if two ions are nearest neighbors by seeing if the Norm of the difference of their locations differ by 1. How should I make these identifications in an easy, programmatic way that Mathematica will understand?

3. Bonus: In an ideal world, I'd know the geometrical locations of all of these ions when the graph is thought of as a polyhedron. Is there something relatively automatic I can do to get at this information? [EDIT: once I have a Graph object I can cast to a Graph3D and use AbsoluteOptions[ graph3D, VertexCoordinates] to get the locations.]

Note that sometimes ions land exactly on the sides of the black triangles (as in the Gi=1;Gj=2 case).

Answer 1

Jason B gave a hint to the tool buckygen. This tool is easily downloaded, unpacked and compiled (by running make in the unpacked folder).

Assuming that the compiled executable of buckygen is contained in the subfolder buckygen-1.0 of folder Downloads of the current user (and maybe also assuming that you run this on a unixoid OS), you can run buckygen and import the results with the following tiny program. It generates all fullerenes with vertex count between m and n. Also some but not all options of buckygen got translated into Mathematica options.

ClearAll[getBucky];
getBucky[m_Integer, n_Integer, OptionsPattern[{
    "Dual" -> False,
    "IPROnly" -> False,
    "Statistics" -> False,
    "Verbose" -> True
    }]
  ] := Module[{file, runstring, path, data, time},
  path = FileNameJoin[{$HomeDirectory, "Downloads", "buckygen-1.0"}];
  file = FileNameJoin[{path, "FullerDump.g6"}];
  runstring = StringJoin[{"! ", FileNameJoin[{path, "buckygen"}],
     " -s",
     If[OptionValue["Dual"], " -d", ""],
     If[OptionValue["IPROnly"], " -I", ""],
     If[OptionValue["Statistics"], " -v", ""],
     If[! OptionValue["Verbose"], " -q", ""],
     " -S", ToString[m], " ",
     ToString[n], " ", file, " 2>&1 &"}];
  Print[Import[runstring, "Text"]];
  Print["Importing ", file, "..."];
  time = AbsoluteTiming[data = Import[file, "Sparse6"];][[1]];
  Print["Time elapsed for import: ", time, "."];
  Flatten[{data}]
  ]

A usage example:

data = getBucky[20,20];
GraphicsGrid[Partition[data, 6, 6, {1, 1}, Graphics[]], ImageSize -> Full]

The dual graphs can be obtained with

data = getBucky[20, 20, "Dual" -> True];
GraphicsGrid[Partition[data, 6, 6, {1, 1}, Graphics[]], ImageSize -> Full]

The graphs are imported as a list of Mathematica Graph objects. (I am using version 11.0.1 and do not know what happens for legacy versions of Mathematica.) For example, you can obtain the AdjacencyMatrix of the first fullerene in the list with AdjacencyMatrix[data[[1]]]. Or you can obtain planar embeddings with PlanarGraph:

GraphicsGrid[Partition[PlanarGraph /@ data, 6, 6, {1, 1}, Graphics[]], ImageSize -> Full]

If the executable is located somewhere else, adjust the variable path in the program to your needs. If you would like to use other options of buckygen, you merely have to adjust the variable runstring to your needs.

What's really bugging me is that Import for Sparse6 is very, very slow. Maybe, someone else knows an improvement.