From e78d453c0bb653605b5f5523697635036adf2713 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 5 Oct 2026 11:35:54 -0400 Subject: [PATCH 1/3] wai-quantum: MPOs with a tensor per site, held sparse Each site carries its own operators between MPO bonds of their own size, stored as its non-zero entries in the order the dense loops visited them, so the uniform chains keep their bits (every pin unchanged). Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/src/quantum_dmrg.rs | 338 +++++++++++++--------------- wai/wai-quantum/src/quantum_tdvp.rs | 26 +-- 2 files changed, 175 insertions(+), 189 deletions(-) diff --git a/wai/wai-quantum/src/quantum_dmrg.rs b/wai/wai-quantum/src/quantum_dmrg.rs index 451c2c2c6..81f49e9f6 100644 --- a/wai/wai-quantum/src/quantum_dmrg.rs +++ b/wai/wai-quantum/src/quantum_dmrg.rs @@ -563,18 +563,34 @@ impl Chain { site[t + 1] = a.scaled(c); site[(t + 1) * w + (w - 1)] = b; } - Mpo { n: self.n, w, site } + let entries: Vec<(usize, usize, Op)> = (0..w * w).filter(|&k| site[k].0.iter().any(|&v| v != 0.0)).map(|k| (k / w, k % w, site[k])).collect(); + let mut left = vec![0.0; w]; + left[0] = 1.0; + let mut right = vec![0.0; w]; + right[w - 1] = 1.0; + Mpo { n: self.n, dims: vec![w; self.n + 1], sites: vec![entries; self.n], left, right } } } -/// A uniform matrix-product operator: `site[a·w + b]` is the one-site -/// operator from state `a` to state `b`; the chain starts in state 0 and -/// ends in state `w − 1`. +/// A matrix-product operator: site `i` maps MPO state `a < dims[i]` to state +/// `b < dims[i + 1]` through the one-site operator `op`, for each non-zero +/// `(a, b, op)` in `sites[i]` (sorted by `(a, b)`). The chain starts in the +/// combination `left` of states and ends in `right`. #[derive(Clone, Debug, PartialEq)] pub struct Mpo { pub n: usize, - pub w: usize, - site: Vec, + /// `dims[i]` is the MPO bond to the left of site `i`; there are `n + 1`. + pub dims: Vec, + pub(crate) sites: Vec>, + pub(crate) left: Vec, + pub(crate) right: Vec, +} + +impl Mpo { + /// The largest MPO bond. + pub fn max_bond(&self) -> usize { + self.dims.iter().copied().max().unwrap_or(1) + } } // --------------------------------------------------------------------------- @@ -892,54 +908,52 @@ pub struct DmrgResult { } /// Contract environments and effective Hamiltonians. Environments are -/// `L[a', w, a]` and `R[b', w, b]`, bra index first. +/// `L[a', w, a]` and `R[b', w, b]`, bra index first; `w` runs over the MPO +/// bond at that cut. pub(crate) struct Engine<'a> { mpo: &'a Mpo, - w: usize, } impl Engine<'_> { pub(crate) fn new(mpo: &Mpo) -> Engine<'_> { - Engine { mpo, w: mpo.w } + Engine { mpo } } - fn op(&self, a: usize, b: usize) -> &Op { - &self.mpo.site[a * self.w + b] + /// The MPO bond to the left of site `i`. + fn w(&self, i: usize) -> usize { + self.mpo.dims[i] } - fn nonzero(&self, a: usize, b: usize) -> bool { - self.op(a, b).0.iter().any(|&v| v != 0.0) + /// Site `i`'s non-zero operators, sorted by `(a, b)`. + fn entries(&self, i: usize) -> &[(usize, usize, Op)] { + &self.mpo.sites[i] } /// The environment at the left edge. pub(crate) fn left_edge(&self) -> Vec { - let mut l = vec![T::ZERO; self.w]; - l[0] = T::ONE; - l + self.mpo.left.iter().map(|&x| T::ONE.scale(x)).collect() } /// The environment at the right edge. pub(crate) fn right_edge(&self) -> Vec { - let mut r = vec![T::ZERO; self.w]; - r[self.w - 1] = T::ONE; - r + self.mpo.right.iter().map(|&x| T::ONE.scale(x)).collect() } - /// `L'[b', w2, b] = Σ A*[a',s',b'] L[a',w1,a] W[w1,w2,s',s] A[a,s,b]`. - pub(crate) fn grow_left(&self, l: &[T], a: &[T], dl: usize, dr: usize) -> Vec { - let w = self.w; + /// Through site `i`: `L'[b', w2, b] = Σ A*[a',s',b'] L[a',w1,a] W[w1,w2,s',s] A[a,s,b]`. + pub(crate) fn grow_left(&self, i: usize, l: &[T], a: &[T], dl: usize, dr: usize) -> Vec { + let (wl, wr) = (self.w(i), self.w(i + 1)); // T[a', w1, s, b] = Σ_a L[a', w1, a] A[a, s, b] - let mut t = vec![T::ZERO; dl * w * 2 * dr]; + let mut t = vec![T::ZERO; dl * wl * 2 * dr]; for ap in 0..dl { - for w1 in 0..w { - let lrow = &l[(ap * w + w1) * dl..(ap * w + w1 + 1) * dl]; + for w1 in 0..wl { + let lrow = &l[(ap * wl + w1) * dl..(ap * wl + w1 + 1) * dl]; for (aa, &lv) in lrow.iter().enumerate() { if lv.is_zero() { continue; } for s in 0..2 { let src = &a[(aa * 2 + s) * dr..(aa * 2 + s + 1) * dr]; - let dst = &mut t[((ap * w + w1) * 2 + s) * dr..((ap * w + w1) * 2 + s + 1) * dr]; + let dst = &mut t[((ap * wl + w1) * 2 + s) * dr..((ap * wl + w1) * 2 + s + 1) * dr]; for (d, &x) in dst.iter_mut().zip(src) { *d = d.add(lv.mul(x)); } @@ -948,32 +962,26 @@ impl Engine<'_> { } } // U[a', w2, s', b] = Σ_{w1, s} W[w1, w2, s', s] T[a', w1, s, b] - let mut u = vec![T::ZERO; dl * w * 2 * dr]; + let mut u = vec![T::ZERO; dl * wr * 2 * dr]; for ap in 0..dl { - for w1 in 0..w { - for w2 in 0..w { - if !self.nonzero(w1, w2) { - continue; - } - let o = self.op(w1, w2).0; - for sp in 0..2 { - for s in 0..2 { - let m = o[sp * 2 + s]; - if m == 0.0 { - continue; - } - let src = ((ap * w + w1) * 2 + s) * dr; - let dst = ((ap * w + w2) * 2 + sp) * dr; - for b in 0..dr { - u[dst + b] = u[dst + b].add(t[src + b].scale(m)); - } + for &(w1, w2, o) in self.entries(i) { + for sp in 0..2 { + for s in 0..2 { + let m = o.0[sp * 2 + s]; + if m == 0.0 { + continue; + } + let src = ((ap * wl + w1) * 2 + s) * dr; + let dst = ((ap * wr + w2) * 2 + sp) * dr; + for b in 0..dr { + u[dst + b] = u[dst + b].add(t[src + b].scale(m)); } } } } } // L'[b', w2, b] = Σ_{a', s'} A*[a', s', b'] U[a', w2, s', b] - let mut out = vec![T::ZERO; dr * w * dr]; + let mut out = vec![T::ZERO; dr * wr * dr]; for ap in 0..dl { for sp in 0..2 { let arow = &a[(ap * 2 + sp) * dr..(ap * 2 + sp + 1) * dr]; @@ -982,9 +990,9 @@ impl Engine<'_> { continue; } let av = av.conj(); - for w2 in 0..w { - let src = &u[((ap * w + w2) * 2 + sp) * dr..((ap * w + w2) * 2 + sp + 1) * dr]; - let dst = &mut out[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; + for w2 in 0..wr { + let src = &u[((ap * wr + w2) * 2 + sp) * dr..((ap * wr + w2) * 2 + sp + 1) * dr]; + let dst = &mut out[(bp * wr + w2) * dr..(bp * wr + w2 + 1) * dr]; for (d, &x) in dst.iter_mut().zip(src) { *d = d.add(av.mul(x)); } @@ -995,52 +1003,46 @@ impl Engine<'_> { out } - /// `R'[a', w1, a] = Σ B*[a',s',b'] W[w1,w2,s',s] R[b',w2,b] B[a,s,b]`. - pub(crate) fn grow_right(&self, r: &[T], bt: &[T], dl: usize, dr: usize) -> Vec { - let w = self.w; + /// Through site `i`: `R'[a', w1, a] = Σ B*[a',s',b'] W[w1,w2,s',s] R[b',w2,b] B[a,s,b]`. + pub(crate) fn grow_right(&self, i: usize, r: &[T], bt: &[T], dl: usize, dr: usize) -> Vec { + let (wl, wr) = (self.w(i), self.w(i + 1)); // T[a, s, b', w2] = Σ_b B[a, s, b] R[b', w2, b] - let mut t = vec![T::ZERO; dl * 2 * dr * w]; + let mut t = vec![T::ZERO; dl * 2 * dr * wr]; for a in 0..dl { for s in 0..2 { let brow = &bt[(a * 2 + s) * dr..(a * 2 + s + 1) * dr]; for bp in 0..dr { - for w2 in 0..w { - let rrow = &r[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; + for w2 in 0..wr { + let rrow = &r[(bp * wr + w2) * dr..(bp * wr + w2 + 1) * dr]; let mut acc = T::ZERO; for (x, y) in brow.iter().zip(rrow) { acc = acc.add(x.mul(*y)); } - t[((a * 2 + s) * dr + bp) * w + w2] = acc; + t[((a * 2 + s) * dr + bp) * wr + w2] = acc; } } } } // U[a, s', b', w1] = Σ_{w2, s} W[w1, w2, s', s] T[a, s, b', w2] - let mut u = vec![T::ZERO; dl * 2 * dr * w]; - for w1 in 0..w { - for w2 in 0..w { - if !self.nonzero(w1, w2) { - continue; - } - let o = self.op(w1, w2).0; - for sp in 0..2 { - for s in 0..2 { - let m = o[sp * 2 + s]; - if m == 0.0 { - continue; - } - for a in 0..dl { - for bp in 0..dr { - let at = ((a * 2 + sp) * dr + bp) * w + w1; - u[at] = u[at].add(t[((a * 2 + s) * dr + bp) * w + w2].scale(m)); - } + let mut u = vec![T::ZERO; dl * 2 * dr * wl]; + for &(w1, w2, o) in self.entries(i) { + for sp in 0..2 { + for s in 0..2 { + let m = o.0[sp * 2 + s]; + if m == 0.0 { + continue; + } + for a in 0..dl { + for bp in 0..dr { + let at = ((a * 2 + sp) * dr + bp) * wl + w1; + u[at] = u[at].add(t[((a * 2 + s) * dr + bp) * wr + w2].scale(m)); } } } } } // R'[a', w1, a] = Σ_{s', b'} B*[a', s', b'] U[a, s', b', w1] - let mut out = vec![T::ZERO; dl * w * dl]; + let mut out = vec![T::ZERO; dl * wl * dl]; for ap in 0..dl { for sp in 0..2 { for bp in 0..dr { @@ -1050,9 +1052,9 @@ impl Engine<'_> { } let bv = bv.conj(); for a in 0..dl { - for w1 in 0..w { - let at = (ap * w + w1) * dl + a; - out[at] = out[at].add(bv.mul(u[((a * 2 + sp) * dr + bp) * w + w1])); + for w1 in 0..wl { + let at = (ap * wl + w1) * dl + a; + out[at] = out[at].add(bv.mul(u[((a * 2 + sp) * dr + bp) * wl + w1])); } } } @@ -1062,19 +1064,20 @@ impl Engine<'_> { } /// Rows `rows` (bra indices `a'`) of `H_eff θ` for the two-site tensor - /// `θ[a, s1, s2, b]`, written to `out` (which holds just those rows). + /// `θ[a, s1, s2, b]` on sites `i, i + 1`, written to `out` (which holds + /// just those rows). #[allow(clippy::too_many_arguments)] - fn apply_rows(&self, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { - let w = self.w; + fn apply_rows(&self, i: usize, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { + let (w0n, w1n, w2n) = (self.w(i), self.w(i + 1), self.w(i + 2)); let first = rows.start; let nrows = rows.len(); // T1[a', w0, s1, s2, b] = Σ_a L[a', w0, a] θ[a, s1, s2, b] let blk = 4 * dr; - let mut t1 = vec![T::ZERO; nrows * w * blk]; + let mut t1 = vec![T::ZERO; nrows * w0n * blk]; for ap in 0..nrows { - for w0 in 0..w { - let lrow = &l[((first + ap) * w + w0) * dl..((first + ap) * w + w0 + 1) * dl]; - let dst = &mut t1[(ap * w + w0) * blk..(ap * w + w0 + 1) * blk]; + for w0 in 0..w0n { + let lrow = &l[((first + ap) * w0n + w0) * dl..((first + ap) * w0n + w0 + 1) * dl]; + let dst = &mut t1[(ap * w0n + w0) * blk..(ap * w0n + w0 + 1) * blk]; for (a, &lv) in lrow.iter().enumerate() { if lv.is_zero() { continue; @@ -1085,53 +1088,41 @@ impl Engine<'_> { } } } - // T2[a', w1, s1', s2, b] = Σ_{w0, s1} W[w0, w1, s1', s1] T1[a', w0, s1, s2, b] + // T2[a', w1, s1', s2, b] = Σ_{w0, s1} W_i[w0, w1, s1', s1] T1[a', w0, s1, s2, b] let half = 2 * dr; - let mut t2 = vec![T::ZERO; nrows * w * blk]; + let mut t2 = vec![T::ZERO; nrows * w1n * blk]; for ap in 0..nrows { - for w0 in 0..w { - for w1 in 0..w { - if !self.nonzero(w0, w1) { - continue; - } - let o = self.op(w0, w1).0; - for s1p in 0..2 { - for s1 in 0..2 { - let m = o[s1p * 2 + s1]; - if m == 0.0 { - continue; - } - let src = (ap * w + w0) * blk + s1 * half; - let dst = (ap * w + w1) * blk + s1p * half; - for k in 0..half { - t2[dst + k] = t2[dst + k].add(t1[src + k].scale(m)); - } + for &(w0, w1, o) in self.entries(i) { + for s1p in 0..2 { + for s1 in 0..2 { + let m = o.0[s1p * 2 + s1]; + if m == 0.0 { + continue; + } + let src = (ap * w0n + w0) * blk + s1 * half; + let dst = (ap * w1n + w1) * blk + s1p * half; + for k in 0..half { + t2[dst + k] = t2[dst + k].add(t1[src + k].scale(m)); } } } } } - // T3[a', s1', w2, s2', b] = Σ_{w1, s2} W[w1, w2, s2', s2] T2[a', w1, s1', s2, b] - let mut t3 = vec![T::ZERO; nrows * 2 * w * 2 * dr]; + // T3[a', s1', w2, s2', b] = Σ_{w1, s2} W_{i+1}[w1, w2, s2', s2] T2[a', w1, s1', s2, b] + let mut t3 = vec![T::ZERO; nrows * 2 * w2n * 2 * dr]; for ap in 0..nrows { - for w1 in 0..w { - for w2 in 0..w { - if !self.nonzero(w1, w2) { - continue; - } - let o = self.op(w1, w2).0; - for s2p in 0..2 { - for s2 in 0..2 { - let m = o[s2p * 2 + s2]; - if m == 0.0 { - continue; - } - for s1p in 0..2 { - let src = (ap * w + w1) * blk + s1p * half + s2 * dr; - let dst = (((ap * 2 + s1p) * w + w2) * 2 + s2p) * dr; - for b in 0..dr { - t3[dst + b] = t3[dst + b].add(t2[src + b].scale(m)); - } + for &(w1, w2, o) in self.entries(i + 1) { + for s2p in 0..2 { + for s2 in 0..2 { + let m = o.0[s2p * 2 + s2]; + if m == 0.0 { + continue; + } + for s1p in 0..2 { + let src = (ap * w1n + w1) * blk + s1p * half + s2 * dr; + let dst = (((ap * 2 + s1p) * w2n + w2) * 2 + s2p) * dr; + for b in 0..dr { + t3[dst + b] = t3[dst + b].add(t2[src + b].scale(m)); } } } @@ -1144,9 +1135,9 @@ impl Engine<'_> { for s2p in 0..2 { for bp in 0..dr { let mut acc = T::ZERO; - for w2 in 0..w { - let trow = &t3[(((ap * 2 + s1p) * w + w2) * 2 + s2p) * dr..(((ap * 2 + s1p) * w + w2) * 2 + s2p + 1) * dr]; - let rrow = &r[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; + for w2 in 0..w2n { + let trow = &t3[(((ap * 2 + s1p) * w2n + w2) * 2 + s2p) * dr..(((ap * 2 + s1p) * w2n + w2) * 2 + s2p + 1) * dr]; + let rrow = &r[(bp * w2n + w2) * dr..(bp * w2n + w2 + 1) * dr]; for (x, y) in trow.iter().zip(rrow) { acc = acc.add(x.mul(*y)); } @@ -1158,20 +1149,21 @@ impl Engine<'_> { } } - /// Rows `rows` of `H_eff C` for the one-site tensor `C[a, s, b]`. + /// Rows `rows` of `H_eff C` for the one-site tensor `C[a, s, b]` on site + /// `i`. #[cfg(feature = "quantum_tdvp")] #[allow(clippy::too_many_arguments)] - fn apply1_rows(&self, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { - let w = self.w; + fn apply1_rows(&self, i: usize, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { + let (w0n, w1n) = (self.w(i), self.w(i + 1)); let first = rows.start; let nrows = rows.len(); // T1[a', w0, s, b] = Σ_a L[a', w0, a] C[a, s, b] let blk = 2 * dr; - let mut t1 = vec![T::ZERO; nrows * w * blk]; + let mut t1 = vec![T::ZERO; nrows * w0n * blk]; for ap in 0..nrows { - for w0 in 0..w { - let lrow = &l[((first + ap) * w + w0) * dl..((first + ap) * w + w0 + 1) * dl]; - let dst = &mut t1[(ap * w + w0) * blk..(ap * w + w0 + 1) * blk]; + for w0 in 0..w0n { + let lrow = &l[((first + ap) * w0n + w0) * dl..((first + ap) * w0n + w0 + 1) * dl]; + let dst = &mut t1[(ap * w0n + w0) * blk..(ap * w0n + w0 + 1) * blk]; for (a, &lv) in lrow.iter().enumerate() { if lv.is_zero() { continue; @@ -1182,26 +1174,20 @@ impl Engine<'_> { } } } - // T2[a', s', w1, b] = Σ_{w0, s} W[w0, w1, s', s] T1[a', w0, s, b] - let mut t2 = vec![T::ZERO; nrows * 2 * w * dr]; + // T2[a', s', w1, b] = Σ_{w0, s} W_i[w0, w1, s', s] T1[a', w0, s, b] + let mut t2 = vec![T::ZERO; nrows * 2 * w1n * dr]; for ap in 0..nrows { - for w0 in 0..w { - for w1 in 0..w { - if !self.nonzero(w0, w1) { - continue; - } - let o = self.op(w0, w1).0; - for sp in 0..2 { - for s in 0..2 { - let m = o[sp * 2 + s]; - if m == 0.0 { - continue; - } - let src = (ap * w + w0) * blk + s * dr; - let dst = ((ap * 2 + sp) * w + w1) * dr; - for b in 0..dr { - t2[dst + b] = t2[dst + b].add(t1[src + b].scale(m)); - } + for &(w0, w1, o) in self.entries(i) { + for sp in 0..2 { + for s in 0..2 { + let m = o.0[sp * 2 + s]; + if m == 0.0 { + continue; + } + let src = (ap * w0n + w0) * blk + s * dr; + let dst = ((ap * 2 + sp) * w1n + w1) * dr; + for b in 0..dr { + t2[dst + b] = t2[dst + b].add(t1[src + b].scale(m)); } } } @@ -1212,9 +1198,9 @@ impl Engine<'_> { for sp in 0..2 { for bp in 0..dr { let mut acc = T::ZERO; - for w1 in 0..w { - let trow = &t2[((ap * 2 + sp) * w + w1) * dr..((ap * 2 + sp) * w + w1 + 1) * dr]; - let rrow = &r[(bp * w + w1) * dr..(bp * w + w1 + 1) * dr]; + for w1 in 0..w1n { + let trow = &t2[((ap * 2 + sp) * w1n + w1) * dr..((ap * 2 + sp) * w1n + w1 + 1) * dr]; + let rrow = &r[(bp * w1n + w1) * dr..(bp * w1n + w1 + 1) * dr]; for (x, y) in trow.iter().zip(rrow) { acc = acc.add(x.mul(*y)); } @@ -1225,15 +1211,15 @@ impl Engine<'_> { } } - /// `out = H_eff θ` for the two-site `θ[a, s1, s2, b]` between `l` (dl) - /// and `r` (dr). Blocks of bra rows run on separate threads; each row is - /// computed in the same order whatever the split, so the result is the - /// same bits. + /// `out = H_eff θ` for the two-site `θ[a, s1, s2, b]` on sites `i, i + 1` + /// between `l` (dl) and `r` (dr). Blocks of bra rows run on separate + /// threads; each row is computed in the same order whatever the split, so + /// the result is the same bits. #[allow(clippy::too_many_arguments)] - pub(crate) fn apply(&self, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { + pub(crate) fn apply(&self, i: usize, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { let threads = Self::split_threads(threads, dl, dr); if threads <= 1 { - self.apply_rows(l, r, theta, out, dl, dr, 0..dl); + self.apply_rows(i, l, r, theta, out, dl, dr, 0..dl); return; } let row = 4 * dr; @@ -1242,19 +1228,19 @@ impl Engine<'_> { for (c, chunk) in out.chunks_mut(per * row).enumerate() { let start = c * per; let end = (start + per).min(dl); - scope.spawn(move || self.apply_rows(l, r, theta, chunk, dl, dr, start..end)); + scope.spawn(move || self.apply_rows(i, l, r, theta, chunk, dl, dr, start..end)); } }); } - /// `out = H_eff C` for the one-site `C[a, s, b]`, threaded like - /// [`Engine::apply`]. + /// `out = H_eff C` for the one-site `C[a, s, b]` on site `i`, threaded + /// like [`Engine::apply`]. #[cfg(feature = "quantum_tdvp")] #[allow(clippy::too_many_arguments)] - pub(crate) fn apply1(&self, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { + pub(crate) fn apply1(&self, i: usize, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { let threads = Self::split_threads(threads, dl, dr); if threads <= 1 { - self.apply1_rows(l, r, c, out, dl, dr, 0..dl); + self.apply1_rows(i, l, r, c, out, dl, dr, 0..dl); return; } let row = 2 * dr; @@ -1263,7 +1249,7 @@ impl Engine<'_> { for (k, chunk) in out.chunks_mut(per * row).enumerate() { let start = k * per; let end = (start + per).min(dl); - scope.spawn(move || self.apply1_rows(l, r, c, chunk, dl, dr, start..end)); + scope.spawn(move || self.apply1_rows(i, l, r, c, chunk, dl, dr, start..end)); } }); } @@ -1390,7 +1376,7 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { let mut rights: Vec> = 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]); } let mut lefts: Vec> = vec![Vec::new(); n + 1]; lefts[0] = eng.left_edge(); @@ -1405,7 +1391,7 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { 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 apply = |x: &[f64], y: &mut [f64]| eng.apply(l, r, x, y, dl, dr, cfg.threads); + let apply = |x: &[f64], y: &mut [f64]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads); let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol); energy = e; let (u, rest, disc, weights) = split(&v, dl * 2, 2 * dr, cfg.max_bond, cfg.cutoff); @@ -1415,14 +1401,14 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { mps.sites[i] = u; mps.sites[i + 1] = rest; 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); } // Right to left: optimise (i, i+1), keep the right part. for i in (0..n - 1).rev() { 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 apply = |x: &[f64], y: &mut [f64]| eng.apply(l, r, x, y, dl, dr, cfg.threads); + let apply = |x: &[f64], y: &mut [f64]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads); let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol); energy = e; // Split on the right: rows of the right factor are eigenvectors @@ -1435,7 +1421,7 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { mps.sites[i + 1] = transpose(&u, 2 * dr, k); mps.sites[i] = transpose(&rest, k, dl * 2); 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); } sweep_energies.push(energy); } @@ -1594,8 +1580,8 @@ mod tests { fn threaded_effective_hamiltonian_is_bit_identical() { // Environments big enough to split across threads. let mpo = Chain::heisenberg(8, 1.0, 0.7, 0.2).mpo(); - let eng = Engine { mpo: &mpo, w: mpo.w }; - let (dl, dr, w) = (56usize, 48usize, mpo.w); + let eng = Engine::new(&mpo); + let (dl, dr, w) = (56usize, 48usize, mpo.dims[2]); let mut s = 3u64; let mut next = move || { s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); @@ -1605,10 +1591,10 @@ mod tests { let r: Vec = (0..dr * w * dr).map(|_| next()).collect(); let theta: Vec = (0..dl * 4 * dr).map(|_| next()).collect(); let mut one = vec![0.0; theta.len()]; - eng.apply(&l, &r, &theta, &mut one, dl, dr, 1); + eng.apply(2, &l, &r, &theta, &mut one, dl, dr, 1); for t in [2, 3, 7] { let mut many = vec![0.0; theta.len()]; - eng.apply(&l, &r, &theta, &mut many, dl, dr, t); + eng.apply(2, &l, &r, &theta, &mut many, dl, dr, t); assert!(one.iter().zip(&many).all(|(a, b)| a.to_bits() == b.to_bits()), "threads={t}"); } } diff --git a/wai/wai-quantum/src/quantum_tdvp.rs b/wai/wai-quantum/src/quantum_tdvp.rs index 74fffe8c6..8feea7c85 100644 --- a/wai/wai-quantum/src/quantum_tdvp.rs +++ b/wai/wai-quantum/src/quantum_tdvp.rs @@ -298,7 +298,7 @@ impl Evolution { let mut rights: Vec> = 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), ¢re, -half, cfg.krylov, cfg.tol) + expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(i, l, r, x, y, dl, k, cfg.threads), ¢re, -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 } @@ -1052,7 +1052,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 = (0..dl * w * dl).map(|_| c()).collect(); @@ -1062,13 +1062,13 @@ mod tests { let r1: Vec = 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}"); } -- 2.54.0 (Apple Git-157) From 9c5a0b0ecc5dcb30711f82ff050e7f0598bb7d76 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 5 Oct 2026 11:39:59 -0400 Subject: [PATCH 2/3] wai-quantum: exact MPOs from sums of operator strings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit OpSum compiles a sum of products of one-site operators into an MPO by the bipartite-graph construction: a minimum vertex cover per bond (Hopcroft-Karp matching, König's theorem) sets the bond, and every coefficient is placed once, so the MPO is exact. dmrg_mpo runs DMRG on any MPO from any start. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/src/quantum_dmrg.rs | 401 +++++++++++++++++++++++++++- 1 file changed, 396 insertions(+), 5 deletions(-) diff --git a/wai/wai-quantum/src/quantum_dmrg.rs b/wai/wai-quantum/src/quantum_dmrg.rs index 81f49e9f6..8eee36b04 100644 --- a/wai/wai-quantum/src/quantum_dmrg.rs +++ b/wai/wai-quantum/src/quantum_dmrg.rs @@ -591,6 +591,317 @@ impl Mpo { pub fn max_bond(&self) -> usize { self.dims.iter().copied().max().unwrap_or(1) } + + /// The operator as a dense `2ⁿ × 2ⁿ` matrix, row-major, site 0 the most + /// significant bit (at most 10 sites). + pub fn to_dense(&self) -> Option> { + let n = self.n; + if n == 0 || n > 10 { + return None; + } + // acc[(r · d + c) · w + a] over the sites so far. + let mut acc = self.left.clone(); + let mut d = 1usize; + for i in 0..n { + let (wl, wr) = (self.dims[i], self.dims[i + 1]); + let nd = d * 2; + let mut next = vec![0.0; nd * nd * wr]; + for r in 0..d { + for c in 0..d { + for &(a, b, o) in &self.sites[i] { + let x = acc[(r * d + c) * wl + a]; + if x == 0.0 { + continue; + } + for sp in 0..2 { + for sq in 0..2 { + let m = o.0[sp * 2 + sq]; + if m != 0.0 { + next[((r * 2 + sp) * nd + c * 2 + sq) * wr + b] += x * m; + } + } + } + } + } + } + acc = next; + d = nd; + } + Some((0..d * d).map(|rc| (0..self.dims[n]).map(|a| acc[rc * self.dims[n] + a] * self.right[a]).sum()).collect()) + } +} + +// --------------------------------------------------------------------------- +// Sums of operator strings, compiled into MPOs +// --------------------------------------------------------------------------- + +impl Op { + /// The matrix product `self · o`. + pub fn matmul(self, o: Op) -> Op { + let (a, b) = (self.0, o.0); + Op([a[0] * b[0] + a[1] * b[2], a[0] * b[1] + a[1] * b[3], a[2] * b[0] + a[3] * b[2], a[2] * b[1] + a[3] * b[3]]) + } +} + +/// A Hamiltonian as a sum of products of one-site operators, compiled into +/// an exact MPO. +/// +/// The construction is the bipartite-graph one. At each bond, a term's left +/// half (an MPO state and the operator on the site just passed) and its right +/// half (the operators still to come) are the two sides of a graph, and the +/// term is an edge. A minimum vertex cover of that graph, from a maximum +/// matching by König's theorem, sets the bond: a covered left half becomes a +/// state and carries its terms' coefficients on; a covered right half takes +/// the sum of the uncovered left halves that reach it, coefficients included. +/// Every coefficient is placed exactly once, so the MPO is exact, with no +/// numerical compression. +#[derive(Clone, Debug, PartialEq)] +pub struct OpSum { + n: usize, + /// One-site operators, interned by their bits. + table: Vec, + /// Each distinct string of operator indices and its summed coefficient. + terms: std::collections::BTreeMap, f64>, +} + +impl OpSum { + pub fn new(n: usize) -> OpSum { + let mut s = OpSum { n, table: Vec::new(), terms: std::collections::BTreeMap::new() }; + s.intern(Op::id()); + s + } + + fn intern(&mut self, op: Op) -> u8 { + let bits = op.0.map(f64::to_bits); + if let Some(k) = self.table.iter().position(|o| o.0.map(f64::to_bits) == bits) { + return k as u8; + } + assert!(self.table.len() < 255, "too many distinct one-site operators"); + self.table.push(op); + (self.table.len() - 1) as u8 + } + + /// Add `coeff · F₁ F₂ ⋯` where each factor is a one-site operator on a + /// site, in operator order. Factors on the same site multiply in that + /// order; sites not named carry the identity. + pub fn add(&mut self, coeff: f64, factors: &[(usize, Op)]) { + if coeff == 0.0 { + return; + } + let mut local = vec![Op::id(); self.n]; + for &(site, op) in factors { + assert!(site < self.n, "site {site} is past the chain"); + local[site] = local[site].matmul(op); + } + if local.iter().any(|o| o.0.iter().all(|&v| v == 0.0)) { + return; + } + let key: Vec = local.into_iter().map(|o| self.intern(o)).collect(); + *self.terms.entry(key).or_insert(0.0) += coeff; + } + + /// The number of distinct operator strings. + pub fn len(&self) -> usize { + self.terms.len() + } + + pub fn is_empty(&self) -> bool { + self.terms.is_empty() + } + + /// Compile into an MPO. + pub fn mpo(&self) -> Mpo { + let n = self.n; + assert!(n >= 1, "an MPO needs a site"); + // Terms at the current bond: (MPO state, coefficient, index into the + // owning string where the right half starts). Strings stay in place. + let strings: Vec<(&Vec, f64)> = self.terms.iter().filter(|e| *e.1 != 0.0).map(|(k, &c)| (k, c)).collect(); + let mut live: Vec<(usize, f64, usize)> = (0..strings.len()).map(|t| (0usize, strings[t].1, t)).collect(); + let mut dims = vec![1usize; n + 1]; + let mut sites: Vec> = Vec::with_capacity(n); + for k in 0..n { + // Left vertices (state, op at k); right vertices (ops after k). + let mut lefts: std::collections::BTreeMap<(usize, u8), usize> = std::collections::BTreeMap::new(); + let mut rights: std::collections::BTreeMap<&[u8], usize> = std::collections::BTreeMap::new(); + for &(state, _, t) in &live { + let s = strings[t].0; + lefts.entry((state, s[k])).or_insert(0); + rights.entry(&s[k + 1..]).or_insert(0); + } + for (i, v) in lefts.values_mut().enumerate() { + *v = i; + } + for (i, v) in rights.values_mut().enumerate() { + *v = i; + } + let (nu, nv) = (lefts.len(), rights.len()); + // Edges with their summed weights, in a fixed order. + let mut edges: std::collections::BTreeMap<(usize, usize), f64> = std::collections::BTreeMap::new(); + for &(state, c, t) in &live { + let s = strings[t].0; + let u = lefts[&(state, s[k])]; + let v = rights[&s[k + 1..]]; + *edges.entry((u, v)).or_insert(0.0) += c; + } + let mut adj = vec![Vec::new(); nu]; + for &(u, v) in edges.keys() { + adj[u].push(v); + } + let (in_u, in_v) = min_vertex_cover(nu, nv, &adj); + // New states: covered left vertices first, then covered right. + let mut state_of_u = vec![usize::MAX; nu]; + let mut state_of_v = vec![usize::MAX; nv]; + let mut next = 0usize; + for u in 0..nu { + if in_u[u] { + state_of_u[u] = next; + next += 1; + } + } + for v in 0..nv { + if in_v[v] { + state_of_v[v] = next; + next += 1; + } + } + let left_keys: Vec<(usize, u8)> = lefts.keys().copied().collect(); + let right_keys: Vec<&[u8]> = rights.keys().copied().collect(); + let mut entries: std::collections::BTreeMap<(usize, usize), Op> = std::collections::BTreeMap::new(); + let mut add_entry = |a: usize, b: usize, op: Op| { + let e = entries.entry((a, b)).or_insert(Op([0.0; 4])); + *e = e.add(op); + }; + for u in 0..nu { + if in_u[u] { + let (a, o) = left_keys[u]; + add_entry(a, state_of_u[u], self.table[o as usize]); + } + } + // Next bond's terms: through a covered left vertex the + // coefficient travels on; into a covered right vertex it is + // absorbed into the entry. + let mut carried: std::collections::BTreeMap<(usize, usize), f64> = std::collections::BTreeMap::new(); + for (&(u, v), &w) in &edges { + if in_u[u] { + *carried.entry((state_of_u[u], v)).or_insert(0.0) += w; + } else { + let (a, o) = left_keys[u]; + add_entry(a, state_of_v[v], self.table[o as usize].scaled(w)); + carried.entry((state_of_v[v], v)).or_insert(1.0); + } + } + // A string's right half is identified by its vertex; find one + // string carrying each so the next bond can read its operators. + let mut owner = vec![usize::MAX; nv]; + for &(_, _, t) in &live { + let v = rights[&strings[t].0[k + 1..]]; + if owner[v] == usize::MAX { + owner[v] = t; + } + } + let _ = right_keys; + live = carried.into_iter().filter(|&(_, c)| c != 0.0).map(|((b, v), c)| (b, c, owner[v])).collect(); + dims[k + 1] = next; + sites.push(entries.into_iter().filter(|(_, o)| o.0.iter().any(|&x| x != 0.0)).map(|((a, b), o)| (a, b, o)).collect()); + } + // After the last site every right half is empty: the remaining + // coefficients form the right boundary. + let mut right = vec![0.0; dims[n]]; + for &(state, c, _) in &live { + right[state] += c; + } + Mpo { n, dims, sites, left: vec![1.0], right } + } +} + +/// A minimum vertex cover of a bipartite graph (`adj[u]` lists the right +/// vertices of left vertex `u`): a maximum matching by Hopcroft–Karp, then +/// König's construction. Returns membership for each side. +fn min_vertex_cover(nu: usize, nv: usize, adj: &[Vec]) -> (Vec, Vec) { + const FREE: usize = usize::MAX; + let mut mu = vec![FREE; nu]; + let mut mv = vec![FREE; nv]; + let mut dist = vec![0usize; nu]; + loop { + // Layer the free left vertices and what alternating paths reach. + let mut queue = std::collections::VecDeque::new(); + for u in 0..nu { + if mu[u] == FREE { + dist[u] = 0; + queue.push_back(u); + } else { + dist[u] = usize::MAX; + } + } + let mut found = false; + while let Some(u) = queue.pop_front() { + for &v in &adj[u] { + let w = mv[v]; + if w == FREE { + found = true; + } else if dist[w] == usize::MAX { + dist[w] = dist[u] + 1; + queue.push_back(w); + } + } + } + if !found { + break; + } + // Augment along vertex-disjoint shortest paths, iteratively. + let mut it = vec![0usize; nu]; + for root in 0..nu { + if mu[root] != FREE { + continue; + } + let mut stack = vec![root]; + while let Some(&u) = stack.last() { + if it[u] == adj[u].len() { + dist[u] = usize::MAX; + stack.pop(); + continue; + } + let v = adj[u][it[u]]; + it[u] += 1; + let w = mv[v]; + if w == FREE { + // Flip the path root → … → u → v. + let mut v = v; + while let Some(u) = stack.pop() { + let prev = mu[u]; + mu[u] = v; + mv[v] = u; + v = prev; + } + break; + } else if dist[w] == dist[u] + 1 { + stack.push(w); + } + } + } + } + // König: Z = vertices reachable from free left vertices by alternating + // paths; the cover is (left ∖ Z) ∪ (right ∩ Z). + let mut zu = vec![false; nu]; + let mut zv = vec![false; nv]; + let mut queue: std::collections::VecDeque = (0..nu).filter(|&u| mu[u] == FREE).collect(); + for &u in &queue { + zu[u] = true; + } + while let Some(u) = queue.pop_front() { + for &v in &adj[u] { + if mu[u] == v || zv[v] { + continue; + } + zv[v] = true; + let w = mv[v]; + if w != FREE && !zu[w] { + zu[w] = true; + queue.push_back(w); + } + } + } + ((0..nu).map(|u| !zu[u]).collect(), zv) } // --------------------------------------------------------------------------- @@ -599,7 +910,7 @@ impl Mpo { /// Lanczos for the lowest eigenpair of a symmetric operator given by /// `apply`, from `start`, with full re-orthogonalisation. -fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, restarts: usize, tol: f64) -> (f64, Vec) { +pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, restarts: usize, tol: f64) -> (f64, Vec) { let dim = start.len(); let norm = |v: &[f64]| v.iter().map(|x| x * x).sum::().sqrt(); let mut v = start.to_vec(); @@ -1365,11 +1676,17 @@ pub(crate) fn right_canonicalise(mps: &mut Mps) { /// Two-site DMRG for the ground state of `chain`. pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { - let n = chain.n; + assert!(chain.n >= 2, "DMRG needs at least two sites"); + dmrg_mpo(&chain.mpo(), Mps::random(chain.n, cfg.max_bond.min(8), cfg.seed), cfg) +} + +/// Two-site DMRG for the ground state of any MPO, from `start`. +pub fn dmrg_mpo(mpo: &Mpo, start: Mps, cfg: &DmrgConfig) -> DmrgResult { + let n = mpo.n; assert!(n >= 2, "DMRG needs at least two sites"); - let mpo = chain.mpo(); - let eng = Engine::new(&mpo); - let mut mps = Mps::random(n, cfg.max_bond.min(8), cfg.seed); + assert_eq!(start.n(), n, "the start state and the MPO differ in length"); + let eng = Engine::new(mpo); + let mut mps = start; // Bring the state into right-canonical form from the right end, and // build the right environments. right_canonicalise(&mut mps); @@ -1599,6 +1916,80 @@ mod tests { } } + #[test] + fn opsum_builds_the_operator_exactly() { + let mut s = 21u64; + let mut next = move || { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (s >> 11) as f64 / 9_007_199_254_740_992.0 + }; + let ops = [Op::x(), Op::z(), Op::plus(), Op::minus(), Op::sz(), Op([0.0, 0.0, 0.0, 1.0])]; + let n = 6; + for round in 0..4 { + let mut sum = OpSum::new(n); + let mut terms: Vec<(f64, Vec<(usize, Op)>)> = Vec::new(); + for _ in 0..40 + 20 * round { + let k = 1 + (next() * 4.0) as usize; + let factors: Vec<(usize, Op)> = (0..k).map(|_| ((next() * n as f64) as usize, ops[(next() * ops.len() as f64) as usize])).collect(); + let c = next() - 0.5; + sum.add(c, &factors); + terms.push((c, factors)); + } + sum.add(0.75, &[]); + terms.push((0.75, Vec::new())); + // Direct: Σ c ⊗ (per-site products). + let dim = 1usize << n; + let mut want = vec![0.0; dim * dim]; + for (c, factors) in &terms { + let mut local = vec![Op::id(); n]; + for &(site, op) in factors { + local[site] = local[site].matmul(op); + } + for r in 0..dim { + for col in 0..dim { + let mut v = *c; + for (i, o) in local.iter().enumerate() { + let (sr, sc) = ((r >> (n - 1 - i)) & 1, (col >> (n - 1 - i)) & 1); + v *= o.0[sr * 2 + sc]; + } + want[r * dim + col] += v; + } + } + } + let mpo = sum.mpo(); + let got = mpo.to_dense().unwrap(); + for (g, w) in got.iter().zip(&want) { + assert!((g - w).abs() < 1e-12, "round {round}"); + } + assert_eq!(mpo.dims[0], 1); + assert_eq!(mpo.dims[n], 1); + } + } + + #[test] + fn opsum_reproduces_the_chains() { + let n = 12; + let (j, jz, h) = (1.0, 0.7, 0.2); + let mut sum = OpSum::new(n); + for i in 0..n - 1 { + sum.add(j / 2.0, &[(i, Op::plus()), (i + 1, Op::minus())]); + sum.add(j / 2.0, &[(i, Op::minus()), (i + 1, Op::plus())]); + sum.add(jz, &[(i, Op::sz()), (i + 1, Op::sz())]); + } + for i in 0..n { + sum.add(-h, &[(i, Op::sz())]); + } + let mpo = sum.mpo(); + // The bond the hand-built automaton uses: one per bond term, plus two. + assert!(mpo.max_bond() <= 5, "{:?}", mpo.dims); + let chain = Chain::heisenberg(n, j, jz, h); + let cfg = DmrgConfig { max_bond: 32, sweeps: 4, ..DmrgConfig::default() }; + let a = dmrg(&chain, &cfg).energy; + let b = dmrg_mpo(&mpo, Mps::random(n, 8, cfg.seed), &cfg).energy; + assert!((a - b).abs() < 1e-9, "{a} vs {b}"); + assert!((a - exact_ground_energy(&chain).unwrap()).abs() < 1e-8); + } + #[test] fn energy_falls_sweep_by_sweep_and_runs_repeat() { let chain = Chain::heisenberg(16, 1.0, 1.0, 0.0); -- 2.54.0 (Apple Git-157) From 90a490b9e7b1bda0c7e690453b0833dbb5905bfe Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 5 Oct 2026 16:15:50 -0400 Subject: [PATCH 3/3] wai-quantum 0.3.36: molecules and neural quantum states, held to exact answers quantum_chem: a molecule's ground state by DMRG from its FCIDUMP integrals. - The second-quantised Hamiltonian lives on spin-orbital sites (alpha and beta side by side) via Jordan-Wigner, and compiles exactly into an MPO through OpSum. A sector penalty costs no bond, because its terms merge with the Coulomb ones. - DMRG starts from Hartree-Fock with complete bases and keeps them complete (cutoff 0). With a positive cutoff it stalled at Hartree-Fock (Brillouin) or at the state's own rank. - Full CI by direct CI over alpha/beta strings is the referee; it matches an independent program to 1e-9. - Results: LiH and H2O reach full CI to 1e-13 in one sweep. N2 (20 spin orbitals) at bond 64 comes within 3.6e-4. quantum_nqs: neural quantum states. - Spin models on any graph, with exact diagonalisation and an OpSum form for DMRG. - A real RBM with Marshall's sign; variational Monte Carlo by Metropolis chains or exact enumeration. - Stochastic reconfiguration, solved in parameter or sample space, which give the same step. - Results: 4x4 Heisenberg comes within 4.6e-4 of exact diagonalisation, and the critical 4x4 Ising model within its error bar. On 6x6 it meets DMRG at 3.6e-3, still descending. DMRG's Lanczos now stops on the Ritz residual estimate, which re-pinned DMRG. New pins chemistry_by_dmrg and neural_quantum_states pass natively and under wasm32-wasip2. The CLI gains chem and nqs, byte-identical under WASI. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 19 +- wai/wai-quantum/README.md | 31 + wai/wai-quantum/examples/chem_dmrg.rs | 33 + wai/wai-quantum/examples/nqs_ground.rs | 57 ++ wai/wai-quantum/src/lib.rs | 6 + wai/wai-quantum/src/quantum_chem.rs | 607 ++++++++++++++++ wai/wai-quantum/src/quantum_dmrg.rs | 128 +++- wai/wai-quantum/src/quantum_nqs.rs | 922 +++++++++++++++++++++++++ wai/wai-quantum/src/quantum_tdvp.rs | 52 +- wai/wai-quantum/src/repro.rs | 6 +- wai/wai-quantum/tests/bit_identity.rs | 45 +- wai/wai-quantum/tests/data/h2.fcidump | 12 + wai/wai-quantum/tests/data/h2o.fcidump | 400 +++++++++++ wai/wai-quantum/tests/data/h6.fcidump | 255 +++++++ wai/wai-quantum/tests/data/lih.fcidump | 194 ++++++ wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 2 +- wai/wai-rs/src/bin/wai_quantum.rs | 99 +++ wai/wai-rs/src/lib.rs | 6 + wai/wai-web/Cargo.lock | 2 +- 22 files changed, 2815 insertions(+), 67 deletions(-) create mode 100644 wai/wai-quantum/examples/chem_dmrg.rs create mode 100644 wai/wai-quantum/examples/nqs_ground.rs create mode 100644 wai/wai-quantum/src/quantum_chem.rs create mode 100644 wai/wai-quantum/src/quantum_nqs.rs create mode 100644 wai/wai-quantum/tests/data/h2.fcidump create mode 100644 wai/wai-quantum/tests/data/h2o.fcidump create mode 100644 wai/wai-quantum/tests/data/h6.fcidump create mode 100644 wai/wai-quantum/tests/data/lih.fcidump diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index e719b744d..ea3a590c1 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/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", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 9368b760c..835fc4e62 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/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", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index e4419e5aa..dac363755 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/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"] diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 5fa68ca94..9544d8bf0 100644 --- a/wai/wai-quantum/README.md +++ b/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 diff --git a/wai/wai-quantum/examples/chem_dmrg.rs b/wai/wai-quantum/examples/chem_dmrg.rs new file mode 100644 index 000000000..33a1472a5 --- /dev/null +++ b/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 = 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)"), + } +} diff --git a/wai/wai-quantum/examples/nqs_ground.rs b/wai/wai-quantum/examples/nqs_ground.rs new file mode 100644 index 000000000..79cf6ca9e --- /dev/null +++ b/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 = 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::() / tail.len() as f64); +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index af3bb2bf9..b95daf6f3 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/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; diff --git a/wai/wai-quantum/src/quantum_chem.rs b/wai/wai-quantum/src/quantum_chem.rs new file mode 100644 index 000000000..9f0805d60 --- /dev/null +++ b/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, + /// `(pq|rs)` in chemists' order, `[((p·n + q)·n + r)·n + s]`. + pub h2: Vec, + /// 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 { + 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 { + 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 { + 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> { + 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 { + 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 = (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 = (0..dim).filter(|x: &usize| x.count_ones() == 2 && (x >> 3 & 1) + (x >> 1 & 1) == 1).collect(); + let k = two.len(); + let block: Vec = 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); + } +} diff --git a/wai/wai-quantum/src/quantum_dmrg.rs b/wai/wai-quantum/src/quantum_dmrg.rs index 8eee36b04..74b0bc162 100644 --- a/wai/wai-quantum/src/quantum_dmrg.rs +++ b/wai/wai-quantum/src/quantum_dmrg.rs @@ -11,9 +11,16 @@ //! are the XXZ Heisenberg chain in a field and the transverse-field Ising //! chain. [`Chain::mpo`] writes it as a matrix-product operator, the finite //! state machine of its terms. +//! - **Any Hamiltonian as an MPO.** [`OpSum`] compiles a sum of products of +//! one-site operators, at any range, into an exact MPO with a tensor per +//! site. At each bond, a minimum vertex cover of the graph joining the +//! terms' left and right halves sets the bond, and every coefficient is +//! placed exactly once. [`dmrg_mpo`] runs DMRG on any MPO from any start. //! - **Two-site DMRG.** [`dmrg`] sweeps a matrix-product state left and right. //! - At each bond it finds the lowest eigenvector of the effective -//! Hamiltonian by Lanczos, with full re-orthogonalisation. +//! Hamiltonian by Lanczos, with full re-orthogonalisation. Lanczos stops +//! as soon as the lowest Ritz pair's residual estimate meets the +//! tolerance. //! - It splits that vector by the eigenvectors of its reduced density //! matrix. It keeps at most `max_bond` of them and drops those below a //! weight cutoff. @@ -39,6 +46,10 @@ //! //! - **The eigensolvers.** They reconstruct random symmetric and Hermitian //! matrices up to 40 × 40 to `10⁻¹²` with orthonormal vectors. +//! - **The MPO compiler.** Random sums of operator strings, up to four +//! factors at any range, match their dense matrices to `10⁻¹²`. The +//! Heisenberg chain rebuilt from its terms comes out at the hand-built +//! automaton's bond of five, with the same DMRG energy. //! - **The free-fermion solution.** It equals exact diagonalisation of the //! Ising chain at three couplings, for 4, 7 and 10 sites. //! - **DMRG against exact diagonalisation.** On Heisenberg chains (isotropic, @@ -62,9 +73,10 @@ //! //! # Honest boundaries //! -//! - **One dimension, nearest neighbours, real Hamiltonians.** Longer-range -//! couplings, complex terms and two-dimensional lattices are not handled. -//! Time evolution is the TDVP module's, on this engine. +//! - **One dimension, two states a site, real Hamiltonians.** Longer-range +//! and two-dimensional couplings go through [`OpSum`], at a bond that +//! grows with their reach; complex terms are not handled. Time evolution is +//! the TDVP module's, and molecules the chemistry module's, on this engine. //! - **No symmetry sectors.** Conserved quantities are not used to block the //! matrices, so a run costs more than a symmetric code would. //! - **Variational.** The energy is an upper bound, converging from above as @@ -941,6 +953,24 @@ pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: if j + 1 == k || bnorm < 1e-13 { break; } + // The lowest Ritz pair's residual is β times the last component + // of its vector; once that meets the tolerance, more vectors + // cannot help. + if j >= 1 { + let m = alpha.len(); + let mut t = vec![0.0; m * m]; + for i in 0..m { + t[i * m + i] = alpha[i]; + if i + 1 < m { + t[i * m + i + 1] = beta[i]; + t[(i + 1) * m + i] = beta[i]; + } + } + let (_, z) = symmetric_eigen(&t, m); + if bnorm * z[(m - 1) * m].abs() < tol { + break; + } + } beta.push(bnorm); basis.push(w.iter().map(|x| x / bnorm).collect()); } @@ -1674,6 +1704,57 @@ pub(crate) fn right_canonicalise(mps: &mut Mps) { } } +/// 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|ψ⟩`. +#[cfg(any(feature = "quantum_tdvp", feature = "quantum_chem"))] +pub(crate) fn complete_bases(mps: &mut Mps, 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> = 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![T::ZERO; cols]; + v[unit] = T::ONE; + // Gram–Schmidt against the rows so far, twice. + for _ in 0..2 { + for r in &rows { + let c = r.iter().zip(&v).fold(T::ZERO, |acc, (x, y)| acc.add(x.conj().mul(*y))); + for (x, y) in v.iter_mut().zip(r) { + *x = x.sub(c.mul(*y)); + } + } + } + let nv = v.iter().map(|x| x.norm2()).sum::().sqrt(); + 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![T::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; + } +} + /// Two-site DMRG for the ground state of `chain`. pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { assert!(chain.n >= 2, "DMRG needs at least two sites"); @@ -1745,6 +1826,45 @@ pub fn dmrg_mpo(mpo: &Mpo, start: Mps, cfg: &DmrgConfig) -> DmrgResult { DmrgResult { energy, sweep_energies, discarded, entropies, mps } } +/// `⟨ψ|H|ψ⟩ / ⟨ψ|ψ⟩` for any state and MPO of the same length, canonical +/// or not. +pub fn expectation(mpo: &Mpo, psi: &Mps) -> f64 { + let n = psi.n(); + assert_eq!(mpo.n, n, "the state and the MPO differ in length"); + let eng = Engine::new(mpo); + let mut l: Vec = eng.left_edge(); + // ⟨ψ|ψ⟩ alongside, through the transfer matrix. + let mut e = vec![T::ONE]; + for i in 0..n { + let (dl, dr) = (psi.dims[i], psi.dims[i + 1]); + let m = &psi.sites[i]; + l = eng.grow_left(i, &l, m, dl, dr); + let mut next = vec![T::ZERO; dr * dr]; + for ap in 0..dl { + for a in 0..dl { + let x = e[ap * dl + a]; + if x.is_zero() { + continue; + } + for s in 0..2 { + for bp in 0..dr { + let bra = m[(ap * 2 + s) * dr + bp].conj().mul(x); + if bra.is_zero() { + continue; + } + for b in 0..dr { + next[bp * dr + b] = next[bp * dr + b].add(bra.mul(m[(a * 2 + s) * dr + b])); + } + } + } + } + } + e = next; + } + let num = l.iter().zip(&mpo.right).fold(T::ZERO, |acc, (x, &r)| acc.add(x.scale(r))); + num.re() / e[0].re() +} + /// `θ[a, s1, s2, b] = Σ_c A[a, s1, c] B[c, s2, b]`. pub(crate) fn two_site(a: &[T], b: &[T], dl: usize, dm: usize, dr: usize) -> Vec { let mut t = vec![T::ZERO; dl * 4 * dr]; diff --git a/wai/wai-quantum/src/quantum_nqs.rs b/wai/wai-quantum/src/quantum_nqs.rs new file mode 100644 index 000000000..8221e5250 --- /dev/null +++ b/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 { + (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, &) 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) -> Option { + 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 = (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, + pub b: Vec, + /// `W[j·n + i]`. + pub w: Vec, + pub sublattice: Option>, +} + +/// 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) -> 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 { + (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, + /// 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, Rng)>, +} + +/// Samples: configurations, their weights (summing to one), local energies +/// and log-derivatives (`params` per row). +struct Batch { + weights: Vec, + energies: Vec, + derivs: Vec, + 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 = (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, logs: &mut [f64], energies: &mut [f64], derivs: &mut [f64]| { + let start = range.start; + for idx in range { + let x = configs[idx]; + let s: Vec = (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 = 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, Rng)| -> (Vec, Vec, 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, Vec, 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::>()) + }) + .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 { + 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 { + { + // (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, k: usize, p: usize, shift: f64) -> Vec { + { + // δ = 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, rng: &mut Rng) -> Vec { + 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 = (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 = 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 = (0..k * p).map(|_| rng.uniform() - 0.5).collect(); + let eps: Vec = (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); + } +} + diff --git a/wai/wai-quantum/src/quantum_tdvp.rs b/wai/wai-quantum/src/quantum_tdvp.rs index 8feea7c85..3263f5e03 100644 --- a/wai/wai-quantum/src/quantum_tdvp.rs +++ b/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; // --------------------------------------------------------------------------- @@ -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, 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> = 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` diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index ee1f70afa..36bf82eb9 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/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; diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index 752ad9ce4..deed66bd0 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/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); +} + diff --git a/wai/wai-quantum/tests/data/h2.fcidump b/wai/wai-quantum/tests/data/h2.fcidump new file mode 100644 index 000000000..25c428564 --- /dev/null +++ b/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 diff --git a/wai/wai-quantum/tests/data/h2o.fcidump b/wai/wai-quantum/tests/data/h2o.fcidump new file mode 100644 index 000000000..c1c7ec989 --- /dev/null +++ b/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 diff --git a/wai/wai-quantum/tests/data/h6.fcidump b/wai/wai-quantum/tests/data/h6.fcidump new file mode 100644 index 000000000..3bd3de1d4 --- /dev/null +++ b/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 diff --git a/wai/wai-quantum/tests/data/lih.fcidump b/wai/wai-quantum/tests/data/lih.fcidump new file mode 100644 index 000000000..425a26b79 --- /dev/null +++ b/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 diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index c89a431b9..c8aadced7 100644 --- a/wai/wai-rs/Cargo.lock +++ b/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", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 45a21bd67..3d73130c7 100644 --- a/wai/wai-rs/Cargo.toml +++ b/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 diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index a4f2c7105..8661fbd65 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/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 = 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. diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index 399e9fe68..86f4df9a9 100644 --- a/wai/wai-rs/src/lib.rs +++ b/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")] diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index fbbb16335..382807f5b 100644 --- a/wai/wai-web/Cargo.lock +++ b/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", -- 2.54.0 (Apple Git-157)