wai/wai-quantum-serve/Cargo.lock

@@ -589,7 +589,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.35"
+version = "0.3.36"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.lock

@@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.35"
+version = "0.3.36"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wai-quantum"
-version = "0.3.35"
+version = "0.3.36"
 edition = "2024"
 rust-version = "1.97"
 license = "Apache-2.0"
@@ -85,6 +85,13 @@ quantum_dmrg = []
 # Time evolution of spin chains by the time-dependent variational principle,
 # held to exact evolution and to the free-fermion solution of a quench.
 quantum_tdvp = ["quantum_dmrg"]
+# Molecules by DMRG: FCIDUMP integrals, the second-quantised Hamiltonian as an
+# exact MPO, and full configuration interaction as the referee.
+quantum_chem = ["quantum_dmrg"]
+# Neural quantum states: a restricted Boltzmann machine trained by variational
+# Monte Carlo with stochastic reconfiguration, held to exact diagonalisation
+# and DMRG on chains and square lattices.
+quantum_nqs = ["quantum_dmrg"]
 # Mixed states and the observables that distinguish quantum theory from a classical
 # account: density matrices, the partial trace, exact Pauli expectations, CHSH.
 quantum_info = ["quantum"]
@@ -124,7 +131,7 @@ quantum_phasor_meter = ["_signing"]
 
 # Everything.
 full = [
-    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "quantum_tdvp", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
+    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "quantum_tdvp", "quantum_chem", "quantum_nqs", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
     "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec",
     "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom",
     "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata",
@@ -135,6 +142,14 @@ full = [
 name = "qio_gset"
 required-features = ["quantum_qio"]
 
+[[example]]
+name = "tdvp_quench"
+required-features = ["quantum_tdvp"]
+
+[[example]]
+name = "chem_dmrg"
+required-features = ["quantum_chem"]
+
 [[example]]
 name = "dmrg_chain"
 required-features = ["quantum_dmrg"]

wai/wai-quantum/README.md

@@ -68,6 +68,8 @@ never creates information.
 | `quantum_uf` | **union-find decoding in near-linear time** over the same matching graph: event-driven cluster growth on integer weights (a radix heap orders the fill events), merges, and peeling. Measured against exact matching on the same shots: within 0.92–1.18× the failures, up to 77× faster at `d = 11` |
 | `quantum_dmrg` | **ground states of spin chains by DMRG**: matrix-product operators for nearest-neighbour chains (XXZ Heisenberg in a field, transverse-field Ising), two-site sweeps with Lanczos and density-matrix truncation, entanglement entropies, a symmetric eigensolver, and two exact referees: full diagonalisation up to 22 sites, and the free-fermion solution of the Ising chain at any length. The critical Ising chain is reproduced to about `10⁻¹¹` per site at 200 sites |
 | `quantum_tdvp` | **time evolution of spin chains by TDVP**: two-site sweeps forward and back, Lanczos exponentials that halve a step they cannot converge, bonds split by singular value decomposition, and every basis completed with null vectors so a product state's small bonds do not hold the evolution back. Measures site and bond observables, entanglement at every bond, energy and the Loschmidt echo. Held to two exact referees: Krylov evolution of all amplitudes up to 20 sites, and the free-fermion solution of an Ising quench at any length. On a 100-site critical quench at bond 64 it stays within `10⁻⁹` of the exact answer until the entanglement fills the bond, and past that its error follows the discarded weight it reports |
+| `quantum_chem` | **molecules by DMRG**: FCIDUMP integrals, the second-quantised Hamiltonian on spin orbitals compiled into an exact MPO (minimum vertex covers, no numerical compression), DMRG from Hartree–Fock with complete bases, and full configuration interaction as the referee. LiH and H₂O reach full CI to `10⁻¹³` |
+| `quantum_nqs` | **neural quantum states**: spin models on any graph, a restricted Boltzmann machine (with Marshall's sign on bipartite lattices), variational Monte Carlo by Metropolis chains or exact enumeration, and stochastic reconfiguration solved in whichever of parameter or sample space is smaller. Held to exact diagonalisation and to DMRG on the same operator sum |
 | `quantum_noisy` | **noisy circuits**: quantum trajectories on the state vector (depolarizing, dephasing, amplitude damping, readout), an exact density-matrix referee up to 12 qubits that reuses the same kernels, linear cross-entropy fidelity (sampled, exact and the circuit's own), and the digital error model. Measures where that model is off: below the true fidelity by a margin that grows with the error rate, not with depth |
 | `quantum_ldpc` | **qLDPC memory at circuit level**: bivariate-bicycle codes from their two polynomials (`k` by rank, logical bases), the depth-8 syndrome cycle with circuit-level noise, its detector error model, and BP-OSD decoding (min-sum with adaptive scaling, ordered statistics with a combination sweep). The error model has the same columns as the source's own decoding matrices, and on identical syndromes the decoder makes the source's decisions shot for shot |
 | `quantum_arch` | **architecture-level estimation**: yoked surface-code storage (row and grid parity checks over patches, hot and cold), magic-state cultivation read off the published simulation curves, CCZ factories, and a clock set by whichever is slowest of lattice surgery, the magic-state supply and the control loop. Assembles whole-machine estimates; its storage search matches the source's footprint tool layout for layout |
@@ -220,6 +222,35 @@ rises with the discarded weight, and the run says so. A chain small enough to
 fit the bond evolves exactly, to rounding, at any time step. Runs repeat bit
 for bit at any thread count, natively and as WebAssembly.
 
+## Molecules
+
+```sh
+wai-quantum chem --fcidump h2o.fcidump --bond 64 --sweeps 2   # against Hartree–Fock and full CI
+```
+
+Quantum chemistry is where quantum computers are most often said to pay off,
+and an energy claim for a molecule is checkable when a classical method
+reaches the same number. `quantum_chem` reads 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 computed in the same
+crate. LiH and water reach full CI to `10⁻¹³` in one sweep. N₂, on 20 spin
+orbitals at bond 64, comes within `4·10⁻⁴` of it, and reports the discarded
+weight that accounts for the difference.
+
+## Neural quantum states
+
+```sh
+wai-quantum nqs --lx 8 --exact --steps 200 --lr 0.1          # every configuration enumerated: no sampling noise
+wai-quantum nqs --model ising --lx 4 --ly 4 --h 3.044 --alpha 1 --steps 400
+```
+
+Neural-network wavefunctions are offered as evidence about systems that no
+exact method reaches. `quantum_nqs` trains one the standard way and holds it
+to exact answers. On a 4 × 4 lattice at the Ising model's critical field, it
+lands within its own error bar of exact diagonalisation. Past exact reach,
+the same Hamiltonian goes to DMRG through its exact MPO, so two independent
+methods meet on one number.
+
 ## Noisy circuits
 
 ```sh

wai/wai-quantum/examples/chem_dmrg.rs

@@ -0,0 +1,33 @@
+//! A molecule's ground state by DMRG from its FCIDUMP integrals, against
+//! Hartree–Fock and, where it fits, full configuration interaction.
+//!
+//! ```text
+//! cargo run --release --example chem_dmrg --features quantum_chem -- FILE BOND SWEEPS [threads]
+//! ```
+use wai_quantum::quantum_chem::{fci, ground_state, hartree_fock_energy, molecular_mpo, read_fcidump, ChemConfig};
+
+fn main() {
+    let a: Vec<String> = std::env::args().collect();
+    let fd = read_fcidump(&std::fs::read_to_string(&a[1]).expect("read the FCIDUMP")).expect("parse the FCIDUMP");
+    let bond: usize = a[2].parse().unwrap();
+    let sweeps: usize = a[3].parse().unwrap();
+    let mut cfg = ChemConfig { max_bond: bond, sweeps, ..ChemConfig::default() };
+    if let Some(t) = a.get(4) {
+        cfg.threads = t.parse().unwrap();
+    }
+    let mpo = molecular_mpo(&fd, &cfg);
+    println!("{}: {} orbitals, {} electrons, MPO bond {}", a[1], fd.norb, fd.nelec, mpo.max_bond());
+    println!("Hartree–Fock  {:.10}", hartree_fock_energy(&fd));
+    let t0 = std::time::Instant::now();
+    let r = ground_state(&fd, &cfg);
+    let secs = t0.elapsed().as_secs_f64();
+    for (i, e) in r.sweep_energies.iter().enumerate() {
+        println!("sweep {:>2}      {e:.10}", i + 1);
+    }
+    println!("discarded     {:.1e}   entropy, mid {:.4}   {secs:.1}s", r.discarded, r.entropies[r.entropies.len() / 2]);
+    let t0 = std::time::Instant::now();
+    match fci(&fd) {
+        Some(f) => println!("full CI       {:.10}   ({} determinants, {:.1}s)   DMRG − full CI {:.2e}", f.energy, f.dim, t0.elapsed().as_secs_f64(), r.energy - f.energy),
+        None => println!("full CI       (too large)"),
+    }
+}

wai/wai-quantum/examples/nqs_ground.rs

@@ -0,0 +1,57 @@
+//! A neural quantum state for a spin model's ground state, by variational
+//! Monte Carlo with stochastic reconfiguration, against exact
+//! diagonalisation (up to 22 spins) or DMRG on the same operator sum.
+//!
+//! ```text
+//! cargo run --release --example nqs_ground --features quantum_nqs -- heisenberg|ising LX LY ALPHA ITERS SAMPLES [h] [dmrg-bond] [threads]
+//! ```
+use wai_quantum::quantum_dmrg::{dmrg_mpo, DmrgConfig, Mps};
+use wai_quantum::quantum_nqs::*;
+
+fn main() {
+    let a: Vec<String> = std::env::args().collect();
+    let (lx, ly): (usize, usize) = (a[2].parse().unwrap(), a[3].parse().unwrap());
+    let alpha: usize = a[4].parse().unwrap();
+    let iters: usize = a[5].parse().unwrap();
+    let samples: usize = a[6].parse().unwrap();
+    let h: f64 = a.get(7).map_or(3.044, |x| x.parse().unwrap());
+    let bond: usize = a.get(8).map_or(128, |x| x.parse().unwrap());
+    let threads: usize = a.get(9).map_or(VmcConfig::default().threads, |x| x.parse().unwrap());
+    let n = lx * ly;
+    let edges = if ly == 1 { chain_edges(n, false) } else { square_edges(lx, ly, false) };
+    let heis = a[1] == "heisenberg";
+    let model = if heis { SpinModel::heisenberg(n, &edges, 1.0) } else { SpinModel::ising(n, &edges, 1.0, h) };
+    let sector = if heis { Some(0) } else { None };
+    let t0 = std::time::Instant::now();
+    let (reference, how) = match model.exact_ground_energy(sector) {
+        Some(e) if n <= 20 => (e, "exact diagonalisation".to_string()),
+        _ => {
+            let mpo = model.opsum().mpo();
+            let r = dmrg_mpo(&mpo, Mps::random(n, 8, 1), &DmrgConfig { max_bond: bond, sweeps: 6, threads, ..DmrgConfig::default() });
+            (r.energy, format!("DMRG, bond {bond}, discarded {:.1e}", r.discarded))
+        }
+    };
+    println!("{} {lx}×{ly}: reference {reference:.8} ({how}, {:.1}s)", a[1], t0.elapsed().as_secs_f64());
+    let mut rbm = Rbm::new(n, alpha, 0.01, 7);
+    if heis {
+        rbm = rbm.with_marshall_sign(checkerboard(lx, ly));
+    }
+    println!("RBM α={alpha}: {} parameters, {samples} samples a step", rbm.params());
+    let cfg = VmcConfig { samples, chains: 16, magnetisation: sector, lr: 0.02, shift: 1e-3, threads, ..VmcConfig::default() };
+    let mut vmc = Vmc::new(model, rbm, cfg);
+    let t0 = std::time::Instant::now();
+    let mut tail = Vec::new();
+    for it in 1..=iters {
+        let e = vmc.step();
+        if it > iters * 4 / 5 {
+            tail.push(e.energy);
+        }
+        if it % (iters / 10).max(1) == 0 {
+            println!("step {it:>5}  E {:.6} ± {:.6}  variance {:.2e}  acceptance {:.2}  relative error {:.2e}  {:.0}s", e.energy, e.error, e.variance, e.acceptance, (e.energy - reference) / reference.abs(), t0.elapsed().as_secs_f64());
+        }
+    }
+    // A long final estimate at the trained parameters.
+    vmc.cfg.samples = samples * 8;
+    let f = vmc.estimate();
+    println!("final  E {:.6} ± {:.6}  variance {:.2e}  relative error {:.2e} (mean of last fifth {:.6})", f.energy, f.error, f.variance, (f.energy - reference) / reference.abs(), tail.iter().sum::<f64>() / tail.len() as f64);
+}

wai/wai-quantum/src/lib.rs

@@ -127,6 +127,12 @@ pub mod quantum_dmrg;
 #[cfg(feature = "quantum_tdvp")]
 pub mod quantum_tdvp;
 
+#[cfg(feature = "quantum_chem")]
+pub mod quantum_chem;
+
+#[cfg(feature = "quantum_nqs")]
+pub mod quantum_nqs;
+
 #[cfg(feature = "quantum_frame")]
 pub mod quantum_match;
 

wai/wai-quantum/src/quantum_chem.rs

@@ -0,0 +1,607 @@
+//! 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);
+    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);
+    }
+}

wai/wai-quantum/src/quantum_nqs.rs

@@ -0,0 +1,922 @@
+//! Neural quantum states — `wai.quantum.nqs`.
+//!
+//! A neural network can stand in for a many-body wavefunction: it gives the
+//! amplitude of each configuration, and variational Monte Carlo trains it
+//! toward the ground state. Results from such networks are offered as
+//! evidence about systems no exact method reaches. This module is the method
+//! in its checkable form, deterministic and held to exact answers.
+//!
+//! - **Models.** [`SpinModel`] is a spin-½ Hamiltonian on any graph, with
+//!   exchange, Ising, and field terms. Presets are Heisenberg and
+//!   transverse-field Ising, on chains and square lattices.
+//! - **The network.** [`Rbm`] is a restricted Boltzmann machine with `α·n`
+//!   hidden units. On a bipartite lattice it can carry Marshall's sign, which
+//!   makes the Heisenberg ground state positive in what remains.
+//! - **Variational Monte Carlo.** [`Vmc`] samples `|ψ|²` by Metropolis
+//!   chains: single flips, or swaps of antiparallel spins at fixed
+//!   magnetisation. Or it enumerates every configuration (at most 24 spins)
+//!   for expectations with no Monte Carlo noise. Each step is stochastic
+//!   reconfiguration, `θ ← θ − η (S + λ)⁻¹ F`, solved in parameter space or
+//!   in sample space, whichever is smaller; the two give the same step.
+//! - **Referees.** [`SpinModel::exact_ground_energy`] is Lanczos on all `2ⁿ`
+//!   amplitudes (at most 22 spins). [`SpinModel::opsum`] gives the model to
+//!   DMRG on its exact MPO, for lattices past exact diagonalisation.
+//!
+//! Chains draw from their own seeded generators, and every sum runs in a
+//! fixed order, so a run is the same bits at any thread count.
+//!
+//! # Checked
+//!
+//! - **The parts.** The model's exact diagonalisation equals the chain
+//!   referees. Its operator sum's MPO equals its dense matrix on a 3 × 3
+//!   lattice with every kind of term. The network's amplitude ratios equal
+//!   ratios of its amplitudes to `10⁻¹²`, Marshall sign included, and its
+//!   log-derivatives agree with central differences. Both spaces take the
+//!   same reconfiguration step.
+//! - **Exact expectations.** On an 8-site Heisenberg chain, 200 steps with
+//!   every configuration enumerated reach the exact energy to `10⁻⁷`.
+//! - **Sampling.** Sampled energies agree with enumerated ones within their
+//!   error bars, and steps repeat bit for bit at any thread count.
+//! - **Against exact answers** (`nqs_ground` example; relative error of the
+//!   final estimate, with its statistical error):
+//!
+//!   | model             | lattice      | reference               | α | steps | relative error          |
+//!   |-------------------|--------------|-------------------------|---|-------|-------------------------|
+//!   | Heisenberg        | chain of 12  | exact diagonalisation   | 2 | 300   | 7.7·10⁻⁵ ± 7·10⁻⁵       |
+//!   | Heisenberg        | 4 × 4 open   | exact diagonalisation   | 2 | 600   | 4.6·10⁻⁴ ± 2·10⁻⁴       |
+//!   | Ising, `h = 3.044`| 4 × 4 open   | exact diagonalisation   | 1 | 400   | −3·10⁻⁶ ± 1·10⁻⁵        |
+//!   | Heisenberg        | 6 × 6 open   | DMRG, bond 256          | 2 | 800   | 3.6·10⁻³ ± 3·10⁻⁴       |
+//!
+//!   Each run used 1000 samples a step. The 6 × 6 DMRG reference is itself
+//!   an upper bound (discarded weight `2·10⁻⁵`). The network was still
+//!   descending when it stopped there.
+//!
+//! # Honest boundaries
+//!
+//! - **A real network.** The sign structure must be known, as Marshall's is
+//!   on bipartite lattices. Frustrated models, whose signs are not known,
+//!   need complex parameters or a learned sign, which are not here.
+//! - **Error bars ignore autocorrelation.** The reported error is the
+//!   standard error of uncorrelated samples. Chains thinned by two sweeps are
+//!   correlated, so the true error is larger.
+//! - **Dense reconfiguration.** `S` (or its sample-space twin) is formed and
+//!   factored in full each step. That bounds the network and sample sizes
+//!   long before iterative solvers would.
+
+use crate::quantum_dmrg::{Op, OpSum, default_threads, lanczos};
+use crate::repro::{exp, ln};
+
+// ---------------------------------------------------------------------------
+// Spin models
+// ---------------------------------------------------------------------------
+
+/// A spin-½ model on any graph, in Pauli matrices:
+/// `H = Σ J⊥ (σˣσˣ + σʸσʸ) + Σ J_z σᶻσᶻ + Σ hₓ σˣ + Σ h_z σᶻ`.
+/// Spin `+1` is `|↑⟩`, the first basis state.
+#[derive(Clone, Debug, PartialEq)]
+pub struct SpinModel {
+    pub n: usize,
+    /// `(i, j, J⊥)`: `J⊥ (σˣᵢσˣⱼ + σʸᵢσʸⱼ)`, which swaps antiparallel spins with
+    /// amplitude `2J⊥`.
+    pub exchange: Vec<(usize, usize, f64)>,
+    /// `(i, j, J_z)`: `J_z σᶻᵢσᶻⱼ`.
+    pub zz: Vec<(usize, usize, f64)>,
+    /// `(i, hₓ)`: `hₓ σˣᵢ`.
+    pub x: Vec<(usize, f64)>,
+    /// `(i, h_z)`: `h_z σᶻᵢ`.
+    pub z: Vec<(usize, f64)>,
+}
+
+/// The bonds of an open (or periodic) chain.
+pub fn chain_edges(n: usize, periodic: bool) -> Vec<(usize, usize)> {
+    let mut e: Vec<(usize, usize)> = (0..n.saturating_sub(1)).map(|i| (i, i + 1)).collect();
+    if periodic && n > 2 {
+        e.push((n - 1, 0));
+    }
+    e
+}
+
+/// The nearest-neighbour bonds of an `lx × ly` square lattice, site
+/// `y·lx + x`, open or periodic in both directions.
+pub fn square_edges(lx: usize, ly: usize, periodic: bool) -> Vec<(usize, usize)> {
+    let mut e = Vec::new();
+    for y in 0..ly {
+        for x in 0..lx {
+            let i = y * lx + x;
+            if x + 1 < lx {
+                e.push((i, i + 1));
+            } else if periodic && lx > 2 {
+                e.push((i, y * lx));
+            }
+            if y + 1 < ly {
+                e.push((i, i + lx));
+            } else if periodic && ly > 2 {
+                e.push((i, x));
+            }
+        }
+    }
+    e
+}
+
+/// The checkerboard sublattice of an `lx × ly` lattice (a chain is `ly = 1`).
+pub fn checkerboard(lx: usize, ly: usize) -> Vec<bool> {
+    (0..lx * ly).map(|i| (i % lx + i / lx).is_multiple_of(2)).collect()
+}
+
+impl SpinModel {
+    /// The Heisenberg model `J Σ Sᵢ·Sⱼ` (with `S = σ/2`) on `edges`.
+    pub fn heisenberg(n: usize, edges: &[(usize, usize)], j: f64) -> SpinModel {
+        SpinModel { n, exchange: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), zz: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), x: Vec::new(), z: Vec::new() }
+    }
+
+    /// The transverse-field Ising model `−J Σ σᶻσᶻ − h Σ σˣ` on `edges`.
+    pub fn ising(n: usize, edges: &[(usize, usize)], j: f64, h: f64) -> SpinModel {
+        SpinModel { n, exchange: Vec::new(), zz: edges.iter().map(|&(a, b)| (a, b, -j)).collect(), x: (0..n).map(|i| (i, -h)).collect(), z: Vec::new() }
+    }
+
+    /// Whether `Σ σᶻ` is conserved (no transverse field).
+    pub fn conserves_magnetisation(&self) -> bool {
+        self.x.is_empty()
+    }
+
+    /// The model as an operator sum, for DMRG on its MPO.
+    pub fn opsum(&self) -> OpSum {
+        let mut s = OpSum::new(self.n);
+        for &(i, j, c) in &self.exchange {
+            // σˣσˣ + σʸσʸ = 2(σ⁺σ⁻ + σ⁻σ⁺).
+            s.add(2.0 * c, &[(i, Op::plus()), (j, Op::minus())]);
+            s.add(2.0 * c, &[(i, Op::minus()), (j, Op::plus())]);
+        }
+        for &(i, j, c) in &self.zz {
+            s.add(c, &[(i, Op::z()), (j, Op::z())]);
+        }
+        for &(i, c) in &self.x {
+            s.add(c, &[(i, Op::x())]);
+        }
+        for &(i, c) in &self.z {
+            s.add(c, &[(i, Op::z())]);
+        }
+        s
+    }
+
+    /// `out = H inp` on all `2ⁿ` amplitudes, site 0 the most significant bit
+    /// (a set bit is `|↓⟩`).
+    fn apply_dense(&self, inp: &[f64], out: &mut [f64]) {
+        let n = self.n;
+        let bit = |i: usize| 1usize << (n - 1 - i);
+        let spin = |x: usize, i: usize| if x & bit(i) == 0 { 1.0 } else { -1.0 };
+        out.iter_mut().for_each(|o| *o = 0.0);
+        for (x, &amp) in inp.iter().enumerate() {
+            if amp == 0.0 {
+                continue;
+            }
+            let mut diag = 0.0;
+            for &(i, j, c) in &self.zz {
+                diag += c * spin(x, i) * spin(x, j);
+            }
+            for &(i, c) in &self.z {
+                diag += c * spin(x, i);
+            }
+            out[x] += diag * amp;
+            for &(i, j, c) in &self.exchange {
+                if spin(x, i) != spin(x, j) {
+                    out[x ^ bit(i) ^ bit(j)] += 2.0 * c * amp;
+                }
+            }
+            for &(i, c) in &self.x {
+                out[x ^ bit(i)] += c * amp;
+            }
+        }
+    }
+
+    /// The exact ground energy by Lanczos on all `2ⁿ` amplitudes (at most 22
+    /// sites), restricted to total `Σ σᶻ = magnetisation` when given (the
+    /// model must conserve it).
+    pub fn exact_ground_energy(&self, magnetisation: Option<i32>) -> Option<f64> {
+        let n = self.n;
+        if n == 0 || n > 22 {
+            return None;
+        }
+        if magnetisation.is_some() && !self.conserves_magnetisation() {
+            return None;
+        }
+        let dim = 1usize << n;
+        let in_sector = |x: usize| magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m);
+        let start: Vec<f64> = (0..dim as u64).map(|x| if in_sector(x as usize) { 1.0 + (x.wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0 } else { 0.0 }).collect();
+        if start.iter().all(|&v| v == 0.0) {
+            return None;
+        }
+        let apply = |inp: &[f64], out: &mut [f64]| self.apply_dense(inp, out);
+        Some(lanczos(&apply, &start, 60, 30, 1e-9).0)
+    }
+}
+
+// ---------------------------------------------------------------------------
+// The restricted Boltzmann machine
+// ---------------------------------------------------------------------------
+
+/// `ln cosh x`, stable for large `|x|`.
+fn lncosh(x: f64) -> f64 {
+    let a = x.abs();
+    a + ln(1.0 + exp(-2.0 * a)) - core::f64::consts::LN_2
+}
+
+fn tanh(x: f64) -> f64 {
+    let e = exp(-2.0 * x.abs());
+    let t = (1.0 - e) / (1.0 + e);
+    if x < 0.0 { -t } else { t }
+}
+
+/// A real restricted Boltzmann machine: `ψ(s) = sign(s) · exp(Σ aᵢsᵢ) ·
+/// Π_j cosh(b_j + Σ W_jᵢ sᵢ)`, with `m = α·n` hidden units. On a bipartite
+/// lattice the optional sign is Marshall's, `(−1)` to the number of up spins
+/// on sublattice A, which makes the Heisenberg ground state positive in the
+/// rest.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Rbm {
+    pub n: usize,
+    pub m: usize,
+    pub a: Vec<f64>,
+    pub b: Vec<f64>,
+    /// `W[j·n + i]`.
+    pub w: Vec<f64>,
+    pub sublattice: Option<Vec<bool>>,
+}
+
+/// A deterministic generator (splitmix64).
+#[derive(Clone, Debug, PartialEq)]
+struct Rng(u64);
+
+impl Rng {
+    fn next(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
+        z ^ (z >> 31)
+    }
+
+    /// Uniform on `[0, 1)`.
+    fn uniform(&mut self) -> f64 {
+        (self.next() >> 11) as f64 / 9_007_199_254_740_992.0
+    }
+
+    fn below(&mut self, k: usize) -> usize {
+        (self.next() % k as u64) as usize
+    }
+}
+
+impl Rbm {
+    /// `α·n` hidden units, weights uniform in `±scale`, biases zero.
+    pub fn new(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm {
+        let m = alpha.max(1) * n;
+        let mut rng = Rng(seed);
+        let w = (0..m * n).map(|_| scale * (2.0 * rng.uniform() - 1.0)).collect();
+        Rbm { n, m, a: vec![0.0; n], b: vec![0.0; m], w, sublattice: None }
+    }
+
+    /// Multiply by Marshall's sign for this sublattice.
+    pub fn with_marshall_sign(mut self, sublattice: Vec<bool>) -> Rbm {
+        assert_eq!(sublattice.len(), self.n);
+        self.sublattice = Some(sublattice);
+        self
+    }
+
+    /// The number of variational parameters.
+    pub fn params(&self) -> usize {
+        self.n + self.m + self.n * self.m
+    }
+
+    fn theta(&self, s: &[i8]) -> Vec<f64> {
+        (0..self.m)
+            .map(|j| {
+                let row = &self.w[j * self.n..(j + 1) * self.n];
+                let mut t = self.b[j];
+                for (wi, &si) in row.iter().zip(s) {
+                    t += wi * f64::from(si);
+                }
+                t
+            })
+            .collect()
+    }
+
+    /// `ln |ψ(s)|` up to a constant.
+    pub fn log_amplitude(&self, s: &[i8]) -> f64 {
+        let mut l: f64 = self.a.iter().zip(s).map(|(a, &si)| a * f64::from(si)).sum();
+        for t in self.theta(s) {
+            l += lncosh(t);
+        }
+        l
+    }
+
+    /// The sign of `ψ(s)`.
+    pub fn sign(&self, s: &[i8]) -> f64 {
+        match &self.sublattice {
+            Some(sub) => {
+                let ups = s.iter().zip(sub).filter(|&(&si, &a)| a && si > 0).count();
+                if ups.is_multiple_of(2) { 1.0 } else { -1.0 }
+            }
+            None => 1.0,
+        }
+    }
+
+    /// `ψ(s′)/ψ(s)` where `s′` is `s` with `flips` flipped, given `θ(s)`.
+    fn ratio(&self, s: &[i8], theta: &[f64], flips: &[usize]) -> f64 {
+        let mut l = 0.0;
+        for &i in flips {
+            l -= 2.0 * self.a[i] * f64::from(s[i]);
+        }
+        for (j, &t) in theta.iter().enumerate() {
+            let mut d = 0.0;
+            for &i in flips {
+                d -= 2.0 * self.w[j * self.n + i] * f64::from(s[i]);
+            }
+            l += lncosh(t + d) - lncosh(t);
+        }
+        let mut r = exp(l);
+        if let Some(sub) = &self.sublattice
+            && flips.iter().filter(|&&i| sub[i]).count() % 2 == 1
+        {
+            r = -r;
+        }
+        r
+    }
+
+    /// `∂ ln ψ / ∂(a, b, W)` at `s`, given `θ(s)`.
+    fn log_derivatives(&self, s: &[i8], theta: &[f64], out: &mut [f64]) {
+        let (n, m) = (self.n, self.m);
+        for (o, &si) in out[..n].iter_mut().zip(s) {
+            *o = f64::from(si);
+        }
+        for (j, &th) in theta.iter().enumerate() {
+            let t = tanh(th);
+            out[n + j] = t;
+            for (o, &si) in out[n + m + j * n..n + m + (j + 1) * n].iter_mut().zip(s) {
+                *o = t * f64::from(si);
+            }
+        }
+    }
+
+    fn shift(&mut self, delta: &[f64], step: f64) {
+        let (n, m) = (self.n, self.m);
+        for (x, d) in self.a.iter_mut().zip(&delta[..n]) {
+            *x -= step * d;
+        }
+        for (x, d) in self.b.iter_mut().zip(&delta[n..n + m]) {
+            *x -= step * d;
+        }
+        for (x, d) in self.w.iter_mut().zip(&delta[n + m..]) {
+            *x -= step * d;
+        }
+    }
+}
+
+/// `E_loc(s) = Σ_{s′} H_{ss′} ψ(s′)/ψ(s)`.
+fn local_energy(model: &SpinModel, rbm: &Rbm, s: &[i8], theta: &[f64]) -> f64 {
+    let mut e = 0.0;
+    for &(i, j, c) in &model.zz {
+        e += c * f64::from(s[i]) * f64::from(s[j]);
+    }
+    for &(i, c) in &model.z {
+        e += c * f64::from(s[i]);
+    }
+    for &(i, j, c) in &model.exchange {
+        if s[i] != s[j] {
+            e += 2.0 * c * rbm.ratio(s, theta, &[i, j]);
+        }
+    }
+    for &(i, c) in &model.x {
+        e += c * rbm.ratio(s, theta, &[i]);
+    }
+    e
+}
+
+// ---------------------------------------------------------------------------
+// Variational Monte Carlo
+// ---------------------------------------------------------------------------
+
+/// Settings for [`Vmc`].
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct VmcConfig {
+    /// Samples per iteration, over all chains.
+    pub samples: usize,
+    pub chains: usize,
+    /// Sweeps (`n` moves each) discarded when a chain starts.
+    pub burn: usize,
+    /// Sweeps between kept samples.
+    pub thin: usize,
+    /// Enumerate every configuration (at most 24 spins) instead of sampling:
+    /// exact expectations, no Monte Carlo noise.
+    pub exact: bool,
+    /// Keep `Σ σᶻ` at this value: moves swap antiparallel spins. The model
+    /// must conserve it.
+    pub magnetisation: Option<i32>,
+    /// Learning rate of stochastic reconfiguration.
+    pub lr: f64,
+    /// Diagonal shift added to the metric `S`.
+    pub shift: f64,
+    pub seed: u64,
+    /// Threads over chains and configurations; the result does not depend on
+    /// it.
+    pub threads: usize,
+}
+
+impl Default for VmcConfig {
+    fn default() -> VmcConfig {
+        VmcConfig { samples: 2000, chains: 8, burn: 50, thin: 2, exact: false, magnetisation: None, lr: 0.05, shift: 1e-3, seed: 1, threads: default_threads() }
+    }
+}
+
+/// One iteration's estimate.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Estimate {
+    /// `⟨E_loc⟩`.
+    pub energy: f64,
+    /// `⟨(E_loc − E)²⟩`: zero exactly at an eigenstate.
+    pub variance: f64,
+    /// The statistical error of `energy`, from the variance and the sample
+    /// count, ignoring autocorrelation (zero when exact).
+    pub error: f64,
+    /// Fraction of Metropolis moves accepted (one when exact).
+    pub acceptance: f64,
+}
+
+/// Variational Monte Carlo with stochastic reconfiguration on an RBM.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Vmc {
+    pub model: SpinModel,
+    pub rbm: Rbm,
+    pub cfg: VmcConfig,
+    chains: Vec<(Vec<i8>, Rng)>,
+}
+
+/// Samples: configurations, their weights (summing to one), local energies
+/// and log-derivatives (`params` per row).
+struct Batch {
+    weights: Vec<f64>,
+    energies: Vec<f64>,
+    derivs: Vec<f64>,
+    acceptance: f64,
+}
+
+impl Vmc {
+    pub fn new(model: SpinModel, rbm: Rbm, cfg: VmcConfig) -> Vmc {
+        assert_eq!(model.n, rbm.n, "the model and the machine differ in size");
+        if cfg.magnetisation.is_some() {
+            assert!(model.conserves_magnetisation(), "the model does not conserve the magnetisation");
+        }
+        let n = model.n;
+        let chains = (0..cfg.chains.max(1))
+            .map(|c| {
+                let mut rng = Rng(cfg.seed ^ (c as u64).wrapping_mul(0xa076_1d64_78bd_642f));
+                let s = start_configuration(n, cfg.magnetisation, &mut rng);
+                (s, rng)
+            })
+            .collect();
+        let mut v = Vmc { model, rbm, cfg, chains };
+        if !cfg.exact {
+            for c in 0..v.chains.len() {
+                let (mut s, mut rng) = v.chains[c].clone();
+                metropolis(&v.rbm, &mut s, &mut rng, cfg.burn * n, cfg.magnetisation.is_some());
+                v.chains[c] = (s, rng);
+            }
+        }
+        v
+    }
+
+    fn batch(&mut self) -> Batch {
+        if self.cfg.exact { self.enumerate() } else { self.sample() }
+    }
+
+    fn enumerate(&self) -> Batch {
+        let n = self.model.n;
+        assert!(n <= 24, "exact enumeration takes at most 24 spins");
+        let configs: Vec<usize> = (0..1usize << n).filter(|&x| self.cfg.magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m)).collect();
+        let p = self.rbm.params();
+        let k = configs.len();
+        let mut logs = vec![0.0; k];
+        let mut energies = vec![0.0; k];
+        let mut derivs = vec![0.0; k * p];
+        let threads = self.cfg.threads.max(1).min(k.div_ceil(64)).max(1);
+        let per = k.div_ceil(threads);
+        let work = |range: core::ops::Range<usize>, logs: &mut [f64], energies: &mut [f64], derivs: &mut [f64]| {
+            let start = range.start;
+            for idx in range {
+                let x = configs[idx];
+                let s: Vec<i8> = (0..n).map(|i| if x >> (n - 1 - i) & 1 == 0 { 1 } else { -1 }).collect();
+                let theta = self.rbm.theta(&s);
+                let r = idx - start;
+                logs[r] = self.rbm.log_amplitude(&s);
+                energies[r] = local_energy(&self.model, &self.rbm, &s, &theta);
+                self.rbm.log_derivatives(&s, &theta, &mut derivs[r * p..(r + 1) * p]);
+            }
+        };
+        if threads <= 1 || cfg!(target_arch = "wasm32") {
+            work(0..k, &mut logs, &mut energies, &mut derivs);
+        } else {
+            std::thread::scope(|scope| {
+                let mut lrest: &mut [f64] = &mut logs;
+                let mut erest: &mut [f64] = &mut energies;
+                let mut drest: &mut [f64] = &mut derivs;
+                for t in 0..threads {
+                    let lo = t * per;
+                    let hi = ((t + 1) * per).min(k);
+                    if lo >= hi {
+                        break;
+                    }
+                    let (l, lr) = lrest.split_at_mut(hi - lo);
+                    let (e, er) = erest.split_at_mut(hi - lo);
+                    let (d, dr) = drest.split_at_mut((hi - lo) * p);
+                    lrest = lr;
+                    erest = er;
+                    drest = dr;
+                    let work = &work;
+                    scope.spawn(move || work(lo..hi, l, e, d));
+                }
+            });
+        }
+        // |ψ|², normalised, in a fixed order.
+        let top = logs.iter().copied().fold(f64::NEG_INFINITY, f64::max);
+        let raw: Vec<f64> = logs.iter().map(|l| exp(2.0 * (l - top))).collect();
+        let z: f64 = raw.iter().sum();
+        Batch { weights: raw.iter().map(|r| r / z).collect(), energies, derivs, acceptance: 1.0 }
+    }
+
+    fn sample(&mut self) -> Batch {
+        let n = self.model.n;
+        let p = self.rbm.params();
+        let nc = self.chains.len();
+        let per_chain = self.cfg.samples.div_ceil(nc).max(1);
+        let exchange = self.cfg.magnetisation.is_some();
+        let thin = self.cfg.thin.max(1) * n;
+        let rbm = &self.rbm;
+        let model = &self.model;
+        let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<f64>, Vec<f64>, f64) {
+            let (s, rng) = chain;
+            let mut energies = Vec::with_capacity(per_chain);
+            let mut derivs = vec![0.0; per_chain * p];
+            let mut accepted = 0usize;
+            for k in 0..per_chain {
+                accepted += metropolis(rbm, s, rng, thin, exchange);
+                let theta = rbm.theta(s);
+                energies.push(local_energy(model, rbm, s, &theta));
+                rbm.log_derivatives(s, &theta, &mut derivs[k * p..(k + 1) * p]);
+            }
+            (energies, derivs, accepted as f64 / (per_chain * thin) as f64)
+        };
+        let threads = self.cfg.threads.max(1).min(nc);
+        let results: Vec<(Vec<f64>, Vec<f64>, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") {
+            self.chains.iter_mut().map(run).collect()
+        } else {
+            let per = nc.div_ceil(threads);
+            std::thread::scope(|scope| {
+                let handles: Vec<_> = self
+                    .chains
+                    .chunks_mut(per)
+                    .map(|group| {
+                        let run = &run;
+                        scope.spawn(move || group.iter_mut().map(run).collect::<Vec<_>>())
+                    })
+                    .collect();
+                handles.into_iter().flat_map(|h| h.join().expect("a sampling thread panicked")).collect()
+            })
+        };
+        let total = per_chain * nc;
+        let mut energies = Vec::with_capacity(total);
+        let mut derivs = Vec::with_capacity(total * p);
+        let mut acc = 0.0;
+        for (e, d, a) in results {
+            energies.extend(e);
+            derivs.extend(d);
+            acc += a;
+        }
+        Batch { weights: vec![1.0 / total as f64; total], energies, derivs, acceptance: acc / nc as f64 }
+    }
+
+    /// Estimate the energy without changing the machine.
+    pub fn estimate(&mut self) -> Estimate {
+        let b = self.batch();
+        summarise(&b, self.cfg.exact)
+    }
+
+    /// One step of stochastic reconfiguration: `θ ← θ − η (S + λ)⁻¹ F`, with
+    /// `S` the covariance of the log-derivatives and `F` their covariance
+    /// with the local energy. It is solved in parameter space or in sample
+    /// space, whichever is smaller; the two give the same step.
+    pub fn step(&mut self) -> Estimate {
+        let b = self.batch();
+        let est = summarise(&b, self.cfg.exact);
+        let delta = sr_direction(&b, self.rbm.params(), self.cfg.shift);
+        self.rbm.shift(&delta, self.cfg.lr);
+        est
+    }
+}
+
+fn summarise(b: &Batch, exact: bool) -> Estimate {
+    let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e).sum();
+    let variance: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * (e - energy) * (e - energy)).sum();
+    let error = if exact { 0.0 } else { (variance / b.energies.len() as f64).sqrt() };
+    Estimate { energy, variance, error, acceptance: b.acceptance }
+}
+
+/// `(S + λ)⁻¹ F` from a batch.
+fn sr_direction(b: &Batch, p: usize, shift: f64) -> Vec<f64> {
+    let k = b.weights.len();
+    let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e).sum();
+    let mut mean = vec![0.0; p];
+    for (r, &w) in b.weights.iter().enumerate() {
+        for (m, &d) in mean.iter_mut().zip(&b.derivs[r * p..(r + 1) * p]) {
+            *m += w * d;
+        }
+    }
+    // Rows √w (O − Ō) and √w (E − Ē).
+    let mut o = vec![0.0; k * p];
+    let mut eps = vec![0.0; k];
+    for r in 0..k {
+        let sw = b.weights[r].sqrt();
+        for c in 0..p {
+            o[r * p + c] = sw * (b.derivs[r * p + c] - mean[c]);
+        }
+        eps[r] = sw * (b.energies[r] - energy);
+    }
+    if p <= k { sr_parameter_space(&o, &eps, k, p, shift) } else { sr_sample_space(&o, eps, k, p, shift) }
+}
+
+/// `(OᵀO + λ)⁻¹ Oᵀε`.
+fn sr_parameter_space(o: &[f64], eps: &[f64], k: usize, p: usize, shift: f64) -> Vec<f64> {
+    {
+        // (OᵀO + λ) δ = Oᵀε.
+        let mut s = vec![0.0; p * p];
+        for r in 0..k {
+            let row = &o[r * p..(r + 1) * p];
+            for i in 0..p {
+                let oi = row[i];
+                if oi == 0.0 {
+                    continue;
+                }
+                for j in i..p {
+                    s[i * p + j] += oi * row[j];
+                }
+            }
+        }
+        for i in 0..p {
+            s[i * p + i] += shift;
+            for j in 0..i {
+                s[i * p + j] = s[j * p + i];
+            }
+        }
+        let mut f = vec![0.0; p];
+        for r in 0..k {
+            for c in 0..p {
+                f[c] += o[r * p + c] * eps[r];
+            }
+        }
+        cholesky_solve(&mut s, &mut f, p);
+        f
+    }
+}
+
+/// `Oᵀ (OOᵀ + λ)⁻¹ ε`.
+fn sr_sample_space(o: &[f64], eps: Vec<f64>, k: usize, p: usize, shift: f64) -> Vec<f64> {
+    {
+        // δ = Oᵀ (OOᵀ + λ)⁻¹ ε.
+        let mut t = vec![0.0; k * k];
+        for i in 0..k {
+            for j in i..k {
+                let mut acc = 0.0;
+                for c in 0..p {
+                    acc += o[i * p + c] * o[j * p + c];
+                }
+                t[i * k + j] = acc;
+                t[j * k + i] = acc;
+            }
+            t[i * k + i] += shift;
+        }
+        let mut x = eps;
+        cholesky_solve(&mut t, &mut x, k);
+        let mut delta = vec![0.0; p];
+        for r in 0..k {
+            for c in 0..p {
+                delta[c] += o[r * p + c] * x[r];
+            }
+        }
+        delta
+    }
+}
+
+/// Solve `A x = b` for symmetric positive-definite `A` (overwritten by its
+/// Cholesky factor), `b` overwritten by `x`.
+#[allow(clippy::needless_range_loop)]
+fn cholesky_solve(a: &mut [f64], b: &mut [f64], n: usize) {
+    for j in 0..n {
+        let mut d = a[j * n + j];
+        for k in 0..j {
+            d -= a[j * n + k] * a[j * n + k];
+        }
+        let d = d.max(1e-300).sqrt();
+        a[j * n + j] = d;
+        for i in j + 1..n {
+            let mut s = a[i * n + j];
+            for k in 0..j {
+                s -= a[i * n + k] * a[j * n + k];
+            }
+            a[i * n + j] = s / d;
+        }
+    }
+    for i in 0..n {
+        let mut s = b[i];
+        for k in 0..i {
+            s -= a[i * n + k] * b[k];
+        }
+        b[i] = s / a[i * n + i];
+    }
+    for i in (0..n).rev() {
+        let mut s = b[i];
+        for k in i + 1..n {
+            s -= a[k * n + i] * b[k];
+        }
+        b[i] = s / a[i * n + i];
+    }
+}
+
+/// A random configuration, at the given magnetisation if any.
+fn start_configuration(n: usize, magnetisation: Option<i32>, rng: &mut Rng) -> Vec<i8> {
+    match magnetisation {
+        None => (0..n).map(|_| if rng.next() & 1 == 0 { 1 } else { -1 }).collect(),
+        Some(m) => {
+            let ups = ((n as i32 + m) / 2).clamp(0, n as i32) as usize;
+            let mut s: Vec<i8> = (0..n).map(|i| if i < ups { 1 } else { -1 }).collect();
+            // Fisher–Yates.
+            for i in (1..n).rev() {
+                let j = rng.below(i + 1);
+                s.swap(i, j);
+            }
+            s
+        }
+    }
+}
+
+/// `moves` Metropolis moves on `s` (single flips, or swaps of antiparallel
+/// spins), accepted with probability `|ψ(s′)/ψ(s)|²`. Returns the number
+/// accepted.
+fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool) -> usize {
+    let n = s.len();
+    let mut theta = rbm.theta(s);
+    let mut accepted = 0;
+    for _ in 0..moves {
+        let flips: Vec<usize> = if exchange {
+            let i = rng.below(n);
+            let j = rng.below(n);
+            if s[i] == s[j] {
+                continue;
+            }
+            vec![i, j]
+        } else {
+            vec![rng.below(n)]
+        };
+        let r = rbm.ratio(s, &theta, &flips);
+        if rng.uniform() < r * r {
+            for &i in &flips {
+                let ds = -2.0 * f64::from(s[i]);
+                for (j, t) in theta.iter_mut().enumerate() {
+                    *t += rbm.w[j * n + i] * ds;
+                }
+                s[i] = -s[i];
+            }
+            accepted += 1;
+        }
+    }
+    accepted
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::quantum_dmrg::{Chain, exact_ground_energy};
+
+    #[test]
+    fn exact_diagonalisation_matches_the_chain_referees() {
+        let n = 10;
+        let e = chain_edges(n, false);
+        let a = SpinModel::heisenberg(n, &e, 1.0).exact_ground_energy(Some(0)).unwrap();
+        let b = exact_ground_energy(&Chain::heisenberg(n, 1.0, 1.0, 0.0)).unwrap();
+        assert!((a - b).abs() < 1e-9, "{a} vs {b}");
+        let a = SpinModel::ising(n, &e, 1.0, 0.7).exact_ground_energy(None).unwrap();
+        let b = exact_ground_energy(&Chain::ising(n, 1.0, 0.7)).unwrap();
+        assert!((a - b).abs() < 1e-9, "{a} vs {b}");
+    }
+
+    #[test]
+    fn the_operator_sum_is_the_model() {
+        // A 3 × 3 lattice with every kind of term.
+        let n = 9;
+        let mut m = SpinModel::heisenberg(n, &square_edges(3, 3, false), 0.8);
+        m.x.push((4, 0.3));
+        m.z.push((2, -0.45));
+        m.zz.push((0, 8, 0.2));
+        let dense = m.opsum().mpo().to_dense().unwrap();
+        let dim = 1usize << n;
+        let mut col = vec![0.0; dim];
+        let mut out = vec![0.0; dim];
+        for c in 0..dim {
+            col.iter_mut().for_each(|x| *x = 0.0);
+            col[c] = 1.0;
+            m.apply_dense(&col, &mut out);
+            for r in 0..dim {
+                assert!((out[r] - dense[r * dim + c]).abs() < 1e-12, "({r}, {c})");
+            }
+        }
+    }
+
+    #[test]
+    fn ratios_and_derivatives_are_the_machines() {
+        let n = 7;
+        let mut rbm = Rbm::new(n, 2, 0.4, 3).with_marshall_sign(checkerboard(n, 1));
+        let mut rng = Rng(5);
+        for v in rbm.a.iter_mut().chain(rbm.b.iter_mut()) {
+            *v = 0.3 * (rng.uniform() - 0.5);
+        }
+        for _ in 0..20 {
+            let s = start_configuration(n, None, &mut rng);
+            let theta = rbm.theta(&s);
+            for flips in [vec![rng.below(n)], vec![0, 3], vec![2, 6]] {
+                let mut t = s.clone();
+                for &i in &flips {
+                    t[i] = -t[i];
+                }
+                let want = rbm.sign(&t) * rbm.sign(&s) * exp(rbm.log_amplitude(&t) - rbm.log_amplitude(&s));
+                let got = rbm.ratio(&s, &theta, &flips);
+                assert!((got - want).abs() < 1e-12 * want.abs().max(1.0), "{got} vs {want}");
+            }
+            // Derivatives by central differences.
+            let p = rbm.params();
+            let mut d = vec![0.0; p];
+            rbm.log_derivatives(&s, &theta, &mut d);
+            let h = 1e-6;
+            for k in [0, n, n + 3, n + rbm.m + 5, p - 1] {
+                let mut e = vec![0.0; p];
+                e[k] = -h;
+                let mut up = rbm.clone();
+                up.shift(&e, 1.0);
+                e[k] = h;
+                let mut dn = rbm.clone();
+                dn.shift(&e, 1.0);
+                let fd = (up.log_amplitude(&s) - dn.log_amplitude(&s)) / (2.0 * h);
+                assert!((fd - d[k]).abs() < 1e-7, "param {k}: {fd} vs {}", d[k]);
+            }
+        }
+    }
+
+    #[test]
+    fn both_spaces_take_the_same_step() {
+        let mut rng = Rng(9);
+        for (k, p) in [(12usize, 5usize), (5, 12)] {
+            let o: Vec<f64> = (0..k * p).map(|_| rng.uniform() - 0.5).collect();
+            let eps: Vec<f64> = (0..k).map(|_| rng.uniform() - 0.5).collect();
+            let a = sr_parameter_space(&o, &eps, k, p, 1e-3);
+            let b = sr_sample_space(&o, eps.clone(), k, p, 1e-3);
+            for (x, y) in a.iter().zip(&b) {
+                assert!((x - y).abs() < 1e-9 * x.abs().max(1.0), "{x} vs {y}");
+            }
+        }
+    }
+
+    #[test]
+    fn exact_vmc_finds_the_ground_state() {
+        let n = 8;
+        let model = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0);
+        let exact = model.exact_ground_energy(Some(0)).unwrap();
+        let rbm = Rbm::new(n, 2, 0.05, 1).with_marshall_sign(checkerboard(n, 1));
+        let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-4, ..VmcConfig::default() };
+        let mut vmc = Vmc::new(model, rbm, cfg);
+        let mut last = vmc.step();
+        for _ in 0..300 {
+            let e = vmc.step();
+            last = e;
+        }
+        let rel = (last.energy - exact) / exact.abs();
+        assert!(last.energy >= exact - 1e-9, "variational: {} vs {exact}", last.energy);
+        assert!(rel < 1e-3, "{} vs {exact} (relative {rel:e})", last.energy);
+        assert!(last.variance < 1e-2, "{}", last.variance);
+    }
+
+    #[test]
+    fn sampling_agrees_with_enumeration_and_repeats() {
+        let n = 8;
+        let model = SpinModel::ising(n, &chain_edges(n, false), 1.0, 1.0);
+        let rbm = Rbm::new(n, 1, 0.3, 4);
+        let exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, ..VmcConfig::default() }).estimate();
+        let cfg = VmcConfig { samples: 8000, chains: 8, threads: 1, ..VmcConfig::default() };
+        let sampled = Vmc::new(model.clone(), rbm.clone(), cfg).estimate();
+        assert!((sampled.energy - exact.energy).abs() < 5.0 * sampled.error, "{sampled:?} vs {exact:?}");
+        // The same bits at any thread count, and again.
+        let mut one = Vmc::new(model.clone(), rbm.clone(), VmcConfig { samples: 400, threads: 1, ..VmcConfig::default() });
+        let mut many = Vmc::new(model, rbm, VmcConfig { samples: 400, threads: 3, ..VmcConfig::default() });
+        for _ in 0..3 {
+            assert_eq!(one.step(), many.step());
+        }
+        assert_eq!(one.rbm, many.rbm);
+        assert_eq!(one.chains, many.chains);
+    }
+}
+

wai/wai-quantum/src/quantum_tdvp.rs

@@ -106,7 +106,7 @@
 //!   times the cost.
 
 use crate::linalg::jacobi_svd_strict;
-use crate::quantum_dmrg::{C, Chain, Engine, Mpo, Mps, Op, Scalar, apply_dense, default_threads, entropy, hermitian_eigen, right_canonicalise, symmetric_eigen, transpose, two_site};
+use crate::quantum_dmrg::{C, Chain, Engine, Mpo, Mps, Op, Scalar, apply_dense, complete_bases, default_threads, entropy, hermitian_eigen, right_canonicalise, symmetric_eigen, transpose, two_site};
 use crate::repro::sin_cos;
 
 // ---------------------------------------------------------------------------
@@ -298,7 +298,7 @@ impl Evolution {
         let mut rights: Vec<Vec<C>> = vec![Vec::new(); n + 1];
         rights[n] = eng.right_edge();
         for i in (1..n).rev() {
-            rights[i] = eng.grow_right(&rights[i + 1], &mps.sites[i], mps.dims[i], mps.dims[i + 1]);
+            rights[i] = eng.grow_right(i, &rights[i + 1], &mps.sites[i], mps.dims[i], mps.dims[i + 1]);
         }
         Evolution { mpo, mps, rights, cfg, steps: 0 }
     }
@@ -327,7 +327,7 @@ impl Evolution {
             let (dl, dm, dr) = (mps.dims[i], mps.dims[i + 1], mps.dims[i + 2]);
             let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr);
             let (l, r) = (&lefts[i], &rights[i + 2]);
-            let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol);
+            let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol);
             let (u, mut rest, disc, weights) = svd_split(&theta, dl * 2, 2 * dr, cfg.max_bond, cfg.cutoff);
             if disc > 0.0 {
                 renormalise(&mut rest, norm(&theta));
@@ -336,10 +336,10 @@ impl Evolution {
             let k = weights.len();
             mps.sites[i] = u;
             mps.dims[i + 1] = k;
-            lefts[i + 1] = eng.grow_left(&lefts[i], &mps.sites[i], dl, k);
+            lefts[i + 1] = eng.grow_left(i, &lefts[i], &mps.sites[i], dl, k);
             mps.sites[i + 1] = if i + 2 < n {
                 let (l, r) = (&lefts[i + 1], &rights[i + 2]);
-                expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(l, r, x, y, k, dr, cfg.threads), &rest, -half, cfg.krylov, cfg.tol)
+                expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(i + 1, l, r, x, y, k, dr, cfg.threads), &rest, -half, cfg.krylov, cfg.tol)
             } else {
                 rest
             };
@@ -348,20 +348,20 @@ impl Evolution {
             let (dl, dm, dr) = (mps.dims[i], mps.dims[i + 1], mps.dims[i + 2]);
             let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr);
             let (l, r) = (&lefts[i], &rights[i + 2]);
-            let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol);
+            let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol);
             let (u, rest, disc, weights) = svd_split(&transpose(&theta, dl * 2, 2 * dr), 2 * dr, dl * 2, cfg.max_bond, cfg.cutoff);
             discarded += disc;
             let k = weights.len();
             mps.sites[i + 1] = transpose(&u, 2 * dr, k);
             mps.dims[i + 1] = k;
-            rights[i + 1] = eng.grow_right(&rights[i + 2], &mps.sites[i + 1], k, dr);
+            rights[i + 1] = eng.grow_right(i + 1, &rights[i + 2], &mps.sites[i + 1], k, dr);
             let mut centre = transpose(&rest, k, dl * 2);
             if disc > 0.0 {
                 renormalise(&mut centre, norm(&theta));
             }
             mps.sites[i] = if i > 0 {
                 let (l, r) = (&lefts[i], &rights[i + 1]);
-                expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(l, r, x, y, dl, k, cfg.threads), &centre, -half, cfg.krylov, cfg.tol)
+                expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(i, l, r, x, y, dl, k, cfg.threads), &centre, -half, cfg.krylov, cfg.tol)
             } else {
                 centre
             };
@@ -393,7 +393,7 @@ impl Evolution {
         let c = &self.mps.sites[0];
         let (dl, dr) = (self.mps.dims[0], self.mps.dims[1]);
         let mut hc = vec![C::ZERO; c.len()];
-        eng.apply1(&eng.left_edge(), &self.rights[1], c, &mut hc, dl, dr, self.cfg.threads);
+        eng.apply1(0, &eng.left_edge(), &self.rights[1], c, &mut hc, dl, dr, self.cfg.threads);
         dotc(c, &hc).re / dotc(c, c).re
     }
 
@@ -492,56 +492,6 @@ impl Evolution {
     }
 }
 
-/// Complete the orthonormal rows of every right-canonical site, from the
-/// right, up to the bond the chain allows or `max_bond`, and give the site to
-/// the left zero columns for them. The state does not change; the tangent
-/// space it spans grows, so a state with small bonds (a product state, say)
-/// is not held back by bases that cannot yet express `H|ψ⟩`.
-fn complete_bases(mps: &mut Mps<C>, max_bond: usize) {
-    let n = mps.n();
-    for i in (1..n).rev() {
-        let (k, dr) = (mps.dims[i], mps.dims[i + 1]);
-        let cols = 2 * dr;
-        let left_space = 1usize.checked_shl(i.min(40) as u32).unwrap_or(usize::MAX);
-        let target = cols.min(left_space).min(max_bond.max(k));
-        if target <= k {
-            continue;
-        }
-        let mut rows: Vec<Vec<C>> = mps.sites[i].chunks(cols).map(|r| r.to_vec()).collect();
-        for unit in 0..cols {
-            if rows.len() == target {
-                break;
-            }
-            let mut v = vec![C::ZERO; cols];
-            v[unit] = C::ONE;
-            // Gram–Schmidt against the rows so far, twice.
-            for _ in 0..2 {
-                for r in &rows {
-                    let c = dotc(r, &v);
-                    for (x, y) in v.iter_mut().zip(r) {
-                        *x = x.sub(c.mul(*y));
-                    }
-                }
-            }
-            let nv = norm(&v);
-            if nv > 0.5 {
-                rows.push(v.iter().map(|x| x.scale(1.0 / nv)).collect());
-            }
-        }
-        let added = rows.len();
-        mps.sites[i] = rows.into_iter().flatten().collect();
-        // The site to the left gains zero columns.
-        let pl = mps.dims[i - 1];
-        let prev = &mps.sites[i - 1];
-        let mut wider = vec![C::ZERO; pl * 2 * added];
-        for r in 0..pl * 2 {
-            wider[r * added..r * added + k].copy_from_slice(&prev[r * k..(r + 1) * k]);
-        }
-        mps.sites[i - 1] = wider;
-        mps.dims[i] = added;
-    }
-}
-
 /// Split `m` (rows × cols) by its singular value decomposition, keeping at
 /// most `max_bond` (and never more than its rank allows) and dropping weight
 /// below `cutoff`. Returns the kept left vectors (rows × k), `u† m`
@@ -1052,7 +1002,7 @@ mod tests {
     fn threaded_complex_effective_hamiltonians_are_bit_identical() {
         let mpo = Chain::heisenberg(8, 1.0, 0.7, 0.2).mpo();
         let eng = Engine::new(&mpo);
-        let (dl, dr, w) = (56usize, 48usize, mpo.w);
+        let (dl, dr, w) = (56usize, 48usize, mpo.dims[2]);
         let mut next = rng(3);
         let mut c = || C { re: next(), im: next() };
         let l: Vec<C> = (0..dl * w * dl).map(|_| c()).collect();
@@ -1062,13 +1012,13 @@ mod tests {
         let r1: Vec<C> = r.clone();
         let mut one = vec![C::ZERO; theta.len()];
         let mut one1 = vec![C::ZERO; site.len()];
-        eng.apply(&l, &r, &theta, &mut one, dl, dr, 1);
-        eng.apply1(&l, &r1, &site, &mut one1, dl, dr, 1);
+        eng.apply(2, &l, &r, &theta, &mut one, dl, dr, 1);
+        eng.apply1(2, &l, &r1, &site, &mut one1, dl, dr, 1);
         for t in [2, 3, 7] {
             let mut many = vec![C::ZERO; theta.len()];
             let mut many1 = vec![C::ZERO; site.len()];
-            eng.apply(&l, &r, &theta, &mut many, dl, dr, t);
-            eng.apply1(&l, &r1, &site, &mut many1, dl, dr, t);
+            eng.apply(2, &l, &r, &theta, &mut many, dl, dr, t);
+            eng.apply1(2, &l, &r1, &site, &mut many1, dl, dr, t);
             assert_eq!(one, many, "threads={t}");
             assert_eq!(one1, many1, "threads={t}");
         }

wai/wai-quantum/src/repro.rs

@@ -101,7 +101,7 @@ pub fn ln(x: f64) -> f64 {
 // Exponential
 // ---------------------------------------------------------------------------
 
-#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))]
+#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
 const EXP_P: [f64; 5] = [
     1.666_666_666_666_660_190_37e-1,
     -2.777_777_777_701_559_338_42e-3,
@@ -113,7 +113,7 @@ const EXP_P: [f64; 5] = [
 /// `eˣ`, bit-identical on every IEEE-754 platform: fdlibm's algorithm
 /// (reduction by `ln 2` in two parts, a rational approximation on
 /// `[−ln2/2, ln2/2]`), using only `+ − × ÷`. Within an ulp of the true value.
-#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))]
+#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
 pub fn exp(x: f64) -> f64 {
     const LN2_HI_E: f64 = 6.931_471_803_691_238_164_90e-1;
     const LN2_LO_E: f64 = 1.908_214_929_270_587_700_02e-10;
@@ -314,7 +314,7 @@ mod tests {
     #[allow(unused_imports)]
     use super::*;
 
-    #[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio"))]
+    #[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
     #[test]
     fn exp_is_within_an_ulp_of_the_platform() {
         let mut s = 0x9e37_79b9_7f4a_7c15u64;

wai/wai-quantum/tests/bit_identity.rs

@@ -26,7 +26,9 @@
     feature = "quantum_uf",
     feature = "quantum_noisy",
     feature = "quantum_dmrg",
-    feature = "quantum_tdvp"
+    feature = "quantum_tdvp",
+    feature = "quantum_chem",
+    feature = "quantum_nqs"
 ))]
 
 use wai_quantum::*;
@@ -473,7 +475,7 @@ fn dmrg_ground_states() {
     let (vals, vecs) = symmetric_eigen(&sym, 6);
     v.extend(vals);
     v.extend(vecs.iter().map(|x| x.abs()));
-    assert_eq!(digest(&v), 0x3b28_8bdd_be3f_48a2);
+    assert_eq!(digest(&v), 0xc0a7_f146_6ad9_b7e8);
 }
 
 #[test]
@@ -508,3 +510,42 @@ fn tdvp_time_evolution() {
     assert_eq!(digest(&v), 0xe1b4_fd28_af23_50cd);
 }
 
+#[test]
+fn chemistry_by_dmrg() {
+    use quantum_chem::*;
+    let mut v = Vec::new();
+    let fd = read_fcidump(include_str!("data/lih.fcidump")).unwrap();
+    v.push(fci(&fd).unwrap().energy);
+    v.push(hartree_fock_energy(&fd));
+    let r = ground_state(&fd, &ChemConfig { max_bond: 10, sweeps: 1, threads: 1, ..ChemConfig::default() });
+    v.extend(r.sweep_energies);
+    v.push(r.discarded);
+    v.extend(r.entropies);
+    let mpo = molecular_mpo(&fd, &ChemConfig::default());
+    v.extend(mpo.dims.iter().map(|&d| d as f64));
+    assert_eq!(digest(&v), 0xe955_a4ff_7eda_9960);
+}
+
+#[test]
+fn neural_quantum_states() {
+    use quantum_nqs::*;
+    let mut v = Vec::new();
+    let n = 8;
+    let model = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0);
+    v.push(model.exact_ground_energy(Some(0)).unwrap());
+    let rbm = Rbm::new(n, 2, 0.05, 1).with_marshall_sign(checkerboard(n, 1));
+    // Exact expectations, then sampled ones, a few steps each.
+    let mut exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, threads: 1, ..VmcConfig::default() });
+    for _ in 0..5 {
+        let e = exact.step();
+        v.extend([e.energy, e.variance]);
+    }
+    let mut sampled = Vmc::new(SpinModel::ising(n, &chain_edges(n, false), 1.0, 0.8), Rbm::new(n, 1, 0.1, 2), VmcConfig { samples: 200, chains: 4, threads: 1, ..VmcConfig::default() });
+    for _ in 0..3 {
+        let e = sampled.step();
+        v.extend([e.energy, e.variance, e.error, e.acceptance]);
+    }
+    v.extend(&sampled.rbm.w);
+    assert_eq!(digest(&v), 0x6d09_94ff_8d98_0dc1);
+}
+

wai/wai-quantum/tests/data/h2.fcidump

@@ -0,0 +1,12 @@
+ &FCI NORB=   2,NELEC= 2,MS2=0,
+  ORBSYM=1,1,
+  ISYM=1,
+ &END
+ 0.6747559268144484    1    1    1    1
+ 0.6637114013508132    1    1    2    2
+ 0.1812104620151968    2    1    2    1
+ 0.6637114013508132    2    2    1    1
+ 0.697651504490461    2    2    2    2
+ -1.253309786645977    1    1  0  0
+ -0.4750688487721783    2    2  0  0
+ 0.7151043390810812  0  0  0  0

wai/wai-quantum/tests/data/h2o.fcidump

@@ -0,0 +1,400 @@
+ &FCI NORB=   7,NELEC=10,MS2=0,
+  ORBSYM=1,1,1,1,1,1,1,
+  ISYM=1,
+ &END
+ 4.74450531202799    1    1    1    1
+ -0.4166568125051141    1    1    2    1
+ 1.004574901946536    1    1    2    2
+ 0.7999304193144103    1    1    3    3
+ -0.18359671151859    1    1    4    1
+ 0.1284623636944898    1    1    4    2
+ 0.9997707458540169    1    1    4    4
+ -1.096680936080461e-15    1    1    5    1
+ -1.234611383062098e-15    1    1    5    4
+ 1.115336271961169    1    1    5    5
+ -0.2379642076436868    1    1    6    1
+ 0.3082974147602019    1    1    6    2
+ 6.912715135774105e-15    1    1    6    3
+ -0.2194080179749168    1    1    6    4
+ -1.555242604178038e-15    1    1    6    5
+ 0.8026929114041939    1    1    6    6
+ 3.33333630322981e-15    1    1    7    1
+ -5.466453597376032e-15    1    1    7    2
+ 0.3624115212778842    1    1    7    3
+ 4.647095366902334e-15    1    1    7    4
+ 1.662523270532883e-15    1    1    7    6
+ 0.8689352741666503    1    1    7    7
+ -0.4166568125051136    2    1    1    1
+ 0.05817684869785823    2    1    2    1
+ -0.01297289116079521    2    1    2    2
+ -0.004406283490518009    2    1    3    3
+ 0.02249791049099267    2    1    4    1
+ -0.009212365591395872    2    1    4    2
+ -0.01356454804023788    2    1    4    4
+ -0.01169439676993425    2    1    5    5
+ 0.03579516654116417    2    1    6    1
+ -0.006646269777667551    2    1    6    2
+ 0.002234961476987691    2    1    6    4
+ -0.006978989180611147    2    1    6    6
+ -0.00750338145384954    2    1    7    3
+ -0.009399378211486308    2    1    7    7
+ 1.004574901946536    2    2    1    1
+ -0.01297289116079522    2    2    2    1
+ 0.7281701524927489    2    2    2    2
+ 0.6451438166645459    2    2    3    3
+ -0.01605176798440034    2    2    4    1
+ -0.004070151432527679    2    2    4    2
+ 0.6756423771903608    2    2    4    4
+ 0.7474056835529891    2    2    5    5
+ -0.000782432844683919    2    2    6    1
+ 0.1429042097541407    2    2    6    2
+ 2.603744830074969e-15    2    2    6    3
+ -0.09546441068766685    2    2    6    4
+ 0.614161718800186    2    2    6    6
+ -2.730885919022776e-15    2    2    7    2
+ 0.1383084739263875    2    2    7    3
+ 1.887395396669715e-15    2    2    7    4
+ 0.6242300746097845    2    2    7    7
+ 0.01099438448594219    3    1    3    1
+ 0.0177649938496416    3    1    3    2
+ 0.003439420521538433    3    1    4    3
+ 0.003149458301461384    3    1    6    3
+ 0.01531468707175935    3    1    7    1
+ 0.0138782501420901    3    1    7    2
+ 0.009634200496180479    3    1    7    4
+ 0.009208022839475303    3    1    7    6
+ 0.0177649938496416    3    2    3    1
+ 0.1444191819226373    3    2    3    2
+ -0.0199811196581604    3    2    4    3
+ -0.04010539607228727    3    2    6    3
+ 1.374212950574065e-15    3    2    6    4
+ 3.480715932451514e-15    3    2    6    6
+ 0.02310286336599912    3    2    7    1
+ 0.0403465306041421    3    2    7    2
+ 0.07710108477473231    3    2    7    4
+ 0.09859556273872483    3    2    7    6
+ -3.599676369507471e-15    3    2    7    7
+ 0.7999304193144107    3    3    1    1
+ -0.004406283490517992    3    3    2    1
+ 0.6451438166645459    3    3    2    2
+ 0.6329857220999804    3    3    3    3
+ -0.006469096930464782    3    3    4    1
+ -0.006903981904082305    3    3    4    2
+ 0.5984897705101471    3    3    4    4
+ 0.6288776966159403    3    3    5    5
+ 0.0002004859018046899    3    3    6    1
+ 0.07587274552656574    3    3    6    2
+ 1.512602221862055e-15    3    3    6    3
+ -0.04325010718698175    3    3    6    4
+ 0.5714432667033956    3    3    6    6
+ -1.657596147047677e-15    3    3    7    2
+ 0.09041020817349989    3    3    7    3
+ 1.156825186122905e-15    3    3    7    6
+ 0.6107416373778662    3    3    7    7
+ -0.1835967115185898    4    1    1    1
+ 0.02249791049099268    4    1    2    1
+ -0.01605176798440034    4    1    2    2
+ -0.006469096930464793    4    1    3    3
+ 0.02777284909384414    4    1    4    1
+ 0.0189209634151068    4    1    4    2
+ 0.01136200052559131    4    1    4    4
+ -0.005118262205227389    4    1    5    5
+ 0.0005620657839054863    4    1    6    1
+ -0.0186515813791622    4    1    6    2
+ 0.002338317445662476    4    1    6    4
+ -0.02118857118656031    4    1    6    6
+ 0.0008237209373926195    4    1    7    3
+ -0.004169127252127583    4    1    7    7
+ 0.1284623636944898    4    2    1    1
+ -0.009212365591395861    4    2    2    1
+ -0.004070151432527629    4    2    2    2
+ -0.006903981904082291    4    2    3    3
+ 0.0189209634151068    4    2    4    1
+ 0.1240556379933925    4    2    4    2
+ 0.1044464677045656    4    2    4    4
+ 0.06880867186643734    4    2    5    5
+ -0.02034480817937473    4    2    6    1
+ -0.02096799027742031    4    2    6    2
+ 1.399183521667672e-15    4    2    6    3
+ -0.03143005920448164    4    2    6    4
+ -0.05857904410649258    4    2    6    6
+ 0.07624760177013633    4    2    7    3
+ 1.48301483276511e-15    4    2    7    6
+ 0.01382975099652568    4    2    7    7
+ 0.003439420521538433    4    3    3    1
+ -0.01998111965816038    4    3    3    2
+ 0.04725737617877469    4    3    4    3
+ 0.0286139802412324    4    3    6    3
+ -1.620106278602005e-15    4    3    6    6
+ 0.004959092277355116    4    3    7    1
+ 0.03407934345368872    4    3    7    2
+ 0.002271311972147127    4    3    7    4
+ -0.04766697950247159    4    3    7    6
+ 1.673096790626979e-15    4    3    7    7
+ 0.9997707458540169    4    4    1    1
+ -0.01356454804023786    4    4    2    1
+ 0.6756423771903608    4    4    2    2
+ 0.5984897705101472    4    4    3    3
+ 0.01136200052559131    4    4    4    1
+ 0.1044464677045657    4    4    4    2
+ 0.7826369265089579    4    4    4    4
+ 0.7288776509205789    4    4    5    5
+ -0.01923646820554455    4    4    6    1
+ 0.08818613371425008    4    4    6    2
+ 3.368914361628211e-15    4    4    6    3
+ -0.1213863402802925    4    4    6    4
+ 0.5490733794041828    4    4    6    6
+ -1.532226487351395e-15    4    4    7    2
+ 0.1600228668704621    4    4    7    3
+ 2.373594622869556e-15    4    4    7    4
+ 0.6082119440506374    4    4    7    7
+ -1.09668093608046e-15    5    1    1    1
+ 0.02604454407916306    5    1    5    1
+ 0.0324620026971651    5    1    5    2
+ 0.01345017628891254    5    1    5    4
+ 0.0157466867492852    5    1    6    5
+ 0.03246200269716511    5    2    5    1
+ 0.1444653280412293    5    2    5    2
+ 0.04690886238904591    5    2    5    4
+ 1.060973968044138e-15    5    2    5    5
+ 0.059106273757118    5    2    6    5
+ -1.00284132890759e-15    5    2    7    5
+ 0.02880002005414797    5    3    5    3
+ 0.02368785864256855    5    3    7    5
+ -1.234611383062099e-15    5    4    1    1
+ 0.01345017628891254    5    4    5    1
+ 0.04690886238904591    5    4    5    2
+ 0.05591041615623339    5    4    5    4
+ 0.001740514519325905    5    4    6    5
+ 1.115336271961169    5    5    1    1
+ -0.01169439676993426    5    5    2    1
+ 0.747405683552989    5    5    2    2
+ 0.6288776966159402    5    5    3    3
+ -0.005118262205227399    5    5    4    1
+ 0.06880867186643731    5    5    4    2
+ 0.7288776509205787    5    5    4    4
+ 1.060973968044138e-15    5    5    5    2
+ 0.8801590933750444    5    5    5    5
+ -0.006208141382350879    5    5    6    1
+ 0.1585774903309637    5    5    6    2
+ 3.654959079105337e-15    5    5    6    3
+ -0.1163101710278052    5    5    6    4
+ 0.588948063505042    5    5    6    6
+ -2.890885345658369e-15    5    5    7    2
+ 0.1898227366397078    5    5    7    3
+ 2.372814355500244e-15    5    5    7    4
+ 0.6249850013591214    5    5    7    7
+ -0.2379642076436866    6    1    1    1
+ 0.0357951665411642    6    1    2    1
+ -0.0007824328446839209    6    1    2    2
+ 0.0002004859018046781    6    1    3    3
+ 0.000562065783905491    6    1    4    1
+ -0.02034480817937472    6    1    4    2
+ -0.01923646820554455    6    1    4    4
+ -0.00620814138235087    6    1    5    5
+ 0.03130931279981761    6    1    6    1
+ 0.006811265496433188    6    1    6    2
+ 0.0002855122873411405    6    1    6    4
+ 0.008408979484920942    6    1    6    6
+ -0.006986551345356871    6    1    7    3
+ -0.005126079207404938    6    1    7    7
+ 0.3082974147602018    6    2    1    1
+ -0.006646269777667536    6    2    2    1
+ 0.1429042097541408    6    2    2    2
+ 0.07587274552656581    6    2    3    3
+ -0.01865158137916219    6    2    4    1
+ -0.02096799027742028    6    2    4    2
+ 0.08818613371425017    6    2    4    4
+ 0.1585774903309637    6    2    5    5
+ 0.006811265496433182    6    2    6    1
+ 0.101890101151612    6    2    6    2
+ 2.150546348926948e-15    6    2    6    3
+ -0.06097304265413169    6    2    6    4
+ 0.09678768146184304    6    2    6    6
+ 0.07673587941703364    6    2    7    3
+ 1.546322779389e-15    6    2    7    4
+ 0.06904640998463688    6    2    7    7
+ 6.927221004821335e-15    6    3    1    1
+ 2.647862510188213e-15    6    3    2    2
+ 0.003149458301461385    6    3    3    1
+ -0.04010539607228727    6    3    3    2
+ 1.501623171712475e-15    6    3    3    3
+ 1.392288592051266e-15    6    3    4    2
+ 0.0286139802412324    6    3    4    3
+ 3.304881771788963e-15    6    3    4    4
+ 3.659688081047736e-15    6    3    5    5
+ 2.1713556148204e-15    6    3    6    2
+ 0.07092212645449124    6    3    6    3
+ -2.128917174079112e-15    6    3    6    4
+ -1.411650380322749e-15    6    3    6    6
+ 0.003863020755633329    6    3    7    1
+ 0.03533060229283891    6    3    7    2
+ 1.803316952931579e-15    6    3    7    3
+ -0.04445451849340859    6    3    7    4
+ -0.06451443741757142    6    3    7    6
+ 4.060671031660444e-15    6    3    7    7
+ -0.2194080179749165    6    4    1    1
+ 0.002234961476987666    6    4    2    1
+ -0.09546441068766688    6    4    2    2
+ 1.38831117260477e-15    6    4    3    2
+ -0.04325010718698176    6    4    3    3
+ 0.002338317445662473    6    4    4    1
+ -0.03143005920448162    6    4    4    2
+ -0.1213863402802925    6    4    4    4
+ -0.1163101710278052    6    4    5    5
+ 0.0002855122873411347    6    4    6    1
+ -0.06097304265413171    6    4    6    2
+ -2.160936453509478e-15    6    4    6    3
+ 0.06874080026687987    6    4    6    4
+ -0.04459738388695341    6    4    6    6
+ 1.483232651420531e-15    6    4    7    2
+ -0.07849156307841754    6    4    7    3
+ -0.04154608497977489    6    4    7    7
+ -1.555242604178036e-15    6    5    1    1
+ 0.01574668674928521    6    5    5    1
+ 0.05910627375711799    6    5    5    2
+ 0.001740514519325887    6    5    5    4
+ 0.03858991540578578    6    5    6    5
+ 0.8026929114041936    6    6    1    1
+ -0.006978989180611163    6    6    2    1
+ 0.6141617188001863    6    6    2    2
+ 3.512113923190604e-15    6    6    3    2
+ 0.5714432667033957    6    6    3    3
+ -0.02118857118656031    6    6    4    1
+ -0.05857904410649251    6    6    4    2
+ -1.631137921943511e-15    6    6    4    3
+ 0.549073379404183    6    6    4    4
+ 0.5889480635050423    6    6    5    5
+ 0.008408979484920894    6    6    6    1
+ 0.09678768146184298    6    6    6    2
+ -1.514935256001914e-15    6    6    6    3
+ -0.04459738388695345    6    6    6    4
+ 0.5971308786426219    6    6    6    6
+ -2.107418731265503e-15    6    6    7    2
+ 0.0379522686064513    6    6    7    3
+ 2.73940564220636e-15    6    6    7    4
+ 3.727210921586574e-15    6    6    7    6
+ 0.5663097768142775    6    6    7    7
+ 3.309414917719409e-15    7    1    1    1
+ 0.01531468707175935    7    1    3    1
+ 0.02310286336599912    7    1    3    2
+ 0.004959092277355118    7    1    4    3
+ 0.003863020755633331    7    1    6    3
+ 0.02139027044584247    7    1    7    1
+ 0.01830755315227746    7    1    7    2
+ 0.01319881251038175    7    1    7    4
+ 0.01219122749270908    7    1    7    6
+ -5.527433074026119e-15    7    2    1    1
+ -2.705752738841816e-15    7    2    2    2
+ 0.0138782501420901    7    2    3    1
+ 0.0403465306041421    7    2    3    2
+ -1.673498562697683e-15    7    2    3    3
+ 0.03407934345368874    7    2    4    3
+ -1.522089726447589e-15    7    2    4    4
+ -2.906553331821465e-15    7    2    5    5
+ 0.03533060229283891    7    2    6    3
+ 1.489865807033632e-15    7    2    6    4
+ -2.150838397427569e-15    7    2    6    6
+ 0.01830755315227745    7    2    7    1
+ 0.06191200289296931    7    2    7    2
+ -2.019554797677553e-15    7    2    7    3
+ 0.01667132443343576    7    2    7    4
+ -0.009945921505355977    7    2    7    6
+ 0.3624115212778843    7    3    1    1
+ -0.007503381453849525    7    3    2    1
+ 0.1383084739263876    7    3    2    2
+ 0.09041020817349986    7    3    3    3
+ 0.0008237209373926215    7    3    4    1
+ 0.07624760177013637    7    3    4    2
+ 0.1600228668704619    7    3    4    4
+ 0.1898227366397079    7    3    5    5
+ -0.006986551345356867    7    3    6    1
+ 0.07673587941703366    7    3    6    2
+ 1.771878334389302e-15    7    3    6    3
+ -0.07849156307841747    7    3    6    4
+ 0.03795226860645124    7    3    6    6
+ -2.014134055633541e-15    7    3    7    2
+ 0.1524934650232258    7    3    7    3
+ 2.380431702900645e-15    7    3    7    4
+ 2.115530786338072e-15    7    3    7    6
+ 0.09320616911421969    7    3    7    7
+ 4.604365715466489e-15    7    4    1    1
+ 1.850926226512467e-15    7    4    2    2
+ 0.009634200496180479    7    4    3    1
+ 0.07710108477473229    7    4    3    2
+ 0.002271311972147138    7    4    4    3
+ 2.372967331185579e-15    7    4    4    4
+ 2.369475964146786e-15    7    4    5    5
+ 1.525381698626244e-15    7    4    6    2
+ -0.04445451849340858    7    4    6    3
+ 2.662393406608827e-15    7    4    6    6
+ 0.01319881251038176    7    4    7    1
+ 0.01667132443343575    7    4    7    2
+ 2.379444600440423e-15    7    4    7    3
+ 0.06897687707695072    7    4    7    4
+ 0.05791860369249648    7    4    7    6
+ -1.182067037200129e-15    7    4    7    7
+ -1.002714729007794e-15    7    5    5    2
+ 0.02368785864256856    7    5    5    3
+ 0.02435072550962004    7    5    7    5
+ 1.575287290602142e-15    7    6    1    1
+ 0.009208022839475308    7    6    3    1
+ 0.09859556273872483    7    6    3    2
+ 1.163632817204229e-15    7    6    3    3
+ 1.528026251876103e-15    7    6    4    2
+ -0.04766697950247158    7    6    4    3
+ -0.06451443741757144    7    6    6    3
+ 3.451306051223984e-15    7    6    6    6
+ 0.01219122749270908    7    6    7    1
+ -0.009945921505355967    7    6    7    2
+ 2.108990168998312e-15    7    6    7    3
+ 0.05791860369249649    7    6    7    4
+ 0.1153064014297518    7    6    7    6
+ -2.90540271318117e-15    7    6    7    7
+ 0.8689352741666504    7    7    1    1
+ -0.009399378211486287    7    7    2    1
+ 0.6242300746097844    7    7    2    2
+ -3.598571497893575e-15    7    7    3    2
+ 0.6107416373778664    7    7    3    3
+ -0.004169127252127587    7    7    4    1
+ 0.01382975099652561    7    7    4    2
+ 1.64838906100722e-15    7    7    4    3
+ 0.6082119440506373    7    7    4    4
+ 0.6249850013591215    7    7    5    5
+ -0.005126079207404943    7    7    6    1
+ 0.06904640998463674    7    7    6    2
+ 3.897057995836122e-15    7    7    6    3
+ -0.04154608497977499    7    7    6    4
+ 0.5663097768142774    7    7    6    6
+ 0.09320616911421972    7    7    7    3
+ -1.215132998625822e-15    7    7    7    4
+ -2.85056499123113e-15    7    7    7    6
+ 0.6195154900510039    7    7    7    7
+ -32.70260435195033    1    1  0  0
+ 0.5581081719507953    2    1  0  0
+ -7.670748541529406    2    2  0  0
+ -1.018772020573078e-15    3    1  0  0
+ 2.834885411551239e-15    3    2  0  0
+ -6.363962604047698    3    3  0  0
+ 0.2351902329540697    4    1  0  0
+ -0.4316849771936422    4    2  0  0
+ 1.681441594146254e-15    4    3  0  0
+ -6.986223540237033    4    4  0  0
+ 1.407073905611408e-15    5    1  0  0
+ -3.377416934408135e-15    5    2  0  0
+ 1.018186100079006e-15    5    3  0  0
+ 6.518590162557252e-15    5    4  0  0
+ -7.45717009854408    5    5  0  0
+ 0.3046055386583093    6    1  0  0
+ -1.381405899634734    6    2  0  0
+ -3.308711234164146e-14    6    3  0  0
+ 1.080199730359876    6    4  0  0
+ 8.129235496684577e-15    6    5  0  0
+ -5.33601461036252    6    6  0  0
+ 2.492743299483893e-14    7    2  0  0
+ -1.709921420260744    7    3  0  0
+ -2.192612893034448e-14    7    4  0  0
+ 2.314841418581213e-15    7    5  0  0
+ -6.770194595807969e-15    7    6  0  0
+ -5.603486827980103    7    7  0  0
+ 9.189533762934902  0  0  0  0

wai/wai-quantum/tests/data/h6.fcidump

@@ -0,0 +1,255 @@
+ &FCI NORB=   6,NELEC= 6,MS2=0,
+  ORBSYM=1,1,1,1,1,1,
+  ISYM=1,
+ &END
+ 0.4295489168856813    1    1    1    1
+ 0.3468506155380016    1    1    2    2
+ -0.07974263467018042    1    1    3    1
+ 0.3643124474081217    1    1    3    3
+ -0.07982511288203843    1    1    4    2
+ 0.3707417649103147    1    1    4    4
+ 0.004536611863310048    1    1    5    1
+ 0.08294888117938162    1    1    5    3
+ 0.3658559707829486    1    1    5    5
+ 0.006079883228768165    1    1    6    2
+ 0.08273202586394117    1    1    6    4
+ 0.458681932122968    1    1    6    6
+ 0.1332007723024807    2    1    2    1
+ 0.1012040659075259    2    1    3    2
+ -0.05122561264726874    2    1    4    1
+ 0.08383365112831728    2    1    4    3
+ 0.04396669132011263    2    1    5    2
+ -0.1047396278519592    2    1    5    4
+ -0.001758103405853429    2    1    6    1
+ 0.05106706248722689    2    1    6    3
+ 0.1364871557940125    2    1    6    5
+ -1.440438479008169e-15    2    1    6    6
+ 0.3468506155380017    2    2    1    1
+ 0.3778345942116631    2    2    2    2
+ 0.02807821092099262    2    2    3    1
+ 0.3466585252364577    2    2    3    3
+ -0.01293997684894165    2    2    4    2
+ 0.3512668956623429    2    2    4    4
+ 0.03643623471367177    2    2    5    1
+ 0.01472231571247107    2    2    5    3
+ -1.196959198423997e-15    2    2    5    4
+ 0.3857483602782906    2    2    5    5
+ 0.036875400248939    2    2    6    2
+ -0.02071351880399119    2    2    6    4
+ 0.3719934864810962    2    2    6    6
+ -0.07974263467018042    3    1    1    1
+ 0.02807821092099261    3    1    2    2
+ 0.1027044849749038    3    1    3    1
+ -0.02107654467081059    3    1    3    3
+ 0.06059029244659133    3    1    4    2
+ -0.02177854546188477    3    1    4    4
+ 0.03338982673125539    3    1    5    1
+ -0.06310854812487629    3    1    5    3
+ 0.0167720419941487    3    1    5    5
+ 0.03153281477382371    3    1    6    2
+ -0.09893780526013846    3    1    6    4
+ -0.08570577708516933    3    1    6    6
+ 0.101204065907526    3    2    2    1
+ 0.1265054827462463    3    2    3    2
+ 0.01519358972525283    3    2    4    1
+ 0.08468268676133547    3    2    4    3
+ 0.001855911283668451    3    2    5    2
+ -0.1200881982509644    3    2    5    4
+ 0.02460147054842003    3    2    6    1
+ -0.008085382690766998    3    2    6    3
+ 0.1067353007850898    3    2    6    5
+ 0.3643124474081217    3    3    1    1
+ 0.3466585252364577    3    3    2    2
+ -0.02107654467081057    3    3    3    1
+ 0.3703455344211578    3    3    3    3
+ -0.002505967824627193    3    3    4    2
+ 0.3646857402777312    3    3    4    4
+ -0.01618227027332034    3    3    5    1
+ 0.01380931478346078    3    3    5    3
+ 0.3620114612301668    3    3    5    5
+ -0.008577808331913114    3    3    6    2
+ 0.02373758610096638    3    3    6    4
+ 0.3929579432040012    3    3    6    6
+ -0.05122561264726874    4    1    2    1
+ 0.01519358972525282    4    1    3    2
+ 0.07932363600022328    4    1    4    1
+ -0.009562021775725035    4    1    4    3
+ -0.05212217080098179    4    1    5    2
+ -0.004601388004221819    4    1    5    4
+ 0.02960125909361224    4    1    6    1
+ -0.07313258390431737    4    1    6    3
+ -0.05166886696266604    4    1    6    5
+ -0.07982511288203847    4    2    1    1
+ -0.01293997684894165    4    2    2    2
+ 0.06059029244659132    4    2    3    1
+ -0.002505967824627203    4    2    3    3
+ 0.0866200805109425    4    2    4    2
+ -0.01457653604125899    4    2    4    4
+ -0.02764284132912115    4    2    5    1
+ -0.08002059568981182    4    2    5    3
+ -0.01915173077818465    4    2    5    5
+ -0.02253604377940627    4    2    6    2
+ -0.06259483081312084    4    2    6    4
+ -0.08733550413212422    4    2    6    6
+ 0.08383365112831728    4    3    2    1
+ 0.08468268676133546    4    3    3    2
+ -0.009562021775725028    4    3    4    1
+ 0.1081289494308392    4    3    4    3
+ -0.03346747921824173    4    3    5    2
+ -0.0858941741526251    4    3    5    4
+ -0.03999894894037705    4    3    6    1
+ 0.01090459091616766    4    3    6    3
+ 0.08942410439787174    4    3    6    5
+ 0.3707417649103147    4    4    1    1
+ 0.3512668956623428    4    4    2    2
+ -0.02177854546188478    4    4    3    1
+ 0.3646857402777312    4    4    3    3
+ -0.01457653604125898    4    4    4    2
+ 0.3789890914489369    4    4    4    4
+ -0.006474192094668946    4    4    5    1
+ 0.01068861382162299    4    4    5    3
+ 0.3703942517777273    4    4    5    5
+ -0.01057032172682465    4    4    6    2
+ 0.02555283264156673    4    4    6    4
+ 0.403341686139579    4    4    6    6
+ 0.004536611863310058    5    1    1    1
+ 0.03643623471367177    5    1    2    2
+ 0.03338982673125539    5    1    3    1
+ -0.0161822702733203    5    1    3    3
+ -0.02764284132912114    5    1    4    2
+ -0.006474192094668924    5    1    4    4
+ 0.05549981522359869    5    1    5    1
+ 0.01992225081782516    5    1    5    3
+ 0.03431870994071447    5    1    5    5
+ 0.05008558336516987    5    1    6    2
+ -0.03075145906197635    5    1    6    4
+ 0.005202992851507899    5    1    6    6
+ 0.04396669132011262    5    2    2    1
+ 0.001855911283668446    5    2    3    2
+ -0.05212217080098179    5    2    4    1
+ -0.0334674792182417    5    2    4    3
+ 0.08556407022984075    5    2    5    2
+ -0.005747342818933686    5    2    5    4
+ 0.03390839675434936    5    2    6    1
+ 0.05157543156832715    5    2    6    3
+ 0.04570023568213209    5    2    6    5
+ 0.08294888117938162    5    3    1    1
+ 0.01472231571247106    5    3    2    2
+ -0.06310854812487629    5    3    3    1
+ 0.01380931478346081    5    3    3    3
+ -0.08002059568981182    5    3    4    2
+ 0.01068861382162304    5    3    4    4
+ 0.01992225081782516    5    3    5    1
+ 0.08623149487954507    5    3    5    3
+ 0.02093226922172396    5    3    5    5
+ 0.02449285467948081    5    3    6    2
+ 0.0649591803143143    5    3    6    4
+ 0.09329288198911408    5    3    6    6
+ -0.1047396278519591    5    4    2    1
+ -1.181682112187318e-15    5    4    2    2
+ -0.1200881982509644    5    4    3    2
+ -0.00460138800422183    5    4    4    1
+ -0.08589417415262512    5    4    4    3
+ -0.005747342818933691    5    4    5    2
+ 0.128982445796214    5    4    5    4
+ -0.02190984245838463    5    4    6    1
+ 0.008331619084783103    5    4    6    3
+ -0.1130168678254847    5    4    6    5
+ 1.140979621303909e-15    5    4    6    6
+ 0.3658559707829486    5    5    1    1
+ 0.3857483602782906    5    5    2    2
+ 0.0167720419941487    5    5    3    1
+ 0.3620114612301669    5    5    3    3
+ -0.01915173077818466    5    5    4    2
+ 0.3703942517777273    5    5    4    4
+ 0.03431870994071444    5    5    5    1
+ 0.02093226922172394    5    5    5    3
+ 0.4126507510401768    5    5    5    5
+ 0.03649148858437951    5    5    6    2
+ -0.01961392259484946    5    5    6    4
+ 0.4024127949974573    5    5    6    6
+ -0.001758103405853435    6    1    2    1
+ 0.02460147054842003    6    1    3    2
+ 0.02960125909361223    6    1    4    1
+ -0.03999894894037701    6    1    4    3
+ 0.03390839675434937    6    1    5    2
+ -0.02190984245838461    6    1    5    4
+ 0.06912553239638672    6    1    6    1
+ -0.02802006549712804    6    1    6    3
+ -0.002076148506383364    6    1    6    5
+ 0.006079883228768178    6    2    1    1
+ 0.036875400248939    6    2    2    2
+ 0.03153281477382371    6    2    3    1
+ -0.00857780833191309    6    2    3    3
+ -0.02253604377940627    6    2    4    2
+ -0.01057032172682461    6    2    4    4
+ 0.05008558336516987    6    2    5    1
+ 0.02449285467948082    6    2    5    3
+ 0.03649148858437955    6    2    5    5
+ 0.05243596890352582    6    2    6    2
+ -0.03137873874068051    6    2    6    4
+ 0.007486653259193091    6    2    6    6
+ 0.0510670624872269    6    3    2    1
+ -0.008085382690766977    6    3    3    2
+ -0.0731325839043174    6    3    4    1
+ 0.01090459091616769    6    3    4    3
+ 0.05157543156832717    6    3    5    2
+ 0.008331619084783087    6    3    5    4
+ -0.02802006549712806    6    3    6    1
+ 0.07772450921252712    6    3    6    3
+ 0.05618663454551731    6    3    6    5
+ -1.245591052033635e-15    6    3    6    6
+ 0.08273202586394117    6    4    1    1
+ -0.02071351880399119    6    4    2    2
+ -0.09893780526013846    6    4    3    1
+ 0.02373758610096634    6    4    3    3
+ -0.06259483081312084    6    4    4    2
+ 0.02555283264156676    6    4    4    4
+ -0.03075145906197634    6    4    5    1
+ 0.06495918031431429    6    4    5    3
+ -0.01961392259484944    6    4    5    5
+ -0.03137873874068051    6    4    6    2
+ 0.1080434259706099    6    4    6    4
+ 0.09526081277536884    6    4    6    6
+ 0.1364871557940125    6    5    2    1
+ 0.1067353007850899    6    5    3    2
+ -0.05166886696266609    6    5    4    1
+ 0.0894241043978718    6    5    4    3
+ 0.04570023568213216    6    5    5    2
+ -0.1130168678254847    6    5    5    4
+ -0.002076148506383372    6    5    6    1
+ 0.05618663454551732    6    5    6    3
+ 0.1546561696289368    6    5    6    5
+ -1.308958149825214e-15    6    5    6    6
+ 0.458681932122968    6    6    1    1
+ -1.39341663207837e-15    6    6    2    1
+ 0.3719934864810962    6    6    2    2
+ -0.08570577708516929    6    6    3    1
+ 0.3929579432040012    6    6    3    3
+ -0.08733550413212425    6    6    4    2
+ 0.4033416861395791    6    6    4    4
+ 0.005202992851507952    6    6    5    1
+ 0.09329288198911405    6    6    5    3
+ 1.068589661201713e-15    6    6    5    4
+ 0.4024127949974575    6    6    5    5
+ 0.007486653259193086    6    6    6    2
+ -1.200428645375951e-15    6    6    6    3
+ 0.09526081277536878    6    6    6    4
+ -1.342299074187513e-15    6    6    6    5
+ 0.5177055407346068    6    6    6    6
+ -2.281751933013601    1    1  0  0
+ -2.040945268831873    2    2  0  0
+ 0.1458668260785023    3    1  0  0
+ -1.88908673174047    3    3  0  0
+ 0.2110592057665517    4    2  0  0
+ -1.998401444325282e-15    4    3  0  0
+ -1.675701879970149    4    4  0  0
+ -0.06418639468496004    5    1  0  0
+ -0.1739059788109305    5    3  0  0
+ -1.38367990976047    5    5  0  0
+ -0.0417230422669391    6    2  0  0
+ 1.644517855226013e-15    6    3  0  0
+ -0.1535423833278287    6    4  0  0
+ 1.006251619321429e-15    6    5  0  0
+ -1.209826611879526    6    6  0  0
+ 4.603841735004002  0  0  0  0

wai/wai-quantum/tests/data/lih.fcidump

@@ -0,0 +1,194 @@
+ &FCI NORB=   6,NELEC= 4,MS2=0,
+  ORBSYM=1,1,1,1,1,1,
+  ISYM=1,
+ &END
+ 1.658551513231926    1    1    1    1
+ -0.1119411348240613    1    1    2    1
+ 0.3673101014750008    1    1    2    2
+ -0.1385319239318432    1    1    3    1
+ 0.01334611115799815    1    1    3    2
+ 0.3956539136685808    1    1    3    3
+ 0.3963189291588849    1    1    4    4
+ 0.3963189291588851    1    1    5    5
+ 0.05263815043175974    1    1    6    1
+ -0.04091412316813224    1    1    6    2
+ 0.01764596781263263    1    1    6    3
+ 0.3617428178975034    1    1    6    6
+ -0.1119411348240612    2    1    1    1
+ 0.01339683430673438    2    1    2    1
+ 0.006258341962824109    2    1    2    2
+ 0.01123036367889496    2    1    3    1
+ -0.003363202997395241    2    1    3    2
+ -0.01106470491241514    2    1    3    3
+ -0.004366866519105297    2    1    4    4
+ -0.0043668665191053    2    1    5    5
+ -0.008878377322386331    2    1    6    1
+ 0.004741251801471238    2    1    6    2
+ -0.00369300693476956    2    1    6    3
+ 0.003316786999488483    2    1    6    6
+ 0.3673101014750008    2    2    1    1
+ 0.006258341962824105    2    2    2    1
+ 0.4876578261514061    2    2    2    2
+ -0.01592569304822697    2    2    3    1
+ -0.04849495258111529    2    2    3    2
+ 0.2237530773670962    2    2    3    3
+ 0.2704181616357398    2    2    4    4
+ 0.27041816163574    2    2    5    5
+ -0.006804882104791949    2    2    6    1
+ 0.1270522788593902    2    2    6    2
+ -0.05134078354716973    2    2    6    3
+ 0.4540419572677863    2    2    6    6
+ -0.1385319239318431    3    1    1    1
+ 0.01123036367889495    3    1    2    1
+ -0.01592569304822697    3    1    2    2
+ 0.02165565463641273    3    1    3    1
+ 0.00017922897593376    3    1    3    2
+ 0.001833245648913372    3    1    3    3
+ -0.004973743953774324    3    1    4    4
+ -0.004973743953774327    3    1    5    5
+ -0.002308669573223772    3    1    6    1
+ 0.0005015830060559282    3    1    6    2
+ 0.004400888432867765    3    1    6    3
+ -0.01133739604976779    3    1    6    6
+ 0.01334611115799818    3    2    1    1
+ -0.00336320299739524    3    2    2    1
+ -0.04849495258111526    3    2    2    2
+ 0.0001792289759337518    3    2    3    1
+ 0.013013969421508    3    2    3    2
+ 0.007418200743174592    3    2    3    3
+ 0.005712903933711085    3    2    4    4
+ 0.005712903933711089    3    2    5    5
+ 0.001669872514413443    3    2    6    1
+ -0.03454100054795995    3    2    6    2
+ 0.009357453432246147    3    2    6    3
+ -0.04329410713641037    3    2    6    6
+ 0.3956539136685808    3    3    1    1
+ -0.01106470491241513    3    3    2    1
+ 0.2237530773670962    3    3    2    2
+ 0.001833245648913375    3    3    3    1
+ 0.007418200743174591    3    3    3    2
+ 0.3379349975894662    3    3    3    3
+ 0.2820037734839136    3    3    4    4
+ 0.2820037734839137    3    3    5    5
+ 0.01040808711633322    3    3    6    1
+ -0.01228417210130269    3    3    6    2
+ 0.03598190448284638    3    3    6    3
+ 0.2414678300084735    3    3    6    6
+ 0.00981794099440618    4    1    4    1
+ 0.007492522187157329    4    1    4    2
+ 0.01025687809344547    4    1    4    3
+ -0.006108198728404846    4    1    6    4
+ 0.00749252218715733    4    2    4    1
+ 0.02345011774315864    4    2    4    2
+ 0.01927261285133701    4    2    4    3
+ -0.01957479577021053    4    2    6    4
+ 0.01025687809344547    4    3    4    1
+ 0.01927261285133701    4    3    4    2
+ 0.04127779388220811    4    3    4    3
+ -0.01373211854998539    4    3    6    4
+ 0.3963189291588849    4    4    1    1
+ -0.004366866519105309    4    4    2    1
+ 0.2704181616357398    4    4    2    2
+ -0.004973743953774335    4    4    3    1
+ 0.005712903933711086    4    4    3    2
+ 0.2820037734839135    4    4    3    3
+ 0.3129455111594089    4    4    4    4
+ 0.279207232132027    4    4    5    5
+ 0.0005730651471392125    4    4    6    1
+ -0.01603690066007866    4    4    6    2
+ 0.002194540624147631    4    4    6    3
+ 0.2681942478263322    4    4    6    6
+ 0.009817940994406185    5    1    5    1
+ 0.007492522187157334    5    1    5    2
+ 0.01025687809344548    5    1    5    3
+ -0.006108198728404849    5    1    6    5
+ 0.007492522187157334    5    2    5    1
+ 0.02345011774315866    5    2    5    2
+ 0.01927261285133702    5    2    5    3
+ -0.01957479577021054    5    2    6    5
+ 0.01025687809344548    5    3    5    1
+ 0.01927261285133702    5    3    5    2
+ 0.04127779388220813    5    3    5    3
+ -0.01373211854998539    5    3    6    5
+ 0.01686913951369098    5    4    5    4
+ 0.3963189291588851    5    5    1    1
+ -0.004366866519105295    5    5    2    1
+ 0.27041816163574    5    5    2    2
+ -0.004973743953774317    5    5    3    1
+ 0.005712903933711086    5    5    3    2
+ 0.2820037734839137    5    5    3    3
+ 0.279207232132027    5    5    4    4
+ 0.3129455111594093    5    5    5    5
+ 0.0005730651471391958    5    5    6    1
+ -0.01603690066007867    5    5    6    2
+ 0.002194540624147631    5    5    6    3
+ 0.2681942478263324    5    5    6    6
+ 0.05263815043175962    6    1    1    1
+ -0.008878377322386329    6    1    2    1
+ -0.006804882104791959    6    1    2    2
+ -0.002308669573223771    6    1    3    1
+ 0.001669872514413437    6    1    3    2
+ 0.01040808711633322    6    1    3    3
+ 0.0005730651471391972    6    1    4    4
+ 0.0005730651471391975    6    1    5    5
+ 0.008491729737162623    6    1    6    1
+ 0.0001275745284848168    6    1    6    2
+ 0.004302208257139307    6    1    6    3
+ -0.00302804462426087    6    1    6    6
+ -0.04091412316813216    6    2    1    1
+ 0.004741251801471238    6    2    2    1
+ 0.1270522788593901    6    2    2    2
+ 0.0005015830060559263    6    2    3    1
+ -0.03454100054795993    6    2    3    2
+ -0.01228417210130271    6    2    3    3
+ -0.01603690066007867    6    2    4    4
+ -0.01603690066007868    6    2    5    5
+ 0.0001275745284848165    6    2    6    1
+ 0.1238723253808021    6    2    6    2
+ -0.03185703217262079    6    2    6    3
+ 0.1345287312734732    6    2    6    6
+ 0.01764596781263265    6    3    1    1
+ -0.003693006934769559    6    3    2    1
+ -0.05134078354716971    6    3    2    2
+ 0.004400888432867778    6    3    3    1
+ 0.009357453432246162    6    3    3    2
+ 0.03598190448284642    6    3    3    3
+ 0.002194540624147617    6    3    4    4
+ 0.002194540624147618    6    3    5    5
+ 0.004302208257139305    6    3    6    1
+ -0.03185703217262077    6    3    6    2
+ 0.02643669328394296    6    3    6    3
+ -0.04405224647857216    6    3    6    6
+ -0.006108198728404846    6    4    4    1
+ -0.01957479577021053    6    4    4    2
+ -0.01373211854998538    6    4    4    3
+ 0.01971345998579963    6    4    6    4
+ -0.00610819872840485    6    5    5    1
+ -0.01957479577021054    6    5    5    2
+ -0.01373211854998539    6    5    5    3
+ 0.01971345998579964    6    5    6    5
+ 0.3617428178975034    6    6    1    1
+ 0.003316786999488494    6    6    2    1
+ 0.4540419572677863    6    6    2    2
+ -0.01133739604976777    6    6    3    1
+ -0.04329410713641037    6    6    3    2
+ 0.2414678300084735    6    6    3    3
+ 0.2681942478263322    6    6    4    4
+ 0.2681942478263323    6    6    5    5
+ -0.003028044624260894    6    6    6    1
+ 0.1345287312734731    6    6    6    2
+ -0.04405224647857225    6    6    6    3
+ 0.4539585121684841    6    6    6    6
+ -4.728421350668932    1    1  0  0
+ 0.1056827931754777    2    1  0  0
+ -1.494577436776578    2    2  0  0
+ 0.1670201071483866    3    1  0  0
+ 0.03303310622847724    3    2  0  0
+ -1.125883391287575    3    3  0  0
+ -1.136267633478075    4    4  0  0
+ -1.136267633478076    5    5  0  0
+ -0.03428713441312929    6    1  0  0
+ -0.05410240910957359    6    2  0  0
+ 0.03053995599123926    6    3  0  0
+ -0.9501040287604586    6    6  0  0
+ 0.9953176380940441  0  0  0  0

wai/wai-rs/Cargo.lock

@@ -3900,7 +3900,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.35"
+version = "0.3.36"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-rs/Cargo.toml

@@ -368,7 +368,7 @@ quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"]
 # qLDPC memory at circuit level: bivariate-bicycle codes, BP-OSD.
 quantum_ldpc = ["quantum_frame", "wai-quantum/quantum_ldpc"]
 # Quantum-inspired optimisation: simulated bifurcation, annealing, exact referee.
-quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg", "wai-quantum/quantum_tdvp"]
+quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg", "wai-quantum/quantum_tdvp", "wai-quantum/quantum_chem", "wai-quantum/quantum_nqs"]
 quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"]
 # Variational quantum machine learning (extensions/quantum-sim § VML,
 # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot

wai/wai-rs/src/bin/wai_quantum.rs

@@ -170,6 +170,11 @@ fn usage() -> ! {
          [--cutoff C] [--exact] [--json]\n  \
          wai-quantum tdvp [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--start plus|neel|up]\n    \
          [--bond D] [--dt DT] [--t T] [--every K] [--cutoff C] [--json]   (time evolution by TDVP)\n\n\
+         MOLECULES (DMRG from FCIDUMP integrals, full CI as referee):\n  \
+         wai-quantum chem --fcidump FILE [--bond D] [--sweeps S] [--penalty P] [--no-fci] [--json]\n\n\
+         NEURAL QUANTUM STATES (RBM, variational Monte Carlo, stochastic reconfiguration):\n  \
+         wai-quantum nqs [--model heisenberg|ising] [--lx L] [--ly L] [--h H] [--alpha A] [--steps N]\n    \
+         [--samples S] [--lr R] [--exact] [--seed S] [--json]\n\n\
          OPTIMISATION (simulated bifurcation, annealing, exact):\n  \
          wai-quantum qio (--gset FILE | --sk N [--instance S]) [--method dsb|bsb|hdsb|hbsb|sa|exact]\n    \
          [--steps N] [--trials T] [--dt D] [--c1 C] [--gamma G] [--beta0 B --beta1 B] [--descend]\n    \
@@ -1116,6 +1121,100 @@ fn main() {
             }
         }
 
+        // A molecule's ground state by DMRG from its integrals, against
+        // Hartree–Fock and full configuration interaction where it fits.
+        "chem" => {
+            use wai::quantum_chem::*;
+            let path = opt(args, "--fcidump", "").unwrap_or_else(|| die("--fcidump FILE is required"));
+            let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}")));
+            let fd = read_fcidump(&text).unwrap_or_else(|e| die(&format!("{path}: {e}")));
+            let cfg = ChemConfig { max_bond: num(args, "--bond", "", 64), sweeps: num(args, "--sweeps", "", 4), penalty: num(args, "--penalty", "", 1.0), ..ChemConfig::default() };
+            let mpo_bond = molecular_mpo(&fd, &cfg).max_bond();
+            let hf = hartree_fock_energy(&fd);
+            let r = ground_state(&fd, &cfg);
+            let full = if flag(args, "--no-fci") { None } else { fci(&fd) };
+            if flag(args, "--json") {
+                let sw: Vec<String> = r.sweep_energies.iter().map(|e| e.to_string()).collect();
+                println!(
+                    "{{\"bond\":{},\"discarded\":{},\"electrons\":{},\"energy\":{},\"fci\":{},\"fci_determinants\":{},\"hartree_fock\":{hf},\"mpo_bond\":{mpo_bond},\"orbitals\":{},\"sweeps\":[{}]}}",
+                    cfg.max_bond,
+                    r.discarded,
+                    fd.nelec,
+                    r.energy,
+                    full.as_ref().map_or("null".into(), |f| f.energy.to_string()),
+                    full.as_ref().map_or("null".into(), |f| f.dim.to_string()),
+                    fd.norb,
+                    sw.join(",")
+                );
+            } else {
+                println!("{path}: {} orbitals, {} electrons, MPO bond {mpo_bond}, DMRG bond ≤ {}", fd.norb, fd.nelec, cfg.max_bond);
+                println!("Hartree–Fock  {hf:.10}");
+                for (i, e) in r.sweep_energies.iter().enumerate() {
+                    println!("sweep {:>2}      {e:.10}", i + 1);
+                }
+                println!("discarded     {:.1e}", r.discarded);
+                match &full {
+                    Some(f) => println!("full CI       {:.10}   ({} determinants; DMRG − full CI {:.2e})", f.energy, f.dim, r.energy - f.energy),
+                    None => println!("full CI       not run"),
+                }
+            }
+        }
+
+        // A neural quantum state for a spin model on a chain or an open
+        // square lattice, against exact diagonalisation where it fits.
+        "nqs" => {
+            use wai::quantum_nqs::*;
+            let lx: usize = num(args, "--lx", "", 10);
+            let ly: usize = num(args, "--ly", "", 1);
+            let n = lx * ly;
+            if !(2..=400).contains(&n) {
+                die("--lx × --ly: 2 to 400 spins");
+            }
+            let kind = opt(args, "--model", "").unwrap_or_else(|| "heisenberg".into());
+            let edges = if ly == 1 { chain_edges(n, false) } else { square_edges(lx, ly, false) };
+            let heis = match kind.as_str() {
+                "heisenberg" => true,
+                "ising" => false,
+                _ => die("--model: heisenberg or ising"),
+            };
+            let model = if heis { SpinModel::heisenberg(n, &edges, 1.0) } else { SpinModel::ising(n, &edges, 1.0, num(args, "--h", "", 1.0)) };
+            let sector = if heis { Some(0) } else { None };
+            let exact = flag(args, "--exact");
+            if exact && n > 20 {
+                die("--exact enumerates every configuration: at most 20 spins");
+            }
+            let mut rbm = Rbm::new(n, num(args, "--alpha", "", 2), 0.01, num(args, "--seed", "", 7));
+            if heis {
+                rbm = rbm.with_marshall_sign(checkerboard(lx, ly));
+            }
+            let cfg = VmcConfig { samples: num(args, "--samples", "", 1000), exact, magnetisation: sector, lr: num(args, "--lr", "", 0.02), seed: num(args, "--seed", "", 7), ..VmcConfig::default() };
+            let steps: usize = num(args, "--steps", "", 300);
+            let reference = if n <= 20 { model.exact_ground_energy(sector) } else { None };
+            let params = rbm.params();
+            let mut vmc = Vmc::new(model, rbm, cfg);
+            let mut last = vmc.estimate();
+            for _ in 0..steps {
+                last = vmc.step();
+            }
+            let fin = vmc.estimate();
+            if flag(args, "--json") {
+                println!(
+                    "{{\"energy\":{},\"error\":{},\"exact\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}",
+                    fin.energy,
+                    fin.error,
+                    reference.map_or("null".into(), |e| e.to_string()),
+                    last.energy,
+                    fin.variance
+                );
+            } else {
+                println!("{kind}, {lx}×{ly}, RBM with {params} parameters, {steps} steps{}", if exact { ", exact expectations" } else { "" });
+                println!("energy        {:.8} ± {:.8}   variance {:.2e}", fin.energy, fin.error, fin.variance);
+                if let Some(e) = reference {
+                    println!("exact         {e:.8}   relative error {:.2e}", (fin.energy - e) / e.abs());
+                }
+            }
+        }
+
         // Quantum-inspired optimisation: Max-Cut (G-set text) or a spin
         // glass, by simulated bifurcation, annealing, or exact enumeration.
         // Max-Cut reports cuts; a spin glass reports energies.

wai/wai-rs/src/lib.rs

@@ -666,6 +666,12 @@ pub use wai_quantum::quantum_dmrg;
 #[cfg(feature = "quantum_qio")]
 pub use wai_quantum::quantum_tdvp;
 
+#[cfg(feature = "quantum_qio")]
+pub use wai_quantum::quantum_chem;
+
+#[cfg(feature = "quantum_qio")]
+pub use wai_quantum::quantum_nqs;
+
 #[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))]
 pub use wai_quantum::gates;
 #[cfg(feature = "quantum_qir")]

wai/wai-web/Cargo.lock

@@ -734,7 +734,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.35"
+version = "0.3.36"
 dependencies = [
  "blake3",
  "ed25519-dalek",