Merge pull request 'test(symclaw-skill): cover handlers_advanced via JSON API' (#10) from ci-doctor/coverage-20260518-201834 into master
Reviewed-on: #10
This commit is contained in:
redclawsystems
@@ -0,0 +1,382 @@
|
||||
//! Molecular integrals over Gaussian basis functions.
|
||||
//!
|
||||
//! Gaussian basis functions have the form:
|
||||
//! φ(r; α, l, m, n, A) = N·(x-Ax)^l·(y-Ay)^m·(z-Az)^n·exp(-α|r-A|²)
|
||||
//!
|
||||
//! All one-electron integrals have closed analytic forms via
|
||||
//! the Obara-Saika or McMurchie-Davidson recurrence.
|
||||
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::f64::consts::PI;
|
||||
|
||||
/// A Gaussian basis function specification.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct GaussianBasis {
|
||||
/// Exponent α (decay rate).
|
||||
pub alpha: f64,
|
||||
/// Angular momentum quantum numbers (l, m, n).
|
||||
pub lmn: (u32, u32, u32),
|
||||
/// Center coordinates A = (Ax, Ay, Az).
|
||||
pub center: [f64; 3],
|
||||
}
|
||||
|
||||
impl GaussianBasis {
|
||||
/// Create an s-type Gaussian (l=m=n=0).
|
||||
#[must_use]
|
||||
pub fn s_type(alpha: f64, center: [f64; 3]) -> Self {
|
||||
Self {
|
||||
alpha,
|
||||
lmn: (0, 0, 0),
|
||||
center,
|
||||
}
|
||||
}
|
||||
|
||||
/// Create a p-type Gaussian (one of l,m,n = 1).
|
||||
#[must_use]
|
||||
pub fn px(alpha: f64, center: [f64; 3]) -> Self {
|
||||
Self {
|
||||
alpha,
|
||||
lmn: (1, 0, 0),
|
||||
center,
|
||||
}
|
||||
}
|
||||
|
||||
/// Angular momentum total L = l + m + n.
|
||||
#[must_use]
|
||||
pub fn angular_momentum(&self) -> u32 {
|
||||
self.lmn.0 + self.lmn.1 + self.lmn.2
|
||||
}
|
||||
|
||||
/// Normalization constant N for this Gaussian.
|
||||
#[must_use]
|
||||
pub fn normalization(&self) -> f64 {
|
||||
let (l, m, n) = self.lmn;
|
||||
let alpha = self.alpha;
|
||||
let lf = l as f64;
|
||||
let mf = m as f64;
|
||||
let nf = n as f64;
|
||||
let prefactor = (2.0 * alpha / PI).powf(0.75);
|
||||
// double_factorial(2k-1): for k=0 → 1 (convention), k=1 → 1, k=2 → 3, ...
|
||||
let df = |k: u32| -> f64 {
|
||||
if k == 0 {
|
||||
1.0
|
||||
} else {
|
||||
double_factorial(2 * k as i32 - 1) as f64
|
||||
}
|
||||
};
|
||||
let lm_norm = ((4.0 * alpha).powf(lf + mf + nf) / (df(l) * df(m) * df(n))).sqrt();
|
||||
prefactor * lm_norm
|
||||
}
|
||||
}
|
||||
|
||||
/// Double factorial: n!! = n·(n-2)·…·1 (or 2 for n=0).
|
||||
/// Returns 1 for n ≤ 0.
|
||||
fn double_factorial(n: i32) -> u64 {
|
||||
if n <= 0 {
|
||||
return 1;
|
||||
}
|
||||
if n > 33 {
|
||||
return u64::MAX;
|
||||
} // prevent overflow; values > 33!! exceed u64
|
||||
let mut result = 1u64;
|
||||
let mut k = n as u64;
|
||||
while k > 1 {
|
||||
result = result.saturating_mul(k);
|
||||
k = k.saturating_sub(2);
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Boys function F_m(x) = ∫₀¹ t^{2m} exp(-x t²) dt.
|
||||
/// Computed via series expansion for small x, asymptotic for large x.
|
||||
#[must_use]
|
||||
pub fn boys_function(m: u32, x: f64) -> f64 {
|
||||
if x < 1e-10 {
|
||||
// F_m(0) = 1 / (2m + 1)
|
||||
return 1.0 / (2 * m + 1) as f64;
|
||||
}
|
||||
if x > 20.0 {
|
||||
// Asymptotic: F_m(x) ≈ (2m-1)!! sqrt(π) / (2^{m+1} x^{m+1/2})
|
||||
let m = m as f64;
|
||||
let df = semi_factorial(m); // (2m-1)!!
|
||||
return df * PI.sqrt() / (2.0_f64.powf(m + 1.0) * x.powf(m + 0.5));
|
||||
}
|
||||
// Series: F_m(x) = e^{-x} Σ_{k=0}^∞ x^k / (2m+2k+1)!!
|
||||
let mut sum = 1.0 / (2 * m + 1) as f64; // k=0 term
|
||||
let mut x_pow = 1.0;
|
||||
let mut term;
|
||||
for k in 1..50usize {
|
||||
x_pow *= x;
|
||||
term = x_pow / double_factorial(2 * m as i32 + 2 * k as i32 + 1) as f64;
|
||||
sum += term;
|
||||
if term.abs() < 1e-15 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
(1.0 / (2 * m + 1) as f64 + sum) * (-x).exp()
|
||||
}
|
||||
|
||||
fn semi_factorial(m: f64) -> f64 {
|
||||
// (2m-1)!! for the asymptotic Boys formula: Γ(m+1/2)/Γ(1/2) = (2m-1)!!/2^m
|
||||
// Use Stirling / direct gamma approximation
|
||||
if m < 0.5 {
|
||||
return 1.0;
|
||||
}
|
||||
// Compute via recursion: (2m-1)!! = (2m-1) * (2m-3)!!
|
||||
let mi = m.round() as usize;
|
||||
if mi == 0 {
|
||||
return 1.0;
|
||||
}
|
||||
let mut result = 1.0_f64;
|
||||
for k in 0..mi {
|
||||
result *= (2 * k + 1) as f64;
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Overlap integral ⟨φ_a|φ_b⟩ between two Gaussian basis functions.
|
||||
///
|
||||
/// For s-type Gaussians: S(α,β,A,B) = (π/(α+β))^{3/2} exp(-αβ/(α+β)|A-B|²)
|
||||
#[must_use]
|
||||
pub fn overlap_integral(a: &GaussianBasis, b: &GaussianBasis) -> f64 {
|
||||
let alpha = a.alpha;
|
||||
let beta = b.alpha;
|
||||
let gamma = alpha + beta;
|
||||
|
||||
// Distance squared |A-B|²
|
||||
let rab2: f64 = a
|
||||
.center
|
||||
.iter()
|
||||
.zip(b.center.iter())
|
||||
.map(|(ai, bi)| (ai - bi).powi(2))
|
||||
.sum();
|
||||
|
||||
let prefactor = (PI / gamma).powf(1.5);
|
||||
let exponential = (-alpha * beta / gamma * rab2).exp();
|
||||
|
||||
// For s-type only (l=m=n=0); for higher angular momentum use recursion
|
||||
let angular_factor = if a.angular_momentum() == 0 && b.angular_momentum() == 0 {
|
||||
1.0
|
||||
} else {
|
||||
// McMurchie-Davidson E coefficient product
|
||||
let ex = mcmd_e(a.lmn.0, b.lmn.0, 0, a.center[0], b.center[0], alpha, beta);
|
||||
let ey = mcmd_e(a.lmn.1, b.lmn.1, 0, a.center[1], b.center[1], alpha, beta);
|
||||
let ez = mcmd_e(a.lmn.2, b.lmn.2, 0, a.center[2], b.center[2], alpha, beta);
|
||||
ex * ey * ez * (PI / gamma).powf(1.5) / prefactor
|
||||
};
|
||||
|
||||
let na = a.normalization();
|
||||
let nb = b.normalization();
|
||||
na * nb * prefactor * exponential * angular_factor
|
||||
}
|
||||
|
||||
/// McMurchie-Davidson E coefficient E^{ij}_{t}: expansion coefficient
|
||||
/// for the product of two 1D Gaussians.
|
||||
///
|
||||
/// Recursion:
|
||||
/// E^{i+1,j}_{t} = 1/(2γ) E^{ij}_{t-1} + (Px-Ax) E^{ij}_t + (t+1) E^{ij}_{t+1}
|
||||
/// E^{i,j+1}_{t} = 1/(2γ) E^{ij}_{t-1} + (Px-Bx) E^{ij}_t + (t+1) E^{ij}_{t+1}
|
||||
/// E^{00}_0 = exp(-μ X_AB²), E^{00}_{t≠0} = 0
|
||||
#[must_use]
|
||||
pub fn mcmd_e(i: u32, j: u32, t: u32, ax: f64, bx: f64, alpha: f64, beta: f64) -> f64 {
|
||||
let gamma = alpha + beta;
|
||||
let px = (alpha * ax + beta * bx) / gamma;
|
||||
let xab = ax - bx;
|
||||
let xpa = px - ax;
|
||||
let xpb = px - bx;
|
||||
let mu = alpha * beta / gamma;
|
||||
|
||||
// Memoize in a small table
|
||||
let max_i = i + j + t + 2;
|
||||
let size = (max_i + 1) as usize;
|
||||
let mut e = vec![vec![vec![0.0; size]; size]; size];
|
||||
e[0][0][0] = (-mu * xab * xab).exp();
|
||||
|
||||
for ii in 0..=(i + j) as usize {
|
||||
for jj in 0..=(i + j - ii as u32) as usize {
|
||||
for tt in 0..=(i + j) as usize {
|
||||
if ii == 0 && jj == 0 && tt == 0 {
|
||||
continue;
|
||||
}
|
||||
// Build up via recursion on j
|
||||
if jj > 0 {
|
||||
let prev = if tt > 0 {
|
||||
e[ii][jj - 1][tt - 1] / (2.0 * gamma)
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
let curr = xpb * e[ii][jj - 1][tt];
|
||||
let next = if tt + 1 < size {
|
||||
(tt as f64 + 1.0) * e[ii][jj - 1][tt + 1]
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
e[ii][jj][tt] = prev + curr + next;
|
||||
} else if ii > 0 {
|
||||
let prev = if tt > 0 {
|
||||
e[ii - 1][jj][tt - 1] / (2.0 * gamma)
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
let curr = xpa * e[ii - 1][jj][tt];
|
||||
let next = if tt + 1 < size {
|
||||
(tt as f64 + 1.0) * e[ii - 1][jj][tt + 1]
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
e[ii][jj][tt] = prev + curr + next;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
e[i as usize][j as usize][t as usize]
|
||||
}
|
||||
|
||||
/// Kinetic energy integral ⟨φ_a|T|φ_b⟩ where T = -½∇².
|
||||
///
|
||||
/// Uses the Obara-Saika relation applied to each Cartesian direction:
|
||||
/// T_x(i,j) = j(j-1)/2 · S(i,j-2) - α_b(2j+1) · S(i,j) + 2α_b² · S(i,j+2)
|
||||
/// Total: T = T_x · S_y · S_z + S_x · T_y · S_z + S_x · S_y · T_z
|
||||
///
|
||||
/// where S(i,j) is the 1D overlap integral between the i-th and j-th components.
|
||||
#[must_use]
|
||||
pub fn kinetic_integral(a: &GaussianBasis, b: &GaussianBasis) -> f64 {
|
||||
let alpha = a.alpha;
|
||||
let beta = b.alpha;
|
||||
let na = a.normalization();
|
||||
let nb = b.normalization();
|
||||
|
||||
// 1D overlap integrals for each Cartesian direction
|
||||
let s1d = |ia: u32, ib: u32, dim: usize| -> f64 {
|
||||
let ax = a.center[dim];
|
||||
let bx = b.center[dim];
|
||||
mcmd_e(ia, ib, 0, ax, bx, alpha, beta) * (PI / (alpha + beta)).sqrt()
|
||||
};
|
||||
|
||||
// 1D kinetic contribution in direction dim
|
||||
let t1d = |ia: u32, ib: u32, dim: usize| -> f64 {
|
||||
let bx = b.center[dim];
|
||||
let ax = a.center[dim];
|
||||
let ax2 = a.center[dim];
|
||||
let bx2 = b.center[dim];
|
||||
|
||||
// T(i,j) = j(j-1)/2 · S(i,j-2) - beta(2j+1) · S(i,j) + 2beta² · S(i,j+2)
|
||||
let term_low = if ib >= 2 {
|
||||
let s = mcmd_e(ia, ib - 2, 0, ax, bx, alpha, beta) * (PI / (alpha + beta)).sqrt();
|
||||
ib as f64 * (ib - 1) as f64 / 2.0 * s
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
|
||||
let term_mid = {
|
||||
let s = mcmd_e(ia, ib, 0, ax2, bx2, alpha, beta) * (PI / (alpha + beta)).sqrt();
|
||||
-beta * (2 * ib + 1) as f64 * s
|
||||
};
|
||||
|
||||
let term_high = {
|
||||
let s = mcmd_e(ia, ib + 2, 0, ax, bx, alpha, beta) * (PI / (alpha + beta)).sqrt();
|
||||
2.0 * beta * beta * s
|
||||
};
|
||||
|
||||
term_low + term_mid + term_high
|
||||
};
|
||||
|
||||
let (la, ma, na_ang) = a.lmn;
|
||||
let (lb, mb, nb_ang) = b.lmn;
|
||||
|
||||
// T = Tx·Sy·Sz + Sx·Ty·Sz + Sx·Sy·Tz
|
||||
let t = t1d(la, lb, 0) * s1d(ma, mb, 1) * s1d(na_ang, nb_ang, 2)
|
||||
+ s1d(la, lb, 0) * t1d(ma, mb, 1) * s1d(na_ang, nb_ang, 2)
|
||||
+ s1d(la, lb, 0) * s1d(ma, mb, 1) * t1d(na_ang, nb_ang, 2);
|
||||
|
||||
na * nb * t
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
|
||||
(a - b).abs() < tol
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn gaussian_s_orbital_normalization() {
|
||||
let g = GaussianBasis::s_type(1.0, [0.0, 0.0, 0.0]);
|
||||
let n = g.normalization();
|
||||
assert!(n > 0.0, "normalization should be positive: {n}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn overlap_ss_same_center() {
|
||||
// Two normalized s-Gaussians at same center: ⟨g|g⟩ = 1
|
||||
let g = GaussianBasis::s_type(1.0, [0.0, 0.0, 0.0]);
|
||||
let s = overlap_integral(&g, &g);
|
||||
assert!(approx_eq(s, 1.0, 1e-6), "⟨g|g⟩ = 1, got {s}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn overlap_ss_far_apart_is_small() {
|
||||
// Two s-Gaussians far apart: overlap → 0
|
||||
let g1 = GaussianBasis::s_type(1.0, [0.0, 0.0, 0.0]);
|
||||
let g2 = GaussianBasis::s_type(1.0, [100.0, 0.0, 0.0]);
|
||||
let s = overlap_integral(&g1, &g2);
|
||||
assert!(s.abs() < 1e-10, "overlap far apart should be ~0, got {s}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn boys_function_at_zero() {
|
||||
// F_0(0) = 1, F_1(0) = 1/3, F_2(0) = 1/5
|
||||
assert!(approx_eq(boys_function(0, 0.0), 1.0, 1e-10));
|
||||
assert!(approx_eq(boys_function(1, 0.0), 1.0 / 3.0, 1e-10));
|
||||
assert!(approx_eq(boys_function(2, 0.0), 1.0 / 5.0, 1e-10));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn boys_function_decreasing() {
|
||||
// F_m(x) should be positive and decreasing in x
|
||||
let f0 = boys_function(0, 1.0);
|
||||
let f1 = boys_function(0, 2.0);
|
||||
assert!(
|
||||
f0 > f1,
|
||||
"Boys function should decrease: F(1)={f0} > F(2)={f1}"
|
||||
);
|
||||
assert!(f0 > 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mcmurchie_davidson_base_case() {
|
||||
// E^{00}_0 = exp(-μ X²) with same center → = 1
|
||||
let e = mcmd_e(0, 0, 0, 0.0, 0.0, 1.0, 1.0);
|
||||
assert!(approx_eq(e, 1.0, 1e-10), "E^00_0(same center) = 1, got {e}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn kinetic_s_orbital() {
|
||||
// ⟨1s|T|1s⟩ = 3α/2 for a normalized s-type Gaussian with exponent α
|
||||
// For α=1.0: expected ≈ 1.5
|
||||
let g = GaussianBasis::s_type(1.0, [0.0, 0.0, 0.0]);
|
||||
let t = kinetic_integral(&g, &g);
|
||||
// Kinetic energy is always positive by definition
|
||||
assert!(
|
||||
t.abs() > 0.0,
|
||||
"kinetic energy integral should be non-zero, got {t}"
|
||||
);
|
||||
// For α=1, the exact value is 3/2 * α = 1.5
|
||||
assert!(
|
||||
(t.abs() - 1.5).abs() < 0.1,
|
||||
"⟨1s|T|1s⟩ ≈ 1.5 for α=1, got {t}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn double_factorial_values() {
|
||||
assert_eq!(double_factorial(-1), 1);
|
||||
assert_eq!(double_factorial(0), 1);
|
||||
assert_eq!(double_factorial(1), 1);
|
||||
assert_eq!(double_factorial(3), 3);
|
||||
assert_eq!(double_factorial(5), 15);
|
||||
assert_eq!(double_factorial(7), 105);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,17 @@
|
||||
//! # symclaw-qchem
|
||||
//!
|
||||
//! Quantum chemistry algebra for SymClaw.
|
||||
//!
|
||||
//! ## Modules
|
||||
//!
|
||||
//! - [`second_quant`] — Fermionic creation/annihilation operators and Wick's theorem
|
||||
//! - [`integrals`] — Gaussian basis functions, overlap/kinetic integrals (McMurchie-Davidson)
|
||||
//! - [`spin`] — Pauli matrices, Clebsch-Gordan coefficients, angular momentum
|
||||
//! - [`vqe`] — UCCSD ansatz, parameter-shift gradients, Hartree-Fock reference
|
||||
#![allow(missing_docs)]
|
||||
#![deny(unsafe_code)]
|
||||
|
||||
pub mod integrals;
|
||||
pub mod second_quant;
|
||||
pub mod spin;
|
||||
pub mod vqe;
|
||||
@@ -0,0 +1,7 @@
|
||||
//! Second quantization: creation/annihilation operators and Wick's theorem.
|
||||
|
||||
pub mod operators;
|
||||
pub mod wick;
|
||||
|
||||
pub use operators::{FermionOp, FermionTerm};
|
||||
pub use wick::{Contraction, NormalOrderedTerm, normal_order, vacuum_expectation, wick_expand};
|
||||
@@ -0,0 +1,238 @@
|
||||
//! Fermionic and bosonic creation/annihilation operators.
|
||||
//!
|
||||
//! Operators are represented symbolically as `Expr` trees using
|
||||
//! `FuncId::Create` and `FuncId::Annihilate`. The orbital index
|
||||
//! is the argument.
|
||||
|
||||
use serde::{Deserialize, Serialize};
|
||||
|
||||
// ── Operator types ────────────────────────────────────────────────
|
||||
|
||||
/// A fermionic operator term: coefficient × product of creation/annihilation ops.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct FermionTerm {
|
||||
/// Rational coefficient numerator.
|
||||
pub coeff_num: i64,
|
||||
/// Rational coefficient denominator.
|
||||
pub coeff_den: i64,
|
||||
/// Operator string in creation/annihilation sequence.
|
||||
/// Each element: (orbital_index, is_creation).
|
||||
pub ops: Vec<(usize, bool)>,
|
||||
}
|
||||
|
||||
impl FermionTerm {
|
||||
#[must_use]
|
||||
pub fn new(coeff_num: i64, coeff_den: i64, ops: Vec<(usize, bool)>) -> Self {
|
||||
Self {
|
||||
coeff_num,
|
||||
coeff_den,
|
||||
ops,
|
||||
}
|
||||
}
|
||||
|
||||
/// Creation operator a†_p.
|
||||
#[must_use]
|
||||
pub fn create(p: usize) -> Self {
|
||||
Self::new(1, 1, vec![(p, true)])
|
||||
}
|
||||
|
||||
/// Annihilation operator a_p.
|
||||
#[must_use]
|
||||
pub fn annihilate(p: usize) -> Self {
|
||||
Self::new(1, 1, vec![(p, false)])
|
||||
}
|
||||
|
||||
/// Coefficient as f64.
|
||||
#[must_use]
|
||||
pub fn coeff_f64(&self) -> f64 {
|
||||
self.coeff_num as f64 / self.coeff_den as f64
|
||||
}
|
||||
|
||||
/// Number of operators in this term.
|
||||
#[must_use]
|
||||
pub fn n_ops(&self) -> usize {
|
||||
self.ops.len()
|
||||
}
|
||||
|
||||
/// True if this is a number-conserving term (equal creation and annihilation ops).
|
||||
#[must_use]
|
||||
pub fn is_number_conserving(&self) -> bool {
|
||||
let creates = self.ops.iter().filter(|&&(_, c)| c).count();
|
||||
let annihilates = self.ops.iter().filter(|&&(_, c)| !c).count();
|
||||
creates == annihilates
|
||||
}
|
||||
}
|
||||
|
||||
impl std::fmt::Display for FermionTerm {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
if self.coeff_den == 1 {
|
||||
write!(f, "{}", self.coeff_num)?;
|
||||
} else {
|
||||
write!(f, "({}/{})", self.coeff_num, self.coeff_den)?;
|
||||
}
|
||||
for &(idx, is_create) in &self.ops {
|
||||
if is_create {
|
||||
write!(f, "·a†_{idx}")?;
|
||||
} else {
|
||||
write!(f, "·a_{idx}")?;
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
|
||||
/// A fermionic operator expression: sum of `FermionTerm`s.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct FermionOp {
|
||||
pub terms: Vec<FermionTerm>,
|
||||
}
|
||||
|
||||
impl FermionOp {
|
||||
/// Zero operator.
|
||||
#[must_use]
|
||||
pub fn zero() -> Self {
|
||||
Self { terms: vec![] }
|
||||
}
|
||||
|
||||
/// Single-term operator.
|
||||
#[must_use]
|
||||
pub fn from_term(t: FermionTerm) -> Self {
|
||||
Self { terms: vec![t] }
|
||||
}
|
||||
|
||||
/// a†_p.
|
||||
#[must_use]
|
||||
pub fn create(p: usize) -> Self {
|
||||
Self::from_term(FermionTerm::create(p))
|
||||
}
|
||||
|
||||
/// a_p.
|
||||
#[must_use]
|
||||
pub fn annihilate(p: usize) -> Self {
|
||||
Self::from_term(FermionTerm::annihilate(p))
|
||||
}
|
||||
|
||||
/// Add two operators.
|
||||
#[must_use]
|
||||
#[allow(clippy::should_implement_trait)]
|
||||
pub fn add(mut self, other: Self) -> Self {
|
||||
self.terms.extend(other.terms);
|
||||
self
|
||||
}
|
||||
|
||||
/// Scale by integer coefficient.
|
||||
#[must_use]
|
||||
pub fn scale(mut self, n: i64) -> Self {
|
||||
for t in &mut self.terms {
|
||||
t.coeff_num *= n;
|
||||
}
|
||||
self
|
||||
}
|
||||
|
||||
/// Multiply (concatenate operator strings).
|
||||
#[must_use]
|
||||
pub fn mul(&self, other: &Self) -> Self {
|
||||
let mut result = Self::zero();
|
||||
for a in &self.terms {
|
||||
for b in &other.terms {
|
||||
let mut ops = a.ops.clone();
|
||||
ops.extend(b.ops.iter().copied());
|
||||
let cn = a.coeff_num * b.coeff_num;
|
||||
let cd = a.coeff_den * b.coeff_den;
|
||||
let g = gcd(cn.abs(), cd.abs());
|
||||
result.terms.push(FermionTerm::new(cn / g, cd / g, ops));
|
||||
}
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Number of terms.
|
||||
#[must_use]
|
||||
pub fn n_terms(&self) -> usize {
|
||||
self.terms.len()
|
||||
}
|
||||
|
||||
/// True if all terms are number-conserving.
|
||||
#[must_use]
|
||||
pub fn is_number_conserving(&self) -> bool {
|
||||
self.terms.iter().all(|t| t.is_number_conserving())
|
||||
}
|
||||
}
|
||||
|
||||
impl std::fmt::Display for FermionOp {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
if self.terms.is_empty() {
|
||||
return write!(f, "0");
|
||||
}
|
||||
let parts: Vec<String> = self.terms.iter().map(|t| format!("{t}")).collect();
|
||||
write!(f, "{}", parts.join(" + "))
|
||||
}
|
||||
}
|
||||
|
||||
fn gcd(a: i64, b: i64) -> i64 {
|
||||
let (mut a, mut b) = (a.abs(), b.abs());
|
||||
while b != 0 {
|
||||
let t = b;
|
||||
b = a % b;
|
||||
a = t;
|
||||
}
|
||||
if a == 0 { 1 } else { a }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn creation_annihilation_display() {
|
||||
let c = FermionOp::create(0);
|
||||
let a = FermionOp::annihilate(0);
|
||||
assert!(format!("{c}").contains("a†_0"));
|
||||
assert!(format!("{a}").contains("a_0"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn number_conserving() {
|
||||
// a†_0 a_0 is number-conserving
|
||||
let op = FermionOp::create(0).mul(&FermionOp::annihilate(0));
|
||||
assert!(op.is_number_conserving());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn not_number_conserving() {
|
||||
// a†_0 alone is not
|
||||
assert!(!FermionOp::create(0).is_number_conserving());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mul_creates_product() {
|
||||
let c0 = FermionOp::create(0);
|
||||
let c1 = FermionOp::create(1);
|
||||
let prod = c0.mul(&c1);
|
||||
assert_eq!(prod.n_terms(), 1);
|
||||
assert_eq!(prod.terms[0].ops.len(), 2);
|
||||
assert_eq!(prod.terms[0].ops[0], (0, true));
|
||||
assert_eq!(prod.terms[0].ops[1], (1, true));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn add_operators() {
|
||||
let c0 = FermionOp::create(0);
|
||||
let c1 = FermionOp::create(1);
|
||||
let sum = c0.add(c1);
|
||||
assert_eq!(sum.n_terms(), 2);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn zero_operator() {
|
||||
let z = FermionOp::zero();
|
||||
assert_eq!(z.n_terms(), 0);
|
||||
assert_eq!(format!("{z}"), "0");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn coeff_arithmetic() {
|
||||
let t = FermionTerm::new(3, 4, vec![(0, true)]);
|
||||
assert!((t.coeff_f64() - 0.75).abs() < 1e-12);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,300 @@
|
||||
//! Wick's theorem: automated normal ordering of fermionic operator products.
|
||||
//!
|
||||
//! A product of creation/annihilation operators can be expressed as a sum of
|
||||
//! normal-ordered products (all a† left of all a) plus contraction terms.
|
||||
//!
|
||||
//! Vacuum contraction rule: ⟨0| a_p a†_q |0⟩ = δ_{pq} (only a then a†)
|
||||
//!
|
||||
//! Implementation: enumerate all subsets of valid contraction pairs via
|
||||
//! bitmask over the operator sequence, compute sign from bubble-sort,
|
||||
//! normal-order the remaining (uncontracted) operators.
|
||||
|
||||
use super::operators::FermionTerm;
|
||||
use serde::{Deserialize, Serialize};
|
||||
|
||||
/// A single contraction pair: ⟨a_{ann} a†_{cre}⟩ = δ_{ann,cre}.
|
||||
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
|
||||
pub struct Contraction {
|
||||
pub ann: usize, // annihilation orbital
|
||||
pub cre: usize, // creation orbital
|
||||
}
|
||||
|
||||
impl Contraction {
|
||||
#[must_use]
|
||||
pub fn new(ann: usize, cre: usize) -> Self {
|
||||
Self { ann, cre }
|
||||
}
|
||||
#[must_use]
|
||||
pub fn delta(&self) -> i64 {
|
||||
if self.ann == self.cre { 1 } else { 0 }
|
||||
}
|
||||
}
|
||||
|
||||
/// One term in the Wick expansion.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct NormalOrderedTerm {
|
||||
pub sign: i64,
|
||||
pub contractions: Vec<Contraction>,
|
||||
/// Remaining uncontracted operators, in normal order (all a† before a).
|
||||
pub ops: Vec<(usize, bool)>,
|
||||
}
|
||||
|
||||
impl NormalOrderedTerm {
|
||||
#[must_use]
|
||||
pub fn contraction_value(&self) -> i64 {
|
||||
self.contractions.iter().map(|c| c.delta()).product()
|
||||
}
|
||||
#[must_use]
|
||||
pub fn is_vacuum_expectation(&self) -> bool {
|
||||
self.ops.is_empty() && self.contraction_value() != 0
|
||||
}
|
||||
}
|
||||
|
||||
/// Wick-expand a sequence of operators `(orbital, is_creation)`.
|
||||
///
|
||||
/// Enumerates all valid contraction subsets (pairing each annihilation op
|
||||
/// with a later creation op) and collects all normal-ordered terms.
|
||||
///
|
||||
/// Time complexity: O(2^n · n) — practical for n ≤ 12.
|
||||
#[must_use]
|
||||
pub fn wick_expand(ops: &[(usize, bool)]) -> Vec<NormalOrderedTerm> {
|
||||
let n = ops.len();
|
||||
if n == 0 {
|
||||
return vec![NormalOrderedTerm {
|
||||
sign: 1,
|
||||
contractions: vec![],
|
||||
ops: vec![],
|
||||
}];
|
||||
}
|
||||
|
||||
// Collect all valid contraction pairs: (i, j) where ops[i] is annihilation,
|
||||
// ops[j] is creation, and i < j.
|
||||
let mut valid_pairs: Vec<(usize, usize)> = Vec::new();
|
||||
for i in 0..n {
|
||||
for j in (i + 1)..n {
|
||||
if !ops[i].1 && ops[j].1 {
|
||||
valid_pairs.push((i, j));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
let n_pairs = valid_pairs.len();
|
||||
let mut results = Vec::new();
|
||||
|
||||
// Enumerate all subsets of contraction pairs via bitmask
|
||||
let n_subsets = 1usize << n_pairs;
|
||||
for mask in 0..n_subsets {
|
||||
// Collect which pairs are contracted in this subset
|
||||
let selected_pairs: Vec<(usize, usize)> = (0..n_pairs)
|
||||
.filter(|&k| mask & (1 << k) != 0)
|
||||
.map(|k| valid_pairs[k])
|
||||
.collect();
|
||||
|
||||
// Check that each operator index appears at most once
|
||||
let mut used = vec![false; n];
|
||||
let mut valid = true;
|
||||
for &(i, j) in &selected_pairs {
|
||||
if used[i] || used[j] {
|
||||
valid = false;
|
||||
break;
|
||||
}
|
||||
used[i] = true;
|
||||
used[j] = true;
|
||||
}
|
||||
if !valid {
|
||||
continue;
|
||||
}
|
||||
|
||||
// Compute the sign: moving contracted ops adjacent, then eliminating pairs.
|
||||
// Sign from anticommutation: each pair (i,j) moved together requires
|
||||
// counting the number of unmatched ops between positions i and j.
|
||||
let sign = contraction_sign(ops, &selected_pairs);
|
||||
|
||||
// Build contraction list
|
||||
let contractions: Vec<Contraction> = selected_pairs
|
||||
.iter()
|
||||
.map(|&(i, j)| Contraction::new(ops[i].0, ops[j].0))
|
||||
.collect();
|
||||
|
||||
// Remaining uncontracted operators (in original order)
|
||||
let remaining: Vec<(usize, bool)> = (0..n).filter(|k| !used[*k]).map(|k| ops[k]).collect();
|
||||
|
||||
// Normal-order the remaining operators (all a† before a)
|
||||
let (ordered, order_sign) = normal_order_ops(&remaining);
|
||||
|
||||
results.push(NormalOrderedTerm {
|
||||
sign: sign * order_sign,
|
||||
contractions,
|
||||
ops: ordered,
|
||||
});
|
||||
}
|
||||
|
||||
results
|
||||
}
|
||||
|
||||
/// Compute the sign from moving contracted pairs out of the operator sequence.
|
||||
///
|
||||
/// For each contraction pair (i, j): moving op[j] to be adjacent to op[i]
|
||||
/// requires passing over all ops between i+1 and j-1 that are NOT already
|
||||
/// contracted in this subset. Each such pass flips the sign.
|
||||
fn contraction_sign(ops: &[(usize, bool)], pairs: &[(usize, usize)]) -> i64 {
|
||||
let n = ops.len();
|
||||
// Build a boolean array: contracted[k] = true if k is in a pair
|
||||
let mut contracted = vec![false; n];
|
||||
for &(i, j) in pairs {
|
||||
contracted[i] = true;
|
||||
contracted[j] = true;
|
||||
}
|
||||
|
||||
let mut sign = 1i64;
|
||||
// Process pairs in order of first index
|
||||
let mut sorted_pairs = pairs.to_vec();
|
||||
sorted_pairs.sort_by_key(|&(i, _)| i);
|
||||
|
||||
// For each pair (i,j), count uncontracted ops strictly between i and j
|
||||
// at the time of contraction (ops already extracted don't count)
|
||||
let mut extracted = vec![false; n];
|
||||
for &(i, j) in &sorted_pairs {
|
||||
// Count ops between i+1 and j-1 that haven't been extracted yet
|
||||
let between = ((i + 1)..j).filter(|&k| !extracted[k]).count();
|
||||
if between % 2 != 0 {
|
||||
sign = -sign;
|
||||
}
|
||||
extracted[i] = true;
|
||||
extracted[j] = true;
|
||||
}
|
||||
sign
|
||||
}
|
||||
|
||||
/// Sort operators into normal order (all a† before a) using bubble sort.
|
||||
/// Returns `(sorted_ops, sign)`.
|
||||
fn normal_order_ops(ops: &[(usize, bool)]) -> (Vec<(usize, bool)>, i64) {
|
||||
let mut sorted = ops.to_vec();
|
||||
let mut sign = 1i64;
|
||||
let n = sorted.len();
|
||||
// Stable bubble sort: move creation ops left
|
||||
for i in 0..n {
|
||||
for j in 0..(n.saturating_sub(i + 1)) {
|
||||
// If annihilation (false) is followed by creation (true): swap
|
||||
if !sorted[j].1 && sorted[j + 1].1 {
|
||||
sorted.swap(j, j + 1);
|
||||
sign = -sign;
|
||||
}
|
||||
}
|
||||
}
|
||||
(sorted, sign)
|
||||
}
|
||||
|
||||
/// Vacuum expectation value ⟨0|ops|0⟩ via Wick's theorem.
|
||||
#[must_use]
|
||||
pub fn vacuum_expectation(ops: &[(usize, bool)]) -> i64 {
|
||||
if !ops.len().is_multiple_of(2) {
|
||||
return 0;
|
||||
}
|
||||
wick_expand(ops)
|
||||
.iter()
|
||||
.filter(|t| t.is_vacuum_expectation())
|
||||
.map(|t| t.sign * t.contraction_value())
|
||||
.sum()
|
||||
}
|
||||
|
||||
/// Normal-order a single `FermionTerm`.
|
||||
#[must_use]
|
||||
pub fn normal_order(term: &FermionTerm) -> Vec<NormalOrderedTerm> {
|
||||
wick_expand(&term.ops)
|
||||
.into_iter()
|
||||
.map(|mut t| {
|
||||
t.sign *= term.coeff_num;
|
||||
t
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn anticommutation_relation_same_index() {
|
||||
// ⟨0|a_0 a†_0|0⟩ = 1
|
||||
let ops = vec![(0, false), (0, true)];
|
||||
assert_eq!(vacuum_expectation(&ops), 1, "⟨0|a_0 a†_0|0⟩ = 1");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn anticommutation_relation_diff_index() {
|
||||
// ⟨0|a_0 a†_1|0⟩ = δ_{0,1} = 0
|
||||
let ops = vec![(0, false), (1, true)];
|
||||
assert_eq!(vacuum_expectation(&ops), 0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn normal_order_single_creation() {
|
||||
let ops = vec![(0usize, true)];
|
||||
let terms = wick_expand(&ops);
|
||||
assert!(!terms.is_empty());
|
||||
assert!(terms.iter().any(|t| t.ops == vec![(0, true)]));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn normal_order_single_product() {
|
||||
// a†_0 a_0 is already normal-ordered
|
||||
let ops = vec![(0, true), (0, false)];
|
||||
let terms = wick_expand(&ops);
|
||||
assert!(!terms.is_empty());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wick_theorem_one_body() {
|
||||
// ⟨0|a†_0 a_0|0⟩ = 0 (no annihilation before creation = no valid contraction)
|
||||
let ops = vec![(0, true), (0, false)];
|
||||
assert_eq!(vacuum_expectation(&ops), 0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wick_theorem_two_body_vacuum() {
|
||||
// ⟨0|a_0 a_1 a†_1 a†_0|0⟩
|
||||
let ops = vec![(0, false), (1, false), (1, true), (0, true)];
|
||||
let vev = vacuum_expectation(&ops);
|
||||
// = ⟨a_0 a†_0⟩⟨a_1 a†_1⟩ sign... should be ±1
|
||||
assert!(vev.abs() <= 1, "two-body VEV should be ±1: {vev}");
|
||||
assert_ne!(
|
||||
vev, 0,
|
||||
"two-body VEV should not be zero for matching indices"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn odd_ops_vev_zero() {
|
||||
let ops = vec![(0, false), (0, true), (1, true)];
|
||||
assert_eq!(vacuum_expectation(&ops), 0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn normal_order_two_body() {
|
||||
let ops = vec![(0, false), (0, true), (1, false), (1, true)];
|
||||
let terms = wick_expand(&ops);
|
||||
assert!(!terms.is_empty());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wick_expand_empty() {
|
||||
let terms = wick_expand(&[]);
|
||||
assert_eq!(terms.len(), 1);
|
||||
assert!(terms[0].ops.is_empty());
|
||||
assert_eq!(terms[0].sign, 1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn normal_order_two_ops_already_normal() {
|
||||
// a†_0 a_1: already normal — sign should be +1
|
||||
let ops = vec![(0, true), (1, false)];
|
||||
let terms = wick_expand(&ops);
|
||||
// No valid contractions (creation before annihilation in original order)
|
||||
// so exactly one term with the ops as-is
|
||||
let no_contraction = terms.iter().find(|t| t.contractions.is_empty());
|
||||
assert!(no_contraction.is_some());
|
||||
let t = no_contraction.unwrap();
|
||||
assert_eq!(t.sign, 1);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,200 @@
|
||||
//! Angular momentum, spin operators, and Clebsch-Gordan coefficients.
|
||||
|
||||
/// Spin-1/2 Pauli matrices as [[complex_re, complex_im]; 2×2].
|
||||
/// Format: `matrix[row][col] = (re, im)`.
|
||||
#[must_use]
|
||||
pub fn pauli_x() -> [[(f64, f64); 2]; 2] {
|
||||
[[(0.0, 0.0), (1.0, 0.0)], [(1.0, 0.0), (0.0, 0.0)]]
|
||||
}
|
||||
#[must_use]
|
||||
pub fn pauli_y() -> [[(f64, f64); 2]; 2] {
|
||||
[[(0.0, 0.0), (0.0, -1.0)], [(0.0, 1.0), (0.0, 0.0)]]
|
||||
}
|
||||
#[must_use]
|
||||
pub fn pauli_z() -> [[(f64, f64); 2]; 2] {
|
||||
[[(1.0, 0.0), (0.0, 0.0)], [(0.0, 0.0), (-1.0, 0.0)]]
|
||||
}
|
||||
|
||||
/// Check [A, B] = AB - BA = iC for Pauli algebra: [σx,σy] = 2i·σz.
|
||||
#[must_use]
|
||||
pub fn pauli_commutator_xy_is_i_sigmaz() -> bool {
|
||||
// [σx, σy] = 2i σz → result should be 2i·σz
|
||||
// σx·σy[0][0] = 0*0 + 1*i = i
|
||||
// σy·σx[0][0] = 0*0 + (-i)*1 = -i
|
||||
// commutator[0][0] = i - (-i) = 2i
|
||||
// 2i·σz[0][0] = 2i*1 = 2i ✓
|
||||
let sx = pauli_x();
|
||||
let sy = pauli_y();
|
||||
|
||||
// Compute σx·σy
|
||||
let xy_00_re = sx[0][0].0 * sy[0][0].0 - sx[0][0].1 * sy[0][0].1 + sx[0][1].0 * sy[1][0].0
|
||||
- sx[0][1].1 * sy[1][0].1;
|
||||
let xy_00_im = sx[0][0].0 * sy[0][0].1
|
||||
+ sx[0][0].1 * sy[0][0].0
|
||||
+ sx[0][1].0 * sy[1][0].1
|
||||
+ sx[0][1].1 * sy[1][0].0;
|
||||
|
||||
// Compute σy·σx
|
||||
let yx_00_re = sy[0][0].0 * sx[0][0].0 - sy[0][0].1 * sx[0][0].1 + sy[0][1].0 * sx[1][0].0
|
||||
- sy[0][1].1 * sx[1][0].1;
|
||||
let yx_00_im = sy[0][0].0 * sx[0][0].1
|
||||
+ sy[0][0].1 * sx[0][0].0
|
||||
+ sy[0][1].0 * sx[1][0].1
|
||||
+ sy[0][1].1 * sx[1][0].0;
|
||||
|
||||
// [σx,σy][0][0] should be 2i (re=0, im=2)
|
||||
let comm_re = xy_00_re - yx_00_re;
|
||||
let comm_im = xy_00_im - yx_00_im;
|
||||
comm_re.abs() < 1e-12 && (comm_im - 2.0).abs() < 1e-12
|
||||
}
|
||||
|
||||
/// Clebsch-Gordan coefficient ⟨j1,m1; j2,m2 | J,M⟩.
|
||||
///
|
||||
/// Arguments are all doubled to avoid half-integers:
|
||||
/// `j1_2` = 2*j1, `m1_2` = 2*m1, etc.
|
||||
///
|
||||
/// Uses the Racah formula.
|
||||
#[must_use]
|
||||
pub fn clebsch_gordan(j1_2: i32, m1_2: i32, j2_2: i32, m2_2: i32, j_2: i32, m_2: i32) -> f64 {
|
||||
// Selection rules
|
||||
if m1_2 + m2_2 != m_2 {
|
||||
return 0.0;
|
||||
}
|
||||
if (j_2 - j1_2 - j2_2).abs() > 0 && j_2 < (j1_2 - j2_2).abs() {
|
||||
return 0.0;
|
||||
}
|
||||
if j_2 < 0 || j_2 > j1_2 + j2_2 {
|
||||
return 0.0;
|
||||
}
|
||||
if m1_2.abs() > j1_2 || m2_2.abs() > j2_2 || m_2.abs() > j_2 {
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
// Convert to actual values for computation
|
||||
let j1 = j1_2 as f64 / 2.0;
|
||||
let m1 = m1_2 as f64 / 2.0;
|
||||
let j2 = j2_2 as f64 / 2.0;
|
||||
let m2 = m2_2 as f64 / 2.0;
|
||||
let j = j_2 as f64 / 2.0;
|
||||
let m = m_2 as f64 / 2.0;
|
||||
|
||||
// Racah formula
|
||||
let delta = delta_factor(j1, j2, j);
|
||||
if delta.abs() < 1e-15 {
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
let prefactor = delta
|
||||
* ((2.0 * j + 1.0)
|
||||
* factorial(j + m)
|
||||
* factorial(j - m)
|
||||
* factorial(j1 + m1)
|
||||
* factorial(j1 - m1)
|
||||
* factorial(j2 + m2)
|
||||
* factorial(j2 - m2))
|
||||
.sqrt();
|
||||
|
||||
// Sum over s
|
||||
let s_min = 0i32;
|
||||
let s_max = 20i32; // sufficient for reasonable j values
|
||||
let mut sum = 0.0;
|
||||
for s in s_min..=s_max {
|
||||
let sf = s as f64;
|
||||
let d1 = j1 + j2 - j - sf;
|
||||
let d2 = j1 - m1 - sf;
|
||||
let d3 = j2 + m2 - sf;
|
||||
let d4 = j - j2 + m1 + sf;
|
||||
let d5 = j - j1 - m2 + sf;
|
||||
if d1 < 0.0 || d2 < 0.0 || d3 < 0.0 || d4 < 0.0 || d5 < 0.0 {
|
||||
continue;
|
||||
}
|
||||
let sign = if s % 2 == 0 { 1.0 } else { -1.0 };
|
||||
let denom = factorial(sf)
|
||||
* factorial(d1)
|
||||
* factorial(d2)
|
||||
* factorial(d3)
|
||||
* factorial(d4)
|
||||
* factorial(d5);
|
||||
if denom.abs() < 1e-15 {
|
||||
continue;
|
||||
}
|
||||
sum += sign / denom;
|
||||
}
|
||||
|
||||
prefactor * sum
|
||||
}
|
||||
|
||||
fn delta_factor(j1: f64, j2: f64, j: f64) -> f64 {
|
||||
(factorial(j1 + j2 - j) * factorial(j1 - j2 + j) * factorial(-j1 + j2 + j)
|
||||
/ factorial(j1 + j2 + j + 1.0))
|
||||
.sqrt()
|
||||
}
|
||||
|
||||
fn factorial(n: f64) -> f64 {
|
||||
if n < 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
let n = n.round() as u64;
|
||||
(1..=n).product::<u64>() as f64
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn approx(a: f64, b: f64) -> bool {
|
||||
(a - b).abs() < 1e-8
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn pauli_matrices_from_spinors() {
|
||||
// σx† = σx (Hermitian)
|
||||
let sx = pauli_x();
|
||||
// For 2×2: Hermitian means sx[i][j] = conj(sx[j][i])
|
||||
// sx[0][1] = (1,0) = conj(sx[1][0]) = conj((1,0)) = (1,0) ✓
|
||||
assert_eq!(sx[0][1].0, sx[1][0].0);
|
||||
assert_eq!(sx[0][1].1, -sx[1][0].1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn angular_momentum_algebra() {
|
||||
// [Jx, Jy] = i Jz → [σx/2, σy/2] = i σz/2 → [σx, σy] = 2i σz
|
||||
assert!(pauli_commutator_xy_is_i_sigmaz(), "[σx,σy] ≠ 2i·σz");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn clebsch_gordan_half_half() {
|
||||
// 1/2 ⊗ 1/2 = 0 ⊕ 1
|
||||
// ⟨1/2,1/2; 1/2,-1/2 | 0,0⟩ = 1/√2
|
||||
let cg = clebsch_gordan(1, 1, 1, -1, 0, 0);
|
||||
assert!(
|
||||
approx(cg, 1.0 / 2.0_f64.sqrt()),
|
||||
"CG(1/2,1/2;1/2,-1/2|0,0) = 1/√2, got {cg}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn clebsch_gordan_triplet_m1() {
|
||||
// ⟨1/2,1/2; 1/2,1/2 | 1,1⟩ = 1
|
||||
let cg = clebsch_gordan(1, 1, 1, 1, 2, 2);
|
||||
assert!(approx(cg, 1.0), "CG(1/2,+1/2;1/2,+1/2|1,1) = 1, got {cg}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn clebsch_gordan_selection_rules() {
|
||||
// m1 + m2 ≠ M → 0
|
||||
let cg = clebsch_gordan(1, 1, 1, 1, 1, -1);
|
||||
assert!(approx(cg, 0.0), "violated selection rule should give 0");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn clebsch_gordan_orthogonality() {
|
||||
// Σ_M |⟨j1,m1; j2,m2|J,M⟩|² should sum correctly
|
||||
// For simplicity: check ⟨1/2,1/2;1/2,-1/2|1,0⟩ = 1/√2
|
||||
let cg = clebsch_gordan(1, 1, 1, -1, 2, 0);
|
||||
assert!(
|
||||
approx(cg, 1.0 / 2.0_f64.sqrt()),
|
||||
"CG(1/2,1/2;1/2,-1/2|1,0) = 1/√2, got {cg}"
|
||||
);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,261 @@
|
||||
//! VQE (Variational Quantum Eigensolver) tools.
|
||||
//!
|
||||
//! Provides:
|
||||
//! - UCCSD ansatz parameter-shift gradient formula
|
||||
//! - Hartree-Fock reference state (symbolic)
|
||||
//! - Parameter-shift rule for quantum gradients
|
||||
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::collections::HashMap;
|
||||
|
||||
/// A parameterized quantum ansatz: list of (gate_type, qubit_indices, parameter_name).
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct Ansatz {
|
||||
/// Layers: each is a list of (gate, qubits, param_name).
|
||||
pub layers: Vec<AnsatzLayer>,
|
||||
/// Current parameter values.
|
||||
pub params: HashMap<String, f64>,
|
||||
/// Number of qubits.
|
||||
pub n_qubits: usize,
|
||||
}
|
||||
|
||||
/// A single ansatz layer.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct AnsatzLayer {
|
||||
pub gates: Vec<AnsatzGate>,
|
||||
}
|
||||
|
||||
/// A parameterized gate in the ansatz.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct AnsatzGate {
|
||||
pub gate_type: String,
|
||||
pub qubits: Vec<usize>,
|
||||
pub param_name: Option<String>,
|
||||
}
|
||||
|
||||
impl Ansatz {
|
||||
/// Create an empty ansatz for n qubits.
|
||||
#[must_use]
|
||||
pub fn new(n_qubits: usize) -> Self {
|
||||
Self {
|
||||
layers: Vec::new(),
|
||||
n_qubits,
|
||||
params: HashMap::new(),
|
||||
}
|
||||
}
|
||||
|
||||
/// Add a Rz(θ) rotation layer.
|
||||
pub fn add_rz_layer(&mut self, qubit: usize, param_name: &str, init_val: f64) -> &mut Self {
|
||||
self.layers.push(AnsatzLayer {
|
||||
gates: vec![AnsatzGate {
|
||||
gate_type: "Rz".into(),
|
||||
qubits: vec![qubit],
|
||||
param_name: Some(param_name.into()),
|
||||
}],
|
||||
});
|
||||
self.params.insert(param_name.into(), init_val);
|
||||
self
|
||||
}
|
||||
|
||||
/// Add a CNOT entangling layer.
|
||||
pub fn add_cnot_layer(&mut self, ctrl: usize, tgt: usize) -> &mut Self {
|
||||
self.layers.push(AnsatzLayer {
|
||||
gates: vec![AnsatzGate {
|
||||
gate_type: "CNOT".into(),
|
||||
qubits: vec![ctrl, tgt],
|
||||
param_name: None,
|
||||
}],
|
||||
});
|
||||
self
|
||||
}
|
||||
|
||||
/// Number of variational parameters.
|
||||
#[must_use]
|
||||
pub fn n_params(&self) -> usize {
|
||||
self.params.len()
|
||||
}
|
||||
|
||||
/// List parameter names.
|
||||
#[must_use]
|
||||
pub fn param_names(&self) -> Vec<&str> {
|
||||
self.params.keys().map(String::as_str).collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// Parameter-shift gradient rule.
|
||||
///
|
||||
/// For a gate G(θ) = exp(-iθ/2 P) where P is a Pauli:
|
||||
/// ∂⟨H⟩/∂θ = [⟨H⟩(θ + π/2) - ⟨H⟩(θ - π/2)] / 2
|
||||
///
|
||||
/// This function returns `(θ_plus, θ_minus, coefficient)` for a parameter at `θ`.
|
||||
#[must_use]
|
||||
pub fn parameter_shift_gradient_rule(theta: f64) -> (f64, f64, f64) {
|
||||
let shift = std::f64::consts::PI / 2.0;
|
||||
(theta + shift, theta - shift, 0.5)
|
||||
}
|
||||
|
||||
/// Compute the numerical gradient of an expectation value using parameter-shift rule.
|
||||
///
|
||||
/// `expectation_fn` takes a parameter map and returns ⟨H⟩.
|
||||
pub fn parameter_shift_gradient<F>(
|
||||
params: &HashMap<String, f64>,
|
||||
param_name: &str,
|
||||
expectation_fn: &F,
|
||||
) -> f64
|
||||
where
|
||||
F: Fn(&HashMap<String, f64>) -> f64,
|
||||
{
|
||||
let theta = *params.get(param_name).unwrap_or(&0.0);
|
||||
let (tp, tm, coeff) = parameter_shift_gradient_rule(theta);
|
||||
|
||||
let mut params_plus = params.clone();
|
||||
params_plus.insert(param_name.into(), tp);
|
||||
let mut params_minus = params.clone();
|
||||
params_minus.insert(param_name.into(), tm);
|
||||
|
||||
coeff * (expectation_fn(¶ms_plus) - expectation_fn(¶ms_minus))
|
||||
}
|
||||
|
||||
/// Hartree-Fock reference state description.
|
||||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||||
pub struct HFReference {
|
||||
/// Number of electrons.
|
||||
pub n_electrons: usize,
|
||||
/// Number of spin-orbitals.
|
||||
pub n_orbitals: usize,
|
||||
/// Occupied orbital indices (0-indexed).
|
||||
pub occupied: Vec<usize>,
|
||||
/// Virtual orbital indices.
|
||||
pub virtual_orbs: Vec<usize>,
|
||||
}
|
||||
|
||||
impl HFReference {
|
||||
/// Create HF reference for `n_electrons` in `n_orbitals` spin-orbitals.
|
||||
#[must_use]
|
||||
pub fn new(n_electrons: usize, n_orbitals: usize) -> Self {
|
||||
assert!(n_electrons <= n_orbitals);
|
||||
let occupied: Vec<usize> = (0..n_electrons).collect();
|
||||
let virtual_orbs: Vec<usize> = (n_electrons..n_orbitals).collect();
|
||||
Self {
|
||||
n_electrons,
|
||||
n_orbitals,
|
||||
occupied,
|
||||
virtual_orbs,
|
||||
}
|
||||
}
|
||||
|
||||
/// Number of occupied orbitals.
|
||||
#[must_use]
|
||||
pub fn n_occ(&self) -> usize {
|
||||
self.occupied.len()
|
||||
}
|
||||
|
||||
/// Number of virtual orbitals.
|
||||
#[must_use]
|
||||
pub fn n_virt(&self) -> usize {
|
||||
self.virtual_orbs.len()
|
||||
}
|
||||
|
||||
/// UCCSD single excitation operators: (i→a) for i in occ, a in virt.
|
||||
#[must_use]
|
||||
pub fn single_excitations(&self) -> Vec<(usize, usize)> {
|
||||
let mut excitations = Vec::new();
|
||||
for &i in &self.occupied {
|
||||
for &a in &self.virtual_orbs {
|
||||
excitations.push((i, a));
|
||||
}
|
||||
}
|
||||
excitations
|
||||
}
|
||||
|
||||
/// UCCSD double excitation operators: (i,j→a,b) for i<j in occ, a<b in virt.
|
||||
#[must_use]
|
||||
pub fn double_excitations(&self) -> Vec<(usize, usize, usize, usize)> {
|
||||
let mut excitations = Vec::new();
|
||||
for (ki, &i) in self.occupied.iter().enumerate() {
|
||||
for &j in &self.occupied[ki + 1..] {
|
||||
for (ka, &a) in self.virtual_orbs.iter().enumerate() {
|
||||
for &b in &self.virtual_orbs[ka + 1..] {
|
||||
excitations.push((i, j, a, b));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
excitations
|
||||
}
|
||||
|
||||
/// Total UCCSD parameters: singles + doubles.
|
||||
#[must_use]
|
||||
pub fn n_uccsd_params(&self) -> usize {
|
||||
self.single_excitations().len() + self.double_excitations().len()
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use std::f64::consts::PI;
|
||||
|
||||
#[test]
|
||||
fn parameter_shift_gradient_rule() {
|
||||
let theta = 0.5;
|
||||
let (tp, tm, coeff) = super::parameter_shift_gradient_rule(theta);
|
||||
assert!((tp - (theta + PI / 2.0)).abs() < 1e-12);
|
||||
assert!((tm - (theta - PI / 2.0)).abs() < 1e-12);
|
||||
assert!((coeff - 0.5).abs() < 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn parameter_shift_numerical_gradient() {
|
||||
// f(θ) = cos(θ), ∂f/∂θ = -sin(θ)
|
||||
let theta = 0.7;
|
||||
let mut params = HashMap::new();
|
||||
params.insert("theta".into(), theta);
|
||||
let grad = parameter_shift_gradient(¶ms, "theta", &|p| p["theta"].cos());
|
||||
let expected = -theta.sin();
|
||||
assert!(
|
||||
(grad - expected).abs() < 1e-6,
|
||||
"gradient {grad} ≠ expected {expected}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hf_reference_state() {
|
||||
let hf = HFReference::new(2, 4);
|
||||
assert_eq!(hf.n_occ(), 2);
|
||||
assert_eq!(hf.n_virt(), 2);
|
||||
assert_eq!(hf.occupied, vec![0, 1]);
|
||||
assert_eq!(hf.virtual_orbs, vec![2, 3]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn uccsd_ansatz_singles() {
|
||||
// 2 occ × 2 virt = 4 singles
|
||||
let hf = HFReference::new(2, 4);
|
||||
assert_eq!(hf.single_excitations().len(), 4);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn uccsd_ansatz_doubles() {
|
||||
// C(2,2) × C(2,2) = 1 double excitation for 2 occ, 2 virt
|
||||
let hf = HFReference::new(2, 4);
|
||||
assert_eq!(hf.double_excitations().len(), 1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn uccsd_parameterised() {
|
||||
let hf = HFReference::new(2, 4);
|
||||
let n = hf.n_uccsd_params();
|
||||
assert_eq!(n, 5, "2e/4o: 4 singles + 1 double = 5 UCCSD params");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ansatz_construction() {
|
||||
let mut ansatz = Ansatz::new(2);
|
||||
ansatz.add_rz_layer(0, "theta0", 0.0);
|
||||
ansatz.add_cnot_layer(0, 1);
|
||||
ansatz.add_rz_layer(1, "theta1", 0.0);
|
||||
assert_eq!(ansatz.n_params(), 2);
|
||||
assert_eq!(ansatz.layers.len(), 3);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user