wai-quantum 0.4.0

//! Molecules by DMRG — `wai.quantum.chem`.
//!
//! Quantum chemistry is where a quantum computer's value is most often
//! claimed, and a claim about a molecule's energy is checkable against a
//! classical method that reaches the same answer. This module takes a
//! molecule's integrals, writes its Hamiltonian as an exact matrix-product
//! operator, and finds the ground state by DMRG, held to full configuration
//! interaction.
//!
//! - **Integrals.** [`read_fcidump`] reads the FCIDUMP text format: one- and
//!   two-electron integrals in an orthonormal orbital basis, the electron
//!   count and spin, and the core energy. Any one of the eight equal orderings
//!   of a two-electron integral fills all eight.
//! - **The Hamiltonian as an MPO.** [`molecular_mpo`] places alpha and beta
//!   of each orbital on neighbouring sites, maps the fermions to them by the
//!   Jordan–Wigner transformation, and compiles every term with
//!   [`OpSum`]. That is an exact construction by minimum vertex covers, with
//!   no numerical compression. A penalty on the wrong electron count or spin
//!   is added, which costs no bond: its terms already appear in the Coulomb
//!   ones.
//! - **DMRG.** [`ground_state`] starts from the Hartree–Fock determinant with
//!   every bond's basis completed by null vectors, and keeps the bases
//!   complete through every split.
//! - **The referee.** [`fci`] is full configuration interaction: Lanczos with
//!   `σ = Hc` formed directly from the integrals over alpha and beta strings.
//!
//! # Checked
//!
//! - **Full CI.** It reproduces an independent program's full-CI energies to
//!   `10⁻⁹` for H₂, LiH, H₂O and a six-atom hydrogen chain (minimal basis),
//!   which also checks the reader.
//! - **The MPO is the Hamiltonian.** For H₂ its dense matrix has the full-CI
//!   energy as its lowest eigenvalue, with the penalty and without it (in the
//!   two-electron, zero-spin block). For H₂, LiH and H₂O its expectation in
//!   the Hartree–Fock determinant is the Hartree–Fock energy, which checks
//!   the fermion signs.
//! - **DMRG against full CI** (`chem_dmrg` example):
//!
//!   | molecule | spin orbitals | MPO bond | DMRG bond | DMRG − full CI |
//!   |----------|---------------|----------|-----------|----------------|
//!   | LiH      | 12            | 74       | 16        | 3·10⁻¹⁴         |
//!   | H₂O      | 14            | 138      | 64        | 3·10⁻¹³         |
//!   | N₂       | 20            | 278      | 64        | 3.6·10⁻⁴        |
//!
//!   LiH and H₂O reach full CI in one sweep. N₂ at bond 64 is limited by the
//!   bond: it converges over sweeps to `3.6·10⁻⁴` above full CI, inside
//!   chemical accuracy, with discarded weight `2.6·10⁻⁵`.
//!
//! # Measured, and fixed
//!
//! DMRG from the Hartree–Fock determinant with a positive weight cutoff did
//! not move at all. A product state's one-dimensional bases confine each
//! two-site update to excitations within its pair, and by Brillouin's theorem
//! those do not lower the energy. Completing the initial bases freed the
//! first sweep: LiH then came within `2·10⁻⁴` of full CI, and H₂O within
//! `4·10⁻²`. After that the cutoff shrank the bases back to the state's own
//! rank, which stalls the sweep the same way. Keeping the bases complete
//! through every split (cutoff zero) brought both to full CI. From a random
//! state, the penalty traps DMRG in a mixture of electron counts; the
//! Hartree–Fock start avoids that.
//!
//! # Honest boundaries
//!
//! - **Spin orbitals, no symmetry blocks.** Each site is one spin orbital,
//!   and conserved particle number and spin are enforced by a penalty, not
//!   used to block the matrices. The cost goes as the MPO bond, which grows
//!   as the square of the orbital count, times the cube of the DMRG bond.
//!   Spatial-orbital sites with conserved quantities would be far cheaper.
//! - **The orbital order is the file's.** Orbitals are not reordered to
//!   shorten the reach of entanglement.
//! - **Real integrals** in the FCIDUMP format only.
//! - **The MPO is exact, not minimal for every structure.** The vertex-cover
//!   construction is optimal for terms with independent coefficients. It
//!   does not find the lower rank of coefficient matrices with structure,
//!   which a numerical compression would.

use crate::quantum_dmrg::{DmrgConfig, DmrgResult, Mpo, Mps, Op, OpSum, complete_bases, dmrg_mpo, lanczos, right_canonicalise};

// ---------------------------------------------------------------------------
// Integrals
// ---------------------------------------------------------------------------

/// A molecular Hamiltonian in an orthonormal orbital basis: one- and
/// two-electron integrals, the electron count and spin, and the core energy.
#[derive(Clone, Debug, PartialEq)]
pub struct Fcidump {
    pub norb: usize,
    pub nelec: usize,
    /// Twice the spin projection, `N_α − N_β`.
    pub ms2: i64,
    /// `h[p][q]`, row-major.
    pub h1: Vec<f64>,
    /// `(pq|rs)` in chemists' order, `[((p·n + q)·n + r)·n + s]`.
    pub h2: Vec<f64>,
    /// Nuclear repulsion plus any frozen-core energy.
    pub ecore: f64,
}

impl Fcidump {
    pub fn h1(&self, p: usize, q: usize) -> f64 {
        self.h1[p * self.norb + q]
    }

    pub fn h2(&self, p: usize, q: usize, r: usize, s: usize) -> f64 {
        let n = self.norb;
        self.h2[((p * n + q) * n + r) * n + s]
    }

    /// Electrons of each spin, `(N_α, N_β)`.
    pub fn electrons(&self) -> (usize, usize) {
        let up = (self.nelec as i64 + self.ms2) / 2;
        (up as usize, self.nelec - up as usize)
    }
}

fn header_int(header: &str, key: &str) -> Option<i64> {
    let at = header.find(key)? + key.len();
    let rest = header[at..].trim_start().strip_prefix('=')?.trim_start();
    let end = rest.find(|c: char| !(c.is_ascii_digit() || c == '-' || c == '+')).unwrap_or(rest.len());
    rest[..end].parse().ok()
}

/// Read the FCIDUMP text format: a namelist header (`NORB`, `NELEC`, `MS2`)
/// closed by `&END` or `/`, then lines `value i j k l` with 1-based indices:
/// two-electron integrals `(ij|kl)` (any one of the eight equal orderings),
/// one-electron integrals with `k = l = 0`, and the core energy with all four
/// zero. Orbital energies (`i > 0`, `j = k = l = 0`) are ignored.
pub fn read_fcidump(text: &str) -> Result<Fcidump, String> {
    let upper = text.to_ascii_uppercase();
    let end = upper.find("&END").or_else(|| upper.find("\n/")).or_else(|| upper.find("/\n")).ok_or("no end to the FCIDUMP header (&END or /)")?;
    let header = &upper[..end];
    let norb = header_int(header, "NORB").ok_or("NORB missing")?;
    let nelec = header_int(header, "NELEC").ok_or("NELEC missing")?;
    let ms2 = header_int(header, "MS2").unwrap_or(0);
    if !(1..=64).contains(&norb) || nelec < 0 || nelec > 2 * norb || (nelec + ms2) % 2 != 0 || ms2.abs() > nelec {
        return Err(format!("inconsistent header: NORB={norb} NELEC={nelec} MS2={ms2}"));
    }
    let n = norb as usize;
    let mut fd = Fcidump { norb: n, nelec: nelec as usize, ms2, h1: vec![0.0; n * n], h2: vec![0.0; n * n * n * n], ecore: 0.0 };
    let body_start = text[end..].find('\n').map_or(text.len(), |k| end + k + 1);
    for (lineno, line) in text[body_start..].lines().enumerate() {
        let toks: Vec<&str> = line.split_whitespace().collect();
        if toks.is_empty() {
            continue;
        }
        if toks.len() != 5 {
            return Err(format!("body line {}: expected `value i j k l`", lineno + 1));
        }
        let v: f64 = toks[0].replace(['D', 'd'], "e").parse().map_err(|_| format!("body line {}: bad value {}", lineno + 1, toks[0]))?;
        let mut idx = [0usize; 4];
        for (slot, t) in idx.iter_mut().zip(&toks[1..]) {
            let x: usize = t.parse().map_err(|_| format!("body line {}: bad index {t}", lineno + 1))?;
            if x > n {
                return Err(format!("body line {}: index {x} past NORB", lineno + 1));
            }
            *slot = x;
        }
        let [i, j, k, l] = idx;
        match (i, j, k, l) {
            (0, 0, 0, 0) => fd.ecore = v,
            (_, 0, 0, 0) => {}
            (i, j, 0, 0) if i > 0 && j > 0 => {
                fd.h1[(i - 1) * n + j - 1] = v;
                fd.h1[(j - 1) * n + i - 1] = v;
            }
            (i, j, k, l) if i > 0 && j > 0 && k > 0 && l > 0 => {
                let (p, q, r, s) = (i - 1, j - 1, k - 1, l - 1);
                for (a, b, c, d) in [(p, q, r, s), (q, p, r, s), (p, q, s, r), (q, p, s, r), (r, s, p, q), (s, r, p, q), (r, s, q, p), (s, r, q, p)] {
                    fd.h2[((a * n + b) * n + c) * n + d] = v;
                }
            }
            _ => return Err(format!("body line {}: indices {i} {j} {k} {l} are not a known kind", lineno + 1)),
        }
    }
    Ok(fd)
}

// ---------------------------------------------------------------------------
// The molecular MPO
// ---------------------------------------------------------------------------

/// The site of spin orbital `(p, σ)`: alpha and beta of each spatial orbital
/// side by side, `2p + σ`.
pub fn site(p: usize, beta: bool) -> usize {
    2 * p + beta as usize
}

/// `|1⟩⟨0|`: creates an electron (index 1 is occupied).
fn cre() -> Op {
    Op([0.0, 0.0, 1.0, 0.0])
}

/// `|0⟩⟨1|`: annihilates one.
fn ann() -> Op {
    Op([0.0, 1.0, 0.0, 0.0])
}

/// `(−1)^n`, the Jordan–Wigner string.
fn parity() -> Op {
    Op([1.0, 0.0, 0.0, -1.0])
}

fn number() -> Op {
    Op([0.0, 0.0, 0.0, 1.0])
}

/// The one-site factors of a fermion operator on site `j` after the
/// Jordan–Wigner transformation: the parity string on every site before it.
fn jordan_wigner(j: usize, create: bool, out: &mut Vec<(usize, Op)>) {
    for k in 0..j {
        out.push((k, parity()));
    }
    out.push((j, if create { cre() } else { ann() }));
}

/// Settings for the molecular Hamiltonian and its DMRG.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct ChemConfig {
    /// Integrals smaller than this in magnitude are left out.
    pub drop: f64,
    /// Weight of `(N̂ − N)² + (2Ŝ_z − 2S_z)²` added to the Hamiltonian. It is
    /// zero on the right electron count and spin and positive elsewhere, so
    /// rounding cannot carry the ground state out of its sector.
    pub penalty: f64,
    pub max_bond: usize,
    /// Full sweeps (left to right and back).
    pub sweeps: usize,
    /// Lanczos residual at which each two-site optimisation stops.
    pub tol: f64,
    /// Threads for the effective Hamiltonian; the result does not depend on
    /// it.
    pub threads: usize,
}

impl Default for ChemConfig {
    fn default() -> ChemConfig {
        ChemConfig { drop: 1e-12, penalty: 1.0, max_bond: 64, sweeps: 4, tol: 1e-8, threads: DmrgConfig::default().threads }
    }
}

/// The second-quantised Hamiltonian
/// `E_core + Σ h_pq a†_pσ a_qσ + ½ Σ (pq|rs) a†_pσ a†_rτ a_sτ a_qσ`, plus the
/// sector penalty, as an operator sum over `2·norb` spin-orbital sites.
pub fn molecular_opsum(fd: &Fcidump, cfg: &ChemConfig) -> OpSum {
    let n = fd.norb;
    let mut sum = OpSum::new(2 * n);
    sum.add(fd.ecore, &[]);
    let mut f = Vec::new();
    for p in 0..n {
        for q in 0..n {
            let h = fd.h1(p, q);
            if h.abs() < cfg.drop {
                continue;
            }
            for beta in [false, true] {
                f.clear();
                jordan_wigner(site(p, beta), true, &mut f);
                jordan_wigner(site(q, beta), false, &mut f);
                sum.add(h, &f);
            }
        }
    }
    for p in 0..n {
        for q in 0..n {
            for r in 0..n {
                for s in 0..n {
                    let v = fd.h2(p, q, r, s);
                    if v.abs() < cfg.drop {
                        continue;
                    }
                    for sigma in [false, true] {
                        for tau in [false, true] {
                            let (ps, qs, rt, st) = (site(p, sigma), site(q, sigma), site(r, tau), site(s, tau));
                            if ps == rt || qs == st {
                                continue;
                            }
                            f.clear();
                            jordan_wigner(ps, true, &mut f);
                            jordan_wigner(rt, true, &mut f);
                            jordan_wigner(st, false, &mut f);
                            jordan_wigner(qs, false, &mut f);
                            sum.add(0.5 * v, &f);
                        }
                    }
                }
            }
        }
    }
    if cfg.penalty != 0.0 {
        // (N̂ − N)² + (2Ŝ_z − m)² = Σ_ij (1 + s_i s_j) n_i n_j
        //   − 2 Σ_j (N + m s_j) n_j + N² + m², with s_j = ±1 for α, β.
        let (ne, m) = (fd.nelec as f64, fd.ms2 as f64);
        let sign = |j: usize| if j.is_multiple_of(2) { 1.0 } else { -1.0 };
        for i in 0..2 * n {
            for j in 0..2 * n {
                let c = 1.0 + sign(i) * sign(j);
                if c != 0.0 {
                    sum.add(cfg.penalty * c, &[(i, number()), (j, number())]);
                }
            }
            sum.add(-2.0 * cfg.penalty * (ne + m * sign(i)), &[(i, number())]);
        }
        sum.add(cfg.penalty * (ne * ne + m * m), &[]);
    }
    sum
}

/// The molecular Hamiltonian as an MPO.
pub fn molecular_mpo(fd: &Fcidump, cfg: &ChemConfig) -> Mpo {
    molecular_opsum(fd, cfg).mpo()
}

/// The Hartree–Fock determinant in the given orbital order: the lowest
/// `N_α` alpha and `N_β` beta orbitals occupied.
pub fn hartree_fock_state(fd: &Fcidump) -> Mps {
    let (na, nb) = fd.electrons();
    let n = 2 * fd.norb;
    let sites = (0..n)
        .map(|j| {
            let (p, beta) = (j / 2, j % 2 == 1);
            let occupied = if beta { p < nb } else { p < na };
            if occupied { vec![0.0, 1.0] } else { vec![1.0, 0.0] }
        })
        .collect();
    Mps { dims: vec![1; n + 1], sites }
}

/// The Hartree–Fock energy of [`hartree_fock_state`]: `E_core + Σ_occ h_ii
/// + ½ Σ_occ,occ [(ii|jj) − δ_σσ' (ij|ji)]`.
pub fn hartree_fock_energy(fd: &Fcidump) -> f64 {
    let (na, nb) = fd.electrons();
    let occ: Vec<(usize, bool)> = (0..na).map(|p| (p, false)).chain((0..nb).map(|p| (p, true))).collect();
    let mut e = fd.ecore;
    for &(i, _) in &occ {
        e += fd.h1(i, i);
    }
    for &(i, si) in &occ {
        for &(j, sj) in &occ {
            e += 0.5 * fd.h2(i, i, j, j);
            if si == sj {
                e -= 0.5 * fd.h2(i, j, j, i);
            }
        }
    }
    e
}

/// The molecule's ground state by two-site DMRG on its MPO.
///
/// DMRG starts from the Hartree–Fock determinant, which has the right
/// electron count and spin, with every bond's basis completed by null vectors
/// up to `max_bond`, and keeps them complete at every split. A product
/// state's one-dimensional bases would otherwise hold the two-site updates to
/// excitations within the pair being updated; by Brillouin's theorem those do
/// not lower the energy, and the sweep would stay at Hartree–Fock. Bases cut
/// back to the state's own rank stall it the same way later on.
pub fn ground_state(fd: &Fcidump, cfg: &ChemConfig) -> DmrgResult {
    let mpo = molecular_mpo(fd, cfg);
    let mut start = hartree_fock_state(fd);
    right_canonicalise(&mut start);
    complete_bases(&mut start, cfg.max_bond);
    // Cutoff zero: every split keeps its full basis up to the bond, null
    // vectors included, so the bases never shrink to what the state uses
    // and the sweep cannot stall on them.
    let dmrg = DmrgConfig { max_bond: cfg.max_bond, cutoff: 0.0, sweeps: cfg.sweeps, tol: cfg.tol, threads: cfg.threads, ..DmrgConfig::default() };
    dmrg_mpo(&mpo, start, &dmrg)
}

// ---------------------------------------------------------------------------
// Full configuration interaction, the referee
// ---------------------------------------------------------------------------

/// All `norb`-bit strings with `k` bits set, ascending.
fn strings(norb: usize, k: usize) -> Vec<u64> {
    let mut out = Vec::new();
    if k > norb {
        return out;
    }
    if k == 0 {
        return vec![0];
    }
    let mut x: u64 = (1u64 << k) - 1;
    let limit: u64 = 1u64 << norb;
    while x < limit {
        out.push(x);
        // Gosper's hack: the next integer with the same number of bits.
        let c = x.isolate_lowest_one();
        let r = x + c;
        x = (((r ^ x) >> 2) / c) | r;
    }
    out
}

/// For each string, every `E_kl |I⟩ = sign |J⟩`: `(k·n + l, J, sign)`.
fn replacements(norb: usize, list: &[u64]) -> Vec<Vec<(usize, usize, f64)>> {
    list.iter()
        .map(|&s| {
            let mut out = Vec::new();
            for l in 0..norb {
                if s >> l & 1 == 0 {
                    continue;
                }
                let below_l = (s & ((1u64 << l) - 1)).count_ones();
                let without = s & !(1u64 << l);
                for k in 0..norb {
                    if without >> k & 1 == 1 {
                        continue;
                    }
                    let below_k = (without & ((1u64 << k) - 1)).count_ones();
                    let t = without | (1u64 << k);
                    let j = list.binary_search(&t).expect("the string list is closed under replacements");
                    let sign = if (below_l + below_k).is_multiple_of(2) { 1.0 } else { -1.0 };
                    out.push((k * norb + l, j, sign));
                }
            }
            out
        })
        .collect()
}

/// The result of [`fci`].
#[derive(Clone, Debug, PartialEq)]
pub struct FciResult {
    pub energy: f64,
    /// Determinants in the `(N_α, N_β)` sector.
    pub dim: usize,
}

/// The exact ground energy in the `(N_α, N_β)` sector by full configuration
/// interaction: Lanczos with `σ = Hc` formed directly from the integrals
/// over alpha and beta strings (at most two million determinants).
pub fn fci(fd: &Fcidump) -> Option<FciResult> {
    let n = fd.norb;
    if n > 32 {
        return None;
    }
    let (na, nb) = fd.electrons();
    let (sa, sb) = (strings(n, na), strings(n, nb));
    let (da, db) = (sa.len(), sb.len());
    let dim = da * db;
    if dim == 0 || dim > 2_000_000 {
        return None;
    }
    let (ra, rb) = (replacements(n, &sa), replacements(n, &sb));
    let nn = n * n;
    // h'_kl = h_kl − ½ Σ_m (km|ml).
    let mut hp = vec![0.0; nn];
    for k in 0..n {
        for l in 0..n {
            let mut s = fd.h1(k, l);
            for m in 0..n {
                s -= 0.5 * fd.h2(k, m, m, l);
            }
            hp[k * n + l] = s;
        }
    }
    let apply = |c: &[f64], sigma: &mut [f64]| {
        // D[kl][I] = Σ_J ⟨I|E_kl|J⟩ c_J, scattered from each J.
        let mut d = vec![0.0; nn * dim];
        for ja in 0..da {
            for jb in 0..db {
                let cj = c[ja * db + jb];
                if cj == 0.0 {
                    continue;
                }
                for &(kl, ia, sg) in &ra[ja] {
                    d[kl * dim + ia * db + jb] += sg * cj;
                }
                for &(kl, ib, sg) in &rb[jb] {
                    d[kl * dim + ja * db + ib] += sg * cj;
                }
            }
        }
        // G[kl][I] = ½ Σ_mn (kl|mn) D[mn][I].
        let mut g = vec![0.0; nn * dim];
        for kl in 0..nn {
            let dst = &mut g[kl * dim..(kl + 1) * dim];
            for mn in 0..nn {
                let v = 0.5 * fd.h2[kl * nn + mn];
                if v == 0.0 {
                    continue;
                }
                for (x, &y) in dst.iter_mut().zip(&d[mn * dim..(mn + 1) * dim]) {
                    *x += v * y;
                }
            }
        }
        // σ = E_core c + Σ_kl h'_kl D_kl + Σ_kl E_kl G_kl.
        for (i, s) in sigma.iter_mut().enumerate() {
            let mut acc = fd.ecore * c[i];
            for kl in 0..nn {
                acc += hp[kl] * d[kl * dim + i];
            }
            *s = acc;
        }
        for ia in 0..da {
            for ib in 0..db {
                let i = ia * db + ib;
                for &(kl, ta, sg) in &ra[ia] {
                    sigma[ta * db + ib] += sg * g[kl * dim + i];
                }
                for &(kl, tb, sg) in &rb[ib] {
                    sigma[ia * db + tb] += sg * g[kl * dim + i];
                }
            }
        }
    };
    // Start from the Hartree–Fock determinant (string 0 of each spin) with
    // a little of every other determinant, so no symmetry hides the ground
    // state.
    let start: Vec<f64> = (0..dim as u64).map(|x| if x == 0 { 1.0 } else { 1e-3 * (1.0 + (x.wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0) }).collect();
    let (energy, _) = lanczos(&apply, &start, 60, 40, 1e-9, false);
    Some(FciResult { energy, dim })
}

#[cfg(test)]
mod tests {
    use super::*;

    const H2: &str = include_str!("../tests/data/h2.fcidump");
    const LIH: &str = include_str!("../tests/data/lih.fcidump");
    const H2O: &str = include_str!("../tests/data/h2o.fcidump");
    const H6: &str = include_str!("../tests/data/h6.fcidump");

    /// Full-CI energies from an independent program for the same integrals.
    const REFERENCE: [(&str, f64); 4] = [("h2", -1.137283834488502), ("lih", -7.882401932290214), ("h2o", -75.01257824109194), ("h6", -3.2360662798923476)];

    fn read(name: &str) -> Fcidump {
        let text = match name {
            "h2" => H2,
            "lih" => LIH,
            "h2o" => H2O,
            "h6" => H6,
            _ => unreachable!(),
        };
        read_fcidump(text).unwrap()
    }

    #[test]
    fn full_ci_reproduces_the_reference_energies() {
        for (name, want) in REFERENCE {
            let r = fci(&read(name)).unwrap();
            assert!((r.energy - want).abs() < 1e-9, "{name}: {} vs {want}", r.energy);
        }
    }

    #[test]
    fn the_mpo_is_the_hamiltonian() {
        // On H₂ (four spin orbitals) the MPO's dense matrix, restricted to
        // two electrons, has the full-CI ground energy as its lowest
        // eigenvalue, and the penalty lifts every other sector.
        let fd = read("h2");
        let mpo = molecular_mpo(&fd, &ChemConfig::default());
        let dense = mpo.to_dense().unwrap();
        let dim = 16;
        let (vals, _) = crate::quantum_dmrg::symmetric_eigen(&dense, dim);
        assert!((vals[0] - REFERENCE[0].1).abs() < 1e-10, "{}", vals[0]);
        // Without the penalty the two-electron block still holds it.
        let bare = molecular_mpo(&fd, &ChemConfig { penalty: 0.0, ..ChemConfig::default() }).to_dense().unwrap();
        let two: Vec<usize> = (0..dim).filter(|x: &usize| x.count_ones() == 2 && (x >> 3 & 1) + (x >> 1 & 1) == 1).collect();
        let k = two.len();
        let block: Vec<f64> = two.iter().flat_map(|&r| two.iter().map(move |&c| (r, c))).map(|(r, c)| bare[r * dim + c]).collect();
        let (bv, _) = crate::quantum_dmrg::symmetric_eigen(&block, k);
        assert!((bv[0] - REFERENCE[0].1).abs() < 1e-10, "{}", bv[0]);
    }

    #[test]
    fn hartree_fock_energy_is_the_determinants() {
        // The first diagonal element of the MPO at the HF determinant.
        for name in ["h2", "lih", "h2o"] {
            let fd = read(name);
            let mpo = molecular_mpo(&fd, &ChemConfig::default());
            let hf = hartree_fock_state(&fd);
            let e = crate::quantum_dmrg::expectation(&mpo, &hf);
            assert!((e - hartree_fock_energy(&fd)).abs() < 1e-10, "{name}: {e} vs {}", hartree_fock_energy(&fd));
        }
    }

    #[test]
    fn dmrg_reaches_full_ci() {
        // H₂'s four spin orbitals fit any bond of four: exact.
        let fd = read("h2");
        let r = ground_state(&fd, &ChemConfig { max_bond: 4, sweeps: 1, ..ChemConfig::default() });
        assert!((r.energy - REFERENCE[0].1).abs() < 1e-10, "h2: {}", r.energy);
        // LiH at bond 12 cannot hold its ground state exactly (bond 16 does):
        // above full CI, never below, by about the discarded weight.
        let fd = read("lih");
        let exact = fci(&fd).unwrap().energy;
        let r = ground_state(&fd, &ChemConfig { max_bond: 12, sweeps: 1, ..ChemConfig::default() });
        assert!(r.energy >= exact - 1e-10, "lih: below full CI: {} vs {exact}", r.energy);
        assert!(r.energy - exact < 1e-5, "lih: {} vs {exact}", r.energy);
        assert!(r.discarded > 0.0);
    }

    #[test]
    fn the_reader_refuses_what_it_cannot_read() {
        assert!(read_fcidump("NORB=2").is_err());
        assert!(read_fcidump(" &FCI NORB=2,NELEC=5,MS2=0,\n &END\n").is_err());
        assert!(read_fcidump(" &FCI NORB=2,NELEC=2,MS2=0,\n &END\n 1.0 3 1 1 1\n").is_err());
        assert!(read_fcidump(" &FCI NORB=2,NELEC=2,MS2=0,\n &END\n 1.0 1 1\n").is_err());
        let fd = read_fcidump(" &FCI NORB=2,NELEC=2,MS2=0,\n &END\n 0.5D0 1 2 1 2\n -1.25 1 1 0 0\n 0.7 0 0 0 0\n").unwrap();
        assert_eq!(fd.h2(1, 0, 2 - 1, 0), 0.5);
        assert_eq!(fd.h2(1, 0, 0, 1), 0.5);
        assert_eq!(fd.h1(0, 0), -1.25);
        assert_eq!(fd.ecore, 0.7);
    }
}