Surface charge density

Real-space surface charge-density heat maps of a general (hkl) slab, built directly from a Wannier Hamiltonian’s real-space H(R). Works on any Wannier tight-binding model — a Tailwater prediction or a DFT-generated Wannier90 Hamiltonian — because the only interchange format it needs is H(R).

Pipeline: re-express H(R) in an integer supercell whose first two lattice vectors lie in the (hkl) plane (an exact, determinant-preserving remap), stack size cells along the surface normal and drop hoppings that leave the slab, integrate \(|\psi|^2\) of the occupied states over the surface BZ to get a per-orbital occupation, then render \(\rho(\mathbf{r}) = \sum_g n_g\,\mathcal{G}(\mathbf{r}-\mathbf{r}_g)\) with the Wannier centres as \(\mathbf{r}_g\).

Quick example

from tailwater import surface_charge_density, load_hr, supercell_self_check

# `model` accepts a tbmodels.Model, an HDF5 path, a Wannier90 *_hr.dat
# (DFT output), or the dict returned by load_hr().
HR_PATH = "outputs/wannier90_hr.hdf5"
MILLER  = (0, 0, 1)     # surface Miller index
SIZE    = 4             # slab thickness in unit cells

# Sanity-gate the general-(hkl) supercell remap (expect ~1e-13 eV).
model = load_hr(HR_PATH)
assert supercell_self_check(model, MILLER) < 1e-8

# Top-view + side cross-section heat maps; everything past `size` is optional.
res = surface_charge_density(
    model, MILLER, SIZE,
    mu=0.0, nk=12, sigma=0.6, tile=3,
    savepath="surface_charge_001.png",
)
rho, top_img, side_img = res["rho"], res["top_img"], res["side_img"]

Image a topological surface state by restricting the occupation to a narrow window around \(E_F\):

surface_charge_density(model, MILLER, SIZE, energy_window=(-0.1, 0.1),
                       savepath="surface_charge_001_tss.png")

A DFT Wannier90 model drops in unchanged — pass the *_hr.dat path directly:

surface_charge_density("path/to/wannier90_hr.dat", (1, 1, 1), 5)

See examples/11_surface_charge_density.py for the full runnable script.

Note

The per-k diagonalisation loop can stall under some OpenBLAS builds. The function limits BLAS threads via threadpoolctl when it is installed; otherwise set OMP_NUM_THREADS=1 in the environment.

API

tailwater.surface_charge.surface_charge_density(model, miller, size, *, mu=0.0, nk=12, sigma=0.6, tile=3, ngrid=320, surface_thickness=None, energy_window=None, cmap='turbo', show=True, savepath=None, title=None)[source]

Render real-space surface charge-density heat maps of a (hkl) slab.

Parameters:

  • model (tbmodels.Model | str | pathlib.Path | dict) – The Wannier Hamiltonian. A tbmodels.Model (e.g. from a Tailwater prediction loaded via tb_model.load), a path to a tbmodels HDF5 or a Wannier90 *_hr.dat (DFT output), or the internal dict from load_hr().

  • miller ((int, int, int)) – Surface Miller index, e.g. (0, 0, 1) or (1, 1, 1).

  • size (int) – Slab thickness in unit cells (cells stacked along the surface normal).

  • mu (float, default 0.0) – Fermi level (eV). Occupied states are those with E < mu (ignored when energy_window is given). Default 0 matches the Tailwater training convention.

  • nk (int, default 12) – Surface-BZ Monkhorst-Pack mesh is nk × nk.

  • sigma (float, default 0.6) – Gaussian radius (Å) used to render each Wannier centre.

  • tile (int, default 3) – Number of in-plane unit-cell repetitions in the top view.

  • ngrid (int, default 320) – Pixels along the long axis of each heat map.

  • surface_thickness (float, optional) – Depth (Å) from the top surface counted as “the surface” for the top view. Defaults to ~0.6 of one layer’s thickness.

  • energy_window ((float, float), optional) – If given, occupy states with emin < E < emax instead of E < mu. Use a small window around E_F (e.g. (-0.1, 0.1)) to image topological surface states.

  • cmap (str, default "turbo") – Matplotlib colormap.

  • show (bool, default True) – Display the figure (plt.show()).

  • savepath (str | pathlib.Path, optional) – If given, save the figure to this path (PNG by default).

  • title (str, optional) – Figure supertitle; defaults to a description built from the inputs.

Returns:

  • dict with keys (rho (per-orbital occupation), pos, depth,)

  • layer, nhat, top_img, side_img, top_extent,

  • side_extent, slab, sup, and fig (the matplotlib Figure,

  • or None when show is False and no figure was created).

tailwater.surface_charge.load_hr(path: str | Path)[source]

Load a Wannier Hamiltonian from disk into the internal model dict.

.hdf5/.h5 -> tbmodels.Model.from_hdf5_file. Anything else (e.g. wannier90_hr.dat) -> from_wannier_files.

Parameters:

path (Union[str, Path])

Return type:

dict

tailwater.surface_charge.supercell_self_check(model: dict, miller: Sequence[int], ntrials: int = 20)[source]

Max band-edge mismatch (eV) between bulk and supercell at the same physical k.

A det-preserving re-basis must reproduce the bulk spectrum exactly; this returns the worst band-edge discrepancy over ntrials random k-points (expect ~1e-13 eV). Useful as a sanity gate before trusting a slab.

Parameters:
Return type:

float