Brillouin zone integration

Table of Contents


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)k(n_{1}, n_{2}, n_{3}) \mathbf k-points parallel to the three primitive reciprocal lattice vectors, G1,G2,G3G_{1}, G_{2}, 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,P3P_{1}, P_{2}, P_{3}. G1,G2,G3G_{1}, G_{2}, G_{3} are the three vectors composing the 3×3 matrix G\mathbf G.

Dividing G1,G2,G3G_{1}, G_{2}, 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 Z 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)f(x)=\sum_{i=0:n-1} C_{i}cos(2\pi ix/L) + S_{i}sin(2\pi ix/L)

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,G3G_{1}, G_{2}, G_{3}). The band energy can be represented in a Fourier series

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

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 a separated by distance 1/G\mathbf{1/\|G\|}. For each of the G1,G2,G3G_{1}, G_{2}, G_{3} directions independently, there are two ways to choose a uniform mesh: it may contain a point that coincides with the origin, or it is offset so that the points straddle the origin. You specify this choice through BZ_BZJOB. Under what circumstances does this matter? For band calculations the main reason is that sometimes a staggered mesh will yield fewer irreducible kk points, even while the total number of points is unchanged. 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,ω)\eta(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}.

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

The Singe 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)dEE_{band} = \int_{-\infty}^{E_{F}} D(E) E dE = \int_{-\infty}^{\infty} f(E)ED(E)dE

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)d3kD(E) = \frac{V}{(2\pi)^{3}}\sum_{n} \int_{BZ} \delta \left(E_{n}(k) - E \right)d^{3} \mathbf{k}

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π)3BZEn(k)d3kE_{band} = \frac{V}{(2\pi)^{3}}\sum_{BZ} E_{n}(k)d^{3}k

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 proprtion 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 QP level at each k\mathbf k.

  1. Tetrahedron integration, enhanced by Peter Bloechl’s improved treatment of the Fermi surface. It divides space into tetrahedra (there are 6 tetra 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. 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 the correct description for finite temperature. But for 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 Bloechl 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/n2 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 Bloechl corrections, the single particle sum is converged to better than 2 μRy at 48 divisions. Without the Bloechl correction term this term convergence is only 100 μRy at 60 divisions. The integration error from the regular tetrahedron method scales as 1/n21/n2; the Bloechl correction knocks out the leading 1/n21/n2 term, as is evident from the Figure. For a simple band the scaling improves to 1/n31/n3.

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, the 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 w: the exact (or 0K) result is obtained only for (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}.

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

δ(x)=limNDN(x),DN(x)=m=0NAmH2m(x)ex2\delta(x) = \lim_{N \to \infty} D_{N}(x), D_{N}(x) = \sum_{m=0}^{N} A_{m}H_{2m}(x)e^{-x^{2}}

Standard gaussian broadening represents δ(x)\delta(x) by the m=0m = 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, right) and the integration quality can be dramatically improved. EbandE_{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 on the left shows how |ΔE\Delta E| for MnBi depends on N, 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). For the coarse mesh no advantage is found by varying N. This is because the kk mesh integration error (i.e. the number of CTC_{T} in Eq.(1) 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 N=0, the approximation to the δ\delta-function dominates ΔE\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 kk-mesh, any adequate (N≥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 N=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 right figure shows the MnBi test case, now fixing w at 0.015 and varying N. (N=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 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 tradeoff: 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 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 sizeable N, but there can be a significant gain.
  3. If you want very hiqh 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 on the right compares Im ε(ω) calculated for Ag, tetrahedron to MP sampling (N=1, W=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 shallowly.

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 left hand figure 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 DN(x)D_{N}(x) oscillates around 0 for N≥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. 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 kk mesh sufficiently.

Generally speaking, for linear response calculations the tetrahedron method is superior. But, be advised:

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 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 30x30x24 divisions. (A check was made with 40x40x32 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’ column. 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


Edit This Page