Metadata-Version: 2.4
Name: hamop
Version: 0.12.0
Summary: One tight-binding Hamiltonian, every observable: bands, DOS, Kubo optical conductivity and NEGF transmission from the same real-space blocks, nonorthogonal bases included
Author-email: "Tanvir M. Mahim" <tanvir.mahim@bracu.ac.bd>
License: Apache-2.0
Project-URL: Homepage, https://tanvir-mahmud-mahim.github.io/software/
Project-URL: Repository, https://github.com/TaN-MM-Org/hamop
Project-URL: Issues, https://github.com/TaN-MM-Org/hamop/issues
Project-URL: Changelog, https://github.com/TaN-MM-Org/hamop/releases
Keywords: tight binding,Kubo-Greenwood,NEGF,Landauer,optical conductivity,quantum transport,LCAO
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.8
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# hamop

[![Tests](https://github.com/TaN-MM-Org/hamop/actions/workflows/ci.yml/badge.svg)](https://github.com/TaN-MM-Org/hamop/actions/workflows/ci.yml)
[![PyPI](https://img.shields.io/pypi/v/hamop?label=PyPI&color=blue&cacheSeconds=3600)](https://pypi.org/project/hamop/)
[![License](https://img.shields.io/badge/License-Apache_2.0-blue.svg)](LICENSE)
[![DOI](https://img.shields.io/badge/DOI-10.5281%2Fzenodo.22311381-blue)](https://doi.org/10.5281/zenodo.22311381)

`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](#a-short-guide-to-the-words-used-here)
- [Install, units and conventions](#install-units-and-conventions)
- [Examples](#examples) (each with the output it prints)
- [What is in the package](#what-is-in-the-package)
- [When it refuses, and why](#when-it-refuses-and-why)
- [How the results are checked](#how-the-results-are-checked)
- [What changed in 0.12.0](#what-changed-in-0120)
- [What changed in 0.11.0](#what-changed-in-0110)
- [Bugs fixed in 0.10.1](#bugs-fixed-in-0101)
- [Limits](#limits)
- [Relation to existing tools](#relation-to-existing-tools)
- [Where it comes from](#where-it-comes-from)
- [Citing, support and license](#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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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)

```python
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

```python
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

```python
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](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](https://www.physics.rutgers.edu/pythtb/) and [pybinding](https://docs.pybinding.site/) build tight-binding models and their spectra, and [Kwant](https://kwant-project.org/) 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](https://doi.org/10.5281/zenodo.22311381),
which always resolves to the latest release; every release is
archived on Zenodo. Citation metadata is in
[CITATION.cff](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](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](LICENSE)).
