A robust algorithm for k-point grid generation and symmetry reduction

Codes that solve the many-body problem using density functional theory (DFT) use uniform grids over the Brillouin zone in order to calculate the total electronic energy, among other material properties. The total electronic energy is calculated by numerically integrating the occupied electronic bands. For metallic systems, there exist surfaces of discontinuities at the boundary between occupied and unoccupied states, collectively known as the Fermi surface. These discontinuities cause the accuracy in the calculation of the total electronic energy to converge extremely slowly and erratically with respect to grid density. This is demonstrated in figure 1 where we compare the convergence of an insulator (silicon) with a metal (aluminum).

Zoom In Zoom Out Reset image size

The poor convergence of the electronic energy means that DFT codes must use extremely dense grids [2, 3] to achieve an accuracy of several meV/atom. To reduce computation time, it is common practice to evaluate eigenvalues at symmetrically equivalent ${\bf{k}}$-points only once. This is the essence of 'symmetry reducing' a ${\bf{k}}$-point grid.

In most DFT codes, even for very dense grids, the setup and symmetry reduction of the grid takes a few seconds at most. Our motivation for an improved algorithm (despite the speed of current routines) is two-fold: 1) enable an automatic grid-generation technique that allows us to scan over thousands of candidate grids, in a few seconds, to find one with the best possible symmetry reduction [4, 5] (in other words, enable a ${\bf{k}}$-point generation method in the same spirit as that of [2] but have the grid generation done on-the-fly [6]), and 2) eliminate (or at least greatly reduce) the probability of incorrect symmetry reduction³ as the result of finite precision errors (the danger of these increases as the density of the integration grid increases).

In this brief report, an algorithm for generating, and subsequently symmetry-reducing, ${\bf{k}}$-point grids is explained. This algorithm builds on concepts such as Hermite Normal Form, Smith Normal Form, and the connection between finite groups and integer matrices. These concepts are briefly explained in the main text; for more details, see the appendix and [8]. The algorithm has been implemented in an open-source code available at https://github.com/msg-byu/kgridGen and incorporated in version 6 of the VASP code [9]. The algorithm has been incorporated into a code for generating generalized regular grids https://github.com/msg-byu/GRkgridgen [6].

As demonstrated in figure 3, every uniform sampling of a reciprocal unit cell can be expressed through the simple integer relationship

$$\begin{eqnarray}&&{\mathbb{R}}={\mathbb{KN}}\end{eqnarray} \tag{ 1 }$$

where ${\mathbb{R}},{\mathbb{K}}$, and ${\mathbb{N}}$ are 3 × 3 matrices; the columns of ${\mathbb{R}}$ are the reciprocal lattice vectors, and the columns of ${\mathbb{K}}$ are the ${\bf{k}}$-point grid generating vectors. Put simply, ${\mathbb{N}}$ describes the integer linear combination of vectors of ${\mathbb{K}}$ that are equivalent to ${\mathbb{R}}$. One obtains Monkhorst–Packgrids (regular grids) when ${\mathbb{N}}$ is an integer, diagonal matrix. More generally, when ${\mathbb{N}}$ is an invertible, integer matrix, one obtains generalized regular grids. Examples of Monkhorst–Packand generalized regular grids are given in figure 2. We use R to refer to the infinite lattice of points defined by integer linear combinations of ${\mathbb{R}}$, and K to refer to the lattice of points defined by ${\mathbb{K}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With no loss of generality, a new basis for the lattice K can be chosen (a different, but equivalent, ${\mathbb{K}}$) so that ${\mathbb{N}}$ is a lower triangular matrix in Hermite normal form (HNF) [11]. (See section II-A of [8] for a brief introduction to HNF.) HNF is a lower-triangular canonical matrix form, where the entries below the diagonal are non-negative and strictly less than the diagonal entry in the same row. Code for converting integer matrices to Hermite Normal Form is available at https://github.com/msg-byu/symlib in the rational_mathemematics module.

A ${\bf{k}}$-point integration grid is the set of points of the lattice K that lie inside one unit cell (one fundamental domain) of the reciprocal lattice R. We refer to this finite subset of K as K_α (See figure 3; black dots are K, dots inside the blue parallelogram comprise K_α.) The number of points that lie within one unit cell of R is given by $| \det ({\mathbb{N}})| =n$.

How then does one generate these n points? If ${\mathbb{N}}$ is in HNF, then the diagonal elements of ${\mathbb{N}}$ are three integers, a, c, and f, such that a · c · f = n. A set of n translationally distinct⁴ points of the lattice K can be generated by taking integer linear combinations of the columns of ${\mathbb{K}}$:

$$\begin{eqnarray}&&\vec{k}=p{\vec{\kappa }}_{1}+q{\vec{\kappa }}_{2}+r{\vec{\kappa }}_{3},\end{eqnarray} \tag{ 2 }$$

where p, q, and r are nonnegative integers such that

$$\begin{eqnarray*}\begin{array}{rcl}0 & \leqslant & p\lt a\\ 0 & \leqslant & q\lt c\\ 0 & \leqslant & r\lt f.\end{array}\end{eqnarray*}$$

The n points generated this way will not generally lie inside the same unit cell, but they can be translated into the same cell by expressing them in 'lattice coordinates' (fractions of the columns of ${\mathbb{R}}$, instead of Cartesian coordinates) and then reducing the coordinates modulo 1 so that they all lie within the range $[0,1)$. This is illustrated by the dashed arrow in figure 4.

Zoom In Zoom Out Reset image size

Expressed as fractions of the lattice vectors of R, these four points are:

$$\begin{eqnarray*}\begin{array}{rcl}{\vec{k}}_{1} & = & (0,0)\\ {\vec{k}}_{2} & = & \left(0,\displaystyle \frac{1}{2}\right)\\ {\vec{k}}_{3} & = & \left(\displaystyle \frac{1}{2},-\displaystyle \frac{1}{4}\right)\mathop{\longrightarrow }\limits^{\mathrm{mod}\ 1}\left(\displaystyle \frac{1}{2},\displaystyle \frac{3}{4}\right)\\ {\vec{k}}_{4} & = & \left(\displaystyle \frac{1}{2},\tfrac{1}{4}\right).\end{array}\end{eqnarray*}$$

Initially, ${\vec{k}}_{3}$ is not in the same unit cell as the other three points; its first coordinate is not between 0 and 1. After reducing the first coordinate modulo 1, ${\vec{k}}_{3}$ moves to an equivalent position in the same unit cell as the other three points.

In summary, this first part of the algorithm generates n translationally distinct points and translates them all into the first unit cell of R. (It is not necessary to translate all the points into the first unit cell, but it is convenient to do so as a first step to translating them into the first Brillouin zone. The translation into the first Brillouin zone is less trivial and is discussed in section 4.)

In many cases the crystal will have some point group symmetries, and these can be exploited to reduce the number of ${\bf{k}}$-points where the energy eigenvalues and corresponding wavefunctions need to be evaluated. The grid is reduced by applying the point group symmetries⁵ of the crystal to each point in the grid. For example, in figure 5, the points connected by green arrows will be mapped onto one another by successive 90° rotations. These four symmetrically equivalent points lie on a 4-fold 'orbit' (as do the points marked by the red arrows). The point at the origin maps onto itself under all symmetry operations and has an orbit of length 1.

Zoom In Zoom Out Reset image size

For a grid containing N_k points and a group (of rotation and reflection symmetries) with N_G operations, a naive algorithm for identifying the symmetrically equivalent points and counting the length of each orbit would be as follows: For each point (${ \mathcal O }({N}_{k})$ loop), compare all rotations of that point (a loop of ${ \mathcal O }(| G| )$) to all other points (another ${ \mathcal O }({N}_{k})$ loop) to find a match; for a total computational complexity of ${ \mathcal O }({N}_{k}^{2}{N}_{G})$. The algorithm is shown in pseudocode in figure 6. N_G will never be larger than 48, but N_k may be as large as 50³ for extremely dense grids, so the ${N}_{k}^{2}$ complexity of this naive approach is undesirable. But using group theory concepts (see the appendix for details), we can construct a hash table for the points that reduces the complexity from ${ \mathcal O }({N}_{k}^{2}{N}_{G}$) to ${ \mathcal O }({N}_{k}{N}_{G})$ by eliminating the k_j loop in algorithm 1. The hash table makes a one-to-one association between the ordinal counter (the index) of each point and its coordinates.

Zoom In Zoom Out Reset image size

The three coefficients p, q, r in equation (2) can be conceptualized as the three 'digits' of a 3-digit mixed-radix number pqr or as the three numerals shown on an odometer with three wheels. The ranges of the values are 0 ≤ p < d₁, 0 ≤ q < d₂, and 0 ≤ r < d₃, where d₁, d₂, d₃ are the 'sizes' of the wheels, or in other words, the base of each digit. Then the mixed-radix number is converted to base 10 by

$$\begin{eqnarray}&&x=p\cdot {d}_{2}\cdot {d}_{1}+q\cdot {d}_{1}+r.\end{eqnarray} \tag{ 3 }$$

The total number of possible readings of the odometer is ${d}_{3}\cdot {d}_{2}\cdot {d}_{1}$. So it must be the case that the number of ${\bf{k}}$-points in the cell is $n={d}_{3}\cdot {d}_{2}\cdot {d}_{1}$. Each reading on the odometer is a distinct point of the n points that are contained in the reciprocal cell. Via equation (3) it is simple to convert a point given in 'lattice coordinates' as (p, q, r) to a base-10 number x. The concept of the hash table is to use this base-10 representation as the index in the hash table. Without the hash table, comparing two points is an ${ \mathcal O }({N}_{k})$ search because one point must be compared to every other point in the list to check for equivalency. But with the hash function, no search is necessary—one simply maps the point (p, q, r) to the index x of the equivalent ${\bf{k}}$-point in the hash table.

It is not generally the case that the coefficients p, q, r for every interior point of the unit cell obey conditions:

$$\begin{eqnarray}&&0\leqslant p\lt {d}_{1},\quad 0\leqslant q\lt {d}_{2},\quad 0\leqslant r\lt {d}_{3}\end{eqnarray} \tag{ 4 }$$

(Figure 4 shows an example where the interior points do not meet these conditions.) These conditions hold only for a certain choice of basis. That basis is found by transforming the matrix ${\mathbb{N}}$ in equation (1) into its Smith Normal Form [11], ${\mathbb{D}}={\mathbb{ANB}}$. By elementary row and column operations, represented by unimodular matrices ${\mathbb{A}}$ and ${\mathbb{B}}$, it is possible to transform ${\mathbb{N}}$ into a diagonal matrix ${\mathbb{D}}$, where each diagonal element divides the ones below it: ${d}_{11}| {d}_{22}| {d}_{33}$, and ${d}_{11}\cdot {d}_{22}\cdot {d}_{33}=n=\left|{\mathbb{N}}\right|$ (the notation $i| j$ means that i is divisible by j). As explained in the appendix (section A.1), when ${\mathbb{N}}$ is expressed in Smith normal form (SNF) and the interior points of the reciprocal cell are expressed as linear combinations of the grid generating vectors ${\mathbb{K}}$, then the coordinates (coefficients) of the interior points will obey equation (4). When these conditions are met, the hashing algorithm discussed above (in particular, equation (3)) becomes possible. This enables the ${ \mathcal O }({N}_{k})$ algorithm, shown in figure 7.

Zoom In Zoom Out Reset image size

For accurate DFT calculations, it is best if the energy eigenvalues (electronic bands) are evaluated at ${\bf{k}}$-points inside the first Brillouin zone, so our algorithm includes a step that finds the translationally equivalent grid points in the Brillouin zone. (In principle, the electronic structure $E({\bf{k}})$ should be the same in every unit cell, but numerically the periodicity of the electronic bands is only approximate, becoming less accurate for ${\bf{k}}$-points in unit cells farther from the origin.)

The first Brillouin zone of the reciprocal lattice is simply the Voronoi cell centered at the origin—all ${\bf{k}}$-points in the first Brillouin zone are closer to the origin than to any other lattice point. Conceptually, an algorithm for translating a ${\bf{k}}$-point of the integration grid into the first zone merely requires one to look at all translationally equivalent 'cousins' of the ${\bf{k}}$-point and select the one closest to the origin. But the number of translationally equivalent cousins is countably infinite, so in practice, the set of cousins must be limited only to those near the origin.

How can we select a set of cousins near the origin that is guaranteed to include the closest cousin? The key idea is illustrated in two-dimensions in figure 8. In three-dimensions, if the basis vectors of the reciprocal unit cell are as short as possible (the so-called Minkowski-reduced basis [12]), then the eight unit cells that all share a vertex at the origin must contain the Brillouin zone. In other words, the boundary of this union of eight cells is guaranteed to circumscribe the first Brilloun zone (i.e., the Voronoi cell containing the origin). A proof of this '8 cells' conjecture is given in the appendix (section A.1). The steps for moving ${\bf{k}}$-points into the Brillouin zone are as follows:

• 1.  
  Minkowski-reduce the reciprocal unit cell [12] (i.e., find the basis with the shortest basis vectors⁶ )
• 2.  
  For each ${\bf{k}}$-point in the reduced grid, find the translation-equivalent cousin in each of the eight unit cells that have a vertex at the origin.
• 3.  
  From these eight cousins, select the one closest to the origin.

Zoom In Zoom Out Reset image size

We have developed an algorithm that i) generates generalized regular ${\bf{k}}$-point grids, ii) reduces the grids by symmetry, and iii) maps the points of the reduced grids into the first Brillouin zone. Whereas the typical algorithm for generating and reducing ${\bf{k}}$-point grids scales quadratically with the number of ${\bf{k}}$-points, this algorithm scales linearly. The improved scaling becomes essential when one generates and symmetry reduces all combinatorically possible generalized regular grids at a given ${\bf{k}}$-point density, in order to select the one with the fewest number of reduced ${\bf{k}}$-points [6].

The algorithm is also useful because it relies primarily on integer-based operations, making it more robust than typical floating point-based algorithms that are prone to finite precision errors. Mapping the grid to the first Brillouin zone is more efficient due to a proof that limits the search for translationally equivalent ${\bf{k}}$-points to the eight unit cells having a vertex at the origin. The algorithm has been incorporated into version 6 of VASP [9].

We are grateful to Martijn Marsmann who found a bug in an early version of the code that spoiled the ${ \mathcal O }({N}_{k})$ scaling and who implemented the algorithm in VASP, version 6. Tim Mueller also contributed by pointing out a simpler approach to generating all the ${\bf{k}}$-points interior to the unit cell of R than we originally proposed. GLWH, WSM, and JJJ are grateful for financial support from the Office of Naval Research (MURI N00014-13-1-0635).

A.1. Proof limiting Brillouin zone location

Given a point x in the space, we will use the term cousin for a point $x^{\prime}$ which differs from x by an element of the lattice—i.e., a coset representative or a lattice-translation equivalent point.

Let R be a basis. Let U_R denote the union of 2^d basis cells around the origin—the set of points which are expressible in terms of the basis R with all coefficients having absolute value < 1. Let V denote the Voronoi cell (Brillouin zone)—the set of all points in the space which are closer to the origin than any other lattice point. Note that U_R depends on the basis R, but V depends only on the lattice. Note also that both U_R and V are convex sets.

We claim (in two and three dimensions) that if R is a Minkowski basis, then $V\subseteq {U}_{R}$. We shall argue by contrapositive— if $V \nsubseteq {U}_{R}$, then the basis is not Minkowski reduced.

If $V \nsubseteq {U}_{R}$ then V must intersect the boundary of U_R, so there exist points on the boundary of U_R which are closer to the origin than to any other lattice points. Equivalently, those points are closer to the origin than any of their cousins.

Note that among the cousins of any point on the boundary of U_R, there is always a closest to the origin. But usually points on the boundary will have closer cousins in the interior. But if $V \nsubseteq {U}_{R}$ there must be points on the boundary which have no closer cousins in the interior of U_R. In other words, there are points (at least one) on the boundary such that all of its cousins in the interior of U_R are farther from the origin.

A.1.1. 2D argument

Let ${\vec{r}}_{1}$ and ${\vec{r}}_{2}$ be basis elements of R. Assuming that $V \nsubseteq {U}_{R}$ there must be a point x on the boundary of U_R whose cousins are all farther from the origin than x.

Without loss of generality (re-label the basis if necessary), we may express one of the bounding edges of U_R as $x={\vec{r}}_{1}+\lambda {\vec{r}}_{2}$ where λ ∈ [−1, 1]. One of its interior cousins is $x^{\prime} =\lambda {\vec{r}}_{2}$, which is illustrated in figure A1. We have (since $x^{\prime}$ must be farther from the origin)

$$\begin{eqnarray*}\begin{array}{rcl}{x}^{2} & \lt & x{{\prime} }^{2}\\ {\left({\vec{r}}_{1}+\lambda {\vec{r}}_{2}\right)}^{2} & \lt & {\left(\lambda {\vec{r}}_{2}\right)}^{2}\\ {\vec{r}}_{1}^{2}+2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}+{\lambda }^{2}{\vec{r}}_{2}^{2} & \lt & {\lambda }^{2}{\vec{r}}_{2}^{2}\\ {\vec{r}}_{1}^{2} & \lt & -2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}\end{array}\end{eqnarray*}$$

Since the expression on the left-hand side is greater than zero, the expression on the right-hand side must be also and taking the absolute value of both sides does not change the inequality:

$$\begin{eqnarray*}&&| {\vec{r}}_{1}^{2}| \lt | -2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}| \Rightarrow | {\vec{r}}_{1}{| }^{2}\lt 2| \lambda | | {\vec{r}}_{1}\cdot {\vec{r}}_{2}| .\end{eqnarray*}$$

Considering the worst case scenario of λ = 1 gives

$$\begin{eqnarray}&&\displaystyle \frac{| {\vec{r}}_{1}| }{2}\lt \displaystyle \frac{| {\vec{r}}_{1}\cdot {\vec{r}}_{2}| }{| {\vec{r}}_{1}| },\end{eqnarray} \tag{ 5 }$$

which violates the condition of a Minkowski basis $| {\vec{r}}_{1}\cdot {\vec{r}}_{2}| /| {\vec{r}}_{1}| \lt | {\vec{r}}_{1}| /2$. The remaining three boundaries are similar to the one just considered, the only differences being permutations of the basis elements ${\vec{r}}_{1}$ and ${\vec{r}}_{2}$ and possibly changes of sign. When applying the same reasoning to the other edges we arrive at the same contradiction. Hence, the points on the boundary of U_R are closer to the origin than interior cousins, $V \nsubseteq {U}_{R}$, only when the basis R is not Minkowski reduced. If R is Minkowski reduced, all points on the bounday of U_R have interior cousins that lie closer to the origin and $V\subseteq {U}_{R}$.

Zoom In Zoom Out Reset image size

A.1.2. 3D argument

Let ${\vec{r}}_{1},{\vec{r}}_{2}$, and ${\vec{r}}_{3}$ be the basis elements of R, and suppose (relabeling the basis vectors if necessary) that $x={\vec{r}}_{1}+\lambda {\vec{r}}_{2}+\delta {\vec{r}}_{3}$ (where λ and δ are elements of [−1, 1]) is a point on the boundary of U_R which is closer to the origin than are its interior cousins.

One of those cousins is a plane through the origin $x^{\prime} =\lambda {\vec{r}}_{2}+\delta {\vec{r}}_{3}$. The boundary and cousin planes are shown in figure A2. Thus

$$\begin{eqnarray*}\begin{array}{rcl}{x}^{2} & \lt & x{{\prime} }^{2}\\ {\left({\vec{r}}_{1}+\lambda {\vec{r}}_{2}+\delta {\vec{r}}_{3}\right)}^{2} & \lt & {\left(\lambda {\vec{r}}_{2}+\delta {\vec{r}}_{3}\right)}^{2}.\end{array}\end{eqnarray*}$$

Simplifying this expression gives

$$\begin{eqnarray}&&{\vec{r}}_{1}^{2}\lt -2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}-2\delta {\vec{r}}_{1}\cdot \end{eqnarray} \tag{ 6 }$$

Since the expression on the left-hand side is greater than zero, the expression on the right-hand side must be also and taking the absolute value of both sides does not change the inequality:

$$\begin{eqnarray}&&| {\vec{r}}_{1}^{2}| \lt | -2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}-2\delta {\vec{r}}_{1}\cdot {\vec{r}}_{3}| \Rightarrow | {\vec{r}}_{1}{| }^{2}\lt | 2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}+2\delta {\vec{r}}_{1}\cdot {\vec{r}}_{3}| \end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

To simplify the right-hand side of equation (7) the triangle inequality is used, making it more likely that the inequality is satisfied:

$$\begin{eqnarray}\begin{array}{rcl}| {\vec{r}}_{1}{| }^{2} & \lt & | 2\lambda {\vec{r}}_{1}\cdot {\vec{r}}_{2}\,+\,2\delta | {\vec{r}}_{1}\cdot {\vec{r}}_{3}| \\ {\vec{r}}_{1}{| }^{2} & \lt & 2| \lambda | | | {\vec{r}}_{1}\cdot {\vec{r}}_{2}| +2| \delta | | {\vec{r}}_{1}\cdot {\vec{r}}_{3}| \end{array}\end{eqnarray} \tag{ 8 }$$

Since the expression in equation (8) does not depend on the sign of λ or δ, we can restrict λ and δ to positive values within [0,1] without loss of generality. Consider now another cousin that lies within U_R and on the same plane as $x^{\prime}$: $x^{\prime\prime} =(\lambda -1){\vec{r}}_{2}+(\delta -1){\vec{r}}_{3}$. Repeating the same process for $x^{\prime}$ with $x^{\prime\prime}$ gives

$$\begin{eqnarray}&&| {\vec{r}}_{1}{| }^{2}\lt 2| \lambda -1| | {\vec{r}}_{1}\cdot {\vec{r}}_{2}| +2| \delta -1| | {\vec{r}}_{1}\cdot {\vec{r}}_{3}| \end{eqnarray} \tag{ 9 }$$

Combining equations (8) and (9) gives

$$\begin{eqnarray}\begin{array}{rcl}| {\vec{r}}_{1}{| }^{2} & \lt & (| \lambda | +| \lambda -1| )| {\vec{r}}_{1}\cdot | {\vec{r}}_{2}| +(| \delta | +| \delta -1| )| {\vec{r}}_{1}\cdot {\vec{r}}_{3}| | \\ {\vec{r}}_{1}| & \lt & \displaystyle \frac{| {\vec{r}}_{1}\cdot | {\vec{r}}_{2}| }{| {\vec{r}}_{1}| }+\displaystyle \frac{| {\vec{r}}_{1}\cdot | {\vec{r}}_{3}| }{| {\vec{r}}_{1}| }\end{array}\end{eqnarray} \tag{ 10 }$$

Assuming $\{{\vec{r}}_{1},{\vec{r}}_{2},{\vec{r}}_{3}\}$ forms a Minkowski basis, and plugging in the largest possible values for all quantities under this assumption on the right-hand side of equation (10) gives the contradiction $| {\vec{r}}_{1}| \lt | {\vec{r}}_{1}|$. The remaining seven bounding planes are similar to the one just considered, the only differences being permutations of the basis elements ${\vec{r}}_{1}$, ${\vec{r}}_{2}$, and ${\vec{r}}_{3}$ and changes of sign. We arrive at the same contradiction when applying the same reasoning to the other planes. Hence, the points on the boundary of U_R are closer to the origin than interior cousins, $V \nsubseteq {U}_{R}$, only when the basis R is not Minkowski reduced. If R is Minkowski reduced, all points on the boundary of U_R have interior cousins that lie closer to the origin and $V\subseteq {U}_{R}$.

A.2. Groups, matrices, and lattices in smith normal form

The discussion below is limited to three-dimensions though the arguments easily generalize to higher dimensions. The purpose of the discussion below is to help the reader make the connection between groups and integer matrices. The Smith Normal Form is a key concept to make this connection.

In this discussion, we show that we can associate a single, finite group with the lattice sites within one tile (i.e., one unit cell) of a superlattice. In our application, this tile is the unit cell of the grid generating vectors and the superlattice is the reciprocal cell. The association between the group and the lattice sites is a homomorphism that maps each lattice site to an element of the group. If two points are translationally equilavent (same site but in two different tiles) they will map to the same element of the group. This homomorphism is the key ingredient to constructing the hash function (see equation (3)) that enables a perfect hash table where points are listed consecutively, from 1 to N. In what follows, we explain in detal how this association is made, i.e., we detail how one finds this homomorphism.

A.2.1. Groups in smith normal form

Begin with the simplest case. Let ${\mathbb{N}}$ be a non-singular 3 × 3 matrix of integers. Its columns represent the basis for a subgroup ${{ \mathcal L }}_{N}$ of the group Z³ (where Z is the set of all integers, and the group operation is addition). The two latttices whose symmetries are represented by these two groups are the 'simple cubic' lattice of all points with all integer coordinates and its superlattice⁷ whose basis is given by the columns of ${\mathbb{N}}$. Since Z³ and its subgroups are Abelian, we know that all the subgroups are normal so there exists a quotient group $G={{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$, and that group is finite.

Note that the cosets which form the elements of that quotient group are simply the distinct translates of the lattice ${{ \mathcal L }}_{N}$ within Z³. In fact, each coset has exactly one representative in each unit cell, so the order of G is equal to the volume of a unit cell (the absolute value of the determinant of ${\mathbb{N}}$). Since the quotient group G is finite, and Abelian, it must be a direct sum of cyclic groups (by the Fundamental theorem of Finite Abelian Groups).

One canonical form for direct sums of groups is called Smith Normal Form, where the direct summands are ordered so that each summand divides the next. In other words, $G\cong {{\bf{Z}}}_{{m}_{1}}\oplus {{\bf{Z}}}_{{m}_{2}}\oplus \cdots \oplus \,{{\bf{Z}}}_{{m}_{k}}$ where ${m}_{1}| {m}_{2}| ...{m}_{k-1}| {m}_{k}$ and (of course) $\prod {m}_{i}=| G|$. Any finite Abelian group can be uniquely written in this form. (Isomorphic groups will yield the same 'invariant factors' m₁, m₂,...,m_k when written in this form.)

Note that, since $G={{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$, there must be a homomorphism from Z³ onto G, having ${{ \mathcal L }}_{N}$ as its kernel. In other words, ${{ \mathcal L }}_{N}=\{p\in {{\bf{Z}}}^{3}\,:\,\psi (p)=0\}$. Our task is to find the direct-sum representation of the quotient group ${{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$, and also to find the homomorphism ψ which maps the points of Z³ onto the group (in such a way that ψ(p) = 0 iff $p\in {{ \mathcal L }}_{N}$). This allows us to work with the elements of the group as proxies for the ${\bf{k}}$-points inside the reciprocal cell.

A.2.2. Matrices in smith normal form

There is a useful connection between the SNF for Abelian groups and the SNF for integer matrices. As the reader may infer, the SNF form of the basis matrix ${\mathbb{N}}$ effectively tells one how to represent the quotient group ${{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$ as a direct sum of cyclic groups in Smith Normal Form, and, as shown in the following section, the row operations used to create the SNF of ${\mathbb{N}}$ give the homomorphism ψ suggested above.

A.2.3. The connection between SNF groups and SNF matrices

In the matrix case, since the operations are elementary row and column operations, we have ${\mathbb{D}}={\mathbb{ANB}}$ where ${\mathbb{A}}$ and ${\mathbb{B}}$ are integer matrices with determinant ±1 representing the accumulated row operations and column operations respectively. The matrix ${\mathbb{D}}$ is completely determined by ${\mathbb{N}}$, but the matrices ${\mathbb{A}}$ and ${\mathbb{B}}$ depend on the algorithm used to arrive at the Smith Normal Form of ${\mathbb{N}}$. A different implementation might yield ${\mathbb{D}}={\mathbb{A}}^{\prime} {\mathbb{NB}}^{\prime}$ (same ${\mathbb{N}}$ and same ${\mathbb{D}}$, but different ${\mathbb{A}}$ and ${\mathbb{B}}$).

Note that, since ${\mathbb{B}}$ represents elementary column operations, the product ${\mathbb{NB}}$ simply represents a change of basis from ${\mathbb{N}}$ to a new basis ${\mathbb{N}}^{\prime} ={\mathbb{NB}}$. In other words, the columns of ${\mathbb{N}}^{\prime}$ are still a basis for ${{ \mathcal L }}_{N}$. But the new basis has the property that ${\mathbb{AN}}^{\prime} ={\mathbb{D}}$. That means that every element $\vec{w}={\mathbb{N}}^{\prime} \vec{z}$ of ${{ \mathcal L }}_{N}$ (where $\vec{z}$ is some element of ${{\mathbb{Z}}}^{3}$) will satisfy the equation ${\mathbb{A}}\vec{w}={\mathbb{D}}\vec{z}=\left(\begin{array}{c}{{\mathbb{D}}}_{11}{z}_{1}\\ {{\mathbb{D}}}_{22}{z}_{2}\\ {{\mathbb{D}}}_{33}{z}_{3}\end{array}\right)$. In other words, ${\mathbb{A}}\vec{w}$ will be a vector whose entries are multiples of the corresponding diagonal entries in ${\mathbb{D}}$.

To put it another way, define $*$ to be the operation that maps $\vec{x}=\left(\begin{array}{c}{x}_{1}\\ {x}_{2}\\ {x}_{3}\end{array}\right)$ in ${{\mathbb{Z}}}^{3}$ to ${\vec{x}}^{* }={\left(\begin{array}{c}{x}_{1}(\mathrm{mod}{{\mathbb{D}}}_{11})\\ {x}_{2}(\mathrm{mod}{{\mathbb{D}}}_{22})\\ {x}_{3}(\mathrm{mod}{{\mathbb{D}}}_{33})\end{array}\right)}^{T}$.

Then we have shown $\vec{w}\in {{ \mathcal L }}_{N}$ iff ${({\mathbb{A}}w)}^{* }=(0,0,0)$ (the zero-element in the group ${G}_{0}={{\mathbb{Z}}}_{{{\mathbb{D}}}_{11}}\oplus {{\mathbb{Z}}}_{{{\mathbb{D}}}_{22}}\oplus {{\mathbb{Z}}}_{{{\mathbb{D}}}_{33}}$.

That suggests we let $\psi {(\vec{w})=({\mathbb{A}}\vec{w})}^{* }$, a homomorphism from Z³ onto the direct-sum G₀. Then, since that homomorphism is easily shown to be onto, and its kernel is ${{ \mathcal L }}_{N}$, we see (by the First Isomorphism theorem of group theory) that ${G}_{0}\cong {{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$, and ψ is precisely the homomorphism we sought.

Thus we have connected the two versions of SNF. The matrix algorithm provides the SNF description of the quotient group by the diagonal entries in ${\mathbb{D}}$, and the transition matrix ${\mathbb{A}}$ provides the homomorphism which maps the parent lattice onto the group.

A.2.3.1. An example

Let ${\mathbb{N}}=\left(\begin{array}{ccc}1 & 2 & -1\\ 1 & 4 & -3\\ 0 & 2 & 4\end{array}\right)$. This describes a lattice ${{ \mathcal L }}_{N}$ which contains the points ${\vec{p}}_{1}=\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right),{\vec{p}}_{2}=\left(\begin{array}{c}2\\ 4\\ 2\end{array}\right)$, and ${\vec{p}}_{3}=\left(\begin{array}{c}-1\\ -3\\ 4\end{array}\right)$, and all the points which are integer linear combinations of those three points. The matrix ${\mathbb{N}}$ has determinant 12, which must be the volume of each lattice tile—and it is also the order of the quotient group ${{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}$.

Using the SNF algorithm to diagonalize this basis matrix, we find ${\mathbb{D}}={\mathbb{ANB}}$ where ${\mathbb{D}}=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 6\end{array}\right)$, with ${\mathbb{A}}=\left(\begin{array}{ccc}0 & 1 & 0\\ 0 & 0 & 1\\ 1 & -1 & -2\end{array}\right)$ and ${\mathbb{B}}=\left(\begin{array}{ccc}1 & 7 & 11\\ 0 & -1 & -2\\ 0 & 1 & 1\end{array}\right)$.

Thus we now know that the quotient group is $G={{\bf{Z}}}^{3}/{{ \mathcal L }}_{N}\cong {{\bf{Z}}}_{1}\oplus {{\bf{Z}}}_{2}\oplus {{\bf{Z}}}_{6}\cong {{\bf{Z}}}_{2}\oplus {{\bf{Z}}}_{6}$.

Further, from the matrix ${\mathbb{A}}$, we may obtain the homomorphism projecting Z³ onto the quotient group, with kernel ${{ \mathcal L }}_{N}$. If $\vec{w}=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)$ then ${\mathbb{A}}\vec{w}=\left(\begin{array}{c}y\\ z\\ x-y-2z\end{array}\right)$ and thus

$$\begin{eqnarray*}\begin{array}{rcl}\psi (\vec{w}) & = & {\left({\mathbb{A}}\vec{w}\right)}^{* }\\ & = & {\left(\begin{array}{c}y(\mathrm{mod}1)\\ z(\mathrm{mod}2)\\ x-y-2z(\mathrm{mod}6)\end{array}\right)}^{T}\\ & = & \left(z\,(\mathrm{mod}\,2),x+5y+4z\,(\mathrm{mod}\,6)\right)\end{array}\end{eqnarray*}$$

(noting that anything mod 1 is zero).

Note that this homomorphism provides a different, but convenient, way to describe the superlattice. Since ${{ \mathcal L }}_{N}$ is the kernel of ψ, it is comprised of the points (x, y, z) ∈ Z³ which satisfy the simultaneous congruences z ≡ 0 (mod 2) and $x+5y+4z\equiv 0$ (mod 6). We note that all three basis points p₁, p₂ and p₃ satisfy these congruences, and thus so will all their integer linear combinations (all points in ${{ \mathcal L }}_{N}$).

A.2.3.2. Algorithmic variation

In the example we computed above, a different application of the SNF matrix algorithm, with the same ${\mathbb{N}}$, might have yielded the same diagonal matrix ${\mathbb{D}}=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 6\end{array}\right)$, but different ${\mathbb{A}}=\left(\begin{array}{ccc}1 & 0 & 0\\ -5 & 3 & 1\\ -2 & 2 & 1\end{array}\right)$ and ${\mathbb{B}}=\left(\begin{array}{ccc}0 & -1 & 2\\ 0 & 0 & 1\\ -1 & -1 & 4\end{array}\right)$, which would change the homomorphism to $(x,y,z)\mapsto (-5x+3y+z$ (mod 2), $-2x+2y+z$ (mod 6)) = $(x+y+z\,(\mathrm{mod}\,2)$, $4x+2y+z$ (mod 6)). The new homomorphism is different, since (1, 0, 1) ↦ (0, 5) now, where previously (1, 2, 3) ↦ (1, 5) (for example), but the kernel is the same. In fact the two homomorphisms are related via an automorphism of the group G.

A.2.4. Non-integer lattices

Now, what about the more complicated situation, where ${\mathbb{N}}$ represents a (possibly HNF) matrix describing the change from some lattice other than the simple integer lattice Z³ to one of its subgroups (superlattice)?

Then we have a basis ${\mathbb{V}}$ and lattice ${{ \mathcal L }}_{V}$ and a basis ${\mathbb{W}}={\mathbb{VN}}$ for a (super) lattice ${{ \mathcal L }}_{W}$. Again, the quotient group $G={{ \mathcal L }}_{V}/{{ \mathcal L }}_{W}$ is Abelian of order $| \det ({\mathbb{N}})|$. Again, G is a direct sum of cyclic groups corresponding to the diagonal entries of ${\mathbb{D}}={\mathbb{ANB}}$ (where ${\mathbb{D}}$ is the SNF of ${\mathbb{N}}$).

The only difference here is that the homomorphism ψ provided by ${\mathbb{A}}$ must depend on the basis ${\mathbb{V}}$ (which might even be irrational). Every point in ${{ \mathcal L }}_{V}$ has the form $\vec{x}={\mathbb{V}}\vec{w}$ where $\vec{w}$ is a column of integers. Then $\psi (\vec{x})={\mathbb{A}}\vec{w}$ (modded by the corresponding entries from ${\mathbb{D}}={\mathbb{ANB}}$). We could write it as $\psi (\vec{x})={\left({{\mathbb{AV}}}^{-1}\vec{x}\right)}^{* }$ (with the entries appropriately modularly reduced and transposed to a horizontal vector).

A.2.5. Example: general (non-integer) lattices

Suppose ${{ \mathcal L }}_{V}$ is the lattice defined by (columns of) the basis matrix ${\mathbb{V}}=\left(\begin{array}{ccc}1 & 1/2 & 0\\ 0 & \sqrt{3}/2 & 0\\ 0 & 0 & 2\end{array}\right)$, and ${{ \mathcal L }}_{W}$ is the subgroup lattice defined by the basis marix ${\mathbb{W}}={\mathbb{VN}}$ where ${\mathbb{N}}=\left(\begin{array}{ccc}4 & 2 & 2\\ 2 & 2 & 2\\ 4 & 0 & 4\end{array}\right)$. In other words, one basis for ${{ \mathcal L }}_{W}$ is given by the columns of ${\mathbb{W}}=\left(\begin{array}{ccc}5 & 3 & 3\\ \sqrt{3} & \sqrt{3} & \sqrt{3}\\ 8 & 0 & 8\end{array}\right)$.

Reducing ${\mathbb{N}}$ to SNF yields

$$\begin{eqnarray*}{\mathbb{D}}=\left(\begin{array}{ccc}2 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 4\end{array}\right)=\left(\begin{array}{ccc}1 & 0 & -1\\ 1 & -1 & -1\\ -6 & 4 & 5\end{array}\right){\mathbb{N}}\left(\begin{array}{ccc}-2 & -3 & -2\\ 2 & 1 & 1\\ 1 & 1 & 1\end{array}\right).\end{eqnarray*}$$

Thus our quotient group is $G={{ \mathcal L }}_{V}/{{ \mathcal L }}_{W}\cong {{\bf{Z}}}_{2}\oplus {{\mathbb{Z}}}_{2}\oplus {{\mathbb{Z}}}_{4}$ and ${\mathbb{A}}=\left(\begin{array}{ccc}1 & 0 & -1\\ 1 & -1 & -1\\ -6 & 4 & 5\end{array}\right)$ so

$$\begin{eqnarray*}{{\mathbb{AV}}}^{-1}=\left(\begin{array}{ccc}1 & -\sqrt{3}/3 & -1/2\\ 1 & -\sqrt{3} & -1/2\\ -6 & 14\sqrt{3}/3 & 5/2\end{array}\right),\end{eqnarray*}$$

which provides our homomorphism $\psi (\vec{x})={\left({{\mathbb{AV}}}^{-1}x\right)}^{* }$ from ${{ \mathcal L }}_{V}$ onto G.

If we let $\vec{x}=\left(\begin{array}{c}2\\ \sqrt{3}\\ 2\end{array}\right)$ which is an element of ${{ \mathcal L }}_{V}$ but not of ${{ \mathcal L }}_{W}$, then ${{\mathbb{AV}}}^{-1}x=\left(\begin{array}{c}0\\ -2\\ 7\end{array}\right)$ and $\psi (\vec{x})=(0,0,3)\in G$ (after reducing the elements modulo 2, 2 and 4 respectively). On the other hand, if we let $\vec{x}=\left(\begin{array}{c}7\\ \sqrt{3}\\ 8\end{array}\right)$, which is an element of ${{ \mathcal L }}_{W}$ (the kernel), then ${{\mathbb{AV}}}^{-1}\vec{x}=\left(\begin{array}{c}2\\ 0\\ -8\end{array}\right)$ and so $\psi (\vec{x})=(0,0,0)$, and $\vec{x}$ is in the group.

By this function ψ, the elements of ${{ \mathcal L }}_{V}$ are all mapped to elements of the group G and, in particular, the elements of ${{ \mathcal L }}_{W}$ are mapped to the zero element of the group. Stated in terms of the cosets, the entire set ${{ \mathcal L }}_{W}$ is mapped to the zero element of the group G, and each of the distinct translates of ${{ \mathcal L }}_{W}$ (within ${{ \mathcal L }}_{V}$) gets mapped to a different element of the group. We might think of this as decorating or labeling the elements of ${{ \mathcal L }}_{V}$ in a periodic manner, using ${{ \mathcal L }}_{W}$ to define the period, and using the elements of the group G as the labels.