Brillouin zone integration

Table of Contents

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


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 (n1,n2,n3)(n_{1}, n_{2}, n_{3}) k-points parallel to the three primitive reciprocal lattice vectors, G1,G2,G3\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}.

The definition of the G\mathbf G vectors is standard: define G=(P1)\mathbf{G=(P^{-1})^\dagger} where P\mathbf P is 3×3 matrix constructed from the three primitive lattice vectors P1,P2,P3\mathbf{P}_{1}, \mathbf{P}_{2}, \mathbf{P}_{3}. G1,G2,G3\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3} are the three vectors composing the 3×3 matrix G\mathbf G.

Dividing G1,G2,G3\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3} respectively into n1,n2,n3n_{1}, n_{2}, n_{3} sections creates a set of (n1×n2×n3)(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 (n1,n2,n3)(n_{1}, n_{2}, n_{3}) through BZ_NKABC. Thus the total number of kk points in the BZ is n1n2n3n_{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 nn uniformly spaced points integrates exactly any function representable by a Fourier series of nn Fourier components or fewer, i.e.

f(x)=i=0:n1Cicos(2πix/L)+Sisin(2πix/L)(1)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 LL.

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 iGi \to G where G\mathbf G is a reciprocal lattice vector (some integer multiple of G1,G2,G3\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3}). The band energy can be represented in a Fourier series

E(k)=TCTe2πikT(2)E(k)=\sum_{T} C_{T}e^{2\pi i k \cdot T} \qquad (2)

k\mathbf k is the analog of xx; it is periodic in the G\mathbf G vectors. Provided that the Fourier coefficients CTC_T decay exponentially with G\mathbf{\|G\|}, integration with a uniform quadrature converges exponentially with the number of k\mathbf k points. CTC_{T} does more or less decay this way, typically. In one-band tight-binding theory CTC_{T} would correspond to a hopping matrix element separated by distance 1/G\mathbf{1/\|G\|}. For each of the G1,G2,G3\mathbf{G}_{1}, \mathbf{G}_{2}, \mathbf{G}_{3} directions independently, there are two natural ways to choose a uniform mesh: a mesh that contain a point that coincides with the origin, or one offset so that points straddle the origin. You specify this choice through BZ_BZJOB. Under what circumstances does this matter? They generate different errors, but the errors are of the same order. There no a priori criterion to that will determine in advance the choice with the smallest integration error. For band calculations the main reason is to choose a staggered mesh is that it will occasionally yield fewer irreducible kk points, even while the total number of points is unchanged.

NH4 eps static

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

One striking example of this is the recently discovered solar cell material, NH4PbI3\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 = n1 = n2 = n3). 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.

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)D(E),

Eband=EFD(E)EdE=f(E)ED(E)dE(3)E_{band} = \int_{-\infty}^{E_{F}} D(E) E dE = \int_{-\infty}^{\infty} f(E)ED(E)dE \qquad (3)

where f(E)f(E) is the Fermi function. D(E)D(E) itself is given by a sum over eigenstates ii as D(E)=iδ(EEi)D(E) = \sum_{i}\delta(E-E_{i}). In a crystal with Bloch states, k\mathbf k is a good quantum number we can write it as an integral over the BZ for each band:

D(E)=V(2π)3nBZδ(En(k)E)d3k(4)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 EbandE_{band} can be written simply as

Eband=V(2π)3noccBZEn(k)d3k(5)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 k\mathbf k, so numerical integration on the uniform mesh enables the band sum to converge with the k\mathbf k mesh in proportion to how quickly the CTC_{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 EFE_{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 k\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. To use the tetrahedron integration scheme, set BZ_TETRA.

  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 ww, 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 NN and ww 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, converges slowly with the k mesh at normal temperatures. To obtain rapid convergence, 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 kBTk_{B}T through BZ_W (in Ry). For a test case that illustrates Fermi-Dirac sampling, try

     fp/test/test.fp copt

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 ΔE\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 kk divisions. The xx axis is scaled so that points are uniformly spaced in 1/n21/n^2 with the origin corresponding to n=n{=}\infty. MnBi forms in the NiAs structure, which is hexagonal with c/ac/a slightly less than 3/2. Accordingly, the number of divisions along the three directions was chosen to be (n,n,2n/3)(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/n21/n^2; the Blöchl correction knocks out the leading 1/n21/n^2 term, as is evident from the Figure. For a simple band the scaling improves to 1/n31/n^3.

MnBi etet

Band crossing error

There are caveats to keep in mind with respect to this method. Sometimes a pair of bands can cross near EFE_{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 k\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 k\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 EbandE_{band}.


As the mesh spacing becomes infinitely fine (n)(n{\to}\infty), the quadrature estimate for integrated property of interest (e.g. EbandE_{band}) becomes independent of nn. But with sampling the converged result at the 1/n01/n \to 0 limit is not exact, but depends on ww: the exact (or 0K) result is obtained only for both (n)(n{\to}\infty) and w0w{\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 nn, and the error in the converged result. These issues are not entirely independent: how rapidly the quadrature convergences with nn 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 EFE_{F}.

MnBi etet

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 NN=1 for various widths ww. When ww 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 w0w{\rightarrow}0, ΔE should vanish with n→∞. Indeed this is the case: if ww=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 ww=5 and 10 mRy cases show behavior intermediate between ww=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 ww 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

δ(x)=limNDN(x),DN(x)=m=0NAmH2m(x)ex2(6)\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)

MP delta

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

MnBi samp2

Choosing any NN>0 can dramatically improve the quality of sampling integration. The figure above shows how |ΔE\Delta E| for MnBi depends on NN, for a relatively coarse kk mesh (n=12)(n{=}12), a very fine one (n=90)(n{=}90), and intermediate one ((n=42)(n{=}42) some intermediate ones. For the coarse mesh no advantage is found by varying NN. This is because the kk mesh integration error (i.e. the number of CTC_{T} in Eq. (2) integrated exactly by the mesh is finite, and essentially equal to n1n2n3n_{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 NN=0, the approximation to the δ\delta-function dominates ΔE\Delta E: increasing NN from 0→1 reduces it by a factor of almost 100! But increasing NN 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 kk-mesh, any adequate (NN≥1) representation of the δ\delta-function will serve equally well. Not insignificantly, ΔE\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 ΔE\Delta E being dominated by δ\delta-function error to the discreteness of the kk mesh doesn’t occur until NN=4. Because the mesh is so fine, ΔE\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 NN. (NN=0 is not shown: ΔE\Delta E is a couple orders of magnitude larger for this w.) In the nn\to\infty limit, increasing the polynomial order NN has an effect somewhat analogous to decreasing ww. The discrepancy with the exact (0K) result decreases, but the quadrature converges more slowly with n.

MnBi samp3

There is trade-off: another parameter is available to tinker with. ΔE\Delta E is now a function of an additional variable: ΔE=ΔE(n1,n2,n3,N,w)\Delta E = \Delta E(n_{1},n_{2},n_{3},N,w). The optimum combination of (n1,n2,n3,N,w)(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 n1:n2:n3n_{1}:n_{2}:n_{3} as closely as possible to the ratios G1:G2:G3G_{1}:G_{2}:G_{3}, where the latter are the lengths primitive reciprocal lattice vectors.
  2. Use at least NN=1 and unless you are looking for energy resolution (see next section) preferably NN=3 or 4 with ww~0.01. There doesn’t seem to be much downside to using a sizable NN, but there can be a significant gain.
  3. If you want very high quality integration, use a fine mesh with ww~0.005 and NN 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 (NN=1, ww=0.02). Calculations for 24 and 48 kk 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.

ag eps

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 NN=1, ww=.01, can pick up the shoulder at ω=0.28 Ry, but there is ringing around the converged result. This is because DN(x)D_{N}(x) oscillates around 0 for NN≥1, as seen in the figure for DN(sampling section). Increasing n to 96 divisions yields Im ε(ω) with definition almost as good as the tetrahedron method.

ag eps

In general, the larger NN you choose, the better definition you will have if your k mesh is sufficiently fine. But the larger NN, the more ‘wiggly’ response function will be when you haven’t converged the kk 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 EFE_{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 kk 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 LSL\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 L10 structure, using a kk 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 LSL \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


[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)

Edit This Page