Metadata-Version: 2.4
Name: hamop
Version: 0.11.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.

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

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

## 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;
- 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;
- 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

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

## 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.
- The package ships no material constants. The only physical constant
  in the code is the Boltzmann constant (CODATA 2018).

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