Brillouin zone integration

Table of Contents

• Documentation
  • Staggered k mesh
  • The Single Particle Sum
• Choices in Numerical Integration Methods
  • Tetrahedron Integration
    • Band crossing error
  • Sampling
    • Methfessel-Paxton sampling
• Comparing Sampling to Tetrahedron Integration
  • Response Functions
  • The Self Consistency Cycle
• References

This page documents how the Questaal suite performs numerical integrations over the Brillouin zone.

Documentation

Numerical integrations over the Brillouin zone (BZ) occur in several contexts, such as the sum of single particle energies. Care must be taken when performing this integration. In practice the first BZ is divided into a uniform mesh of $(n_{1}, n_{2}, n_{3})$ k-points parallel to the three primitive reciprocal lattice vectors, $\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}$.

Note: Questaal works with dimensionless reciprocal lattice vectors. $\mathbf{Q}_{1}, \mathbf{Q}_{2}, \mathbf{Q}_{3}$, with $\mathbf{G}_{i}{=}\mathbf{Q}_{i}{\times}2\pi/a$. ($a$ is the lattice constant in atomic units, read as ALAT in the input file).

The definition of the $\mathbf Q$ vectors is standard: construct a 3×3 matrix $P$ from the three (dimensionless) primitive lattice vectors $\mathbf{P}_{1}, \mathbf{P}_{2}, \mathbf{P}_{3}$ (the actual lattice is $a{\times}P$). Similarly let $Q$ be the three vectors composing the 3×3 matrix made from $\mathbf{Q}_{1}, \mathbf{Q}_{2}, \mathbf{Q}_{3}$. Then ${Q=(P^{-1})^\dagger}$ and $G=Q{\times}a/2\pi$.

Dividing $\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}$ respectively into $n_{1}, n_{2}, n_{3}$ sections creates a set of $(n_{1} {\times} n_{2} {\times} n_{3})$ microcells all of the same shape and size. The BZ integration is performed numerically by summing over each of the microcells. You specify $(n_{1}, n_{2}, n_{3})$ through BZ_NKABC. Thus the total number of $k$ points in the BZ is $n_{1}n_{2}n_{3}$. If symmetry operations apply, this number may be reduced.

A uniform mesh is the optimal general purpose quadrature for analytic, periodic functions, in the following sense. Consider the 1D case; the 3D case is exactly analogous. A numerical quadrature from a set of $n$ uniformly spaced points integrates exactly any function representable by a Fourier series of $n$ Fourier components or fewer, i.e. $$f(x) = \sum_{i=0:n-1} C_{i}\cos(2\pi ix/L) + S_{i}\sin(2\pi ix/L) \qquad (1)$$

for a function periodic in $L$.

Thus when integrating analytic functions, such as the single particle sum for the energy bands of an insulator, the uniform mesh is a good general purpose quadrature. In order to return the discussion to the three dimensional BZ, make the substitution $i \to G$ where $\mathbf G$ is a reciprocal lattice vector (some integer multiple of $\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}$). The band energy can be represented in a Fourier series $$E(k)=\sum_{T} C_{T}e^{2\pi i k \cdot T} \qquad (2)$$

$\mathbf k$ is the analog of $x$; it is periodic in the $\mathbf G$ vectors. Provided that the Fourier coefficients $C_T$ decay exponentially with $\mathbf{\|G\|}$, integration with a uniform quadrature converges exponentially with the number of $\mathbf k$ points. $C_{T}$ does more or less decay this way, typically. In one-band tight-binding theory $C_{T}$ would correspond to a hopping matrix element separated by distance $\mathbf{1/\|G\|}$.

Staggered k mesh

These expressions apply to any uniform mesh, but there is no requirement that the mesh include the origin, $\mathbf{k}{=}0$. In other words the mesh can be translated by an arbitrary amount in each of the three directions. By default the mesh is selected to include $\mathbf{k}{=}0$, but there is another natural choice, that applies to each of the three independent vectors $\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}$: a mesh that contain a point at the origin, or one offset so that points straddle the origin so one point midway between two mesh points falls at the origin. You specify this choice through BZ_BZJOB. Note that BZ_BZJOB is a tag for three numbers, one for each of the primitive reciprocal lattice vectors.

Under what circumstances does this matter? The default and staggered mesh generate different errors, but the errors are of the same order. There no a priori criterion to that will determine in advance which mesh will have the smallest integration error. For band calculations the main reason to choose a staggered mesh is that it may yield fewer irreducible $k$ points, even while the total number of points is unchanged. It also may yield more irreducible points. If computation time is an issue, you can check which yields the fewer number of points.

In the somewhat different context of linear response calculations, the response function as $\mathbf{q}{\to}0$ is difficult to treat. The GW code has an option to supply a staggered $k$ for calculating $\epsilon(r,r',\omega)$. Sometimes an offset mesh can significantly improve matters.

One striking example of this is the recently discovered solar cell material, $\mathsf{NH_4PbI_3}$. The figure above shows the static dielectric constant (more precisely $\epsilon_\infty$), calculated in the Quasiparticle Self-Consistent GW approximation as a function of n (n = n₁ = n₂ = n₃). The abscissa is scaled so that points are uniformly spaced in 1/n with the origin corresponding to n→∞. Calculations with both the normal mesh (which includes the Γ point) and an offset mesh are shown. Extrapolating either curve to n→∞ results in an estimate for ε_∞ of 4.3. Note how much more rapidly $\epsilon_\infty$ converges with the number of k divisions using a staggered mesh.

)

See also the description of this figure in the GWinput file documentation.

Note: Meshes generated by BZ_NKABC and BZ_BZJOB are identical to meshes generated by the Monkhorst Pack scheme, with suitably chosen parameters. But Questaal’s specification is simpler and more transparent.

The Single Particle Sum

Numerical integration over the BZ is necessary for band calculations, to obtain charge densities and the sum of single particle energies for total energy. The single-particle sum is the energy-weighted integral of the density-of-states $D(E)$, $$E_{band} = \int_{-\infty}^{E_{F}} D(E) E dE = \int_{-\infty}^{\infty} f(E)ED(E)dE \qquad (3)$$

where $f(E)$ is the Fermi function. $D(E)$ itself is given by a sum over eigenstates $i$ as $D(E) = \sum_{i}\delta(E-E_{i})$. In a crystal with Bloch states, $\mathbf k$ is a good quantum number we can write it as an integral over the BZ for each band: $$D(E) = \frac{V}{(2\pi)^{3}}\sum_{n} \int_{BZ} \delta \left(E_{n}(k) - E \right)d^{3} \mathbf{k} \qquad (4)$$

When the material is an insulator, sum over occupied states means sum over all the filled bands. Then $E_{band}$ can be written simply as $$E_{band} = \frac{V}{(2\pi)^{3}}\sum_n^{occ} \int_{BZ} E_{n}(k)d^{3}k \qquad (5)$$

Each band is an analytic, periodic function of $\mathbf k$, so numerical integration on the uniform mesh enables the band sum to converge with the $\mathbf k$ mesh in proportion to how quickly the $C_{T}$ of the Fourier series above decay.

For a metal the situation is more complicated. There is an abrupt truncation of the occupied part of the band at the Fermi level $E_{F}$. At 0K the behavior is non-analytic; at normal finite temperatures the non-analyticity is only a little softened.

Choices in Numerical Integration Methods

The band codes have three integration schemes that determine integration weights for every Quasi-Particle (QP) level at each $\mathbf k$.

1. Tetrahedron integration, enhanced by Peter Blöchl’s improved treatment of the Fermi surface [2]. It divides space into tetrahedra (there are 6 tetrahedra for each parallelepiped) and interpolates the energy linearly between the points of the tetrahedron. The Dirac δ-function $\delta [E_{n}(\mathbf{k}) - E]$ can be integrated over the tetrahedron by expressing it in the form $1/[\delta (E_{n}(\mathbf{k}) - E + i0]$ and approximating $E_{n}(k)$ by a linear function of k in the tetrahedron.
   To use the tetrahedron integration scheme, set BZ_TETRA=t.

2. Gaussian-broadened sampling integration, as extended by Methfessel and Paxton [1]. The basic idea of Gaussian broadening is that a quasiparticle level, which is a $\delta$-function of the energy in band theory, is smeared out into a Gaussian function. This makes properties analytic again. Methfessel and Paxton extended definition of the $\delta$-function from the traditional choice of a simple Gaussian of width $w$, to the product a Gaussian and a Hermite polynomial of order 2N. The resulting form of $\delta$-function has a richer shape: it is still peaked at the quasiparticle level and decays as a Gaussian, but it oscillates as it decays. The generalized Gaussian can rather dramatically improve convergence, as we show below. To use this method, turn off the tetrahedron method (BZ_TETRA=0) and specify $N$ and $w$ with BZ_N and BZ_W.

3. Fermi-Dirac-broadened sampling integration. This is similar to Gaussian broadening, but uses the (derivative of) the Fermi distribution for gaussian. It is the physically correct description for finite temperature, but it converges slowly with the k mesh at normal temperatures. To obtain rapidconvergence, artificially high temperatures are often used. To employ FD sampling, turn off the tetrahedron method (BZ_TETRA=0) and set BZ_N=−1. This tells the integrator to broaden the QP level with the (derivative of a) Fermi Dirac function instead of a Gaussian function. Specify the electron temperature $k_{B}T$ through BZ_W (in Ry). For a test case that illustrates Fermi-Dirac sampling, try

    fp/test/test.fp copt
   

   There is also a tutorial that uses this method to calculate the specific heat in Al.

Tetrahedron Integration

For self-consistent LDA calculations, Blöchl-enhanced tetrahedron integration is usually the method of choice for small metallic systems. The Figure on the left shows the power of the tetrahedron integration, especially when enhanced with the Blöchl weights. It shows the integration error $\Delta E$ in the Harris-Foulkes total energy for a fixed density (equivalent to the convergence of the single-particle sum) in MnBi, as a function of the number of $k$ divisions. The $x$ axis is scaled so that points are uniformly spaced in $1/n^2$ with the origin corresponding to $n{=}\infty$. MnBi forms in the NiAs structure, which is hexagonal with $c/a$ slightly less than 3/2. Accordingly, the number of divisions along the three directions was chosen to be $(n, n, 2 n/3)$, for an approximately uniform mesh spacing in the three directions. With Blöchl corrections, the single particle sum is converged to better than 2 μRy at 48 divisions. Without the Blöchl correction term this term convergence is only 100 μRy at 60 divisions. The integration error from the regular tetrahedron method scales as $1/n^2$; the Blöchl correction knocks out the leading $1/n^2$ term, as is evident from the Figure. For a simple band the scaling improves to $1/n^3$.

Band crossing error

There are caveats to keep in mind with respect to this method. Sometimes a pair of bands can cross near $E_{F}$. The tetrahedron integrator, which must interpolate a single band between the four corners of a tetrahedron, expects the band to behave in an analytic way. If two bands cross inside the tetrahedron, both the lower band and upper band are non-analytic (in reality both are analytic but the integrator isn’t smart enough to detect that there is a band crossing and switch the levels). When this occurs usually the integrator yields a non-integral number of electrons. It will detect this and print out a warning message:

 (warning): non-integral number of electrons --- possible band crossing at E_f

When you get this message usually (not always!) it is because of a band crossing problem. As you proceed to self-consistency it may go away. If not, you need to modify your $\mathbf k$ mesh. The larger the system with a denser mesh of bands, the more serious this issue becomes.

Another drawback with the tetrahedron scheme is that all the eigenvalues must be known before the weight for any one can be determined. This is because integrals are done over tetrahedra, each of which involves four $\mathbf k$ points. Thus unless some trickery is used, two band passes are required: one pass to get all the eigenvalues to determine integration weights, and a second pass to accumulate quantities that depend on eigenfunctions, such as the charge density. The FP code has a “mixed” scheme available: the charge density is calculated by sampling while the band sum is calculated by the tetrahedron method. This avoids a double band pass, and is effective because, thanks to the variational principle, density can be evaluated to less precision than $E_{band}$.

Sampling

As the mesh spacing becomes infinitely fine $(n{\to}\infty)$, the quadrature estimate for integrated property of interest (e.g. $E_{band}$) becomes independent of $n$. But with sampling the converged result at the $1/n \to 0$ limit is not exact, but depends on $w$: the exact (or 0K) result is obtained only for both $(n{\to}\infty)$ and $w{\to}0$. This should be immediately obvious on physical grounds: the band energy should depend e.g. on temperature.

Thus with sampling there are two independent issues to contend with: how rapidly the result converges with $n$, and the error in the converged result. These issues are not entirely independent: how rapidly the quadrature convergences with $n$ also depends on w: the larger w, the more rapid the convergence. This is also easy to understand: the greater the broadening, the smoother and more analytic the behavior at $E_{F}$.

The figure above illustrates these points. The BZ integration in the Harris-Foulkes energy is shown for MnBi — the same system used to demonstrate tetrahedron integration. In this figure integrations are carried out with Methfessel-Paxton (MP) sampling described below using $N$=1 for various widths $w$. When $w$ is 0.015 (rather large), the energy has converged to within 1μRy of the n→∞ limit for n≥42. But, it converges to a result differing from the exact one by 40μRy. As $w{\rightarrow}0$, ΔE should vanish with n→∞. Indeed this is the case: if $w$=2 mRy the integral is converged to well within 1μRy. But the approach to convergence is much slower: only for n≥90 does the integral stabilize to within 1μRy. The $w$=5 and 10 mRy cases show behavior intermediate between $w$=2 and 15 mRy. The converged integration errors are 2.6μRy and 3.5μRy, respectively.

Methfessel-Paxton sampling

If the $\delta$-function is broadened with a Gaussian in the standard manner, good results are obtained only when $w$ is very small, and the situation would seem to be hopeless. What saves the day is Methfessel and Paxton’s [1] generalization of the broadening function. It starts from the representation of the exact $\delta$-function as bilinear combinations of Hermite polynomials and a Gaussian $$\delta(x) = \lim_{N \to \infty} D_{N}(x), D_{N}(x) = \sum_{m=0}^{N} A_{m}H_{2m}(x)e^{-x^{2}} \qquad (6)$$

Standard Gaussian broadening represents $\delta(x)$ by the $m{=}0$ term alone. But by truncating the series at some low order $N$ larger than 0, better representations of the $\delta$-function are possible (figure above) and the integration quality can be dramatically improved. $E_{band}$ can be well converged without requiring an excessively small $w$.

Choosing any $N$>0 can dramatically improve the quality of sampling integration. The figure above shows how |$\Delta E$| for MnBi depends on $N$, for a relatively coarse $k$ mesh $(n{=}12)$, a very fine one $(n{=}90)$, and intermediate one ($(n{=}42)$ some intermediate ones. For the coarse mesh no advantage is found by varying $N$. This is because the $k$ mesh integration error (i.e. the number of $C_{T}$ in Eq. (2) integrated exactly by the mesh is finite, and essentially equal to $n_{1}n_{2}n_{3}$) is larger than the error originating from approximations to the $\delta$-function.

For the intermediate mesh (n=42) the situation is a bit different. If $N$=0, the approximation to the $\delta$-function dominates $\Delta E$: increasing $N$ from 0→1 reduces it by a factor of almost 100! But increasing $N$ above 1 has no additional effect: the error remains fixed in the 3-4μRy range. It has already become dominated by the fineness of the $k$-mesh, any adequate ($N$≥1) representation of the $\delta$-function will serve equally well. Not insignificantly, $\Delta E$ from the tetrahedron method is in the same range (4μRy).

For the very fine mesh (n=90) much the same story appears, only now that the transition from $\Delta E$ being dominated by $\delta$-function error to the discreteness of the $k$ mesh doesn’t occur until $N$=4. Because the mesh is so fine, $\Delta E$ is tiny, of order 0.1μRy.

MP makes it possible to perform quality integrations with sampling. The figure below shows the MnBi test case, now fixing w at 0.015 and varying $N$. ($N$=0 is not shown: $\Delta E$ is a couple orders of magnitude larger for this w.) In the $n\to\infty$ limit, increasing the polynomial order $N$ has an effect somewhat analogous to decreasing $w$. The discrepancy with the exact (0K) result decreases, but the quadrature converges more slowly with n.

There is trade-off: another parameter is available to tinker with. $\Delta E$ is now a function of an additional variable: $\Delta E = \Delta E(n_{1},n_{2},n_{3},N,w)$. The optimum combination of $(n_{1},n_{2},n_{3},N,w)$ depends on the precision you need and computer time available. While there is no universal prescription that can be worked out in advance, the following guidelines seem to work fairly well:

1. Choose the ratios $n_{1}:n_{2}:n_{3}$ as closely as possible to the ratios $G_{1}:G_{2}:G_{3}$, where the latter are the lengths primitive reciprocal lattice vectors.
2. Use at least $N$=1 and unless you are looking for energy resolution (see next section) preferably $N$=3 or 4 with $w$~0.01. There doesn’t seem to be much downside to using a sizable $N$, but there can be a significant gain.
3. If you want very high quality integration, use a fine mesh with $w$~0.005 and $N$ around 3 or 4.

Comparing Sampling to Tetrahedron Integration

Which method should you choose? The answer depends strongly on context.

Response Functions

Consider first the calculation of quantities that involve transitions between two states, such as the joint density-of-states or dielectric function. The Figure below compares Im ε(ω) calculated for Ag, tetrahedron to MP sampling ($N$=1, $w$=0.02). Calculations for 24 and 48 $k$ divisions are shown. The sharp shoulder at 0.2 Ry marks the onset of transitions from the Ag d states to conduction band states. This shoulder is seen in experiment, at a slightly higher energy; this is because the LDA puts the Ag d states too shallow.

In the tetrahedron case, doubling the mesh affects Im ε(ω) only slightly, showing that it is already well converged at 24 divisions. Sampling with 24 divisions shows Im ε(ω) oscillating around the converged result. With 48 divisions, the oscillations are much weaker, but you can still see Im ε(ω) dip below 0 just before the onset of the shoulder, and round off the first peak around ω=0.28 Ry.

It is possible to get good definition in Im ε(ω) via MP sampling, by increasing N or reducing W, or some judicious combination of the two. However, many more k-points are required. The figure below shows that Im ε(ω), calculated with $N$=1, $w$=.01, can pick up the shoulder at ω=0.28 Ry, but there is ringing around the converged result. This is because $D_{N}(x)$ oscillates around 0 for $N$≥1, as seen in the figure for D_N(sampling section). Increasing n to 96 divisions yields Im ε(ω) with definition almost as good as the tetrahedron method.

In general, the larger $N$ you choose, the better definition you will have if your k mesh is sufficiently fine. But the larger $N$, the more ‘wiggly’ response function will be when you haven’t converged the $k$ mesh sufficiently.

Generally speaking, for linear response calculations the tetrahedron method is superior. But, see Note, next.

Note: the existing implementations of tetrahedron may have a problem for tetrahedra that contain multiple bands crossing $E_{F}$. We have not fully understood the problem yet, but have discovered that for some metals, (e.g. a four-atom supercell of Ag), a spurious peak in Im ε(ω) can appear for ω→0. A peak actually should exist there, but it originates from the Drude term which is not properly accounted for in this method. So if the peak appears, it does so for spurious reasons.

The Self Consistency Cycle

BZ integration is also required to make the ordinary density-of-states (either total, or resolved by orbital) for the self-consistency cycle or to calculate the single-particle sum for total energy. Here also tetrahedron integration is generally superior, and is usually preferred unless problems arise, e.g. issues with band crossings.

Nevertheless there are circumstances where sampling is a better choice, even apart from the band crossing issue. To calculate the shear modulus of Al, tetrahedron integration requires a very fine $k$ mesh. But doing integrations with Fermi Dirac statistics at a rather high temperature, rapidly convergent results can be obtained. To check convergence you can use several temperatures and extrapolate to 0K.

Another example is the magnetic anisotropy. As with the Al shear constant, it is not the absolute value of the energy (mostly single particle sum) that matters, but its variation with some perturbation such as a shear or the addition $L\cdot S$ coupling to the hamiltonian. In such cases the error of the energy difference can be smaller, and more stable, with a suitably chosen sampling method. The table below shows the magnetic anisotropy, in Ry, for CoPt in the L1₀ structure, using a $k$ mesh of 30×30×24 divisions. (A check was made with 40×40×32 divisions using sampling, and essentially the same results were found.) Calculations were carried out in the LDA with both tetrahedron and sampling methods, either by evaluating the change in the Harris-Foulkes total energy, or by calculating the energy through coupling constant integration of $L \cdot S$. Calculations were both of the non-self-consistent variety (frozen density at some suitable starting point) and also self-consistent. While all the numbers are fairly similar, there is more variation in the tetrahedron results. Note especially the ehf and cc-int columns. They should give the same answer in all cases.

                                                 ehf(tet)   cc-int(tet)  ehf(samp)  cc-int(samp)
non self-consistent, rho from L.S=0             -0.000084   -0.000077   -0.000078   -0.000077
non self-consistent, rho from full L.S || z     -0.000076   -0.000069   -0.000075   -0.000075
self consistent with 8 symmetry operations      -0.000077   -0.000070   -0.000078   -0.000078
self consistent with no symmetry                -0.000078               -0.000077   -0.000076

References

[1] M. Methfessel and A.T. Paxton, Phys. Rev. B, 40, 3616 (1989)

[2] P.E. Blöchl, O. Jepsen and O.K. Andersen, Phys. Rev. B 49, 16223 (1994)