k-point meshes

Periodic SCF replaces the molecular sum-over-occupied-orbitals with an integral over the first Brillouin zone (BZ). Practically that means picking a finite mesh of \(\mathbf{k}\)-points, solving the SCF self-consistently at each, and integrating to get the total energy. This page is the reference for how to choose, build, and interpret that mesh — Monkhorst-Pack vs Γ-centered, IBZ reduction, mesh convergence, band paths.

For the Bravais-lattice + reciprocal-lattice background — what the BZ actually looks like for FCC, BCC, hexagonal, etc., and where Γ, X, K, L, M sit — see crystal lattices.

Why we sample the BZ

For a periodic crystal the Bloch states \(\psi_{n\mathbf{k}}(\mathbf{r})\) are labeled by a band index \(n\) and a continuous wavevector \(\mathbf{k}\) in the first BZ. The total energy per unit cell is

\[ E = \frac{1}{V_\text{BZ}} \int_\text{BZ} \sum_{n \text{ occ}} \varepsilon_{n\mathbf{k}} \, d^3\mathbf{k} \]

(plus exchange-correlation, electron-electron, nuclear, etc. terms). We approximate this integral by a finite weighted sum over a mesh of points:

\[ E \approx \sum_i w_i \sum_{n \text{ occ}} \varepsilon_{n,\mathbf{k}_i} \quad \text{with} \quad \sum_i w_i = 1. \]

The choice of mesh — its size, its centring, whether it exploits crystal symmetry — controls how accurate the integral is. All properties (energy, forces, band gap, density of states) inherit that accuracy.

The Monkhorst-Pack mesh

The standard choice. Monkhorst and Pack (1976) proposed sampling on a uniform grid in fractional reciprocal-lattice coordinates:

\[ \mathbf{k}_{p,q,r} = \frac{2p - n_1 - 1}{2 n_1}\,\mathbf{b}_1 + \frac{2q - n_2 - 1}{2 n_2}\,\mathbf{b}_2 + \frac{2r - n_3 - 1}{2 n_3}\,\mathbf{b}_3 \]

for \(p = 1\ldots n_1\), \(q = 1\ldots n_2\), \(r = 1\ldots n_3\). The mesh has \(n_1 n_2 n_3\) points, each with weight \(1/(n_1 n_2 n_3)\) before any symmetry reduction.

Two common conventions for the origin of the mesh:

  • Centered (above) — point \((p,q,r) = ((n_1+1)/2, \ldots)\) sits at \(\mathbf{k} = 0\) only when all three \(n_i\) are odd. This is vibe-qc’s default, and it’s the original 1976 Monkhorst-Pack specification.

  • Γ-centered — shifted so \(\mathbf{k} = 0\) is always on the mesh. Useful for hexagonal lattices where the symmetry of the Γ point matters; some VASP-style workflows default to this.

In vibe-qc, the user-facing entry point is the :class:vibeqc.KPoints class (Phase K1 — first-class k-point ergonomics matching what users coming from VASP / CRYSTAL / Quantum ESPRESSO expect):

import numpy as np
import vibeqc as vq

# 8×8×8 Monkhorst-Pack mesh on a cubic crystal.
sysp = vq.PeriodicSystem(
    dim=3,
    lattice=8.0 * np.eye(3),
    unit_cell=[vq.Atom(11, [0, 0, 0])],
)
kp = vq.KPoints.monkhorst_pack(sysp, mesh=[8, 8, 8])

print(kp.n_kpoints)        # 512 (full mesh)
print(kp.kind)             # "monkhorst_pack"
print(kp.shift)            # (1, 1, 1) — classical-MP auto-shift
                           # for even meshes; (0, 0, 0) for odd
print(kp.is_symmetry_reduced)   # False — full mesh, no IBZ folding yet

Other common builders:

# Γ-centred mesh (k=0 is always on the grid):
kp = vq.KPoints.gamma_centred(sysp, [8, 8, 8])

# Γ-only (single k-point at the BZ origin):
kp = vq.KPoints.gamma(sysp)

# Custom shift on top of the MP base:
kp = vq.KPoints.monkhorst_pack(sysp, [8, 8, 8], shift=(0, 0, 1))

# From an explicit list of fractional k-points (Phase K4):
import numpy as np
kp = vq.KPoints.from_list(
    sysp,
    np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5], [0.5, 0.0, 0.0]]),
    weights=np.array([1, 8, 6]),   # auto-normalised to sum to 1
)

The mesh size is 3 integers, even for 1D and 2D systems — vibe-qc auto-clamps non-periodic axes to 1 k-point regardless of the input:

# 1D polymer: only mesh[0] matters; mesh[1] and mesh[2] silently → 1.
kp = vq.KPoints.monkhorst_pack(sysp_1d, [12, 4, 4])
print(kp.n_kpoints)      # 12, not 192

Tip

Back-compat. The old bare-function vq.monkhorst_pack(sysp, mesh) still ships and returns a :class:vibeqc.BlochKMesh directly (the lower-level type that periodic SCF drivers consume). KPoints.to_bloch_kmesh() adapts new-style KPoints to that format for any driver that hasn’t been updated to take KPoints natively. New code should prefer the KPoints builder; old code keeps working.

Density-based auto-meshes

For convenience — pick a mesh size from a single physically-motivated parameter rather than guessing a 3-integer mesh by hand. Three conventions ship (Phase K5):

# AFLOW convention: ~1000 k-points per reciprocal atom (Curtarolo 2012)
kp = vq.KPoints.from_kppra(sysp, n_kpts_per_atom=1000)

# Materials Project / ASE convention: target k-spacing in 2π/Å
kp = vq.KPoints.from_kspacing(sysp, dk_2pi_per_A=0.20)

# VASP "Auto" convention: a single length parameter
kp = vq.KPoints.auto(sysp, length_A=30.0)

Each builder takes an optional metallic=True flag that bumps the default density and warns if smearing isn’t enabled — metals need both fine k-meshes and Fermi-Dirac smearing to converge.

Equivalent k-points and IBZ reduction

The crystal’s space-group symmetry maps points in the BZ onto each other: \(\mathbf{k}\) and \(R\mathbf{k}\) (for \(R\) a point-group rotation) give identical energies and properties. Including both is redundant. The irreducible Brillouin zone (IBZ) is the smallest piece of the BZ such that no two points within it are symmetry-related; the rest is recoverable by applying the symmetry operations.

For an FCC lattice with cubic point-group \(O_h\) (48 operations), an 8×8×8 mesh has 512 points but only 35 in the IBZ — a 14× speedup on the SCF if you exploit it.

In vibe-qc you opt in via the KPoints.symmetry_reduce() builder method (Phase K2 — IBZ reduction via spglib):

from vibeqc import attach_symmetry

attach_symmetry(sysp)             # spglib detects the space group
kp = vq.KPoints.monkhorst_pack(sysp, [8, 8, 8])
kp_irr = kp.symmetry_reduce()
print(kp_irr.n_kpoints)           # ~35 for FCC
print(kp_irr.is_symmetry_reduced) # True
print(sum(kp_irr.weights))        # still 1.0 — symmetry-equivalent
                                   # points fold into boosted weights

The returned KPoints instance carries the multiplicity per point through the weights array; downstream periodic-SCF drivers see a smaller k-list with the right BZ-integration weighting.

Warning

Hex / trigonal cells (SG 143–194) refuse non-zero MP shifts — the classical (½,½,½) offset breaks the three-fold symmetry, so symmetry_reduce() raises an actionable error pointing at KPoints.gamma_centred(...) (the safe choice for those spacegroups).

Note

Caveat at v0.5: the multi-k SCF drivers (run_rhf_periodic_scf, run_rks_periodic_scf, run_uhf_periodic_multi_k_ewald3d, run_uks_periodic_multi_k_ewald3d) currently use the full mesh, not the IBZ. The IBZ reduction is a mesh-construction helper — useful for k-path / DOS / property post-processing, not yet plumbed through the SCF inner loop. SCF-side IBZ requires the D(R) transformation matrices for basis-function pairs under each symmetry operation; that’s tracked as a performance phase in the roadmap. Performance-wise, the SCF cost is therefore O(N_full) until the IBZ wiring lands.

Choosing the mesh size

The minimum mesh that gives a converged total energy depends on the system. The right approach is always: start with a small mesh, double, recompute, stop when the energy stops moving.

Rough starting points:

System type

Suggested starting mesh

Finite cluster in periodic box (molecular limit)

\([1, 1, 1]\)

1D polymer

\([6, 1, 1]\)

2D slab / surface

\([6, 6, 1]\)

3D wide-gap insulator (NaCl, MgO)

\([4, 4, 4]\)

3D narrow-gap semiconductor (Si, GaAs)

\([8, 8, 8]\)

3D metal

\([12, 12, 12]\) + smearing (see roadmap)

Anisotropic cells want anisotropic meshes. A surface slab with \(c \gg a\) doesn’t need many k-points perpendicular to the surface (the bands there are flat). HCP with \(c/a \approx 1.6\) wants roughly equal sampling in-plane and along \(c\).

# Slab: lots of in-plane k, very few perpendicular
kp = vq.KPoints.monkhorst_pack(sysp_slab, [8, 8, 1])

# HCP: mesh density along c roughly matches in-plane
kp = vq.KPoints.monkhorst_pack(sysp_hcp, [6, 6, 4])

The general rule of thumb is roughly equal density along each direction in reciprocal space: \(n_i \cdot |\mathbf{a}_i|\) should be roughly equal across \(i = 1, 2, 3\). So a long real-space direction wants fewer k-points along that direction.

Convergence study workflow

import vibeqc as vq

results = {}
for n in [2, 4, 6, 8, 10, 12]:
    kp = vq.KPoints.monkhorst_pack(sysp, [n, n, n])
    result = vq.run_rhf_periodic_scf(sysp, basis, opts,
                                      kmesh=kp.to_bloch_kmesh())
    results[n] = result.energy
    print(f"  n={n}:  E = {result.energy:.6f} Ha   |Δ| from n-2: "
          f"{abs(result.energy - results.get(n-2, result.energy)) * 1e3:.3f} mHa")

# Stop doubling when |Δ| drops below your target — typically 1 mHa
# for energies, 1 µHa for forces (forces converge faster than
# energy in the k-mesh limit).

Tutorial 12 walks through a complete convergence study on a 1D H-chain.

Tip

For metals (or near-metals), the total energy is discontinuous in the k-mesh — the Fermi level can drop below or rise above discrete eigenvalues as the mesh changes. The fix is smearing (Fermi-Dirac, Gaussian, or Methfessel-Paxton): occupations become fractional near the Fermi level, smoothing the integral. Vibe-qc’s Fermi-Dirac smearing landed in Phase C1b; without it, k-mesh convergence on metals is essentially impossible. Insulators don’t need smearing — the gap means the occupations are 0 or 2 by definition.

Bands: k-paths through the BZ

For band-structure plots you want a path of k-points along the edges of the IBZ, not a uniform mesh. The standard convention (Setyawan-Curtarolo 2010) defines a path per Bravais lattice:

Lattice

Standard path

Cubic (cP)

Γ → X → M → Γ → R → X → M → R

FCC (cF)

Γ → X → W → K → Γ → L → U → W → L → K

BCC (cI)

Γ → H → N → Γ → P → H → P → N

Hexagonal (hP)

Γ → M → K → Γ → A → L → H → A → L → M → K → H

Tetragonal (tP)

Γ → X → M → Γ → Z → R → A → Z → X → R → M → A

Each segment is sampled with N points (typically 30-60 per leg); the SCF result is post-processed to give \(\varepsilon_n(\mathbf{k})\) along the path. This is non-self-consistent — the density from the converged uniform-mesh SCF is fixed; only the band energies are re-evaluated at the path k-points.

In vibe-qc, the lattice-aware path is auto-detected via seekpath using the HPKOT 2017 convention (Phase K3). One call:

# Auto-detected path for the system's Bravais lattice (HPKOT/Hinuma):
kp_path = vq.KPoints.band_path(sysp, scheme="auto",
                                points_per_segment=30)
print(kp_path.kind)        # "band_path"
print(kp_path.n_kpoints)   # ~210 for FCC's standard 7-leg path

# Explicit override — pass your own segments:
kp_path = vq.KPoints.band_path(
    sysp,
    scheme="manual",
    segments=[("Γ", "X"), ("X", "W"), ("W", "K"), ("K", "Γ")],
    points_per_segment=40,
)

# Convert to the band-energy evaluation format:
bands = vq.compute_bands(sysp, basis, scf_result,
                          kp_path.to_kpath())
# bands.shape = (n_kpath_points, n_bands)

The auto-detection step calls spglib + seekpath under the hood: spglib identifies the Bravais lattice, seekpath returns the canonical path for that lattice in the HPKOT convention. Manual override is available for cases where you want a non-standard path or a specific high-symmetry segment.

See tutorial 12: band structure of an H-chain for the complete plotting workflow, and examples/periodic/input-h-chain-bands.py for the runnable script.

Density of states

The density of states (DOS) is

\[ g(\varepsilon) = \sum_n \int_\text{BZ} \delta(\varepsilon - \varepsilon_{n\mathbf{k}}) \, \frac{d^3\mathbf{k}}{V_\text{BZ}}. \]

Like the energy, this needs a finite k-mesh, but the right mesh for DOS is typically finer than for energy. Energy integrates the occupied states, where \(\varepsilon < \varepsilon_F\) is well-defined; DOS evaluates the integrand at every \(\varepsilon\), including in regions where the bands are flat (high DOS) — the mesh needs to resolve those features.

Practical rule: double the mesh size you used for energy convergence to get a smooth DOS. A 16×16×16 mesh on Si gives a visually clean DOS where 8×8×8 was sufficient for the total energy.

vibe-qc’s DOS API:

from vibeqc import compute_dos

kp_dense = vq.KPoints.monkhorst_pack(sysp, [16, 16, 16])
dos = vq.compute_dos(sysp, basis, scf_result,
                     kmesh=kp_dense.to_bloch_kmesh(),
                     sigma=0.05)  # eV broadening
# dos.energies, dos.dos

See the band-structure user guide for the full DOS + projected DOS (PDOS) workflow.

Smearing for metals

When the occupation function \(f(\varepsilon)\) is sharp (a step at \(\varepsilon_F\)), the BZ integral has discontinuities every time a band crosses the Fermi level — a finite mesh can’t resolve them cleanly. The fix is to broaden the step:

  • Fermi-Dirac: \(f = (1 + e^{(\varepsilon - \mu)/k_B T})^{-1}\). Has a thermodynamic interpretation (electronic temperature \(T\)); use this for finite-T studies.

  • Gaussian: \(f\) ≈ erf-shaped. No physical interpretation, but converges faster than Fermi-Dirac at the same broadening width.

  • Methfessel-Paxton: a hierarchy of corrections that approach the zero-T limit as the order increases. Common in plane-wave DFT.

Vibe-qc currently exposes Fermi-Dirac smearing through smearing_temperature_K on the periodic SCF options (Phase C1b). Gaussian and Methfessel-Paxton are tracked on the roadmap.

opts = vq.PeriodicSCFOptions()
opts.smearing_temperature_K = 300.0    # 300 K Fermi-Dirac
kp = vq.KPoints.monkhorst_pack(sysp, [12, 12, 12])
result = vq.run_rks_periodic_scf(sysp, basis, opts,
                                  kmesh=kp.to_bloch_kmesh())
# Reported entropy contribution: result.entropy_contribution

Worked example: NaCl rocksalt

Putting it together — a converged SCF on rocksalt NaCl with a properly-sized k-mesh:

import numpy as np
import vibeqc as vq
from vibeqc import attach_symmetry

# 1. Build the lattice + basis (FCC primitive + 2-atom basis).
a = 10.65   # bohr (NaCl: a = 5.64 Å)
sysp = vq.PeriodicSystem(
    dim=3,
    lattice=(a / 2) * np.array([
        [0, 1, 1], [1, 0, 1], [1, 1, 0],
    ]),
    unit_cell=[
        vq.Atom(11, [0.0, 0.0, 0.0]),
        vq.Atom(17, [0.5, 0.5, 0.5]),
    ],
)

# 2. Symmetry analysis — informs the k-mesh size.
attach_symmetry(sysp)

# 3. Convergence-screened k-mesh; 6×6×6 is fine for NaCl
#    (wide-gap insulator). symmetry_reduce() builds the IBZ via
#    spglib (Phase K2).
kp = vq.KPoints.monkhorst_pack(sysp, [6, 6, 6]).symmetry_reduce()
print(f"Full mesh: {6**3} points; IBZ: {kp.n_kpoints} points "
      f"({6**3 / kp.n_kpoints:.1f}× speedup if SCF supports it)")

# 4. Converged SCF.
basis = vq.BasisSet(sysp.unit_cell_molecule(), "pob-tzvp")
opts = vq.PeriodicSCFOptions()
opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D
result = vq.run_rks_periodic_scf(sysp, basis, opts,
                                 kmesh=kp.to_bloch_kmesh(),
                                 omega=0.5, spacing_bohr=0.5,
                                 functional="PBE")
print(f"E(NaCl) = {result.energy:.6f} Ha/cell  "
      f"({result.n_iter} iters, converged={result.converged})")

For more end-to-end periodic examples — including the Madelung- constant cross-validation that exercises the k=0 limit of this machinery — see the bundled example outputs.

References

  • H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations,” Phys. Rev. B 13, 5188 (1976) — defines the uniform-grid scheme that bears their name.

  • W. Setyawan and S. Curtarolo, “High-throughput electronic band structure calculations: Challenges and tools,” Comput. Mater. Sci. 49, 299 (2010) — defines the standard band-path conventions per Bravais lattice.

  • M. Methfessel and A. T. Paxton, “High-precision sampling for Brillouin-zone integration in metals,” Phys. Rev. B 40, 3616 (1989) — the corrected-Gaussian smearing scheme.

  • N. Marzari, D. Vanderbilt, and others, “Thermal contraction and disordering of the Al(110) surface,” Phys. Rev. Lett. 82, 3296 (1999) — discussion of finite-T smearing for metallic surfaces.

  • The VASP wiki KPOINTS page is the de-facto reference for plane-wave DFT k-mesh conventions; most of what’s there carries over to Gaussian-basis periodic SCF, modulo “smooth-cutoff” considerations that are PW-specific.