How to build a neighbor table for the Hexagonal lattice

I will take this opportunity to showcase the abilities of IGraph/M for lattice generation and mesh / graph / matrix conversions. IGraph/M thrives on user feedback, so if you find it useful, please take some time to write a few comments about your experience. It will help me improve the package.

Non-periodic case

You can directly generate a (non-periodic) lattice with IGraph/M.

<< IGraphM`

mesh = IGLatticeMesh["Hexagonal", Polygon@CirclePoints[3, 6], 
  MeshCellLabel -> {2 -> "Index"}]

The second argument of IGLatticeMesh may be a region. This region will be filled with cells. In this case, I chose a big hexagon to be filled with small hexagonal cells.

The cell adjacency matrix:

am = IGMeshCellAdjacencyMatrix[mesh, 2]

"2" means 2-dimensional cells, i.e. little hexagons. "1" would mean edges and "0" points.

MatrixPlot[am]

If you need the graph,

graph = IGMeshCellAdjacencyGraph[mesh, 2, 
  VertexCoordinates -> Automatic]

Notice that this is actually a triangular connectivity, which could also be generated directly (in some shapes) with IGTriangularLattice. Demo:

{IGTriangularLattice[4], IGTriangularLattice[{3, 5}]}

We could have use IGLatticeMesh too:

IGLatticeMesh["Triangular", {3, 3}]

Let's get the point-to-point connectivity now (instead of the cell-to-cell one):

IGMeshCellAdjacencyGraph[%, 0]

Periodic case

Now let us do the periodic case.

We start with a hex lattice arranged in an $n\times m$ grid.

{n, m} = {5, 6};
mesh = IGLatticeMesh["Hexagonal", {n, m}, MeshCellLabel -> {2 -> "Index"}]

Convert it to a graph. This time I will not preserve the vertex coordinates so that we can get a clearer layout after we make the lattice periodic.

graph = IGMeshCellAdjacencyGraph[mesh, 2, VertexLabels -> "Name"];
graph = VertexReplace[graph, {2, i_} :> i]

I have also converted the vertex names, which were of the form {2, index} (2 indicating two-dimensional mesh cells) to simply index.

We add the extra edges needed for periodic boundary conditions.

extraEdges = DeleteDuplicates@Flatten@Table[
    {If[Mod[i, m] == 0, {i <-> i - m + 1, i <-> Mod[i - 2 m + 1, m n, 1]}, {}],
     If[i <= m, {i <-> i + m n - m, i <-> Mod[i + m n - m + 1, m n, 1]}, {}]},
    {i, m n}
    ]

pgraph = EdgeAdd[graph, extraEdges]

Then we can get (or plot) the graph's adjacency matrix.

IGAdjacencyMatrixPlot[pgraph]

am = AdjacencyMatrix[pgraph]

Extra visualization: here's the graph in 3D with {m,n} = {10,20}:

(* remember to re-create graph and extraEdges after setting {m,n} *)
pgraph = Graph3D[EdgeAdd[graph, extraEdges], VertexLabels -> None]

An alternative solution for the periodic case

The adjacency relations of hexagonal cells form a triangular lattice. There is a function in IGraph/M for directly generating a triangular lattice graph, and it has an option to make it periodic:

IGTriangularLattice[{5, 10}]

IGTriangularLattice[{5, 10}, "Periodic" -> True]

Then you can just get the adjacency matrix again.

Note that the {m,n} syntax in IGLatticeMesh and IGTriangularLattice do not have the exact same meaning—pay attention to the difference if you mix these approaches! The vertex labelling will also be different. Presumably, at some point you will want to use the visualization of the hex lattice mesh to plot your results. Thus it is useful to be able to map back to mesh cell indices.

---

Update: How to do this for the connectivity of vertices in a hexagonal graph?

OP is asking how to do this if the vertices of the graph are the vertices (not faces) of the hexagonal mesh.

The simplest way is to use the same method as above, but start with the dual lattice of the hexagonal one, i.e. a triangular lattice.

IGLatticeMesh["Triangular", {4, 5}]

IGMeshCellAdjacencyGraph[triMesh, 2, VertexCoordinates -> Automatic]

We can also do it directly with the vertices of a hexagonal lattice, but it is a bit more trouble because of those two hanging out points you can see in the graph above.

Let us start by creating the graph directly from a hexagonal mesh.

{n, m} = {4, 5};
graph = IGMeshGraph[
  IGLatticeMesh["Hexagonal", {n, m}],
  VertexShapeFunction -> "Name",
  PerformanceGoal -> "Quality"
 ]

Now we need to add periodicity. This time, I am not going to add extra edges to connect the left and right, top and bottom of the lattice. If we simply repeat this partial lattice in both directions to see which node would need to be connected to which other one, we will immediately see that it is not enough to add connections. It would also be necessary to add two new vertices (red dots in the illustration below).

We are going to merge corresponding vertices at bottom and top, left and right of the lattice. The formulas for correspondences are easy to figure out by making drawings like the one above. For convenience, we will use VertexReplace instead of VertexContract.

bottom = Range[m + 1, 2 n (m + 1), m + 1];
repl1 = Thread[bottom + m -> bottom]
(* {11 -> 6, 17 -> 12, 23 -> 18, 29 -> 24, 35 -> 30, 41 -> 36, 47 -> 42, 53 -> 48} *)

left = Range[1, 2 m];
repl2 = Thread[left + 2 n (m + 1) -> left]
(* {49 -> 1, 50 -> 2, 51 -> 3, 52 -> 4, 53 -> 5, 54 -> 6, 55 -> 7, 56 -> 8, 57 -> 9, 58 -> 10} *)

If you look carefully at the replacement lists, you will notice that we are not done yet. I kept the output for this specific size of lattice so you can see that vertex 53 is replaced by 48 in the top -> bottom replacement and the same vertex 53 is replaced by 5 in the right -> left replacement. This creates an inconsistency. To get the correct result, we also need to merge 5 and 48 in a third step.

repl3 = {2 n (m + 1) -> m}
(* {48 -> 5} *)

The replacement lists must be applied successively and in the correct order, rather than concurrently, because of repeated treatment of the same vertices. We use Fold for this.

pgraph = SimpleGraph@Fold[VertexReplace, graph, {repl1, repl2, repl3}]

In version 11.3, the vertex coordinates are lost in this process. Let us re-add them so we can see the result better, and we can verify that it is correct.

coord = AssociationThread[VertexList[graph], GraphEmbedding[graph]];

pgraph = Graph[pgraph,
 VertexCoordinates -> Normal@KeyTake[coord, VertexList[pgraph]],
 VertexShapeFunction -> "Name", PerformanceGoal -> "Quality"
 ]

Notice that with this layout, 5 and 46 are the two vertices that would have been missing if we naively repeat the lattice in every direction and try to add edges (instead of contracting vertices).

I was still not completely confident about the result. As you can see from the necessity for repl3, it is easy to make mistakes. Thus, let us make further checks. We expect the result to be vertex transitive. That means that for any two vertices, the graph has a symmetry which transforms them into each other. Loosely speaking, all vertices look the same, they cannot be distinguished based on their position in the graph (at least not without a reference point).

IGraph/M has a function for this.

IGVertexTransitiveQ[pgraph]
(* True *)

Are all edges interchangeable as well? That is not the case. Clearly, we have three categories of edges, running in three different directions in the geometrically laid out lattice.

To show this, let us make a function that categorizes edges based whether they may be transformed into each other by any graph automorphisms.

edgeCategory[graph_] := With[{lg = LineGraph[graph]},
  IGPartitionsToMembership[lg]@
   GroupOrbits@PermutationGroup@IGBlissAutomorphismGroup[lg]
  ]

This function returns a category number for each edge, in the same order as EdgeList.

We can use these numbers for colouring:

Graph[pgraph, EdgeStyle -> Thick] // 
 IGEdgeMap[ColorData[100], EdgeStyle -> edgeCategory]

Again, everything looks good. Every vertex is incident to three edges of distinct categories, and there are precisely three categories.

pgraph has the symmetries we expect for an infinite hexagonal lattice.

Just for fun, here's a force directed layout for a $12\times 16$ periodic lattice.

Edges

This problem can be handled elegantly and efficiently using ListCorrelate.

ntab[r_, x_, pad_: "Cyclic"] := (
   2^Partition[Range[r x - 1, 0, -1], x]
    // ListCorrelate[1 - IdentityMatrix[3], #, 2, pad] &
    // IntegerDigits[Join @@ #, 2, r x] &
 )

With and without wrapping:

ntab[5, 6] // MatrixPlot

ntab[5, 6, 0] // MatrixPlot

Vertices

Regarding your comment to Szabolcs the same methods can be applied to a vertex graph.

Now much faster and using far less memory via SparseArray rather than powers-of-two.

ntabV[n_, m_] :=
  Module[{r = m + 1, x = 2 (n + 1), a, k},
    a = Partition[Hold /@ Range[r x], r];
    k[1] = {{1, 1}, {1, 0}};
    k[2] = {{0, 1}, {1, 1}};
    ListCorrelate[k[#], a[[# ;; ;; 2]], 2 (-1)^#] & /@ {2, 1}
     // Thread[{Level[Riffle @@ #, {-1}]}] &
     // SparseArray[Automatic, {r x, r x}, 0,
         {1, {Range[0, 3 r x, 3], #}, ConstantArray[1, 3 r x]}] &
  ]

ntabV[12, 16] // AdjacencyGraph

Borrowing Szabolcs's illustration we can label vertices like this:

From the adjacency table output of ntabV we can create a matching explicit list:

ntabV[4, 5]["AdjacencyLists"];
Thread[Range[0, 59] -> (% - 1)]

{0 -> {6, 11, 54}, 1 -> {6, 7, 55}, 2 -> {7, 8, 56}, 3 -> {8, 9, 57}, 
 4 -> {9, 10, 58}, 5 -> {10, 11, 59}, 6 -> {0, 1, 12}, 7 -> {1, 2, 13}, 
 8 -> {2, 3, 14}, 9 -> {3, 4, 15}, 10 -> {4, 5, 16}, 11 -> {0, 5, 17}, 
 12 -> {6, 18, 23}, 13 -> {7, 18, 19}, 14 -> {8, 19, 20}, 15 -> {9, 20, 21}, 
 16 -> {10, 21, 22}, 17 -> {11, 22, 23}, 18 -> {12, 13, 24}, 19 -> {13, 14, 25}, 
 20 -> {14, 15, 26}, 21 -> {15, 16, 27}, 22 -> {16, 17, 28}, 23 -> {12, 17, 29}, 
 24 -> {18, 30, 35}, 25 -> {19, 30, 31}, 26 -> {20, 31, 32}, 27 -> {21, 32, 33}, 
 28 -> {22, 33, 34}, 29 -> {23, 34, 35}, 30 -> {24, 25, 36}, 31 -> {25, 26, 37}, 
 32 -> {26, 27, 38}, 33 -> {27, 28, 39}, 34 -> {28, 29, 40}, 35 -> {24, 29, 41}, 
 36 -> {30, 42, 47}, 37 -> {31, 42, 43}, 38 -> {32, 43, 44}, 39 -> {33, 44, 45}, 
 40 -> {34, 45, 46}, 41 -> {35, 46, 47}, 42 -> {36, 37, 48}, 43 -> {37, 38, 49}, 
 44 -> {38, 39, 50}, 45 -> {39, 40, 51}, 46 -> {40, 41, 52}, 47 -> {36, 41, 53}, 
 48 -> {42, 54, 59}, 49 -> {43, 54, 55}, 50 -> {44, 55, 56}, 51 -> {45, 56, 57}, 
 52 -> {46, 57, 58}, 53 -> {47, 58, 59}, 54 -> {0, 48, 49}, 55 -> {1, 49, 50}, 
 56 -> {2, 50, 51}, 57 -> {3, 51, 52}, 58 -> {4, 52, 53}, 59 -> {5, 48, 53}}

Observe that the seams wrap, e.g.

 2  -> { 7,  8, 56}
 24 -> {18, 30, 35}
 35 -> {24, 29, 41}
 59 -> { 5, 48, 53}

---

Explanation

Szabolcs implied that I need to explain this better. Let's start with a drawing of the lattice:

hex = {Polygon@CirclePoints[#, {1.1, 90 °}, 6], Yellow, Text[i++, #]} &;
i = 1;
Array[hex@{2 #2 + #, -Sqrt[3] #} &, {5, 6}] // Graphics

We can see that this is a skewed rectangular matrix, equivalent to:

MatrixForm[m = Partition[HoldForm /@ Range[5*6], 6]]

$\left( \begin{array}{cccccc} 1 & 2 & 3 & 4 & 5 & 6 \\ 7 & 8 & 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 & 17 & 18 \\ 19 & 20 & 21 & 22 & 23 & 24 \\ 25 & 26 & 27 & 28 & 29 & 30 \\ \end{array} \right)$

In that rectangular matrix the neighbors of 8 are {1, 2, 3, 7, 9, 13, 14, 15} but in this skewed version 1 and 15 are too far. We eliminate these, as well as the center, using zeros in the "mask" that is our convolution kernel:

k = {{0,  1,  1},            (* visually skewed to make the hexagon apparent *)
       {1,  0,  1},
         {1,  1,  0}};

We apply this kernel using ListCorrelate:

nt = ListCorrelate[k, m, 2];

The neighbors for 8:

nt[[2, 2]]

2 + 3 + 7 + 9 + 13 + 14

Note also that by default ListCorrelate wraps around:

nt[[1, 1]]  (* neighbors of 1 *)

2 + 6 + 7 + 12 + 25 + 26

(ListCorrelate has a fourth parameter that controls padding; if 0 is specified it effectively does not wrap.)

In this example I used HoldForm so that the integers would not sum. In the complete function ntab I used powers of two to allow them to sum, then get the binary output you wanted using IntegerDigits, e.g.

{2, 6, 7, 12, 25, 26};
2^(% - 1)
Total[%]
IntegerDigits[%, 2, 5*6] // Reverse
Position[%, 1] // Flatten

{2, 32, 64, 2048, 16777216, 33554432}

50333794

{0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, \
0, 0, 0}

{2, 6, 7, 12, 25, 26}

The vertex graph case is similar but I needed two different kernels for the "up" and "down" triads.