integral 0.1.0

Native-Rust Gaussian integrals for quantum mechanics (driver + public API).
Documentation
//! Property and layout tests for the two-electron (ERI) builder.
//!
//! Pure-Rust, no external library. These tests assert the
//! permutational symmetry and the documented block layout/strides — the latter
//! matters now that ERI blocks are *not* index-symmetric the way the Phase-1
//! operators were.

use integral::{Basis, Shell};
use integral_math::am::n_cart;

/// Four shells of different angular momenta on different centers, so every
/// (ij|kl) block is genuinely distinct (no accidental symmetry to hide a
/// transposition bug).
fn mixed_basis() -> Basis {
    Basis::new(vec![
        Shell::new(0, [0.0, 0.0, 0.0], vec![1.4, 0.5], vec![0.6, 0.5]).unwrap(), // s
        Shell::new(1, [0.6, -0.3, 0.2], vec![0.9], vec![1.0]).unwrap(),          // p
        Shell::new(2, [-0.4, 0.7, -0.1], vec![1.1], vec![1.0]).unwrap(),         // d
        Shell::new(3, [0.2, 0.5, 0.8], vec![0.7], vec![1.0]).unwrap(),           // f
    ])
}

fn nao(b: &Basis) -> usize {
    b.nao_cart()
}

#[test]
fn ssss_is_positive_and_finite() {
    // A single s shell: (00|00) must be a finite positive number.
    let basis = Basis::new(vec![
        Shell::new(0, [0.0, 0.0, 0.0], vec![0.8], vec![1.0]).unwrap()
    ]);
    let eri = basis.eri();
    assert_eq!(eri.len(), 1);
    assert!(eri[0].is_finite() && eri[0] > 0.0, "(ss|ss) = {}", eri[0]);
}

#[test]
fn eight_fold_permutational_symmetry() {
    // Genuine assertion: the tensor is built block-by-block with no
    // symmetrization, so equality of the 8 permutations is a real check.
    let basis = mixed_basis();
    let n = nao(&basis);
    let eri = basis.eri();
    let at = |i: usize, j: usize, k: usize, l: usize| eri[((i * n + j) * n + k) * n + l];

    let mut worst = 0.0_f64;
    for i in 0..n {
        for j in 0..n {
            for k in 0..n {
                for l in 0..n {
                    let base = at(i, j, k, l);
                    let perms = [
                        at(j, i, k, l),
                        at(i, j, l, k),
                        at(j, i, l, k),
                        at(k, l, i, j),
                        at(l, k, i, j),
                        at(k, l, j, i),
                        at(l, k, j, i),
                    ];
                    for &p in &perms {
                        let d = (base - p).abs();
                        worst = worst.max(d / base.abs().max(1e-300));
                        assert!(
                            d <= 1e-11 * base.abs().max(1e-12),
                            "8-fold symmetry broken at ({i}{j}|{k}{l}): {base} vs {p}"
                        );
                    }
                }
            }
        }
    }
    assert!(worst < 1e-11, "worst symmetry residual {worst:e}");
}

#[test]
fn block_layout_matches_dense_strides() {
    // `eri_block(i,j,k,l)` must equal the dense tensor slice at the four shell
    // offsets under the documented row-major (a,b,c,d) layout. Use all-distinct,
    // different-L shells so a stride/transposition bug cannot stay hidden.
    let basis = mixed_basis();
    let n = nao(&basis);
    let dense = basis.eri();
    let shells = basis.shells();
    let offs: Vec<usize> = {
        let mut acc = 0;
        let mut v = Vec::new();
        for s in shells {
            v.push(acc);
            acc += s.n_cart();
        }
        v
    };

    for (qi, qj, qk, ql) in [(1usize, 0usize, 2usize, 3usize), (0, 1, 3, 2), (2, 3, 1, 0)] {
        let block = basis.eri_block(qi, qj, qk, ql);
        let (nb, nc, nd) = (
            n_cart(shells[qj].l()),
            n_cart(shells[qk].l()),
            n_cart(shells[ql].l()),
        );
        let (na2,) = (n_cart(shells[qi].l()),);
        assert_eq!(block.len(), na2 * nb * nc * nd);
        for a in 0..na2 {
            for b in 0..nb {
                for c in 0..nc {
                    for d in 0..nd {
                        let src = ((a * nb + b) * nc + c) * nd + d;
                        let (i, j, k, l) = (offs[qi] + a, offs[qj] + b, offs[qk] + c, offs[ql] + d);
                        let dval = dense[((i * n + j) * n + k) * n + l];
                        assert!(
                            (block[src] - dval).abs() <= 1e-12 * dval.abs().max(1e-12),
                            "layout mismatch at block {src} ({i}{j}|{k}{l}): {} vs {}",
                            block[src],
                            dval
                        );
                    }
                }
            }
        }
    }
}