Exporting models: sparse .npz, HDF5, _hr.dat, pybinding, PythTB, Kwant

By default the API returns the tight-binding Hamiltonian in sparse form — a wannier90_hr.npz (a tailwater.SparseHR) that stores only the non-zero hoppings, so it is O(N) in memory and file size rather than O(N²). tailwater.tw_api_call() handles the two size regimes for you:

• Small systems (< 30 atoms) — the .npz is automatically converted to a dense tbmodels HDF5 on your machine and returned under paths["hdf5"] (the .npz is kept too, under paths["npz"]), so every dense recipe on this page — tb_model.load("...hdf5") onward — works unchanged.

• Large systems — the result stays sparse (paths["npz"]), and a note is printed pointing at the conversions here. For big cells you should keep it sparse (see Why maintaining sparsity matters for large systems).

Force either format with tw_api_call(..., output_format="hdf5") or "sparse".

Importantly, the conversion functions below accept either a sparse .npz / SparseHR or a dense HDF5 / tbmodels.Model, so the same one-liner works whatever you’re holding.

The sparse .npz format (SparseHR)

Load a .npz with tailwater.SparseHR:

from tailwater import SparseHR

shr = SparseHR.load("wannier90_hr.npz")
shr.num_wann      # number of Wannier orbitals (the H(k) matrix dimension)
shr.nnz           # number of stored hoppings

Internally the .npz holds the Hamiltonian as a COO sparse list of hoppings plus the on-site diagonal and (optionally) the geometry — a dense [num_wann, num_wann] matrix is never formed:

on_site	real on-site energy per orbital, shape [num_wann]
rows, cols	orbital indices i, j of each stored hopping
Rs	lattice vector R per hopping, shape [nnz, 3]
vals	complex hopping amplitude H_ij(R)
cell	3×3 lattice vectors (Å), when geometry was recovered
positions	per-orbital Cartesian positions, when available

Only the forward half of each ±R pair is stored (the Hermitian conjugate at -R is implied) and the R = 0 diagonal lives in on_site. That is what makes it O(N): a .npz for a 14,000-orbital moiré cell is a few MB, where the dense HDF5 would be tens of GB.

Convert the .npz to any format (one call, auto-detecting the input)

The top-level converters accept a SparseHR / .npz path or a tbmodels.Model / .hdf5 / _hr.dat and dispatch automatically:

from tailwater import (as_tbmodels, to_hdf5, to_hr_dat,
                       to_pb, to_pythtb, to_kwant)

npz = "wannier90_hr.npz"

model      = as_tbmodels(npz)               # tbmodels.Model (dense)
to_hdf5(npz,   "wannier90_hr.hdf5")         # tbmodels HDF5
to_hr_dat(npz, "wannier90_hr.dat")          # Wannier90 _hr.dat
pb_lattice = to_pb(npz)                      # pybinding.Lattice
py_model   = to_pythtb(npz)                 # pythtb model
syst, lat  = to_kwant(npz)                  # kwant (Builder, lattice)

# the identical calls work on a dense HDF5 / _hr.dat, too:
to_hr_dat("wannier90_hr.hdf5", "wannier90_hr.dat")

as_tbmodels(npz) is the bridge to the rest of this page: once you have the tbmodels.Model (or the auto-converted HDF5 for a small system), every pybinding / PythTB / Kwant / WannierBerri recipe below applies.

Note

_hr.dat and HDF5 are dense on-disk formats — size grows as num_R · num_wann². They are guarded for very large systems; pass max_wann= to SparseHR.to_hr_dat() / SparseHR.to_hdf5() to override the guard if you really intend to write a huge file.

Staying sparse: pybinding, Kwant, and built-in solvers

For large systems you usually do not want to densify at all. SparseHR builds pybinding and Kwant models straight from the COO list (no dense matrix is ever formed) and carries its own sparse solvers:

from tailwater import SparseHR

shr = SparseHR.load("wannier90_hr.npz")

# --- built-in sparse spectra (large num_wann OK) ---
Hk = shr.Hk([0.0, 0.0, 0.0])                   # scipy sparse H(k) at Γ
w  = shr.eigsh_near_fermi([0, 0, 0], e_fermi=0.0, num=20)  # 20 states near E_F
Rd = shr.hr_dict()                             # {R: scipy.sparse.csr_matrix}

# --- hand the sparse model to pybinding / Kwant (built from the COO) ---
pb_lattice = shr.to_pb()                       # pybinding.Lattice
syst, lat  = shr.to_kwant()                    # kwant (Builder, lattice)

eigsh_near_fermi uses a shift-invert sparse eigensolver, so you can get just the handful of bands nearest the Fermi level for a Hamiltonian far larger than a dense H(k) could hold. hr_dict returns the H(R) blocks as scipy sparse matrices to feed your own KPM / Green’s-function / transport code. to_pb and to_kwant are sparse-native and scale to large num_wann — pybinding and Kwant are themselves sparse solvers, so the whole pipeline stays O(N).

Why maintaining sparsity matters for large systems

A dense Hamiltonian stores every [num_wann, num_wann] block for every lattice vector R: memory and file size scale as O(num_R · num_wann²). But a physical tight-binding Hamiltonian is sparse — each orbital hops only to a bounded number of neighbours — so the number of non-zeros scales just as O(num_wann). For large cells the difference is decisive:

• A twisted-bilayer moiré cell with ~14,000 orbitals is a few MB as a .npz but tens of GB as a dense _hr.dat / HDF5 — the dense form often cannot be written, let alone loaded into RAM.

• Diagonalising a dense H(k) is O(num_wann³) and needs the entire matrix resident; the sparse shift-invert path (SparseHR.eigsh_near_fermi()) touches only the non-zeros and returns just the near-Fermi bands you ask for.

• pybinding and Kwant are sparse-native, so converting via SparseHR.to_pb() / SparseHR.to_kwant() keeps memory and compute O(N) end to end.

Rule of thumb: for small systems, let tw_api_call convert to HDF5 and use the dense recipes below. For large systems, keep the ``.npz`` sparse — write _hr.dat / HDF5 only if an external tool demands it, and prefer the sparse pybinding / Kwant / built-in solvers above.

Working with the dense tbmodels.Model

The sections below operate on a dense tbmodels.Model — the HDF5 you get for a small system (tb_model.load("wannier90_hr.hdf5")), or as_tbmodels("wannier90_hr.npz") for a converted sparse model. From a dense model there are four common downstream needs, each a one-liner:

1. Write the model to a Wannier90-style _hr.dat file — so it can be consumed by external tools (Z2Pack, WannierTools, downstream DFT pipelines, custom analysis scripts that expect the plain-text Wannier90 format).

2. Convert the model to a pybinding Lattice — so you can use pybinding’s solvers, KPM routines, eigenvalue plotters, and transport tools on top of the Tailwater Hamiltonian.

3. Convert the model to a PythTB tb_model — so you can use PythTB’s band-path helpers, slab/wire builders, Berry-phase / Wannier-charge-centre routines, and the body of literature that targets PythTB.

4. Convert the model to a Kwant Builder — so you can use Kwant’s transport machinery (leads, smatrix, greens_function), the wraparound trick for bulk H(k), or any of Kwant’s sample-builder utilities.

Writing an _hr.dat file

Two equivalent entry points:

from tailwater import tb_model, write_hr_output

model = tb_model.load("wannier90_hr.hdf5")

# Option A — tbmodels' native method, attached to the loaded model:
model.to_hr_file("wannier90_hr.dat")

# Option B — Tailwater's thin wrapper, useful when the format is a
# runtime choice rather than hard-coded:
write_hr_output(model, "wannier90_hr.dat", fmt="hr_dat")
write_hr_output(model, "wannier90_hr.hdf5", fmt="hdf5")     # re-emit HDF5

Both produce the standard Wannier90 column layout:

Rx Ry Rz   i  j   Re(H)   Im(H)

with the unit-cell weights written in the same convention Wannier90 uses, so the result drops into any tool that already reads wannier90_hr.dat. The HDF5 round-trip is bit-identical.

No optional dependencies are needed for the hr export — it’s pure tbmodels under the hood.

Converting to a pybinding Lattice

Every model returned by tailwater.tb_model.load() carries a .to_pb() instance method:

from tailwater import tb_model, k_cart_from_frac
import pybinding as pb

model = tb_model.load("wannier90_hr.hdf5")
lat   = model.to_pb()

# Build a pybinding model and sample H(k) at any fractional k:
pmod = pb.Model(lat, pb.translational_symmetry())
pmod.set_wave_vector(k_cart_from_frac([0.0, 0.0, 0.0], model.uc))    # Gamma
eig  = np.sort(np.linalg.eigvalsh(pmod.hamiltonian.todense()))

These eigenvalues match np.sort(np.linalg.eigvalsh(model.hamilton([0,0,0]))) to float32 precision (~1e-6 eV).

Important

Pybinding’s pb.Model.set_wave_vector() expects k in Cartesian (rad/length) — not fractional. Use tailwater.k_cart_from_frac() to convert from the fractional convention tbmodels.Model.hamilton uses. Passing fractional k directly to set_wave_vector is the most common source of “the pybinding bands don’t match the tbmodels bands” reports.

What .to_pb() does internally:

• Reads the on-site energies off the diagonal of the R = (0, 0, 0) hop block (doubling them, see “Convention notes” below).

• Adds one Sublattice per Wannier orbital, with the position converted from fractional → Cartesian via pos_cart = pos_frac @ LM.

• Iterates the hopping dict and adds each (R, i, j) to lat.add_one_hopping. pybinding implies the Hermitian conjugate automatically, so duplicates returned by tbmodels are silently skipped.

Computing a full band structure

import numpy as np
from tailwater import tb_model, k_cart_from_frac
import pybinding as pb

model = tb_model.load("wannier90_hr.hdf5")
pmod  = pb.Model(model.to_pb(), pb.translational_symmetry())

# Gamma -> M -> K -> Gamma on the (k_x, k_y) plane
k_frac_path = np.array([
    [0.000, 0.000, 0],
    [0.500, 0.000, 0],
    [0.333, 0.333, 0],
    [0.000, 0.000, 0],
])

bands = []
for kf in k_frac_path:
    pmod.set_wave_vector(k_cart_from_frac(kf, model.uc))
    bands.append(np.sort(np.linalg.eigvalsh(pmod.hamiltonian.todense())))
bands = np.array(bands)                # shape (Npts, num_wann)

Comparing this against [np.sort(np.linalg.eigvalsh(model.hamilton(kf))) for kf in k_frac_path] gives matching curves to ~1e-6 eV — both routes diagonalise the same Hamiltonian, just with different per-orbital phase conventions on the eigenvectors.

Convention notes

These are the two conventions to_pb quietly handles for you; the CHANGELOG covers the same story for upgraders from 0.4.2 or earlier.

• On-site doubling. tbmodels.Model.hamilton(k) constructs the Hamiltonian via Σ_R stored[R] e^{ikR} followed by H += H.c. This H.c. step supplies the missing -R half for R ≠ 0, but at R=(0,0,0) it doubles the stored block. tbmodels therefore stores exactly half the physical on-site block. to_pb multiplies hop[(0,0,0)] by 2 before feeding it to pybinding so pybinding’s H(k) recovers the full physical on-site contribution.

• Position basis. tbmodels.Model.pos is in fractional coordinates; pybinding’s add_one_sublattice expects Cartesian. to_pb does the conversion pos_cart = pos_frac @ LM so the resulting pybinding Lattice has physically meaningful positions (its Brillouin-zone routines, real-space LDOS plotters, etc. all see the right geometry).

Overriding the lattice vectors

API-produced HDF5 files don’t always carry the unit cell as model.uc (it’s deliberately left as None so the same model works under any choice of lattice convention). .to_pb() falls back to the identity np.eye(3) in that case, which is fine for algebraic work but wrong for any plot in physical k-coordinates.

Pass the real lattice explicitly:

import numpy as np

a, c = 4.143, 28.636            # Bi2Se3 hexagonal in Å
LM = np.array([
    [ a,             0,  0],
    [-a/2,  a*np.sqrt(3)/2, 0],
    [ 0,             0,  c],
])

lat = model.to_pb(lattice_vectors=LM)

After this, the pybinding brillouin_zone() returns the correct hexagonal BZ and any band plot is in physical Å^-1.

Filtering tiny hops

For exploratory analysis or visualisation you may want to drop the smallest hops to keep lat lightweight:

lat = model.to_pb(hop_threshold=1e-6)         # drop |H_ij(R)| < 1e-6 eV

Keep hop_threshold low — it’s there to clean up sparse-storage zeros, not to act as a physical cutoff. The right place to set the physical threshold is the inference / hr-build step (default 0.01 eV), not here.

Converting to a PythTB tb_model

Every model returned by tailwater.tb_model.load() also carries a .to_pythtb() instance method:

from tailwater import tb_model

model    = tb_model.load("wannier90_hr.hdf5")
py_model = model.to_pythtb()

# Sample H(k) at any fractional k:
eig_gamma = py_model.solve_one([0.0, 0.0, 0.0])      # Γ
eig_m     = py_model.solve_one([0.5, 0.0, 0.0])      # M (hex zone)

These eigenvalues match np.linalg.eigvalsh(model.hamilton([0,0,0])) to float64 precision (~5×10⁻¹⁴ eV) — a much tighter agreement than the pybinding path, which is float32 (~10⁻⁶ eV).

The PythTB path is generally the easier of the two:

• PythTB takes orbital positions in fractional coordinates, the same convention tbmodels.Model.pos uses — no Cartesian conversion needed.

• PythTB’s solve_one(k) accepts fractional k directly — no analogue of tailwater.k_cart_from_frac() is needed.

• PythTB ships rich first-class helpers for band paths, slabs/wires (cut_piece), supercells (make_supercell), and Wannier-centre / Berry-phase analyses.

Computing a band structure with PythTB’s built-in helper

import numpy as np
from tailwater import tb_model
import matplotlib.pyplot as plt

model    = tb_model.load("wannier90_hr.hdf5")
py_model = model.to_pythtb()

# PythTB does the path interpolation for you:
k_path, k_dist, k_node = py_model.k_path(
    [[0,0,0], [0.5,0,0], [0.333,0.333,0], [0,0,0]],
    nk=101, report=False,
)
bands = py_model.solve_all(k_path)                    # (num_wann, nk)

fig, ax = plt.subplots()
for band in bands:
    ax.plot(k_dist, band, lw=0.7, color="k")
ax.set_xticks(k_node, [r"$\Gamma$", "M", "K", r"$\Gamma$"])
ax.set_ylabel("E (eV)")
fig.savefig("bands_pythtb.png", dpi=150)

Slabs and wires

PythTB’s cut_piece makes a 1D / 2D slab from the bulk model. For a 6-layer Bi₂Se₃ slab terminated along the c-axis:

py_slab = py_model.cut_piece(num=6, fin_dir=2, glue_edgs=False)
print(py_slab.get_num_orbitals())                 # 6 * 124 = 744

The resulting model is 2D-periodic (in-plane) and 0D along the surface-normal direction — solve it the same way:

eig_2d = py_slab.solve_one([0.0, 0.0])            # Γ of the surface BZ

Both model.to_pythtb() and model.to_pb() produce models with identical bulk Hamiltonians; pick whichever ecosystem (PythTB, pybinding, or both) fits your downstream analysis.

Converting to a Kwant Builder

Every model returned by tailwater.tb_model.load() also carries a .to_kwant() instance method:

import numpy as np, kwant
from tailwater import tb_model

model   = tb_model.load("wannier90_hr.hdf5")
builder = model.to_kwant()              # kwant.Builder, 3D periodic

# For bulk H(k) sampling, wrap-around and finalise:
syst = kwant.wraparound.wraparound(builder).finalized()

# Kwant's wraparound takes 2π·k_frac as k_x, k_y, k_z (the
# per-cell Bloch phase) — NOT Cartesian rad/length like pybinding.
k_frac = np.array([0.5, 0.0, 0.0])
phase  = 2 * np.pi * k_frac
H      = syst.hamiltonian_submatrix(
    params=dict(k_x=phase[0], k_y=phase[1], k_z=phase[2]),
)
eig    = np.sort(np.linalg.eigvalsh(H))

This matches np.linalg.eigvalsh(model.hamilton([0.5, 0, 0])) to float64 precision (~10⁻¹³ eV) — Kwant is double-precision internally.

Important

Kwant’s kwant.wraparound is the one place in this ecosystem where k is neither fractional nor Cartesian rad/length. The k_x/k_y/k_z parameters are 2π · k_frac — the per-cell Bloch phase, independent of cell size. (Use tailwater.k_cart_from_frac() only for pybinding’s set_wave_vector.) Passing a Cartesian k to Kwant is the most common source of “the Kwant bands don’t match the tbmodels bands” reports.

What .to_kwant() returns

It returns an unfinalised kwant.Builder with a 3D kwant.TranslationalSymmetry, one site per Wannier orbital. From there the user has two natural next steps:

• Bulk H(k) — wrap the Builder via kwant.wraparound.wraparound (as in the snippet above) and finalise.

• Transport / scattering — build a finite scattering region on top of the bulk Builder, attach leads via kwant.Builder(kwant.TranslationalSymmetry(...)) for each lead, and call kwant.smatrix / kwant.greens_function.

  A complete end-to-end recipe (quantum wire of a Bi₂Se₃ hr-model + two semi-infinite leads + G(E) sweep) lives at examples/05_kwant_scattering.py. The full Kwant tutorial is at https://kwant-project.org/doc/.

  One Kwant API quirk to know: bulk.symmetry.periods is a custom Kwant array class that doesn’t support direct integer indexing. Cast to a plain NumPy array first:

  periods = np.asarray(bulk.symmetry.periods)        # (3, 3)
  sym_x   = kwant.TranslationalSymmetry(periods[0])  # transport axis
  

Computing a band structure with Kwant

import numpy as np, kwant
import matplotlib.pyplot as plt
from tailwater import tb_model

model = tb_model.load("wannier90_hr.hdf5")
syst  = kwant.wraparound.wraparound(model.to_kwant()).finalized()

# Gamma -> M -> K -> Gamma  (linearly sampled fractional path)
nodes = np.array([[0,0,0], [0.5,0,0], [0.333,0.333,0], [0,0,0]])
path  = np.vstack([np.linspace(nodes[i], nodes[i+1], 33)
                   for i in range(len(nodes) - 1)])

bands = []
for k_frac in path:
    phase = 2 * np.pi * k_frac
    H = syst.hamiltonian_submatrix(
        params=dict(k_x=phase[0], k_y=phase[1], k_z=phase[2]),
    )
    bands.append(np.sort(np.linalg.eigvalsh(H)))
bands = np.array(bands)                       # (Npts, num_wann)

fig, ax = plt.subplots()
ax.plot(bands, lw=0.7, color="k")
ax.set_ylabel("E (eV)")
fig.savefig("bands_kwant.png", dpi=150)

The bands are identical to those from model.hamilton(k_frac) or model.to_pythtb().solve_one(k_frac) at every k — Kwant just gives you the rest of the transport ecosystem for free.

Using the model with WannierBerri

WannierBerri (Berry-curvature integrator for transport quantities, optical responses, Wannier-charge centres, etc.) reads tbmodels.Model directly — no model.to_X() method needed. Just pass the loaded model:

import wannierberri as wb
from tailwater import tb_model

model  = tb_model.load("wannier90_hr.hdf5")
sys_wb = wb.system.System_R.from_tbmodels(model, berry=True)

grid = wb.Grid(sys_wb, NK=(8, 8, 8), NKFFT=(4, 4, 4))

Efermi = np.linspace(-3, 3, 61)
result = wb.run(
    sys_wb, grid=grid,
    calculators={
        "dos":   wb.calculators.static.DOS           (Efermi=Efermi, tetra=True),
        "ohmic": wb.calculators.static.Ohmic_FermiSea(Efermi=Efermi),
        "ahc":   wb.calculators.static.AHC           (Efermi=Efermi),
    },
    parallel=False, symmetrize=False, dump_results=False,
)

berry=True tells WannierBerri to construct the position matrix elements <0,n|r|R,m> from model.pos — these are what every Berry-curvature-derived calculator (AHC, optical conductivity, Wannier charge centres) needs.

A complete worked recipe with DOS, longitudinal conductivity, and AHC vs. Fermi energy is at examples/06_wannierberri_conductivity.py.

Spin Hall conductivity (SHC and friends) need an extra ingredient: the spin matrix elements <n,0 | S^alpha | m,R>. WannierBerri’s from_tbmodels(..., spin=True) is documented but doesn’t actually populate these — tbmodels carries only the Hamiltonian. Tailwater synthesises them from the known sigma_z-eigenstate structure of the 18-orbital Wannier basis via tailwater.wb_system_with_spin():

import numpy as np, wannierberri as wb
from tailwater import tb_model, wb_system_with_spin

model = tb_model.load("wannier90_hr.hdf5")
sys   = wb_system_with_spin(model)        # SS_R populated; berry=True too

Efermi = np.linspace(-2.0, 2.0, 41)
grid   = wb.Grid(sys, NK=(8, 8, 8), NKFFT=(4, 4, 4))
result = wb.run(
    sys, grid=grid,
    calculators={
        "shc": wb.calculators.static.SHC(
            Efermi=Efermi,
            kwargs_formula={"spin_current_type": "simple"},
        ),
    },
    parallel=False, symmetrize=False, dump_results=False,
)
shc = np.asarray(result.results["shc"].data)   # (Nef, 3, 3, 3) in (ℏ/e)·S/cm

The function infers the sigma_z eigenstate doublets from the model’s geometric structure by default — it walks the atomic positions to group Wannier functions by atom, then pairs consecutive Kramers partners and verifies their on-site energies match. A worked end-to-end recipe with Bi₂Se₃ lives at examples/07_spin_hall_conductivity.py; the in-gap plateau in \(\sigma^{z}_{xy}\) is the topological signature.

WannierBerri brings its own optional dependencies — install them explicitly when needed:

pip install wannierberri numba   # numba: tetrahedron integration

Round-trip: HDF5 → pybinding → HDF5

A pybinding Lattice is not directly serialisable to HDF5, but you don’t need to round-trip through pybinding — just keep the original tbmodels.Model alongside:

from tailwater import tb_model, write_hr_output

model = tb_model.load("wannier90_hr.hdf5")
lat   = model.to_pb(lattice_vectors=LM)

# ... do pybinding work with `lat` ...

# When you want to persist or share the Hamiltonian, write the
# original tbmodels model — same content, two interchangeable
# serialisations:
write_hr_output(model, "Bi2Se3.hdf5",  fmt="hdf5")
write_hr_output(model, "Bi2Se3_hr.dat", fmt="hr_dat")

Building a model from raw head predictions (advanced)

The standard workflow uses the HDF5 the API ships back. If you instead need to assemble the tight-binding model on the client side from the raw dense head outputs — e.g. to experiment with a different hop_threshold — use build_hr_model_fast():

from tailwater import build_hr_model_fast, write_hr_output

hr_model = build_hr_model_fast(
    edge_pred   = edge_pred,           # [num_edges, 18, 18, 2] from API/heads
    onsite_pred = onsite_pred,         # [num_atoms, 18, 18, 2]
    gdata       = gdata,               # PyG Data the model consumed
    LM          = lattice_matrix,      # 3x3 real lattice (Å)
    atoms       = [(sym, xyz), ...],   # Cartesian per-atom positions
    hop_threshold = 0.01,              # drop |H_ij(R)| <= this (eV)
)

write_hr_output(hr_model, "wannier90_hr.dat", fmt="hr_dat")

build_hr_model_fast is byte-identical to build_hr_model at ~100-300× the speed; prefer it unless you’re debugging a build discrepancy. Both require pip install pybinding-dev — see Installation.

API reference

Sparse Hamiltonian + format-detecting converters

class tailwater.SparseHR(num_wann, on_site, rows, cols, Rs, vals, cell=None, positions=None)

Bases: object

Self-contained sparse Wannier H(R): COO hops + on-site + (optional) geometry (primitive lattice cell and per-orbital positions). Carries everything needed to reload the inference output straight into scipy / tbmodels / pybinding / Kwant — no crystal graph or GNN required.

Construct one with load() (from an API-produced .npz); the per-hop arrays are the forward-only half of each +/-R pair (the H.c. is implied), with the R=0 diagonal carried separately in on_site.

num_wann

on_site

rows

cols

Rs

vals

cell

positions

property nnz

Number of stored (forward-only) hops.

save(path)

Write the full sparse model (hops + on-site + geometry) to .npz.

classmethod load(path)

Reload a SparseHR from an API wannier90_hr.npz (or one written by save()) — no graph/model needed.

Hk(k)

Sparse H(k) matching tbmodels.hamilton (convention 2). Each stored hop (val, i, j, R) contributes e^{2*pi*i*k.R}*val at (i, j) and its conjugate at (j, i); on_site fills the diagonal. Vectorised -> O(nnz).

eigvals_grid(kpts, dense=True)

Dense eigenvalues (sorted ascending) at each fractional k in kpts. Returns an array of shape [len(kpts), num_wann].

eigsh_near_fermi(k, e_fermi=0.0, num=40)

The num eigenvalues nearest e_fermi at fractional k, via shift-invert sparse diagonalisation — usable for num_wann far beyond what a dense H(k) can hold.

hr_dict()

Real-space H(R) as {R_tuple: scipy.sparse.csr_matrix} — the natural sparse in-memory form for large systems (feed to your own solver, KPM, pybinding, Kwant, …). O(nnz).

to_tbmodels(uc=None)

Convert to a tbmodels.Model (dense per-R; small/medium systems). uc defaults to the stored cell. From the returned model use .to_hr_file() (Wannier90 _hr.dat) or .to_hdf5_file(), or the tbmodels solvers directly.

to_pb()

Load directly into a pybinding pb.Lattice (one sublattice per orbital), built straight from the COO — no dense matrix. Mirrors tb_model.load(hdf5).to_pb() for the dense path.

to_pybinding()

Deprecated alias for to_pb() (kept for older sparse-API code).

to_kwant()

Build a bulk Kwant (Builder, lattice) (single site carrying all num_wann orbitals as matrix-valued blocks), straight from the COO — scales to large num_wann. Mirrors tb_model.load(hdf5).to_kwant().

to_hr_dat(path, uc=None, max_wann=4000)

Write a Wannier90 _hr.dat (via tbmodels). DENSE format: file size ~ num_R * num_wann**2, so it is guarded to small/medium systems. For large systems keep it sparse (hr_dict/Hk/KPM). Pass max_wann to override the guard.

to_hdf5(path, uc=None, max_wann=4000)

Write a tbmodels HDF5 (dense per-R). Same dense-format guard as to_hr_dat(). This is the format the tw_api_call auto-conversion produces for small systems, so downstream code that expects wannier90_hr.hdf5 keeps working unchanged.

tailwater.convert.as_tbmodels(src, uc=None)

Return src as a tbmodels.Model (the sparse -> dense bridge).

A dense input is returned unchanged; a sparse input is densified (uc defaults to its stored cell).

tailwater.convert.to_hdf5(src, path, **kw)

Write src to a tbmodels HDF5 (auto-detects sparse/dense).

DENSE on-disk format. For a sparse input the max_wann guard applies (pass max_wann= to override).

tailwater.convert.to_hr_dat(src, path, **kw)

Write src to a Wannier90 _hr.dat (auto-detects sparse/dense).

DENSE on-disk format (~ num_R * num_wann**2). For a sparse input the max_wann guard applies (pass max_wann= to override); prefer keeping large systems sparse.

tailwater.convert.to_pb(src, lattice_vectors=None, hop_threshold=1e-12)

Convert src to a pybinding pb.Lattice (auto-detects sparse/dense).

For a sparse input the lattice is built straight from the COO (no dense matrix); lattice_vectors is honored only for the dense path.

tailwater.convert.to_pythtb(src, lattice_vectors=None, hop_threshold=1e-12)

Convert src to a PythTB model (auto-detects sparse/dense).

PythTB is inherently dense, so a sparse input is densified via tbmodels first (small/medium systems only).

tailwater.convert.to_kwant(src, lattice_vectors=None, hop_threshold=1e-12)

Convert src to a Kwant kwant.Builder (auto-detects sparse/dense).

For a sparse input the builder carries matrix-valued blocks built straight from the COO (scales to large num_wann) and returns (Builder, lattice); lattice_vectors is honored only for the dense path.

Dense (tbmodels) converters

tailwater.client._to_pb_method(self, lattice_vectors=None, hop_threshold: float = 1e-12)

Convert this tbmodels.Model into a pybinding.Lattice.

Bound as model.to_pb on instances returned by tb_model.load(). After conversion, eigenvalues of model.hamilton(k_frac) and the pybinding model match to float32 precision (~1e-6 eV) at every k.

To plug pybinding into a band-structure calculation, pair the returned lattice with the companion helper k_cart_from_frac() — pybinding expects set_wave_vector(k_cart) in rad/length:

from tailwater import tb_model, k_cart_from_frac
import pybinding as pb

model = tb_model.load("wannier90_hr.hdf5")
lat   = model.to_pb()
pmod  = pb.Model(lat, pb.translational_symmetry())

for k_frac in k_path:
    pmod.set_wave_vector(k_cart_from_frac(k_frac, model.uc))
    bands.append(np.linalg.eigvalsh(pmod.hamiltonian.todense()))

Parameters:

• lattice_vectors (array-like (3, 3), optional) – Real-space lattice vectors as rows. If None, uses self.uc (the unit cell the tbmodels.Model carried when it was loaded); if that’s also None, falls back to np.eye(3).

• hop_threshold (float, default 1e-12) – Skip hops with |val| <= hop_threshold. Keep this low — we’re only filtering exact-zero entries from sparse hop storage; the band-relevant threshold should have been applied upstream when the HDF5 was first written.

Returns:

• pb.Lattice – with the same sublattices, on-site energies, lattice vectors, and hops as self, producing the same H(k) eigenvalues at every k.

• Conventions

• ———–

• **On-site doubling.* tbmodels’ Hamiltonian construction sums*

• stored[R] * exp(i k . R) over R, then adds its Hermitian

• conjugate to symmetrise. That second step supplies the missing

• minus-R half for nonzero R, but at R=0 it doubles the stored block

• on top of itself. tbmodels therefore stores half the user-supplied

• on-site value at R=0, and the round-trip Hamiltonian matches the

• physical Hamiltonian. Pybinding has no such double-up step, so

• we feed it twice the stored R=0 block, restoring the physical

• contribution.

• **Position basis.** tbmodels stores self.pos in fractional

• coordinates. Pybinding expects positions in Cartesian. We convert

• pos_cart = pos_frac @ LM so the resulting lattice’s

• Brillouin-zone and real-space geometry routines are physically

• meaningful. Eigenvalues are invariant under the per-orbital phase

• change induced by this choice — only the eigenvectors get rephased.

• **Hop duplicates.** For nonzero R, both (R, i, j) and

• (R, j, i) entries of the stored hop matrix are added explicitly

• to pybinding; the H.c. of each pybinding add automatically supplies

• the matching minus-R contribution, so the full Hamiltonian is

• reconstructed. For R = 0, the auto-implied H.c. of a given

• add_one_hopping call lands at the transposed indices —

• pybinding rejects the explicit second add as a duplicate. We catch

• that rejection silently.

tailwater.client._to_pythtb_method(self, lattice_vectors=None, hop_threshold: float = 1e-12)

Convert this tbmodels.Model into a pythtb.tb_model.

Bound as model.to_pythtb on instances returned by tb_model.load(). After conversion, eigenvalues of model.hamilton(k_frac) and py_model.solve_one(k_frac) match to float64 precision (~1e-12 eV) at every k.

The PythTB conversion is simpler than the pybinding one in two ways:

• PythTB takes orb (orbital positions) in fractional coordinates, the same convention as tbmodels.Model.pos, so no Cartesian conversion is needed.

• PythTB’s solve_one(k) accepts fractional k directly, so no analogue of k_cart_from_frac() is needed:

from tailwater import tb_model

model    = tb_model.load("wannier90_hr.hdf5")
py_model = model.to_pythtb()

# Sample H(k) at any fractional k:
eig = py_model.solve_one([0.0, 0.0, 0.0])    # Γ
eig = py_model.solve_one([0.5, 0.0, 0.0])    # M

# ...or use PythTB's built-in band-path helpers:
k_path, k_dist, k_node = py_model.k_path(
    [[0,0,0], [0.5,0,0], [0.333,0.333,0], [0,0,0]],
    nk=101, report=False,
)
bands = py_model.solve_all(k_path)            # (num_wann, nk)

Parameters:

• lattice_vectors (array-like (3, 3), optional) – Real-space lattice vectors as rows. If None, uses self.uc; if that’s also None, falls back to np.eye(3).

• hop_threshold (float, default 1e-12) – Skip hops with |val| <= hop_threshold. Keep this low — only filters exact-zero entries from sparse hop storage.

Returns:

• pythtb.tb_model – A 3D periodic tight-binding model with the same Hamiltonian as self. For slabs / wires, call PythTB’s .cut_piece(num, fin_dir) on the returned model.

• Conventions

• ———–

• The same on-site doubling story as _to_pb_method() applies

• here — PythTB counts each (0,0,0) entry once (like pybinding

• does), while tbmodels effectively counts it twice via its

• H += H.c. symmetrisation. To match, we multiply

• hop[(0,0,0)] by 2 before feeding PythTB.

• PythTB defaults to ``allow_conjugate_pair=False`` (only one of)

• each H.c. pair may be added explicitly. At R=(0,0,0) we therefore

• only add the upper-triangle off-diagonals; PythTB fills in the

• lower triangle via its automatic conjugate. For R != 0 we add

• every nonzero entry of hop[+R]; PythTB implies the

• corresponding (-R, j, i) contribution from each.

tailwater.client._to_kwant_method(self, lattice_vectors=None, hop_threshold: float = 1e-12)

Convert this tbmodels.Model into a kwant.Builder.

Bound as model.to_kwant on instances returned by tb_model.load(). Returns an unfinalised kwant.Builder with a 3D kwant.TranslationalSymmetry so callers can either:

• Finalise immediately for bulk H(k) sampling:

  import numpy as np, kwant
  from tailwater import tb_model
  
  model = tb_model.load("wannier90_hr.hdf5")
  syst  = kwant.wraparound.wraparound(model.to_kwant()).finalized()
  
  # Kwant's wraparound takes 2π·k_frac (per-cell Bloch phase),
  # NOT Cartesian rad/length.
  k_frac = [0.5, 0.0, 0.0]
  phase  = 2 * np.pi * np.asarray(k_frac)
  H      = syst.hamiltonian_submatrix(
      params=dict(k_x=phase[0], k_y=phase[1], k_z=phase[2]),
  )
  eigs   = np.sort(np.linalg.eigvalsh(H))
  

• Attach leads / build a finite scattering region on top of the bulk Builder for transport calculations:

  bulk = model.to_kwant()
  # Cut a finite slab, add leads, attach to bulk, etc.
  # See the Kwant tutorial: https://kwant-project.org/doc/
  

The returned eigenvalues match np.linalg.eigvalsh(model.hamilton(k_frac)) to ~float64 precision (~1e-12 eV) at every k.

Parameters:

• lattice_vectors (array-like (3, 3), optional) – Real-space lattice vectors as rows. If None, uses self.uc; if that’s also None, falls back to np.eye(3).

• hop_threshold (float, default 1e-12) – Skip hops with |val| <= hop_threshold. Keep this low — only filters exact-zero entries from sparse hop storage.

Returns:

• kwant.Builder – A 3D-periodic Builder with one site per Wannier orbital. Sublattices are accessible in the same order as model.pos via the Builder’s lattice.sublattices (the lattice object is the first argument of the TranslationalSymmetry stored on the Builder).

• Conventions

• ———–

• The same on-site doubling story as _to_pb_method() and

• _to_pythtb_method() applies — Kwant counts each (0,0,0)

• entry once, while tbmodels effectively counts it twice via its

• H += H.c. symmetrisation. To match, we multiply

• hop[(0,0,0)] by 2 before feeding Kwant.

• Positions are Cartesian (Kwant’s convention; converted from

• fractional via pos_cart = pos_frac @ LM).

• .. warning:: – The k-parameters that kwant.wraparound exposes (k_x, k_y, k_z) are not Cartesian rad/length like pybinding’s. They are the per-cell Bloch phase (i.e. 2π · k_frac), independent of the physical cell size. To sample H(k) at the same fractional k tbmodels uses, pass k_x = 2π · k_frac[0], etc. — see the worked example above. This is the most common source of “the Kwant bands don’t match the tbmodels bands” reports.

tailwater.client.k_cart_from_frac(k_frac, lattice_vectors)

Convert a fractional k-point to Cartesian (rad/length) for pybinding.

Pybinding’s set_wave_vector(k) expects k in rad/length — i.e. in the basis of the Cartesian reciprocal-lattice vectors b_i, not the fractional k_i Wannier90 and tbmodels use by default. The conversion is:

k_cart = 2π · inv(LM) @ k_frac

where LM has the real-space lattice vectors as rows.

Parameters:

• k_frac (array-like, shape (3,) or (N, 3)) – Fractional k (or batch of k-points), in the same units tbmodels.Model.hamilton(k) expects.

• lattice_vectors (array-like, shape (3, 3)) – Real-space lattice vectors as rows (e.g. model.uc).

Returns:

Cartesian k in rad/length, ready for pb.Model.set_wave_vector.

Return type:

np.ndarray of shape (3,) or (N, 3)

Example

import numpy as np, pybinding as pb
from tailwater import tb_model, k_cart_from_frac

model = tb_model.load("wannier90_hr.hdf5")
lat   = model.to_pb()
pmod  = pb.Model(lat, pb.translational_symmetry())

# Sample H(k) at Gamma → M (Bi2Se3) on a fractional path:
k_path_frac = np.array([[0,0,0], [0.5, 0, 0]])
bands = []
for kf in k_path_frac:
    pmod.set_wave_vector(k_cart_from_frac(kf, model.uc))
    bands.append(np.sort(np.linalg.eigvalsh(pmod.hamiltonian.todense())))

tailwater.hr_export.write_hr_output(hr_model: Model, out_path: str, fmt: str = 'hdf5')

Persist a tbmodels.Model. fmt is “hdf5” or “hr_dat”.

Parameters:

• hr_model (Model)

• out_path (str)

• fmt (str)

Return type:

str

tailwater.hr_export.build_hr_model(edge_pred, onsite_pred, gdata, LM, atoms: List[Tuple[str, List[float]]], hop_threshold: float = 0.01)

Build a tbmodels.Model from the model’s dense predictions.

Parameters:

• edge_pred (torch.Tensor or ndarray, shape [num_edges, 18, 18, 2]) – (or anything reshape-compatible). Real/imag in the last dim. Self-loop entries are overwritten by onsite_pred internally — caller may pass the raw head output unchanged.

• onsite_pred (torch.Tensor or ndarray, shape [num_atoms, 18, 18, 2].)

• gdata (PyG Data with edge_index, edge_vectors, inv_data,) – node_features. Same object the model consumed.

• LM (3x3 lattice matrix (rows = lattice vectors, Å). Passed) – to tbmodels.Model as uc.

• atoms ([(symbol, [x, y, z]), ...]. Per-atom Cartesian) – positions; used as the sublattice positions of the tbmodels.Model so the per-atom orbitals carry the same geometric labels they have in the structure.

• hop_threshold (drop hops with |val| <= this (eV). Default 0.01.)

Returns:

hr_model

Return type:

tbmodels.Model populated with on-site energies and hops.

tailwater.hr_export.build_hr_model_fast(edge_pred, onsite_pred, gdata, LM, atoms: List[Tuple[str, List[float]]], hop_threshold: float = 0.01)

Vectorized equivalent of build_hr_model — same output, much faster.

The reference build_hr_model loops l × atm1 × s1o × atm2 × s2o in pure Python, which is ~num_R * N^2 * 18^2 inner iterations and crosses into multi-minute territory for ~50-atom inputs. This function keeps the EXACT same iteration order (so tbmodels’ first-add-wins semantics and the model’s non-Hermitian inter-atom predictions produce the same final recorded values) but lifts the threshold check, the active-mask filter, the magnitude computation, the value gather, and the orbital-index lookup out of Python and into NumPy. The surviving per-hop Python loop only visits hops that actually got recorded.

Order preservation is what makes this drop-in safe: np.nonzero on a [num_R, N, 18, N, 18] mask returns indices in C (row-major) order, which is identical to the reference’s outer-to-inner nested-loop ordering of those five axes.

Approximate speedup on a 50-atom material: ~100-300x for the hop insertion phase. Output (the tbmodels.Model) is byte-identical to build_hr_model’s output as long as pb.Lattice / tbmodels.Model.add_hop are deterministic — both are.

Parameters:

• atoms (List[Tuple[str, List[float]]])

• hop_threshold (float)

Return type:

Model