k-Point Grids and Convergence

Learning Objectives

The previous seven sections of this chapter built up reciprocal space — the reciprocal lattice (3.2), the Brillouin zone (3.3), high-symmetry paths (3.4), Bloch's theorem (3.5), the structure factor (3.6), and systematic absences (3.7). Now we put it to work. Every electronic-structure code, VASP included, replaces continuous Brillouin-zone integrals by discrete sums over a finite mesh of k-points. Choosing that mesh is the single most consequential numerical decision in any DFT calculation. By the end of this section you should be able to:

1. Explain why integration over the BZ is unavoidable, and what quantity is actually being averaged.
2. Translate the integral $\frac{V}{(2\pi)^3}\int_{\text{BZ}} f(\mathbf{k})\,d\mathbf{k}$ into the discrete sum $\sum_{\mathbf{k}} w_{\mathbf{k}}\,f(\mathbf{k})$ and know what the weights $w_{\mathbf{k}}$ mean.
3. State the inverse-volume rule: $(\text{real cell side}) \times (\text{k-mesh density}) \;=\; \text{constant}$. A 2×2×2 supercell needs half the k-mesh of the primitive cell per axis.
4. Recognise that convergence is not always monotonic — band crossings, Fermi surfaces, and grid-edge effects can make the energy oscillate before it settles.
5. Read VASP's KPOINTS file and understand what each line of the simplest automatic form actually requests.
6. Run a clean convergence test in Python or VASP and decide when the answer is "converged enough" for the question you are asking.

One-line preview: a k-point grid is a Riemann-sum discretisation of the Brillouin zone. The denser it is, the closer the sum gets to the true integral — but the right density depends on the smoothness of the integrand and the size of your real-space cell.

Why Do We Sample at All?

Bloch's theorem (Section 3.5) tells us that every single-particle eigenstate of a periodic crystal can be labelled by a wavevector $\mathbf{k}$ living in the first Brillouin zone and a band index $n$. The total energy of the electrons, the charge density, the density of states, the optical spectrum — every observable in this chapter and the next — is built from a Brillouin-zone average of band quantities. For total energy at zero temperature:

where $f(\varepsilon)$ is the Fermi–Dirac occupation (effectively a step function at $T = 0$ for insulators, smeared by temperature for metals). The factor $V_{\text{cell}}/(2\pi)^3$ converts the integral into a per-cell quantity using the result we proved in Section 3.2 — the reciprocal-cell volume is $(2\pi)^3/V_{\text{cell}}$.

For an infinite crystal this is a clean continuous integral. For a computer it is something we must approximate, because the integrand $\varepsilon_{n\mathbf{k}}$ is the eigenvalue of the Kohn–Sham Hamiltonian at $\mathbf{k}$ — every evaluation costs an entire diagonalisation. We pick a finite set of $\mathbf{k}$-points, diagonalise once at each, and add the results with weights.

From an Integral to a Weighted Sum

Replace the integral by a Riemann-style sum:

Three ideas are bundled into that line.

1. The mesh. A discrete subset of points in the BZ. The most popular choice (Section 3.9) is a regular $N_1\times N_2\times N_3$ grid — Monkhorst–Pack — but other meshes (Γ-only, Chadi–Cohen, hybrid) are common.
2. The weights. $w_{\mathbf{k}} = (V_{\text{BZ}}/N_{\text{tot}}) / V_{\text{BZ}} = 1/N_{\text{tot}}$ for an unfolded uniform grid. After symmetry reduction (Section 3.4 taught us what symmetries do to k), points related by a symmetry operation are merged into one representative whose weight equals its orbit size. So in a VASP IBZKPT file, weights of 1, 4, 6, 8, … encode "how many full-BZ partners does this irreducible point stand in for?".
3. The convergence claim. If $f(\mathbf{k})$ is smooth (band-gap material, insulator) the sum approaches the integral with error that falls off rapidly in the mesh density. If $f(\mathbf{k})$ has a kink — a Fermi surface, a Dirac cone — convergence is much slower and requires either denser meshes or special tricks (tetrahedron method, Gaussian smearing).

Read the discrete formula as "average f over the mesh, weighting each point by how much of the BZ it stands for". The mesh is your ruler; the weights are the tick spacing.

The Inverse-Volume Rule

How dense should the mesh be? Section 3.2 already gave the answer in compressed form: $V_{\text{BZ}} = (2\pi)^3 / V_{\text{cell}}$. Spelt out for a uniform grid:

The thing that matters for accuracy is $\Delta k_i$, not $N_i$. A 2×2×2 supercell has each side $a_i$ twice as long, so its reciprocal vector $|\mathbf{b}_i|$ is twice as short, and a mesh with $N_i = 4$ on the supercell achieves the same $\Delta k_i$ as $N_i = 8$ on the primitive cell. Doubling the cell halves the required N.

VASP exposes this directly: the KSPACING tag in INCAR sets $\Delta k$ in Å⁻¹; VASP then picks the smallest $N_i$ that satisfies it. A typical metallic calculation uses KSPACING = 0.2 (Å⁻¹); a typical insulator uses 0.3–0.4 Å⁻¹. We will see why in a moment.

Interactive — k-Grid in a 2D Brillouin Zone

Slide the controls below to change the grid density $N$ on a 2-D square Brillouin zone (a tight-binding analogue of a Cu (001) surface, with $\varepsilon(\mathbf{k}) = -2t\,(\cos k_x a + \cos k_y a)$). The colour map is the band; cyan dots are filled k-points ($\varepsilon < 0$), pink dots are empty. The right panel reports the current sum estimate, the "exact" value computed on a 256×256 reference grid, and the absolute error. The log–log convergence curve below shows how the error falls as you slide $N$ upward.

Three things to play with:

• Slide N up. Watch the error fall. Notice it does not fall monotonically — it overshoots and undershoots. The orange ring on the convergence plot is your current $N$.
• Toggle Γ-centred shift. Without the half-step, the unshifted grid hits the BZ edges (and the Fermi-surface line at $k_x = \pm\pi/a$) directly. Edge points are ambiguous: do they belong to this BZ or the next? Convergence slows down. The shift dodges this by sampling only the interior of each mesh cell.
• Compare even and odd N. Near band-crossing lines, an even-$N$ grid that straddles the crossing can be lucky or unlucky. The oscillation in the curve is real — it is the metallic Fermi-surface signature of the integrand having a kink at $\varepsilon = 0$.

How Convergence Actually Behaves

Three regimes show up in practice. Knowing which one you are in is half the battle.

1. Insulators and large-gap semiconductors

Smooth integrand, no occupation discontinuity. Errors fall like a polynomial — for a 3-D Monkhorst–Pack grid, like $1/N^4$ for many ground-state quantities. A coarse 4×4×4 mesh on the primitive Si cell (8 atoms in the conventional cube) is already excellent; 6×6×6 is overkill for most properties.

2. Metals

The integrand has a step at the Fermi surface; the analogue of taking the Riemann sum of a step function. Bare uniform grids converge as slowly as $1/N$. The cure is smearing — replace the abrupt step $f(\varepsilon)$ by a smoothed version, recover smoothness, and pay a tiny accuracy tax. VASP supports several flavours:

• ISMEAR = -5 — tetrahedron with Blöchl correction. Best for static properties of metals; needs a Γ-centred mesh and at least 4 k-points in the irreducible BZ.
• ISMEAR = 1, 2 — Methfessel–Paxton of order 1 or 2. Excellent for total energies and forces in metals at modest mesh density.
• ISMEAR = 0, SIGMA = 0.05 — Gaussian smearing. Robust default for ionic relaxations of metals.

3. Narrow-gap, doped, or magnetic systems

Treat as metal-like for k-point convergence. Doped semiconductors have a Fermi level inside a band; magnetic systems with spin polarisation have Fermi surfaces in each spin channel. Both want the dense mesh + smearing combo, not the sparse-mesh insulator workflow.

Choosing a Smearing — Gaussian, MP, Tetrahedron

We just hand-waved the metal cure as "smearing". That word hides four distinct numerical schemes — VASP's ISMEAR tag — that look superficially similar but differ in which order in σ they cancel. Picking the wrong one is one of the most common mistakes in plane-wave DFT. This block lays the four side by side, derives the why for each, and ends with a one-glance decision rule.

The integral we are approximating

For a metal, the total electron number per cell is

with $\Theta$ the step function. The same kind of integral, with extra factors, gives the total energy and the density. The integrand jumps at the Fermi surface, so a finite Riemann sum on a uniform mesh has unavoidable noise unless we replace $\Theta$ by a smooth approximation $f_\sigma(\mu - \varepsilon)$ of width $\sigma$. The four VASP options are four different choices of $f_\sigma$ — and a fifth, the tetrahedron method, throws the smooth-occupation idea away entirely and integrates a piecewise-linear interpolation of the bands instead.

The four kernels — what each one is

The kernels, plotted side by side

Toggle each scheme on or off and slide $\sigma$. Three things should be visible at a glance:

• All curves cross $f = \tfrac{1}{2}$ at the Fermi level — that is the definition of $\mu$.
• As $\sigma \to 0$, every kernel converges to the grey step. The rate of convergence is what differs.
• Methfessel–Paxton overshoots 1 just below $\mu$ and undershoots 0 just above. The area of those two wiggles is exactly equal — and that cancellation is what kills the order-1 error in $\sigma$. Gaussian and Fermi–Dirac do not have this cancellation.

Why MP's wiggles cancel — one paragraph

For a smooth band structure, the total energy can be Taylor-expanded around its $\sigma = 0$ value. Methfessel and Paxton (PRB 40, 3616) chose the Hermite-correction coefficients precisely so that the coefficients of $\sigma, \sigma^2, \ldots, \sigma^{2N}$ in that expansion are zero. The first surviving error is $\mathcal{O}(\sigma^{2N+1})$. So N = 1 means the answer is correct through $\sigma^2$, and you can run with $\sigma \sim 0.2$ eV without paying for it — provided the mesh is dense enough that the surviving piecewise wobble averages out. That last clause is why MP requires a denser k-mesh than Gaussian: under-sampled, the wiggles do not cancel and you can get a wrong sign on the force.

Why tetrahedron is exact at zero σ

The tetrahedron method (ISMEAR = −5) does not smear at all. It treats the Brillouin zone as a tiling of tiny tetrahedra (8 per Monkhorst–Pack cube) with one band eigenvalue at each vertex, linearly interpolates the bands inside each tetrahedron, and then integrates the interpolation analytically. The Blöchl correction adjusts the integration weights to remove the leading curvature error, lifting accuracy from $\mathcal{O}(\Delta k^2)$ to $\mathcal{O}(\Delta k^4)$. The result is the most accurate static-energy / static-DOS method on a fixed k-mesh — but it depends on neighbouring k-points sharing a connected tetrahedron, which fails for sparse meshes (hence the "at least 4 irreducible k-points" rule), Γ-only sampling, and non-Γ-centred MP grids.

Anisotropic Cells Need Anisotropic Grids

The inverse-volume rule applies per axis. A wurtzite GaN cell with $a = 3.19$ Å and $c = 5.19$ Å has $|\mathbf{b}_3|$ shorter than $|\mathbf{b}_1|$ by the ratio $a/c \approx 0.61$. To match $\Delta k$ on all three axes you want an $N_1 \times N_1 \times N_3$ grid with $N_3 \approx 0.61\,N_1$. A common pragmatic choice: 8×8×5 instead of 8×8×8.

Slab calculations push this even further. A slab cell with vacuum along $\mathbf{a}_3$ often has $c \sim 25$ Å of vacuum, making $|\mathbf{b}_3|$ tiny. Sampling along that axis is wasteful — most authors set $N_3 = 1$ for slabs (a single Γ-line of k-points) and use the in-plane $N_1 \times N_2$ as for the bulk surface unit cell.

Symmetry — The Irreducible Wedge

Section 3.4 noticed that the cubic BZ has 48 symmetry operations; the hexagonal BZ has 24. Apply these to the full mesh and most points fall on top of each other. VASP automatically reduces the full $N_1 N_2 N_3$ grid to the irreducible Brillouin zone (IBZ): the smallest wedge from which the rest of the BZ can be regenerated by symmetry.

For an FCC primitive cell with an 8×8×8 mesh:

• Full BZ: 512 points.
• Irreducible BZ for the 48-element $O_h$ symmetry: 60 points (printed at the top of IBZKPT).
• Cost saving: ~8.5×. The diagonalisation runs only over the irreducible set; the BZ averages are reconstructed using the orbit-size weights $w_{\mathbf{k}}$.

Two everyday consequences:

1. If you break the symmetry — say, by introducing a single Mn dopant in a CdSe supercell — the irreducible set grows. A 4×4×4 mesh that gave you 8 irreducible k-points in pristine CdSe might give 32 in the doped supercell. Plan your CPU budget accordingly.
2. A Γ-centred grid keeps every point group operation as a strict symmetry of the mesh. A non-Γ-centred (offset) grid can break inversion symmetry and double the irreducible set. For magnetic and spin-orbit calculations, Γ-centred is almost always the right choice.

VASP — The KPOINTS File at a Glance

VASP's simplest k-point input — and the one you will write most often — is the automatic form. Five lines, no surprises:

1Automatic mesh                         ! line 1: comment, free-text
20                                       ! line 2: 0 → automatic generation
3Gamma                                    ! line 3: Gamma  (Γ-centred) or Monkhorst-Pack
48 8 8                                    ! line 4: subdivisions N1 N2 N3
50 0 0                                    ! line 5: optional shift (s1 s2 s3)

1. Line 1. A free-form comment. Use it to document the mesh density (e.g. KSPACING ≈ 0.20 Å⁻¹) so your future self can remember what was tested.
2. Line 2. Number of k-points to read explicitly. 0 means "generate automatically from the next three lines". A positive integer would be followed by a manual list — used for band-structure paths (Section 3.10).
3. Line 3. Mesh type. Gamma places one point exactly at $\mathbf{k} = 0$ and is the safe default for hexagonal cells, slabs, magnetic systems, and tetrahedron method. Monkhorst-Pack shifts the grid by half a step so no point sits at Γ — better for some metallic cubic systems but breaks inversion-related savings in hexagonal lattices. When in doubt, choose Gamma.
4. Line 4. Subdivisions $N_1, N_2, N_3$ along $\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3$. This is what you change in a convergence sweep.
5. Line 5. Additional shift in fractional reciprocal coordinates. Almost always zero. Used in rare cases where you need to hit a specific high-symmetry point exactly.

Running a Convergence Test in Python

Before launching a 10-day VASP job on the right mesh, test the same idea in Python with a toy band. The 2-D tight-binding model from the interactive above is perfect: it has a Fermi surface (so it behaves like a metal), a closed-form band, and runs in microseconds. Click any line for the variable values and the why-each-argument story.

[[-0.589, -0.196, +0.196, +0.589],
 [-0.589, -0.196, +0.196, +0.589],
 [-0.589, -0.196, +0.196, +0.589],
 [-0.589, -0.196, +0.196, +0.589]]

[[-0.589, -0.589, -0.589, -0.589],
 [-0.196, -0.196, -0.196, -0.196],
 [+0.196, +0.196, +0.196, +0.196],
 [+0.589, +0.589, +0.589, +0.589]]

[[+2.828,  0.000,  0.000, +2.828],
 [ 0.000, -2.828, -2.828,  0.000],
 [ 0.000, -2.828, -2.828,  0.000],
 [+2.828,  0.000,  0.000, +2.828]]

[[F, F, F, F],
 [F, T, T, F],
 [F, T, T, F],
 [F, F, F, F]]

1import numpy as np
2
3def band_energy(kx, ky, t=1.0, a=4.0):
4    return -2.0 * t * (np.cos(kx * a) + np.cos(ky * a))
5
6def average_filled(N, a=4.0):
7    shifts = (np.arange(N) + 0.5) / N - 0.5
8    grid = 2 * np.pi / a * shifts
9    kx, ky = np.meshgrid(grid, grid)
10    eps = band_energy(kx, ky)
11    filled = eps < -1e-10  # tolerance dodges floating-point ε≈0 on the BZ edge
12    return eps[filled].mean()
13
14for N in [3, 4, 5, 6, 8, 12, 16, 24, 32]:
15    print(f"{N:3d}x{N:3d}  E = {average_filled(N):+.6f} eV")

Run it and look at the output. The energy lurches from $-1.60$ at $N=3$ down to $-2.83$ at $N=4$, back up to $-1.61$ at $N=5$, then creeps monotonically toward the converged value $-1.6339$ eV from below. That dramatic odd/even swing at low $N$ is the metallic Fermi-surface signature in action — exactly what we predicted from the integrand having a kink at $\varepsilon = 0$. The fix in a real VASP calculation is ISMEAR = 1, SIGMA = 0.1 — smear the step out and convergence becomes monotonic and fast.

Common Pitfalls

Summary

• Brillouin-zone integrals are unavoidable: every observable in electronic-structure theory is a BZ average of band quantities. We replace the integral by a weighted sum over a finite mesh.
• The mesh density that matters is $\Delta k_i = |\mathbf{b}_i|/N_i$, not $N_i$. Doubling the real-space cell halves the required N — the inverse-volume rule.
• Insulators converge fast (polynomial in $1/N$); metals converge slowly because the integrand has a Fermi-surface step. Smearing (ISMEAR) cures the metal case.
• The interactive 2-D viz makes the oscillation visible: even a textbook square-lattice Fermi surface produces a non-monotonic convergence curve until enough density is reached.
• Anisotropic cells need anisotropic grids. Slabs need $N_3 = 1$.
• Symmetry shrinks the mesh dramatically — VASP's irreducible BZ can be 8× smaller than the full mesh on cubic cells. Breaking symmetry with dopants undoes that saving.
• The simplest VASP KPOINTS file is five lines; KSPACING in INCAR is even simpler. Always run a convergence sweep on the property you care about, not just on total energy.
• The Python sweep above mirrors the VASP sweep one-to-one. If a toy band converges at $N \sim 30$, expect the same for a real metal at the same $\Delta k$.

Coming next: Section 3.9 — Monkhorst–Pack Grids — where we make the "regular grid" idea precise, derive the optimal half-step shift you saw above, and learn why Gamma and Monkhorst-Pack can give different irreducible point counts for the same cell.