Harrison Tight-Binding (universal parameters)

Building a solid's electronic band structure normally means DFT or an SCF calculation -- this module skips both. Harrison's tight-binding model builds the Hamiltonian from two universal ingredients only: tabulated free-atom orbital energies, and one bond-scaling law with the same four coefficients for every element pair. No fitting per material -- the tradeoff for that convenience is accuracy (see the Accuracy section below).

Step 1. A real Si-Si bond, eigenvalues

import numpy as np
import dense_evolution as de

H = de.sp3_dimer_hamiltonian('Si', 'Si', bond_length_angstrom=2.35)
np.round(np.linalg.eigvalsh(H), 3)

array([-16.626, -12.25 ,  -9.847,  -7.638,  -7.638,  -5.402,  -5.402,  -1.418])

sp3_dimer_hamiltonian(element_a, element_b, bond_length_angstrom) builds an 8x8 Hamiltonian -- one s and three p orbitals per atom, two atoms -- for a real Si-Si bond at Si's real bulk bond length (2.35 Angstrom). The eigenvalues are the dimer's molecular-orbital energies in eV, built entirely from Harrison's tabulated atomic term values and universal eta bond-scaling coefficients, no per-material fitting anywhere.

Step 2. The full periodic crystal, at one k-point

H_k = de.zincblende_hamiltonian((0, 0, 0), 'Si', 'Si', lattice_constant_angstrom=5.43)
eigvals = np.linalg.eigvalsh(H_k)
float(eigvals[4] - eigvals[3])

3.665861271362627

zincblende_hamiltonian(k, cation, anion, lattice_constant_angstrom) sums the same bond physics with Bloch phases over the 4 nearest-neighbor bonds of a real zincblende/ diamond crystal -- k=(0,0,0) is the Gamma point, 5.43 Angstrom is silicon's real lattice constant. With 4 valence electrons per Si atom (8 total, 4 filled bands), the gap between the 4th and 5th eigenvalue at Gamma comes out 3.67 eV -- real silicon's actual band structure is more subtle than this single number (see Accuracy below).

---

Details

Accuracy: this is a universal model -- one parameter table for every material, no per-material fitting, no d-orbitals. That buys zero setup cost per new material at the price of accuracy: gaps typically come out ~2-3x off from experiment, and for indirect-gap materials (silicon included) this simple Gamma-only reading can altogether miss where the true conduction-band minimum sits -- see vhd_tb's own Step 2 for the real Si example (Gamma-only gives the wrong picture; scanning along a k-path finds the true, off-Gamma minimum). Validation numbers against real experimental gaps (GaAs, Si, Ge) are tracked in Dense-Evolution-Discovery.

Source: Walter A. Harrison, Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond (W. H. Freeman, 1980; Dover, 1989, ISBN 0-486-66021-4). Atomic term values (HARRISON_ELEMENTS) and universal eta coefficients (HARRISON_ETA) are transcribed from that book's Solid State Table, cross-checked against jarvist/HarrisonSolidStateTable.jl, an independent Julia implementation of the same table.

The bond-scaling law: off-diagonal hopping matrix elements follow V = eta * hbar^2/(m_e * d^2) (hbar^2/m_e = 7.62 eV*Angstrom^2), with eta_ss_sigma = -1.40, eta_sp_sigma = +1.84, eta_pp_sigma = +3.24, eta_pp_pi = -0.81 -- the same four numbers for every element pair, only the bond length d and each element's own atomic term values change between materials. Full sp3 Slater-Koster (1954) matrix elements for a bond of direction cosines (l, m, n) are in the module source.

harrison_tb

Harrison empirical tight-binding parameters -- builds an sp3 tight-binding Hamiltonian for a cluster of real atoms directly from published atomic term values and the universal bond-scaling law, with no SCF/DFT and no external quantum-chemistry dependency (PySCF, OpenFermion).

Source: Walter A. Harrison, "Electronic Structure and the Properties of Solids" (Dover reprint). Atomic term values (ELEMENTS) and the universal eta coefficients (ETA) below are transcribed from that book's Solid State Table, cross-checked against the numbers in jarvist/HarrisonSolidStateTable.jl (github.com/jarvist/HarrisonSolidStateTable.jl), a Julia implementation of the same table. Only elements with both an s and a p term value in that table are included here (the "simple atom" sp3 entries Harrison uses for tetrahedral semiconductors); d-block entries are left out since Harrison's own table flags them "not well checked".

Sign convention: ELEMENTS stores true orbital energies (negative, eV -- bound states below vacuum), i.e. the negative of the magnitudes printed in Harrison's table.

Two-center bond integrals follow Harrison's universal scaling law: V_ll'm = eta_ll'm * hbar^2 / (m_e * d^2) with d the bond length in Angstrom and hbar^2/m_e = 7.62 eV*Angstrom^2 (the standard constant quoted alongside this law). eta_ssσ/spσ/ppσ/ppπ are dimensionless and materials-independent -- the same four numbers apply to every element pair.

Slater-Koster sp3 matrix elements (sp3_bond_block) follow the standard 1954 Slater-Koster table for an (s, px, py, pz) basis.

zincblende_hamiltonian builds the periodic Bloch Hamiltonian for a two-atom zinc-blende basis (nearest-neighbor sp3, 4 bonds per atom) -- validated against real GaAs (a=5.6533 Angstrom): computed direct gap at Gamma is 2.91 eV vs. the experimental 1.42 eV, roughly 2x too large. This is a known, documented limitation of Harrison's universal (materials-independent) parameter set on polar/ionic compound semiconductors -- not a bug here -- since it uses no per-material fitting and omits d-orbitals. Useful as a fast, dependency-free qualitative estimate; not a substitute for this project's DFT-derived GaAs parameters where quantitative accuracy matters.

hopping_integral

hopping_integral(eta, d_angstrom)

V_ll'm = eta * hbar^2/(m_e d^2) [eV], for bond length d in Angstrom.

Source code in dense_evolution/solvers/harrison_tb.py

def hopping_integral(eta, d_angstrom):
    """V_ll'm = eta * hbar^2/(m_e d^2) [eV], for bond length d in Angstrom."""
    if d_angstrom <= 0:
        raise ValueError(f"bond length must be positive, got {d_angstrom}")
    return eta * HBAR2_OVER_M_EV_ANG2 / d_angstrom ** 2

sp3_bond_block

sp3_bond_block(l, m, n, d_angstrom, eta=ETA)

4x4 Slater-Koster hopping block for a bond from atom A to atom B along direction cosines (l, m, n) (unit vector, l^2+m^2+n^2 = 1) and bond length d_angstrom.

Basis order: s, px, py, pz. Standard Slater-Koster (1954) table.

Source code in dense_evolution/solvers/harrison_tb.py

def sp3_bond_block(l, m, n, d_angstrom, eta=ETA):
    """
    4x4 Slater-Koster hopping block <A, {s,px,py,pz}| H |B, {s,px,py,pz}>
    for a bond from atom A to atom B along direction cosines (l, m, n)
    (unit vector, l^2+m^2+n^2 = 1) and bond length d_angstrom.

    Basis order: s, px, py, pz. Standard Slater-Koster (1954) table.
    """
    norm = l * l + m * m + n * n
    if not np.isclose(norm, 1.0, atol=1e-6):
        raise ValueError(f"(l, m, n) must be a unit vector, got norm={norm}")

    Vssσ = hopping_integral(eta['ss_sigma'], d_angstrom)
    Vspσ = hopping_integral(eta['sp_sigma'], d_angstrom)
    Vppσ = hopping_integral(eta['pp_sigma'], d_angstrom)
    Vppπ = hopping_integral(eta['pp_pi'], d_angstrom)

    block = np.zeros((4, 4), dtype=np.complex128)
    block[0, 0] = Vssσ
    block[0, 1], block[0, 2], block[0, 3] = l * Vspσ, m * Vspσ, n * Vspσ
    block[1, 0], block[2, 0], block[3, 0] = -l * Vspσ, -m * Vspσ, -n * Vspσ

    block[1, 1] = l * l * Vppσ + (1 - l * l) * Vppπ
    block[2, 2] = m * m * Vppσ + (1 - m * m) * Vppπ
    block[3, 3] = n * n * Vppσ + (1 - n * n) * Vppπ

    block[1, 2] = block[2, 1] = l * m * (Vppσ - Vppπ)
    block[2, 3] = block[3, 2] = m * n * (Vppσ - Vppπ)
    block[1, 3] = block[3, 1] = l * n * (Vppσ - Vppπ)
    return block

sp3_dimer_hamiltonian

sp3_dimer_hamiltonian(
    element_a,
    element_b,
    bond_length_angstrom,
    direction=(0.0, 0.0, 1.0),
    eta=ETA,
)

8x8 sp3 tight-binding Hamiltonian for a 2-atom A-B cluster (one bond), basis order [A:s,px,py,pz, B:s,px,py,pz]. On-site blocks are each atom's diagonal (eps_s, eps_p, eps_p, eps_p); the A-B off-diagonal block is sp3_bond_block along direction (unit vector, default: bond along z), Hermitian-conjugated into the B-A block.

This is a minimal, directly checkable sanity case (bonding/antibonding sp3 splitting), not a periodic-solid band structure.

Source code in dense_evolution/solvers/harrison_tb.py

def sp3_dimer_hamiltonian(element_a, element_b, bond_length_angstrom,
                           direction=(0.0, 0.0, 1.0), eta=ETA):
    """
    8x8 sp3 tight-binding Hamiltonian for a 2-atom A-B cluster (one bond),
    basis order [A:s,px,py,pz, B:s,px,py,pz]. On-site blocks are each
    atom's diagonal (eps_s, eps_p, eps_p, eps_p); the A-B off-diagonal
    block is sp3_bond_block along `direction` (unit vector, default: bond
    along z), Hermitian-conjugated into the B-A block.

    This is a minimal, directly checkable sanity case (bonding/antibonding
    sp3 splitting), not a periodic-solid band structure.
    """
    for name in (element_a, element_b):
        if name not in ELEMENTS:
            raise ValueError(f"no Harrison sp term values for element {name!r}; "
                              f"available: {sorted(ELEMENTS)}")
    l, m, n = direction
    a, b = ELEMENTS[element_a], ELEMENTS[element_b]

    H = np.zeros((8, 8), dtype=np.complex128)
    H[0, 0] = a['eps_s']
    H[1, 1] = H[2, 2] = H[3, 3] = a['eps_p']
    H[4, 4] = b['eps_s']
    H[5, 5] = H[6, 6] = H[7, 7] = b['eps_p']

    hop = sp3_bond_block(l, m, n, bond_length_angstrom, eta=eta)
    H[0:4, 4:8] = hop
    H[4:8, 0:4] = hop.conj().T
    return H

zincblende_hamiltonian

zincblende_hamiltonian(
    k, cation, anion, lattice_constant_angstrom, eta=ETA
)

8x8 Bloch Hamiltonian for a zinc-blende crystal's two-atom basis (cation at (0,0,0), anion at (1/4,1/4,1/4) of the conventional cubic cell), sp3 nearest-neighbor tight-binding, at crystal momentum k (Cartesian, 1/Angstrom -- e.g. Gamma=(0,0,0)).

Basis order [cation:s,px,py,pz, anion:s,px,py,pz]. The four cation->anion nearest-neighbor bonds are the standard zinc-blende tetrahedral set (lattice_constant/4)*(1,1,1), (1,-1,-1), (-1,1,-1), (-1,-1,1); each contributes sp3_bond_block(...) weighted by its Bloch phase exp(i k . d), summed into the off-diagonal block.

Source code in dense_evolution/solvers/harrison_tb.py

def zincblende_hamiltonian(k, cation, anion, lattice_constant_angstrom, eta=ETA):
    """
    8x8 Bloch Hamiltonian for a zinc-blende crystal's two-atom basis
    (cation at (0,0,0), anion at (1/4,1/4,1/4) of the conventional cubic
    cell), sp3 nearest-neighbor tight-binding, at crystal momentum k
    (Cartesian, 1/Angstrom -- e.g. Gamma=(0,0,0)).

    Basis order [cation:s,px,py,pz, anion:s,px,py,pz]. The four
    cation->anion nearest-neighbor bonds are the standard zinc-blende
    tetrahedral set (lattice_constant/4)*(1,1,1), (1,-1,-1), (-1,1,-1),
    (-1,-1,1); each contributes sp3_bond_block(...) weighted by its
    Bloch phase exp(i k . d), summed into the off-diagonal block.
    """
    for name in (cation, anion):
        if name not in ELEMENTS:
            raise ValueError(f"no Harrison sp term values for element {name!r}; "
                              f"available: {sorted(ELEMENTS)}")
    k = np.asarray(k, dtype=float)
    d_vectors = (lattice_constant_angstrom / 4) * np.array([
        [1, 1, 1], [1, -1, -1], [-1, 1, -1], [-1, -1, 1],
    ], dtype=float)
    bond_length = np.linalg.norm(d_vectors[0])

    T = np.zeros((4, 4), dtype=np.complex128)
    for d in d_vectors:
        l, m, n = d / bond_length
        phase = np.exp(1j * np.dot(k, d))
        T += phase * sp3_bond_block(l, m, n, bond_length, eta=eta)

    c, a = ELEMENTS[cation], ELEMENTS[anion]
    H = np.zeros((8, 8), dtype=np.complex128)
    H[0, 0] = c['eps_s']
    H[1, 1] = H[2, 2] = H[3, 3] = c['eps_p']
    H[4, 4] = a['eps_s']
    H[5, 5] = H[6, 6] = H[7, 7] = a['eps_p']
    H[0:4, 4:8] = T
    H[4:8, 0:4] = T.conj().T
    return H

---

See also: vhd_tb for material-specific fitted parameters when the universal table's ~2-3x gap error isn't good enough.