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, a test of whether the data need a richer model, the chemical potential and quantum capacitance that a measured carrier density implies, and the transmittance of a 2D sample that an optical measurement would see;
• the conductivity, the thermopower (Seebeck coefficient) and the electronic heat conductivity of a bulk crystal or a gated 2D sheet, in the relaxation-time (Boltzmann) picture.

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)
• Tutorials (three complete analyses, from data to a result with error bars)
• What is in the package
• When it refuses, and why
• How the results are checked
• What changed in 0.14.0
• What changed in 0.13.0
• 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). Since 0.14.0 a model can also repeat in fewer directions than its space has, such as a ribbon that repeats along x but has a finite width in y (example 14).
• 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.
• Relaxation time tau, Boltzmann transport -- in a bulk crystal or a large sheet an electron scatters, on average, after a time tau. In the Boltzmann picture the conductivity, the Seebeck coefficient and the heat conductivity of the electrons follow from the band velocities and tau. The size of tau drops out of the Seebeck coefficient, but how tau changes with energy does not (tutorial 3).
• 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).
• Bulk transport (boltzmann_transport) is in SI units: the conductivity in S/m for a 3D model, in S (per square) for a 2D model and in S m for a chain or a ribbon; the Seebeck coefficient in V/K; tau in seconds.
• 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.14.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.

13. Conductivity and thermopower of a doped crystal

import numpy as np
from hamop import graphene, boltzmann_transport

g = graphene(t=-2.7, a=2.46)
mu = np.array([0.3, 0.4, -0.3])               # chemical potentials, eV
r = boltzmann_transport(g, mu, T=300.0, mesh=300, tau=1e-14)   # tau in s

kB = 8.617333262e-5                           # Boltzmann constant, eV/K
for m, sig, S, L in zip(mu, r["sigma"], r["seebeck"], r["lorenz"]):
    cone = -(np.pi ** 2 / 3) * kB * (kB * 300.0) / m   # Dirac cone, V/K
    print(f"mu = {m:+.1f} eV: sigma = {sig[0, 0] * 1e3:.3f} mS, "
          f"S = {S[0, 0] * 1e6:+.2f} microvolt/K (cone {cone * 1e6:+.2f}), "
          f"Lorenz ratio {L[0]:.3e} V^2/K^2")

mu = +0.3 eV: sigma = 0.353 mS, S = -24.43 microvolt/K (cone -24.43), Lorenz ratio 2.380e-08 V^2/K^2
mu = +0.4 eV: sigma = 0.471 mS, S = -18.35 microvolt/K (cone -18.32), Lorenz ratio 2.407e-08 V^2/K^2
mu = -0.3 eV: sigma = 0.353 mS, S = +24.43 microvolt/K (cone +24.43), Lorenz ratio 2.380e-08 V^2/K^2

boltzmann_transport gives the conductivity, Seebeck coefficient and electronic heat conductivity of a bulk crystal or a large 2D sheet in the relaxation-time approximation of the Boltzmann equation, the picture used by the BoltzTraP code (Madsen and Singh, Comput. Phys. Commun. 175, 67 (2006)). Here graphene is doped to 0.3 and 0.4 eV (electrons) and to -0.3 eV (holes), with a constant relaxation time of 10 fs. Three checks against the Dirac cone, all exact for a cone at this doping: the conductivity equals (e^2/h)(2 mu tau / hbar) (0.353 and 0.471 mS); the Seebeck coefficient equals -(pi^2/3)(k_B/e)(kT/mu) (to the printed digits at 0.3 eV, within 0.2 % at 0.4 eV, where the 300 x 300 mesh limits the accuracy), and it changes sign with the carriers; and the Lorenz ratio is below (pi^2/3)(k_B/e)^2 = 2.443e-8 V^2/K^2 by the factor 1 - (pi^2/3)(kT/mu)^2. A constant tau cancels in S and in the Lorenz ratio; tau may also be a function of energy (tutorial 3).

14. A ribbon: periodic in one direction of the plane

import numpy as np
from hamop import TightBindingModel, bands, dos

# A strip of W = 3 rows of the square lattice (spacing 1 Angstrom,
# hopping t = -1 eV), repeated along x only: one lattice vector in a plane.
W, t = 3, -1.0
strip = TightBindingModel(positions=[[0.0, float(y)] for y in range(W)],
                          norb=1, cell=[[1.0, 0.0]])
for y in range(W):
    strip.add_hop(y, y, (1,), [[t]])                 # along the strip
    if y + 1 < W:
        strip.add_hop(y, y + 1, (0,), [[t]])         # across it

k = np.array([[0.0, 0.0], [np.pi / 2, 0.0]])         # Cartesian k, 1/Angstrom
exact = [sorted(2 * t * np.cos(kx) + 2 * t * np.cos(p * np.pi / (W + 1))
                for p in range(1, W + 1)) for kx in k[:, 0]]
print("bands agree with the closed form to 1e-12 eV:",
      np.abs(bands(strip, k) - exact).max() < 1e-12)
print("cell length:", strip.cell_volume, "Angstrom; reciprocal vector:",
      np.round(strip.reciprocal_vectors, 4))
rho = dos(strip, [-2.0, 0.0, 2.0], mesh=600, eta=0.05)
print("DOS at -2, 0, 2 eV:", np.round(rho, 4), "states/eV per cell")

bands agree with the closed form to 1e-12 eV: True
cell length: 1.0 Angstrom; reciprocal vector: [[6.2832 0.    ]]
DOS at -2, 0, 2 eV: [0.7805 0.6105 0.7805] states/eV per cell

A model that repeats in fewer directions than its space has (a ribbon, a nanotube, a chain lying in a plane) is built with fewer lattice vectors than Cartesian components, here one vector in two dimensions. The bands of this strip are the standing waves across it, -2|t| cos(kx) - 2|t| cos(p pi / (W + 1)), p = 1..W. Wave vectors stay Cartesian: monkhorst_pack places them along the lattice vector, and reciprocal_vectors gives the matching reciprocal vector. The velocity across the strip enters the optics (sigma_optical with direction=1 gives the transitions between the standing waves), while nothing conducts across it in boltzmann_transport. Before 0.14.0 such a model was accepted, but every function given mesh= failed with a NumPy error.

Tutorials

Three scripts in examples/ carry a complete analysis from data to a result with error bars. Each uses only hamop, makes its own synthetic "measured" data with a fixed seed (replace them with yours), runs in seconds and prints the result below. The test suite runs all three and checks every printed number. Since 0.14.0 the scripts are also in the source distribution, so the check runs from an unpacked sdist too.

Gated graphene: from a Hall density to the absorption edge

examples/tutorial_gated_graphene.py starts from eight Hall readings of the electron density of a gated graphene sheet. It averages them, turns the density into the chemical potential mu with chemical_potential_from_density (the interval comes from solving for the density minus and plus its error), reports the quantum capacitance, and computes the optical conductivity and the transmittance relative to the bare substrate around the Pauli-blocking edge at 2 mu. Below the edge the transitions are blocked because their final states are already filled; above it the absorption returns to the universal value. The numbers in brackets are the conductivity at the two ends of the mu interval.

Hall density: n = (5.980 +/- 0.047) x 1e12 cm^-2
chemical potential: mu = 0.2448 eV, interval [0.2438, 0.2458] eV
Pauli-blocking edge: 2 mu = 0.490 +/- 0.002 eV
quantum capacitance: 7.57 microfarad/cm^2
photon (eV)  Re sigma / (e^2/4 hbar)   T/T_bare (substrate n = 1.5)
   0.30      0.028 [0.027, 0.029]           0.9995
   0.45      0.326 [0.318, 0.334]           0.9940
   0.60      0.888 [0.884, 0.892]           0.9839
   0.90      1.014 [1.014, 1.014]           0.9817

Band fit with a model check

examples/tutorial_band_fit_model_check.py fits 20 measured band energies of a chain. The nearest-neighbor model fails its p-value test. compare_models asks whether a second-neighbor hopping t2 is needed and answers yes, with a p-value of 2e-53. The extended model fits well (p = 0.93) and recovers the true parameters within their error bars. The last lines show why the check matters: at k = pi, a point the measurement did not reach, the nearest-neighbor model predicts the band top 0.17 eV too low, about 14 times its own error bar, while the extended model is within its error bar of the truth.

step 1, nearest-neighbor fit: chi2 = 245.8 for 18 dof, p = 5.8e-42
step 2, is t2 needed? Delta chi2 = 236.5 for 1 dof, p = 2.3e-53
step 3, extended fit: chi2 = 9.3 for 17 dof, p = 0.93
   t  = -0.9991 +/- 0.0049 eV   (truth -1)
   e0 = +0.1966 +/- 0.0067 eV   (truth +0.2)
   t2 = +0.0749 +/- 0.0049 eV   (truth +0.08)
step 4, quantities not measured directly:
   band width 4|t|      = 3.9962 +/- 0.0194 eV   (truth 4)
   band top, k = pi     = 2.3446 +/- 0.0156 eV   (truth 2.36)
   same, nearest-neighbor model only: 2.1886 +/- 0.0118 eV

Thermopower: what scatters the electrons?

examples/tutorial_thermopower_scattering.py takes six gate densities of a graphene sheet, each with a Hall density (2 % error) and a Seebeck coefficient S (error 1.5 microvolt/K) at 300 K. In the Boltzmann picture S does not depend on the size of the relaxation time, but it does depend on how the relaxation time changes with energy. The script models that change as tau ~ |E|^p and finds p from the data: it turns each density into mu with chemical_potential_from_density, computes S(mu) with boltzmann_transport, adds the spread of S over each mu interval to the error of that point, and fits p by weighted least squares. The data were made with p = 1, for which the conductivity of the Dirac cone grows linearly with the density. The fit recovers p = 1.00 +/- 0.02 with an acceptable chi-squared, and the constant-tau model (p = 0), which predicts about half the measured thermopower, is ruled out.

n (1e12/cm^2)   mu (eV)   S measured     S fit   S, constant tau
   2.00      0.1364    -94.0 +/- 1.5    -94.9        -50.3
   3.08      0.1727    -81.3 +/- 1.5    -80.2        -43.1
   4.10      0.2008    -67.8 +/- 1.5    -68.2        -35.3
   4.95      0.2218    -66.4 +/- 1.5    -63.9        -33.2
   6.46      0.2546    -52.5 +/- 1.5    -55.1        -28.4
   7.92      0.2827    -50.8 +/- 1.5    -50.8        -26.0
(S in microvolt/K; negative: electrons)
scattering exponent p = 1.00 +/- 0.02   (tau ~ |E|^p)
fit: chi2 = 5.8 for 5 dof, p-value 0.33
constant tau (p = 0): chi2 = 2779 for 6 dof, p-value 10^-597

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; reciprocal_vectors (new in 0.14.0) gives the reciprocal lattice vectors. The cell may have fewer lattice vectors than Cartesian dimensions (a ribbon, example 14); lattice vectors that are linearly dependent are refused.
• 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. Where bands are degenerate at a k-point it sums over the degenerate pairs, so the result does not depend on which eigenvectors the solver returns.
• 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.
• compare_models (tutorial 2): fits the same data with a simple model and with an extended model that contains it, and says whether the extra parameters are needed. With measurement errors it uses the likelihood-ratio (chi-squared) test (Wilks, Ann. Math. Statist. 9, 60 (1938)); without them, the F test of linear regression. Both are exact when the band energies are linear in the parameters. An optional embed function, which places the simple parameters inside the extended model, lets it verify that the models are nested.

Carrier density and gating

• carrier_density -- the net density of added electrons (negative for holes) at a chemical potential, per cm^2 for a 2D model (per cm in 1D, per cm^3 in 3D), counted from a neutral filling you give.
• chemical_potential_from_density (tutorial 1) -- the inverse: the chemical potential that holds a measured density (from the Hall effect or from the gate voltage), with the interval that the density's error bar maps to, dn/dmu and the quantum capacitance e^2 dn/dmu in F/cm^2. For graphene these follow the Dirac-cone forms n = mu^2 / (pi (hbar v_F)^2) and C_Q = 2 e^2 |mu| / (pi (hbar v_F)^2). Graphene's quantum capacitance was measured by Xia et al., Nature Nanotechnology 4, 505 (2009).

Bulk transport (new in 0.14.0)

• boltzmann_transport -- the conductivity tensor, the Seebeck tensor, the electronic thermal conductivity at zero current and the Lorenz ratio of a periodic model, in the relaxation-time approximation of the Boltzmann equation (Madsen and Singh, Comput. Phys. Commun. 175, 67 (2006)), for one chemical potential or many at once (example 13, tutorial 3). The relaxation time is a constant or a function of the state's energy. The band velocities are the same velocity matrix elements as in the Kubo functions, so overlap models are handled, and where bands are degenerate the sum runs over the degenerate pairs, as in drude_weight. With a constant relaxation time the conductivity equals the Drude weight times tau / (pi hbar).

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;
• kpm_dos or kpm_sigma is given bounds that do not enclose the spectrum (the Chebyshev moments then grow without limit; before 0.13.0 the result was silently wrong, a DOS of -4.7e7 in the test);
• principal_layers gets a layer width that is not a positive length;
• compare_models is given an extended model without more parameters, too few data for it, an embed that does not reproduce the simple model, or an extended fit that ends worse than the simple one (the models are then not nested, or the fit stopped in another minimum);
• chemical_potential_from_density is asked for a density the bands cannot hold, at T = 0, or on a k-mesh so coarse for the temperature that one state carries more than 5 % of dn/dmu (then mu would jump from one mesh state to the next instead of following the density); carrier_density and its inverse refuse finite models;
• boltzmann_transport gets a finite model, T <= 0, a relaxation time that is not positive and finite for every state, a mu with no state within 60 kT (far outside the bands or deep inside a wide gap), a mu at which nothing conducts (a flat band), or a k-mesh so coarse for the temperature that one state carries more than 5 % of the conductivity;
• symmetry_fold is given a mesh that an operation maps off itself (since 0.14.0; for example the four-fold rotation on an (8, 4) mesh, which before gave wrong weights), or a model with fewer lattice vectors than dimensions (also refused by find_point_group);
• TightBindingModel is given lattice vectors that are linearly dependent or more numerous than the Cartesian dimensions (since 0.14.0; before, the first k-grid failed with a NumPy error);
• bands is asked for a model whose number of bands changes with k (overlap directions dropped at some k-points), which cannot form one array; dos, fermi_level, band_edges, carrier_count and carrier_density handle such models.

How the results are checked

253 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. At graphene's Dirac point, where the two bands are degenerate, it equals the eigenvector-free form Tr[(dH/dk)^2] (-df/de) to 1e-12 in both directions, and on a mesh through K and K' D_xx = D_yy to 1e-12, as the hexagonal symmetry requires.
• The routines that loop over k-points (bands, dos, fermi_level, carrier_count, sigma_optical, sigma_tensor, drude_weight) agree with the 0.12.1 one-k-point-at-a-time code, kept in the tests as the reference, to 1e-12 relative for six models (orthogonal, overlapping, complex, spinful, with dipole blocks, and one whose overlap matrix drops a state at every k-point).
• 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 block versions of berry_curvature (and so chern_number), band_curvatures, berry_dipole and quantum_weight (0.14.0) agree with the 0.13.0 point-by-point code, kept in the tests as the reference: plaquette fluxes to 1e-10 for six models in both frames (Haldane in both phases, with overlap, Kane-Mele, all bands, and a Hofstadter cell with three filled bands; measured differences are below 1e-14), the energies at the plaquette centers to 1e-12 and the curvature dipole (method="fermi") to 1e-10 relative, and the quantum-weight tensor to 1e-9 relative (its finite differences amplify rounding; measured 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. A chain of spacing 1.42 Angstrom cut with width 1.42 gives 40 one-site layers (width 2.84: 20 two-site layers), and the impurity chain built that way matches its closed form to 1e-5.
• 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. The in-place Chebyshev recurrence equals the 0.12.1 one to 1e-12 (exact trace) and bit for bit (stochastic trace); bounds that cut the spectrum are refused. 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. Operations given as a generator are checked like a list. On an uneven mesh, (8, 4) or (6, 3), a group that moves grid points off the grid is refused, and the two axis mirrors of the square lattice, which fit, reproduce the full-grid DOS to 1e-12. The array versions of the orbit fold and of the time-reversal fold (0.14.0) give the points, order and weights of the 0.13.0 loops (weights exactly, k-points to 1e-13).
• A model whose overlap matrix drops a state at one k-point only: fermi_level returns the root of carrier_count to 1e-9, band_edges the extremes of a k-point-by-k-point loop exactly, and carrier_density and its inverse are consistent to 1e-8.
• 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 %.
• compare_models: for the chain with and without a second-neighbor hopping (energies linear in the parameters) the chi-squared and F statistics equal the least-squares closed forms computed with numpy.linalg.lstsq to 1e-6 relative, and the p-values the chi-squared and F tail functions to 1e-6; across 60 seeded simulated experiments without t2 the statistic averages 1 (within 0.45; it came out 1.25) and p < 0.1 occurs 1 to 13 times (11); with t2 = 0.1 eV the p-value is below 1e-10 and the fitted t2 is within 4 error bars of the truth.
• carrier_density: the chain matches the exact integral of its density of states (adaptive quadrature) to 1e-9 at five chemical potentials, and graphene matches its closed-form dispersion summed on the same grid to 1e-10, with electron-hole symmetry to 1e-10. The inverse returns the chemical potential to 1e-9, its interval ends reproduce the density minus and plus its error to 1e-6 of the error, dn/dmu equals a finite difference of the density to 1e-6 and the chain's closed form to 1e-3. For graphene at 0.2 and 0.4 eV (300 K) the density matches the Dirac cone (by quadrature) and the quantum capacitance matches 2 e^2 mu / (pi (hbar v_F)^2) within 1 %; the lattice bands deviate from the cone by 0.1 % and 0.4 %. The elementary charge equals SciPy's exact SI value.

Bulk transport and ribbons (0.14.0)

• boltzmann_transport for the chain matches an adaptive quadrature of the closed-form dispersion (no eigenvectors, no mesh) to 1e-9 for the conductivity, the Seebeck coefficient and the Lorenz ratio, with a constant and with an energy-dependent relaxation time (measured 1e-13). Graphene at +/-0.3 eV and 300 K matches the Dirac-cone closed forms within 1 % (measured deviations below 2.5e-3): the conductivity e^2 tau v_F^2 N(mu) / 2, the Seebeck coefficient -(pi^2/3)(k_B/e)(kT/mu), the Lorenz ratio (pi^2/3)(k_B/e)^2 [1 - (pi^2/3)(kT/mu)^2], and, for tau ~ |E|, the Seebeck coefficient of a conductivity that grows as E^2. With a constant relaxation time the conductivity equals the Drude weight times tau / (pi hbar) to 1e-10, at the Dirac point too; rotating an anisotropic crystal rotates all three tensors to 1e-9; a constant relaxation time scales the conductivity and drops out of S and the Lorenz ratio to 1e-12. The Planck constant equals SciPy's exact SI value.
• A chain lying at 30 degrees in a plane gives the DOS, chemical potential and carrier density of linear_chain to 1e-12, and its Drude weight and Boltzmann conductivity projected on x and y are the chain's times cos^2 and sin^2 of the angle. A strip of the square lattice matches its closed-form bands to 1e-12 and, built a second way as a 2D model with a vacuum cell, gives the same DOS, optics and Boltzmann results to 1e-12.

What changed in 0.14.0

New

• boltzmann_transport: the conductivity, Seebeck coefficient and electronic heat conductivity of a bulk crystal or a gated sheet, with a constant or an energy-dependent relaxation time, for many chemical potentials at once (example 13). Tutorial 3 uses it to find how the relaxation time changes with energy from thermopower data.
• Models periodic in fewer directions than their space, such as ribbons (example 14). TightBindingModel.reciprocal_vectors gives the reciprocal lattice vectors.
• The tutorial scripts are part of the source distribution, so their check runs from an unpacked sdist as well.

Faster

• chern_number, berry_curvature, band_curvatures, berry_dipole and quantum_weight assemble and diagonalize all k-points at once. Measured with one thread (median of 3 runs; the lower of two measurements on a shared machine): chern_number for the Haldane model on a 48 x 48 mesh 98 times faster, berry_dipole on a 48 x 48 mesh 140 times (method="grad") and 160 times (method="fermi"), quantum_weight on an 18 x 18 mesh 120 times. The results agree with the old loops to rounding (1e-12 or better).
• symmetry_fold and monkhorst_pack(..., time_reversal=True) fold the grid with integer array operations: 10 times faster for graphene's 12 operations on a 60 x 60 mesh, 17 times for the time-reversal fold of a 300 x 300 mesh. The k-points and weights are those of 0.13.0 (the weights exactly, the k-points to rounding).

Fixed

• symmetry_fold with the operations given as a generator used the generator up before the spectral check, so a non-symmetry was folded without being checked. Graphene with one strained bond and the three-fold rotation gave a DOS of 0.2636 instead of 0.2393 at -3 eV (10 % off). It is now refused, as it was for a list.
• symmetry_fold on an uneven mesh, such as (8, 4), mapped the grid indices instead of the grid points, so an operation that moves points off the grid was not noticed and the weights were wrong: the square lattice's DOS at 0.5 eV came out 0.0422 instead of 0.0729, and graphene's at 0.5 eV on a (6, 3) mesh 0.0000 instead of 0.0737. Such a fold is now refused, and a group that fits the mesh gives the full-grid result. Even meshes are unchanged.
• fermi_level, band_edges and carrier_density failed with a NumPy error ("inhomogeneous shape") for a model whose overlap matrix drops a state at some k-points only; dos and carrier_count handled it. They now count the states one by one. bands, which must return one array, refuses such a model with an explanation.
• A model with fewer lattice vectors than dimensions (a ribbon) was accepted, but monkhorst_pack, cell_volume and every function given mesh= failed with a NumPy error. It is now supported (see New); lattice vectors that cannot form a lattice are refused when the model is built.

What changed in 0.13.0

New

• carrier_density and chemical_potential_from_density connect a measured carrier density (Hall effect, gate voltage) with the chemical potential that the optics and transport functions take, with an error interval, dn/dmu and the quantum capacitance (tutorial 1).
• compare_models tests whether measured band energies need a richer model, by the likelihood-ratio or F test for nested models (tutorial 2).
• Two tutorials in examples/, each a complete analysis with error bars; the tests run them and check their output.

Faster

• bands, dos, fermi_level, band_edges, carrier_count, sigma_optical, sigma_tensor and drude_weight now assemble and diagonalize the Bloch matrices for many k-points at once instead of one at a time. The results agree with the old loop to 1e-12 (the order of floating-point sums differs). Measured with one thread (median of 5 runs; the lower of two measurements on a shared machine): sigma_optical for graphene on a 120 x 120 mesh 140 times faster, sigma_tensor on a 48 x 48 mesh 270 times, dos 26 times, bands at 2000 k-points 100 times, carrier_count 110 times. Small cells gain the most; for large cells the matrix algebra dominates either way. The whole test suite runs in under half the time.
• kpm_dos updates its Chebyshev vectors in place: 2.5 times faster for a 1000-site chain with the exact trace, same result to 1e-15.

Fixed

• drude_weight at a k-point where bands are degenerate (a Dirac point on the mesh) summed only the diagonal of whatever eigenvectors the solver returned inside the degenerate pair, so the result was arbitrary there. For graphene at mu = 0 with only K on the grid it gave 3068.8 along x but 0.0 along y; on a 30 x 30 mesh D_xx and D_yy came out 3.538 and 3.281 (7 % apart, where symmetry makes them equal). It now sums over the degenerate pairs, the zero-frequency limit of the Kubo formula: 3068.8 in both directions, and 6.820 for both on the mesh. Away from degeneracies nothing changes (1e-12).
• principal_layers assigned a site that sits exactly on a layer boundary by floating-point division, which can land just below the integer (15 x 1.42 / 1.42 = 14.999999999999998). A chain of spacing 1.42 Angstrom cut with width 1.42 was refused as "layers 14 and 16 couple", and with width 2.84 one layer got three sites and the next one. Sites within 1e-9 layer widths of a boundary now count on it.

Clearer errors and documentation

• kpm_dos and kpm_sigma refuse bounds that do not enclose the spectrum (before: silently wrong, a DOS of -4.7e7 for a 200-site chain with bounds (-1.9, 1.9)); principal_layers refuses a layer width that is not a positive length.
• sigma_tensor and monkhorst_pack now say that a grid folded by time reversal is not valid for the Hall component: gapped graphene gives sigma_xy(0) = -0.912 on such a grid instead of 0. Use the full grid for directions=(0, 1).

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.
• carrier_density and chemical_potential_from_density assume rigid bands: doping or gating only moves mu and does not change the Hamiltonian (no self-consistent electrostatics, no screening), and the density is counted from the neutral filling you state. The quantum capacitance is e^2 dn/dmu of the model's bands, without disorder broadening.
• compare_models assumes independent Gaussian errors and nested models. Its p-value is exact for energies linear in the parameters and holds for small errors otherwise; it does not apply when the simple model sits on the edge of the extended model's allowed range (a parameter that must be positive, tested at zero).
• sigma_tensor with a != b (the Hall component) needs the full k-grid; a grid folded by time reversal gives a wrong result.
• 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.
• boltzmann_transport uses the relaxation-time approximation: one relaxation time per state that depends on the state's energy only (not on its direction or on k), no scattering back into the state (no vertex corrections), and only the electrons' heat conduction (no phonons, no phonon drag). It needs hbar / tau small against kT and the band spacing, and it does not give the Hall coefficient. On a model with fewer lattice vectors than dimensions the Seebeck tensor is taken in the conducting directions only.
• Models with fewer lattice vectors than dimensions (ribbons) work with the bands, DOS, optics, carrier and Boltzmann functions, berry_phase, the time-reversal fold and Fourier interpolation. berry_curvature, chern_number, band_curvatures, berry_dipole, quantum_weight and magnetic_supercell need a 2D model with two lattice vectors, and symmetry_fold and find_point_group need as many lattice vectors as dimensions; they refuse a ribbon. The optical conductivity of a ribbon is per unit length (unit e^2/(4 hbar) times Angstrom), as for a 1D chain.
• The package ships no material constants. The only physical constants in the code are the Boltzmann constant, the fine-structure constant, (since 0.13.0, for the quantum capacitance) the elementary charge and (since 0.14.0, for the Boltzmann coefficients) the Planck constant, all CODATA 2018 values (Tiesinga, Mohr, Newell and Taylor, Rev. Mod. Phys. 93, 025010 (2021)); the elementary charge and the Planck constant are exact in the SI.

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).