Monkhorst–Pack Grids

Learning Objectives

Section 3.8 made the case that any property of a crystal is, deep down, an integral over the Brillouin zone, and that we approximate that integral by sampling at a discrete set of k-points. This section answers the obvious next question: which discrete points? By the end you should be able to:

1. Write down the Monkhorst–Pack formula $u_r = (2r - N - 1) / (2N)$ and explain why this particular choice is "the most uniform" uniform mesh.
2. Decide between a Γ-centred and an MP-shifted mesh from the parity of $N$ and the symmetry of the lattice — and know why hexagonal cells demand the Γ-centred form.
3. Use the point group of the crystal to fold an $N_1 \times N_2 \times N_3$ mesh into a much smaller irreducible wedge, with multiplicities (weights) that sum back to the original count.
4. Read a real convergence curve $E_{\text{tot}}(N) - E_\infty$ and recognise why metals oscillate while insulators decay smoothly.
5. Apply the k·a ≈ const rule of thumb to set a consistent mesh density across anisotropic cells and supercells of different size.
6. Read, write, and debug a VASP KPOINTS file in automatic mode, and predict what its IBZKPT counterpart will contain.

One-line preview: Monkhorst and Pack's 1976 prescription is the unique uniform mesh whose centroid lies exactly at Γ for every direction independently — and whose half-step shift makes even $N$ avoid Γ on purpose, which turns out to be exactly what metallic Brillouin-zone integrals need.

Why a Grid? The Integral Hiding in Every DFT Run

Every observable VASP prints — total energy, Fermi level, density of states, dielectric tensor, magnetic moment — has the same skeleton:

The sum runs over bands $n$; the integral runs over the first Brillouin zone; $f_{n\mathbf{k}}$ is the Fermi–Dirac (or smearing-broadened) occupation. We never compute this integral analytically — the integrand is the output of a self-consistent Kohn–Sham diagonalisation that we do not have in closed form. The only practical option is quadrature: replace the integral by a weighted sum,

Three knobs control accuracy: (i) where the sample points $\mathbf{k}_i$ sit, (ii) how many of them you take, and (iii) how the weights $w_i$ are assigned. Ideally the points cover the BZ with the same density everywhere — anything biased towards one corner would systematically over- or under-count parts of the electronic structure. So we want a uniform mesh.

The Monkhorst–Pack Recipe

The 1976 paper (Phys. Rev. B 13, 5188) writes down a one-line formula for the k-points along each reciprocal-lattice direction:

That single line is the entire scheme. The k-points themselves are

Read the formula: the increment between consecutive $u_r$ is exactly $1/N$, so the points are evenly spaced. The smallest is $u_1 = (1-N)/(2N) = -1/2 + 1/(2N)$ and the largest is $u_N = (N-1)/(2N) = 1/2 - 1/(2N)$. The whole list sits symmetrically around zero — the centroid is exactly at Γ. That symmetry is the design principle Monkhorst and Pack engineered into the formula.

The Γ-centred alternative — what VASP calls Gamma — uses $u_r = (r-1)/N$, equivalent to shifting MP by $+1/(2N)$. That moves the mesh so that $u_1 = 0$ always — Γ is in. Algebraically the two schemes differ by a constant translation; physically they are very different when symmetry is involved (next section).

Even vs Odd: The Γ-Point Question

The single most important practical fact about MP grids is the parity rule:

Why does anyone want to avoid Γ? Because in a metal the partial occupancies $f_{n\mathbf{k}}$ have a step discontinuity at the Fermi surface. If the mesh hits exactly the Γ-point — which often coincides with a band extremum — you sample the integrand at a kink, and the quadrature error becomes first-order in $1/N$ instead of second-order. Even-N MP grids deliberately step around Γ, which restores the higher convergence rate. This is the technical reason the original paper argued for the half-step shift.

Hexagonal lattices are the famous exception. The C₃ rotation of the basal plane maps an even-N MP mesh onto a different mesh — the Brillouin-zone symmetry is broken by the sampling. The remedy is to force $u_1 = 0$ via Γ-centring, which keeps the mesh symmetric under all operations of the point group.

Interactive — Watch the Grid Fold

The widget below puts the formula in your hands. The grey square is the first Brillouin zone of a 2-D square lattice ($a_1 \perp a_2$, $|a_1| = |a_2| = a$); its symmetry group is C_4v. Slide $N_1$ and $N_2$, toggle the centring, and read the irreducible-point table. Each colour in the dot field is one irreducible representative — every dot of that colour is a member of its star.

Five things to do with the widget:

1. Set $N_1 = N_2 = 4$, MP centring. Note that no dot sits at Γ. Switch to Γ-centred — a dot appears at the origin and the entire mesh rigidly slides up-right by $1/(2N)$.
2. Stay at $4 \times 4$ MP. Count the irreducible points in the table — there are 3. Multiplicities are 4, 8, 4 → sum 16, matching $N_1 \times N_2$. Each irreducible k-point carries the "weight of its star."
3. Crank to $8 \times 8$. The dot field densifies; the irreducible count grows roughly as $N^2 / 8$ — the C_4v reduction factor. For most production work, the irreducible mesh is what dominates wall time, not the full mesh.
4. Try $N_1 = 5, N_2 = 5$ Γ-centred and MP. Now MP also includes Γ — that's the odd-N case. The irreducible counts differ slightly because the on-axis points sit on mirror lines and pick up extra symmetry.
5. Push $N_1 = 8, N_2 = 2$. The mesh is dense in $k_1$ but coarse in $k_2$. This is the wrong choice for a square BZ — but the right one if the real cell is squashed (anisotropic reciprocal vectors), which we discuss in the "k·a" rule below.

Symmetry, the Irreducible Wedge, and Weights

Look at the shaded triangle in the widget — the lower-right wedge with vertices Γ, X = $(\tfrac{1}{2}, 0)$, M = $(\tfrac{1}{2}, \tfrac{1}{2})$. That is the irreducible Brillouin zone (IBZ) of C_4v: one-eighth of the full BZ, the smallest region you can rotate or reflect to cover everything. Two k-points are physically equivalent iff the point group of the crystal sends one to the other, so we never need to compute eigenvalues at both — we work in the IBZ and then weight by multiplicity:

where the weight $w_{\mathbf{k}^*}$ is simply the number of full-BZ k-points that fold onto $\mathbf{k}^*$. For a generic interior point of the IBZ that number is $|G|$ (the order of the point group). Points on a mirror plane have weight $|G|/2$; points on the intersection of two mirrors carry $|G|/4$; the corners (Γ, X, M for C_4v) sit on the full stabiliser and have weight 1 each — their star has only one member.

Hover over a row in the table beside the widget — the corresponding star lights up in the dot field. Notice that the brighter, dark-edged dot in each colour is the canonical representative VASP would print to IBZKPT; the faint dots are its symmetry partners.

Wall-time scales with the irreducible count, so the symmetry reduction is a first-class optimisation, not a bookkeeping detail. For an FCC Brillouin zone (point group O_h, order 48) and a $16 \times 16 \times 16$ mesh, VASP reduces 4096 full-BZ points to roughly 145 irreducible — a 28× speed-up, for free, every SCF iteration.

Convergence: Metals vs Semiconductors

How dense is dense enough? You answer that question with a convergence test: run a series of identical calculations differing only in $N$, and watch an observable settle. The widget below shows representative shapes for a metal (FCC Al) and a wide-gap semiconductor (wurtzite CdSe), with the y-axis on a log scale so you can see four decades of error at once.

Three lessons that show up immediately:

1. Semiconductors are easy. The CdSe curve drops as a near-perfect exponential, hitting the 1 meV/atom bar by $N \approx 5$. Smooth integrand → exponential convergence. This is why every published CdSe study uses meshes around $6 \times 6 \times 4$.
2. Metals are hard. Aluminium oscillates with a period of ~2 in $N$ and only really sits below 1 meV/atom past $N \approx 16$. The oscillation is the partially-filled band edge moving across sample points; its amplitude is exactly what a Methfessel–Paxton or Marzari–Vanderbilt smearing tames (next chapter). Without smearing you would need $N \sim 30$ for a metal — and even then convergence is only first-order.
3. Density, not absolute count. The same $N = 8$ mesh on a small primitive cell is far denser in real-space-equivalent units than on a 2×2×2 supercell. We turn this observation into a quantitative rule next.

Anisotropic Cells and the k·a = const Rule

The reciprocal-space spacing along axis $i$ is

For a uniform sampling density we want $\Delta k_i$ the same in every direction, i.e.

The constant $\ell$ has units of length and is often quoted as the "k-mesh density" in Å. A common rule of thumb is $\ell \approx 30\,\text{Å}$ for a well-converged semiconductor and $\ell \approx 60\,\text{Å}$ for a metal.

For slabs and isolated molecules in vacuum ($|\mathbf{a}_3| \to \infty$), $|\mathbf{b}_3| \to 0$ — a single k-point along that axis is more than enough. $N_3 = 1$ is the right answer; spending more is pure waste.

Reading and Writing a VASP KPOINTS File

The KPOINTS file is the smallest file in a VASP calculation but it controls a remarkable amount. In automatic mode it is exactly five lines. Click any line below to see what VASP actually does with it.

1K-Points for Mn-doped CdSe primitive cell
20
3Monkhorst-Pack
48 8 6
50 0 0

What VASP writes to IBZKPT

After parsing the KPOINTS above, VASP writes the irreducible list to the IBZKPT file. For a wurtzite cell with C_6v symmetry and an $8 \times 8 \times 6$ mesh, the top of IBZKPT looks roughly like this (numbers schematic):

1Automatically generated mesh
2      36
3Reciprocal lattice
4    0.00000000000000    0.00000000000000    0.00000000000000             1
5    0.12500000000000    0.00000000000000    0.00000000000000             6
6    0.25000000000000    0.00000000000000    0.00000000000000             6
7    0.37500000000000    0.00000000000000    0.00000000000000             6
8    0.50000000000000    0.00000000000000    0.00000000000000             3
9    0.12500000000000    0.12500000000000    0.00000000000000             6
10    0.25000000000000    0.12500000000000    0.00000000000000           12
11    ...
12                                                       sum            384

The four columns are $(k_1, k_2, k_3)$ in fractional coordinates of the reciprocal basis (not Cartesian!) and the integer multiplicity$w_{\mathbf{k}}$. The sum of the multiplicities equals $N_1 \times N_2 \times N_3 = 384$ — a sanity check you should make every time. If the sum disagrees, VASP has applied the wrong symmetry, usually because ISYM in INCAR was set wrong or the POSCAR symmetry was broken by a typo.

Common Pitfalls

Summary

• Every observable in a DFT calculation is a Brillouin-zone integral. We approximate it with a uniform mesh of k-points and a quadrature rule.
• The Monkhorst–Pack formula $u_r = (2r - N - 1)/(2N)$ is the unique uniform mesh whose centroid sits at Γ for any $N$. The Γ-centred variant $u_r = (r-1)/N$ shifts the same mesh by $1/(2N)$ so that $u_1 = 0$.
• Even $N$ with MP misses Γ — exactly what you want for metals. Hexagonal cells must use Γ-centring to preserve C₃ symmetry.
• Point-group symmetry folds the $N_1 \times N_2 \times N_3$ mesh into a much smaller irreducible wedge; multiplicities (weights) recover the full sum. Wall time scales with the irreducible count, not the full count.
• Convergence behaviour distinguishes metals (oscillating, slow) from insulators (smooth, exponential). The right convergence target is 1 meV/atom on the quantity you care about, not just the total energy.
• The k·a ≈ ℓ rule sets a consistent k-mesh density across anisotropic cells and supercells. Use $\ell \approx 30\,\text{Å}$ for insulators, $60\,\text{Å}$ for metals, $N_3 = 1$ for slabs and molecules.
• A VASP KPOINTS file in automatic mode is exactly five lines: title, 0, scheme name, three integers, three-integer shift. The IBZKPT file VASP writes back tells you the irreducible set with weights — the sum of weights must equal $N_1 N_2 N_3$.

Coming next: Section 3.10 — Reciprocal Space in VASP — where we wire today's KPOINTS theory into a complete first-principles workflow: KSPACING, ISYM, the band-structure line mode, and how to interpret the eigenvalue tables in OUTCAR for the Mn:CdSe quantum-dot supercell that anchors the rest of the book.