Fermi alignment for semiconductors and insulators

For any material with a band gap — semiconductors, insulators, topological insulators — it’s natural to anchor the energy scale so the valence band maximum (VBM) sits at zero. The 0.4.0 release adds two helpers for exactly this:

  • compute_band_edges() — measures where the VBM, CBM, and gap fall on a uniform k-mesh.

  • align_to_vbm() — returns a new model with on-site energies shifted so VBM = 0.

After align_to_vbm, every downstream calculator (BulkDOS, SurfaceSpectralDensity, SurfaceGreensFunction, FermiArcMap, BandStructure) automatically sees the band edge at zero — no per-calculator flags needed.

Why this matters

The Tailwater training data is Fermi-shifted so the DFT-chosen \(E_F\) sits at 0 eV in every training sample. For non-metals this puts the band gap straddling \(E = 0\): the VBM lands just below zero and the CBM just above. That’s a reasonable convention but inconvenient when you want band-edge plots, since the natural physical reference is the band edge itself, not a mid-gap point.

align_to_vbm re-anchors a single model so the band edge IS the zero. Every eigenvalue at every k is shifted by exactly the same constant — gap-preserving, Hermiticity-preserving.

Quick start

For a non-metal like Bi2Se3 (~0.18 eV bulk gap):

from tailwater import tb_model, align_to_vbm, BulkDOS, bulk_band_structure

model = tb_model.load("wannier90_hr.hdf5")
model = align_to_vbm(model)               # VBM is now exactly at E = 0

# All downstream calculators inherit the new zero — no extra args needed:
dos = BulkDOS(model, energies=(-3, 3), k_mesh=(8, 8, 8)).run()
dos.figure.savefig("bulk_dos.png")

fig = bulk_band_structure(
    model,
    k_points=[[0, 0, 0], [0.5, 0.5, 0], [0, 0, 0]],
    k_labels=[r"$\Gamma$", "M", r"$\Gamma$"],
    e_range=(-3, 3),
)
fig.savefig("bands.png")

Inspecting the band edges first

If you’d like to see the gap before deciding to align:

from tailwater import compute_band_edges

edges = compute_band_edges(model)
# -> {"vbm": -0.139, "cbm": 0.040, "gap": 0.179, "is_metal": False}

print(f"VBM = {edges['vbm']:+.4f} eV    CBM = {edges['cbm']:+.4f} eV")
print(f"gap = {edges['gap']:.4f} eV    metal = {edges['is_metal']}")

Detection runs on a uniform 4×4×4 k-mesh by default; pass k_mesh=(8, 8, 8) (or any tuple) for a denser grid.

How VBM is identified

The heuristic is intentionally simple. Assuming the model’s existing zero is somewhere inside the gap (the Tailwater training convention), compute_band_edges diagonalises \(H(\mathbf{k})\) on the k-mesh, collects every eigenvalue across every k-point, then reports:

  • vbm = max(eigs < 0) — the negative eigenvalue closest to zero,

  • cbm = min(eigs > 0) — the positive eigenvalue closest to zero,

  • gap = cbm - vbm,

  • is_metal = gap <= 0.

This avoids needing to know how many bands are occupied — a number that’s awkward to source for Wannier-projected models, especially under spin–orbit coupling.

Metals

For a metal (bands crossing \(E = 0\)), the heuristic still returns numbers but is_metal is True and the VBM/gap aren’t physically meaningful. By default, align_to_vbm recognises this and emits a RuntimeWarning and returns the unshifted model so downstream code keeps working:

>>> model = tb_model.load("some_metal.hdf5")
>>> aligned = align_to_vbm(model)
RuntimeWarning: align_to_vbm: no clean gap around E=0 on the (4, 4, 4)
k-mesh (vbm=-0.001, cbm=+0.002, gap=0.003). Consistent with a metal
(or a non-metal whose current zero isn't in the gap). Returning
unshifted model.

Override via the if_metal argument:

if_metal=...

Behavior

"warn"

(default) emit RuntimeWarning and return unshifted model.

"raise"

raise RuntimeError — fail loudly.

"skip"

silently return the unshifted model — useful in batch processing where you’d rather log + continue without warnings.

Overriding the auto-detection

If you already know the Fermi level (from a DFT calculation, a self-consistent integration, or just preference for a different reference like the CBM or mid-gap), pass it explicitly:

# Put your chosen E_F at the new zero — shift every eigenvalue by -fermi_level:
aligned = align_to_vbm(model, fermi_level=-0.123)

This bypasses both the k-mesh detection and the metal check — it simply adds a constant offset to all on-site energies.

Worked patterns

Bulk DOS, VBM-aligned:

from tailwater import tb_model, align_to_vbm, BulkDOS

model = align_to_vbm(tb_model.load("wannier90_hr.hdf5"))
dos = BulkDOS(model, energies=(-3, 3), k_mesh=(12, 12, 12)).run()
dos.figure.savefig("bulk_dos.png")

Surface Green’s function near the band edge:

import numpy as np
from tailwater import tb_model, align_to_vbm, SurfaceGreensFunction

model = align_to_vbm(tb_model.load("wannier90_hr.hdf5"))
sgf = SurfaceGreensFunction(
    model,
    surface=np.eye(3),
    energies=np.linspace(-0.5, +0.5, 201),     # ±0.5 eV around the VBM
    k_path=[[0, 0.5, 0], [0, 0, 0], [0.333, 0.333, 0]],
    k_labels=["M", r"$\Gamma$", "K"],
    n_jobs=-1,                                 # see :doc:`performance`
).run()
sgf.figure_top.savefig("surface_top.png")

Just inspecting, no shift:

edges = compute_band_edges(model, k_mesh=(8, 8, 8))
if edges["is_metal"]:
    print("metal — skipping VBM alignment")
else:
    print(f"VBM = {edges['vbm']:+.4f} eV   "
          f"CBM = {edges['cbm']:+.4f} eV   "
          f"gap = {edges['gap']:.4f} eV")

Implementation note

The on-site shift is applied to the (0, 0, 0) block of the tbmodels HopDict. Because tbmodels includes both the stored +R matrix and its Hermitian conjugate when building \(H(\mathbf{k})\), the R=(0,0,0) block contributes to \(H(\mathbf{k})\) twice. align_to_vbm therefore adds \(\frac{1}{2}\, \text{shift} \cdot I\) to hop[(0,0,0)] so the eigenvalue shift comes out to exactly the requested value (this is verified numerically — the per-band shift std is at the 1e-14 level across random k-points).

API reference

The canonical entries live in Post-processing: bands, DOS, surface states, Fermi arcs. Reproduced here for convenience (:no-index: to avoid double-indexing):

tailwater.wannier_wizard.compute_band_edges(model_or_path: str | Model, k_mesh: Tuple[int, int, int] = (4, 4, 4))[source]

Locate VBM / CBM / gap on a uniform Monkhorst-Pack k-mesh.

Assumes the model’s current zero of energy is roughly E_F (the Tailwater training convention). Diagonalizes H(k) on every k in a k_mesh[0] x k_mesh[1] x k_mesh[2] uniform grid in fractional reciprocal coordinates, then takes:

  • VBM = the negative eigenvalue closest to zero (max of e < 0)

  • CBM = the positive eigenvalue closest to zero (min of e > 0)

  • gap = CBM - VBM

  • is_metal = (gap <= 0) — i.e. bands overlap E=0 across the mesh

Parameters:

  • model_or_path (str | tbmodels.Model) – Path to the HDF5 hr-model the API produced, or a tbmodels.Model already in memory.

  • k_mesh ((int, int, int)) – Grid density. Default (4, 4, 4) — denser meshes catch the VBM/CBM at off-symmetry k more accurately at small extra cost.

Return type:

dict

Returns:

  • dict with keys {“vbm”, “cbm”, “gap”, “is_metal”} (floats, except

  • is_metal which is bool; vbm/cbm/gap may be None in degenerate cases

  • where the spectrum has no eigenvalues on one side of zero).

tailwater.wannier_wizard.align_to_vbm(model_or_path: str | Model, k_mesh: Tuple[int, int, int] = (4, 4, 4), fermi_level: float | None = None, if_metal: str = 'warn')[source]

Return a NEW model with its on-site energies shifted so VBM = 0.

This re-anchors the energy scale so the band-edge sits exactly at zero — the natural reference for plotting / DOS / surface-state computations on semiconductors and insulators, instead of whatever DFT-chosen E_F the training data was referenced against.

Parameters:

  • model_or_path (str | tbmodels.Model) – Path to the HDF5 hr-model, or a tbmodels.Model in memory.

  • k_mesh ((int, int, int)) – k-mesh used to auto-detect the VBM (only consulted if fermi_level is None). Default (4, 4, 4).

  • fermi_level (float, optional) – If supplied, bypass auto-detection and shift on-site energies by -fermi_level (i.e. put your chosen Fermi value at the new zero). Useful if you already know E_F from a DFT calculation.

  • if_metal ({"warn", "raise", "skip"}) – What to do when compute_band_edges reports the spectrum has no clean gap around E=0 (signature of a metal, or of a non-metal whose current zero isn’t in the gap). Default "warn": emits a RuntimeWarning and returns the unshifted model so downstream code still runs. "raise" errors out; "skip" silently returns the unshifted model.

Returns:

A deep copy of the input model with its (0,0,0) hop block adjusted by shift * I so every eigenvalue at every k is offset by shift = -VBM. The input model is not mutated.

Return type:

tbmodels.Model