# Nonuniform grids for Brillouin zone integration and interpolation

Siyu Chen<EMAIL_ADDRESS>TCM Group, Cavendish Laboratory, University of
Cambridge,J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom  Pascal T.
Salzbrenner  Department of Materials Science and Metallurgy, University of
Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom
Bartomeu Monserrat<EMAIL_ADDRESS>TCM Group, Cavendish Laboratory,
University of Cambridge,J. J. Thomson Avenue, Cambridge CB3 0HE, United
Kingdom  Department of Materials Science and Metallurgy, University of
Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom

###### Abstract

We present two developments for the numerical integration of a function over
the Brillouin zone. First, we introduce a nonuniform grid, which we refer to
as the Farey grid, that generalizes regular grids. Second, we introduce
symmetry-adapted Voronoi tessellation, a general technique to assign weights
to the points in an arbitrary grid. Combining these two developments, we
propose a strategy to perform Brillouin zone integration and interpolation
that provides a significant computational advantage compared to the usual
approach based on regular uniform grids. We demonstrate our methodology in the
context of first principles calculations with the study of Kohn anomalies in
the phonon dispersions of graphene and MgB2, and in the evaluation of the
electron-phonon driven renormalization of the band gaps of diamond and
bismuthene. In the phonon calculations, we find speedups by a factor of $3$ to
$4$ when using density functional perturbation theory, and by a factor of $6$
to $7$ when using finite differences in conjunction with supercells. As a
result, the computational expense between density functional perturbation
theory and finite differences becomes comparable. For electron-phonon coupling
calculations we find even larger speedups. Finally, we also demonstrate that
the Farey grid can be expressed as a combination of the widely used regular
grids, which should facilitate the adoption of this methodology.

††preprint: APS/123-QED

##  I Introduction

Many macroscopic properties of crystalline materials are approximated by
averaging the contributions from every unit cell in an infinite crystal.
Invoking Bloch’s theorem, this real-space average can be equivalently
evaluated as an integral of a periodic function of the wavevector $\bm{k}$ in
the Brillouin zone of reciprocal space. In practice, the integral is
approximated with a sum over a discrete set of points $\\{\bm{k}_{i}\\}$ and
corresponding weights $\\{w_{i}\\}$, and several schemes have been proposed to
generate special points to perform these integrals at small computational cost
[1, 2, 3, 4, 5, 6, 7, 8], with the most popular being that of Monkhorst and
Pack [4]. These grids are a cornerstone of first principles calculations, and
they have been used to evaluate Brillouin zone averages of properties related
to (quasi-)particles including electrons [9], phonons [10], excitons [11], and
magnons [12], as well as their interactions [13, 14, 15].

One common feature of special points grids is that the $\bm{k}$-points are
uniformly distributed over the Brillouin zone and the corresponding weights
$w$ can be easily obtained using the point group symmetries of the system.
These regular grids have been tremendously successful in many applications,
and a key example is their central role in density functional theory where
they are used to evaluate the self-consistent potential. However, some
quantities exhibit sharp variations over the Brillouin zone, and the required
grid sizes can become computationally prohibitive. Examples include Kohn
anomalies in phonon dispersions [16, 17] and electron-phonon coupling matrix
elements [18, 19, 20], which often necessitate sophisticated interpolation
methods [21, 14, 22, 23] to achieve the required grid sizes.

In this paper, we propose two independent developments that generalize regular
$\bm{k}$-point grids and weights. First, we introduce a nonuniform grid, the
Farey grid, which contains regular grids as a subset of its $\bm{k}$-points.
Second, we introduce symmetry-adapted Voronoi tessellation as a general
strategy to calculate the weights $w$ of an arbitrary grid in the Brillouin
zone and subject to translational and point group symmetries.

Farey grids provide two advantages compared to regular grids. The first
advantage is that the grids used in convergence tests, which are discarded
when using regular grids, are all incorporated into the Farey grid. This
typically implies that convergence is reached using fewer $\bm{k}$-points,
accelerating calculations. The second advantage is restricted to calculations
in which the Brillouin zone is sampled exploiting commensurate real space
supercells, for example when calculating phonons using finite differences. In
this case, converged Farey grids typically include $\bm{k}$-points that
correspond to smaller commensurate real space supercells, again accelerating
calculations.

Symmetry-adapted Voronoi tessellation provides a general strategy to assign
weights to any $\bm{k}$-point grid. The weights reduce to the standard weights
for uniform grids, but are also applicable to nonuniform grids. Although we
use symmetry-adapted Voronoi tessellation in conjunction with Farey grids in
this paper, other applications could include stochastic grids or physically-
inspired tailored irregular grids.

We find that combining Farey grids with symmetry-adapted Voronoi tesselation
dramatically reduces the computational cost of first principles calculations,
and a computational cost reduction is observed irrespective of whether the
Brillouin zone grids are accessed directly or using commensurate real space
supercells. For example, using density functional perturbation theory, which
directly samples the Brillouin zone grid, the convergence of the Kohn
anomalies in the phonon dispersions of graphene and MgB2 is accelerated by
factors of three and four, respectively. The same calculations using finite
differences with real space supercells are accelerated by factors of six and
seven, respectively. The larger acceleration observed in the latter implies
that overall computational costs become similar between the two methods, and
in the case of graphene, the finite displacement method using supercells is
computationally cheaper than density functional perturbation theory. Similar
computational cost gains are observed in the evaluation of electron-phonon
interactions in diamond and bismuthene.

We also highlight that, while the results reported for graphene, MgB2, and
diamond have been reported elsewhere and are here used as a proof-of-
principle, bismuthene has not been studied before in this context. One
important conclusion revealed by our results is that bismuthene, a known
topological insulator at low temperature, will remain a topological insulator
up to room temperature.

The paper is organized as follows. We first introduce the nonuniform Farey
grid in Sec. II and compare it with the widely used uniform regular grids. In
Sec. III, we describe the assignment of weights to points on an arbitrary grid
by means of symmetry-adapted Voronoi tessellation. Section IV extends the
ideas of Farey grids and symmetry-adapted Voronoi tesselation from BZ
integration to BZ interpolation. The proposed methods are illustrated in Sec.
V with two examples: the calculation of Kohn anomalies in graphene and MgB2,
and the calculation of electron-phonon driven band gap renormalization in
diamond and bismuthene. Finally, a short summary and conclusions are presented
in Sec. VI.

##  II Brillouin zone grids

Let $f(\bm{k})$ be a scalar function defined in $\tau$-dimensional reciprocal
space and with the same symmetry as the reciprocal lattice of the crystal. A
general integral of this function over the Brillouin zone (BZ) reads

$I=\frac{1}{\Omega}\int_{\mathrm{BZ}}f(\bm{k})\;\mathrm{d}^{\tau}\bm{k},$ (1)

where $\Omega$ is the $\tau$-dimensional volume of the BZ. According to the
definition of the Riemann integral, Equation (1) can always be approximated by
a weighted sum over a discrete set of points in the BZ:

$I\approx\frac{1}{W}\sum_{\bm{k}_{i}\in\mathrm{BZ}}f\left(\bm{k}_{i}\right)w_{i},$
(2)

where $W=\sum_{i}w_{i}$, and $w_{i}$ is the weight for the $i$-th point. We
henceforth refer to the discrete set $\\{\bm{k}_{i}\\}$ as a grid.

Figure 1: (color online). Workflow of the convergence test to obtain the
converged sampling density using the regular grid (left, red panel) and the
Farey grid (right, green panel). Figure 2: (color online). (a-c) 2D Farey
grids of order $l=\text{(a) }8,\text{(b) }9\text{ and }\text{(c) }10$, where
the red dots mark the sampling points of the $l\times l$ regular grid, and the
black dots mark the “easy” sampling points missed by it. (d) The graph of
$J(\bm{k})$ in 2D for $s=1$, where green represents small and red large
values.

###  II.1 Regular Brillouin zone grids

One possible implementation of Eq. (2) is to find an optimal set
$\\{(\bm{k}_{i},w_{i})\\}$ such that as many terms as possible in the Fourier
expansion of $f(\bm{k})$ are integrated exactly, and this is the strategy
followed by the various special points methods [1, 2, 3, 4, 5, 6, 7, 8]. It is
worth noting that the logic behind these methods is essentially an analogy of
Gauss-Legendre quadrature, where Legendre polynomials are replaced by
trigonometric polynomials in order to exploit the translational symmetry of
$f(\bm{k})$.

A regular grid is indeed the optimal grid in this case in terms of maximizing
the number of Fourier components that are exactly integrated given a fixed
sampling density. This can be rationalized by considering the one-dimensional
(1D) case in which the integration evaluated from $n$ equidistant points
yields an exact result for any function representable by a Fourier series of
$n$ Fourier components or fewer:

$f(k)=\sum_{i=0}^{n-1}\left[C_{i}\cos\left(\frac{2\pi
ik}{G}\right)+S_{i}\sin\left(\frac{2\pi ik}{G}\right)\right],$ (3)

where $C_{i}$ and $S_{i}$ are the Fourier coefficients for the $i$-th Fourier
components and $G$ is the period of $f(k)$. The conclusion holds for higher
dimensions and also holds when the Fourier series is symmetrized according to
the point group of the crystal.

In summary, the special points methods (which could be called Fourier
quadrature methods) invariably lead to regular grids in which the sampling
points are uniformly distributed.

###  II.2 Computational cost of regular Brillouin zone grids

In a practical BZ integration using a regular grid, the converged sampling
density is not known in advance. The general convergence workflow is depicted
in the left panel of Fig. 1: sequential calculations are performed with
increasing linear grid size $l$ until the integral agrees within some pre-
defined tolerance between successive grid sizes. In this process, all
intermediate calculations necessary to reach the converged grid size are
discarded.

Figures 2(a-c) depict the regular grids considered in a typical convergence
workflow. The red dots in Fig. 2(a) correspond to a $8\times 8$ grid, and the
black dots correspond to all the discarded grids of sizes between $1\times 1$
and $7\times 7$. Figure 2(b) depicts the subsequent step in the convergence
workflow, with the red dots corresponding to a $9\times 9$ grid, and the black
dots to all the discarded grids of sizes between $1\times 1$ and $8\times 8$.
Figure 2(c) corresponds to the next step, when the sampled grid is of size
$10\times 10$.

In calculations where the BZ grids are accessed directly, the computational
cost of every $\bm{k}$-point is equivalent. By contrast, in calculations where
the BZ grids are accessed using commensurate real space supercells, the
smallest supercell commensurate with $\tau$-dimensional wavevector
$\bm{k}=\left({u_{1}}/{v_{1}},\cdots,{u_{\tau}}/{v_{\tau}}\right)$, written in
reduced fractional coordinates, has a size of
$S_{\text{min}}=\mathrm{lcm}\left(v_{1},\cdots,v_{\tau}\right)$ [24]. Based on
these insights, we define a reciprocal computational cost function on the
$\tau$-dimensional BZ as follows:

$J(\bm{k})=\begin{cases}\frac{1}{|S_{\text{min}}|^{s}}&\text{ if all entries
of }\bm{k}\text{ are rational }\\\ \;0^{s}&\text{ if }\text{any entry of
}\bm{k}\text{ is irrational}\end{cases}.$ (4)

In other words, $J(\bm{k})$ is defined such that the computational cost needed
to evaluate the $\bm{k}$-point is proportional to $J(\bm{k})^{-1}$. The
parameter $s$ characterizes the scaling of the solver used. For example, in
calculations where the BZ grids are accessed directly, $s=0$ because each
$\bm{k}$-point has an equivalent cost. As another example, in calculations
where the BZ grids are accessed using commensurate real space supercells, and
in which density functional theory (DFT) is used, then $s=3$ because the
computational cost of DFT scales cubically with system size. The 2D version of
$J(\bm{k})$ with $s=1$ is visualized in Fig. 2(d), and it exhibits a fractal
structure dominated by the lines with rational slopes. Conceptually, Figure
2(d) shows how computationally expensive it is to sample a particular point in
the BZ.

Figure 2 conveys two important insights. The first insight is that,
irrespective of the method used to sample the BZ, a typical convergence
workflow (Fig. 1, left panel) requires many calculations that are ultimately
discarded, as shown by the black dots in Figs. 2(a-c). The second insight is
that, when the BZ grid is sampled using commensurate real space supercells,
the computational cost of accessing a given $\bm{k}$-point varies
significantly across the BZ, as shown in Fig. 2(d). Both insights suggest that
the use of regular grids with standard convergence workflows is
computationally wasteful.

###  II.3 Farey Brillouin zone grids

Unlike the widely used special points methods and associated regular grids [1,
2, 3, 4, 5, 6, 7, 8], whose objective is to maximize the number of Fourier
components of $f(\bm{k})$ that are exactly integrated given a fixed sampling
density, our starting point is to maximize the sampling density one can
achieve with a fixed computational cost.

####  II.3.1 Farey Brillouin zone grids in 1D

We start with the 1D case to illustrate the key ideas. In 1D, the reciprocal
computational cost function defined in Eq. (4) becomes:

$J(k)=\begin{dcases}\frac{1}{|q|^{s}}&\text{ if
}|k|=\frac{|p|}{|q|},\;p,q\in\mathbb{Z},\;\gcd(p,q)=1\\\ \;0^{s}&\text{ if
}k\text{ is irrational }\end{dcases}.$ (5)

Note that $-{1}/{2}<k\leq{1}/{2}$ is represented in units of the reciprocal
lattice and $|q|$ corresponds to the size of the smallest supercell
commensurate with $k$. Also remember that $s$ characterizes the scaling of the
solver used, and that overall $J(k)$ is defined such that the computational
cost needed to evaluate the $k$-point is proportional to $J(k)^{-1}$.

It is worth noting that when $s=1$, $J(k)$ reduces to the Thomae’s function, a
famous pathological function with an intriguing fractal structure. Thomae’s
function is nowhere differentiable, but it has been shown that its local
maxima are located strictly at every rational number. Therefore, one can
easily deduce that the positions of the top $l\in\mathbb{N}^{*}$ maxima of
$J(k)$ for any $s>0$ are at:

$\displaystyle\mathscr{F}\left(l\right)=$ (6)
$\displaystyle\left\\{\frac{p}{q}\;\middle|\;q\in\mathbb{N^{*}}\leq
l,\;-\frac{q}{2}<p\in\mathbb{Z}\leq\frac{q}{2},\;\gcd(p,q)=1\right\\}.$

It is clear that the most cost-effective scheme to sample the BZ is to follow
the points provided by $\mathscr{F}\left(l\right)$. We henceforth refer to
this kind of grid generated according to $\mathscr{F}\left(l\right)$ as the 1D
Farey grid of order $l$, as the grid is isomorphic to the Farey sequence of
order $l$ modulo $1$ [25] when its sampling points are arranged in ascending
order.

The 1D Farey grid has two noteworthy properties inherited from the Farey
sequence [26]:

  * •

For any subinterval $I=\left(\alpha,\beta\right]\subseteq(-1/2,1/2]$,
$\left|I\cap\mathscr{F}\left(l\right)\right|/\left|\mathscr{F}\left(l\right)\right|=\left(\beta-\alpha\right)+\mathscr{O}\left(l^{-1}\log
l\right)$.

  * •

$\mathscr{F}\left(l\right)\supsetneq\mathscr{F}\left(l-1\right)$ and
$\left|\mathscr{F}\left(l\right)-\mathscr{F}\left(l-1\right)\right|=\phi(l)$,
where $\phi(l)$ is Euler’s totient function.

The first property states that, as $l$ increases, the 1D Farey grid tends to
become uniform at the same speed as $\mathscr{O}\left(l^{-1}\log l\right)$
approaches zero. This asymptotically uniform nature guarantees that the 1D
Farey grid is an almost unbiased sampling grid. It follows that for a
sufficiently large sampling density, the error for the BZ integration
calculated with the Farey grid ought to be close to that estimated by Chadi
and Cohen for the regular grid [2]. For the same reason, interpolation between
different points on a sufficiently dense Farey grid should also perform as
well as that on a regular grid with the same sampling density. Overall, just
like for regular grids, by deploying successively finer Farey grids one can
achieve increasingly precise BZ integration and interpolation.

The second property states the 1D Farey grid has a compact nested hierarchy. A
Farey grid of order $l$ contains all sampling points of any lower-order grid,
all the way to the Farey grid of order $1$ (that is, the $\Gamma$-only grid).
More specifically, increasing the order $l$ of the Farey grid by $1$, then the
new Farey 1D grid or order $l+1$ will have $\phi(l+1)$ more sampling points
than the order $l$ grid. The average growth rate of Euler’s totient function
is known to be approximately constant, namely ${6}/{\pi^{2}}$. From this
linear growth, one can estimate that a 1D Farey grid of order $l$ samples
approximately $\sum_{i=1}^{l}{6i}/{\pi^{2}}=\mathscr{O}\left(l^{2}\right)$
points in the BZ and the corresponding computational cost is
$\sum_{i=1}^{l}i^{s}=\mathscr{O}\left(l^{s+1}\right)$. By comparison, a
regular grid would require a computational effort of
$\mathscr{O}\left(l^{2s}\right)$ to achieve the same sampling density.

As a result of the second property, the Farey grid resolves the computational
waste associated with the discarded points in a typical convergence workflow
using regular grids (Fig. 1, left panel). The true computational cost of a
calculation that converges with a regular grid of size $l$ requires the
addition of an extra computational cost of $\mathscr{O}\left(l^{2s+1}\right)$.
In contrast, due to the nested hierarchical structure of the Farey grid, the
corresponding convergence workflow does not add any computational cost (Fig.
1, right panel), as the sampled results from the lower-order Farey grids are
recycled. Therefore, even considering a practical convergence workflow, the
total computational cost of the Farey grid is still
$\mathscr{O}\left(l^{s+1}\right)$.

In this sense, sampling with the Farey grid always outperforms sampling with
the regular grid. Consider an example 1D system whose converged results
require semi-uniformly sampling $128$ points. In this case, one needs to carry
out $128$ loops of the regular grid convergence test in order to determine
this result, while only $20$ loops of the Farey grid convergence test will be
necessary (as $\sum_{j=1}^{20}\phi(j)=128$). Importantly, this conclusion
holds for any method used to sample the BZ. For methods in which the grid
points are accessed directly [$s=0$ in the context of Eq. (5)], the
computational gains arise due to overall sampling fewer $k$-points. For
methods in which the grid points are accessed using commensurate real space
supercells [$s>0$ in the context of Eq. (5)], the advantage of the Farey grid
becomes more pronounced as convergence is achieved with smaller supercells.
The advantage of the Farey grid compounds as the scaling $s$ of the solver
becomes larger.

####  II.3.2 Farey Brillouin zone grids in arbitrary dimensions

We define the $\tau$-dimensional Farey grid of linear order $l$ as follows:

$\displaystyle\mathcal{K}^{(\tau)}_{\text{far}}\left(l\right)=$ (7)
$\displaystyle\bigg{\\{}\left(\frac{u_{1}}{v_{1}},\cdots,\frac{u_{\tau}}{v_{\tau}}\right)\;\bigg{|}\;\frac{u_{i}}{v_{i}}\in\mathscr{F}\left(l\right),\mathrm{lcm}\left(v_{1},\cdots,v_{\tau}\right)\leq
l\bigg{\\}}.$

In this notation,
$\mathcal{K}^{(1)}_{\text{far}}\left(l\right)=\mathscr{F}\left(l\right)$.
Also, in the context of supercell calculations, the condition
$\mathrm{lcm}\left(v_{1},\cdots,v_{n}\right)\leq l$ filters out all sampling
points whose corresponding smallest commensurate supercells have sizes greater
than $l$ such that the sampling density is, as desired, strictly maximized at
a fixed computational cost. It is worth noting that this condition preserves
the asymptotic uniformity and nested hierarchy which characterized the 1D
Farey grids.

Going back to Figs. 2(a-c), they clearly demonstrate how a Farey grid resolves
the wasteful calculations associated with regular grids. For a 2D Farey grid
of order $8$, all points (red and black dots) in Fig. 2(a) are included in the
evaluation of the BZ integral. By contrast, a regular grid requires
calculations on all points, but only the points indicated with red dots are
included in the final BZ integral. Figure 2(b) depicts the corresponding
situation for a Farey grid of linear order $9$, and Fig. 2(c) for a Farey grid
of linear order $10$.

The main difference between the higher-dimensional Farey grid compared to 1D
is in the growth rate of the grid size. Figure 3 depicts
$|\mathcal{K}^{(\tau)}_{\text{far}}\left(l\right)-\mathcal{K}^{(\tau)}_{\text{far}}\left(l-1\right)|^{1/\tau}$,
which characterizes the number of new sampling points introduced upon
increasing the Farey grid linear order, in $\tau=1,2\text{ and }3$ dimensions.
It is clear that irrespective of the dimensionality,
$|\mathcal{K}^{(\tau)}_{\text{far}}\left(l\right)-\mathcal{K}^{(\tau)}_{\text{far}}\left(l-1\right)|^{1/\tau}$,
on average, grows linearly. The local fluctuations around the curves are
reduced in higher dimensions, with the 3D case approaching a smooth curve with
an average slope close to $1$. Overall, a $\tau$-dimensional Farey grid of
order $l$ has $\mathscr{O}\left(l^{\tau}\right)$ more sampling points than one
of order $l-1$, and it contains $\mathscr{O}\left(l^{\tau+1}\right)$ points in
total.

Figure 3: (color online).
$|\mathcal{K}^{(\tau)}_{\text{far}}\left(l\right)-\mathcal{K}^{(\tau)}_{\text{far}}\left(l-1\right)|^{1/\tau}$,
which characterizes how many new sampling points are introduced when one
increments the $\tau$-dimensional Farey grid, as a function of $l$.

####  II.3.3 Relation between Farey and regular grids

A Farey grid of order $l$ can be treated as a union of regular grids with
their linear size increasing step-by-step from $\lfloor l/2\rfloor+1$ to $l$:

$\mathcal{K}^{(\tau)}_{\text{far}}\left(l\right)=\bigcup_{j=\lfloor
l/2\rfloor+1}^{l}\mathcal{K}^{(\tau)}_{\text{reg}}\left(j\right).$ (8)

In this expression,

$\displaystyle\mathcal{K}^{(\tau)}_{\text{reg}}\left(l\right)=\left\\{\left(\frac{u_{1}}{v_{1}},\cdots,\frac{u_{3}}{v_{3}}\right)\;\middle|\;\frac{u_{i}}{v_{i}}\in\mathscr{A}\left(l\right)\right\\}$
(9)

denotes a $\tau$-dimensional $\Gamma$-centered regular grid generated by a
finite set of the arithmetic sequence

$\displaystyle\mathscr{A}\left(l\right)=$
$\displaystyle\left\\{\frac{j}{l}\;\middle|\;-\frac{l}{2}<j\in\mathbb{Z}\leq\frac{l}{2}\right\\}$
(10)

with a common difference of $1/l$.

It follows from Eq. (8) that the Farey grid can be constructed by combining
the results (insofar as they exist) from consecutive regular grids, and this
implies that software implementing the regular grid, including all major DFT
codes, can be easily modified to adopt the Farey grid. Indeed, this is the
approach taken in the calculations reported in Sec. V of this work.

We note that, in practice, it may be more natural to start with a tentative
regular grid of linear size $a\in\mathbb{N^{*}}$ and combine it with
increasingly large regular grids, step by step, until a linear size
$b\in\mathbb{N^{*}}$ :

$\mathcal{K}^{(\tau)}_{\text{far}}\left(a,b\right)=\bigcup_{j=a}^{b}\mathcal{K}^{(\tau)}_{\text{reg}}\left(j\right).$
(11)

If $b<2(a-1)$, this corresponds to a “truncated” Farey grid. This kind of
Farey grid misses a small number of “easy” points whose fractional coordinates
have a component of $p/q\in\mathscr{H}_{\text{miss}}\left(a,b\right)$, where

$\displaystyle\mathscr{H}_{\text{miss}}\left(a,b\right)=$ (12)
$\displaystyle\left\\{\frac{p}{q}\in\mathscr{F}\left(b,0\right)\;\middle|\;\forall\;a\leq
n\in\mathbb{N^{*}}\leq b\;\gcd(np,q)\neq q\right\\},$

but it can still be very useful in practical calculations.

##  III Symmetry-Adapted Voronoi Tessellations

Figure 4: (color online). SAVTs of the oblique lattice IBZ for the regular
grids (a) $\mathcal{K}^{(2)}_{\text{reg}}\left(3\right)$, (b)
$\mathcal{K}^{(2)}_{\text{reg}}\left(4\right)$, (c)
$\mathcal{K}^{(2)}_{\text{reg}}\left(5\right)$ and (d) the Farey grid
$\mathcal{K}^{(2)}_{\text{far}}\left(5\right)=\mathcal{K}^{(2)}_{\text{reg}}\left(3\right)\cup\mathcal{K}^{(2)}_{\text{reg}}\left(4\right)\cup\mathcal{K}^{(2)}_{\text{reg}}\left(5\right)$.
The grid points in the IBZ of interest (colored) are marked by ‘x’. The color
distinguishes the Voronoi cells generated by the different regular grids.

In this section we consider the assignment of weights to an arbitrary
Brillouin zone grid.

###  III.1 Weights for a regular grid

All the points of a regular grid carry exactly the same weight. This is
reasonable, since by definition the grid points are evenly spaced in each
independent direction, so there is no reason to think of any sampling point as
more important than any other.

Regular grids allow for an efficient symmetry reduction since a regular grid
spanned by the reciprocal basis vectors of the BZ is essentially a scaled
version of the reciprocal lattice. This means that all the points of a regular
grid can be mapped to the same irreducible BZ (IBZ) by some point group
symmetry operation of the crystal. As a result, the BZ integration can be
restricted to the IBZ,

$I\approx\frac{1}{W}\sum_{\bm{k}_{i}\in\mathrm{IBZ}}f\left(\bm{k}_{i}\right)w_{i},$
(13)

by recasting the weights as the symmetry multiplicity of the sampling points
in the IBZ. In other words, the weight of a grid point is equal to the
cardinality of its orbit in the crystal point group:

$w_{i}=\left|\mathrm{Orb}\left(\bm{k}_{i}\right)\right|.$ (14)

###  III.2 Weights for a general grid

One may expect that a general grid (for example, the nonuniform Farey grid
introduced in Sec. II.3) can also be reduced by symmetry, as in the case of
regular grids. To address this question, we introduce symmetry-adapted Voronoi
tessellation (SAVT).

Consider first plain Voronoi tessellation (VT) in Euclidean space. Given a
$\tau$-dimensional space $\mathbb{R}^{\tau}$ and a set of distinct points
$\mathcal{K}=\left\\{\bm{k}_{i}\right\\}$, for each point $\bm{k}_{i}$, the
corresponding Voronoi cell is given by

$C_{i}(\mathcal{K})=\left\\{\bm{k}\in\mathbb{R}^{\tau}\mid\;\left\|\bm{k}-\bm{k}_{i}\right\|<\left\|\bm{k}-\bm{k}_{j}\right\|\;\forall\;i\neq
j\right\\},$ (15)

where $\|\dots\|$ denotes the Euclidean distance. That is, the Voronoi cell of
a given point includes all space closer to it than to any other point in the
set.

To incorporate the crystal point group symmetry $\mathcal{G}$ and the
periodicity of the reciprocal lattice $\mathbf{G}$, we change the space metric
to extend Eq. (15) to define the SAVT of a bounded domain $\Pi=\text{IBZ}$ as
follows:

$C_{i}(\mathcal{K})=\left\\{\bm{k}\in\Pi\mid\;{}_{\mathcal{G}}\left\|\bm{k}-\bm{k}_{i}\right\|_{\mathbf{G}}<{}_{\mathcal{G}}\left\|\bm{k}-\bm{k}_{j}\right\|_{\mathbf{G}}\;\forall\;i\neq
j\right\\},$ (16)

where

${}_{\mathcal{G}}\|\bm{k}-\bm{k}_{i}\|_{\mathbf{G}}=\min_{g\in\mathcal{G},\;\eta\in\mathbf{G}}\left(\|\bm{k}-g\bm{k}_{i}-\eta\|\right).$
(17)

One can easily show that for $i\neq j$, $C_{i}\cap C_{j}=\emptyset$, so the
set $\left\\{C_{i}\right\\}$ defined above is indeed a tessellation of $\Pi$.
It is worth further noting that:

  * •

There are two types of Voronoi cells in the SAVT, inner cells which are
entirely inside the domain and marginal cells which touch the boundary of the
domain. The inner cell is the same as the plain Voronoi cell, while the
marginal cell can either be composed of just one singly connected polyhedron
or several visually disconnected polyhedra.

  * •

The metric function defined in Eq. (17) guarantees that the SAVT of the IBZ
conserves both the point group and translational symmetries of the crystal.
Therefore, the corresponding Voronoi diagram allows for a symmetry reduction
similar to that used in regular grids. The weight of a sampling point is just
equal to the volume of its corresponding symmetry-adapted Voronoi cell:

$w_{i}=\left|C_{i}\right|.$ (18)

This is one of the main conclusions of this work.

Equation (18) can be rationalized by considering the 2D example of the SAVTs
of the oblique lattice IBZ shown in Fig. 4. From (a) to (d), the SATVs are
constructed for the regular grids
$\mathcal{K}^{(2)}_{\text{reg}}\left(3\right)$,
$\mathcal{K}^{(2)}_{\text{reg}}\left(4\right)$ and
$\mathcal{K}^{(2)}_{\text{reg}}\left(5\right)$ as well as the Farey grid
$\mathcal{K}^{(2)}_{\text{far}}\left(5\right)$, which is equivalent to their
union. For the regular grids, the SAVT of the IBZ generates a number of
Voronoi cells with exactly the shape of the BZ, except for some marginal
Voronoi cells that are cut off by the IBZ boundary. It is important to realize
that some marginal cells appear not to be associated with any grid point, but
they are in fact connected to some other marginal cell that visually possesses
a point. These “point-less” marginal cells, from another perspective, can be
treated as parts of a Voronoi cell associated with a symmetry image of the
point in an adjacent IBZ which intersects with the IBZ of interest. The volume
associated with such a point is the sum of all symmetry-connected volumes
within the IBZ of interest. This is the main difference between SAVTs and
plain VTs. In the special case when the IBZ boundary edges coincide with
mirror planes, SAVTs and VTs become equivalent because the IBZ boundary
coincides with the boundaries of all marginal cells. In 2D space, except for
the oblique lattice, the IBZ boundaries of the other four Bravais lattices are
mirror-like. In 3D space, four out of the fourteen Bravais lattices also
possess this kind of mirror-like IBZ boundaries. These are the orthorhombic,
tetragonal, rhombohedral, and hexagonal lattices.

It is also worth noting that there is an intimate connection between the
geometry of the Voronoi cells and the symmetry of the associated points. In
order to show this, we consider the little group $\mathcal{L}_{\bm{k}}$ of
$\mathcal{G}$ for a given $\bm{k}$-point and the action of the symmetry
operator $g\in\mathcal{L}_{\bm{k}}$ on the corresponding Voronoi cell. For an
inner cell whose associated point always has a trivial little group containing
the identity operator only, it will, by definition, be simply mapped to itself
through $g$. However, for a marginal cell whose associated point is located
exactly on the IBZ boundary, thus having a higher symmetry, it will be mapped
in its entirety to an adjacent IBZ by $g$ other than the identity operator. In
the usual scheme for the regular grid, the weight of the boundary point is
generally reduced compared to that of an inner point by
$1/|\mathcal{L}_{\bm{k}}|$. In our geometric picture, because the symmetry
operations $g\in\mathcal{L}_{\bm{k}}$ map points on the IBZ boundary to
themselves, images created by the operation of these symmetries can never
contribute to the volume of the marginal cell within the IBZ of interest. The
volume of marginal cells associated with points on the IBZ boundary is
therefore also reduced by $1/|\mathcal{L}_{\bm{k}}|$ compared to the volume of
an inner cell – that is, its weight, determined from the volume of its
marginal cell, is reduced to exactly the same extent as its weight determined
by symmetry. This is why the geometric and the symmetry pictures are
equivalent when one considers the regular grid. For the more general case,
only the geometric approach developed here works.

We emphasize that using the SAVTs to reduce the sampling points to the IBZ and
obtain the corresponding weights is not restricted to Farey grids, but is also
suitable for any other nonuniform grid. Examples of these might include
stochastic grids or manually customized grids which sample some specific
regions of the BZ heavily but the remaining regions in a light-weight manner.

Finally, we point out the relation between the SAVT method and the currently
used weighting methods for the special points regular grids [1, 2, 3, 4, 5, 6,
7, 8]. As pointed out in Sec. II.1, the logic behind these special points
grids is analogous to Gaussian quadrature. The various methods are based on
the choice of a small set $\\{(\bm{k}_{i},w_{i})\\}$ such that as many terms
as possible in the symmetrized Fourier expansion of $f(\bm{k})$ are integrated
exactly. For this reason, these methods work effectively when the Fourier
coefficients converge rapidly to zero, that is, when $f(\bm{k})$ is
sufficiently smooth. When $f(\bm{k})$ is a rapidly varying function, a better
approach would be to expand $f(\bm{k})$ on a symmetrized basis that is
localized in $\bm{k}$-space

$\delta_{i}(\bm{k}-\bm{k_{i}})=\begin{dcases}\frac{1}{\left|C_{i}\right|}&\text{
if }k\in C_{i}\\\ 0&\text{ if }k\notin C_{i}\end{dcases}.$ (19)

This is precisely what the SAVT method does.

To some extent, the SAVT method can be treated as a generalization of the
special point method combined with the tetrahedron method, which was proposed
to resolve the BZ integration of metals where $f(\bm{k})$ usually varies
sharply near the Fermi energy [27, 28, 29]. The method divides the BZ into
tetrahedra in such a way that the edges of the tetrahedra have minimal lengths
for the given geometry of the Bravais lattice (in order to make the tetrahedra
as compact as possible). In fact, this is very close to the idea of Delaunay
triangulations. The relation becomes even more clear when one notes that the
dual graph of the Delaunay triangulation of a discrete point set is exactly
the Voronoi tessellation of the same point set.

##  IV Natural neighbor interpolation in the Billouin zone

Rather than the integral of $f(\bm{k})$, some calculations require the value
of $f$ at some particular point $\bm{k}^{*}$. The construction of
$f(\bm{k}^{*})$ from the knowledge of $f$ on a $\bm{k}$-point grid is known as
BZ interpolation. In this section we describe a strategy to perform BZ
interpolation starting from a nonuniform Farey grid.

We can recast BZ interpolation as an instance of BZ integration using the
following integral in terms of the Dirac delta function:

$\displaystyle f(\bm{k}^{*})$
$\displaystyle=\int_{\mathrm{BZ}}f(\bm{k})\delta(\bm{k}-\bm{k}^{*})\;\mathrm{d}^{\tau}\bm{k}.$
(20)

BZ interpolation is the name given to evaluating this integral as a weighted
sum over a discrete $\bm{k}$-point grid, which results in the interpolation of
$f$ from the grid to $\bm{k}^{*}$. Due to the locality of the Dirac delta
function, the weighted sum for calculating the above integral only needs to
involve $\bm{k}$-points that are located close to the target point
$\bm{k}^{*}$. Based on this, Sibson proposed the so-called natural neighbor
interpolation scheme [30] to interpolate a function defined in
$\mathbb{R}^{\tau}$:

$\mathcal{N}[f(\bm{k}^{*})|f\left(\mathcal{K}\right)]=\frac{1}{W}\sum_{i=1}^{|\mathcal{K}|}f\left(\bm{k}_{i}\right)w_{i}^{*}.$
(21)

In this expression, $\mathcal{N}[f(\bm{k}^{*})|f\left(\mathcal{K}\right)]$ is
the natural neighbor interpolation of $f$ at $\bm{k}^{*}$ given the values of
$f$ at the grid points $\bm{k}_{i}\in\mathcal{K}$, and $W=\sum_{i}w_{i}^{*}$
is the normalization factor. According to Sibson’s definition, the weight
elements $w^{*}_{i}$ are equal to the volume of each of the surrounding
Voronoi cells which is “stolen” when inserting $\bm{k}^{*}$ into the original
Voronoi tessellation, that is

$w^{*}_{i}=\left|C_{i}(\mathcal{K}\cup\bm{k}^{*})\cap
C_{i}(\mathcal{K})\right|.$ (22)

To perform BZ interpolation starting from arbitrary grids, we extend this
interpolation strategy to SAVT. As we are interested in functions $f$ that
have the periodicity of the reciprocal lattice, we replace
$f\left(\bm{k}_{i}\right)$ in Eq. (21) by a series of Fourier interpolated
values at $\bm{k}^{*}$, that is

$\mathcal{N}[f(\bm{k}^{*})|f\left(\mathcal{K}\right)]=\frac{1}{W}\sum_{i=1}^{|\mathcal{K}|}\mathcal{F}[f(\bm{k}^{*})|f\left(\mathcal{K}_{i}\right)]w_{i}^{*},$
(23)

where $\mathcal{F}[f(\bm{k}^{*})|f\left(\mathcal{K}_{i}\right)]$ denotes the
Fourier interpolation of $f$ at $\bm{k}^{*}$ given the values of $f$ at the
grid points $\mathcal{K}_{i}$. $\mathcal{K}_{i}\subseteq\mathcal{K}$ is the
finest regular grid which contains the $i$-th point
$\bm{k}_{i}\in\mathcal{K}$, and $w^{*}_{i}$ is calculated using SAVT. It is
straightforward to see that when $\mathcal{K}$ is itself a regular grid,
Equation (23) reduces to Fourier interpolation, and when $\mathcal{K}$ is a
Farey grid, the decomposition of the grid into a series of regular grids is
always feasible according to Eq (11). As a result, the Fourier interpolation
required in Eq. (23) can always be performed. Henceforth we refer to the
interpolation based on Eq. (23) as natural neighbor Fourier interpolation.

In particular, when $\bm{k}^{*}$ coincides exactly with some grid point
$\bm{k}_{i}$, the natural neighbor Fourier interpolated value of $f$ is
guaranteed to be equal to the actual value, that is
$\mathcal{N}[f(\bm{k}^{*}=\bm{k}_{i})|f\left(\mathcal{K}\right)]=f\left(\bm{k}_{i}\right)$.
A limitation of this strategy is that
$\mathcal{N}[f(\bm{k}^{*})|f\left(\mathcal{K}\right)]$ is not differentiable.
In other words, the function obtained by natural neighbor Fourier
interpolation is smooth everywhere in the BZ except at the locations of the
grid points, where it is only piece-wise continuous. To obtain a smooth
interpolation of $f$ everywhere, we propose to first map the results evaluated
on the nonuniform grid $\mathcal{K}$ to a target regular grid
$\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}$ using natural neighbor Fourier
interpolation, and then to obtain an interpolation that is smooth everywhere
in the BZ through plain Fourier interpolation, as summarized below:

$\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}\xrightarrow{{k}^{\text{tar}}_{i}\in\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}}f(\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}})=\\{\mathcal{N}[f(\bm{k}^{\text{tar}}_{i})|f\left(\mathcal{K}\right)]\\}\xrightarrow{\bm{k}^{*}\in\text{BZ}}f\left(\bm{k}^{*}\right)=\mathcal{F}[f(\bm{k}^{*}|f(\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}})].$ (24)

The size of the target regular grid is, in principle, a free parameter.
However, we suggest the following two options for practical calculations when
the nonuniform grid is the Farey grid $\mathcal{K}=\mathcal{K}_{\text{far}}$:
(i) $|\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}|\approx|\mathcal{K}_{\text{far}}|$, that
is, the regular grid is chosen such that the sampling densities of the two
grids are almost equal; (ii)
$\mathcal{K}_{\text{far}}\subseteq\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}$, that is, choosing the smallest
$\mathcal{K}_{\begin{subarray}{c}\scalebox{0.6}{reg}\\\
\scalebox{0.6}{tar}\end{subarray}}$ which contains all grid points of
$\mathcal{K}_{\text{far}}$. As discussed in Sec. II.3, the Farey grid is an
asymptotically uniform grid, so option (i) is usually a sufficiently good
choice to access the interpolated function $f$. Indeed, we use option (i) for
the calculations reported in Sec. V below, in which we match the sampling
density along high symmetry lines of interest in the BZ. Option (ii)
guarantees that the interpolated curve strictly passes through all the data
points on the starting grid, but its evaluation can become computationally
costly when the size of the Farey grid is large.

##  V Farey grids and symmetry-adapted Voronoi tessellation in action

We now describe how to use Farey grids and SAVT to perform first principles
calculations in practice. We focus on two examples that have traditionally
proven challenging for calculations based on regular grids: (i) the
calculation of phonon dispersions with Kohn anomalies in metallic and
semimetallic systems, and (ii) the calculation of electron-phonon coupling
(EPC).

Unless otherwise indicated, first-principles calculations are carried out at
the DFT level [31, 32] with the Vienna ab initio Simulation Package (vasp)
[9]. The interaction between ions and valence electrons is modeled with
pseudopotentials based on the projector-augmented wave [33, 34] method. Some
of the calculations, clearly identified below, are performed using the castep
plane-wave DFT code [35] using norm-conserving pseudopotentials [36].

The phonon and EPC calculations are evaluated within the harmonic
approximation, and most are preformed using the finite differences (FD) method
[37] in conjunction with nondiagonal supercells [24]. Some phonon calculations
are also performed using linear response in the context of density functional
perturbation theory (DFPT) [38].

To follow standard notation, in this section we use $\bm{q}$ instead of
$\bm{k}$ to label points in the BZ related to phonons.

###  V.1 Kohn anomalies in the phonon dispersion

Phonons are a critical concept for understanding lattice vibration in solids.
Accurate phonon dispersions and densities of states are necessary to evaluate
properties of solids including thermal expansion, heat capacity, vibrational
entropy, and free energy, and also form the starting point for calculations of
quantities such as thermal conductivity and electronic conductivity. The
fundamental physical quantity for accessing phonons in solids is the so-called
dynamical matrix:

$D_{\alpha\alpha^{\prime}}(\bm{q})=\frac{1}{N\sqrt{m_{\alpha}m_{\alpha^{\prime}}}}\sum_{p\,p^{\prime}}C_{\alpha\alpha^{\prime}}\left(\bm{\mathbf{R}}_{p},\bm{\mathbf{R}}_{p^{\prime}}\right)\mathrm{e}^{i\bm{q}\cdot\left(\bm{\mathbf{R}}_{p}-\bm{\mathbf{R}}_{p^{\prime}}\right)},$
(25)

whose eigenvalues and eigenvectors correspond to the squared phonon
frequencies $\omega^{2}_{\nu\bm{q}}$ and phonon polarization vectors
$\bm{\xi}_{\nu\bm{q}}$, respectively. In Eq. (25), $N$ is the number of
primitive cells contained in a supercell subject to Born–von Karman periodic
boundary conditions, $m_{\alpha}$ denotes the mass of the $\alpha$-th atom in
the primitive cell, and
$C_{\alpha\alpha^{\prime}}\left(\bm{\mathbf{R}}_{p},\bm{\mathbf{R}}_{p^{\prime}}\right)$
is the force constant matrix, which is defined as the derivative of the force
on atom $\alpha$ in the primitive cell indicated by the Bravais lattice vector
$\bm{\mathbf{R}}_{p}$ with respect to displacement of atom $\alpha^{\prime}$
in the primitive cell indicated by $\bm{\mathbf{R}}_{p^{\prime}}$.

First principles DFT calculations of phonons can be performed using two
different methodologies. The first method is based on linear response theory
in the context of DFPT, which is a method that provides direct access to the
dynamical matrix on a grid of $\bm{q}$-points, often called the coarse
$\bm{q}$-point grid. The second method is based on finite differences and it
constructs the dynamical matrix on a grid of $\bm{q}$-points using
commensurate real space supercells. In both cases, the dynamical matrix
$D_{\alpha\alpha^{\prime}}(\bm{q})$ is exact at the $\bm{q}$-points on the
grid, and it can be approximated at other $\bm{q}$-points in the BZ using
Fourier interpolation. The size of the $\bm{q}$-point grid becomes a
convergence parameter as the dynamical matrix may deviate significantly from
its exact value when the $\bm{q}$-point grid is not dense enough.

The changes driven by the displacement of an atom from its equilibrium
position are usually screened by electrons, which implies that the matrix of
force constants decays with a rapid power law as a function of distance. As a
result, relatively small regular $\bm{q}$-point grids have been very
successful at converging phonon calculations. An exception to this rule-of-
thumb is provided by Kohn anomalies, which arise at nesting vectors of the
Fermi surface of metallic and semimetallic systems, and whose convergence
typically requires large regular $\bm{q}$-point grids that capture the
relevant nesting vectors. In the following, we evaluate the phonon dispersions
of graphene and MgB2 which host Kohn anomalies.

####  V.1.1 Computational details

For graphene, the exchange-correlation functional is treated in the local
density approximation (LDA) [39]. A vacuum layer with a thickness of $15$ Å is
used in the calculation to avoid periodic image interactions along the
direction perpendicular to the plane. Precise forces are obtained when the
energy is converged to within the strict criterion of $10^{-4}$ meV per cell.
This is achieved with a plane-wave energy cutoff of $500$ eV and by sampling
the electronic BZ using a $24\times 24$ $\Gamma$-centered $\bm{k}$-point grid.

For MgB2, the exchange-correlation functional is treated in the generalized
gradient approximation (GGA) parametrized by Perdew-Burke-Ernzerhoff (PBE)
[40]. The same energy convergence criterion as that in graphene is achieved
with a plane-wave energy cutoff of $500$ eV and by sampling the electronic BZ
using an $18\times 18\times 12$ $\Gamma$-centered $\bm{k}$-point grid.

####  V.1.2 Phonon dispersion and Kohn anomaly of graphene

Figure 5(a) shows the phonon dispersion of graphene along the boundary of its
IBZ, calculated using regular grids of different sizes as well as the Farey
grid combining them. Both acoustic and optical branches rapidly converge with
respect to the size of the $\bm{q}$-point grid in most regions of the high-
symmetry path, except for the window in the vicinity of the K point
highlighted in yellow shade and magnified in Fig. 5(b). The dispersion of the
highest frequency optical branch near the K point is very flat when the
dynamical matrix is sampled on the $4\times 4$, $5\times 5$ and $7\times 7$
regular $\bm{q}$-point grids. In contrast, it shows a strong
$\bm{q}$-dependence when the dynamical matrix is instead sampled on the
$6\times 6$ regular $\bm{q}$-point grid.

The fractional coordinates of the K point are $(1/3,1/3)$, which implies that
of all the regular grids considered, only the $6\times 6$ grid contains the
point $\bm{q}=\mathrm{K}$. This implies that the dynamical matrix at the K
point is only exact for this grid, and is obtained through Fourier
interpolation for the other considered grid sizes. In general, the dynamical
matrix at the K point will be exact for any regular $\bm{q}$-point grid of
size $3l\times 3l$ for any positive integer $l$, and interpolated for any
other grid size. Since Fourier interpolation does not tend to produce peaks or
pits that are not already represented by the explicit input data points, it is
expected that the regular grids in Fig. 5(a) that do not include the K point
miss the abrupt change of phonon frequency around that point.

This frequency drop is due to the abrupt change in the screening of lattice
vibrations by conduction electrons in metals and semi-metallic systems, which
is known as the Kohn anomaly [41]. For graphene, this Kohn anomaly has been
experimentally observed [42, 43] and theoretically characterized [16].
However, the vast majority of phonon calculations in the literature report a
flat dispersion near the K point [44, 45, 46, 47, 48, 49, 50, 51, 52]. Among
the above works, Ref. [47] uses the largest regular $\bm{q}$-point grid of
size $16\times 16$, but the Fourier interpolated frequency at the K point is
still significantly overestimated compared to the exact value. In our
simulations, we used regular grids up to a size of $20\times 20$, but the
frequency at the K point was still substantially overestimated. These results
demonstrate that it is computationally challenging to converge the phonon
frequency of the upper branch at the K point with respect to the regular grid
size.

However, as we have already noted earlier, it is not actually computationally
expensive to capture the exact phonon frequencies at the K point – a $3\times
3$ $\bm{q}$-grid can capture them exactly. As a result, in this example the
Farey grid has a clear advantage compared to regular grids: while regular
grids alternately sample the K point or not as the grid size grows, Farey
grids always guarantee the K point is sampled exactly when the grid order
exceeds $l=3$.

Figure 5: (color online). (a) Phonon dispersion of graphene along the high-
symmetry path
$\Gamma\rightarrow\mathrm{M}\rightarrow\mathrm{K}\rightarrow\Gamma$,
calculated with regular grids of different sizes (colored dotted line:
$l\times l$ where $l$ ranges from $4$ to $7$) and the Farey grid of order $7$
(black solid line: $\mathcal{K}^{(2)}_{\text{far}}\left(7\right)$). The shaded
region highlights the Kohn anomaly located near the K point, whose details are
magnified and shown in Figure (b). The phonon dispersion calculated with the
Farey grid is also compared to that calculated with even finer regular grids
(colored dashed line: $3l\times 3l$ where $l$ ranges from $4$ to $7$) and
shown in Figure (c).

In addition to the phonon frequency itself, the slope of the phonon dispersion
curve also has an important physical meaning: it is proportional to the phonon
momentum derivative of the dielectric function of the solid. Furthermore,
Piscanec and co-workers showed that, in graphene, the slope of the phonon
dispersion at the K point, where the Kohn anomaly is located, is proportional
to the square of the EPC strength [16]. The dynamical matrix of graphene is
nonanalytic at the K point, which leads to a nonzero slope. This nonanalytic
behavior cannot be captured with standard phonon calculations and Fourier
interpolation, as the dynamical matrix is assumed to be analytically dependent
on $\bm{q}$. In fact, if the dynamical matrix is analytical, symmetry enforces
the slope at the K point to always be zero in graphene. This means we cannot
directly access the slope at the K point, and instead we consider the slope at
points very close to the K point. To be precise, we use
$\bm{q}=\mathrm{K}+\delta\bm{q}$, where $|\delta\bm{q}|=1/10$ of the
$\mathrm{M}-\mathrm{K}$ distance. As shown in Fig. 5(c), We find that at least
an $18\times 18$ $\bm{q}$-point grid is necessary to converge this slope with
a regular grid. In contrast, the Farey grid of order 7,
$\mathcal{K}^{(2)}_{\text{far}}\left(7\right)$, which has the same sampling
density along the high-symmetry path
$\Gamma\rightarrow\mathrm{M}\rightarrow\mathrm{K}\rightarrow\Gamma$, leads to
the same level of convergence at a much smaller computational cost.

####  V.1.3 Phonon dispersion and Kohn anomalies of MgB2

We carry out a similar phonon calculation for the 3D system MgB2 to further
verify the superior performance of the Farey grid. The crystal structure of
MgB2 consists of hexagonal graphene-like planes of boron atoms intercalated
vertically with magnesium atoms located in the centers of the hexagons. This
simple binary compound has an exceptionally high superconducting transition
temperature $T_{\text{c}}=39$ K [53]. Another important feature is that the
phonons of the second highest frequency on the high-symmetry segment from
$\Gamma$ to A couple strongly to the electrons [54], which leads to a Kohn
anomaly as discussed below. These phonons correspond to a transverse mode with
$E_{2g}$ symmetry propagating along the vertical direction with the atomic
displacement parallel to the horizontal plane.

Figure 6 shows the phonon dispersion of MgB2 calculated with the Farey grid
$\mathcal{K}^{(3)}_{\text{far}}\left(6\right)$, as well as with the regular
grids whose sizes range from $4\times 4\times 4$ to $7\times 7\times 7$ for
comparison. We observe Kohn anomalies in the $E_{2g}$ optical branch near the
$\Gamma$ point (along the $\Gamma-\mathrm{K}$ and $\Gamma-\mathrm{M}$
directions) both in the phonon dispersions calculated with the Farey grid and
the regular grids whose sizes exceed $6\times 6\times 6$, which is in
agreement with previously reported results [55, 56, 57, 58]. In addition, we
also find Kohn anomalies on the first- and second-highest frequency phonon
branches along the $\Gamma-\mathrm{A}$ direction, and at $|\bm{q}|\approx
9/11$ of the $\Gamma-\mathrm{A}$ distance (highlighted in yellow shade and
magnified in Fig. 6(b)).

These Kohn anomalies are relatively weak compared to those around the $\Gamma$
point, but their presence has also been observed in inelastic x-ray scattering
experiments of the phonon dispersion [59]. From a computational point of view,
early first-principles calculations usually adopted regular $\bm{q}$-point
grids with a maximum size of $6\times 6\times 6$, and as a result overlooked
these Kohn anomalies [60, 55, 59, 61, 56, 57, 58]. The Kohn anomalies along
the $\Gamma-\mathrm{A}$ direction have only recently been correctly captured,
requiring a regular grid of size $12\times 12\times 12$ or larger [62, 63]. We
highlight that our Farey grid $\mathcal{K}^{(3)}_{\text{far}}\left(6\right)$
has the same sampling density on the $\Gamma-\mathrm{A}$ high-symmetry segment
and is able to capture the Kohn anomalies at a lower computational cost.
Finally, we note that Calandra and co-workers proposed a Wannier interpolation
scheme of the dynamical matrix which also successfully describes the Kohn
anomalies along the $\Gamma-\mathrm{A}$ direction [17]. They also show that
the presence of Kohn anomalies on the $E_{2g}$ branch in MgB2 can affect the
position of the main peak in the Eliashberg function, which is a crucial
quantity for calculating the superconducting transition temperature.

Figure 6: (color online). (a) Phonon dispersion of MgB2 along the high-
symmetry path
$\Gamma\rightarrow\mathrm{M}\rightarrow\mathrm{K}\rightarrow\Gamma\rightarrow\mathrm{A}$,
calculated with regular grids of different sizes (colored dotted line:
$l\times l\times l$ where $l$ ranges from $4$ to $7$) and the Farey grid of
order $6$ (black solid line: $\mathcal{K}^{(3)}_{\text{far}}\left(6\right)$).
The shaded region highlights the Kohn anomaly located near the A point, whose
details are magnified and shown in Figure (b). The left panel indicates the
lower and the right the higher mode.

####  V.1.4 Comparison between DFPT and FD with Farey grids

Table 1:  The computational expense of FD and DFPT calculations using regular and Farey grids. Timings are given as fractions of the time taken for the regular grid DFPT calculations. $\mathcal{K}_{\text{reg}}^{(\tau)}$ and $\mathcal{K}_{\text{far}}^{(\tau)}$ are the regular and Farey $\bm{q}$-grids necessary to obtain convergence of the Kohn anomaly. |  |  | FD | DFPT
---|---|---|---|---
| $\mathcal{K}^{(\tau)}_{\text{reg}}$ | $\mathcal{K}^{(\tau)}_{\text{far}}$ | t${}_{\text{reg}}$ | t${}_{\text{far}}$ | t${}_{\text{reg}}$ | t${}_{\text{far}}$
Graphene | $\mathcal{K}^{(2)}_{\text{reg}}(18)$ | $\mathcal{K}^{(2)}_{\text{far}}(7)$ | $1.456$ | $0.254$ | $1$ | $0.294$
MgB2 | $\mathcal{K}^{(3)}_{\text{reg}}(12)$ | $\mathcal{K}^{(3)}_{\text{far}}(6)$ | $3.894$ | $0.560$ | $1$ | $0.233$

The calculation of phonons using DFPT directly samples the $\bm{q}$-point grid
to construct the dynamical matrix, which can be accomplished with calculations
on a primitive cell. By contrast, the calculation of phonons using FD samples
the $\bm{q}$-point grid by constructing commensurate real space supercells. As
a result, DFPT phonon calculations are computationally cheaper than FD phonon
calculations when using regular grids. This implies that if both methods are
available, DFPT is typically the method of choice. The traditional advantage
of FD phonon calculations is that the formalism is simpler, which means that
it is available in some contexts in which DFPT phonon calculations are not, a
recent example being the calculations of phonons using DFT combined with
dynamical mean-field theory to study strongly correlated materials [64]. As
discussed in Sec. II, Farey grids provide a computational speedup for both
DFPT and FD, but the latter benefit from a larger speedup. In this section, we
show how the use of Farey grids leads to comparable computational expense
between DFPT and FD for the calculation of phonons, with examples in which FD
is computationally advantageous compared to DFPT.

In contrast to the rest of this work, we use the castep plane-wave DFT code
[35] for the comparisons carried out in this section. The reason for this is
that its implementation of DFPT is more efficient than that of the vasp
package, enabling a fairer comparison. We use norm-conserving
pseudopotentials, and target a converged electronic self-consistent field
energy per atom below $1$ meV.

For graphene, this is achieved with an energy cutoff of $1250$ eV, a
Monkhorst-Pack grid with adjacent points separated by $2\pi\times 0.03$ Å-1,
and a vacuum spacing of $9$ Å between graphene layers. For MgB2, the
corresponding quantities are a $900$ eV cutoff and a $\bm{k}$-point spacing of
$2\pi\times 0.03$ Å-1.

Table 1 compares the computational time required for different methods to
converge the phonon calculations of graphene and MgB2 to be able to capture
the various Kohn anomalies, with the reference cost set to $1$ for DFPT
calculations using a regular grid, which is the current state-of-the-art. As
expected, the computational cost of FD calculations using a regular grid is
significantly higher than the corresponding cost of DFPT calculations, with a
slowdown by a factor of $1.5$ for graphene and of $3.9$ for MgB2.

The results shown in Table 1 demonstrate that the use of Farey grids provides
a notable speedup compared to the use of regular grids. For DFPT, we find a
computational speedup by a factor of $3.4$ for graphene and of $4.3$ for MgB2.
This is due to a reduction in the number of $\bm{q}$-points necessary to
obtain a converged result with Farey grids, highlighting the extensive
computational savings enabled by an appropriate choice of the BZ sampling
points. With FD, the speedup is even larger, since, in addition to the number
of $\bm{q}$-points, the size of the largest supercells required in the
calculations is also reduced. Explicitly, the computational speedup compared
to regular-grid FD is $5.7$ for graphene and $7.0$ for MgB2.

It should be noted that the timings in Table 1 are taken only for the
converged grids in each case; as explained in Sec. II, a typical convergence
workflow for regular grids requires a significant amount of calculations with
insufficiently dense grids, which are ultimately discarded. This implies that
a realistic regular grid calculation requires more time than accounted for in
the timings presented here. By contrast, the convergence workflow for the
Farey grid does not discard any points, so that the timings presented here
correspond to a realistic calculation. As a result, the speedups obtained with
Farey grids compared to regular grids in Table 1 are a lower bound to the
expected speedups in a standard calculation that includes a convergence test.

It is remarkable that for graphene, the overall fastest method for a converged
phonon calculation is FD combined with Farey grids. The likely origin of this
is a favorable scaling of the supercell sizes in $2$D compared to 3D, combined
with the reduction of the number of $\bm{k}$-points required to sample the
electronic BZ of a supercell compared to a primitive cell. In this sense, we
can conclude that Farey grids unlock the full potential of the FD method. This
will prove particularly useful in calculations using electronic structure
methods for which DFPT is unavailable. This has so far been a bottleneck in
describing lattice dynamics with theories beyond semi-local DFT, such as meta-
GGA exchange-correlation functional DFT or even DFT combined with dynamical
mean-field theory, as DFPT requires mathematics specific to a given theory
[65]. FD, on the other hand, is completely agnostic to the underlying theory,
as long as it can calculate the forces on an atom. Our results significantly
broaden the range of applicability of FD, demonstrating that FD can become
competitive with DFPT. We expect this to open the doors to many new
applications.

###  V.2 Bandgap renormalization driven by electron-phonon coupling

Electron-phonon coupling dominates a variety of phenomena, from transport to
superconductivity. An area of recent interest is the temperature dependence of
electronic structures driven by EPC [66, 67, 68, 69, 70, 71], and we use this
as an example of calculating electron-phonon interactions from first
principles.

The temperature dependence of an electronic energy level $E_{n\bm{k}}$ can be
understood starting from the Allen-Heine-Cardona theory [72, 73, 74], which
states that

$E_{n\bm{k}}(T)=E_{n\bm{k}}^{\text{static}}+\frac{1}{2}\sum_{\nu\bm{q}}\frac{1}{2\omega_{\nu\bm{q}}}\frac{\partial^{2}E_{n\bm{k}}}{\partial\mu_{\nu\bm{q}}^{2}}\left[1+2n_{\mathrm{B}}\left(\omega_{\nu\bm{q}},T\right)\right],$
(26)

where $E_{n\bm{k}}^{\text{static}}$ is the static electronic eigenenergy,
$\mu_{\nu\bm{q}}$ is the phonon amplitude that one can directly obtain from
the phonon polarization $\bm{\xi}_{\nu\bm{q}}$ [75], and
$n_{\mathrm{B}}\left(\omega_{\nu\bm{q}},T\right)=\left(\mathrm{e}^{\omega_{\nu\bm{q}}/k_{\mathrm{B}}T}-1\right)^{-1}$
is the Bose-Einstein factor. It is worth noting that there is a difference
between the static energy band $E_{n\bm{k}}^{\text{static}}$ and the zero-
temperature energy band $E_{n\bm{k}}(T=0)$, arising from zero-point quantum
fluctuations and typically referred to as the zero-point correction.

The numerical implementation of Eq. (26) requires the solution of the lattice
dynamics of the system, and the subsequent evaluation of the second
derivatives of the electronic eigenenergy of interest with respect to all
possible phonon displacements. In practice, the calculations are performed on
a $\bm{q}$-point grid of the BZ. However, this procedure often requires a very
fine $\bm{q}$-point grid whose corresponding computational cost can become
prohibitive, particularly in the context of FD where the grid is sampled
through commensurate real space supercells. To understand the requirement of a
fine $\bm{q}$-point grid, we invoke the Hellmann-Feynman theorem to split the
second derivative in Eq. (26) into three terms:

$\displaystyle\frac{\partial^{2}E_{n\bm{k}}}{\partial\mu_{\nu\bm{q}}^{2}}$
$\displaystyle=\left\langle\varphi_{n\bm{k}}\left|\frac{\partial^{2}\mathcal{H}_{\mathrm{el}}}{\partial\mu_{\nu\bm{q}}^{2}}\right|\varphi_{n\bm{k}}\right\rangle$
(27)
$\displaystyle+\left\langle\varphi_{n\bm{k}}\left|\frac{\partial\mathcal{H}_{\mathrm{el}}}{\partial\mu_{\nu\bm{q}}}\right|\frac{\partial\varphi_{n\bm{k}}}{\partial\mu_{\nu\bm{q}}}\right\rangle+\left\langle\frac{\partial\varphi_{n\bm{k}}}{\partial\mu_{\nu\bm{q}}}\left|\frac{\partial\mathcal{H}_{\mathrm{el}}}{\partial\mu_{\nu\bm{q}}}\right|\varphi_{n\bm{k}}\right\rangle.$

The derivative of the state appearing in Eq. (27) can be expanded in the
complete (unperturbed) basis $\\{|\varphi_{n\bm{k}}\rangle\\}$ according to
perturbation theory:

$\left|\frac{\partial\varphi_{n\bm{k}}}{\partial\mu_{\nu\bm{q}}}\right\rangle=\sum_{\tiny{\begin{array}[]{c}(m,\bm{k^{\prime}})\\\
\neq(n,\bm{k})\end{array}}}\left|\varphi_{m\bm{k}^{\prime}}\right\rangle\frac{\left\langle\varphi_{m\bm{k}^{\prime}}\left|\frac{\partial\mathcal{H}_{\mathrm{el}}}{\partial\mu_{\nu\bm{q}}}\right|\varphi_{n\bm{k}}\right\rangle}{E_{n\bm{k}}^{(0)}-E_{m\bm{k}^{\prime}}^{(0)}}.$
(28)

Substituting Eq. (27) and Eq. (28) into Eq. (26) yields

$E_{n\bm{k}}(T)=E_{n\bm{k}}^{\text{static}}+\frac{1}{N}\sum_{\nu\bm{q}}\left[\sum_{m}\frac{\left|g_{mn\nu}(\bm{k},\bm{q})\right|^{2}}{E_{n\bm{k}}^{(0)}-E_{m\bm{k}+\bm{q}}^{(0)}}+g_{nn\nu\nu}^{\mathrm{DW}}(\bm{k},\bm{q},-\bm{q})\right]\left(1+2n_{\mathrm{B}}\left(\omega_{\nu\bm{q}},T\right)\right),$
(29)

where $g_{mn\nu}(\bm{k},\bm{q})$ and
$g_{nn\nu\nu}^{\mathrm{DW}}(\bm{k},\bm{q},-\bm{q})$ are the EPC matrix
elements of first and second order respectively, defined as in Ref [65]. It is
worth noting that the sum excludes $\bm{q}=\bm{0}$ for the first-order term
when $m=n$, but the denominator $E_{n\bm{k}}^{(0)}-E_{m\bm{k+q}}^{(0)}$ can
still tend to $0$ for some $\bm{q}\neq\bm{0}$. When this happens, Equation
(26) becomes an improper integral over the BZ: the integrand diverges at some
$\bm{q}$-point(s), leading to a sharp peak in the EPC. A converged calculation
requires that this singular behavior is correctly captured, which can only be
accomplished if the distance between two adjacent sampling points is less than
the width of these peaks. This implies that fine sampling grids become
necessary.

In the following, we evaluate the temperature dependence of the band gaps of
diamond and bismuthene.

####  V.2.1 Computational details

For diamond, the exchange-correlation functional is treated in the LDA [39].
The energy is converged to within the criterion of $10^{-4}$ meV per cell with
a plane-wave energy cutoff of $500$ eV and by sampling the electronic BZ using
an $11\times 11\times 11$ $\Gamma$-centered $\bm{k}$-point grid.

For bismuthene, the exchange-correlation functional is modeled using the PBE
GGA parametrization [40]. A vacuum layer with an even larger thickness than
that in the graphene calculation, $20$ Å, is used given the buckling structure
in bismuthene. The same energy convergence criterion as that in diamond is
achieved with a plane-wave energy cutoff of $300$ eV and by sampling the
electronic BZ using an $11\times 11$ $\Gamma$-centered $\bm{k}$-point grid.
The spin-orbital interaction is included in the calculations via a
perturbation to the scalar relativistic Hamiltonian [76].

####  V.2.2 Temperature dependence of the indirect band gap of diamond

Diamond has been the testbed for the development of computational tools to
study the electron-phonon driven temperature dependence of band structures
[66, 77, 78, 69, 79, 24]. It is a particularly useful test case for our method
in $3$D, owing to the fact that it is possible to obtain converged results for
its thermal (indirect) band gap using a regular grid with both DFPT [78, 69]
and FD [79, 24]. Therefore, we can carry out a direct comparison between the
predictions from the regular and Farey grids.

Figure 7: Zero point renormalization of the thermal band gap of diamond as a
function of the scale of the regular (black squares) and Farey (red circles)
$\bm{q}$-point grids. Inset: static band structure of diamond along the high-
symmetry path
$W\rightarrow\Gamma\rightarrow\mathrm{L}\rightarrow\mathrm{U}\rightarrow\mathrm{X}\rightarrow\Gamma$.

The inset of Fig. 7 shows the static band structure of diamond, highlighting
its wide thermal band gap, $E_{\mathrm{gap}}(\mu_{\nu\bm{q}}=0)=4.165$ eV.
This result, in a manner typical of semi-local DFT, underestimates the
experimental value of $5.450$ eV [80]. The valence band maximum is located at
the $\Gamma$ point and the conduction band minimum lies on the high-symmetry
line between $\Gamma$ and $\mathrm{X}$. Figure 7 also shows the
renormalization of this band gap at $0$ K due to nuclear zero-point motion,
calculated with regular grids, $\mathcal{K}^{(3)}_{\text{reg}}\left(l\right)$,
and Farey grids, $\mathcal{K}^{(3)}_{\text{far}}\left(l\right)$, of different
sizes: $l$ in both cases ranges from $1$ to $20$. The Farey grid of order $l$
is constructed by combining the EPC calculations of all regular grids with
sizes up to $l$.

The results obtained from the Farey grid exhibit a smoother convergence
compared to those obtained with the regular grid. The zero-point
renormalization is converged to $-325\pm 5$ meV when using the Farey grid
$\mathcal{K}^{(3)}_{\text{far}}\left(6\right)$. This uncertainty is an order
of magnitude smaller than the change to the zero-point renormalization due to
the choice of pseudopotential for carbon [81]. Using the regular grid, a
$12\times 12\times 12$ grid is needed to achieve the same level of
convergence. It is worth noting that the grid size of the Farey grid here
($|\mathcal{K}^{(3)}_{\text{far}}\left(6\right)|=396$) is less than a quarter
of that of the regular grid
($|\mathcal{K}^{(3)}_{\text{reg}}\left(12\right)|=1728$), and the actual
computational cost of the Farey grid is approximately a sixth that of the
regular grid. This computational cost reduction comes not only from reducing
the number of required sampling points, but also from the fact that the
largest supercells needed for the Farey grid calculations are only 6 times
larger than the primitive cell, whereas we require supercells that are 12
times larger than the primitive cell in the calculation with a $12\times
12\times 12$ regular grid. It is also clear that the regular and Farey grids
converge to the same value as $l$ increases.

Overall, the converged result from the first-principles calculations is very
close to the experimental measurement of $-340$ meV [82]. The small difference
compared to the experimental value is likely to be due to the use of the LDA
to describe the exchange-correlation energy. In principle, this can be
overcome with more expensive calculations using hybrid exchange-correlation
functionals or by invoking the $GW$ approximation from many-body perturbation
theory. Recently, Karsai and co-workers evaluated the renormalization of the
diamond thermal band gap using the $G_{0}W_{0}$ approximation, finding that
$G_{0}W_{0}$ further increases the magnitude of the zero-point renormalization
by 17 meV [70]. Adding this value to our converged DFT result, one obtains a
theoretical value very close to the experimental result.

We also note that thermal expansion can also renormalize the electronic energy
levels of materials, but in the case of diamond the effect is minimal compared
to the changes driven by EPC, so it can be safely neglected [68].

####  V.2.3 Topologically nontrivial band gap of bismuthene

Bismuthene is a 2D topological insulator with a topologically nontrivial band
gap at the BZ center due to its intrinsic strong spin-orbital coupling [83,
84, 85, 86, 87]. Motivated by the search for topological insulators which can
work at ambient temperature ($T=300$ K), here we investigate the EPC-driven
renormalization of this topological band gap. The top panel of Fig. 8 shows
the static band structure of bismuthene near the Fermi level with the spin-
orbital coupling interaction taken into account. The conduction band has
conical shape in the zone center, while the highest valence band exhibits a
double peak around the $\Gamma$-point, with a trough exactly at it. This
results in a sizable static $\Gamma$-point band gap of
$E_{\mathrm{gap}}(\mu_{\nu\bm{q}}=0)=572$ meV.

Figure 8: (color online). (top panel) Static band structure of bismuthene
(where the Fermi level is set to zero and the dashed line indicates the
valence band energy at the $\Gamma$-point) and (middle panel) the EPC strength
along the high-symmetry path
$\mathrm{M}\rightarrow\Gamma\rightarrow\mathrm{K}$ (where the symbols on the
line mark the $\bm{q}$-points that can be explicitly assessed by the regular
$\bm{q}$-point grids of different sizes). The bottom panel highlights the EPC
strength ‘seen’ by regular $\bm{q}$-point grids of sizes $16\times 16$,
$17\times 17$ and $18\times 18$. Figure 9: (color online). Renormalization of
the topologically nontrivial band gap of bismuthene at room temperature
($T=300$ K) as a function of the scale of the regular (black squares) and the
truncated Farey (red circles) $\bm{q}$-point grids.

Figure 9 shows the renormalization of this band gap at $300$ K calculated with
a regular and a truncated Farey $\bm{q}$-point grids. Its value converges to
$18.8\pm 0.05$ meV using the truncated Farey $\bm{q}$-point grid generated by
$\mathcal{K}^{(2)}_{\text{far}}\left(10,18\right)$. This stands in stark
contrast with the results calculated with the regular $\bm{q}$-point grids,
which oscillate so strongly as to prevent even qualitatively consistent
conclusions. For instance, while the calculation with the $14\times 14$
regular $\bm{q}$-point grid suggests the band gap is almost independent of
temperature, the calculation with the $16\times 16$ ($18\times 18$) regular
$\bm{q}$-point grid implies that it would be significantly lower (higher) at
room temperature. Even at a grid size of $21\times 21$, the largest regular
grid considered here, no clear convergence trend can be observed.

To obtain more insight into these difficulties with converging the EPC
calculations in bismuthene using regular grids, we define the EPC amplitude in
the context of the Allen-Heine-Cardona formula [Eq. (26)] as

$S^{\mathrm{el-
ph}}_{\nu\bm{q}}=\frac{\partial^{2}E_{\mathrm{gap}}}{\partial\mu_{\nu\bm{q}}^{2}},$
(30)

and the $\bm{q}$-resolved EPC amplitude:

$S^{\mathrm{el-ph}}(\bm{q})=\sum_{\nu}S^{\mathrm{el-ph}}_{\nu\bm{q}}.$ (31)

The middle panel of Fig. 8 shows the $\bm{q}$-resolved EPC strength along the
high-symmetry path $\mathrm{M}-\Gamma-\mathrm{K}$, where it can be clearly
seen that the curve exhibits strong peaks at around $1/4$ of the
$\Gamma-\mathrm{M}$ distance and $1/5$ of the $\Gamma-\mathrm{K}$ distance.
This agrees with the perturbation picture provided by Eq. (29): the double-
peak shape of the valence band of bismuthene favors processes where a
$\Gamma$-point electron scatters low-frequency phonons with wavevectors
$|\bm{q}|\approx 1/4$ of the $\Gamma-\mathrm{M}$ distance or $|\bm{q}|\approx
1/5$ of the $\Gamma-\mathrm{K}$ distance. These phonons scatter the electron
to another state at a similar energy (as shown in the top panel), but with
different momentum corresponding to that contributed by the phonon. More
generally, we find that the EPC amplitude is almost isotropic around $\Gamma$,
which implies that the band gap renormalization is dominated by low-frequency
phonons with wavevectors $|\bm{q}|\approx b/8$, where $b$ is the lattice
constant of the bismuthene reciprocal lattice.

The bottom panel of Fig. 8 shows the EPC strength ‘seen’ by regular
$\bm{q}$-point grids of a range of different sizes, which explains the
qualitatively inconsistent results observed in Fig. 9. The $18\times 18$
uniform $\bm{q}$-point grid overestimates the EPC amplitude of the system as
it only samples the peak points of the curve, while the $16\times 16$ uniform
$\bm{q}$-point grid underestimates it due to only sampling the valley points.
As shown in the middle panel of Fig. 8, to capture this peak, it is necessary
to use a regular grid that at least contains
$\mathrm{lcm}(16,17,18)\times\mathrm{lcm}(16,17,18)$ points, that is, a
$2448\times 2448$ regular $\bm{q}$-point grid. This grid size is beyond the
capabilities of modern high-performance computing. By comparison, the Farey
grid leads to converged results at a moderate computational cost.

It should also be noted that, although the sharpest peaks in the middle panel
of Fig. 8 agree very well with the expectation of Eq. (29), there are two
smaller peaks in the EPC amplitude appearing at around $|\bm{q}|\approx 1/2$
of the $\Gamma-\mathrm{M}$ distance and $|\bm{q}|\approx 3/8$ of the
$\Gamma-\mathrm{K}$ distance. These cannot be explained by the perturbation
picture, and we infer that this behavior is mostly likely subject to details
of the EPC matrix element and high-order EPC interaction terms.

More generally, the temperature dependence of the band gap of bismuthene has
not been reported elsewhere. Our results show that bismuthene will have a
topological band of almost $600$ meV at $300$ K, making it an ideal material
to study topological phenomena at room temperature. We also note that our
calculations do not include thermal expansion, but the known dependence of the
band gap of bismuthene on the lattice parameter [88, 89] implies that the
inclusion of thermal expansion should not change our conclusions.

###  V.3 Discussion and outlook

First principles calculations of the electronic structure of materials have
reached a high level of sophistication, enabling the description of electrons
moving in a static lattice with high accuracy. Examples of this are hybrid
functional DFT, which accurately captures the exchange-correlation energy [90,
91], the $GW$ method [92, 93], which includes many-body effects, and the DFT+U
method [94, 95] and dynamical mean-field theory [96, 97], which treat strong
correlation effects. These developments have laid a solid foundation for
studying the coupling between phonons and electrons in a wide variety of
materials from first principles. In addition, with advances in the study of
other quasiparticles, such as the Heisenberg ferromagnet equation for magnons
[98, 99] and the Bethe-Salpeter equation for excitons [100, 101], one can
investigate phenomena involving magnon-phonon and exciton-phonon coupling. In
most cases, the lattice dynamics and related properties predicted by the
methods mentioned in this section can only be calculated using the FD method.
By reducing the computational cost associated with phonon calculations
dramatically, the Farey grid and the SAVT technique proposed in this paper
pave the way for realistic implementations using these state-of-the-art
methods to carry out electron-phonon (or magnon-phonon or exciton-phonon)
coupling calculations.

##  VI Conclusion

We have introduced two developments to carry out BZ integration. First, a
nonuniform grid, called the Farey grid, which is computationally advantageous
compared to the widely used regular uniform grids. Second, the symmetry-
adapted Voronoi tessellation technique to evaluate the weights necessary for a
numerical BZ integration using an arbitrary grid. We have demonstrated that
these developments provide significant computational speedups in the
calculation of phonons and EPC, irrespective of whether DFPT or finite
difference methods are used. However, the most dramatic advantages are
observed for finite difference methods, which have traditionally been slower
than DFPT, but that now become competitive from a computational expense
perspective. Specifically, we apply our methodology to the study of Kohn
anomalies in graphene and MgB2, and to the study of electron-phonon driven
band gap renormalizations in diamond and bismuthene. The results for
bismuthene reveal that this material should have a sizeable topological band
gap at room temperature, making it an ideal platform to study room temperature
topological phenomena.

A practical insight of our work is that the Farey grid can be constructed by a
combination of differently-sized regular grids. This implies that it should be
relatively straightforward to adapt modern first principles software, which is
typically based on regular grids, to incorporate the Farey grid.

An important conclusion from our work is that, in many problems of physical
interest, the regular grid is far from the most efficient sampling method: it
samples regions in the BZ where the integrand varies slowly just as densely as
those where the integrand varies rapidly. Ideally, however, one would use a
few points to represent the former, while focusing more computational power on
accurately representing the latter. This is a possibility that opens up as a
result of the symmetry-adapted Voronoi tessellation technique. Although we
have only used it in conjunction with the nonuniform Farey grid, it is
generally applicable to any grid (uniform or not) that one may conceive to
evaluate the BZ integration. Thus, it provides a solid platform on which to
design the most appropriate nonuniform grid for any given task, with the
prospect of saving valuable computational time and allowing for new frontiers
in computational condensed matter physics to be charted.

###### Acknowledgements.

S.C. acknowledges financial support from the Cambridge Trust and from the
Winton Programme for the Physics of Sustainability. P.T.S. gratefully
acknowledges funding from the Department of Materials Science and Metallurgy
at the University of Cambridge and from a Trinity Hall Research Studentship.
P.T.S. also thanks Chris J. Pickard for helpful discussions and support. B.M.
acknowledges support from a UKRI Future Leaders Fellowship (Grant No.
MR/V023926/1), from the Gianna Angelopoulos Programme for Science, Technology,
and Innovation, and from the Winton Programme for the Physics of
Sustainability. The calculations in this work have been performed using
resources provided by the Cambridge Tier-2 system (operated by the University
of Cambridge Research Computing Service and funded by EPSRC [EP/P020259/1]),
as well as by the UK Materials and Molecular Modelling Hub (partially funded
by EPSRC [EP/P020194]), Thomas, and by the UK National Supercomputing Service,
ARCHER. Access to Thomas and ARCHER was obtained via the UKCP consortium and
funded by EPSRC [EP/P022561/1].

## References

  * Baldereschi [1973] A. Baldereschi, Mean-value point in the Brillouin zone, Phys. Rev. B 7, 5212 (1973).
  * Chadi and Cohen [1973] D. J. Chadi and M. L. Cohen, Special points in the Brillouin zone, Phys. Rev. B 8, 5747 (1973).
  * Cunningham [1974] S. L. Cunningham, Special points in the two-dimensional Brillouin zone, Phys. Rev. B 10, 4988 (1974).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Fehlner and Vosko [1977] W. R. Fehlner and S. H. Vosko, Special points for the evaluation of form factors and fourier series representations in hexagonal crystals, Can. J. Phys. 55, 2041 (1977).
  * Froyen [1989] S. Froyen, Brillouin-zone integration by fourier quadrature: Special points for superlattice and supercell calculations, Phys. Rev. B 39, 3168 (1989).
  * Moreno and Soler [1992] J. Moreno and J. M. Soler, Optimal meshes for integrals in real- and reciprocal-space unit cells, Phys. Rev. B 45, 13891 (1992).
  * Wisesa  _et al._ [2016] P. Wisesa, K. A. McGill, and T. Mueller, Efficient generation of generalized monkhorst-pack grids through the use of informatics, Phys. Rev. B 93, 155109 (2016).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996).
  * Alfè [2009] D. Alfè, Phon: A program to calculate phonons using the small displacement method, Comput. Phys. Commun. 180, 2622 (2009).
  * Marini  _et al._ [2009] A. Marini, C. Hogan, M. Grüning, and D. Varsano, yambo: An ab initio tool for excited state calculations, Comput. Phys. Commun. 180, 1392 (2009).
  * Tancogne-Dejean  _et al._ [2020] N. Tancogne-Dejean, F. G. Eich, and A. Rubio, Time-dependent magnons from first principles, J. Chem. Theory Comput. 16, 1007 (2020).
  * Sabiryanov and Jaswal [1999] R. F. Sabiryanov and S. S. Jaswal, Magnons and magnon-phonon interactions in iron, Phys. Rev. Lett. 83, 2062 (1999).
  * Poncé  _et al._ [2016] S. Poncé, E. Margine, C. Verdi, and F. Giustino, EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized wannier functions, Comput. Phys. Commun. 209, 116 (2016).
  * Antonius and Louie [2022] G. Antonius and S. G. Louie, Theory of exciton-phonon coupling, Phys. Rev. B 105, 085111 (2022).
  * Piscanec  _et al._ [2004] S. Piscanec, M. Lazzeri, F. Mauri, A. C. Ferrari, and J. Robertson, Kohn anomalies and electron-phonon interactions in graphite, Phys. Rev. Lett. 93, 185503 (2004).
  * Calandra  _et al._ [2010] M. Calandra, G. Profeta, and F. Mauri, Adiabatic and nonadiabatic phonon dispersion in a wannier function approach, Phys. Rev. B 82, 165111 (2010).
  * Giustino  _et al._ [2007] F. Giustino, M. L. Cohen, and S. G. Louie, Electron-phonon interaction using wannier functions, Phys. Rev. B 76, 165108 (2007).
  * Sjakste  _et al._ [2015] J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Wannier interpolation of the electron-phonon matrix elements in polar semiconductors: Polar-optical coupling in gaas, Phys. Rev. B 92, 054307 (2015).
  * Verdi and Giustino [2015] C. Verdi and F. Giustino, Fröhlich electron-phonon vertex from first principles, Phys. Rev. Lett. 115, 176401 (2015).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications, Rev. Mod. Phys. 84, 1419 (2012).
  * Zhou  _et al._ [2021] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Perturbo: A software package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics, Comput. Phys. Commun. 264, 107970 (2021).
  * Cepellotti  _et al._ [2021] A. Cepellotti, J. Coulter, A. Johansson, N. S. Fedorova, and B. Kozinsky, Phoebe: a collection of phonon and electron boltzmann equation solvers, arXiv:2111.14999 (2021).
  * Lloyd-Williams and Monserrat [2015] J. H. Lloyd-Williams and B. Monserrat, Lattice dynamics and electron-phonon coupling calculations using nondiagonal supercells, Phys. Rev. B 92, 184301 (2015).
  * Niven  _et al._ [1991] I. Niven, H. S. Zuckerman, and H. L. Montgomery,  _An introduction to the theory of numbers_ (John Wiley & Sons, 1991) pp. 297–300.
  * Ledoan [2018] A. Ledoan, The discrepancy of farey series, Acta Math. Hungar 156, 465 (2018).
  * Kleinman [1983] L. Kleinman, Error in the tetrahedron integration scheme, Phys. Rev. B 28, 1139 (1983).
  * Jepsen and Andersen [1984] O. Jepsen and O. K. Andersen, No error in the tetrahedron integration scheme, Phys. Rev. B 29, 5965 (1984).
  * Blöchl  _et al._ [1994] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B 49, 16223 (1994).
  * Sibson [1981] R. Sibson, A brief description of natural neighbour interpolation, Interpreting multivariate data (1981).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Clark  _et al._ [2005] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson, and M. C. Payne, First principles methods using CASTEP, Z. Kristallogr. Cryst. Mater. 220, 567 (2005).
  * Hamann  _et al._ [1979] D. R. Hamann, M. Schlüter, and C. Chiang, Norm-conserving pseudopotentials, Phys. Rev. Lett. 43, 1494 (1979).
  * Parlinski  _et al._ [1997] K. Parlinski, Z. Q. Li, and Y. Kawazoe, First-principles determination of the soft mode in cubic ${\mathrm{zro}}_{2}$, Phys. Rev. Lett. 78, 4063 (1997).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Kohn [1959] W. Kohn, Image of the fermi surface in the vibration spectrum of a metal, Phys. Rev. Lett. 2, 393 (1959).
  * Maultzsch  _et al._ [2004] J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phonon dispersion in graphite, Phys. Rev. Lett. 92, 075501 (2004).
  * Mohr  _et al._ [2007] M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Phonon dispersion of graphite by inelastic X-ray scattering, Phys. Rev. B 76, 035439 (2007).
  * Sánchez-Portal  _et al._ [1999] D. Sánchez-Portal, E. Artacho, J. M. Soler, A. Rubio, and P. Ordejón, Ab initio structural, elastic, and vibrational properties of carbon nanotubes, Phys. Rev. B 59, 12678 (1999).
  * Dubay and Kresse [2003] O. Dubay and G. Kresse, Accurate density functional calculations for the phonon dispersion relations of graphite layer and carbon nanotubes, Phys. Rev. B 67, 035401 (2003).
  * Yanagisawa  _et al._ [2005] H. Yanagisawa, T. Tanaka, Y. Ishida, M. Matsue, E. Rokuta, S. Otani, and C. Oshima, Analysis of phonons in graphene sheets by means of hreels measurement and ab initio calculation, Surf. Interface Anal. 37, 133 (2005).
  * Mounet and Marzari [2005] N. Mounet and N. Marzari, First-principles determination of the structural, vibrational and thermodynamic properties of diamond, graphite, and derivatives, Phys. Rev. B 71, 205214 (2005).
  * Zhang  _et al._ [2011] H. Zhang, G. Lee, and K. Cho, Thermal transport in graphene and effects of vacancy defects, Phys. Rev. B 84, 115460 (2011).
  * Marquina  _et al._ [2013] J. Marquina, C. Power, J. Gonzalez, and J.-M. Broto, Ab initio study of the electronic and vibrational properties of 1-nm-diameter single-walled nanotubes, Adv. Mater. Phys. Chem. 03, 178 (2013).
  * Slotman  _et al._ [2014] G. J. Slotman, G. A. de Wijs, A. Fasolino, and M. I. Katsnelson, Phonons and electron-phonon coupling in graphene-h-BN heterostructures, Ann. Phys. 526, 381 (2014).
  * Gu and Yang [2015] X. Gu and R. Yang, First-principles prediction of phononic thermal conductivity of silicene: A comparison with graphene, J. Appl. Phys. 117, 025102 (2015).
  * Diery  _et al._ [2018] W. Diery, E. A. Moujaes, and R. Nunes, Nature of localized phonon modes of tilt grain boundaries in graphene, Carbon 140, 250 (2018).
  * Nagamatsu  _et al._ [2001] J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu, Superconductivity at 39 K in magnesium diboride, Nature 410, 63 (2001).
  * Kong  _et al._ [2001] Y. Kong, O. V. Dolgov, O. Jepsen, and O. K. Andersen, Electron-phonon interaction in the normal and superconducting states of MgB2, Phys. Rev. B 64, 020501(R) (2001).
  * Bohnen  _et al._ [2001] K.-P. Bohnen, R. Heid, and B. Renker, Phonon dispersion and electron-phonon coupling in MgB2 and AlB2, Phys. Rev. Lett. 86, 5771 (2001).
  * d’Astuto  _et al._ [2007] M. d’Astuto, M. Calandra, S. Reich, A. Shukla, M. Lazzeri, F. Mauri, J. Karpinski, N. D. Zhigadlo, A. Bossak, and M. Krisch, Weak anharmonic effects in $\mathrm{Mg}{\mathrm{b}}_{2}$: A comparative inelastic X-ray scattering and raman study, Phys. Rev. B 75, 174508 (2007).
  * Alarco  _et al._ [2014] J. A. Alarco, A. Chou, P. C. Talbot, and I. D. Mackinnon, Phonon modes of MgB2: super-lattice structures and spectral response, Phys. Chem. Chem. Phys. 16, 24443 (2014).
  * Mackinnon  _et al._ [2017] I. D. Mackinnon, P. C. Talbot, and J. A. Alarco, Phonon dispersion anomalies and superconductivity in metal substituted MgB2, Comput. Mater. Sci. 130, 191 (2017).
  * Shukla  _et al._ [2003] A. Shukla, M. Calandra, M. d’Astuto, M. Lazzeri, F. Mauri, C. Bellin, M. Krisch, J. Karpinski, S. M. Kazakov, J. Jun, D. Daghero, and K. Parlinski, Phonon dispersion and lifetimes in MgB2, Phys. Rev. Lett. 90, 095506 (2003).
  * Yildirim  _et al._ [2001] T. Yildirim, O. Gülseren, J. W. Lynn, C. M. Brown, T. J. Udovic, Q. Huang, N. Rogado, K. A. Regan, M. A. Hayward, J. S. Slusky, T. He, M. K. Haas, P. Khalifah, K. Inumaru, and R. J. Cava, Giant anharmonicity and nonlinear electron-phonon coupling in MgB2: A combined first-principles calculation and neutron scattering study, Phys. Rev. Lett. 87, 037001 (2001).
  * Baron  _et al._ [2004] A. Q. R. Baron, H. Uchiyama, Y. Tanaka, S. Tsutsui, D. Ishikawa, S. Lee, R. Heid, K.-P. Bohnen, S. Tajima, and T. Ishikawa, Kohn anomaly in MgB2 by inelastic X-ray scattering, Phys. Rev. Lett. 92, 197004 (2004).
  * Pešić  _et al._ [2019] J. Pešić, I. Popov, A. Šolajić, V. Damljanović, K. Hingerl, M. Belić, and R. Gajić, Ab initio study of the electronic, vibrational, and mechanical properties of the magnesium diboride monolayer, Condens. Matter 4, 37 (2019).
  * Novko  _et al._ [2020] D. Novko, F. Caruso, C. Draxl, and E. Cappelluti, Ultrafast hot phonon dynamics in ${\mathrm{mgb}}_{2}$ driven by anisotropic electron-phonon coupling, Phys. Rev. Lett. 124, 077001 (2020).
  * Koçer  _et al._ [2020] C. P. Koçer, K. Haule, G. L. Pascut, and B. Monserrat, Efficient lattice dynamics calculations for correlated materials with $\mathrm{DFT}+\mathrm{DMFT}$, Phys. Rev. B 102, 245104 (2020).
  * Giustino [2017] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 015003 (2017).
  * Giustino  _et al._ [2010] F. Giustino, S. G. Louie, and M. L. Cohen, Electron-phonon renormalization of the direct band gap of diamond, Phys. Rev. Lett. 105, 265501 (2010).
  * Bhosale  _et al._ [2012] J. Bhosale, A. K. Ramdas, A. Burger, A. Muñoz, A. H. Romero, M. Cardona, R. Lauck, and R. K. Kremer, Temperature dependence of band gaps in semiconductors: Electron-phonon interaction, Phys. Rev. B 86, 195208 (2012).
  * Monserrat  _et al._ [2013] B. Monserrat, N. D. Drummond, and R. J. Needs, Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress, Phys. Rev. B 87, 144302 (2013).
  * Poncé  _et al._ [2015] S. Poncé, Y. Gillet, J. Laflamme Janssen, A. Marini, M. Verstraete, and X. Gonze, Temperature dependence of the electronic structure of semiconductors and insulators, J. Chem. Phys. 143, 102813 (2015).
  * Karsai  _et al._ [2018] F. Karsai, M. Engel, E. Flage-Larsen, and G. Kresse, Electron–phonon coupling in semiconductors within the GW approximation, New J. Phys. 20, 123008 (2018).
  * Zhang  _et al._ [2020] Y. Zhang, Z. Wang, J. Xi, and J. Yang, Temperature-dependent band gaps in several semiconductors: from the role of electron–phonon renormalization, J. Condens. Matter Phys. 32, 475503 (2020).
  * Allen and Heine [1976] P. Allen and V. Heine, Theory of the temperature dependence of electronic band structures, J. Phys. C: Solid State Phys. 9, 2305 (1976).
  * Allen and Cardona [1981] P. B. Allen and M. Cardona, Theory of the temperature dependence of the direct gap of germanium, Phys. Rev. B 23, 1495 (1981).
  * Allen and Cardona [1983] P. B. Allen and M. Cardona, Temperature dependence of the direct gap of Si and Ge, Phys. Rev. B 27, 4760 (1983).
  * Monserrat [2018] B. Monserrat, Electron–phonon coupling from finite differences, J. Condens. Matter Phys. 30, 083001 (2018).
  * Koelling and Harmon [1977] D. D. Koelling and B. N. Harmon, A technique for relativistic spin-polarised calculations, J. Phys. C: Solid State Phys. 10, 3107 (1977).
  * Zacharias and Giustino [2016] M. Zacharias and F. Giustino, One-shot calculation of temperature-dependent optical spectra and phonon-induced band-gap renormalization, Phys. Rev. B 94, 075125 (2016).
  * Poncé  _et al._ [2014] S. Poncé, G. Antonius, Y. Gillet, P. Boulanger, J. Laflamme Janssen, A. Marini, M. Côté, and X. Gonze, Temperature dependence of electronic eigenenergies in the adiabatic harmonic approximation, Phys. Rev. B 90, 214304 (2014).
  * Monserrat and Needs [2014] B. Monserrat and R. J. Needs, Comparing electron-phonon coupling strength in diamond, silicon, and silicon carbide: First-principles study, Phys. Rev. B 89, 214304 (2014).
  * O’Donnell and Chen [1991] K. P. O’Donnell and X. Chen, Temperature dependence of semiconductor band gaps, Appl. Phys. Lett. 58, 2924 (1991).
  * Poncé  _et al._ [2014] S. Poncé, G. Antonius, P. Boulanger, E. Cannuccia, A. Marini, M. Côté, and X. Gonze, Verification of first-principles codes: Comparison of total energies, phonon frequencies, electron–phonon coupling and zero-point motion correction to the gap between ABINIT and QE/Yambo, Comput. Mater. Sci. 83, 341 (2014).
  * Cardona [2001] M. Cardona, Renormalization of the optical response of semiconductors by electron–phonon interaction, Phys. Status Solidi A 188, 1209 (2001).
  * Murakami [2006] S. Murakami, Quantum spin hall effect and enhanced magnetic response by spin-orbit coupling, Phys. Rev. Lett. 97, 236805 (2006).
  * Koroteev  _et al._ [2008] Y. M. Koroteev, G. Bihlmayer, E. V. Chulkov, and S. Blügel, First-principles investigation of structural and electronic properties of ultrathin Bi films, Phys. Rev. B 77, 045428 (2008).
  * Wada  _et al._ [2011] M. Wada, S. Murakami, F. Freimuth, and G. Bihlmayer, Localized edge states in two-dimensional topological insulators: Ultrathin Bi films, Phys. Rev. B 83, 121310(R) (2011).
  * Liu  _et al._ [2011] Z. Liu, C.-X. Liu, Y.-S. Wu, W.-H. Duan, F. Liu, and J. Wu, Stable nontrivial ${Z}_{2}$ topology in ultrathin Bi (111) films: A first-principles study, Phys. Rev. Lett. 107, 136805 (2011).
  * Reis  _et al._ [2017] F. Reis, G. Li, L. Dudy, M. Bauernfeind, S. Glass, W. Hanke, R. Thomale, J. Schäfer, and R. Claessen, Bismuthene on a SiC substrate: A candidate for a high-temperature quantum spin hall material, Science 357, 287 (2017).
  * Chen  _et al._ [2013] L. Chen, Z. F. Wang, and F. Liu, Robustness of two-dimensional topological insulator states in bilayer bismuth against strain and electrical field, Phys. Rev. B 87, 235420 (2013).
  * Liu  _et al._ [2017] M.-Y. Liu, Y. Huang, Q.-Y. Chen, Z.-Y. Li, C. Cao, and Y. He, Strain and electric field tunable electronic structure of buckled bismuthene, RSC Adv. 7, 39546 (2017).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * Krukau  _et al._ [2006] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, J. Chem. Phys. 125, 224106 (2006).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * Aryasetiawan and Gunnarsson [1998] F. Aryasetiawan and O. Gunnarsson, The gw method, Rep. Prog. Phys. 61, 237 (1998).
  * Anisimov and Gunnarsson [1991] V. I. Anisimov and O. Gunnarsson, Density-functional calculation of effective coulomb interactions in metals, Phys. Rev. B 43, 7570 (1991).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong interactions: Orbital ordering in mott-hubbard insulators, Phys. Rev. B 52, R5467 (1995).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
  * Gebauer and Baroni [2000] R. Gebauer and S. Baroni, Magnons in real materials from density-functional theory, Phys. Rev. B 61, R6459 (2000).
  * He  _et al._ [2021] X. He, N. Helbig, M. J. Verstraete, and E. Bousquet, Tb2J: A python package for computing magnetic interaction parameters, Comput. Phys. Commun. 264, 107938 (2021).
  * Salpeter and Bethe [1951] E. E. Salpeter and H. A. Bethe, A relativistic equation for bound-state problems, Phys. Rev. 84, 1232 (1951).
  * Rohlfing and Louie [1998] M. Rohlfing and S. G. Louie, Electron-hole excitations in semiconductors and insulators, Phys. Rev. Lett. 81, 2312 (1998).