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

We develop an algorithm for i) computing generalized regular k -point grids, ii) reducing the grids to their symmetrically distinct points, and iii) mapping the reduced grid points into the Brillouin zone. The algorithm exploits the connection between integer matrices and finite groups to achieve a computational complexity that is linear with the number of k -points. The favorable scaling means that, at a given k -point density, all possible commensurate grids can be generated (as suggested by Moreno and Soler) and quickly reduced to identify the grid with the fewest symmetrically unique k -points. These optimal grids provide significant speed-up compared to Monkhorst–Pack k -point grids; they have better symmetry reduction resulting in fewer irreducible k -points at a given grid density. The integer nature of this new reduction algorithm also simplifies issues with finite precision in current implementations. The algorithm is available as open source software.


I. INTRODUCTION
Codes that simulate materials using density functional theory (DFT) create 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 regions of the 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 Fig. 1 where we compare the convergence of an insulator with a metal.
The poor convergence of the electronic energy means that DFT codes must generate extremely dense grids for a target accuracy of several meV/atom 2 .It is common practice to increase computational efficiency by using the crystal's symmetry to reduce the number of electronic structure evaluations: grid points that are symmetrically equivalent have the same eigenenergies, so only one has to be evaluated.
In this brief report, an algorithm for generating, and subsequently symmetry-reducing, k-point grids is explained.This algorithm relies heavily on concepts used in the derivative structure enumeration algorithm of the enumlib code. 4The algorithm has been implemented in an open-source code available at https://github.com/msg-byu/kgridGen.The algorithm has been incorporated in version 6 of the VASP code. 5n 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) 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 1.Total energy error vs. k-point density for the cases of silicon and aluminum.Silicon does not have a Fermi surface so there is no discontinuity in the occupied bands; convergence is super-exponential.(See the discussion of example 1 in Ref. 3.) In contrast, the total energy of aluminum converges very slowly, and the convergence is quite erratic.For typical target accuracies in the total energy, around 10 −3 eV/atom, metals require 10-50 times more k-points than semiconductors.
of the integration grid increases), and 2) 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 (in other words, to enable a k-point generation method in the same spirit as that of Ref. 1 but for which the grid generation can be done on-the-fly).
Example of the integer relationship between reciprocal lattice vectors R and the grid generating vectors K.In the picture, the grid generating vectors, κ1 and κ2, the columns of K, define a lattice of points, four of which are inside the unit cell of R. Note that in the most general case, the relationship between the two lattices, N need not be diagonal (as it is for Monkhorst-Pack 6 k-point grids.)

II. GENERATING GRIDS
Every uniform sampling of a Brillouin zone can be expressed through the simple integer relationship In Eq. ( 1), R, K, and N are 3 × 3 matrices; the columns of R are the reciprocal lattice vectors, the columns of K are the k-point grid generating vectors, and N is a nonsingular integer matrix.Simply, N describes the integer linear combination of vectors of K that is equivalent to R. One example of the concept, illustrated in two dimensions, is shown in Fig. 2. We refer to the infinite lattice of points defined by integer linear combinations of R as R, and the lattice of points defined by K as the lattice K.
With no loss of generality, a new basis for the grid K can be chosen (a different, but equivalent, K) so that N is a lower triangular matrix in Hermite normal form.(HNF; see Sec.II-A of Ref. 4 for a brief introduction to Hermite normal form.HNF is a lower-triangular canonical matrix form, where the non-zero, off-diagonal terms are nonnegative and less than the diagonal entry in the same row.Code for generating the HNF [and SNF, see below], is available at https://github.com/msg-byu/symlib.) The 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 Fig. 2; 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(N)| = n.How then does one generate these n points?If N is in HNF, then the diagonal elements of N are three integers (call them a, c, and f ) such that a 3. An example of generating the points of K (black lattice) that lie within one unit cell (blue parallelogram) of the lattice R (blue lattice).The lattice K is generated by the basis { κ1, κ2} (columns of K).The four points of Kα are generated by k = m1 κ1 + m2 κ2, where 0 ≤ m1 < 2, 0 ≤ m2 < 2. Note that the upper limits of m1 and m2 are the diagonals of N when it is expressed in HNF.
be generated by taking integer linear combinations of the basis of K (columns of the matrix K): where p, q, and r are nonnegative integers such that 0 ≤ p < a and κ 1 , κ 2 , and κ 3 are the columns of K.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 R, instead of Cartesian coordinates) and then modding the coordinates by 1 so that they all lie within the range [0, 1).This is illustrated by the dashed arrow in Fig. 3. Expressed as fractions of the lattice vectors of R, these four points are: Initially, k 3 is not in the same unit cell as the other three points; its first coordinate is not between 0 and r 1 r 2 FIG. 4.An example of symmetry reducing a grid.The reciprocal unit cell is a square.This example assumes that the wavefunctions have square symmetry as well (the D4 group, 8 operations).The example grid is a 3 × 3 sampling of the reciprocal unit cell.The point at (0, 0) is not equivalent to any of the other eight points.There are two sets of equivalent points, each set with 4 points in the orbit, connected by red and green arrows, respectively.The points marked by red arrows are equivalent under horizontal, diagonal, and vertical reflections about the center of the square.The green-marked points are equivalent by 90 • rotations.Thus the nine points are reduced (or "folded") into three symmetrically distinct points.
1.After modding the first coordinate by 1, 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 first unit cell.It is, however, good practice to translate all the points into the first Brillouin zone, and a method for doing this is discussed in Sec.IV.)

III. SYMMETRY REDUCTION OF THE GRID
In many cases the crystal will have some point group symmetries, and these can be exploited to reduce the number of k-points where the energy eigenvalues and corresponding wavefunctions need to be evaluated.The grid is reduced by applying the symmetry operations of the crystal 8 to each point in the grid.For example, in Fig. 4, 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.
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 (O(N k ) loop), compare all rotations of that point (a loop of O(|G|)) to all other points (another O(N k ) loop) to find a match; for a total computational complexity of O(N 2 k N G ) (where N G is the number of rotation and reflection symmetries).In pseudocode, this algorithm would be the following: uniqueCount ++ add k i to list of unique points N G will never be larger than 48, but N k may be as large as 50 3 for extremely dense grids, so the N 2 k 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 ) by eliminating the k j loop in the above algorithm.The hash table makes a one-to-one association between the ordinal counter (i.e., the index) for each point and its coordinates.
The three coefficients p, q, r in Eq. ( 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 1 , 0 ≤ q < d 2 , and 0 ≤ r < d 3 , where d 1 , d 2 , d 3 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 as (3) The total number of possible readings of the odometer is Each reading on the odometer is a distinct point of the n points that are contained in the reciprocal cell.Via Eq. (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 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, mapping (p, q, r) to x only requires a single comparison.
It is not generally the case that the coefficients p, q, r for every interior point of the unit cell obey conditions: (Fig. 3 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 N in Eq. ( 1 Algorithm 2 #point and all its symmetry #equivalent points has already been #indexed At the end of the algorithm, the array Wt will be a list of the weights for each symmetrically unique point, and the index of each unique point in K will be stored in the array First.

IV. MOVING POINTS INTO THE FIRST BRILLOUIN ZONE
For accurate DFT calculations, it is best if the energy eigenvalues (electronic bands) are evaluated at k-points inside the first Brillouin zone, so our algorithm includes a FIG. 5.In 2D, the Brillion zone, shown in blue, is a subset of the union of 4 basis cells around the origin, shown in red, when the basis vectors are chosen to be as short as possible (the so-called Minkowski basis).If the basis is not Minkowski reduced, regions of the Brillouin zone may lie outside the union of the 4 basis cells, which is depicted by the cell in green.
step that finds the translationally equivalent grid points in the Brillouin zone.(In principle, the electronic structure should be the same in every unit cell, but numerically the periodicity of the electronic bands is only approximate, becoming less accurate for 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 k-points in the first Brillouin zone are closer to the origin than to any other lattice point.Conceptually, an algorithm for translating a k-point of the integration grid into the first zone merely requires one to look at all translationally equivalent "cousins" of the k-point and select the one closest to the origin.The number of translationally equivalent cousins is countably infinite, of course, 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 Fig. 5.In threedimensions, if the basis vectors of the reciprocal unit cell are as short as possible (the so-called Minkowski-reduced basis), then the eight unit cells that all share a vertex at the origin must contain every point that is closer to the origin than any other lattice point.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.The steps for moving k-points into the Brillouin zone are as follows: 1. Minkowski-reduce the reciprocal unit cell 9 (i.e., find the basis with the shortest basis vectors 10 ) 2. For each 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.

V. CONCLUSION
We have developed an algorithm for generating generalized regular k-point grids, symmetry reducing those grids, and mapping each point of the reduced grid into the first Brillouin zone.This algorithm scales linearly with the total number of k-points.Mapping the grid to the Brillouin zone benefits from a proof that limits the search for translationally equivalent k-points to the eight unit cells having a vertex at the origin.For dense grids, the application of the algorithm can reduce time and memory requirements for DFT calculations.The algorithm is also useful because it relies primarily on integerbased operations, making it more robust than typical floating point-based algorithms that are prone to finite precision errors.This algorithm has been incorporated into version 6 of VASP. 5

VI. ACKNOWLEDGMENTS
We are grateful to Martijn Marsmann who found a bug in an early version of the code that spoiled the 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 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).

VII. APPENDIX A. Proof limiting Brillouin zone location
Given a point x in the space, we will use the term cousin for a point x which differs from x by an element of the lattice-i.e., a coset representative or a latticetranslation equivalent point.Brillouin 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 to 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 ⊆ U R .We shall argue by contrapositive-if V U R , then the basis is not Minkowski reduced.
If V 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 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.

2D argument
Let r 1 and r 2 be basis elements of R. Assuming that V 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 = r 1 + λ r 2 where λ ∈ [−1, 1].One of its interior cousins is x = λ r 2 , which is illustrated in Fig. 6.We have (since x must be farther from the origin) 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: Considering the worst case scenario of λ = 1 gives The remaining three boundaries are similar to the one just considered, the only differences being permutations of the basis elements r 1 and 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 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 ⊆ U R .

3D argument
Let r 1 , r 2 , and r 3 be the basis elements of R, and suppose (relabeling the basis vectors if necessary) that x = r 1 + λ r 2 + δ 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 = λ r 2 + δ r 3 .The boundary and cousin planes are shown in Fig. 7. Thus

Simplifying this expression gives
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: 7) Since using the triangle inequality makes it more likely that the inequality in Eq. 7 is true, we can use it to simplify the right-hand side: Since the expression in Eqn. 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 : x = (λ − 1) r 2 + (δ − 1) r 3 .Repeating the same process for x with x gives Combining Eqns.8 and 9 gives Assuming { r 1 , r 2 , 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 Eqn. 10 gives the contradiction | r 1 | < | r 1 |.The remaining seven bounding planes are similar to the one just considered, the only differences being permutations of the basis elements r 1 , r 2 , and 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 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 ⊆ U R .

B. 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 Eq. 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.

Groups in Smith Normal Form
Begin with the simplest case.Let N be a non-singular 3 × 3 matrix of integers.Its columns represent the basis for a subgroup L N of the group Z 3 (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 11 whose basis is given by the columns of N. Since Z 3 and its subgroups are Abelian, we know that all the subgroups are normal so there exists a quotient group G = Z 3 /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 L N within Z 3 .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 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, Any finite Abelian group can be uniquely written in this form.(Isomorphic groups will yield the same "invariant factors" m 1 , m 2 ,. . .,m k when written in this form.) Note that, since G = Z 3 /L N , there must be a homomorphism from Z 3 onto G, having L N as its kernel.In other words, L N = {p ∈ Z 3 : ψ(p) = 0}.Our task is to find the direct-sum representation of the quotient group Z 3 /L N , and also to find the homomorphism ψ which maps the points of Z 3 onto the group (in such a way that ψ(p) = 0 iff p ∈ L N ).This allows us to work with the elements of the group as proxies for the k-points inside the reciprocal cell.

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 N effectively tells one how to represent the quotient group Z 3 /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 N give the homomorphism ψ suggested above.To put it another way, define * to be the operation Then we have shown w ∈ L N iff (Aw) * = (0, 0, 0) (the zero-element in the group That suggests we let ψ( w) = (A w) * , a homomorphism from Z 3 onto the direct-sum G 0 .Then, since that homomorphism is easily shown to be onto, and its kernel is L N , we see (by the First Isomorphism Theorem of group theory) that G 0 ∼ = Z 3 /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 D, and the transition matrix A provides the homomorphism which maps the parent lattice onto the group.
Thus we now know that the quotient group is Further, from the matrix A, we may obtain the homomorphism projecting Z 3 onto the quotient group, with Note that this homomorphism provides a different, but convenient, way to describe the superlattice.Since L N is the kernel of ψ, it is comprised of the points (x, y, z) ∈ Z 3 which satisfy the simultaneous congruences z ≡ 0 (mod 2) and x + 5y + 4z ≡ 0 (mod 6).We note that all three basis points p 1 , p 2 and p 3 satisfy these congruences, and thus so will all their integer linear combinations (all points in L N ).z (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.

Non-integer lattices
Now, what about the more complicated situation, where N represents a (possibly HNF) matrix describing the change from some lattice other than the simple integer lattice Z 3 to one of its subgroups (superlattice)?
Then we have a basis V and lattice L V and a basis W = VN for a (super) lattice L W . Again, the quotient group G = L V /L W is Abelian of order |det(N)|.Again, G is a direct sum of cyclic groups corresponding to the diagonal entries of D = ANB (where D is the SNF of N).
The only difference here is that the homomorphism ψ provided by A must depend on the basis V (which might even be irrational).Every point in L V has the form x = V w where w is a column of integers.Then ψ( x) = A w (modded by the corresponding entries from D = ANB).We could write it as ψ( x) = AV −1 x * (with the entries appropriately modularly reduced and transposed to a horizontal vector).
FIG.1.Total energy error vs. k-point density for the cases of silicon and aluminum.Silicon does not have a Fermi surface so there is no discontinuity in the occupied bands; convergence is super-exponential.(See the discussion of example 1 in Ref. 3.) In contrast, the total energy of aluminum converges very slowly, and the convergence is quite erratic.For typical target accuracies in the total energy, around 10 −3 eV/atom, metals require 10-50 times more k-points than semiconductors.

FIG. 6 .
FIG.6.Each point along the boundary of UR has at least one interior cousin closer to the origin when R is Minkowski reduced.For the points along the edge in red, these interior cousins are the points along the dashed red line.

FIG. 7 .
FIG. 7.Each point along the boundary of UR, the edges of which are shown in black, has at least one interior cousin closer to the origin when R is Minkowski reduced.For the points on the bounding plane in red, the interior cousins are the points on the plane in blue.(The origin is contained in the blue plane.)

3.
The connection between SNF groups and SNF matricesIn the matrix case, since the operations are elementary row and column operations, we have D = ANB where A and B are integer matrices with determinant ±1 representing the accumulated row operations and column operations respectively.The matrix D is completely determined by N, but the matrices A and B depend on the algorithm used to arrive at the Smith Normal Form of N. A different implementation might yield D = A NB (same N and same D, but different A and B).Note that, since B represents elementary column operations, the product NB simply represents a change of basis from N to a new basis N = NB.In other words, the columns of N are still a basis for L N .But the new basis has the property that AN = D.That means that every element w = N z of L N (where z is some element of Z 3 ) will satisfy the equation A w = D z = w will be a vector whose entries are multiples of the corresponding diagonal entries in D.
all the points which are integer linear combinations of those three points.The matrix N has determinant 12, which must be the volume of each lattice tile-and it is also the order of the quotient group Z 3 /L N .Using the SNF algorithm to diagonalize this basis matrix, we find D = ANB where D = mod 2), x + 5y + 4z (mod 6) (noting that anything mod 1 is zero).
Algorithmic variation.In the example we computed above, a different application of the SNF matrix algorithm, with the same N, might have yielded the same diagonal matrix D = would change the homomorphism to (x, y, z) → −5x + 3y + z (mod 2), −2x + 2y + z (mod 6) = x + y + L W is the subgroup lattice defined by the basis marix W = VN where N = other words, one basis for L W is given by the columns of W = Thus our quotient group isG = L V /L W ∼ = Z 2 ⊕Z 2 ⊕Z 4 ) into its Smith Normal Form, D = ANB.By elementary row and column operations, represented by unimodular matrices A and B, it is possible to transform N into a diagonal matrix D, where each diagonal element divides the ones below it: d 11 |d 22 |d 33 , and d 11 • d 22 • d 33 = n = |N| (the notation i|j means that i is divisible by j).As explained in the appendix (Sec.VII), when 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 K, then the coordinates (coefficients) of the interior points will obey Eq. 4. When these conditions are met, the hashing algorithm discussed above (in particular, Eq. 3) becomes possible.This enables the O(N k ) algorithm.