Goal: get silicon's actual band gap — 0.61 eV, not the 2.28 eV tutorial 02 read off an SCF grid — and understand why the SCF grid could never have given it.
An SCF calculation samples the Brillouin zone wherever the integrals need sampling. It converges the density, and the density is an integral over occupied states; a Monkhorst-Pack grid is a quadrature rule chosen for that job. Nothing about it guarantees a point sits at the valence-band maximum or the conduction-band minimum.
For silicon those two extrema are at different places, and one of them is at no rational grid point you would pick by accident: the VBM is at Γ, and the CBM is about 83% of the way from Γ to X. Tutorial 02's 2×2×2 grid held two inequivalent points, neither of them near that minimum, so the smallest gap it could see was 2.28 eV.
The fix is to decouple the two jobs. Converge the density on a grid, then stop iterating and diagonalize the resulting fixed Hamiltonian at whatever k-points you actually care about. That second step is a non-self-consistent field (NSCF) calculation, and it is cheap: one diagonalization per k-point, no loop.
The SCF step is an ordinary calculation — the only requirements are that it be well converged (tutorial 03) and that you remember its outdir and prefix, because that is how the NSCF step finds it. Write scf.yml:
Now nscf.yml. Same material, same method — a different Hamiltonian would make the reloaded potential meaningless — but calculation: nscf, no kgrid, and an explicit kpoints_list of [kx, ky, kz, weight] rows. The coordinates are crystal — fractions of the reciprocal lattice vectors, the same convention a generated kgrid uses. These four are Γ, the point 83% along Γ–X, X itself, and L:
That first line is the whole mechanism. The converged effective potential — the .save/ contents tutorial 02 listed — is read back instead of being recomputed, and the four k-points are diagonalized against it.
The NSCF run overwrites <prefix>_results.json. Same outdir, same prefix, same filename. The SCF energies in it are replaced by NSCF eigenvalues at the new k-points. Copy the SCF results aside first if you need both.
Reading the eigenvalues out — silicon has 8 valence electrons, so bands 0–3 are occupied and band 4 is the first empty one:
0.607 eV, from four k-points. Every one of those direct gaps — the vertical distance at a single k-point — is above 2.5 eV. The fundamental gap is small only because the maximum and the minimum are at different k, which is precisely the information a Monkhorst-Pack grid throws away.
You can also use more bands in the NSCF step than the SCF used (number_of_bands: 16 works against a 12-band SCF): the extra conduction bands never entered the density, so they cost one diagonalization and nothing else.
Four points prove the mechanism; a plot needs the path filled in. That bookkeeping — looking up high-symmetry coordinates for the lattice, interpolating segments, tracking cumulative distance for the x-axis — is what q3.workflows.band_structure does. It needs the Python package (tutorial 05 covers the install):
It runs both steps for you and returns eigenvalues_ry, kpoints, kpath_distances (the plot's x-axis) and fermi_energy_ry. On the same analysis as above the 151-point path gives:
| Quantity | Value |
|---|---|
| Indirect gap | 0.607 eV |
| Valence band maximum | Γ |
| Conduction minimum | (0.4167, 0, 0.4167) — 83% of Γ–X |
| Direct gap at Γ | 2.547 eV |
The same 0.607 eV as the four hand-picked points, which is the reassuring outcome: the path did not find anything lower, so those four points were the right four. Silicon's experimental gap is 1.17 eV — PBE underestimates it, as semilocal functionals do, and 0.6 eV is the correct PBE answer. Getting 1.17 eV out of this calculation would mean something was wrong, not right.
Two committed manifests do not run through the binary directly. examples/band_structure--*/nc.yml and examples/density_of_states--*/nc.yml declare calculation: band_structure and calculation: density_of_states, which the YAML parser accepts but the flavor factory does not implement — ./build/q3 on either exits with unsupported calculation type. They are inputs for the Python workflows above, which rewrite the calculation type to self_consistent_field before invoking the binary.
The DOS asks a different question: not where the states are in k, but how many there are at each energy. It comes from the same NSCF machinery — a dense grid instead of a path, then a broadened histogram of the eigenvalues.
The check worth running on any DOS — integrate it up to the Fermi level and see whether you get your electrons back:
Eight, as it must be for two silicon atoms at Z=4. This catches a broken normalization, a wrong Fermi level, or a k-point weight bug in one line — and it is the reason to trust anything else the DOS tells you.
Do not read the band gap off this plot. With broadening_ry: 0.01 — 0.14 eV — the Gaussian tails from the band edges fill the gap in, and the width you measure depends entirely on the threshold you call "zero": 0.27 eV at 2% of peak height, 0.45 eV at 5%. Neither is the 0.607 eV that the band path gives. Broadening is a plotting parameter; the gap is not a plotting parameter. Get gaps from eigenvalues on a path, and use the DOS for what it is good at — where the weight is, how it is distributed, whether there is a peak at the Fermi level.
05 — Python and ASE: installing the package used above, and what else it drives — relaxation, phonons, and the ASE calculator.
1.9.8 