diff --git a/Cargo.lock b/Cargo.lock index 3aeeede49..2af0cdbc0 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -228,6 +228,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "816825ba1726f300e5558f701bf3be2d2d2d2b6bdd0e648f338366c9ce12798a" dependencies = [ "approx-derive", + "num-complex", "num-traits", ] @@ -450,7 +451,7 @@ dependencies = [ "hoomd-rand", "hoomd-simulation", "hoomd-spatial", - "hoomd-vector", + "hoomd-vector 1.1.0", "log", "parquet", "parquet_derive", @@ -2816,8 +2817,8 @@ dependencies = [ "hoomd-microstate", "hoomd-simulation", "hoomd-spatial", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "parquet", "parquet_derive", "rand 0.10.1", @@ -3633,7 +3634,7 @@ dependencies = [ "hoomd-manifold", "hoomd-microstate", "hoomd-simulation", - "hoomd-vector", + "hoomd-vector 1.1.0", "itertools", "wayland-sys", ] @@ -3656,10 +3657,10 @@ dependencies = [ "arrayvec", "assert2", "divan", - "hoomd-linear-algebra", + "hoomd-linear-algebra 1.1.0", "hoomd-manifold", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "itertools", "rand 0.10.1", "robust", @@ -3676,7 +3677,7 @@ dependencies = [ "anyhow", "assert2", "clap", - "hoomd-vector", + "hoomd-vector 1.1.0", "itertools", "memmap2", "tempfile", @@ -3695,8 +3696,8 @@ dependencies = [ "hoomd-geometry", "hoomd-microstate", "hoomd-spatial", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "rand 0.10.1", "rstest", "serde", @@ -3716,14 +3717,24 @@ dependencies = [ "serde_with", ] +[[package]] +name = "hoomd-linear-algebra" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "71cdc3ecd5a9301198375544d900e83603effd4608adb7e2fe50cf8a4cec2d46" +dependencies = [ + "serde", + "serde_with", +] + [[package]] name = "hoomd-manifold" version = "1.1.0" dependencies = [ "approxim", "divan", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "num", "paste", "rand 0.10.1", @@ -3747,8 +3758,8 @@ dependencies = [ "hoomd-microstate", "hoomd-simulation", "hoomd-spatial", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "itertools", "log", "rand 0.10.1", @@ -3774,8 +3785,8 @@ dependencies = [ "hoomd-manifold", "hoomd-rand", "hoomd-spatial", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "itertools", "log", "rand 0.10.1", @@ -3786,6 +3797,19 @@ dependencies = [ "thiserror 2.0.18", ] +[[package]] +name = "hoomd-order" +version = "1.1.0" +dependencies = [ + "approxim", + "divan", + "hoomd-vector 1.1.0 (registry+https://github.com/rust-lang/crates.io-index)", + "num-complex", + "rand 0.10.1", + "rstest", + "sphrs", +] + [[package]] name = "hoomd-rand" version = "1.1.0" @@ -3814,8 +3838,8 @@ name = "hoomd-spatial" version = "1.1.0" dependencies = [ "assert2", - "hoomd-utility", - "hoomd-vector", + "hoomd-utility 1.1.0", + "hoomd-vector 1.1.0", "log", "rand 0.10.1", "rstest", @@ -3837,6 +3861,18 @@ dependencies = [ "thiserror 2.0.18", ] +[[package]] +name = "hoomd-utility" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9174f96aa110e76dd39c1d61cea1d414ba59d5f0ad3c86a58a4aeb82e552201b" +dependencies = [ + "parquet", + "parquet_derive", + "serde", + "thiserror 2.0.18", +] + [[package]] name = "hoomd-vector" version = "1.1.0" @@ -3845,9 +3881,9 @@ dependencies = [ "approxim", "assert2", "divan", - "hoomd-linear-algebra", + "hoomd-linear-algebra 1.1.0", "hoomd-rand", - "hoomd-utility", + "hoomd-utility 1.1.0", "paste", "rand 0.10.1", "rand_distr 0.6.0", @@ -3857,6 +3893,22 @@ dependencies = [ "thiserror 2.0.18", ] +[[package]] +name = "hoomd-vector" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fbb457e90d8c763336856a2d3b0e8c69219e96b71be92ebafd4a843f8bfd993d" +dependencies = [ + "approxim", + "hoomd-linear-algebra 1.1.0 (registry+https://github.com/rust-lang/crates.io-index)", + "hoomd-utility 1.1.0 (registry+https://github.com/rust-lang/crates.io-index)", + "rand 0.10.1", + "rand_distr 0.6.0", + "serde", + "serde_with", + "thiserror 2.0.18", +] + [[package]] name = "iana-time-zone" version = "0.1.65" @@ -3869,7 +3921,7 @@ dependencies = [ "js-sys", "log", "wasm-bindgen", - "windows-core 0.58.0", + "windows-core 0.61.2", ] [[package]] @@ -6108,6 +6160,17 @@ dependencies = [ "serde", ] +[[package]] +name = "sphrs" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1908b37a136e29c075b9e46581fb4a1b9ff7feaacc5cf6100fabb9421927da1" +dependencies = [ + "num", + "num-complex", + "num-traits", +] + [[package]] name = "spin" version = "0.10.0" diff --git a/Cargo.toml b/Cargo.toml index 28129e7f7..ebf4189de 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -18,7 +18,7 @@ members = ["benchmarks", "hoomd-vector", "tools/build-wasm-doc-examples", "hoomd-rand" - ] + , "hoomd-order"] # Include only core hoomd simulation crates by default so that `cargo build`, # `cargo test` and similar commands operate only on the core code by default. @@ -74,6 +74,7 @@ assert2 = "=0.4.0" divan = "=0.1.21" paste = "=1.0.15" rand_sfc = "=0.2.0" +sphrs = "=0.2.2" rstest = "=0.26.1" tempfile = "=3.27.0" strum = "=0.28.0" @@ -97,9 +98,10 @@ parquet_derive = { version = "58.0.0", default-features=false } # Dependencies exposed via Public APIs. anyhow = "1.0.100" -approxim = "0.6.6" +approxim = {version = "0.6.6", features = ["num-complex"]} arrayvec = { version = "0.7.6", features = ["serde"] } num = { version = "0.4.3", features = ["serde"] } +num-complex = "0.4" rand = { version = "0.10.0", default-features=false, features = ["std", "std_rng"] } serde = { version = "1.0.228", features = ["derive"]} thiserror = "2.0.12" diff --git a/hoomd-order/Cargo.toml b/hoomd-order/Cargo.toml new file mode 100644 index 000000000..5a59fc39a --- /dev/null +++ b/hoomd-order/Cargo.toml @@ -0,0 +1,30 @@ +[package] +name = "hoomd-order" +description = "Evaluate Spherical Harmonics at a given Cartesian point." +version.workspace = true +edition.workspace = true +rust-version.workspace = true +homepage.workspace = true +documentation.workspace = true +repository.workspace = true +license.workspace = true +keywords.workspace = true +categories.workspace = true + +[dependencies] +hoomd-vector = "1.1.0" +num-complex.workspace = true + +[dev-dependencies] +approxim.workspace = true +divan.workspace = true +rand.workspace = true +rstest.workspace = true +sphrs.workspace = true + +[lints] +workspace = true + +[[bench]] +name = "harmonics" +harness = false diff --git a/hoomd-order/README.md b/hoomd-order/README.md new file mode 100644 index 000000000..1b0239649 --- /dev/null +++ b/hoomd-order/README.md @@ -0,0 +1,5 @@ +# hoomd-order + +Order parameter analysis library for [hoomd-rs]. + +[hoomd-rs]: https://github.com/glotzerlab/hoomd-rs diff --git a/hoomd-order/benches/harmonics.rs b/hoomd-order/benches/harmonics.rs new file mode 100644 index 000000000..f3cb427c8 --- /dev/null +++ b/hoomd-order/benches/harmonics.rs @@ -0,0 +1,34 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! ... +use hoomd_order::SphericalHarmonic; +use hoomd_vector::{Cartesian, InnerProduct}; +use rand::{Rng, RngExt, SeedableRng, rngs::StdRng}; + +/// Generate a random point on the unit sphere. +#[inline] +fn random_unit_point(rng: &mut R) -> Cartesian<3> { + let (x, y, z) = (rng.random(), rng.random(), rng.random()); + Cartesian::from([x, y, z]) +} +fn main() { + divan::main(); +} + +/// Measure per-point performance at each l. +#[divan::bench( + consts = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], + sample_count = 1_000, + sample_size = 10_000, +)] +fn recurrence(bencher: divan::Bencher<'_, '_>) { + let sh = SphericalHarmonic::::new(); + let mut rng = StdRng::seed_from_u64(1); + bencher + .with_inputs(|| { + let point = random_unit_point::(&mut rng); + (point.to_unit().expect("non-zero point")).0 + }) + .bench_local_values(|xyz| sh.eval(xyz)); +} diff --git a/hoomd-order/src/lib.rs b/hoomd-order/src/lib.rs new file mode 100644 index 000000000..3fc4cbf4f --- /dev/null +++ b/hoomd-order/src/lib.rs @@ -0,0 +1,540 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! Tools for evaluating complex spherical harmonics in rust. +//! +//! This library uses a recurrence relation to evaluate spherical harmonics of a +//! particular azimuthal quantum number `l` and all positive magnetic quantum numbers +//! `m=0..=L`. The approach taken is much faster than more general recurrences, which +//! typically attempt to evaluate all values of `l` up to the target. When computing +//! Steinhardt order parameters or similar algorithms, this code is much faster than +//! alternatives, with competetive numerical stability even out to large values of `l`. +//! +//! # Example +//! ``` +//! use hoomd_order::SphericalHarmonic; +//! use num_complex::Complex64; +//! use hoomd_vector::{Cartesian, InnerProduct}; +//! use approxim::assert_abs_diff_eq; +//! use std::f64::consts::PI; +//! +//! # fn main() -> Result<(), Box> { +//! // Initialize the SphericalHarmonic container, which can be reused +//! // to compute Y_6^m at a large number of points. +//! let y_6 = SphericalHarmonic::<6>::new(); +//! +//! // Values of m in 0..=L are returned as a HarmonicOutput container, which behaves +//! // like a [f64; L+1] array. +//! let (point, _) = Cartesian::<3>::from([1.0; 3]).to_unit()?; +//! let sh = y_6.eval(point); +//! assert_eq!(sh.len(), 6+1); +//! +//! // Zonal harmonic (m=0) is always purely real +//! assert_eq!(sh[0].im, 0.0); +//! +//! // Y_6^0 = sqrt(13/(4pi)) * P_6(1/sqrt(3)) = sqrt(13/(4pi)) * 2/9 +//! let expected_m0 = 2.0 * f64::sqrt(13.0 / (4.0 * PI)) / 9.0; +//! assert_abs_diff_eq!(sh[0].re, expected_m0, epsilon = 1e-15); +//! +//! /// Implement the Steinhardt order parameter q6. +//! fn q6(bonds: &[Cartesian<3>]) -> f64 { +//! let mut accum = [Complex64::ZERO; 7]; +//! let y6 = SphericalHarmonic::<6>::new(); +//! +//! for &bond in bonds { +//! let (unit_bond, _) = bond.to_unit().expect("Bond has zero distance!"); +//! let qlmi = y6.eval(unit_bond); +//! for m in 0..7 { accum[m] += qlmi[m]; } +//! } +//! +//! // We multiply the `m>0` components by two to account for `-m` contributions. +//! let sum_sq = accum[0].norm_sqr() +//! + 2.0 * accum[1..].iter().map(Complex64::norm_sqr).sum::(); +//! let n = bonds.len() as f64; +//! +//! (4.0 * PI / 13.0 * sum_sq).sqrt() / n +//! } +//! +//! // FCC nearest neighbors: permutations of (±1, ±1, 0) +//! let fcc_bonds: Vec> = [ +//! [-1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], +//! [-1.0, 0.0, -1.0], [-1.0, 0.0, 1.0], [1.0, 0.0, -1.0], [1.0, 0.0, 1.0], +//! [ 0.0, -1.0, -1.0], [ 0.0, -1.0, 1.0], [0.0, 1.0, -1.0], [0.0, 1.0, 1.0], +//! ].map(Cartesian::<3>::from).to_vec(); +//! assert_abs_diff_eq!(q6(&fcc_bonds), 0.57452416, epsilon = 1e-6); +//! # Ok(()) +//! # } +//! ``` + +use hoomd_vector::{Cartesian, Unit}; +use num_complex::Complex64; +use std::{ + f64::consts::{FRAC_1_SQRT_2, PI, SQRT_2}, + fmt, + ops::Index, +}; + +/// Precomputed coefficients for evaluating complex spherical harmonics of degree L. +/// +/// Once a [`SphericalHarmonic`] has been created with [`new`](Self::new), +/// [`eval`](Self::eval) can be called to rapidly evaluate the harmonic at a set of +/// points in three-dimensional space. +/// +/// ``` +/// use hoomd_order::SphericalHarmonic; +/// use hoomd_vector::{Cartesian, InnerProduct}; +/// +/// # fn main() -> Result<(), Box> { +/// let sh = SphericalHarmonic::<2>::new(); +/// let (point, _) = Cartesian::<3>::from([0.0, 0.0, 1.0]).to_unit()?; +/// let out = sh.eval(point); +/// +/// // m=0 (zonal harmonic) is always purely real +/// assert_eq!(out[0].im, 0.0); +/// assert_eq!(out.len(), 3); +/// # Ok(()) +/// # } +/// ``` +#[derive(Copy, Clone, Debug)] +pub struct SphericalHarmonic { + /// Initial value for the recurrence. + normalized_recurrence_seed: f64, + /// Coefficient of the `z * h[m]` term in the Legendre recurrence. + z_coeff: [f64; L], + /// Coefficient of the `rxy2 * h[m+1]` term in the Legendre recurrence. + rxy_coeff: [f64; L], +} + +impl SphericalHarmonic { + /// Precompute L-dependent coefficients. + #[must_use] + #[inline] + pub fn new() -> Self { + let normalized_recurrence_seed = { + let mut r = 1.0; + for k in 1..=L { + r *= (2 * k - 1) as f64 / (2 * k) as f64; + } + f64::sqrt((2 * L + 1) as f64 * r / (2.0 * PI)) * FRAC_1_SQRT_2 + }; + + let mut z_coeff = [0.0; L]; + let mut rxy_coeff = [0.0; L]; + + let sqrt_2l = f64::sqrt(2.0 * L as f64); + let mut carry = sqrt_2l; + + for m in (1..L).rev() { + let denom = f64::sqrt(((L - m) * (L + m + 1)) as f64); + z_coeff[m] = 2.0 * (m + 1) as f64 / denom; + rxy_coeff[m] = carry / denom; + carry = denom; + } + + // m=0 step: √2 fused into coefficients + if L > 0 { + let denom_0 = f64::sqrt((2 * L * (L + 1)) as f64); + z_coeff[0] = 2.0 * SQRT_2 / denom_0; + rxy_coeff[0] = carry * SQRT_2 / denom_0; + } + + Self { + normalized_recurrence_seed, + z_coeff, + rxy_coeff, + } + } + + /// Evaluate `Y_L^m` for m = 0..=L at spherical coordinates `(theta, phi)`. + /// + /// `theta` is the polar angle (from the z-axis), `phi` is the azimuthal angle. + /// + /// ``` + /// use approxim::assert_abs_diff_eq; + /// use hoomd_order::SphericalHarmonic; + /// use hoomd_vector::{Cartesian, InnerProduct, Unit}; + /// + /// # fn main() -> Result<(), Box> { + /// let sh = SphericalHarmonic::<2>::new(); + /// + /// let (x, y, z) = (0.6, 0.8, 0.0); + /// let (cart_point, _) = Cartesian::<3>::from([x, y, z]).to_unit()?; + /// let cartesian_result = sh.eval(cart_point); + /// + /// let (theta, phi) = (f64::acos(z), f64::atan2(y, x)); + /// let spherical_result = sh.eval_spherical(theta, phi); + /// + /// assert_abs_diff_eq!( + /// cartesian_result[0], + /// spherical_result[0], + /// epsilon = 1e-15 + /// ); + /// assert_abs_diff_eq!( + /// cartesian_result[1], + /// spherical_result[1], + /// epsilon = 1e-15 + /// ); + /// # Ok(()) + /// # } + /// ``` + #[must_use] + #[inline] + pub fn eval_spherical(&self, theta: f64, phi: f64) -> HarmonicOutput { + let (sin_theta, cos_theta) = theta.sin_cos(); + let (sin_phi, cos_phi) = phi.sin_cos(); + self.eval_unchecked([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]) + } + + /// Evaluate `Y_L^m` for m = 0..=L at a point on the unit sphere. + #[must_use] + #[inline] + pub fn eval(&self, point: Unit>) -> HarmonicOutput { + self.eval_unchecked(point.get().coordinates) + } + + /// Evaluate `Y_l^m` for a point *assumed* to be on the unit sphere. + #[must_use] + #[inline] + fn eval_unchecked(&self, point: [f64; 3]) -> HarmonicOutput { + let [x, y, z] = point; + let rxy2 = x * x + y * y; + + let h_0; + let mut h = [0.0; L]; + + if L == 0 { + h_0 = f64::sqrt(1.0 / (4.0 * PI)); + } else { + h[L - 1] = self.normalized_recurrence_seed; + let mut h_plus1 = 0.0; + + for m in (1..L).rev() { + h[m - 1] = self.z_coeff[m] * z * h[m] - rxy2 * self.rxy_coeff[m] * h_plus1; + h_plus1 = h[m]; + } + + h_0 = self.z_coeff[0] * z * h[0] - rxy2 * self.rxy_coeff[0] * h_plus1; + } + + let mut result = [Complex64::ZERO; L]; + + if L > 0 { + let mut cm = x; + let mut sm = y; + result[0] = Complex64::new(h[0] * cm, h[0] * sm); + + for m in 1..L { + let prev_cm = cm; + let prev_sm = sm; + cm = prev_cm * x - prev_sm * y; + sm = prev_cm * y + prev_sm * x; + result[m] = Complex64::new(h[m] * cm, h[m] * sm); + } + } + + HarmonicOutput { + m0: Complex64::new(h_0, 0.0), + mp: result, + } + } +} + +impl Default for SphericalHarmonic { + #[inline] + fn default() -> Self { + Self::new() + } +} + +/// Complex spherical harmonics `Y_L^m` for a single degree L. +/// +/// Index with `[m]` to access `Y_L^m` for m = 0..=L. +/// The m = 0 term is always purely real. +#[derive(Copy, Clone, Debug, PartialEq)] +pub struct HarmonicOutput { + /// `Y_L^0` (zonal harmonic, always real). + m0: Complex64, + /// `Y_L^m` for m = 1..=L, stored at index m − 1. + mp: [Complex64; L], +} + +impl Index for HarmonicOutput { + type Output = Complex64; + + #[inline] + fn index(&self, index: usize) -> &Complex64 { + match index { + 0 => &self.m0, + n => &self.mp[n - 1], + } + } +} + +impl IntoIterator for HarmonicOutput { + type Item = Complex64; + type IntoIter = + std::iter::Chain, std::array::IntoIter>; + + #[inline] + fn into_iter(self) -> Self::IntoIter { + std::iter::once(self.m0).chain(self.mp) + } +} + +impl fmt::Display for HarmonicOutput { + #[inline] + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + writeln!(f, "[")?; + for m in 0..=L { + writeln!(f, " {:+.12}{:+.12}i, // m={m}", self[m].re, self[m].im)?; + } + write!(f, "]") + } +} + +impl HarmonicOutput { + /// Get the harmonic degree L of the container. + #[inline] + #[must_use] + pub const fn l(&self) -> usize { + L + } + /// Iterate over all `L + 1` values, starting with `Y_L^0`. + /// + /// ``` + /// use approxim::assert_abs_diff_eq; + /// use hoomd_order::SphericalHarmonic; + /// use hoomd_vector::{Cartesian, InnerProduct}; + /// use num_complex::Complex64; + /// + /// # fn main() -> Result<(), Box> { + /// let (point, _) = Cartesian::<3>::from([0.0, 0.0, 1.0]).to_unit()?; + /// let out = SphericalHarmonic::<4>::new().eval(point); + /// + /// // Build Y_4^m for m = -4..=4: `Y_l^{-m} = (-1)^m · conj(Y_l^m)`. + /// let full: Vec = (1..=4) + /// .rev() + /// .map(|m| (-1.0f64).powi(m) * out[m as usize].conj()) + /// .chain(out.iter()) + /// .collect(); + /// assert_eq!(full.len(), 9); + /// assert_abs_diff_eq!(full[0].conj(), full[2 * out.l()]); + /// # Ok(()) + /// # } + /// ``` + #[inline] + pub fn iter(&self) -> impl Iterator + '_ { + std::iter::once(self.m0).chain(self.mp.iter().copied()) + } + /// The length of the container, equal to `L + 1`. + #[inline] + #[must_use] + pub const fn len(&self) -> usize { + L + 1 + } + /// Check if the container is empty. This will always be false. + #[inline] + #[must_use] + pub const fn is_empty(&self) -> bool { + false + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approxim::assert_abs_diff_eq; + use hoomd_vector::InnerProduct; + use rstest::rstest; + use std::marker::PhantomData; + + type Degree = PhantomData<[(); L]>; + fn degree() -> Degree { + Degree::default() + } + + fn unit(arr: [f64; 3]) -> Result>, hoomd_vector::Error> { + Ok(Cartesian::from(arr).to_unit()?.0) + } + + #[test] + fn l0() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<0>::new(); + let out = sh.eval(unit([0.0, 0.0, 1.0])?); + let expected = 1.0 / (2.0 * f64::sqrt(PI)); + assert_abs_diff_eq!(out[0], Complex64::new(expected, 0.0f64), epsilon = 1e-12); + assert_eq!(out.mp.len(), 0); + Ok(()) + } + + #[test] + fn l1_north_pole() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<1>::new(); + let out = sh.eval(unit([0.0, 0.0, 1.0])?); + let c = f64::sqrt(3.0 / (4.0 * PI)); + assert_abs_diff_eq!(out[0], Complex64::new(c, 0.0), epsilon = 1e-12); + assert_abs_diff_eq!(out[1], Complex64::ZERO, epsilon = 1e-12); + Ok(()) + } + + #[test] + fn l1_x_axis() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<1>::new(); + let out = sh.eval(unit([1.0, 0.0, 0.0])?); + let c = f64::sqrt(3.0 / (8.0 * PI)); + assert_abs_diff_eq!(out[0], Complex64::ZERO, epsilon = 1e-12); + assert_abs_diff_eq!(out[1], Complex64::new(c, 0.0), epsilon = 1e-12); + Ok(()) + } + + #[test] + fn l1_y_axis() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<1>::new(); + let out = sh.eval(unit([0.0, 1.0, 0.0])?); + let c = f64::sqrt(3.0 / (8.0 * PI)); + assert_abs_diff_eq!(out[0], Complex64::ZERO, epsilon = 1e-12); + assert_abs_diff_eq!(out[1], Complex64::new(0.0, c), epsilon = 1e-12); + Ok(()) + } + + #[test] + fn l2_finite() -> Result<(), hoomd_vector::Error> { + let inv3 = 3.0_f64.sqrt().recip(); + let sh = SphericalHarmonic::<2>::new(); + let out = sh.eval(unit([inv3, inv3, inv3])?); + assert_eq!(out.mp.len(), 2); + assert!(out.m0.re.is_finite()); + assert!(out.m0.im.is_finite()); + for v in &out.mp { + assert!(v.re.is_finite()); + assert!(v.im.is_finite()); + } + Ok(()) + } + + #[test] + fn into_iter_matches_index() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<4>::new(); + let reference = sh.eval(unit([0.6, 0.3, 0.4])?); + let out = sh.eval(unit([0.6, 0.3, 0.4])?); + for (m, val) in out.into_iter().enumerate() { + assert_abs_diff_eq!(val, reference[m], epsilon = 1e-15); + } + Ok(()) + } + + #[test] + fn iter_matches_index() -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::<4>::new(); + let out = sh.eval(unit([0.6, 0.3, 0.4])?); + let values: Vec<_> = out.iter().collect(); + assert_eq!(values.len(), 5); + for m in 0..=4 { + assert_abs_diff_eq!(values[m], out[m], epsilon = 1e-15); + } + Ok(()) + } + + /// Validate against sphrs via `Y_l^m` = (`S_l^{+m`} + i·S_l^{-m}) / √2. + fn check_against_sphrs(point: [f64; 3]) -> Result<(), hoomd_vector::Error> { + use sphrs::{Coordinates, RealSH, SHEval}; + let l = i64::try_from(L).expect("L would overflow i64"); + + let sh = SphericalHarmonic::::new(); + let out = sh.eval(unit(point)?); + let [x, y, z] = point; + let coords = Coordinates::cartesian(x, y, z); + + let expected_m0: f64 = RealSH::Spherical.eval(l, 0, &coords); + assert_abs_diff_eq!(out[0], Complex64::new(expected_m0, 0.0), epsilon = 1e-8); + + for m in 1..=L { + let m_i64 = i64::try_from(m).expect("m would overflow i64"); + let s_pos: f64 = RealSH::Spherical.eval(l, m_i64, &coords); + let s_neg: f64 = RealSH::Spherical.eval(l, -m_i64, &coords); + assert_abs_diff_eq!( + out[m], + Complex64::new(s_pos * FRAC_1_SQRT_2, s_neg * FRAC_1_SQRT_2), + epsilon = 1e-8 + ); + } + Ok(()) + } + + #[rstest] + #[expect( + clippy::used_underscore_binding, + reason = "Required for const generic parameterization." + )] + fn sphrs_test( + #[values( + degree::<0>(), + degree::<1>(), + degree::<2>(), + degree::<3>(), + degree::<4>(), + degree::<5>(), + degree::<6>(), + degree::<7>(), + degree::<8>(), + degree::<9>(), + degree::<10>() + // Values of L>10 overflow sphrs's factorial implementation + )] + _d: Degree, + #[values( + [0.0, 0.0, 1.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [3.0_f64.sqrt().recip(); 3], + [0.6_f64.sin() * 0.3_f64.cos(), 0.6_f64.sin() * 0.3_f64.sin(), 0.6_f64.cos()], + )] + point: [f64; 3], + ) { + check_against_sphrs::(point).unwrap(); + } + + /// Completeness: |`Y_l^0|²` + 2·Σ_{m=1}^l |`Y_l^m|²` = (2l+1) / (4π). + fn check_completeness(point: [f64; 3]) -> Result<(), hoomd_vector::Error> { + let sh = SphericalHarmonic::::new(); + let out = sh.eval(unit(point)?); + let mut sum = out[0].norm_sqr(); + for m in 1..=L { + sum += 2.0 * out[m].norm_sqr(); + } + let expected = (2 * L + 1) as f64 / (4.0 * PI); + assert_abs_diff_eq!(sum, expected, epsilon = 1e-10); + Ok(()) + } + + #[rstest] + #[expect( + clippy::used_underscore_binding, + reason = "Required for const generic parameterization." + )] + fn completeness_test( + #[values( + degree::<0>(), degree::<1>(), degree::<2>(), degree::<3>(), + degree::<4>(), degree::<5>(), degree::<6>(), degree::<7>(), + degree::<8>(), degree::<9>(), degree::<10>(), degree::<11>(), + degree::<12>(), degree::<13>(), degree::<14>(), degree::<15>(), + degree::<16>(), degree::<17>(), degree::<18>(), degree::<19>(), + degree::<20>(), degree::<21>(), degree::<22>(), degree::<23>(), + degree::<24>(), degree::<25>(), degree::<26>(), degree::<27>(), + degree::<28>(), degree::<29>(), degree::<30>(), degree::<31>(), + degree::<32>(), degree::<33>(), degree::<34>(), degree::<35>(), + degree::<36>(), degree::<37>(), degree::<38>(), degree::<39>(), + degree::<40>(), degree::<41>(), degree::<42>(), degree::<43>(), + degree::<44>(), degree::<45>(), degree::<46>(), degree::<47>(), + degree::<48>(), degree::<49>(), degree::<50>(), + )] + _d: Degree, + ) { + let point = [ + 0.7_f64.sin() * 0.3_f64.cos(), + 0.7_f64.sin() * 0.3_f64.sin(), + 0.7_f64.cos(), + ]; + check_completeness::(point).unwrap(); + } +}