hamop

hamop is a Python package for computing what electrons do in a material or molecule described by a tight-binding model: a short list of orbitals and the energies that let an electron hop between them. You build the model once. From that one model the package computes:

• the allowed electron energies (the band structure) and how many states sit at each energy (the density of states);
• how strongly the material absorbs light at each photon energy (the optical conductivity);
• how easily electrons pass through a small device placed between two wires (the transmission), and from it the conductance, current and thermovoltage a measurement at a given temperature would report;
• topological numbers such as the Chern number, and related geometric quantities of the electron states;
• for your own measurements: model parameters fitted to measured band energies with error bars, error bars for anything computed from them, and the transmittance of a 2D sample that an optical measurement would see.

The point of the package is that all of these come from the same matrices. When the optics and the energy levels of one model are computed by different codes, small differences in conventions (how the velocity is defined, how overlapping orbitals are treated, how peaks are broadened) can make the results disagree, for example an absorption edge that no longer sits at the band gap. Here the spectrum and the optics are computed from the same Bloch matrices through the same eigensolver, the topology functions work from the same hopping blocks, and principal_layers cuts the same model into the layer blocks that the transport functions take.

The package also works when neighbouring orbitals overlap (a nonorthogonal basis, as in LCAO calculations, which build electron states from atomic orbitals), and it can read Hamiltonians exported by Wannier90, a widely used program that turns first-principles (quantum-mechanical, parameter-free) calculations into tight-binding models. It ships no material data: the numbers of a real material come from your own calculation or measurement.

Contents

• A short guide to the words used here
• Install, units and conventions
• Examples (each with the output it prints)
• What is in the package
• When it refuses, and why
• How the results are checked
• What changed in 0.12.0
• What changed in 0.11.0
• Bugs fixed in 0.10.1
• Limits
• Relation to existing tools
• Where it comes from
• Citing, support and license

A short guide to the words used here

• Site, orbital -- an atom (site) carries one or more orbitals, the states an electron can occupy there.
• Hopping -- the energy that couples an orbital on one site to an orbital on another. You give each coupling once, in one direction; the reverse direction is added for you.
• Hermitian -- a matrix equal to its own conjugate transpose. An energy matrix must be Hermitian so that its energies are real; the reverse hopping added for you is the "Hermitian partner" of the one you give.
• Overlap, nonorthogonal basis -- in some models neighbouring orbitals are not independent: they overlap by an amount s. The overlap matrix is called S. If you give no overlap, S is the identity (an orthogonal basis).
• Periodic model, unit cell -- a crystal repeats one unit cell. The package handles periodicity in any number of directions, and also finite models (a molecule, a flake, a device) with no repetition (cell=None).
• Wave vector k, Brillouin zone, k-mesh -- in a crystal each electron state carries a wave vector k. All distinct k fill one cell of "k-space" called the Brillouin zone. Many results are sums over an even grid of k points, the k-mesh (mesh= in the functions).
• Bloch Hamiltonian H(k) -- the model's energy matrix at one k. Its eigenvalues are the band energies at that k.
• Band structure, band gap -- the band energies as k varies. A gap is an energy range with no states.
• Density of states (DOS) -- how many states lie near each energy.
• Chemical potential mu, filling -- states below mu are occupied (smoothly, at temperature T); the filling is the number of occupied states per cell.
• Optical conductivity -- how strongly light of photon energy hbar omega is absorbed. It is computed with the Kubo-Greenwood formula (the standard weak-field response formula) and returned in units of e^2/(4 hbar). That unit is the measured absorption plateau of graphene, so graphene gives about 1.
• Broadening eta -- a small energy width given to each sharp line, so that sums over a finite k-mesh give smooth curves.
• Drude weight -- the part of the conductivity carried by freely moving electrons in a metal (the "intraband" part).
• Transmission T(E) -- the number of electron channels that pass through a device at energy E; in the Landauer picture the conductance is this number times e^2/h per spin. The device sits between two semi-infinite leads (wires), and is cut into principal layers that couple only to their nearest neighbours. The method is the nonequilibrium Green function (NEGF) method.
• Conductance, Seebeck coefficient, Lorenz ratio -- at a finite temperature T a measurement averages T(E) over a window of a few kT around the chemical potential (k is Boltzmann's constant). The conductance is that average times the number of spin directions, in units of e^2/h. The Seebeck coefficient is the voltage that appears per kelvin of temperature difference across the device. The Lorenz ratio is the heat conductance carried by the electrons divided by G T; for a smooth T(E) it is close to (pi^2/3)(k/e)^2 = 2.443e-8 V^2/K^2.
• Green function -- the matrix (zS - H)^-1 at the complex energy z = E + i eta. Its imaginary part gives the local density of states (the DOS on each orbital, LDOS), and the transmission is built from it.
• Self-energy -- a term added to a layer's Hamiltonian that describes what the rest of the system (a lead, disorder, a probe) does to it.
• Zeeman term, spin-orbit coupling -- energy terms that act on the electron's spin: the Zeeman term comes from a magnetic field and splits each level by 2B (B in eV); spin-orbit coupling comes from the electron's own motion.
• Berry phase, Chern number -- properties of how the electron states change as k moves. The Chern number is a whole number that cannot change unless a band gap closes. It sets the quantized Hall conductivity C e^2/h.
• Quantum metric -- how fast the occupied electron states change as k moves; the Berry curvature is its partner.
• Kernel polynomial method (KPM) -- a way to get the DOS and the conductivity of very large finite models without computing every eigenvalue.

Install, units and conventions

pip install hamop

It needs Python 3.9 or newer, NumPy 1.22 or newer and SciPy 1.8 or newer, and nothing else.

Units and conventions, stated once:

• Energies in eV, positions in Angstrom, k in 1/Angstrom, all Cartesian. Temperatures T are in kelvin, except berry_dipole, which takes kT in energy units.
• Each directed hopping block is added once; its Hermitian partner is implied.
• The optical conductivity is the real sheet conductivity in units of e^2/(4 hbar). Spin degeneracy enters as an explicit factor spin (default 2). For a model made spinful with with_spin, pass spin=1. For a finite model the result is the conductivity times the system area; divide by your geometric area.
• dos counts states per unit cell per eV without spin; fermi_level takes the filling without spin; carrier_count includes spin.
• The velocity operator uses the usual tight-binding position convention: the position operator is diagonal at the sites. Extra on-site position blocks (model.set_dipole) add the intra-atomic term on top.
• Magnetic flux for with_peierls is given per square Angstrom, in units of the flux quantum; magnetic_supercell takes the flux p/q per unit cell.
• Temperature T = 0 means the exact step occupation; a level exactly at mu counts as half occupied (the value the Fermi function has at mu for every T > 0).
• Transport at finite temperature: conductance in units of e^2/h (multiply by e^2/h = 3.874e-5 S for siemens), current in units of e^2/h times 1 V (multiply by the same number for amperes), Seebeck coefficient in V/K, Lorenz ratio in V^2/K^2. spin=2 counts both spin directions, as in the optics.

Examples

Each example below runs as written, and the output shown is what it printed with hamop 0.12.0. The model parameters are illustrative values, not fitted to any material, unless the text says otherwise.

1. Graphene: bands and the absorption plateau

import numpy as np
from hamop import graphene, bands, sigma_optical

g = graphene(t=-2.7, a=2.46)      # hopping t in eV, lattice constant a in Angstrom

b = 2 * np.pi * np.linalg.inv(g.cell).T   # reciprocal lattice vectors (rows)
K = (2 * b[0] + b[1]) / 3                 # a corner of the Brillouin zone
e_gamma = bands(g, [np.zeros(2)])[0]
e_K = bands(g, [K])[0]
print(f"band energies at the zone centre: {e_gamma[0]:+.4f} {e_gamma[1]:+.4f} eV")
print("gap at the zone corner K below 1e-9 eV:", abs(e_K[1] - e_K[0]) < 1e-9)

omega = np.array([1.0, 1.3])              # photon energies, eV
sigma = sigma_optical(g, omega, mu=0.0, mesh=120, eta=0.12, T=10.0)
for w, s in zip(omega, sigma):
    print(f"photon {w:.1f} eV: sigma = {s:.3f} x e^2/(4 hbar)")

band energies at the zone centre: -8.1000 +8.1000 eV
gap at the zone corner K below 1e-9 eV: True
photon 1.0 eV: sigma = 1.020 x e^2/(4 hbar)
photon 1.3 eV: sigma = 1.032 x e^2/(4 hbar)

The two bands sit at -/+ 3|t| at the zone centre and touch at the corner K (the "Dirac point"). The absorption is close to 1 in units of e^2/(4 hbar), the universal optical conductivity of graphene (Kuzmenko et al., Phys. Rev. Lett. 100, 117401 (2008)). The test suite checks these two photon energies to within 5 %. The script examples/graphene_universal_conductivity.py prints the same quantity from 0.6 to 1.8 eV.

2. Build your own model: a dimerized chain

import numpy as np
from hamop import TightBindingModel, band_edges, berry_phase

def ssh_chain(t1, t2, a=2.0):
    """Two sites per cell, at 0 and a/2; hopping t1 inside the cell,
    t2 to the next cell."""
    m = TightBindingModel(positions=[[0.0], [0.5 * a]], norb=1, cell=[[a]])
    m.add_hop(0, 1, (0,), [[t1]])        # site 0 -> site 1, same cell
    m.add_hop(1, 0, (1,), [[t2]])        # site 1 -> site 0, next cell
    return m

for t1, t2 in [(-1.0, -0.6), (-0.6, -1.0)]:
    m = ssh_chain(t1, t2)
    vbm, cbm, gap = band_edges(m, mu=0.0, mesh=2001)
    b = 2 * np.pi / 2.0                              # reciprocal period
    loop = [[b * i / 60] for i in range(60)]         # one pass through k
    zak = berry_phase(m, loop, n_occ=1)
    print(f"t1={t1:+.1f} t2={t2:+.1f}: gap {gap:.4f} eV, "
          f"Zak phase {zak / np.pi:+.4f} pi")

t1=-1.0 t2=-0.6: gap 0.8000 eV, Zak phase +0.0000 pi
t1=-0.6 t2=-1.0: gap 0.8000 eV, Zak phase -1.0000 pi

This is the Su-Schrieffer-Heeger (SSH) chain. Both orderings of the bonds have the same gap, 2 | |t1| - |t2| |. They differ in the Berry phase of the filled band (here called the Zak phase): the two values differ by pi. Which of the two is 0 depends on where the unit cell starts; the difference of pi does not.

3. Topology: Chern number and the quantized Hall conductivity

import numpy as np
from hamop import haldane, chern_number, sigma_tensor

for mass in (0.0, 0.9):
    h = haldane(t1=-1.0, t2=0.1, phi=np.pi / 2, m_ab=mass)
    C = chern_number(h, mesh=18)
    sxy = sigma_tensor(h, [0.0], mu=0.0, directions=(0, 1), mesh=48,
                       T=10.0, eta=1e-4, spin=1)[0]
    hall = sxy.real * np.pi / 2       # e^2/h is 2/pi in the unit e^2/(4 hbar)
    print(f"mass {mass}: Chern number {round(C, 6) + 0.0:+.6f}, "
          f"sigma_xy(0) = {round(hall, 6) + 0.0:+.6f} e^2/h")

mass 0.0: Chern number +1.000000, sigma_xy(0) = +1.000000 e^2/h
mass 0.9: Chern number +0.000000, sigma_xy(0) = +0.000000 e^2/h

The Haldane model (Haldane, Phys. Rev. Lett. 61, 2015 (1988)) has a topological phase (Chern number 1) and, when the sublattice energy difference m_ab is large enough, an ordinary one (Chern number 0). Two separate calculations agree: the Chern number from the band states, and the Hall conductivity at zero frequency from the Kubo formula. Their agreement is the TKNN relation (Thouless, Kohmoto, Nightingale and den Nijs, Phys. Rev. Lett. 49, 405 (1982)).

4. Transport: one impurity in a chain

import numpy as np
from hamop import chain_lead_blocks, transmission

t, eps = -1.0, 0.8                     # hopping and impurity energy, eV
H00, H01 = chain_lead_blocks(t=t)      # one lead layer and its coupling
layers = [H00, H00 + eps, H00]         # three device layers, impurity in the middle
coup = [H01, H01]                      # couplings between device layers

E = np.array([-1.2, 0.3, 1.5, 2.5])
T = transmission(E, layers, coup, H00, H01, eta=1e-8)
exact = np.where(np.abs(E) < 2, (4 - E**2) / ((4 - E**2) + eps**2), 0.0)
for e, a, b in zip(E, T, exact):
    print(f"E = {e:+.1f} eV: T = {a:.6f}   closed form {b:.6f}")

E = -1.2 eV: T = 0.800000   closed form 0.800000
E = +0.3 eV: T = 0.859341   closed form 0.859341
E = +1.5 eV: T = 0.732218   closed form 0.732218
E = +2.5 eV: T = 0.000000   closed form 0.000000

A clean chain passes one channel (T = 1) at energies inside its band, between -2|t| and +2|t|, and none outside. A single impurity of energy eps reflects part of the wave; the closed form for this chain is T = (4t^2 - E^2) / ((4t^2 - E^2) + eps^2). At 2.5 eV the energy is outside the band, so nothing passes.

5. Overlapping orbitals: moving the energy zero changes nothing

from hamop import linear_chain, drude_weight

# A chain whose neighbouring orbitals overlap (s = 0.2), half filled.
m1 = linear_chain(t=-1.0, e0=0.0, s=0.2)
D1 = drude_weight(m1, mu=0.0, mesh=800, T=100.0)

# The same chain with the energy zero moved by c = 5 eV.  With overlap
# that shift is H -> H + c S, so the hopping changes by c * s as well.
c = 5.0
m2 = linear_chain(t=-1.0 + c * 0.2, e0=c, s=0.2)
D2 = drude_weight(m2, mu=c, mesh=800, T=100.0)
print(f"Drude weight, original zero: {D1:.8f}")
print(f"Drude weight, shifted zero:  {D2:.8f}")
print("difference below 1e-10:", abs(D1 - D2) < 1e-10)

Drude weight, original zero: 16.01338983
Drude weight, shifted zero:  16.01338983
difference below 1e-10: True

Where you put the zero of energy is a free choice, so no physical result may depend on it. With overlapping orbitals, that only holds if the velocity includes the overlap term v = dH/dk - (E_n + E_m)/2 dS/dk, which the package uses for the Drude weight, the optical conductivity and the conductivity tensor.

6. Fit a model to measured bands, and plan the measurement first

import numpy as np
from hamop import linear_chain, band_information, fit_bands

def build(theta):                      # theta = (t, e0), both in eV
    return linear_chain(t=theta[0], e0=theta[1])

k = np.linspace(0.1, 3.0, 12)[:, None] # 12 measured k-points (1/Angstrom)
sigma = 0.03                           # measurement error, eV

plan = band_information(build, [-1.0, 0.2], k, sigmas=sigma)
print("identifiable:", plan["identifiable"])
print("predicted error bars (eV):", np.round(plan["sigma"], 5))

# Illustrative "measured" data: the chain with t = -1.0, e0 = 0.2,
# plus seeded random noise of size sigma.
rng = np.random.default_rng(13)
measured = 0.2 + 2 * (-1.0) * np.cos(k[:, 0]) + sigma * rng.standard_normal(12)
fit = fit_bands(build, [-1.3, 0.5], k, measured, sigmas=sigma)
print("fitted (t, e0):", np.round(fit.theta, 4))
print("error bars:", np.round(fit.sigma, 5))
print(f"chi-squared {fit.chi2:.2f} for {fit.chi2_dof} degrees of freedom")

identifiable: True
predicted error bars (eV): [0.0061  0.00866]
fitted (t, e0): [-1.0002  0.2087]
error bars: [0.0061  0.00866]
chi-squared 18.87 for 10 degrees of freedom

build is your own parameterization: any function that turns a list of numbers into a model. Before measuring, band_information says whether the planned k-points can tell the parameters apart and how large the error bars will be; design_kpoints picks the most useful k-points from a candidate list. fit_bands then fits the measured energies and reports error bars from the standard least-squares covariance.

7. Quantum geometry of the occupied band

import numpy as np
from hamop import haldane, quantum_geometric_tensor, quantum_weight

m = haldane(t1=-1.0, t2=0.1)                 # lattice constant 1 Angstrom
b = 2 * np.pi * np.linalg.inv(m.cell).T
K = (2 * b[0] + b[1]) / 3
Q = quantum_geometric_tensor(m, K)           # at the zone corner
g = Q.real                                   # quantum metric, Angstrom^2
omega = 2 * Q[0, 1].imag                     # Berry curvature, Angstrom^2
print(f"tr g = {np.trace(g):.4f}, 2 sqrt(det g) = "
      f"{2 * np.sqrt(np.linalg.det(g)):.4f}, |Omega| = {abs(omega):.4f}")

out = quantum_weight(m, mesh=18)             # integrated over the zone
print(f"tr K = {out['tr_K']:.4f}, Chern number from the same data = {out['chern']:.4f}")

tr g = 1.3889, 2 sqrt(det g) = 1.3889, |Omega| = 1.3889
tr K = 1.1827, Chern number from the same data = 1.0000

At every k, tr g >= 2 sqrt(det g) >= |Omega| must hold; at this point all three are equal to four decimals. Integrated over the zone the same chain gives tr K >= |C| (Onishi and Fu, PRX 14, 011052 (2024)): here 1.1827 against a Chern number of 1.

8. Reading a Wannier90 file

import os, tempfile
import numpy as np
import hamop

# Write graphene's real-space Hamiltonian in the Wannier90 hr.dat format,
# read it back, and compare with the built-in graphene model.
t, a = -2.7, 2.46
H_R = {(0, 0): [[0, t], [t, 0]],
       (1, 0): [[0, 0], [t, 0]], (-1, 0): [[0, t], [0, 0]],
       (0, 1): [[0, 0], [t, 0]], (0, -1): [[0, t], [0, 0]]}
path = os.path.join(tempfile.mkdtemp(), "graphene_hr.dat")
hamop.save_wannier90_hr(path, H_R)

cell = a * np.array([[1.0, 0.0], [0.5, np.sqrt(3) / 2]])
centres = np.array([[0.0, 0.0], (cell[0] + cell[1]) / 3])
m = hamop.from_wannier90(path, cell, centers=centres)

k = np.random.default_rng(0).uniform(-2, 2, (20, 2))
diff = np.abs(hamop.bands(m, k) - hamop.bands(hamop.graphene(t, a), k)).max()
print("bands agree to 1e-12 eV at 20 random k-points:", diff < 1e-12)

bands agree to 1e-12 eV at 20 random k-points: True

A seedname_hr.dat file (Mostofi et al., Comput. Phys. Commun. 185, 2309 (2014)) holds the real-space Hamiltonian but not the lattice vectors or the orbital centres. You supply the lattice vectors (cell). The centres are optional: band energies do not depend on them, but optical matrix elements do, so give the true centres when you need optics.

9. Transport at a finite temperature: conductance, current, thermovoltage

import numpy as np
from hamop import (chain_lead_blocks, transmission, landauer_conductance,
                   landauer_current, thermoelectric)

t, eps = -1.0, 0.8                     # the impurity chain of example 4
H00, H01 = chain_lead_blocks(t=t)
layers, coup = [H00, H00 + eps, H00], [H01, H01]
E = np.arange(-1.0, 1.6, 0.002)        # energy grid, eV
T_E = transmission(E, layers, coup, H00, H01, eta=1e-8)

mu = 0.3
exact = 2 * (4 - mu**2) / ((4 - mu**2) + eps**2)   # 2 T(mu), spin 2
for temp in (30.0, 300.0):
    G = landauer_conductance(E, T_E, mu, T=temp)
    print(f"T = {temp:5.1f} K: G = {G:.6f} e^2/h   (2 T(mu) = {exact:.6f})")

r = thermoelectric(E, T_E, mu, T=300.0)
print(f"Seebeck coefficient: {r['S'] * 1e6:+.4f} microvolt/K")
print(f"Lorenz ratio: {r['lorenz']:.4e} V^2/K^2")

I = landauer_current(E, T_E, mu_L=mu + 0.05, mu_R=mu - 0.05, T=300.0)
print(f"current at 0.1 V bias: {I:.6f} x e^2/h x 1 V")

T =  30.0 K: G = 1.718680 e^2/h   (2 T(mu) = 1.718681)
T = 300.0 K: G = 1.718534 e^2/h   (2 T(mu) = 1.718681)
Seebeck coefficient: +0.1589 microvolt/K
Lorenz ratio: 2.4423e-08 V^2/K^2
current at 0.1 V bias: 0.171848 x e^2/h x 1 V

transmission gives T(E) at the energies you choose; the three functions used here turn it into what an experiment at temperature T reports. At 30 K the conductance agrees with 2 T(mu) of the closed form to about 1e-6; at 300 K the thermal average over the curved T(E) lowers it slightly. T(E) falls with energy at mu = 0.3 eV, so the Seebeck coefficient is positive (hole-like); the rough low-temperature estimate -(pi^2/3)(k^2 T/e) d ln T/dE (the Mott formula) gives +0.158 microvolt/K here. The energy grid must reach 30 kT beyond the chemical potentials and resolve kT; otherwise the functions refuse and say how to fix the grid.

10. Reading a Wannier90 tb.dat file (cell and centres included)

import os, tempfile
import numpy as np
import hamop

# A seedname_tb.dat in the Wannier90 layout: date line, lattice vectors
# (Angstrom), num_wann, nrpts, degeneracies, then H(R) and r(R) blocks.
# Here: one orbital centred at x = 0.4 Angstrom in a chain of period 2.0,
# on-site energy 0.1 eV and hopping -1.0 eV to both neighbours.
text = """ written by hand
  2.0 0.0 0.0
  0.0 10.0 0.0
  0.0 0.0 10.0
 1
 3
    1    1    1

   -1    0    0
    1    1   -0.10000000E+01  0.00000000E+00
    0    0    0
    1    1    0.10000000E+00  0.00000000E+00
    1    0    0
    1    1   -0.10000000E+01  0.00000000E+00

   -1    0    0
    1    1    0.00000000E+00  0.0  0.0  0.0  0.0  0.0
    0    0    0
    1    1    0.40000000E+00  0.0  0.0  0.0  0.0  0.0
    1    0    0
    1    1    0.00000000E+00  0.0  0.0  0.0  0.0  0.0
"""
path = os.path.join(tempfile.mkdtemp(), "chain_tb.dat")
with open(path, "w") as fh:
    fh.write(text)

m = hamop.from_wannier90_tb(path)       # cell and centre come from the file
print("periodic directions:", m.cell.shape[0], " cell:", m.cell.ravel(),
      " centre:", m.positions.ravel())
k = np.linspace(-1.5, 1.5, 7)[:, None]
ref = hamop.linear_chain(t=-1.0, e0=0.1, a=2.0)
diff = np.abs(hamop.bands(m, k) - hamop.bands(ref, k)).max()
print("bands agree with linear_chain to 1e-12 eV:", diff < 1e-12)

periodic directions: 1  cell: [2.]  centre: [0.4]
bands agree with linear_chain to 1e-12 eV: True

seedname_tb.dat (written by Wannier90 with write_tb = true) holds the lattice vectors, the Hamiltonian H(R) and the position matrix elements <m0|r|nR> in one file, so from_wannier90_tb needs no cell and no centres from you. The orbital positions are taken from the diagonal elements <m0|r|m0>; the off-diagonal ones are read by load_wannier90_tb but not used (see Limits). Directions in which no hopping occurs are dropped, so this chain imports as a 1D model.

11. Error bars for what you did not measure, and a check of the model

import numpy as np
from hamop import linear_chain, bands, fit_bands, propagate_uncertainty

def build(theta):                      # theta = (t, e0), both in eV
    return linear_chain(t=theta[0], e0=theta[1])

k = np.linspace(0.1, 3.0, 12)[:, None] # the measurement of example 6
sigma = 0.03
rng = np.random.default_rng(13)
measured = 0.2 + 2 * (-1.0) * np.cos(k[:, 0]) + sigma * rng.standard_normal(12)
fit = fit_bands(build, [-1.3, 0.5], k, measured, sigmas=sigma)
print(f"p-value: {fit.p_value:.3f}")

# Error bars of quantities that were not measured directly.
width = propagate_uncertainty(fit, lambda th: 4 * abs(th[0]))
top = propagate_uncertainty(fit, lambda th: bands(build(th), [[np.pi]])[0, 0])
print(f"band width 4|t|: {width['value']:.4f} +/- {width['sigma']:.4f} eV")
print(f"band top at k = pi: {top['value']:.4f} +/- {top['sigma']:.4f} eV")

# The same chain with a second-neighbor hopping of 0.1 eV that the
# fitted model does not have: the p-value flags the missing physics.
wrong = measured + 2 * 0.1 * np.cos(2 * k[:, 0])
bad = fit_bands(build, [-1.3, 0.5], k, wrong, sigmas=sigma)
print(f"chi-squared {bad.chi2:.1f} for {bad.chi2_dof} degrees of freedom, "
      f"p-value {bad.p_value:.1e}")

p-value: 0.042
band width 4|t|: 4.0010 +/- 0.0244 eV
band top at k = pi: 2.2092 +/- 0.0151 eV
chi-squared 246.3 for 10 degrees of freedom, p-value 3.2e-47

propagate_uncertainty turns the fit's covariance into the error bar of any number you compute from the parameters: here the band width and the band energy at k = pi, a point the measurement did not reach. By default it propagates to first order, and it first checks that the quantity is close to linear over one standard deviation of the parameters. If it is not (for example at a band crossing, or where the quantity has a minimum), the call is refused and method="monte_carlo" samples the parameters instead.

fit.p_value is the chance of a chi-squared at least this large if the model is right and the stated errors are correct. The first fit gives 0.042: these data were generated from the model, and a value this low happens by chance in about 4 % of such experiments. The second data set contains physics the model lacks, and its p-value of 3e-47 says so clearly. A small p-value means either the errors are underestimated or the model is incomplete; the fit alone cannot tell which.

12. A 2D sample in an optical measurement

import numpy as np
from hamop import (graphene, sigma_optical, sheet_optics,
                   sheet_conductivity_from_transmittance)

g = graphene(t=-2.7, a=2.46)
omega = np.array([1.0, 1.3])                # photon energies, eV
sigma = sigma_optical(g, omega, mu=0.0, mesh=120, eta=0.12, T=10.0)

free = sheet_optics(sigma)                  # free-standing sheet in air
on_glass = sheet_optics(sigma, n_in=1.0, n_out=1.5)   # on a substrate, index 1.5
for w, T, A, T_rel in zip(omega, free["T"], free["A"], on_glass["T_rel"]):
    print(f"photon {w:.1f} eV: free-standing T = {T:.4f}, A = {A:.4f}; "
          f"on the substrate T/T_bare = {T_rel:.4f}")

# The other way round: a measured T/T_bare of 0.9805 +/- 0.0010 on the
# same substrate gives the sheet conductivity in the package's unit.
m = sheet_conductivity_from_transmittance(0.9805, n_in=1.0, n_out=1.5,
                                          T_rel_err=0.0010)
print(f"measured Re sigma = {m['sigma']:.3f} +/- {m['sigma_err']:.3f} "
      "x e^2/(4 hbar)")

photon 1.0 eV: free-standing T = 0.9770, A = 0.0229; on the substrate T/T_bare = 0.9815
photon 1.3 eV: free-standing T = 0.9768, A = 0.0231; on the substrate T/T_bare = 0.9813
measured Re sigma = 1.079 +/- 0.056 x e^2/(4 hbar)

An optical experiment on a 2D material reports the transmittance T (the fraction of light that passes), often as the ratio to the bare substrate, T/T_bare. sheet_optics computes it, together with the reflectance R and the absorptance A, from the sheet conductivity of a 2D model at normal incidence. The free-standing graphene model absorbs about 2.3 %, the value pi alpha measured by Nair et al., Science 320, 1308 (2008) (alpha is the fine-structure constant). sheet_conductivity_from_transmittance goes the other way, from a measured T/T_bare and its error to Re sigma with an error bar, so a measured spectrum can be compared with sigma_optical point by point. The substrate index 1.5 and the measured value 0.9805 are illustrative numbers, not data for a real sample.

What is in the package

Each function's docstring (help(hamop.sigma_optical), for example) gives its inputs, units and conventions.

The model

• TightBindingModel(positions, norb, cell=None) -- sites, orbitals per site and lattice vectors. add_hop(i, j, image, H_block, S_block=None) adds a hopping (and optional overlap) block; set_dipole(i, X) adds on-site position blocks; bloch(k) and bloch_derivative(k, direction) return H(k), S(k) and their exact k-derivatives; monkhorst_pack(mesh, time_reversal=False) returns an even k-grid with weights.
• gen_eigh(H, S) -- eigenvalues of H c = E S c. It drops overlap directions below a threshold ("canonical orthogonalization"; Szabo and Ostlund, Modern Quantum Chemistry, sec. 3.4.5), so a nearly redundant set of orbitals does not produce huge wrong eigenvalues.
• Ready-made models: linear_chain, two_site (a two-atom molecule), graphene, ssh, haldane, and chain_lead_blocks (lead blocks of a chain, for transport).

Bands and density of states

• bands, k_path (a path through k-space for a band plot), dos (Gaussian-broadened DOS), fermi_level (chemical potential for a given filling), band_edges (valence-band top, conduction-band bottom and gap around mu).

Optics

• sigma_optical -- the real optical conductivity (Gaussian or Lorentzian broadening).
• sigma_tensor -- the complex conductivity tensor component sigma_ab(omega) between bands, including the Hall component sigma_xy.
• drude_weight -- the intraband (Drude) weight.
• carrier_count -- occupied states per cell, spin included.

Topology and band geometry

• berry_phase (along a closed loop in k-space), berry_curvature and chern_number (on a 2D k-grid, by the lattice method of Fukui, Hatsugai and Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005)). With overlap, two conventions are offered (frame="lowdin" or "atomic"); a sparse solver (solver="sparse") serves large cells with few occupied bands.
• chern_marker -- the local Chern marker of a finite 2D model (Bianco and Resta, Phys. Rev. B 84, 241106(R) (2011)).
• berry_dipole, band_curvatures, fermi_occupation -- the Berry curvature dipole behind the nonlinear Hall effect (Sodemann and Fu, Phys. Rev. Lett. 115, 216806 (2015)), computed in two independent ways (method="grad" or "fermi"), the per-band curvature, and the Fermi occupation factor.
• quantum_geometric_tensor, quantum_metric, quantum_weight -- the quantum metric and Berry curvature of the occupied states (Provost and Vallee, Commun. Math. Phys. 76, 289 (1980); Peotta and Torma, Nat. Commun. 6, 8944 (2015); Onishi and Fu, PRX 14, 011052 (2024)).

Spin and magnetic fields

• with_spin -- a spinful copy of a model; PAULI -- the Pauli matrices, for adding Zeeman or spin-orbit terms; kane_mele -- the Kane-Mele model (Kane and Mele, Phys. Rev. Lett. 95, 226801 (2005)).
• with_peierls -- a uniform magnetic field on a finite model, by Peierls substitution: each hopping gets a phase set by the magnetic flux (Peierls, Z. Phys. 80, 763 (1933)).
• magnetic_supercell -- a periodic 2D model at rational flux p/q per cell: the unit cell is repeated q times so that the field fits it (a Hofstadter magnetic cell).

Transport

• transmission -- T(E) by a recursive sweep through the layers; transmission_direct and transmission_sparse compute the same quantity by dense and by sparse inversion. sancho_rubio gives the surface Green function of a lead (Lopez Sancho, Lopez Sancho and Rubio, J. Phys. F 15, 851 (1985)); every result is checked against the lead's own Dyson equation, and when the decimation fails that check it is recomputed from the lead's decaying modes (Lee and Joannopoulos, Phys. Rev. B 23, 4988 and 4997 (1981)).
• principal_layers -- cuts a finite model into layers along an axis and checks that no coupling skips a layer.
• buttiker_transmission -- one dephasing probe, an imaginary contact that scrambles the phase of the electrons passing it but draws no net current (Buttiker, Phys. Rev. B 33, 3020 (1986)); multiprobe_transmission -- probes on many layers (D'Amato and Pastawski, Phys. Rev. B 41, 7411 (1990)); scba_transmission -- averaged on-site disorder, solved self-consistently (the self-consistent Born approximation, SCBA).
• device_greens, device_ldos, bond_currents -- the device Green function, the local density of states per orbital, and the map of currents between orbitals (Paulsson and Brandbyge, Phys. Rev. B 76, 115117 (2007)).
• transmission, transmission_direct and transmission_sparse also accept your own self-energy per layer (sigma_int).
• landauer_conductance, landauer_current, thermoelectric -- the conductance, the current at a finite bias, and the Seebeck coefficient, electronic heat conductance and Lorenz ratio at a finite temperature, from T(E) sampled on an energy grid (example 9; formulas from Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press (1995), ch. 2, and Sivan and Imry, Phys. Rev. B 33, 551 (1986)).

Large systems

• bloch_sparse, bloch_derivative_sparse -- the same matrices as bloch and bloch_derivative, stored sparse.
• lowest_bands -- only the lowest few eigenvalues (Lanczos method, an iterative solver that never forms the full eigenvector set).
• kpm_dos, kpm_sigma -- DOS and optical conductivity of large finite models by the kernel polynomial method (Weisse, Wellein, Alvermann and Fehske, Rev. Mod. Phys. 78, 275 (2006)).

Fewer k-points

• monkhorst_pack(..., time_reversal=True) pairs k with -k.
• symmetry_fold folds a k-grid by a point group you supply; find_point_group finds that group. Both check every operation on the model first.

Interpolation and file exchange

• fourier_interpolation, FourierInterpolator -- sample H(k) on a grid, transform to real space, and evaluate at any k.
• from_wannier90, load_wannier90_hr, save_wannier90_hr -- read and write the Wannier90 seedname_hr.dat format.
• from_wannier90_tb, load_wannier90_tb -- read the Wannier90 seedname_tb.dat format, which also carries the lattice vectors and the position matrix elements (example 10).

Fitting to measurements

• fit_bands, BandFit, band_information, design_kpoints -- example 6; design_kpoints uses greedy D-optimal selection, which adds, one at a time, the k-point that most shrinks the joint uncertainty of the parameters (Pukelsheim, Optimal Design of Experiments, SIAM (2006)).
• propagate_uncertainty and BandFit.p_value (example 11): the error bar of any quantity computed from fitted parameters, by first-order propagation through the fit covariance (refused when the quantity is not close to linear over one standard deviation) or by Monte Carlo sampling; and the chi-squared p-value of a fit, which says whether the data are consistent with the model and the stated errors.

Comparing with optical measurements

• sheet_optics, sheet_conductivity_from_transmittance (example 12): the transmittance, reflectance and absorptance of a 2D sheet at normal incidence between two transparent media, from the sheet conductivity of a 2D model, and the inverse, Re sigma with an error bar from a measured transmittance relative to the bare substrate.

When it refuses, and why

hamop raises an error with an explanation, instead of returning a doubtful number, when:

• a periodic model is given to dos, fermi_level, band_edges, sigma_optical, sigma_tensor, drude_weight or carrier_count without mesh or kpts;
• a function for finite models (with_peierls, principal_layers, kpm_dos, kpm_sigma, chern_marker) gets a periodic one, or a function that needs a k-grid (monkhorst_pack, symmetry_fold, find_point_group, fourier_interpolation) gets a finite one; berry_curvature, chern_number, band_curvatures, berry_dipole, quantum_weight and magnetic_supercell need a periodic 2D model;
• a function that does not handle overlap gets a model with overlap: kpm_sigma, chern_marker, the quantum-geometry functions, band_curvatures / berry_dipole, and solver="sparse" in the Berry functions;
• monkhorst_pack(time_reversal=True) is asked of a model with complex hopping blocks, where pairing k with -k is not guaranteed;
• an operation given to symmetry_fold does not map the lattice to itself, changes the eigenvalues at a random test k, or moves grid points off the grid;
• magnetic_supercell cannot represent the flux in its gauge; the message names an equivalent flux that works;
• principal_layers finds a coupling that skips a layer (the layers are thinner than the hopping range);
• the quantum geometry is asked at a k where the occupied and empty bands are closer than gap_min (the geometry is undefined there);
• a Berry link between neighbouring k-points vanishes ("refine the grid"), or the overlap matrix is not positive definite in the lowdin frame;
• a Wannier90 file is malformed (wrong counts, duplicate or missing elements, a short degeneracy list, an orbital index out of range), its Hamiltonian is not Hermitian, or the cell does not match the file; a tb.dat file lacks its position section, lists different R vectors in its two sections, or has a cell that cannot be cut to the number of periodic directions without changing its geometry;
• explicit k-point weights do not match the k-points in number, are negative, or do not add up to 1;
• the lead surface Green function (sancho_rubio, used by every transport function) fails its Dyson-equation check by both routes (decimation and decaying modes); this always happens at eta = 0 inside a lead band;
• the energy grid given to landauer_conductance, landauer_current or thermoelectric does not reach 30 kT beyond the chemical potentials, has a step larger than kT there, is not increasing, or T <= 0; thermoelectric also refuses when nothing is transmitted in the Fermi window;
• a planned or measured set of band energies cannot tell the fit parameters apart, or there are fewer measurements than parameters;
• propagate_uncertainty (linear method) is given a quantity that is not close to linear over one standard deviation of the parameters;
• sheet_optics gets a negative Re sigma (gain) or a refractive index that is not real and positive; sheet_conductivity_from_transmittance gets a transmittance ratio above what a passive sheet allows (by more than three times its stated error, or at all when no error is given);
• monkhorst_pack(time_reversal=True) is asked of a model with complex dipole blocks (set_dipole), where the optics at k and -k differ;
• a site index is negative or out of range, or a lattice image has the wrong number of components;
• a block has the wrong shape, an on-site or dipole block is not Hermitian, or a finite model is given a non-zero lattice image;
• a broadening eta is not positive (dos, sigma_optical; negative in sigma_tensor), a dephasing gamma or disorder strength W2 is negative, or n_occ exceeds the number of bands;
• photon energies given to sigma_optical or kpm_sigma are zero or negative (sigma_tensor accepts omega = 0, the static Hall conductivity of example 3), or the lineshape is unknown;
• the temperature is negative; drude_weight also refuses T = 0, and berry_dipole(method="fermi") needs kT > 0;
• lowest_bands is asked for the full spectrum, a filling or mu lies outside the spectrum, or the SCBA iteration does not converge.

How the results are checked

190 automated tests run on every change, on Python 3.9, 3.10, 3.11, 3.12, 3.13 and 3.14, and once more on Python 3.10 with the oldest NumPy (1.22.0) and SciPy (1.8.0) the package allows. The numerical tests compare the package with a closed-form result, a symmetry, or a second calculation done a different way; none compares against a number stored from an earlier run. The main checks, with the tolerances the tests use:

Bands, DOS, optics

• The chain's bands equal e0 + 2t cos ka to 1e-12; with overlap s they equal 2t cos ka / (1 + 2s cos ka) to 1e-12.
• The chain's DOS at the band centre equals 1/(2 pi |t|) to 5e-4, and the DOS integrates to the orbital count to 1e-6.
• Graphene's bands touch at K to 1e-9 and sit at -/+ 3|t| at the zone centre to 1e-12. Its optical conductivity at 1.0 and 1.3 eV is within 5 % of e^2/(4 hbar).
• At T = 0 a level exactly at mu is half occupied: the optics then equal the T = 1 K result to 1e-12 and half the result with mu in the gap, and carrier_count gives exactly 1 (spin 2 times 1/2). Wrong-length, negative or unnormalized k-point weights are refused by dos, sigma_optical, drude_weight, carrier_count and fermi_level.
• The two-site molecule absorbs at 2|t| (within 2e-3 eV) with the hand-derived peak height to 1e-3 relative; the Lorentzian lineshape also gives its hand-derived peak height to 1e-3.
• Shifting the energy zero with overlap (H -> H + cS) changes the conductivity, the tensor and the Drude weight by less than 1e-10 (relative, for the Drude weight).
• The Drude weight of the half-filled chain equals 8 spin |t| a to 1e-3 relative, and is below 1e-12 for an empty or full band.
• An on-site transition is exactly dark without a dipole block; with one, its peak matches the hand-derived value to 1e-3 relative, in the dense route. The same holds in a nonorthogonal basis (two uncoupled atoms give exactly twice one atom).

Topology and geometry

• The Haldane model's Chern number is 1 in the topological phase and 0 in the trivial phase, to 1e-12; it flips sign with the flux; all bands together give 0. The same integers come out with overlap and in both frames; the sparse solver gives 2 for two stacked copies, like the dense one.
• sigma_xy(0) equals the Chern number times e^2/h to 1e-6, sign included, and the tensor is antisymmetric to 1e-12.
• The SSH chain's Zak phases are 0 or pi to 1e-9 and differ by pi, with and without overlap.
• The Kane-Mele model equals two Haldane copies to 1e-12, and its gap at K is 6 sqrt(3) lambda_so to 1e-9.
• The Chern marker of a 10 x 10 flake sums to zero to 1e-8, and its bulk average is within 0.05 of the Chern number.
• The Berry curvature dipole is below 1e-10 with inversion symmetry and below 1e-12 for completely filled bands. For gapped graphene it at least halves when the mesh doubles, and one strained bond makes it more than 10 times larger. With one mirror line, the component along the mirror is below 1e-8 times the other one. The two integral forms agree within 10 %.
• The quantum metric and curvature match the two-band closed forms to 1e-4 relative; tr g >= 2 sqrt(det g) >= |Omega| holds at 25 random points per phase; the integrated Chern number is within 0.02 of the lattice one; tr K >= |C| holds in both phases.

Magnetic fields and spin

• A flux-threaded ring matches 2t cos((2 pi j + Theta)/N) to 1e-12; two gauges give the same spectrum to 1e-12; the flux through one square is exact to 1e-12; the spectrum repeats after one flux quantum to 1e-12.
• The lowest Landau level of the square lattice sits at -4|t| + hbar omega_c / 2 within 3 %.
• At zero flux the magnetic supercell equals folded bands to 1e-12; at half a flux quantum the square lattice matches +/- 2|t| sqrt(cos^2 kx + cos^2 ky) to 1e-12; at 1/3 the lowest band's Chern number matches sigma_xy to 1e-4.
• Spin doubling gives doublets to 1e-12, and a Zeeman term splits them by 2B to 1e-12. With overlap, the Zeeman-split bands equal (2t cos ka -/+ B)/(1 + 2s cos ka) to 1e-12.

Transport

• The lead surface Green function matches its closed form to 1e-6.
• A clean chain transmits 1 inside the band (to 1e-4) and 0 outside (to 1e-8); with overlap it still transmits 1 inside the band; two uncoupled chains transmit 2.
• One impurity matches the closed form to 1e-5.
• The recursive sweep equals direct inversion to 1e-10 (1e-12 with self-energies, and with overlap at eta = 0.05); sparse equals dense to 1e-12.
• With overlap, the lead surface Green function matches its closed form to 1e-12 at eta = 0.05, the Green function of the closed device (leads detached) equals (zS - H)^-1 to 1e-12, and i(G - G^dag) = G(GamL + GamR + 2 eta S)G^dag holds to 1e-12.
• The Buttiker probe at zero coupling equals the coherent result to 1e-12 and matches the single-site closed form to 1e-8; a probe on every layer gives a resistance linear in length (R^2 > 0.9999); current is conserved to 1e-12.
• SCBA with zero disorder equals the coherent result to 1e-12; its fixed point is converged below 1e-10 with Im Sigma <= 0; the central layer of a long chain matches the bulk SCBA equation to 1e-3.
• The closed-device LDOS equals the eigenvalue sum to 1e-13 (and to 1e-12 with overlap); bond currents obey Kirchhoff's law to 1e-7 and every cut carries T to 1e-6.
• principal_layers reproduces hand-built blocks exactly and refuses layers that are too thin.
• sancho_rubio refuses at eta = 0 inside the band; at eta = 0 outside the band it matches the real closed form to 1e-12, and at eta = 1e-12 inside the band the complex one to 1e-10. At the band centre and other lead resonances (E = 0, 8.9e-16, 1e-10, 1.0 eV) with eta = 1e-8 it matches the closed form at the same complex energy to 1e-8 for the chain, the chain with overlap and a two-site layer; the impurity chain's T(E) there matches its closed form to 1e-7, and the 30 K conductance on a grid containing E = 0 matches the one from the closed-form T(E) to 1e-7.
• Finite temperature: for a constant T(E) = tau the conductance is spin tau to 1e-10, the Seebeck coefficient is below 1e-12 V/K and the Lorenz ratio equals (pi^2/3)(k/e)^2 to 1e-9 relative at 30, 300 and 900 K (exact for constant T(E)); for a linear T(E) the Seebeck coefficient equals the Mott form and the Lorenz ratio its closed form, both to 1e-9 relative (exact for linear T(E)); the current through a constant T(E) is spin tau V to 1e-10 and odd in the bias. For the impurity chain at 30 K the conductance matches 2 T(mu) of the closed form to 1e-5 and I/V at a 0.1 mV bias matches G to 1e-6 relative; a band-edge step matches 2 [f(-2|t|) - f(2|t|)] to 1e-4 on a 1e-5 eV grid.

Large systems, k-grids, interpolation, files, fitting

• Sparse and dense matrices are identical; Lanczos matches the open chain's closed form to 1e-8 (with overlap too).
• KPM DOS: band-centre value within 1 %, integral within 0.5 %, zero outside the band with overlap. KPM conductivity: molecular line weight within 3 % of the closed form, and within 2 % of the dense route on a dimerized chain.
• The time-reversal fold reproduces DOS and conductivity to 1e-12 and the Drude weight to 1e-10; the six-fold fold of graphene reproduces the DOS to 1e-12 and the chemical potential to 1e-9, and wrong operations are refused. find_point_group finds 12, 8 and 2 operations for the hexagonal, square and chain lattices.
• Fourier interpolation reproduces the bands to 1e-12 for five models and flags a too-coarse grid.
• A hand-written graphene hr.dat reproduces graphene() bands to 1e-12 and its optical conductivity to 1e-10; a hand-written graphene tb.dat (Wannier90 number formats, a degeneracy-2 shell, centres off the plane) reproduces the bands to 1e-12 and both in-plane optical conductivities to 1e-7 relative (the file stores 8 significant digits), while moving the centres to the origin changes the optics by more than 1 %; a save/load round trip keeps each block to 1e-9 (the file stores 10 decimals); four kinds of corrupted file (too few lines, a duplicate element, a short degeneracy list, a non-Hermitian Hamiltonian) and a mismatched cell are refused.
• On the chain the fit covariance matches its closed form to 1e-6 relative; 300 seeded simulated experiments match the reported error bars within 15 %; a design that cannot separate the parameters is refused.
• Propagated error bars of band energies (linear in the parameters) match the closed form g^T C g to 1e-6 relative; across 200 seeded simulated experiments the scatter of the fitted band width matches its propagated error bar within 15 %, and the p-values of these correct-model fits fall below 0.1 in 4 % to 16 % of cases (10 % expected). The p-value equals the closed-form chi-squared tail for 2 and 10 degrees of freedom to 1e-12. The Monte Carlo method agrees with the linear one within 5 %; for a quantity with zero slope the linear method is refused and the Monte Carlo mean and spread match the chi-squared closed form within 6 % and 10 %.
• The sheet formulas match an independent thin-film (Airy) calculation of a film of vanishing thickness carrying the same sheet conductance to 1e-6, for five combinations of conductivity and refractive indices, with T + R + A = 1 to 1e-14. The fine-structure constant agrees with SciPy's CODATA value to 1e-8; free-standing graphene absorbs 2.3 %, as reported by Nair et al. (2008); the inverse recovers the film's conductivity to 1e-5, and its error bar matches 20000 seeded noisy measurements within 3 %.

What changed in 0.12.0

New

• propagate_uncertainty gives the error bar of any quantity computed from parameters fitted with fit_bands, and BandFit.p_value says whether the measured energies are consistent with the model (example 11).
• sheet_optics and sheet_conductivity_from_transmittance connect the sheet conductivity of a 2D model with an optical transmittance measurement, in both directions (example 12). They use the fine-structure constant (CODATA 2018).

Silent wrong answers that are now errors or fixed

• monkhorst_pack(time_reversal=True) checked only the hopping and overlap blocks. With a complex dipole block (set_dipole) the bands stay even in k but the optical matrix elements do not, so the folded grid gave a wrong optical conductivity: 1.949 instead of 2.096 at 1 eV (7 % low; 12 % at 2 eV) for the test model. Such a model is now refused, with a pointer to symmetry_fold for eigenvalue-only quantities.
• A negative site index, as in add_hop(0, -1, ...), passed the shape check but was placed outside the matrix, so the hopping vanished without an error: a two-site molecule had levels (0, 0) instead of (-1, +1). set_dipole had the same problem. Negative and out-of-range site indices are now refused.
• save_wannier90_hr with an empty comment, or one containing a line break, wrote a file that load_wannier90_hr could not read back. The header is now always one non-blank line.

Clearer errors

• Inputs that used to fail later with an unrelated message, or to give a meaningless number, are refused up front: a lattice image with the wrong number of components, norb < 1, positions and lattice vectors of different dimension, a mesh of the wrong length, a broadening eta that is not positive, a negative dephasing gamma or disorder W2, an SCBA mixing outside (0, 1], a Landauer window that is not positive, n_occ larger than the number of bands, and an R-vector with more than three components.
• sigma_optical accepts a single photon energy.
• Public functions now carry type hints.

What changed in 0.11.0

New

• Finite-temperature transport from T(E): landauer_conductance, landauer_current and thermoelectric (example 9).
• from_wannier90_tb reads a Wannier90 seedname_tb.dat file, taking the lattice vectors and the orbital centres from the file itself (example 10).

Silent wrong answers that are now errors or fixed

• sancho_rubio could return a wrong lead surface Green function without warning, in two ways. At eta = 0 inside a band the decimation never converges and the last iterate was returned: the chain lead at E = 0.5 gave 0.924 instead of 0.25 - 0.968i, and transmission gave T = 0 inside the band; this now raises an error. And when E equals an eigenvalue of the lead layer (the band centre E = 0 of the chain) the decimation loses up to log10(1/eta) digits and can settle on a wrong answer: at eta = 1e-8 it gave -6.7e7 i instead of -i, the impurity chain's T(0) came out 1.4e-15 instead of 0.862, and at the default eta = 1e-6 T(0) was off by 5e-6. Every result is now checked against the lead's Dyson equation and, if the check fails, recomputed from the lead's decaying modes.
• At T = 0, a level exactly at mu gave 0/0 and NaN in sigma_optical, sigma_tensor, carrier_count and friends. It is now half occupied (for the two-level test atom, carrier_count goes from NaN to 1.0).
• Explicit k-point weights were not checked. A shorter weight list silently dropped k-points (the chain's DOS at 2 eV on a 10-point grid came out 0.0 instead of 0.798 with only 3 weights), and weights summing to 2 doubled the result. Both are now refused.

Bugs fixed in 0.10.1

0.10.1 fixed six problems.

• Overlap doubled by on-site terms. In a model with overlap, every on-site block added without its own S_block added another identity to S. Adding an on-site energy or a Zeeman term in a second add_hop call (as the with_spin docstring suggests) therefore changed S and gave wrong bands. Now each site's on-site overlap is the identity unless you give it explicitly.
• Transport with overlap at finite eta. The NEGF functions built the reverse coupling block as the conjugate transpose of zS - H, which is conj(z) S^dag - H^dag rather than z S^dag - H^dag. The error is of order eta * S: negligible at the default eta = 1e-6 for transmission, but visible in device_greens and device_ldos at the larger eta often used for a local DOS.
• device_ldos ignored the inter-layer overlap (coup_S). With overlap the LDOS is -Im (G S)_ii / pi (the Mulliken convention, which shares the overlap between the orbitals involved), so leaving out part of S made the orbital LDOS not add up to the total DOS.
• with_spin dropped dipole blocks, so a spinful copy lost the on-site transitions set with set_dipole.
• carrier_count without mesh or kpts summed a periodic model at k = 0 only (2.0 instead of 1.0 for the half-filled chain). It now refuses, like the other functions.
• Temperature. drude_weight returned NaN at T = 0 and a negative weight at T < 0; the other functions accepted T < 0. Both are now refused.

Earlier README and CHANGELOG text also overstated some checks; the 0.10.1 entry of CHANGELOG.md lists the corrections. The full history is in the CHANGELOG.

Limits

These are deliberate choices, not oversights:

• The velocity uses the site-diagonal position approximation; the intra-atomic part enters only through the dipole blocks you supply.
• sigma_optical sums transitions between bands only (pairs closer than 1 meV are left out); add drude_weight for metals.
• symmetry_fold is valid for results that depend on eigenvalues only (DOS, chemical potential, band edges, carrier count), not for direction-dependent ones such as sigma_xx or the Drude weight.
• The SCBA is elastic (disorder only) and has no vertex corrections (the next-order correction to the disorder-averaged conductance); inelastic electron-phonon scattering is not implemented. Electron-electron interaction enters only as a self-energy you supply; there is no mean-field or GW self-consistency.
• bond_currents works in an orthogonal basis only (it takes no overlap arguments).
• kpm_sigma covers orthogonal, finite models and the longitudinal response only. A KPM Hall conductivity of a finite system is not offered: in this position convention it is zero for any bounded system (Im Tr[PxQy] = 0, itself a test), so chern_marker is the finite-system alternative.
• fourier_interpolation is the Fourier step of Wannier interpolation; it does not build maximally localized Wannier functions (Marzari and Vanderbilt, Phys. Rev. B 56, 12847 (1997), is not implemented).
• The quantum-geometry functions and the Berry curvature dipole handle orthogonal bases only.
• landauer_conductance, landauer_current and thermoelectric take T(E) as given: a current at finite bias uses the same T(E) for every bias unless you recompute it (no self-consistent potential drop), and the heat conductance is the electrons' part only (no phonons). A step in T(E) inside the Fermi window (a band edge) converges only linearly in the grid step.
• from_wannier90_tb uses the diagonal position elements (the centres) and ignores the off-diagonal ones, in line with the site-diagonal position approximation above. Both Wannier90 readers divide every element by the degeneracy listed in the file; files written with write_ndegen_applied = true have not been tested.
• The lead surface Green function needs eta > 0 inside a lead band. Its check (the Dyson-equation residual below 1e-10) bounds the error only up to the conditioning of that equation, which grows as eta shrinks; at eta = 1e-8 the tests find errors below 1e-8, smaller than the shift the broadening itself causes.
• propagate_uncertainty and BandFit.p_value assume, like the fit's covariance, independent Gaussian measurement errors. The p-value tests the model and the stated errors together; it cannot say which of the two is wrong.
• sheet_optics treats normal incidence on a sheet much thinner than the wavelength, between transparent media. The back face of a substrate and multiple reflections inside it are not included (they largely cancel in T/T_bare). The sheet conductivity unit e^2/(4 hbar) is a sheet conductance only for a model periodic in two directions. A transmittance alone cannot separate Re sigma from Im sigma; the inverse assumes the Im sigma you give (0 by default), which enters only at second order.
• The package ships no material constants. The only physical constants in the code are the Boltzmann constant and the fine-structure constant, both CODATA 2018 values (Tiesinga, Mohr, Newell and Taylor, Rev. Mod. Phys. 93, 025010 (2021)).

Relation to existing tools

Excellent tools cover parts of this space: PythTB and pybinding build tight-binding models and their spectra, and Kwant is the standard for quantum transport. hamop does not replace any of them, and for their core use cases they are more capable. Its niche is the combination they leave open: overlap matrices handled across the observables listed above, optics computed from the same Bloch matrices as the spectrum, and a small NumPy/SciPy-only core checked against closed forms.

Where it comes from

Methodological basis:

"Learning the quantum Hamiltonian of defective monolayer MoS2 reveals collective vacancy brightness decoupled from defect count"; code for the paper: https://github.com/Tanvir-Mahmud-Mahim/mos2-vacancy-optics

That study computes the optics, the electronic structure and the transport of vacancy-disordered MoS2 supercells from one density-functional Hamiltonian, so that a defect configuration's optical and electronic signatures are consistent, and its conclusions depend on that consistency. This package is the general-purpose engine distilled from that pipeline: the same observables for any Hamiltonian a user supplies, with the material-specific machinery (DFT extraction, machine-learned Hamiltonians, MoS2 structures) left in the paper repository.

Citing, support and license

If hamop helps your work, please cite it with the concept DOI 10.5281/zenodo.22311381, which always resolves to the latest release; every release is archived on Zenodo. Citation metadata is in CITATION.cff.

The package is written and maintained by Tanvir Mahmud Mahim (Department of Electrical and Electronic Engineering, BRAC University), who reviews every change and takes the final decision on scope and releases. There is no separate governance body; design questions are discussed in the open in issues and pull requests, and the standing rule of CONTRIBUTING.md binds the maintainer exactly as it binds contributors: a change that touches physics arrives with a test, and a constant arrives with its source.

Support runs through the issue tracker at https://github.com/TaN-MM-Org/hamop/issues. Usage questions are welcome there alongside bug reports; a docstring that left a unit or a sign convention unclear is treated as a documentation bug, not as user error. The maintainer aims to respond within a week.

While the version is below 1.0 the API may still move between minor versions; such changes are called out in the release notes.

Licensed under Apache-2.0 (see LICENSE).