Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

Across this chapter you have built the entire conceptual machinery of reciprocal space — the duality condition, the Brillouin zone, Bloch's theorem, structure factors, k-point grids, and the Monkhorst–Pack prescription. This closing section answers a single, very practical question: where does each of those concepts actually live inside a VASP calculation?

By the end of the section you should be able to:

Trace each of the five canonical VASP reciprocal-space artefacts — POSCAR, KPOINTS, IBZKPT, EIGENVAL, and DOSCAR — back to the exact mathematical object from this chapter that it represents.
Read and write all four KPOINTS modes (automatic mesh, line-mode for band paths, explicit list, and generalised regular grid) and pick the right one for SCF, DOS, and band-structure calculations.
Interpret VASP's IBZKPT file — including the integer weights — and verify that the irreducible point count matches what symmetry predicts.
Use the rule of thumb $n_i \,|\mathbf{b}_i| \;\gtrsim\; \ell$ (with $\ell$ in Å⁻¹) to pick a converged Monkhorst–Pack mesh from the lattice alone.
Build a one-screen Python helper that converts any VASP POSCAR into a converged KPOINTS file — the same script you will reuse throughout Chapter 6 for the Mn:CdSe simulation.

One-line preview: every reciprocal-space concept in Chapter 3 has a textual analogue in a VASP file. POSCAR carries the real lattice and so implicitly carries $\mathbf{b}_i$ ; KPOINTS picks $\mathbf{k}$ ; IBZKPT records the symmetry reduction; EIGENVAL stores $\varepsilon_n(\mathbf{k})$ ; DOSCAR is the Brillouin-zone integral of those eigenvalues.

The Reciprocal-Space Tour of a VASP Run

Most VASP first-timers think of the program as a black box that consumes INCAR, POSCAR, POTCAR, and KPOINTS, and emits OUTCAR. From a reciprocal-space perspective, however, you can read the entire workflow as one Brillouin-zone integral taken apart into discrete steps. The next table is the index for the rest of this section.

VASP file	What lives inside	Chapter 3 concept	Where covered
POSCAR	Direct lattice rows aᵢ in Å	Real basis. Implicitly fixes B = 2π·(A⁻¹)ᵀ.	§3.2
OUTCAR (header)	Both A and B printed (B without 2π)	Reciprocal lattice as a verifiable object.	§3.2
KPOINTS	Mesh density, centring, optional shifts; or a band path	Where in the BZ to sample. Monkhorst–Pack rule.	§3.8 — §3.9
IBZKPT	Irreducible k-points and integer weights	Symmetry-reduced BZ. Wedge of the point group.	Chapter 2 + §3.3
EIGENVAL	ε_n(k) for each band n at each irreducible k	Bloch eigenvalues — the values being summed.	§3.5
DOSCAR	Density of states n(ε) and partial DOS	Brillouin-zone integral of δ(ε − ε_n(k)).	Chapter 5
vasprun.xml	Structured XML containing all of the above	The one file post-processing tools really read.	Chapter 6

The mental model

Each row of the table is the same $\mathbf{k}$ -information re-expressed in a different form: as a lattice (POSCAR), as a sampling rule (KPOINTS), as a list (IBZKPT), as eigenvalues (EIGENVAL), as an integrated spectrum (DOSCAR). When you read a VASP run, you are walking down this column.

POSCAR → A → B — The Implicit Reciprocal Lattice

VASP never asks you for $\mathbf{b}_i$ explicitly. The program builds them on the fly from your POSCAR in the second line of subroutine latgen.f90 using exactly the closed-form recipe from §3.2: $\mathbf{B} = 2\pi\,(\mathbf{A}^{-1})^{\!\top}$ . The first thing OUTCAR prints back is the result, divided by $2\pi$ , in units of Å⁻¹.

📝text

1direct lattice vectors                 reciprocal lattice vectors
2   0.000000000  3.038500000  3.038500000   -0.164554  0.164554  0.164554
3   3.038500000  0.000000000  3.038500000    0.164554 -0.164554  0.164554
4   3.038500000  3.038500000  0.000000000    0.164554  0.164554 -0.164554
5
6  length of vectors
7       4.298  4.298  4.298     0.285  0.285  0.285

Multiply each printed reciprocal row by $2\pi$ and you recover the textbook $\mathbf{b}_i = (2\pi/a)(\pm 1, \pm 1, \pm 1)$ for FCC CdSe with $a = 6.077\,\text{Å}$ . The columns match $(2\pi/a) \approx 1.034\,\text{Å}^{-1}$ exactly.

A 2-second VASP sanity check

Open OUTCAR and grep for reciprocal lattice vectors. The three rows should have the same $|\mathbf{b}_i|$ as your hand calculation to better than 1 part in 10⁵. If they don't, your POSCAR scale factor is wrong — the single most common typo new VASP users make.

Anatomy of the KPOINTS File

After the lattice itself, the KPOINTS file is the single most consequential reciprocal-space input you provide. It lives in a strict 5-line ASCII format, but those five lines encode four very different sampling strategies. Every working VASP practitioner can recognise all four on sight.

Mode 1 — Automatic Monkhorst–Pack mesh

The default, used for SCF and DOS calculations. The integer on line 2 is zero — the cue that lines 3–5 below describe a generator rather than an explicit list.

📝text

1Automatic mesh
20
3Gamma
412 12 12
50  0  0

Reading top to bottom:

Line 1. Free-form comment.
Line 2. Number of k-points. 0 means "generate them automatically".
Line 3. Centring keyword. Gamma (Γ-centred) places one point exactly on Γ. Monkhorst-pack shifts the entire mesh by half a step so Γ falls between points. For hexagonal lattices Γ-centred is mandatory; for cubic it is the safe default.
Line 4. Mesh subdivisions $(n_1, n_2, n_3)$ along $\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3$ .
Line 5. Additional fractional shift $(s_1, s_2, s_3)$ . Almost always 0 0 0; non-zero shifts are an advanced trick to break a symmetry on purpose.

Mode 2 — Line mode for band structures

For plotting $\varepsilon_n(\mathbf{k})$ along the high-symmetry path of §3.4, switch to line mode. Line 2 now sets how many k-points to interpolate per segment. Each consecutive pair of non-blank coordinate lines defines one segment.

📝text

1FCC band path  Γ-X-W-L-Γ-K
220
3Line
4Reciprocal
50.000  0.000  0.000   ! Γ
60.500  0.000  0.500   ! X
7
80.500  0.000  0.500   ! X
90.500  0.250  0.750   ! W
10
110.500  0.250  0.750   ! W
120.500  0.500  0.500   ! L
13
140.500  0.500  0.500   ! L
150.000  0.000  0.000   ! Γ
16
170.000  0.000  0.000   ! Γ
180.375  0.375  0.750   ! K

Read this as: from Γ to X interpolate 20 points; from X to W another 20; and so on along the standard Setyawan–Curtarolo path Γ–X–W–L–Γ–K. The keyword Reciprocal on line 4 tells VASP that the coordinates are fractional in units of $\mathbf{b}_i$ — never confuse this with Cartesian, which expects Å⁻¹.

Modes 3 and 4 — Explicit lists and generalised grids

Less common, but you will eventually meet them. Explicit-list mode (line 2 set to a positive integer $N$ ) feeds VASP $N$ hand-picked k-points with weights — the format VASP itself uses to write IBZKPT. Generalised-grid mode (line 4 contains three reciprocal basis rows instead of three integers) lets you sample on a non-cubic-aligned grid; it is essentially how phonopy and related tools talk to VASP.

Choosing centring on a non-cubic system

For any lattice whose $\mathbf{b}_i$ are not orthogonal — hexagonal, rhombohedral, monoclinic, triclinic — Monkhorst–Pack centring can place k-points off the symmetry axes, breaking degeneracies that the crystal actually has. Γ-centred is the correct default. The performance penalty is zero; the reduction in surprise is enormous.

KSPACING — the One-Line INCAR Alternative

The KPOINTS file is the classical interface, but it has one annoying feature: every time the cell shape changes (relaxation, equation of state, scan over lattice constant) the right mesh changes too, and you have to remember to regenerate KPOINTS. VASP's answer is the KSPACING tag, an INCAR option that builds the mesh on the fly from the current cell.

Set it once, in INCAR:

📝text

1# INCAR fragment
2KSPACING = 0.20      ! target spacing between k-points, in 1/Angstrom
3KGAMMA   = .TRUE.    ! Γ-centred mesh (mandatory for hex; recommended for cubic)

VASP then generates a Γ-centred Monkhorst–Pack mesh whose subdivisions along each reciprocal axis are

n_i \;=\; \max\!\left(1,\; \Bigl\lceil \frac{|\mathbf{b}_i|}{\text{KSPACING}} \Bigr\rceil\right).

That is exactly the inverse-volume rule from §3.8 expressed as a single number, and it carries the same physical meaning: keep the spacing between adjacent k-points constant in 1/Å, so a stretched cell (smaller $|\mathbf{b}_i|$ ) automatically gets a sparser mesh and a compressed cell automatically gets a denser one. If you write KSPACING = 0.20 in INCAR and provide no KPOINTS file at all, VASP picks the mesh on its own; if you do provide a KPOINTS file, that file wins.

KSPACING (Å⁻¹)	Roughly equivalent to ℓ in §3.8	When to use
0.40	ℓ ≈ 16 Å⁻¹	Insulator/semiconductor relaxation; coarse
0.30	ℓ ≈ 21 Å⁻¹	Insulator total energy; safe default
0.20	ℓ ≈ 31 Å⁻¹	Semiconductor DOS, defects, large supercells
0.15	ℓ ≈ 42 Å⁻¹	Metals, magnetism, optical properties
0.10	ℓ ≈ 63 Å⁻¹	Fermi-surface integrals, GW, hybrids

Why two interfaces for the same thing?

The KPOINTS file gives you total control — line mode for band structures, explicit lists, generalised grids, custom shifts. KSPACING gives you total convenience — set one number, never regenerate. The community split is roughly: production throughput (high-throughput databases, AiiDA workflows, relaxation campaigns) uses KSPACING; per-system one-offs (a band structure for a paper) uses an explicit KPOINTS file. Use whichever matches the task; both produce identical numbers if you match $n_i = \lceil |\mathbf{b}_i|/\text{KSPACING}\rceil$ .

KSPACING and band structures do not mix

Line-mode k-points cannot be auto-generated from KSPACING — it is an SCF-only convenience. If you want a band plot, you switch to an explicit KPOINTS file with Line-mode set, regardless of what is in INCAR. (This is the recipe in the next section.)

The IBZKPT File — Symmetry in Plain Text

After VASP reads KPOINTS it expands the mesh, applies every operation of the crystal's point group, and discards symmetry-equivalent points. The result is dumped into IBZKPT:

📝text

1Automatically generated mesh
2       72
3Reciprocal lattice
4   0.00000000  0.00000000  0.00000000      1
5   0.08333333  0.00000000  0.00000000      8
6   0.16666667  0.00000000  0.00000000      8
7   0.25000000  0.00000000  0.00000000      8
8   0.33333333  0.00000000  0.00000000      8
9   ...
10 Tetrahedra
11       6912    0.0000010851
12        1     1     2     3   ...

The four columns mean exactly:

Lines 1–2. Header. Line 2 is the count of irreducible k-points the run will actually evaluate. For a 12×12×12 mesh on FCC, with full Oh symmetry, this number is 72.
Line 3. Coordinate type — almost always Reciprocal (fractional in $\mathbf{b}_i$ ).
Lines 4 onward. Three fractional coordinates $(k_1, k_2, k_3)$ followed by an integer weight equal to the size of that point's symmetry orbit.

The weights add up to $n_1 \, n_2 \, n_3$ — for a 12×12×12 mesh, 1 728. That is the integer identity which guarantees nothing has been double-counted when VASP later evaluates a BZ-average like

\langle f \rangle \;=\; \frac{1}{n_1 n_2 n_3}\sum_{i \in \text{IBZ}} w_i \, f(\mathbf{k}_i).

The sum runs only over the irreducible set, weighted by $w_i$ — that is the entire reason IBZKPT exists. VASP also writes a Tetrahedra block at the bottom; it is the connectivity list used by the Blöchl tetrahedron method (we will meet it again in Chapter 5 when computing DOS).

Reading the speedup

Divide the total mesh count $n_1 n_2 n_3$ by the irreducible count and you get the symmetry speedup. For 12×12×12 on FCC: 1 728 / 72 = 24× fewer Hamiltonian diagonalisations than a brute-force calculation. This factor scales like the order of the point group divided by 2 (the 2 is time-reversal symmetry, which VASP applies automatically when ISYM ≥ 1).

Interactive — Monkhorst–Pack Mesh in the BZ

The picture below is the geometric content of IBZKPT: the FCC Brillouin zone (truncated octahedron, §3.3) with every MP grid point folded into it. Amber dots are the irreducible representatives — the only ones VASP actually solves. Slate dots are the symmetry-equivalent partners that VASP knows about but never diagonalises. Slide the mesh density and watch the irreducible count grow far slower than the total — exactly the speedup the previous Tip just described.

Loading 3D visualisation…

Three things to try while you play with the slider:

Set $n = 4$ with Γ-centring. You should read64 total, 8 irreducible, ×8 reduction in the side panel — the textbook FCC numbers.
Toggle Γ-centring off (Monkhorst shift). Watch the irreducible count change for odd meshes — the shift can accidentally place mesh points on a symmetry-violating site and the reduction factor drops.
Push $n$ to 8. The grid now has 512 points but the irreducible set is still well under 50. This is the asymptotic regime where a denser mesh costs almost nothing — and the regime where you should run convergence tests in.

What this viewer is not

The reduction is computed with the full cubic point group (Oh, 48 elements). For a real Mn:CdSe supercell the substitutional Mn breaks Oh down to a smaller subgroup, and the irreducible count grows accordingly. VASP handles that automatically through ISYM; this viewer is the unbroken-symmetry baseline.

Reciprocal-Space Outputs — EIGENVAL, DOSCAR, vasprun.xml

Once the SCF cycle finishes (Chapter 4), VASP emits three reciprocal-space artefacts. None of them is hard to read by hand once you know the format.

EIGENVAL — Bloch eigenvalues per k-point

EIGENVAL pairs each irreducible k-point (with its weight) to the column of band energies $\varepsilon_n(\mathbf{k})$ . Skip the first six header lines and you find a repeating block: a blank line, the four-column k-point/weight row, then one line per band.

📝text

11   72  ...
2                                                    <-- header lines
3
4  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.4629630E-03
5   1     -10.4321
6   2      -8.1245
7   3      -5.9871
8   ...
9
10  0.8333333E-01  0.0000000E+00  0.0000000E+00  0.3703704E-02
11   1     -10.2917
12   2      -7.9034
13   ...

For the 12×12×12 CdSe mesh you would see 72 such blocks back to back. Plotting $\varepsilon_n(\mathbf{k})$ against the cumulative path length gives the band structure we will analyse in Chapter 5.

DOSCAR — Brillouin-zone integral

DOSCAR is the BZ integral of those eigenvalues:

n(\varepsilon) \;=\; \sum_n \frac{1}{(2\pi)^3}\int_{\text{BZ}} \delta\!\left(\varepsilon - \varepsilon_n(\mathbf{k})\right) d^3 k.

Numerically this is the weighted sum over the irreducible set, broadened either by tetrahedra (the connectivity list at the bottom of IBZKPT) or by a Gaussian smearing of width SIGMA. The reciprocal-space bookkeeping is the same; only the broadening kernel changes.

vasprun.xml — the one file post-processors actually read

Tools like pymatgen, ase, sumo, and VASPKIT never look at OUTCAR; they parse vasprun.xml. The XML contains everything in the plain-text files plus the lattice, the k-mesh, and per-orbital projections — all under structured tags. If you ever automate post-processing, this is the file you target.

The Band-Structure Workflow — SCF → Non-SCF → Line Mode

Section 3.4 told us where to draw a band path; §3.10 has so far told us how to fill a single KPOINTS file. But a band-structure plot in VASP is the result of two VASP runs, not one — and the file that travels between them is the single most-misused file in this whole pipeline. The recipe is short, mechanical, and where 90 % of beginner band plots silently go wrong.

Why two runs?

The Hamiltonian $\hat{H} = -\tfrac{\hbar^2}{2m}\nabla^2 + V_{\text{eff}}[\rho]$ depends on the self-consistent density $\rho(\mathbf{r})$ , and $\rho$ in turn depends on the eigenstates through $\rho = \sum_{n,\mathbf{k}} f_{n\mathbf{k}}\,|\psi_{n\mathbf{k}}|^2$ . The only honest way to find $\rho$ is to integrate over a uniform sample of the Brillouin zone — the dense Monkhorst–Pack mesh of §3.9. The line-mode path of §3.4, with all its k-points stacked along Γ–X–W–L–Γ–K, samples a single one-dimensional curve through the BZ and is useless for that integral. Run it self-consistently and you would converge to a density that thinks the Brillouin zone is a string instead of a volume.

The fix is the two-step pattern below, used by every band-structure tool and tutorial since the 1990s:

SCF run on a dense MP mesh. Converge the density $\rho$ properly. Save it to CHGCAR.
Non-SCF run on the line-mode path. Read CHGCAR, build $\hat{H}$ from it once, and diagonalise at each path k-point without ever touching $\rho$ again.

The INCAR tag that does the magic is ICHARG = 11: tens digit 1 means "read CHGCAR", ones digit 1 means "do not update it". Pair that with LCHARG = .FALSE. in the band run (so VASP does not overwrite the carefully-converged CHGCAR that you just copied in) and you have the entire pattern in three lines.

The orchestration script — every line annotated

The bash script below runs the full workflow for FCC CdSe, end to end. Click any line to read the meaning of that command, the role of every flag, and what file it produces or consumes. The two critical lines — LCHARG = .TRUE. in the SCF INCAR and ICHARG = 11 in the band INCAR — are highlighted in their cards.

Band-structure workflow — SCF → CHGCAR → line-mode bands

⚡run_bands.sh

Explanation(13)

Code(75)

1#!/usr/bin/env bash

Shebang line. Tells the operating system to run the rest of this file with /usr/bin/env bash. Using env keeps the script portable: it picks up bash from the user's PATH instead of hard-coding /bin/bash, so the same file works on Linux, macOS, and most HPC clusters.

EXECUTION STATE

#! = Two-byte magic that the kernel inspects when you exec a file. The interpreter named after #! is invoked with the script's path as its argument.

/usr/bin/env bash = Resolves bash via PATH. Recommended over /bin/bash because some systems install bash elsewhere (Apple Silicon Homebrew puts it in /opt/homebrew/bin, for example).

4Step 1 — mkdir -p band-cdse/scf band-cdse/bands

Create both run directories in one call. The -p flag means: create parent directories as needed and do not error if they already exist. We separate SCF and band runs into different folders so VASP outputs (OUTCAR, vasprun.xml, etc.) do not overwrite each other and so we can re-run either step independently.

EXECUTION STATE

📚 mkdir = POSIX command to create directories. Without -p, mkdir errors if the parent does not exist or the directory already exists.

⬇ arg: -p = 'parents' flag. Creates missing parent directories silently and does not fail on existing directories. Idempotent — safe to re-run.

⬇ arg: band-cdse/scf, band-cdse/bands = Two paths. mkdir -p creates band-cdse first, then both subdirectories. After this line: band-cdse/, band-cdse/scf/, band-cdse/bands/ all exist.

5cp POSCAR POTCAR band-cdse/scf/

Stage the SCF inputs. POSCAR (lattice + atomic positions) and POTCAR (PAW pseudopotentials) are the two inputs that never change between SCF and band runs. We copy rather than symlink so the run directory is self-contained — useful when the cluster runs on a different filesystem than the staging area.

EXECUTION STATE

📚 cp = POSIX file copy. cp src1 src2 ... dst/ copies each source into the destination directory.

⬇ arg: POSCAR = Atomic structure file. Header has the lattice matrix in Å and a list of fractional or Cartesian coordinates. VASP rebuilds the reciprocal lattice B from this on every run (§3.10 'POSCAR → A → B').

⬇ arg: POTCAR = PAW pseudopotential bundle, concatenated species by species in the same order as POSCAR's species list. Provides ENMAX (default ENCUT) and the projector functions.

⬇ arg: band-cdse/scf/ = Destination directory (trailing slash optional but recommended — it forces cp to treat the target as a directory and not as a rename).

6cd band-cdse/scf

Move into the SCF directory. Subsequent cat redirections (next two blocks) write INCAR and KPOINTS here, and mpirun launches VASP from this directory.

EXECUTION STATE

📚 cd = Change working directory. Affects this shell only — cd in a script does not affect the parent shell.

9Step 2 — cat > INCAR <<'EOF' ... EOF (SCF INCAR)

Heredoc redirection. The text between <<'EOF' and the closing EOF is written verbatim to the file INCAR. Single-quoting the EOF disables shell variable expansion inside the heredoc — what you see is exactly what lands on disk. The eight lines that follow are the entire INCAR for the SCF run.

EXECUTION STATE

📚 cat > FILE <<'EOF' ... EOF = Bash heredoc + redirection. cat reads from stdin, > redirects stdout to FILE, and the heredoc supplies stdin from the lines until the matching EOF marker.

PREC = Accurate = Sets a coherent set of grid densities and projection precision. Use Accurate (not Normal) whenever total-energy differences matter, which they always do for bands.

ENCUT = 350 = Plane-wave cutoff in eV. The basis includes every G with |k+G|² · ℏ²/2m ≤ ENCUT (§3.5, §4.9). 350 eV is well above the largest ENMAX in the POTCAR for Cd and Se. Identical ENCUT in SCF and band runs is mandatory — otherwise the basis sizes mismatch and VASP refuses to read the CHGCAR.

ISMEAR = 0, SIGMA = 0.05 = Gaussian smearing of width 0.05 eV. Safe semiconductor default — keeps the integration well-behaved without smearing across the band gap. Do not use ISMEAR = 1 or 2 (Methfessel–Paxton) here; those are designed for metals and produce a small, non-physical electron count when applied to a gapped material.

EDIFF = 1E-6 = Self-consistency tolerance on the total energy in eV. 1E-6 is more than tight enough for bands — looser tolerances leave noise in the high-ε bands.

LWAVE = .FALSE. = Do not write WAVECAR. The band run reads CHGCAR, not WAVECAR, so saving a 200-MB WAVECAR is wasted disk for this workflow.

LCHARG = .TRUE. ⚑ critical = Write CHGCAR. This is the entire point of the SCF: produce a converged self-consistent density that the band run will read. Forgetting LCHARG = .TRUE. is the #1 reason a band-structure workflow silently fails.

19Step 3 — cat > KPOINTS (SCF mesh)

Standard 5-line KPOINTS file (Mode 1 from §3.10). A 12×12×12 Γ-centred Monkhorst–Pack mesh on the FCC primitive cell folds, by symmetry (Oh + time reversal), to 72 irreducible k-points — the IBZKPT file confirms this after the run.

EXECUTION STATE

Line 1: 'SCF mesh' = Free-form comment, ignored by VASP.

Line 2: 0 = Number of k-points = 0 ⇒ generate them automatically from the rest of the file.

Line 3: Gamma = Γ-centred mesh. Mandatory for non-orthogonal reciprocal axes; safe default for cubic. Avoid Monkhorst-pack on non-cubic to keep symmetry-protected degeneracies intact.

Line 4: 12 12 12 = Mesh subdivisions n₁ n₂ n₃ along b₁, b₂, b₃. With |b_i| ≈ 1.79 Å⁻¹ for FCC CdSe, this gives a k-spacing of 0.15 Å⁻¹ — tight enough for total-energy bands.

Line 5: 0 0 0 = Optional fractional shift; almost always zeros.

28Step 4 — mpirun -np 32 vasp_std (SCF)

Launch VASP with 32 MPI ranks. After ~50 SCF iterations VASP writes CHGCAR (the binary self-consistent density), CONTCAR (= POSCAR if no relaxation), OSZICAR (energy convergence log), and OUTCAR (verbose log). The next step needs only CHGCAR.

EXECUTION STATE

📚 mpirun -np N = Launch N MPI processes of the executable that follows. On Slurm clusters you would normally use srun --ntasks=N instead and let the scheduler decide the allocation.

vasp_std = Standard VASP binary (Γ-only and non-collinear builds are vasp_gam and vasp_ncl). For FCC CdSe with several non-equivalent k-points, the standard build is the right pick.

⬆ result: CHGCAR = Binary file with the converged charge density on the FFT grid. Size ≈ 2 NGX·NGY·NGZ doubles. This is the carrier of self-consistent information into the band run.

31Step 5 — cd ../bands

Hop sideways to the band-structure directory. We deliberately do not run the band step in the SCF directory because the band run will overwrite EIGENVAL, OUTCAR, and friends — keeping them in different folders preserves the SCF outputs for inspection.

32cp ../scf/POSCAR ../scf/POTCAR ../scf/CHGCAR .

Copy the three files the band run needs from the SCF directory into the current directory. POSCAR and POTCAR define the system; CHGCAR carries the converged density. The lattice in POSCAR must be identical to the SCF (it is — we copied the same file) so the FFT grid lines up bit-for-bit.

EXECUTION STATE

⬇ arg: ../scf/CHGCAR ⚑ critical = The bridge between the two runs. ICHARG = 11 in the next INCAR will read this file and freeze it. If you forget to copy it, ICHARG = 11 silently falls back to a superposition of atomic densities — and your bands will be wrong, often subtly, but always wrong.

⬇ arg: . (dot) = Destination — the current directory (band-cdse/bands). Equivalent to ./ — three files copied here.

35Step 6 — cat > INCAR <<'EOF' ... EOF (band INCAR)

Same heredoc trick, different INCAR. The two changed-from-SCF lines are ICHARG = 11 (read CHGCAR, do not update) and the LWAVE/LCHARG flags (we do not need to write either; the bands themselves are written to EIGENVAL). NBANDS and LORBIT are added for the band-plotting workflow.

EXECUTION STATE

PREC, ENCUT, ISMEAR, SIGMA = Identical to the SCF — must be identical, otherwise the basis size implied by ENCUT changes and CHGCAR cannot be read.

ICHARG = 11 ⚑ critical = Tens digit 1 → read CHGCAR; ones digit 1 → do not update it. So this is exactly 'frozen self-consistent density'. With ICHARG = 1 (no tens digit) VASP would still read CHGCAR but then overwrite it during a fresh SCF — which we do not want for a band run.

LWAVE = .FALSE. = WAVECAR not needed for plotting bands; saves disk and I/O.

LCHARG = .FALSE. = Do not overwrite the CHGCAR we just copied in.

NBANDS = 24 = Number of bands to compute at every k-point in the path. Must include enough conduction states to see the band gap and beyond. For CdSe with 18 valence electrons in the primitive cell, 9 occupied bands + 15 conduction bands = 24 is comfortable.

LORBIT = 11 = Write site- and orbital-projected weights into PROCAR. Lets a post-processor (sumo, pymatgen, vaspkit) colour the bands by orbital character — Cd 5s, Se 4p, etc. — which is the difference between a useful band plot and a wall of black lines.

47Step 7 — cat > KPOINTS <<'EOF' (line mode)

Line-mode KPOINTS (Mode 2 from §3.10). Each pair of non-blank coordinate lines defines a segment of the high-symmetry path; blank lines separate segments. The integer on line 2 (40) is the number of k-points to interpolate per segment — finer = smoother band plot, coarser = faster.

EXECUTION STATE

Line 2: 40 = Points per segment. With 5 segments (Γ-X, X-W, W-L, L-Γ, Γ-K), 5 × 40 = 200 k-points total. Each k-point is one diagonalisation of the (frozen-density) Hamiltonian — typically a few seconds, so the band run is much cheaper than the SCF.

Line 3: Line-mode = Switches the mesh generator off and the segment interpolator on. Anything other than the prefix L (case-insensitive) leaves you in Mode 1 — a frequent silent typo.

Line 4: Reciprocal = Coordinates are fractional in units of b₁, b₂, b₃. The alternative is Cartesian (units of 1/Å) — which would require multiplying every coordinate by 2π/a. Using Reciprocal makes the coordinates portable across cell sizes.

The path Γ-X-W-L-Γ-K = Setyawan & Curtarolo, Comp. Mat. Sci. 49 (2010): canonical FCC path that visits one representative point per orbit of the irreducible BZ. Reading the (h k ℓ) of each point: Γ = (0,0,0); X = (½,0,½); W = (½,¼,¾); L = (½,½,½); K = (⅜,⅜,¾).

70Step 8 — mpirun -np 32 vasp_std (band run)

Launch VASP again, this time on the band-mode KPOINTS. Because ICHARG = 11 fixes the density, no SCF cycle runs — VASP simply diagonalises H at each of the 200 k-points and writes the eigenvalues into EIGENVAL. Wall time is typically 5–20 minutes for a small cell, much shorter than the SCF.

EXECUTION STATE

⬆ result: EIGENVAL = Plain-text file: header (number of electrons, nk, NBANDS), then for each k-point the fractional coords + weight, then NBANDS eigenvalues + occupations. This is the file the band-plot tool will read.

⬆ result: PROCAR = Plain-text orbital projections, one block per k-point per band. Required for orbital-coloured band plots.

⬆ result: vasprun.xml = Structured XML containing both EIGENVAL and PROCAR data plus the lattice and k-mesh metadata. The file every modern post-processor reads first.

73Step 9 — sumo-bandplot --filename vasprun.xml -o bands.png

Render ε_n(k) along the path. sumo (Scanlon Materials Theory Group) is a Python tool that parses vasprun.xml, identifies the path segments from the line-mode KPOINTS, sets the energy zero at the valence-band maximum (or Fermi level for metals), and produces a publication-quality plot. Equivalents: pymatgen.electronic_structure.plotter, vaspkit option 211.

EXECUTION STATE

📚 sumo-bandplot = Command-line entry point of the sumo Python package. Install via pip install sumo. Generates SVG, PDF, or PNG and writes a band-data CSV alongside.

⬇ arg: --filename vasprun.xml = Source file. sumo also accepts EIGENVAL but vasprun.xml is preferred — it has the lattice + k-path encoded in the same file, so axis tick labels (Γ, X, W, L, K) are automatic.

⬇ arg: -o bands.png = Output file. Extension picks the format: bands.png for slides, bands.pdf for a paper, bands.svg for vector editing.

62 lines without explanation

1#!/usr/bin/env bash
2# Band structure for FCC CdSe — two VASP runs, one density.
3
4# 1. Two-directory layout: SCF first, bands second
5mkdir -p band-cdse/scf band-cdse/bands
6cp POSCAR POTCAR band-cdse/scf/
7cd band-cdse/scf
8
9# 2. SCF INCAR — converge the density on a dense MP mesh
10cat > INCAR <<'EOF'
11PREC   = Accurate
12ENCUT  = 350
13ISMEAR = 0
14SIGMA  = 0.05
15EDIFF  = 1E-6
16LWAVE  = .FALSE.
17LCHARG = .TRUE.
18EOF
19
20# 3. Dense Γ-centred KPOINTS for the SCF
21cat > KPOINTS <<'EOF'
22SCF mesh
230
24Gamma
2512 12 12
260 0 0
27EOF
28
29# 4. Run SCF — writes CHGCAR (the converged charge density)
30mpirun -np 32 vasp_std
31
32# 5. Hand the density to the band-structure run
33cd ../bands
34cp ../scf/POSCAR ../scf/POTCAR ../scf/CHGCAR .
35
36# 6. Non-SCF INCAR — read CHGCAR, do not touch it
37cat > INCAR <<'EOF'
38PREC   = Accurate
39ENCUT  = 350
40ISMEAR = 0
41SIGMA  = 0.05
42ICHARG = 11
43LWAVE  = .FALSE.
44LCHARG = .FALSE.
45NBANDS = 24
46LORBIT = 11
47EOF
48
49# 7. KPOINTS in line mode along the FCC path Γ-X-W-L-Γ-K
50cat > KPOINTS <<'EOF'
51FCC band path
5240
53Line-mode
54Reciprocal
550.000 0.000 0.000   ! Gamma
560.500 0.000 0.500   ! X
57
580.500 0.000 0.500   ! X
590.500 0.250 0.750   ! W
60
610.500 0.250 0.750   ! W
620.500 0.500 0.500   ! L
63
640.500 0.500 0.500   ! L
650.000 0.000 0.000   ! Gamma
66
670.000 0.000 0.000   ! Gamma
680.375 0.375 0.750   ! K
69EOF
70
71# 8. Run non-SCF — writes EIGENVAL with bands along the path
72mpirun -np 32 vasp_std
73
74# 9. Plot ε_n(k)
75sumo-bandplot --filename vasprun.xml -o bands.png

The four ways this workflow silently fails

Forgetting LCHARG = .TRUE. in the SCF. CHGCAR is never written; the band run quietly falls back to a superposition of atomic densities and produces bands that look roughly right but mis-place every avoided crossing.
Different ENCUT in the two runs. The plane-wave basis size is set by ENCUT, and CHGCAR is stored on a grid sized to match. VASP refuses to read a CHGCAR whose grid does not match the current ENCUT — but the error message is short and easy to miss in a 30-MB OUTCAR.
Using KSPACING in the band INCAR. KSPACING auto-generates a Monkhorst–Pack mesh and ignores the line-mode KPOINTS file entirely. You will get a dense uniform sampling and no band path. Either remove KSPACING or set it higher than the path will ever reach.
ICHARG = 1 instead of 11. ICHARG = 1 reads CHGCAR and runs SCF on top of it on the line-mode path — which converges to a wrong density (string instead of volume) and gives wrong bands. Always 11 for the band run.

From this point onward

Every band-structure plot you see for the rest of this book — the gap engineering of §5.3, the orbital-projected bands of §5.4, the Mn:CdSe defect bands of Chapter 7 — is generated by exactly this two-run pattern. Internalise it once and you will read those plots as the natural output of one converged density and one frozen-density diagonalisation.

End-to-End Recipe — POSCAR to KPOINTS in Python

Putting the chapter together: given any POSCAR, build a converged KPOINTS file in one screen of Python. The script below does it for the FCC CdSe primitive cell and prints the diagnostics you should always check before launching a VASP run. Click any line to see the values flowing through it.

Build a converged KPOINTS file from a POSCAR

🐍build_kpoints.py

Explanation(26)

Code(36)

1import numpy as np

NumPy is the workhorse for every numeric step here. We will use it for matrix inversion (np.linalg.inv), transposition (.T), vector norms (np.linalg.norm), and ceiling rounding (np.ceil). All of this runs as compiled C code under the hood, so even on a triclinic supercell with 200 atoms the lattice maths takes microseconds.

EXECUTION STATE

numpy = Numerical Python library — provides ndarray (multi-dimensional array), linear algebra, broadcasting, statistical functions, and random number generation

as np = Universal alias. Means we write np.linalg.inv() instead of numpy.linalg.inv() throughout. Every VASP-related Python script you will read uses this same alias.

2from pathlib import Path

pathlib.Path is the modern Python way to handle filesystem paths. It works on Windows, macOS, and Linux without changing the code. We use it here to read the POSCAR and write the KPOINTS file in two short lines.

EXECUTION STATE

pathlib = Standard-library module (Python 3.4+) for filesystem paths

Path = Class with convenient methods like .read_text(), .write_text(), .exists(), .glob(). Path('POSCAR').read_text() returns the entire file as one string — no boilerplate open()/close().

4def read_lattice(poscar_path) → np.ndarray

Parses a VASP POSCAR header and returns the 3×3 direct lattice matrix in Å. We deliberately ignore species, counts, and atomic coordinates — for k-point generation only the lattice matters.

EXECUTION STATE

⬇ input: poscar_path = String path to a POSCAR file. Example: 'POSCAR' (current directory) or '/data/CdSe/POSCAR'.

POSCAR layout (first 5 lines) =

Line 1: comment
Line 2: scale factor (a₀)
Line 3: lattice vector a₁
Line 4: lattice vector a₂
Line 5: lattice vector a₃

⬆ returns = np.ndarray of shape (3, 3). Row i is aᵢ in Å. Always fully expanded (scale already multiplied through).

5Docstring

One-line summary of the function purpose and units. Crucial in VASP work — many bugs come from mixing Å and Bohr or forgetting the scale factor on line 2.

6text = Path(poscar_path).read_text().splitlines()

Loads the POSCAR file and splits it into a list of lines (no trailing newlines). For a 100-line file this allocates one Python list of 100 strings.

EXECUTION STATE

📚 Path(...).read_text() = Reads the entire file as a single string. Decodes using the system default encoding (UTF-8 on every modern system). For POSCAR this is fine — POSCARs are pure ASCII.

📚 .splitlines() = String method: splits on \n, \r, \r\n. Handles Windows CRLF transparently. Example: 'a\nb\r\nc'.splitlines() → ['a', 'b', 'c'].

⬆ text = List[str], one entry per line of POSCAR

→ example for CdSe primitive POSCAR = [ 'CdSe primitive', ' 6.0770000000000000', ' 0.0000000000000000 0.5000000000000000 0.5000000000000000', ' 0.5000000000000000 0.0000000000000000 0.5000000000000000', ' 0.5000000000000000 0.5000000000000000 0.0000000000000000', ... ]

7scale = float(text[1].strip())

Reads line 2 of the POSCAR (index 1 since Python is 0-indexed). This is the scale factor — every lattice row will be multiplied by it before being returned.

EXECUTION STATE

text[1] = ' 6.0770000000000000' — for CdSe zincblende, a = 6.077 Å

📚 .strip() = Removes leading/trailing whitespace (spaces, tabs, newlines). Vital because POSCARs from different programs have inconsistent indentation.

📚 float() = Parses '6.0770000' → 6.077. Throws ValueError on bad input — better than silent wrong answers.

⬆ scale = 6.077 (Å)

→ POSCAR convention = If scale > 0: it is the lattice constant in Å, multiplying the 3 vectors below. If scale < 0: |scale| is the cell volume in Å³ (rarer convention). VASP accepts both.

8rows = [list(map(float, text[2 + i].split()[:3])) for i in range(3)]

List comprehension that pulls lines 3, 4, 5 (indices 2, 3, 4), splits each on whitespace, takes the first three numbers, and converts them to floats.

EXECUTION STATE

📚 text[2 + i] = i=0 → line 3 (a₁), i=1 → line 4 (a₂), i=2 → line 5 (a₃). The three direct lattice rows.

📚 .split() = No-argument form — splits on any whitespace (spaces, tabs) and discards empty entries. ' 0.0 0.5 0.5 '.split() → ['0.0', '0.5', '0.5'].

📚 [:3] = Slice keeping the first 3 elements. Defensive — some POSCAR variants append selective-dynamics flags to lattice rows.

📚 map(float, ...) = Lazy iterator that calls float() on each string. Wrapped in list() to materialise.

⬆ rows =

[
  [0.0, 0.5, 0.5],
  [0.5, 0.0, 0.5],
  [0.5, 0.5, 0.0],
]
The FCC primitive lattice in scale-factor units (still has to be multiplied by 6.077).

9return scale * np.array(rows)

Multiplies the unit-scale lattice rows by the lattice constant and wraps everything in a NumPy array. The result is the 3×3 matrix A in Å.

EXECUTION STATE

📚 np.array(rows) = Converts a Python list-of-lists into an ndarray of shape (3, 3). Required for matrix algebra — Python lists do not support np.linalg.inv.

📚 scale * ndarray = NumPy broadcasting: the scalar 6.077 multiplies every element. Element-wise, runs in C, no Python loop.

⬆ return: A (3×3, Å) =

        a_x        a_y        a_z
a₁   0.00000   3.03850   3.03850
a₂   3.03850   0.00000   3.03850
a₃   3.03850   3.03850   0.00000

→ check: |aᵢ| = √(0² + 3.0385² + 3.0385²) = 3.0385·√2 ≈ 4.297 Å. Matches the OUTCAR readout in §3.2.

11def reciprocal(A) → np.ndarray

Computes the physicist-convention reciprocal lattice in one line. This is the matrix form of the cross-product recipe from §3.2 — equivalent but faster.

EXECUTION STATE

⬇ input: A (3×3, Å) = Direct lattice matrix from read_lattice(). Rows are aᵢ.

⬆ returns = B (3×3, Å⁻¹). Rows are bᵢ. Satisfies bᵢ·aⱼ = 2π δᵢⱼ exactly.

12Docstring — physics convention

Reminds the reader (and your future self) which convention is in force. VASP's OUTCAR drops the 2π. Solid-state textbooks include it. Mixing them is the most common bug in plane-wave codes.

13return 2.0 * np.pi * np.linalg.inv(A).T

Three operations chained: invert A, transpose, multiply by 2π. The result is exactly the reciprocal basis from §3.2's cross-product formula.

EXECUTION STATE

📚 np.linalg.inv(A) = Computes A⁻¹ via LU decomposition. Cost: O(n³) but here n=3 — about a microsecond. Throws LinAlgError if A is singular (which would mean a degenerate cell).

→ A⁻¹ for FCC primitive =

         x       y       z
row 1  -1/a   +1/a   +1/a
row 2  +1/a   -1/a   +1/a
row 3  +1/a   +1/a   -1/a
For a = 6.077: each non-zero entry = ±0.1646 Å⁻¹.

📚 .T = Matrix transpose. For a symmetric pattern like FCC primitive A⁻¹ this is a no-op, but for low-symmetry cells (monoclinic, triclinic) it matters — the transpose is what makes Bᵢⱼ·Aⱼₖ = 2π δᵢₖ.

→ why transpose? = By definition, A·B^T = 2π·I (rows of A times rows of B). That means B = 2π·(A⁻¹)^T. Skip the transpose and you accidentally compute (B^T·A) = 2π·I which is the same on a symmetric cell but wrong on a tilted one.

📚 2.0 * np.pi * = The physics convention multiplier. Drop it and you reproduce VASP's OUTCAR display (which omits 2π).

⬆ return: B (3×3, Å⁻¹) =

        bₓ        b_y        b_z
b₁   -1.0339    1.0339    1.0339
b₂    1.0339   -1.0339    1.0339
b₃    1.0339    1.0339   -1.0339

→ magnitude = |bᵢ| = (2π/a)·√3 = 1.0339·1.7321 ≈ 1.791 Å⁻¹.

15def mp_density(B, length=20.0) → (n_x, n_y, n_z)

Picks the Monkhorst–Pack mesh density along each reciprocal direction so the product nᵢ·|bᵢ| reaches a target 'length' in Å⁻¹. The 'length × |b|' rule is the de-facto industry standard for choosing KPOINTS densities (see Hinuma et al. 2017 and the Materials Project workflow).

EXECUTION STATE

⬇ input: B (3×3, Å⁻¹) = Reciprocal lattice from reciprocal()

⬇ input: length (Å⁻¹) = Convergence-target length. Default 20.0. Rule of thumb: 15-25 for routine DFT, 30-50 for tight DOS or band-edge calculations, 40+ for elastic constants.

⬆ returns = Tuple (n_x, n_y, n_z) of positive integers — the three numbers VASP wants on line 4 of KPOINTS.

16Docstring — convergence rule

Documents the decision rule. n·|b| ≥ length means you sample at least 'length' Å⁻¹ along that reciprocal direction. Larger cells get smaller |b| and therefore smaller n — exactly the inverse-volume law from §3.2.

17b_lens = np.linalg.norm(B, axis=1)

Computes the length of each reciprocal vector. axis=1 means 'reduce along columns within each row' — that is, for each row of B compute its Euclidean norm.

EXECUTION STATE

📚 np.linalg.norm = Vector/matrix norm. With axis=None it reduces the whole array to a scalar. With axis=1 on a (3,3) matrix it returns a (3,) vector — one norm per row.

⬇ arg: axis=1 = Tells NumPy to take norms row-wise. axis=0 would be column-wise. axis=None would give a single Frobenius norm.

⬆ b_lens = [1.7910, 1.7910, 1.7910] Å⁻¹

→ all equal? why? = Because FCC primitive is highly symmetric — all three reciprocal vectors are images of each other under the cubic point group. For monoclinic / triclinic cells the three numbers will differ.

18return tuple(int(np.ceil(length / bl)) for bl in b_lens)

Generator expression that turns the three lengths into three integer mesh sizes. Each nᵢ is the smallest integer with nᵢ·|bᵢ| ≥ length.

EXECUTION STATE

📚 np.ceil = Element-wise ceiling. np.ceil(11.16) → 12.0 (still a float). We wrap in int() to convert to a Python int suitable for KPOINTS.

→ arithmetic for CdSe = length / |b₁| = 20.0 / 1.7910 ≈ 11.166 ceil(11.166) = 12 Same for b₂, b₃ (cubic symmetry).

📚 tuple(...) = Materialises the generator into an immutable tuple. Tuples are the conventional return type for fixed-length VASP arguments — it signals 'three numbers, do not append'.

⬆ return: (n_x, n_y, n_z) = (12, 12, 12)

20def write_kpoints(n, path='KPOINTS')

Writes the actual KPOINTS file VASP will read. The 5-line format is the Γ-centred automatic mesh — the safest default for nearly every cubic and hexagonal system.

EXECUTION STATE

⬇ input: n = (n_x, n_y, n_z) from mp_density(). Example: (12, 12, 12).

⬇ input: path = Output filename. Defaults to 'KPOINTS' (VASP's expected name).

21Docstring — Γ-centred

Notes the centring choice. 'Gamma' (Γ-centred) puts one mesh point exactly on the zone centre. 'Monkhorst' (capital M) shifts the whole mesh so Γ is between points. Γ-centred is preferred for hexagonal/trigonal lattices and is the safer default — it never accidentally breaks a symmetry by placing two points on the same Wyckoff site.

22Path(path).write_text(...)

One-shot file write. Path.write_text creates the file (overwriting if it exists), writes the string, and closes the handle — no boilerplate.

EXECUTION STATE

📚 Path.write_text = Writes the entire string. Equivalent to: with open(path, 'w') as f: f.write(text). Encoding defaults to UTF-8, perfect for ASCII KPOINTS.

→ resulting KPOINTS file = Automatic mesh 0 Gamma 12 12 12 0 0 0

→ line-by-line meaning = Line 1: free-form comment (anything you want) Line 2: 0 means automatic generation (no manual list) Line 3: Gamma or Monkhorst-Pack centring keyword Line 4: nx ny nz mesh subdivisions Line 5: optional shift in fractional reciprocal coords (here zero)

28f-strings + escaped newlines

We use Python f-strings for n[0], n[1], n[2] interpolation and \n inside the literal for newlines. The double backslash in our source becomes a single backslash + n which Python then interprets as the newline character.

EXECUTION STATE

→ 5-line layout = After write_text completes, the file on disk is exactly 5 lines, the format VASP's k-point parser expects.

31A = read_lattice('POSCAR')

Top-level call. Loads the lattice from the current directory's POSCAR file.

EXECUTION STATE

⬆ A (3×3, Å) for CdSe primitive =

        a_x       a_y       a_z
a₁   0.0000   3.0385   3.0385
a₂   3.0385   0.0000   3.0385
a₃   3.0385   3.0385   0.0000

32B = reciprocal(A)

Computes the reciprocal lattice. Same matrix you would get by hand from the cross-product formulae in §3.2.

EXECUTION STATE

⬆ B (3×3, Å⁻¹) =

        b_x        b_y        b_z
b₁   -1.0339    1.0339    1.0339
b₂    1.0339   -1.0339    1.0339
b₃    1.0339    1.0339   -1.0339

33n = mp_density(B, length=20.0)

Decides the mesh. We pass length=20.0 — a moderate-tight density appropriate for an SCF calculation that will feed a band structure. Tighten to 25-30 for DOS plots that need smooth peaks.

EXECUTION STATE

⬆ n = (12, 12, 12)

→ why 12? = ceil(20.0 / 1.791) = ceil(11.166) = 12. The same calculation in pymatgen and ASE produces the same number (within their own length conventions).

34write_kpoints(n)

Writes the file. Side effect only — returns None. After this line, ./KPOINTS exists on disk and VASP can read it.

35print |a_i| diagnostic

Sanity check — prints each direct lattice length so you can verify against the POSCAR header you printed by eye.

EXECUTION STATE

→ expected output = |a_i| = [4.297 4.297 4.297] angstrom

36print |b_i| diagnostic

Reciprocal lengths. Multiplied by 12 (the chosen n) they give 12·1.791 ≈ 21.5 Å⁻¹, comfortably above the 20 Å⁻¹ target.

EXECUTION STATE

→ expected output = |b_i| = [1.791 1.791 1.791] 1/angstrom

37print mesh summary

Prints the chosen mesh. After this line you can visually check that VASP receives the same numbers you intended.

EXECUTION STATE

→ expected output = Mesh : (12, 12, 12)

10 lines without explanation

1import numpy as np
2from pathlib import Path
3
4def read_lattice(poscar_path: str) -> np.ndarray:
5    """Return the 3x3 direct lattice matrix in angstroms."""
6    text = Path(poscar_path).read_text().splitlines()
7    scale = float(text[1].strip())
8    rows = [list(map(float, text[2 + i].split()[:3])) for i in range(3)]
9    return scale * np.array(rows)
10
11def reciprocal(A: np.ndarray) -> np.ndarray:
12    """Physics convention: B[i] dot A[j] = 2*pi delta_ij."""
13    return 2.0 * np.pi * np.linalg.inv(A).T
14
15def mp_density(B: np.ndarray, length: float = 20.0) -> tuple[int, int, int]:
16    """Pick n_i so that n_i * |b_i| >= length (in inverse angstroms)."""
17    b_lens = np.linalg.norm(B, axis=1)
18    return tuple(int(np.ceil(length / bl)) for bl in b_lens)
19
20def write_kpoints(n: tuple[int, int, int], path: str = "KPOINTS") -> None:
21    """Gamma-centred Monkhorst-Pack mesh, no shift."""
22    Path(path).write_text(
23        "Automatic mesh\n"
24        "0\n"
25        "Gamma\n"
26        f"{n[0]} {n[1]} {n[2]}\n"
27        "0  0  0\n"
28    )
29
30A = read_lattice("POSCAR")
31B = reciprocal(A)
32n = mp_density(B, length=20.0)
33write_kpoints(n)
34print(f"|a_i| = {np.linalg.norm(A, axis=1).round(3)} angstrom")
35print(f"|b_i| = {np.linalg.norm(B, axis=1).round(3)} 1/angstrom")
36print(f"Mesh : {n}")

Running the script in the same directory as the CdSe POSCAR prints:

📝text

1|a_i| = [4.297 4.297 4.297] angstrom
2|b_i| = [1.791 1.791 1.791] 1/angstrom
3Mesh : (12, 12, 12)

and writes the five-line KPOINTS file from Mode 1 above. Three things worth taking with you from this listing:

The reciprocal lattice is one line of NumPy. B = 2*np.pi*np.linalg.inv(A).T is the entire content of §3.2.
Convergence is one rule of thumb. $n_i = \lceil \ell / |\mathbf{b}_i| \rceil$ covers the 90 % case. For the harder 10 % (defects, surfaces, metals), do an explicit $n$ -vs-energy convergence sweep — it is the first project of Chapter 6.
The KPOINTS format is unforgiving. One stray space, a missing line, or the wrong centring keyword and VASP fails with a cryptic error. Always read the file you wrote with cat KPOINTS before launching.

Practical Recipe for the Mn:CdSe Target

The end-game of this book is a Mn-doped CdSe quantum-dot calculation. Reciprocal-space settings change as the cell grows. The table below is a cheat sheet you will return to in Chapter 6.

Calculation stage	Cell	Mesh density (ℓ)	Resulting KPOINTS
Bulk SCF (primitive)	FCC primitive, a = 6.077 Å	ℓ = 20 Å⁻¹ — moderate	12 × 12 × 12 Γ-centred
Bulk DOS / band structure	Same primitive	ℓ = 30 Å⁻¹ — tight	16 × 16 × 16 SCF, then line mode
2×2×2 supercell (clean)	Conventional 64-atom cube, a ≈ 12.15 Å	ℓ = 20 Å⁻¹	6 × 6 × 6 Γ-centred
2×2×2 + 1 Mn substitution	Same supercell	ℓ = 20 Å⁻¹	4 × 4 × 4 Γ-centred (lower symmetry)
3×3×3 supercell (Mn isolation)	216-atom cube, a ≈ 18.23 Å	ℓ = 18 Å⁻¹	3 × 3 × 3 Γ-centred
Optical absorption tensor	Same as the SCF cell	ℓ = 40 Å⁻¹	Doubled in each direction

Why does the mesh get coarser as the cell grows?

Because real-space cell size and reciprocal-space mesh density are inversely related (§3.2). Doubling the supercell halves $|\mathbf{b}_i|$ , so half as many subdivisions reach the same $\ell$ target. The total number of irreducible k-points falls roughly like $1/N_{\text{cell}}^2$ for a fixed convergence target — which is the reason supercell calculations are tractable in practice.

Common Pitfalls

Pitfall	Symptom	Fix
Wrong centring for hexagonal cell	VASP warning: 'k-mesh breaks symmetry of crystal'.	Use the Gamma keyword on line 3 of KPOINTS for any non-cubic system.
Confusing Reciprocal vs Cartesian in line mode	Band structure looks plausible but the high-symmetry labels are wrong.	Use Reciprocal coordinates for line mode unless you have a very specific reason.
Missing the 2π in OUTCAR comparison	Hand-derived \|b_i\| differs from OUTCAR by exactly 2π.	OUTCAR prints b_i / (2π). Multiply by 2π before comparing to a textbook formula.
Mesh too coarse for a metal	Total energy oscillates as the mesh increases.	Metals need finer meshes (ℓ ≥ 50) and Methfessel–Paxton smearing. Insulators tolerate ℓ ≈ 20.
Re-using bulk KPOINTS for a defect supercell	Calculation runs for a week. Wrong physics.	Scale the mesh inversely with supercell size. The Python recipe above does this automatically.
Lower symmetry than expected	IBZKPT has many more irreducible points than a textbook table predicts.	Symmetry-breaking impurities (Mn substitution) reduce the point group. Check the SYMMETRY block in OUTCAR.

Summary

Every reciprocal-space concept of Chapter 3 has a textual analogue in a VASP file. POSCAR carries the real basis (and so implicitly $\mathbf{B}$ ); KPOINTS picks $\mathbf{k}$ ; IBZKPT records the symmetry reduction; EIGENVAL stores $\varepsilon_n(\mathbf{k})$ ; DOSCAR is the BZ integral.
OUTCAR's reciprocal-lattice header is $\mathbf{B} / (2\pi)$ — the crystallographer convention. Multiply by $2\pi$ to recover the physics-textbook $\mathbf{b}_i$ .
KPOINTS has four modes — automatic, line, explicit, generalised. SCF and DOS use automatic; band structures use line mode along the Setyawan–Curtarolo path.
IBZKPT integer weights add up to $n_1 n_2 n_3$ . The symmetry speedup over a brute-force calculation is the order of the point group divided by 2 — typically ×24 for cubic systems.
The convergence rule of thumb $n_i \,|\mathbf{b}_i| \geq \ell$ with $\ell \in [15, 30]\,\text{Å}^{-1}$ covers most insulator and semiconductor calculations. Metals need $\ell \geq 50$ .
One screen of NumPy converts any POSCAR into a converged KPOINTS file. The script you built here returns in Chapter 6 as the first step of every Mn:CdSe simulation.

Section 3.10 Core Insight

"A VASP run is the ASCII shadow of a Brillouin-zone integral. POSCAR fixes the lattice, KPOINTS picks the integration mesh, IBZKPT carves out the symmetry wedge, and EIGENVAL/DOSCAR is the integrand itself, sampled and summed."

Coming next: Chapter 4 — Quantum Mechanics for Solids — where we open the box VASP has been hiding behind: why does each k-point produce a column of eigenvalues in the first place? Free electrons, nearly-free electrons, tight binding, and the Kohn–Sham equations all wait on the other side.