From 04fc6c7c837e43b139ce2fca723ae500a9caadee Mon Sep 17 00:00:00 2001 From: mthran Date: Thu, 4 Jun 2026 11:44:03 -0400 Subject: [PATCH 01/16] add spherical_vec_cell --- Cargo.lock | 1 + hoomd-spatial/Cargo.toml | 1 + hoomd-spatial/src/lib.rs | 2 + hoomd-spatial/src/spherical_vec_cell.rs | 641 ++++++++++++++++++++++++ 4 files changed, 645 insertions(+) create mode 100644 hoomd-spatial/src/spherical_vec_cell.rs diff --git a/Cargo.lock b/Cargo.lock index af766146b..ad7a100ba 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3725,6 +3725,7 @@ name = "hoomd-spatial" version = "1.1.0" dependencies = [ "assert2", + "hoomd-manifold", "hoomd-utility", "hoomd-vector", "log", diff --git a/hoomd-spatial/Cargo.toml b/hoomd-spatial/Cargo.toml index d3d6bec1d..19e25042a 100644 --- a/hoomd-spatial/Cargo.toml +++ b/hoomd-spatial/Cargo.toml @@ -23,6 +23,7 @@ rustc-hash.workspace = true serde.workspace = true serde_with.workspace = true +hoomd-manifold.workspace = true hoomd-vector.workspace = true hoomd-utility.workspace = true diff --git a/hoomd-spatial/src/lib.rs b/hoomd-spatial/src/lib.rs index 1e5f48e31..5e7754e4b 100644 --- a/hoomd-spatial/src/lib.rs +++ b/hoomd-spatial/src/lib.rs @@ -41,10 +41,12 @@ use hoomd_utility::valid::PositiveReal; mod all_pairs; mod hash_cell; +mod spherical_vec_cell; mod vec_cell; pub use all_pairs::AllPairs; pub use hash_cell::{HashCell, HashCellBuilder}; +pub use spherical_vec_cell::{SphericalVecCell, SphericalVecCellBuilder}; pub use vec_cell::{VecCell, VecCellBuilder}; /// Allow incremental updates to points in the spatial data. diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs new file mode 100644 index 000000000..d1bfa8270 --- /dev/null +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -0,0 +1,641 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +#![expect( + clippy::cast_possible_truncation, + reason = "the necessary conversions are necessary and have been checked" +)] +#![expect( + clippy::cast_sign_loss, + reason = "the necessary conversions are necessary and have been checked" +)] + +//! Implement `VecCell` + +use serde::{Deserialize, Serialize}; +use serde_with::serde_as; +use std::{array, cmp::Eq, fmt, f64::consts::PI, hash::Hash, iter, marker::PhantomData, mem}; + +use log::trace; +use rustc_hash::FxHashMap; + +use hoomd_manifold::Spherical; +use hoomd_utility::valid::PositiveReal; + +use super::{PointUpdate, PointsNearBall, WithSearchRadius}; + +use crate::{IndexFromPosition, hash_cell::CellIndex}; + +/// TODO: docs +#[serde_as] +#[derive(Clone, Debug, Serialize, Deserialize)] +pub struct SphericalVecCell +where + K: Eq + Hash, +{ + /// The Euclidean width of each cell. + cell_width: PositiveReal, + + /// A map from cell indices to cell contents. + keys_map: Vec>, + + /// A map from particle indices to cell indices. + cell_index: FxHashMap>, + + /// The shape of `keys_map` is `(half_extent * 2 + 1).powi(D)`. + half_extent: u32, + + /// Pre-computed stencils. + #[serde_as(as = "Vec>")] + stencils: Vec>, +} + +/// Iterate over cell indices in row major order. +struct CellIndexIterator { + /// The current cell index. + cell_index: [i64; D], + /// The half extent of the cube. + half_extent: u32, +} + +impl CellIndexIterator { + /// Iterate over a cube. + /// + /// The cube extends from `[-half_extent, -half_extent, ..., -half_extent]` to + /// `[half_extent, half_extent, ..., half_extent]`. + fn cube(half_extent: u32) -> Self { + let mut cell_index = [-i64::from(half_extent); D]; + cell_index[D - 1] -= 1; + Self { + cell_index, + half_extent, + } + } + + /// Increment the cell index. + #[inline] + fn increment_cell_index(&mut self) -> Option<[i64; D]> { + self.cell_index[D - 1] += 1; + + for i in (0..D).rev() { + if self.cell_index[i] > self.half_extent.into() { + if i == 0 { + return None; + } + + self.cell_index[i] = -(i64::from(self.half_extent)); + self.cell_index[i - 1] += 1; + } + } + + Some(self.cell_index) + } +} + +impl Iterator for CellIndexIterator { + type Item = [i64; D]; + + fn next(&mut self) -> Option { + self.increment_cell_index() + } +} + +/// Generate the stencil for a given radius. +fn generate_stencil(radius: u32) -> Vec<[i64; D]> { + assert!(radius >= 1, "cell list must have a minimum radius of 1"); + + let mut result = CellIndexIterator::cube(radius).collect::>(); + result.sort_by(|a, b| { + let r_a = a.iter().map(|x| x.pow(2)); + let r_b = b.iter().map(|x| x.pow(2)); + r_a.cmp(r_b) + }); + result +} + +/// Generate the stencils up to a given radius. +pub(crate) fn generate_all_stencils(max_radius: u32) -> Vec> { + assert!(max_radius >= 1, "cell list must have a minimum radius of 1"); + let mut result = Vec::new(); + + // A cell list will never use radius=0, so stencil[i] stores the stencil needed for + // radius i+1. + for radius in 0..max_radius { + result.push(generate_stencil(radius + 1)); + } + result +} + +/// TODO: docs +pub struct SphericalVecCellBuilder { + /// Most commonly used search radius, in Euclidean metric. + nominal_search_radius: PositiveReal, + + /// Largest possible search radius, in Euclidean metric. + maximum_search_radius: f64, + + /// Track the key type. + phantom_key: PhantomData, +} + +impl SphericalVecCellBuilder +where + K: Copy + Eq + Hash, +{ + /// TODO + #[inline] + #[must_use] + pub fn nominal_search_radius(mut self, nominal_search_radius: PositiveReal) -> Self { + self.nominal_search_radius = nominal_search_radius; + self + } + + /// TODO + #[inline] + #[must_use] + pub fn maximum_search_radius(mut self, maximum_search_radius: f64) -> Self { + self.maximum_search_radius = maximum_search_radius; + self + } + + /// TODO + #[inline] + #[must_use] + pub fn build(self) -> SphericalVecCell { + let maximum_stencil_radius = + (self.maximum_search_radius / self.nominal_search_radius.get()).ceil() as u32; + let half_extent: u32 = 1; + + SphericalVecCell { + cell_width: self.nominal_search_radius, + keys_map: iter::repeat_n(Vec::new(), (half_extent * 2 + 1).pow(D as u32) as usize) + .collect(), + cell_index: FxHashMap::default(), + half_extent, + stencils: generate_all_stencils(maximum_stencil_radius.max(1)), + } + } +} + +impl Default for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// TODO + #[inline] + fn default() -> Self { + Self::builder().build() + } +} + +impl WithSearchRadius for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// TODO + #[inline] + fn with_search_radius(spherical_radius: PositiveReal) -> Self { + let euclidean_radius = (spherical_radius.get().clamp(0.0, PI/2.0)).sin(); + Self::builder() + .nominal_search_radius( + PositiveReal::try_from(euclidean_radius) + .expect("positive number given previous check") + ) + .maximum_search_radius(euclidean_radius) + .build() + } +} + +impl SphericalVecCell +where + K: Eq + Hash, +{ + /// Compute the cell index given a `Spherical` position in space. + #[inline] + fn cell_index_from_position(&self, position: &Spherical) -> [i64; D] { + std::array::from_fn(|j| (position.coordinates()[j] / self.cell_width.get()).floor() as i64) + } + + /// Compute the vector index from a cell index + /// + /// Returns `None` when the index is out of bounds. + #[inline] + fn map_index_from_cell(half_extent: u32, cell_index: &[i64; D]) -> Option { + assert!(D > 1); + + let mut vec_index: usize = 0; + let mut width = 1; + + for i in (0..D).rev() { + let needed_extent = cell_index[i].unsigned_abs(); + if needed_extent > u64::from(half_extent) { + return None; + } + let v: usize = (cell_index[i] + i64::from(half_extent)) + .try_into() + .expect("cell index should be in bounds"); + + vec_index += v * width; + width *= (half_extent * 2 + 1) as usize; + } + Some(vec_index) + } + + /// Get the keys in a given cell index + #[inline] + fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { + let index = Self::map_index_from_cell(self.half_extent, cell_index) + .expect("cell_index should be in bounds"); + &self.keys_map[index] + } +} + +impl SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Construct a `VecCell` builder. + /// + /// Use the builder to set any or all parameters and construct a [`VecCell`]. + /// + /// # Example + /// + /// ``` + /// use hoomd_spatial::VecCell; + /// + /// # fn main() -> Result<(), Box> { + /// let vec_cell = VecCell::::builder() + /// .nominal_search_radius(2.5.try_into()?) + /// .maximum_search_radius(7.5) + /// .build(); + /// # Ok(()) + /// # } + /// ``` + #[expect( + clippy::missing_panics_doc, + reason = "hard-coded constant will never panic" + )] + #[inline] + #[must_use] + pub fn builder() -> SphericalVecCellBuilder { + SphericalVecCellBuilder { + nominal_search_radius: 1.0 + .try_into() + .expect("hard-coded constant is a positive real"), + maximum_search_radius: 1.0, + phantom_key: PhantomData, + } + } + + /// Remove excess capacity from dynamically allocated arrays. + /// + /// At this time, `shrink_to_fit` only reduces the memory utilized by the + /// cell contents. It does not shrink the `D`-dimensional storage to match + /// the range spanned by points currently in the data structure. + #[inline] + pub fn shrink_to_fit(&mut self) { + for keys in &mut self.keys_map { + keys.shrink_to_fit(); + } + self.keys_map.shrink_to_fit(); + self.cell_index.shrink_to_fit(); + } + + /// Double the number of cells stored along each axis until it includes the target. + fn expand_to(&mut self, target: u32) { + if self.half_extent >= target { + return; + } + + let mut new_half_extent = self.half_extent.min(1) * 2; + + while new_half_extent < target { + new_half_extent *= 2; + } + + trace!("Expanding to {}^{} cells", new_half_extent * 2 + 1, D); + + let mut new_keys_map: Vec> = + iter::repeat_n(Vec::new(), (new_half_extent * 2 + 1).pow(D as u32) as usize).collect(); + let old_half_extent = self.half_extent; + let old_keys_map = &mut self.keys_map; + + for old_cell_index in CellIndexIterator::cube(old_half_extent) { + let old_vec_index = Self::map_index_from_cell(old_half_extent, &old_cell_index) + .expect("cell_index should be consistent with keys_map"); + let new_vec_index = Self::map_index_from_cell(new_half_extent, &old_cell_index) + .expect("old_cell_index should be inside the new keys_map"); + new_keys_map[new_vec_index] = mem::take(&mut old_keys_map[old_vec_index]); + } + + self.half_extent = new_half_extent; + self.keys_map = new_keys_map; + } +} + +impl PointUpdate, K> for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Insert or update a point identified by a key. + /// + /// TODO + #[inline] + fn insert(&mut self, key: K, position: Spherical) { + let cell_index = self.cell_index_from_position(&position); + let old_cell_index = self.cell_index.insert(key, CellIndex(cell_index)); + let map_index = + Self::map_index_from_cell(self.half_extent, &cell_index).unwrap_or_else(|| { + let max_half_extent = cell_index + .iter() + .map(|x| x.unsigned_abs()) + .reduce(u64::max) + .expect("D should be greater than 1"); + self.expand_to( + max_half_extent + .try_into() + .expect("max extent cannot exceed u32::MAX"), + ); + Self::map_index_from_cell(self.half_extent, &cell_index) + .expect("cell_index should be in the expanded VecCell") + }); + + // This checks if old_cell_index is None or if it is different from the new cell index. + if old_cell_index != Some(CellIndex(cell_index)) { + // Add the particle index to the new cell index vector. + self.keys_map[map_index].push(key); + + if let Some(old_cell_index) = old_cell_index { + // If the particle was in a different cell, we need to remove it from the old cell. + let old_map_index = Self::map_index_from_cell(self.half_extent, &old_cell_index.0) + .expect("cell_index and keys_map should agree"); + let old_keys = &mut self.keys_map[old_map_index]; + if let Some(pos) = old_keys.iter().position(|x| *x == key) { + old_keys.swap_remove(pos); + } + } + } + } + + /// Remove the point with the given key. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, VecCell}; + /// + /// let mut vec_cell = VecCell::default(); + /// vec_cell.insert(0, [1.25, 2.5].into()); + /// + /// vec_cell.remove(&0) + /// ``` + #[inline] + fn remove(&mut self, key: &K) { + let cell_index = self.cell_index.remove(key); + if let Some(cell_index) = cell_index { + let map_index = Self::map_index_from_cell(self.half_extent, &cell_index.0); + if let Some(map_index) = map_index { + let keys = &mut self.keys_map[map_index]; + if let Some(idx) = keys.iter().position(|x| x == key) { + keys.swap_remove(idx); + } + } + } + } + + /// Get the number of points in the spatial data structure. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, VecCell}; + /// + /// let mut vec_cell = VecCell::default(); + /// vec_cell.insert(0, [1.25, 2.5].into()); + /// + /// assert_eq!(vec_cell.len(), 1) + /// ``` + #[inline] + fn len(&self) -> usize { + self.cell_index.len() + } + + /// Test if the spatial data structure is empty. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, VecCell}; + /// + /// let mut vec_cell = VecCell::::default(); + /// + /// assert!(vec_cell.is_empty()); + /// ``` + #[inline] + fn is_empty(&self) -> bool { + self.cell_index.is_empty() + } + + /// Test if the spatial data structure contains a key. + /// ``` + /// use hoomd_spatial::{PointUpdate, VecCell}; + /// + /// let mut vec_cell = VecCell::default(); + /// vec_cell.insert(0, [1.25, 2.5].into()); + /// + /// assert!(vec_cell.contains_key(&0)); + /// ``` + #[inline] + fn contains_key(&self, key: &K) -> bool { + self.cell_index.contains_key(key) + } + + /// Remove all points. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, VecCell}; + /// + /// let mut vec_cell = VecCell::default(); + /// vec_cell.insert(0, [1.25, 2.5].into()); + /// + /// vec_cell.clear(); + /// ``` + #[inline] + fn clear(&mut self) { + self.cell_index.clear(); + for keys in &mut self.keys_map { + keys.clear(); + } + } +} + +/// Iterate over keys in the cell list around a given center cell. +struct PointsIterator<'a, K, const D: usize> +where + K: Eq + Hash, +{ + /// Keys of the current cell iteration (None if the cell is empty) + keys: Option<&'a Vec>, + + /// The cell list we are iterating in. + cell_list: &'a SphericalVecCell, + + /// Current location of the iteration in the cell. + index_in_current_cell: usize, + + /// Current location of the iteration in the stencil. + current_stencil: usize, + + /// Cell offsets to iterate over. + stencil: &'a [[i64; D]], + + /// The cell at the center of the iteration. + center: [i64; D], +} + +impl Iterator for PointsIterator<'_, K, D> +where + K: Copy + Eq + Hash, +{ + type Item = K; + + #[inline] + fn next(&mut self) -> Option { + loop { + if let Some(keys) = self.keys + && self.index_in_current_cell < keys.len() + { + let last_index = self.index_in_current_cell; + self.index_in_current_cell += 1; + return Some(keys[last_index]); + } + + self.index_in_current_cell = 0; + self.current_stencil += 1; + + if self.current_stencil >= self.stencil.len() { + return None; + } + + let cell_index = + array::from_fn(|i| self.center[i] + self.stencil[self.current_stencil][i]); + let map_index = + SphericalVecCell::::map_index_from_cell(self.cell_list.half_extent, &cell_index); + self.keys = map_index.map(|index| &self.cell_list.keys_map[index]); + } + } +} + +impl PointsNearBall, K> for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Find all the points that *might* be in the given ball. + /// + /// `points_near_ball` will iterate over all points in the given ball *and + /// possibly others as well*. [`VecCell`] may iterate over the points in + /// any order. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, PointsNearBall, VecCell}; + /// + /// let mut vec_cell = VecCell::default(); + /// vec_cell.insert(0, [1.25, 0.0].into()); + /// vec_cell.insert(1, [3.25, 0.75].into()); + /// vec_cell.insert(2, [-10.0, 12.0].into()); + /// + /// for key in vec_cell.points_near_ball(&[2.0, 0.0].into(), 1.0) { + /// println!("{key}"); + /// } + /// ``` + /// Prints (in any order): + /// ```text + /// 0 + /// 1 + /// ``` + /// + /// # Panics + /// + /// Panics when `radius` is larger than the *maximum search radius* + /// provided at construction, rounded up to the nearest integer multiple + /// of the *nominal search radius*. + #[inline] + fn points_near_ball(&self, position: &Spherical, radius: f64) -> impl Iterator { + let stencil_index = (radius / self.cell_width.get()).ceil() as usize - 1; + assert!( + stencil_index < self.stencils.len(), + "search radius must be less than or equal to the maximum search radius" + ); + + let center = self.cell_index_from_position(position); + let stencil = &self.stencils[stencil_index]; + let map_index = Self::map_index_from_cell( + self.half_extent, + &array::from_fn(|i| center[i] + stencil[0][i]), + ); + + PointsIterator { + keys: map_index.map(|index| &self.keys_map[index]), + cell_list: self, + index_in_current_cell: 0, + current_stencil: 0, + stencil, + center, + } + } +} + +impl fmt::Display for SphericalVecCell +where + K: Eq + Hash, +{ + /// Summarize the contents of the cell list. + /// + /// This is a slow operation. It is meant to be printed to logs only + /// occasionally, such as at the end of a benchmark or simulation. + /// + /// # Example + /// + /// ``` + /// use hoomd_spatial::VecCell; + /// use log::info; + /// + /// let vec_cell = VecCell::::default(); + /// + /// info!("{vec_cell}"); + /// ``` + #[allow( + clippy::missing_inline_in_public_items, + reason = "no need to inline display" + )] + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let largest_cell_size = self.keys_map.iter().map(Vec::len).fold(0, usize::max); + + writeln!(f, "VecCell:")?; + writeln!( + f, + "- {} total cells with {} cells on each side.", + self.keys_map.len(), + self.half_extent * 2 + 1 + )?; + writeln!(f, "- {} points.", self.cell_index.len())?; + writeln!( + f, + "- Nominal, maximum search radii: {}, {}", + self.cell_width, + self.cell_width.get() * self.stencils.len() as f64 + )?; + write!(f, "- Largest cell length: {largest_cell_size}") + } +} + +impl IndexFromPosition> for SphericalVecCell +where + K: Eq + Hash, +{ + type Location = [i64; D]; + + #[inline] + fn location_from_position(&self, position: &Spherical) -> Self::Location { + self.cell_index_from_position(position) + } +} \ No newline at end of file From a1f11be1e38fbc310b419c8eb933c03b58ece5f9 Mon Sep 17 00:00:00 2001 From: mthran Date: Thu, 4 Jun 2026 16:52:12 -0400 Subject: [PATCH 02/16] implement checkerboard for spherical --- hoomd-mc/src/hypercuboid.rs | 55 +++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index 88ecef8d0..f62fc7953 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -19,6 +19,7 @@ use serde_with::serde_as; use std::array; use hoomd_geometry::shape::Hypercuboid; +use hoomd_manifold::Spherical; use hoomd_microstate::boundary::{Closed, Periodic}; use hoomd_utility::valid::PositiveReal; use hoomd_vector::Cartesian; @@ -136,6 +137,60 @@ impl Checkerboard> for HypercuboidCheckerboard { } } +impl Checkerboard> for HypercuboidCheckerboard { + #[inline] + fn point_to_space_index(&self, point: &Spherical) -> Option { + let p = *point.point() - self.origin; + let mut space_multi_index: [i64; N] = + array::from_fn(|i| (p.coordinates[i] / self.space_width[i]).floor() as i64); + + for (index, shape, periodic) in izip!(&mut space_multi_index, self.shape, self.periodic) { + // The origin is in the lower left corner of the box and may be up + // to one space width to the left of simulation boundary. Therefore, + // negative indices are out of bounds (and should never been seen + // for wrapped inputs). + if *index < 0 { + return None; + } + + if periodic { + // In periodic boundaries, the checkerboard spaces end before + // the right side. The space at the rightmost edge is identical + // to space 0 to make the checkerboard coloring commensurate + // with the periodic boundary conditions. + if *index as usize == shape { + *index = 0; + } + if *index as usize > shape { + return None; + } + } else if *index as usize >= shape { + // In non-periodic boundaries, the checkerboard extends one full + // space to the right of the simulation boundary so that when + // it is shifted to the left it will still cover the boundary. + // Therefore, any points outside the checkerboard shape are + // invalid. + return None; + } + } + + Some(Self::multi_index_to_index( + array::from_fn(|i| space_multi_index[i] as usize), + self.shape, + )) + } + + #[inline] + fn space_indices_by_color(&self) -> &[Vec] { + &self.space_indices_by_color + } + + #[inline] + fn num_spaces(&self) -> usize { + self.shape.iter().product() + } +} + impl HypercuboidCheckerboard { /// Collapse a multi-dimensional index to a single value in `[0, num_spaces]`. #[inline] From 2aad1bc7acbf053a3d7658535c4d261f73158b49 Mon Sep 17 00:00:00 2001 From: mthran Date: Thu, 4 Jun 2026 17:00:15 -0400 Subject: [PATCH 03/16] construct_space_indices_by_color for HypercuboidCheckerboard<4> --- hoomd-mc/src/hypercuboid.rs | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index f62fc7953..23600b84b 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -310,6 +310,38 @@ impl HypercuboidCheckerboard { return result; } + if N == 4 { + for offset_l in 0..=1 { + for offset_k in 0..=1 { + for offset_j in 0..=1 { + for offset_i in 0..=1 { + let mut space_indices = Vec::new(); + let mut multi_index = [0; N]; + + for l in 0..shape[0] / 2 { + multi_index[0] = 2 * l + offset_l; + for k in 0..shape[1] / 2 { + multi_index[1] = 2 * k + offset_k; + for j in 0..shape[2] / 2 { + multi_index[2] = 2 * j + offset_j; + for i in 0..shape[3] / 2 { + multi_index[3] = 2 * i + offset_i; + space_indices + .push(Self::multi_index_to_index(multi_index, shape)); + } + } + } + } + + result.push(space_indices); + } + } + } + } + + return result; + } + todo!("Implement a general method"); } From 5d9790746b5d37090097e1b2ad5f823c77474619 Mon Sep 17 00:00:00 2001 From: mthran Date: Fri, 5 Jun 2026 12:31:26 -0400 Subject: [PATCH 04/16] switch norm check to norm squared check in from_cartesian_coordinates for Spherical --- hoomd-manifold/src/sphere.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 98c0a9100..41f40b70d 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -53,7 +53,7 @@ impl Spherical { #[inline] #[must_use] pub fn from_cartesian_coordinates(point: Cartesian) -> Spherical { - let rad = point.norm(); + let rad = point.norm_squared(); assert_relative_eq!(rad, 1.0_f64, epsilon = 1e-6); Spherical { point } } From 915a96f117e208fa369a7fd4abb32bc96c516562 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Sun, 7 Jun 2026 17:32:03 -0400 Subject: [PATCH 05/16] test code for SphericalVecCell --- Cargo.lock | 1 + hoomd-mc/src/hypercuboid.rs | 6 +- hoomd-spatial/Cargo.toml | 1 + hoomd-spatial/src/spherical_vec_cell.rs | 553 +++++++++++++++++++++--- 4 files changed, 501 insertions(+), 60 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index ad7a100ba..e599cc796 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3724,6 +3724,7 @@ dependencies = [ name = "hoomd-spatial" version = "1.1.0" dependencies = [ + "approxim", "assert2", "hoomd-manifold", "hoomd-utility", diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index 23600b84b..00d4f50e1 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -326,8 +326,10 @@ impl HypercuboidCheckerboard { multi_index[2] = 2 * j + offset_j; for i in 0..shape[3] / 2 { multi_index[3] = 2 * i + offset_i; - space_indices - .push(Self::multi_index_to_index(multi_index, shape)); + space_indices.push(Self::multi_index_to_index( + multi_index, + shape, + )); } } } diff --git a/hoomd-spatial/Cargo.toml b/hoomd-spatial/Cargo.toml index 19e25042a..2a3c61eb7 100644 --- a/hoomd-spatial/Cargo.toml +++ b/hoomd-spatial/Cargo.toml @@ -28,6 +28,7 @@ hoomd-vector.workspace = true hoomd-utility.workspace = true [dev-dependencies] +approxim.workspace = true assert2.workspace = true rand.workspace = true rstest.workspace = true diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index d1bfa8270..323833599 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -10,11 +10,11 @@ reason = "the necessary conversions are necessary and have been checked" )] -//! Implement `VecCell` +//! Implement `SphericalVecCell` use serde::{Deserialize, Serialize}; use serde_with::serde_as; -use std::{array, cmp::Eq, fmt, f64::consts::PI, hash::Hash, iter, marker::PhantomData, mem}; +use std::{array, cmp::Eq, f64::consts::PI, fmt, hash::Hash, iter, marker::PhantomData, mem}; use log::trace; use rustc_hash::FxHashMap; @@ -126,9 +126,25 @@ pub(crate) fn generate_all_stencils(max_radius: u32) -> Vec Result<(), Box> { +/// let two_sphere_vec_cell = SphericalVecCell::::builder() +/// .spherical_nominal_search_radius((PI / 12.0).try_into()?) +/// .spherical_maximum_search_radius(PI / 4.0) +/// .build(); +/// # Ok(()) +/// # } +/// ``` pub struct SphericalVecCellBuilder { - /// Most commonly used search radius, in Euclidean metric. + /// Most commonly used search radius, in Euclidean metruic. nominal_search_radius: PositiveReal, /// Largest possible search radius, in Euclidean metric. @@ -142,23 +158,65 @@ impl SphericalVecCellBuilder where K: Copy + Eq + Hash, { - /// TODO + /// Choose the search radius from a given `Spherical` search radius. + /// + /// The nominal search radius is the edge length of each hypercube in + /// `SphericalVecCell`. Note that `Spherical` geodesic distances are always larger + /// than the euclidean "line-of-sight" distance. + /// + /// # Panics + /// Method will panic if `spherical_nominal_search_radius` is larger than $`\pi/2`$. #[inline] #[must_use] - pub fn nominal_search_radius(mut self, nominal_search_radius: PositiveReal) -> Self { - self.nominal_search_radius = nominal_search_radius; + pub fn spherical_nominal_search_radius( + mut self, + spherical_nominal_search_radius: PositiveReal, + ) -> Self { + self.nominal_search_radius = (spherical_nominal_search_radius.get().clamp(0.0, PI / 2.0)) + .sin() + .try_into() + .expect("clamp ensures number is positive"); self } - /// TODO + /// Choose the largest search radius from a given `Spherical` search radius. + /// + /// As in `VecCell`, the maximum radius is rounded up to the nearest + /// integer multiple of the nominal search radius. [`SphericalVecCell`] + /// will panic when asked to search for points within a radius larger + /// than the maximum. #[inline] #[must_use] - pub fn maximum_search_radius(mut self, maximum_search_radius: f64) -> Self { - self.maximum_search_radius = maximum_search_radius; + pub fn spherical_maximum_search_radius(mut self, spherical_maximum_search_radius: f64) -> Self { + self.maximum_search_radius = (spherical_maximum_search_radius.clamp(0.0, PI / 2.0)).sin(); self } - /// TODO + /// Choose the search radius from a given Euclidean search radius. + #[inline] + #[must_use] + pub fn euclidean_nominal_search_radius( + mut self, + euclidean_nominal_search_radius: PositiveReal, + ) -> Self { + self.nominal_search_radius = euclidean_nominal_search_radius; + self + } + + /// Choose the largest search radius from a given Euclidean search radius. + /// + /// As in `VecCell`, the maximum radius is rounded up to the nearest + /// integer multiple of the nominal search radius. [`SphericalVecCell`] + /// will panic when asked to search for points within a radius larger + /// than the maximum. + #[inline] + #[must_use] + pub fn euclidean_maximum_search_radius(mut self, euclidean_maximum_search_radius: f64) -> Self { + self.maximum_search_radius = euclidean_maximum_search_radius; + self + } + + /// Construct the [`SphericalVecCell`] with the chosen parameters. #[inline] #[must_use] pub fn build(self) -> SphericalVecCell { @@ -192,16 +250,34 @@ impl WithSearchRadius for SphericalVecCell where K: Copy + Eq + Hash, { - /// TODO + /// Construct a [`SphericalVecCell`] with a given search radius in the + /// `Spherical` metric. + /// + /// Both the nominal and maximum search radii are set to `spherical_radius` + /// + /// # Example + /// + /// ``` + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{SphericalVecCell, WithSearchRadius}; + /// use std::f64::consts::PI; + /// + /// # fn main() -> Result<(), Box> { + /// let two_sphere_vec_cell = SphericalVecCell::::with_search_radius( + /// (PI / 2.0).try_into()?, + /// ); + /// # Ok(()) + /// # } + /// ``` #[inline] - fn with_search_radius(spherical_radius: PositiveReal) -> Self { - let euclidean_radius = (spherical_radius.get().clamp(0.0, PI/2.0)).sin(); + fn with_search_radius(radius: PositiveReal) -> Self { + let euclidean_radius = (radius.get().clamp(0.0, PI / 2.0)).sin(); Self::builder() - .nominal_search_radius( + .euclidean_nominal_search_radius( PositiveReal::try_from(euclidean_radius) - .expect("positive number given previous check") + .expect("positive number given previous check"), ) - .maximum_search_radius(euclidean_radius) + .euclidean_maximum_search_radius(euclidean_radius) .build() } } @@ -242,6 +318,7 @@ where } /// Get the keys in a given cell index + #[cfg(test)] #[inline] fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { let index = Self::map_index_from_cell(self.half_extent, cell_index) @@ -254,19 +331,20 @@ impl SphericalVecCell where K: Copy + Eq + Hash, { - /// Construct a `VecCell` builder. + /// Construct a `SphericalVecCell` builder. /// - /// Use the builder to set any or all parameters and construct a [`VecCell`]. + /// Use the builder to set any or all parameters and construct a [`SphericalVecCell`]. /// /// # Example /// /// ``` - /// use hoomd_spatial::VecCell; + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::SphericalVecCell; /// /// # fn main() -> Result<(), Box> { - /// let vec_cell = VecCell::::builder() - /// .nominal_search_radius(2.5.try_into()?) - /// .maximum_search_radius(7.5) + /// let two_sphere_vec_cell = SphericalVecCell::::builder() + /// .euclidean_nominal_search_radius(0.2.try_into()?) + /// .euclidean_maximum_search_radius(0.4) /// .build(); /// # Ok(()) /// # } @@ -339,7 +417,23 @@ where { /// Insert or update a point identified by a key. /// - /// TODO + /// # Example + /// ``` + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// use std::f64::consts::PI; + /// + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// + /// spherical_vec_cell.insert( + /// 0, + /// Spherical::<4>::from_polar_coordinates( + /// PI / 4.0, + /// PI / 2.0, + /// 3.0 * PI / 2.0, + /// ), + /// ); + /// ``` #[inline] fn insert(&mut self, key: K, position: Spherical) { let cell_index = self.cell_index_from_position(&position); @@ -381,12 +475,17 @@ where /// /// # Example /// ``` - /// use hoomd_spatial::{PointUpdate, VecCell}; + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// use std::f64::consts::PI; /// - /// let mut vec_cell = VecCell::default(); - /// vec_cell.insert(0, [1.25, 2.5].into()); + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// spherical_vec_cell.insert( + /// 0, + /// Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 4.0, PI), + /// ); /// - /// vec_cell.remove(&0) + /// spherical_vec_cell.remove(&0) /// ``` #[inline] fn remove(&mut self, key: &K) { @@ -406,12 +505,17 @@ where /// /// # Example /// ``` - /// use hoomd_spatial::{PointUpdate, VecCell}; + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// use std::f64::consts::PI; /// - /// let mut vec_cell = VecCell::default(); - /// vec_cell.insert(0, [1.25, 2.5].into()); + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// spherical_vec_cell + /// .insert(0, Spherical::<3>::from_polar_coordinates(0.0, 0.0)); + /// spherical_vec_cell + /// .insert(1, Spherical::<3>::from_polar_coordinates(PI / 4.0, 0.0)); /// - /// assert_eq!(vec_cell.len(), 1) + /// assert_eq!(spherical_vec_cell.len(), 2) /// ``` #[inline] fn len(&self) -> usize { @@ -422,11 +526,11 @@ where /// /// # Example /// ``` - /// use hoomd_spatial::{PointUpdate, VecCell}; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; /// - /// let mut vec_cell = VecCell::::default(); + /// let mut spherical_vec_cell = SphericalVecCell::::default(); /// - /// assert!(vec_cell.is_empty()); + /// assert!(spherical_vec_cell.is_empty()); /// ``` #[inline] fn is_empty(&self) -> bool { @@ -435,12 +539,17 @@ where /// Test if the spatial data structure contains a key. /// ``` - /// use hoomd_spatial::{PointUpdate, VecCell}; + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// use std::f64::consts::PI; /// - /// let mut vec_cell = VecCell::default(); - /// vec_cell.insert(0, [1.25, 2.5].into()); + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// spherical_vec_cell.insert( + /// 0, + /// Spherical::<3>::from_polar_coordinates(3.0 * PI / 5.0, 0.0), + /// ); /// - /// assert!(vec_cell.contains_key(&0)); + /// assert!(spherical_vec_cell.contains_key(&0)); /// ``` #[inline] fn contains_key(&self, key: &K) -> bool { @@ -451,12 +560,17 @@ where /// /// # Example /// ``` - /// use hoomd_spatial::{PointUpdate, VecCell}; + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// use std::f64::consts::PI; /// - /// let mut vec_cell = VecCell::default(); - /// vec_cell.insert(0, [1.25, 2.5].into()); + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// spherical_vec_cell.insert( + /// 0, + /// Spherical::<4>::from_polar_coordinates(PI / 4.0, PI / 2.0, 0.0), + /// ); /// - /// vec_cell.clear(); + /// spherical_vec_cell.clear(); /// ``` #[inline] fn clear(&mut self) { @@ -517,8 +631,10 @@ where let cell_index = array::from_fn(|i| self.center[i] + self.stencil[self.current_stencil][i]); - let map_index = - SphericalVecCell::::map_index_from_cell(self.cell_list.half_extent, &cell_index); + let map_index = SphericalVecCell::::map_index_from_cell( + self.cell_list.half_extent, + &cell_index, + ); self.keys = map_index.map(|index| &self.cell_list.keys_map[index]); } } @@ -528,22 +644,33 @@ impl PointsNearBall, K> for SphericalVecCell::default(); + /// spherical_vec_cell + /// .insert(0, Spherical::<3>::from_polar_coordinates(PI / 12.0, 0.0)); + /// spherical_vec_cell + /// .insert(1, Spherical::<3>::from_polar_coordinates(PI / 12.0, PI)); + /// spherical_vec_cell.insert( + /// 2, + /// Spherical::<3>::from_polar_coordinates(2.0 * PI / 3.0, 0.0), + /// ); /// - /// for key in vec_cell.points_near_ball(&[2.0, 0.0].into(), 1.0) { + /// for key in spherical_vec_cell.points_near_ball( + /// &Spherical::<3>::from_cartesian_coordinates([0.0, 0.0, 1.0].into()), + /// PI / 4.0, + /// ) { /// println!("{key}"); /// } /// ``` @@ -560,7 +687,9 @@ where /// of the *nominal search radius*. #[inline] fn points_near_ball(&self, position: &Spherical, radius: f64) -> impl Iterator { - let stencil_index = (radius / self.cell_width.get()).ceil() as usize - 1; + // convert spherical distance to Euclidean distance + let euclidean_radius = radius.asin(); + let stencil_index = (euclidean_radius / self.cell_width.get()).ceil() as usize - 1; assert!( stencil_index < self.stencils.len(), "search radius must be less than or equal to the maximum search radius" @@ -596,12 +725,12 @@ where /// # Example /// /// ``` - /// use hoomd_spatial::VecCell; + /// use hoomd_spatial::SphericalVecCell; /// use log::info; /// - /// let vec_cell = VecCell::::default(); + /// let spherical_vec_cell = SphericalVecCell::::default(); /// - /// info!("{vec_cell}"); + /// info!("{spherical_vec_cell}"); /// ``` #[allow( clippy::missing_inline_in_public_items, @@ -638,4 +767,312 @@ where fn location_from_position(&self, position: &Spherical) -> Self::Location { self.cell_index_from_position(position) } -} \ No newline at end of file +} + +#[cfg(test)] +mod tests { + use super::*; + use hoomd_manifold::{Spherical, SphericalDisk}; + use rand::{ + RngExt, SeedableRng, + distr::{Distribution, Uniform}, + rngs::StdRng, + }; + use rstest::*; + + use approxim::assert_relative_eq; + use assert2::{assert, check}; + use std::f64::consts::PI; + + #[test] + fn two_sphere_cell_index() { + let spherical_cell_list = SphericalVecCell::::builder() + .euclidean_nominal_search_radius(0.5.try_into().expect("hard-coded positive number")) + .build(); + check!( + spherical_cell_list + .cell_index_from_position(&Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0)) + == [2, 0, 0] + ); + check!( + spherical_cell_list.cell_index_from_position(&Spherical::<3>::from_polar_coordinates( + PI / 2.0, + PI / 2.0 + )) == [0, 2, 0] + ); + check!( + spherical_cell_list + .cell_index_from_position(&Spherical::<3>::from_polar_coordinates(0.0, 0.0)) + == [0, 0, 2] + ); + } + + #[test] + fn three_sphere_cell_index() { + let spherical_cell_list = SphericalVecCell::::builder() + .euclidean_nominal_search_radius(0.5.try_into().expect("hard-coded positive number")) + .build(); + check!( + spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( + PI / 2.0, + 0.0, + 0.0 + )) == [2, 0, 0, 0] + ); + check!( + spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( + PI / 2.0, + PI / 2.0, + 0.0 + )) == [0, 2, 0, 0] + ); + check!( + spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( + PI / 2.0, + PI / 2.0, + PI / 2.0 + )) == [0, 0, 2, 0] + ); + check!( + spherical_cell_list + .cell_index_from_position(&Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0)) + == [0, 0, 0, 2] + ); + } + + #[test] + fn two_sphere_geodesic_search_radius() { + let spherical_cell_list = SphericalVecCell::::with_search_radius( + (PI / 2.0).try_into().expect("hard-coded positive number"), + ); + assert!(spherical_cell_list.cell_width.get() == 1.0); + + let spherical_cell_list = SphericalVecCell::::with_search_radius( + (PI / 4.0).try_into().expect("hard-coded positive number"), + ); + assert_relative_eq!( + spherical_cell_list.cell_width.get(), + (0.5_f64).sqrt(), + epsilon = 1e-12 + ); + } + + #[test] + fn three_sphere_geodesic_search_radius() { + let spherical_cell_list = SphericalVecCell::::with_search_radius( + (PI / 2.0).try_into().expect("hard-coded positive number"), + ); + assert!(spherical_cell_list.cell_width.get() == 1.0); + + let spherical_cell_list = SphericalVecCell::::with_search_radius( + (PI / 4.0).try_into().expect("hard-coded positive number"), + ); + assert_relative_eq!( + spherical_cell_list.cell_width.get(), + (0.5_f64).sqrt(), + epsilon = 1e-12 + ); + } + + #[test] + fn two_sphere_insert() { + let mut sphere_cell_list = SphericalVecCell::::default(); + + sphere_cell_list.insert(0, Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0)); + sphere_cell_list.insert( + 1, + Spherical::<3>::from_polar_coordinates(PI / 2.0 - 0.1, PI), + ); + sphere_cell_list.insert(2, Spherical::<3>::from_polar_coordinates(PI / 4.0, 0.0)); + sphere_cell_list.insert( + 3, + Spherical::<3>::from_polar_coordinates(PI / 2.0, 3.0 * PI / 2.0 + 0.001), + ); + sphere_cell_list.insert( + 4, + Spherical::<3>::from_polar_coordinates(PI / 4.0, 3.0 * PI / 2.0 + 0.001), + ); + + assert!(sphere_cell_list.cell_index.get(&0) == Some(&CellIndex([1, 0, 0]))); + assert!(sphere_cell_list.cell_index.get(&1) == Some(&CellIndex([-1, 0, 0]))); + assert!(sphere_cell_list.cell_index.get(&2) == Some(&CellIndex([0, 0, 0]))); + assert!(sphere_cell_list.cell_index.get(&3) == Some(&CellIndex([0, -1, 0]))); + assert!(sphere_cell_list.cell_index.get(&4) == Some(&CellIndex([0, -1, 0]))); + + let keys = sphere_cell_list.get_keys(&[0, 0, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&2)); + + let keys = sphere_cell_list.get_keys(&[1, 0, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&0)); + + let keys = sphere_cell_list.get_keys(&[-1, 0, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&1)); + + let keys = sphere_cell_list.get_keys(&[0, -1, 0]); + assert!(keys.len() == 2); + check!(keys.contains(&3)); + check!(keys.contains(&4)); + } + + #[test] + fn three_sphere_insert() { + let mut sphere_cell_list = SphericalVecCell::::default(); + + sphere_cell_list.insert( + 0, + Spherical::<4>::from_polar_coordinates(PI / 2.0 - 0.1, 0.0, 0.0), + ); + sphere_cell_list.insert( + 1, + Spherical::<4>::from_polar_coordinates( + PI / 2.0 - 0.1, + PI / 2.0, + 3.0 * PI / 2.0 + 0.001, + ), + ); + sphere_cell_list.insert( + 2, + Spherical::<4>::from_polar_coordinates(PI / 4.0, PI / 4.0, 3.0 * PI / 2.0), + ); + sphere_cell_list.insert( + 3, + Spherical::<4>::from_polar_coordinates(PI / 2.0, 3.0 * PI / 2.0 + 0.001, 0.0), + ); + sphere_cell_list.insert( + 4, + Spherical::<4>::from_polar_coordinates(PI / 4.0, 3.0 * PI / 2.0 + 0.001, 0.0), + ); + + assert!(sphere_cell_list.cell_index.get(&0) == Some(&CellIndex([0, 0, 0, 0]))); + assert!(sphere_cell_list.cell_index.get(&1) == Some(&CellIndex([0, 0, -1, 0]))); + assert!(sphere_cell_list.cell_index.get(&2) == Some(&CellIndex([0, -1, -1, 0]))); + assert!(sphere_cell_list.cell_index.get(&3) == Some(&CellIndex([-1, 0, 0, 0]))); + assert!(sphere_cell_list.cell_index.get(&4) == Some(&CellIndex([-1, 0, 0, 0]))); + + let keys = sphere_cell_list.get_keys(&[0, 0, 0, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&0)); + + let keys = sphere_cell_list.get_keys(&[-1, 0, 0, 0]); + assert!(keys.len() == 2); + check!(keys.contains(&3)); + check!(keys.contains(&4)); + + let keys = sphere_cell_list.get_keys(&[0, 0, -1, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&1)); + + let keys = sphere_cell_list.get_keys(&[0, -1, -1, 0]); + assert!(keys.len() == 1); + check!(keys.contains(&2)); + } + + #[rstest] + fn two_sphere_consistency() { + const N_STEPS: usize = 10_000; + let mut rng = StdRng::seed_from_u64(0); + let mut reference = FxHashMap::default(); + + let cell_width = 0.2; + let mut cell_list = SphericalVecCell::::builder() + .euclidean_nominal_search_radius( + cell_width + .try_into() + .expect("hard-coded cell with should be positive"), + ) + .build(); + let position_distribution = SphericalDisk { + disk_radius: (PI / 4.0).try_into().expect("hard-coded positive numbers"), + point: Spherical::<3>::from_polar_coordinates(0.0, 0.0), + }; + let key_distribution = + Uniform::new(0, N_STEPS / 4).expect("hardcoded distribution should be valid"); + + for _ in 0..N_STEPS { + // Add more keys than removing + if rng.random_bool(0.7) { + let position: Spherical<3> = position_distribution.sample(&mut rng); + let key = key_distribution.sample(&mut rng); + + cell_list.insert(key, position); + reference.insert(key, cell_list.cell_index_from_position(&position)); + } else { + let key = key_distribution.sample(&mut rng); + cell_list.remove(&key); + reference.remove(&key); + } + } + + // Validate that cell_index contains the expected keys and that + // keys_map is consistent. + assert!(cell_list.cell_index.len() == reference.len()); + for (reference_key, reference_value) in reference.drain() { + let value = cell_list.cell_index.get(&reference_key); + check!(value == Some(&CellIndex(reference_value))); + + let keys = cell_list.get_keys(&reference_value); + check!(keys.contains(&reference_key)); + } + + // Ensure that there are no extra values in keys_map. + let total = cell_list.keys_map.iter().map(Vec::len).sum(); + check!(cell_list.cell_index.len() == total); + check!(total > 300); + } + + #[rstest] + fn three_sphere_consistency() { + const N_STEPS: usize = 10_000; + let mut rng = StdRng::seed_from_u64(0); + let mut reference = FxHashMap::default(); + + let cell_width = 0.2; + let mut cell_list = SphericalVecCell::::builder() + .euclidean_nominal_search_radius( + cell_width + .try_into() + .expect("hard-coded cell with should be positive"), + ) + .build(); + let position_distribution = SphericalDisk { + disk_radius: (PI / 4.0).try_into().expect("hard-coded positive numbers"), + point: Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0), + }; + let key_distribution = + Uniform::new(0, N_STEPS / 4).expect("hardcoded distribution should be valid"); + + for _ in 0..N_STEPS { + // Add more keys than removing + if rng.random_bool(0.7) { + let position: Spherical<4> = position_distribution.sample(&mut rng); + let key = key_distribution.sample(&mut rng); + + cell_list.insert(key, position); + reference.insert(key, cell_list.cell_index_from_position(&position)); + } else { + let key = key_distribution.sample(&mut rng); + cell_list.remove(&key); + reference.remove(&key); + } + } + + // Validate that cell_index contains the expected keys and that + // keys_map is consistent. + assert!(cell_list.cell_index.len() == reference.len()); + for (reference_key, reference_value) in reference.drain() { + let value = cell_list.cell_index.get(&reference_key); + check!(value == Some(&CellIndex(reference_value))); + + let keys = cell_list.get_keys(&reference_value); + check!(keys.contains(&reference_key)); + } + + // Ensure that there are no extra values in keys_map. + let total = cell_list.keys_map.iter().map(Vec::len).sum(); + check!(cell_list.cell_index.len() == total); + check!(total > 300); + } +} From 3f2e0d3700542f3cc1ef0675d3c755a954dc583a Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Sun, 7 Jun 2026 19:16:38 -0400 Subject: [PATCH 06/16] docs for SphericalVecCell --- hoomd-mc/src/hypercuboid.rs | 17 +++-------------- hoomd-spatial/src/spherical_vec_cell.rs | 25 +++++++++++++++++++++++-- 2 files changed, 26 insertions(+), 16 deletions(-) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index 00d4f50e1..e495950b6 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -144,7 +144,7 @@ impl Checkerboard> for HypercuboidCheckerboard { let mut space_multi_index: [i64; N] = array::from_fn(|i| (p.coordinates[i] / self.space_width[i]).floor() as i64); - for (index, shape, periodic) in izip!(&mut space_multi_index, self.shape, self.periodic) { + for (index, shape) in izip!(&mut space_multi_index, self.shape) { // The origin is in the lower left corner of the box and may be up // to one space width to the left of simulation boundary. Therefore, // negative indices are out of bounds (and should never been seen @@ -153,18 +153,7 @@ impl Checkerboard> for HypercuboidCheckerboard { return None; } - if periodic { - // In periodic boundaries, the checkerboard spaces end before - // the right side. The space at the rightmost edge is identical - // to space 0 to make the checkerboard coloring commensurate - // with the periodic boundary conditions. - if *index as usize == shape { - *index = 0; - } - if *index as usize > shape { - return None; - } - } else if *index as usize >= shape { + if *index as usize >= shape { // In non-periodic boundaries, the checkerboard extends one full // space to the right of the simulation boundary so that when // it is shifted to the left it will still cover the boundary. @@ -253,7 +242,7 @@ impl HypercuboidCheckerboard { /// Partition the space indices by color. #[expect( clippy::todo, - reason = "there are no known use-cases for parallel 4D, 5D, ... simulations at this time" + reason = "there are no known use-cases for parallel 5D, 6D, ... simulations at this time" )] fn construct_space_indices_by_color(shape: [usize; N]) -> Vec> { for width in shape { diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index 323833599..e07fb4597 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -26,7 +26,18 @@ use super::{PointUpdate, PointsNearBall, WithSearchRadius}; use crate::{IndexFromPosition, hash_cell::CellIndex}; -/// TODO: docs +/// Implement [`VecCell`] for [`Spherical`] bodies. +/// +/// `Spherical` bodies exist on the surface of an $`(N-1)`$-sphere embedded +/// in $`N`$-dimensional Cartesian space. As such, methods for `VecCell` +/// can be adapted to work for `Spherical`. +/// +/// `SphericalVecCell` only differs from `VecCell` in that the user may choose +/// to use either the `Spherical` geodesic distance or the Euclidean +/// "line-of-site" distance when specifying nominal and maximum search radii. +/// +/// [`VecCell`]: hoomd_spatial::VecCell; +/// [`Spherical`]: hoomd_manifold::Spherical; #[serde_as] #[derive(Clone, Debug, Serialize, Deserialize)] pub struct SphericalVecCell @@ -239,7 +250,17 @@ impl Default for SphericalVecCell where K: Copy + Eq + Hash, { - /// TODO + /// Construct a default [`SphericalVecCell`]. + /// + /// The default sets both the nominal and maximum search to 1.0. + /// + /// # Example + /// + /// ``` + /// use hoomd_spatial::SphericalVecCell; + /// + /// let two_sphere_vec_cell = SphericalVecCell::::default(); + /// ``` #[inline] fn default() -> Self { Self::builder().build() From f1f27a75fb1a8df66cfbbf459f6fac765ff9824d Mon Sep 17 00:00:00 2001 From: mthran Date: Mon, 8 Jun 2026 14:06:44 -0400 Subject: [PATCH 07/16] add new `ClosedSpherical` boundary` --- hoomd-mc/src/hypercuboid.rs | 34 ++++++++++- hoomd-microstate/src/boundary.rs | 2 + .../src/boundary/closed_spherical.rs | 58 +++++++++++++++++++ hoomd-spatial/src/spherical_vec_cell.rs | 2 +- 4 files changed, 94 insertions(+), 2 deletions(-) create mode 100644 hoomd-microstate/src/boundary/closed_spherical.rs diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index e495950b6..a7f28f045 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -20,7 +20,7 @@ use std::array; use hoomd_geometry::shape::Hypercuboid; use hoomd_manifold::Spherical; -use hoomd_microstate::boundary::{Closed, Periodic}; +use hoomd_microstate::boundary::{Closed, ClosedSpherical, Periodic}; use hoomd_utility::valid::PositiveReal; use hoomd_vector::Cartesian; @@ -469,6 +469,38 @@ impl Cover> for Periodic> { } } +impl Cover> for ClosedSpherical { + type Checkerboard = HypercuboidCheckerboard; + #[inline] + fn cover( + &self, + rng: &mut R, + interaction_range: PositiveReal, + ) -> Self::Checkerboard { + HypercuboidCheckerboard::new( + rng, + interaction_range, + [2.0.try_into().expect("hard-coded positive number"); N], + [false; N], + ) + } + + #[inline] + fn cover_into( + &self, + checkerboard: &mut Self::Checkerboard, + rng: &mut R, + interaction_range: PositiveReal, + ) { + checkerboard.update( + rng, + interaction_range, + [2.0.try_into().expect("hard-coded positive number"); N], + [false; N], + ); + } +} + #[cfg(test)] mod tests { use assert2::{assert, check}; diff --git a/hoomd-microstate/src/boundary.rs b/hoomd-microstate/src/boundary.rs index 6b09d365b..bcfe5589e 100644 --- a/hoomd-microstate/src/boundary.rs +++ b/hoomd-microstate/src/boundary.rs @@ -24,10 +24,12 @@ use arrayvec::ArrayVec; use thiserror::Error; mod closed; +mod closed_spherical; mod open; mod periodic; pub use closed::Closed; +pub use closed_spherical::ClosedSpherical; pub use open::Open; pub use periodic::Periodic; diff --git a/hoomd-microstate/src/boundary/closed_spherical.rs b/hoomd-microstate/src/boundary/closed_spherical.rs new file mode 100644 index 000000000..e8a3ec096 --- /dev/null +++ b/hoomd-microstate/src/boundary/closed_spherical.rs @@ -0,0 +1,58 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! Implement `ClosedSpherical` + +use arrayvec::ArrayVec; + +use super::{Error, GenerateGhosts, MAX_GHOSTS, Wrap}; +use crate::property::Position; +use hoomd_manifold::Spherical; + +/// [`ClosedSpherical`] implements a hypercubic box enclosing a unit-radius +/// $`(N-1)`$-sphere. Use `ClosedSpherical` alongside [`SphericalVecCell`] to +/// implement [`ParallelSweep`] for [`Spherical`] bodies. `ClosedSpherical` is +/// otherwise functionally identical to using [`Open`] for [`Spherical`] +/// simulations. +/// +/// # Example +/// ``` +/// use hoomd_microstate::boundary::ClosedSpherical; +/// +/// # fn main() -> Result<(), Box> { +/// let closed_spherical: ClosedSpherical<3> = ClosedSpherical {}; +/// # Ok(()) +/// # } +/// ``` +/// +/// Similar to [`Closed`], `ClosedSpherical` does not wrap bodies and sites, +/// nor does it generate ghost sites. +/// +/// [`SphericalVecCell`]: hoomd_spatial::SphericalVecCell; +/// [`ParallelSweep`]: hoomd_mc::ParallelSweep; +/// [`Spherical`]: hoomd_manifold::Spherical; +pub struct ClosedSpherical {} + +impl Wrap for ClosedSpherical +where + BS: Position>, +{ + #[inline] + fn wrap(&self, properties: BS) -> Result { + Ok(properties) + } +} + +impl GenerateGhosts for ClosedSpherical +where + S: Default, +{ + #[inline] + fn maximum_interaction_range(&self) -> f64 { + std::f64::consts::PI + } + #[inline] + fn generate_ghosts(&self, _site_properties: &S) -> ArrayVec { + ArrayVec::new() + } +} diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index e07fb4597..843fce0d7 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -259,7 +259,7 @@ where /// ``` /// use hoomd_spatial::SphericalVecCell; /// - /// let two_sphere_vec_cell = SphericalVecCell::::default(); + /// let two_sphere_vec_cell = SphericalVecCell::::default(); /// ``` #[inline] fn default() -> Self { From af1f41e1d12d2ce6e2b4eb0bc82c91b663e7ec19 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Thu, 11 Jun 2026 12:47:49 -0400 Subject: [PATCH 08/16] remove duplicate functions from VecCell --- hoomd-spatial/src/spherical_vec_cell.rs | 104 +----------------------- hoomd-spatial/src/vec_cell.rs | 22 ++--- 2 files changed, 13 insertions(+), 113 deletions(-) diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index 843fce0d7..3cf640ac8 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -24,7 +24,7 @@ use hoomd_utility::valid::PositiveReal; use super::{PointUpdate, PointsNearBall, WithSearchRadius}; -use crate::{IndexFromPosition, hash_cell::CellIndex}; +use crate::{IndexFromPosition, hash_cell::CellIndex, vec_cell::{CellIndexIterator, generate_all_stencils, PointsIterator}}; /// Implement [`VecCell`] for [`Spherical`] bodies. /// @@ -61,82 +61,6 @@ where stencils: Vec>, } -/// Iterate over cell indices in row major order. -struct CellIndexIterator { - /// The current cell index. - cell_index: [i64; D], - /// The half extent of the cube. - half_extent: u32, -} - -impl CellIndexIterator { - /// Iterate over a cube. - /// - /// The cube extends from `[-half_extent, -half_extent, ..., -half_extent]` to - /// `[half_extent, half_extent, ..., half_extent]`. - fn cube(half_extent: u32) -> Self { - let mut cell_index = [-i64::from(half_extent); D]; - cell_index[D - 1] -= 1; - Self { - cell_index, - half_extent, - } - } - - /// Increment the cell index. - #[inline] - fn increment_cell_index(&mut self) -> Option<[i64; D]> { - self.cell_index[D - 1] += 1; - - for i in (0..D).rev() { - if self.cell_index[i] > self.half_extent.into() { - if i == 0 { - return None; - } - - self.cell_index[i] = -(i64::from(self.half_extent)); - self.cell_index[i - 1] += 1; - } - } - - Some(self.cell_index) - } -} - -impl Iterator for CellIndexIterator { - type Item = [i64; D]; - - fn next(&mut self) -> Option { - self.increment_cell_index() - } -} - -/// Generate the stencil for a given radius. -fn generate_stencil(radius: u32) -> Vec<[i64; D]> { - assert!(radius >= 1, "cell list must have a minimum radius of 1"); - - let mut result = CellIndexIterator::cube(radius).collect::>(); - result.sort_by(|a, b| { - let r_a = a.iter().map(|x| x.pow(2)); - let r_b = b.iter().map(|x| x.pow(2)); - r_a.cmp(r_b) - }); - result -} - -/// Generate the stencils up to a given radius. -pub(crate) fn generate_all_stencils(max_radius: u32) -> Vec> { - assert!(max_radius >= 1, "cell list must have a minimum radius of 1"); - let mut result = Vec::new(); - - // A cell list will never use radius=0, so stencil[i] stores the stencil needed for - // radius i+1. - for radius in 0..max_radius { - result.push(generate_stencil(radius + 1)); - } - result -} - /// Construct a [`SphericalVecCell`] with given parameters. /// /// # Example @@ -602,31 +526,7 @@ where } } -/// Iterate over keys in the cell list around a given center cell. -struct PointsIterator<'a, K, const D: usize> -where - K: Eq + Hash, -{ - /// Keys of the current cell iteration (None if the cell is empty) - keys: Option<&'a Vec>, - - /// The cell list we are iterating in. - cell_list: &'a SphericalVecCell, - - /// Current location of the iteration in the cell. - index_in_current_cell: usize, - - /// Current location of the iteration in the stencil. - current_stencil: usize, - - /// Cell offsets to iterate over. - stencil: &'a [[i64; D]], - - /// The cell at the center of the iteration. - center: [i64; D], -} - -impl Iterator for PointsIterator<'_, K, D> +impl Iterator for PointsIterator<'_, K, D, SphericalVecCell> where K: Copy + Eq + Hash, { diff --git a/hoomd-spatial/src/vec_cell.rs b/hoomd-spatial/src/vec_cell.rs index 72f4a8ea7..11cdf1ed3 100644 --- a/hoomd-spatial/src/vec_cell.rs +++ b/hoomd-spatial/src/vec_cell.rs @@ -114,7 +114,7 @@ where } /// Iterate over cell indices in row major order. -struct CellIndexIterator { +pub(crate) struct CellIndexIterator { /// The current cell index. cell_index: [i64; D], /// The half extent of the cube. @@ -126,7 +126,7 @@ impl CellIndexIterator { /// /// The cube extends from `[-half_extent, -half_extent, ..., -half_extent]` to /// `[half_extent, half_extent, ..., half_extent]`. - fn cube(half_extent: u32) -> Self { + pub(crate) fn cube(half_extent: u32) -> Self { let mut cell_index = [-i64::from(half_extent); D]; cell_index[D - 1] -= 1; Self { @@ -164,7 +164,7 @@ impl Iterator for CellIndexIterator { } /// Generate the stencil for a given radius. -fn generate_stencil(radius: u32) -> Vec<[i64; D]> { +pub fn generate_stencil(radius: u32) -> Vec<[i64; D]> { assert!(radius >= 1, "cell list must have a minimum radius of 1"); let mut result = CellIndexIterator::cube(radius).collect::>(); @@ -617,30 +617,30 @@ where } /// Iterate over keys in the cell list around a given center cell. -struct PointsIterator<'a, K, const D: usize> +pub(crate) struct PointsIterator<'a, K, const D: usize, X> where K: Eq + Hash, { /// Keys of the current cell iteration (None if the cell is empty) - keys: Option<&'a Vec>, + pub(crate) keys: Option<&'a Vec>, /// The cell list we are iterating in. - cell_list: &'a VecCell, + pub(crate) cell_list: &'a X, //VecCell, /// Current location of the iteration in the cell. - index_in_current_cell: usize, + pub(crate) index_in_current_cell: usize, /// Current location of the iteration in the stencil. - current_stencil: usize, + pub(crate) current_stencil: usize, /// Cell offsets to iterate over. - stencil: &'a [[i64; D]], + pub(crate) stencil: &'a [[i64; D]], /// The cell at the center of the iteration. - center: [i64; D], + pub(crate) center: [i64; D], } -impl Iterator for PointsIterator<'_, K, D> +impl Iterator for PointsIterator<'_, K, D, VecCell> where K: Copy + Eq + Hash, { From eafd4ea4bb16bebd95fa48a59b15378d5f5abcc8 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Thu, 11 Jun 2026 13:05:58 -0400 Subject: [PATCH 09/16] docs --- hoomd-microstate/src/boundary/closed_spherical.rs | 5 ++++- hoomd-spatial/src/spherical_vec_cell.rs | 8 ++++++-- hoomd-spatial/src/vec_cell.rs | 4 ++-- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/hoomd-microstate/src/boundary/closed_spherical.rs b/hoomd-microstate/src/boundary/closed_spherical.rs index e8a3ec096..ddcf4c4c9 100644 --- a/hoomd-microstate/src/boundary/closed_spherical.rs +++ b/hoomd-microstate/src/boundary/closed_spherical.rs @@ -26,7 +26,10 @@ use hoomd_manifold::Spherical; /// ``` /// /// Similar to [`Closed`], `ClosedSpherical` does not wrap bodies and sites, -/// nor does it generate ghost sites. +/// nor does it generate ghost sites. However, unlike [`Closed`], +/// `ClosedSpherical` does not implement `IsPointInside`, `MapPoint`, +/// `Scale`, or `Volume`, as these methods are unrelated to the boundary +/// conditions for `Spherical` simulations. /// /// [`SphericalVecCell`]: hoomd_spatial::SphericalVecCell; /// [`ParallelSweep`]: hoomd_mc::ParallelSweep; diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index 3cf640ac8..5433a4764 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -24,7 +24,11 @@ use hoomd_utility::valid::PositiveReal; use super::{PointUpdate, PointsNearBall, WithSearchRadius}; -use crate::{IndexFromPosition, hash_cell::CellIndex, vec_cell::{CellIndexIterator, generate_all_stencils, PointsIterator}}; +use crate::{ + IndexFromPosition, + hash_cell::CellIndex, + vec_cell::{CellIndexIterator, PointsIterator, generate_all_stencils}, +}; /// Implement [`VecCell`] for [`Spherical`] bodies. /// @@ -526,7 +530,7 @@ where } } -impl Iterator for PointsIterator<'_, K, D, SphericalVecCell> +impl Iterator for PointsIterator<'_, K, D, SphericalVecCell> where K: Copy + Eq + Hash, { diff --git a/hoomd-spatial/src/vec_cell.rs b/hoomd-spatial/src/vec_cell.rs index 11cdf1ed3..4877f24ae 100644 --- a/hoomd-spatial/src/vec_cell.rs +++ b/hoomd-spatial/src/vec_cell.rs @@ -625,7 +625,7 @@ where pub(crate) keys: Option<&'a Vec>, /// The cell list we are iterating in. - pub(crate) cell_list: &'a X, //VecCell, + pub(crate) cell_list: &'a X, // VecCell, /// Current location of the iteration in the cell. pub(crate) index_in_current_cell: usize, @@ -640,7 +640,7 @@ where pub(crate) center: [i64; D], } -impl Iterator for PointsIterator<'_, K, D, VecCell> +impl Iterator for PointsIterator<'_, K, D, VecCell> where K: Copy + Eq + Hash, { From 03b5b95874fff0644e5c78d11910b4158cc25b7c Mon Sep 17 00:00:00 2001 From: mthran Date: Mon, 15 Jun 2026 16:40:09 -0400 Subject: [PATCH 10/16] add serde derive for ClosedSpherical --- hoomd-microstate/src/boundary/closed_spherical.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hoomd-microstate/src/boundary/closed_spherical.rs b/hoomd-microstate/src/boundary/closed_spherical.rs index ddcf4c4c9..65730b43d 100644 --- a/hoomd-microstate/src/boundary/closed_spherical.rs +++ b/hoomd-microstate/src/boundary/closed_spherical.rs @@ -4,6 +4,7 @@ //! Implement `ClosedSpherical` use arrayvec::ArrayVec; +use serde::{Deserialize, Serialize}; use super::{Error, GenerateGhosts, MAX_GHOSTS, Wrap}; use crate::property::Position; @@ -34,6 +35,7 @@ use hoomd_manifold::Spherical; /// [`SphericalVecCell`]: hoomd_spatial::SphericalVecCell; /// [`ParallelSweep`]: hoomd_mc::ParallelSweep; /// [`Spherical`]: hoomd_manifold::Spherical; +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] pub struct ClosedSpherical {} impl Wrap for ClosedSpherical From 07e1cc4174c7161383ab161f21ae7f518d88b246 Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Wed, 24 Jun 2026 15:42:42 -0400 Subject: [PATCH 11/16] cargo fmt --- hoomd-mc/src/hypercuboid.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index a0026d8b0..e6fb28143 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -268,7 +268,6 @@ impl HypercuboidCheckerboard { }) .collect(); - (0..num_colors) .map(|color| { let offset = Self::index_to_multi_index(color, [2; N]); From 8d9a10b64f6bb84c4976484d7d407a1c8a184ede Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Wed, 24 Jun 2026 15:54:51 -0400 Subject: [PATCH 12/16] Avoid code duplication for HypercuboidCheckerboard --- hoomd-mc/src/hypercuboid.rs | 32 +++----------------------------- 1 file changed, 3 insertions(+), 29 deletions(-) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index e6fb28143..c868fbfe8 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -140,43 +140,17 @@ impl Checkerboard> for HypercuboidCheckerboard { impl Checkerboard> for HypercuboidCheckerboard { #[inline] fn point_to_space_index(&self, point: &Spherical) -> Option { - let p = *point.point() - self.origin; - let mut space_multi_index: [i64; N] = - array::from_fn(|i| (p.coordinates[i] / self.space_width[i]).floor() as i64); - - for (index, shape) in izip!(&mut space_multi_index, self.shape) { - // The origin is in the lower left corner of the box and may be up - // to one space width to the left of simulation boundary. Therefore, - // negative indices are out of bounds (and should never been seen - // for wrapped inputs). - if *index < 0 { - return None; - } - - if *index as usize >= shape { - // In non-periodic boundaries, the checkerboard extends one full - // space to the right of the simulation boundary so that when - // it is shifted to the left it will still cover the boundary. - // Therefore, any points outside the checkerboard shape are - // invalid. - return None; - } - } - - Some(Self::multi_index_to_index( - array::from_fn(|i| space_multi_index[i] as usize), - self.shape, - )) + >>::point_to_space_index(self, &Cartesian::from(*point.coordinates())) } #[inline] fn space_indices_by_color(&self) -> &[Vec] { - &self.space_indices_by_color + >>::space_indices_by_color(self) } #[inline] fn num_spaces(&self) -> usize { - self.shape.iter().product() + >>::num_spaces(self) } } From 3d42bff586abf2bfa751a2c24a606c5cb1cd41ea Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Wed, 24 Jun 2026 16:20:28 -0400 Subject: [PATCH 13/16] cargo fmt --- hoomd-mc/src/hypercuboid.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index c868fbfe8..dff841abd 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -140,7 +140,10 @@ impl Checkerboard> for HypercuboidCheckerboard { impl Checkerboard> for HypercuboidCheckerboard { #[inline] fn point_to_space_index(&self, point: &Spherical) -> Option { - >>::point_to_space_index(self, &Cartesian::from(*point.coordinates())) + >>::point_to_space_index( + self, + &Cartesian::from(*point.coordinates()), + ) } #[inline] From c5efd6c1a231cd23507183a248cf7a899b79dd6f Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 26 Jun 2026 14:44:45 -0400 Subject: [PATCH 14/16] avoid duplicate code in spherical_vec_cell.rs --- hoomd-spatial/src/spherical_vec_cell.rs | 346 +++++------------------- hoomd-spatial/src/vec_cell.rs | 22 +- 2 files changed, 82 insertions(+), 286 deletions(-) diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index 5433a4764..211d262fa 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -14,10 +14,7 @@ use serde::{Deserialize, Serialize}; use serde_with::serde_as; -use std::{array, cmp::Eq, f64::consts::PI, fmt, hash::Hash, iter, marker::PhantomData, mem}; - -use log::trace; -use rustc_hash::FxHashMap; +use std::{array, cmp::Eq, f64::consts::PI, fmt, hash::Hash}; use hoomd_manifold::Spherical; use hoomd_utility::valid::PositiveReal; @@ -26,8 +23,7 @@ use super::{PointUpdate, PointsNearBall, WithSearchRadius}; use crate::{ IndexFromPosition, - hash_cell::CellIndex, - vec_cell::{CellIndexIterator, PointsIterator, generate_all_stencils}, + vec_cell::{PointsIterator, VecCell, VecCellBuilder}, }; /// Implement [`VecCell`] for [`Spherical`] bodies. @@ -44,26 +40,7 @@ use crate::{ /// [`Spherical`]: hoomd_manifold::Spherical; #[serde_as] #[derive(Clone, Debug, Serialize, Deserialize)] -pub struct SphericalVecCell -where - K: Eq + Hash, -{ - /// The Euclidean width of each cell. - cell_width: PositiveReal, - - /// A map from cell indices to cell contents. - keys_map: Vec>, - - /// A map from particle indices to cell indices. - cell_index: FxHashMap>, - - /// The shape of `keys_map` is `(half_extent * 2 + 1).powi(D)`. - half_extent: u32, - - /// Pre-computed stencils. - #[serde_as(as = "Vec>")] - stencils: Vec>, -} +pub struct SphericalVecCell(pub VecCell); /// Construct a [`SphericalVecCell`] with given parameters. /// @@ -82,16 +59,7 @@ where /// # Ok(()) /// # } /// ``` -pub struct SphericalVecCellBuilder { - /// Most commonly used search radius, in Euclidean metruic. - nominal_search_radius: PositiveReal, - - /// Largest possible search radius, in Euclidean metric. - maximum_search_radius: f64, - - /// Track the key type. - phantom_key: PhantomData, -} +pub struct SphericalVecCellBuilder(pub VecCellBuilder); impl SphericalVecCellBuilder where @@ -111,7 +79,7 @@ where mut self, spherical_nominal_search_radius: PositiveReal, ) -> Self { - self.nominal_search_radius = (spherical_nominal_search_radius.get().clamp(0.0, PI / 2.0)) + self.0.nominal_search_radius = (spherical_nominal_search_radius.get().clamp(0.0, PI / 2.0)) .sin() .try_into() .expect("clamp ensures number is positive"); @@ -127,7 +95,7 @@ where #[inline] #[must_use] pub fn spherical_maximum_search_radius(mut self, spherical_maximum_search_radius: f64) -> Self { - self.maximum_search_radius = (spherical_maximum_search_radius.clamp(0.0, PI / 2.0)).sin(); + self.0.maximum_search_radius = (spherical_maximum_search_radius.clamp(0.0, PI / 2.0)).sin(); self } @@ -138,7 +106,7 @@ where mut self, euclidean_nominal_search_radius: PositiveReal, ) -> Self { - self.nominal_search_radius = euclidean_nominal_search_radius; + self.0.nominal_search_radius = euclidean_nominal_search_radius; self } @@ -151,7 +119,7 @@ where #[inline] #[must_use] pub fn euclidean_maximum_search_radius(mut self, euclidean_maximum_search_radius: f64) -> Self { - self.maximum_search_radius = euclidean_maximum_search_radius; + self.0.maximum_search_radius = euclidean_maximum_search_radius; self } @@ -159,18 +127,7 @@ where #[inline] #[must_use] pub fn build(self) -> SphericalVecCell { - let maximum_stencil_radius = - (self.maximum_search_radius / self.nominal_search_radius.get()).ceil() as u32; - let half_extent: u32 = 1; - - SphericalVecCell { - cell_width: self.nominal_search_radius, - keys_map: iter::repeat_n(Vec::new(), (half_extent * 2 + 1).pow(D as u32) as usize) - .collect(), - cell_index: FxHashMap::default(), - half_extent, - stencils: generate_all_stencils(maximum_stencil_radius.max(1)), - } + SphericalVecCell(VecCellBuilder::::build(self.0)) } } @@ -235,44 +192,11 @@ impl SphericalVecCell where K: Eq + Hash, { - /// Compute the cell index given a `Spherical` position in space. - #[inline] - fn cell_index_from_position(&self, position: &Spherical) -> [i64; D] { - std::array::from_fn(|j| (position.coordinates()[j] / self.cell_width.get()).floor() as i64) - } - - /// Compute the vector index from a cell index - /// - /// Returns `None` when the index is out of bounds. - #[inline] - fn map_index_from_cell(half_extent: u32, cell_index: &[i64; D]) -> Option { - assert!(D > 1); - - let mut vec_index: usize = 0; - let mut width = 1; - - for i in (0..D).rev() { - let needed_extent = cell_index[i].unsigned_abs(); - if needed_extent > u64::from(half_extent) { - return None; - } - let v: usize = (cell_index[i] + i64::from(half_extent)) - .try_into() - .expect("cell index should be in bounds"); - - vec_index += v * width; - width *= (half_extent * 2 + 1) as usize; - } - Some(vec_index) - } - /// Get the keys in a given cell index #[cfg(test)] #[inline] fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { - let index = Self::map_index_from_cell(self.half_extent, cell_index) - .expect("cell_index should be in bounds"); - &self.keys_map[index] + VecCell::::get_keys(&self.0, cell_index) } } @@ -305,13 +229,7 @@ where #[inline] #[must_use] pub fn builder() -> SphericalVecCellBuilder { - SphericalVecCellBuilder { - nominal_search_radius: 1.0 - .try_into() - .expect("hard-coded constant is a positive real"), - maximum_search_radius: 1.0, - phantom_key: PhantomData, - } + SphericalVecCellBuilder(VecCell::::builder()) } /// Remove excess capacity from dynamically allocated arrays. @@ -321,42 +239,7 @@ where /// the range spanned by points currently in the data structure. #[inline] pub fn shrink_to_fit(&mut self) { - for keys in &mut self.keys_map { - keys.shrink_to_fit(); - } - self.keys_map.shrink_to_fit(); - self.cell_index.shrink_to_fit(); - } - - /// Double the number of cells stored along each axis until it includes the target. - fn expand_to(&mut self, target: u32) { - if self.half_extent >= target { - return; - } - - let mut new_half_extent = self.half_extent.min(1) * 2; - - while new_half_extent < target { - new_half_extent *= 2; - } - - trace!("Expanding to {}^{} cells", new_half_extent * 2 + 1, D); - - let mut new_keys_map: Vec> = - iter::repeat_n(Vec::new(), (new_half_extent * 2 + 1).pow(D as u32) as usize).collect(); - let old_half_extent = self.half_extent; - let old_keys_map = &mut self.keys_map; - - for old_cell_index in CellIndexIterator::cube(old_half_extent) { - let old_vec_index = Self::map_index_from_cell(old_half_extent, &old_cell_index) - .expect("cell_index should be consistent with keys_map"); - let new_vec_index = Self::map_index_from_cell(new_half_extent, &old_cell_index) - .expect("old_cell_index should be inside the new keys_map"); - new_keys_map[new_vec_index] = mem::take(&mut old_keys_map[old_vec_index]); - } - - self.half_extent = new_half_extent; - self.keys_map = new_keys_map; + VecCell::::shrink_to_fit(&mut self.0); } } @@ -385,39 +268,7 @@ where /// ``` #[inline] fn insert(&mut self, key: K, position: Spherical) { - let cell_index = self.cell_index_from_position(&position); - let old_cell_index = self.cell_index.insert(key, CellIndex(cell_index)); - let map_index = - Self::map_index_from_cell(self.half_extent, &cell_index).unwrap_or_else(|| { - let max_half_extent = cell_index - .iter() - .map(|x| x.unsigned_abs()) - .reduce(u64::max) - .expect("D should be greater than 1"); - self.expand_to( - max_half_extent - .try_into() - .expect("max extent cannot exceed u32::MAX"), - ); - Self::map_index_from_cell(self.half_extent, &cell_index) - .expect("cell_index should be in the expanded VecCell") - }); - - // This checks if old_cell_index is None or if it is different from the new cell index. - if old_cell_index != Some(CellIndex(cell_index)) { - // Add the particle index to the new cell index vector. - self.keys_map[map_index].push(key); - - if let Some(old_cell_index) = old_cell_index { - // If the particle was in a different cell, we need to remove it from the old cell. - let old_map_index = Self::map_index_from_cell(self.half_extent, &old_cell_index.0) - .expect("cell_index and keys_map should agree"); - let old_keys = &mut self.keys_map[old_map_index]; - if let Some(pos) = old_keys.iter().position(|x| *x == key) { - old_keys.swap_remove(pos); - } - } - } + VecCell::::insert(&mut self.0, key, *position.point()) } /// Remove the point with the given key. @@ -438,16 +289,7 @@ where /// ``` #[inline] fn remove(&mut self, key: &K) { - let cell_index = self.cell_index.remove(key); - if let Some(cell_index) = cell_index { - let map_index = Self::map_index_from_cell(self.half_extent, &cell_index.0); - if let Some(map_index) = map_index { - let keys = &mut self.keys_map[map_index]; - if let Some(idx) = keys.iter().position(|x| x == key) { - keys.swap_remove(idx); - } - } - } + VecCell::::remove(&mut self.0, key); } /// Get the number of points in the spatial data structure. @@ -468,7 +310,7 @@ where /// ``` #[inline] fn len(&self) -> usize { - self.cell_index.len() + self.0.cell_index.len() } /// Test if the spatial data structure is empty. @@ -483,7 +325,7 @@ where /// ``` #[inline] fn is_empty(&self) -> bool { - self.cell_index.is_empty() + self.0.cell_index.is_empty() } /// Test if the spatial data structure contains a key. @@ -502,7 +344,7 @@ where /// ``` #[inline] fn contains_key(&self, key: &K) -> bool { - self.cell_index.contains_key(key) + self.0.cell_index.contains_key(key) } /// Remove all points. @@ -523,10 +365,7 @@ where /// ``` #[inline] fn clear(&mut self) { - self.cell_index.clear(); - for keys in &mut self.keys_map { - keys.clear(); - } + self.0.clear() } } @@ -556,11 +395,9 @@ where let cell_index = array::from_fn(|i| self.center[i] + self.stencil[self.current_stencil][i]); - let map_index = SphericalVecCell::::map_index_from_cell( - self.cell_list.half_extent, - &cell_index, - ); - self.keys = map_index.map(|index| &self.cell_list.keys_map[index]); + let map_index = + VecCell::::map_index_from_cell(self.cell_list.0.half_extent, &cell_index); + self.keys = map_index.map(|index| &self.cell_list.0.keys_map[index]); } } } @@ -614,27 +451,7 @@ where fn points_near_ball(&self, position: &Spherical, radius: f64) -> impl Iterator { // convert spherical distance to Euclidean distance let euclidean_radius = radius.asin(); - let stencil_index = (euclidean_radius / self.cell_width.get()).ceil() as usize - 1; - assert!( - stencil_index < self.stencils.len(), - "search radius must be less than or equal to the maximum search radius" - ); - - let center = self.cell_index_from_position(position); - let stencil = &self.stencils[stencil_index]; - let map_index = Self::map_index_from_cell( - self.half_extent, - &array::from_fn(|i| center[i] + stencil[0][i]), - ); - - PointsIterator { - keys: map_index.map(|index| &self.keys_map[index]), - cell_list: self, - index_in_current_cell: 0, - current_stencil: 0, - stencil, - center, - } + VecCell::::points_near_ball(&self.0, position.point(), euclidean_radius) } } @@ -662,23 +479,7 @@ where reason = "no need to inline display" )] fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let largest_cell_size = self.keys_map.iter().map(Vec::len).fold(0, usize::max); - - writeln!(f, "VecCell:")?; - writeln!( - f, - "- {} total cells with {} cells on each side.", - self.keys_map.len(), - self.half_extent * 2 + 1 - )?; - writeln!(f, "- {} points.", self.cell_index.len())?; - writeln!( - f, - "- Nominal, maximum search radii: {}, {}", - self.cell_width, - self.cell_width.get() * self.stencils.len() as f64 - )?; - write!(f, "- Largest cell length: {largest_cell_size}") + VecCell::::fmt(&self.0, f) } } @@ -690,13 +491,14 @@ where #[inline] fn location_from_position(&self, position: &Spherical) -> Self::Location { - self.cell_index_from_position(position) + self.0.cell_index_from_position(position.point()) } } #[cfg(test)] mod tests { use super::*; + use crate::hash_cell::CellIndex; use hoomd_manifold::{Spherical, SphericalDisk}; use rand::{ RngExt, SeedableRng, @@ -707,6 +509,7 @@ mod tests { use approxim::assert_relative_eq; use assert2::{assert, check}; + use rustc_hash::FxHashMap; use std::f64::consts::PI; #[test] @@ -715,20 +518,19 @@ mod tests { .euclidean_nominal_search_radius(0.5.try_into().expect("hard-coded positive number")) .build(); check!( - spherical_cell_list - .cell_index_from_position(&Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0)) - == [2, 0, 0] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0).point() + ) == [2, 0, 0] ); check!( - spherical_cell_list.cell_index_from_position(&Spherical::<3>::from_polar_coordinates( - PI / 2.0, - PI / 2.0 - )) == [0, 2, 0] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<3>::from_polar_coordinates(PI / 2.0, PI / 2.0).point() + ) == [0, 2, 0] ); check!( - spherical_cell_list - .cell_index_from_position(&Spherical::<3>::from_polar_coordinates(0.0, 0.0)) - == [0, 0, 2] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<3>::from_polar_coordinates(0.0, 0.0).point() + ) == [0, 0, 2] ); } @@ -738,30 +540,24 @@ mod tests { .euclidean_nominal_search_radius(0.5.try_into().expect("hard-coded positive number")) .build(); check!( - spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( - PI / 2.0, - 0.0, - 0.0 - )) == [2, 0, 0, 0] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<4>::from_polar_coordinates(PI / 2.0, 0.0, 0.0).point() + ) == [2, 0, 0, 0] ); check!( - spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( - PI / 2.0, - PI / 2.0, - 0.0 - )) == [0, 2, 0, 0] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, 0.0).point() + ) == [0, 2, 0, 0] ); check!( - spherical_cell_list.cell_index_from_position(&Spherical::<4>::from_polar_coordinates( - PI / 2.0, - PI / 2.0, - PI / 2.0 - )) == [0, 0, 2, 0] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, PI / 2.0).point() + ) == [0, 0, 2, 0] ); check!( - spherical_cell_list - .cell_index_from_position(&Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0)) - == [0, 0, 0, 2] + spherical_cell_list.0.cell_index_from_position( + &Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0).point() + ) == [0, 0, 0, 2] ); } @@ -770,13 +566,13 @@ mod tests { let spherical_cell_list = SphericalVecCell::::with_search_radius( (PI / 2.0).try_into().expect("hard-coded positive number"), ); - assert!(spherical_cell_list.cell_width.get() == 1.0); + assert!(spherical_cell_list.0.cell_width.get() == 1.0); let spherical_cell_list = SphericalVecCell::::with_search_radius( (PI / 4.0).try_into().expect("hard-coded positive number"), ); assert_relative_eq!( - spherical_cell_list.cell_width.get(), + spherical_cell_list.0.cell_width.get(), (0.5_f64).sqrt(), epsilon = 1e-12 ); @@ -787,13 +583,13 @@ mod tests { let spherical_cell_list = SphericalVecCell::::with_search_radius( (PI / 2.0).try_into().expect("hard-coded positive number"), ); - assert!(spherical_cell_list.cell_width.get() == 1.0); + assert!(spherical_cell_list.0.cell_width.get() == 1.0); let spherical_cell_list = SphericalVecCell::::with_search_radius( (PI / 4.0).try_into().expect("hard-coded positive number"), ); assert_relative_eq!( - spherical_cell_list.cell_width.get(), + spherical_cell_list.0.cell_width.get(), (0.5_f64).sqrt(), epsilon = 1e-12 ); @@ -818,11 +614,11 @@ mod tests { Spherical::<3>::from_polar_coordinates(PI / 4.0, 3.0 * PI / 2.0 + 0.001), ); - assert!(sphere_cell_list.cell_index.get(&0) == Some(&CellIndex([1, 0, 0]))); - assert!(sphere_cell_list.cell_index.get(&1) == Some(&CellIndex([-1, 0, 0]))); - assert!(sphere_cell_list.cell_index.get(&2) == Some(&CellIndex([0, 0, 0]))); - assert!(sphere_cell_list.cell_index.get(&3) == Some(&CellIndex([0, -1, 0]))); - assert!(sphere_cell_list.cell_index.get(&4) == Some(&CellIndex([0, -1, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&0) == Some(&CellIndex([1, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&1) == Some(&CellIndex([-1, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&2) == Some(&CellIndex([0, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&3) == Some(&CellIndex([0, -1, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&4) == Some(&CellIndex([0, -1, 0]))); let keys = sphere_cell_list.get_keys(&[0, 0, 0]); assert!(keys.len() == 1); @@ -871,11 +667,11 @@ mod tests { Spherical::<4>::from_polar_coordinates(PI / 4.0, 3.0 * PI / 2.0 + 0.001, 0.0), ); - assert!(sphere_cell_list.cell_index.get(&0) == Some(&CellIndex([0, 0, 0, 0]))); - assert!(sphere_cell_list.cell_index.get(&1) == Some(&CellIndex([0, 0, -1, 0]))); - assert!(sphere_cell_list.cell_index.get(&2) == Some(&CellIndex([0, -1, -1, 0]))); - assert!(sphere_cell_list.cell_index.get(&3) == Some(&CellIndex([-1, 0, 0, 0]))); - assert!(sphere_cell_list.cell_index.get(&4) == Some(&CellIndex([-1, 0, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&0) == Some(&CellIndex([0, 0, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&1) == Some(&CellIndex([0, 0, -1, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&2) == Some(&CellIndex([0, -1, -1, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&3) == Some(&CellIndex([-1, 0, 0, 0]))); + assert!(sphere_cell_list.0.cell_index.get(&4) == Some(&CellIndex([-1, 0, 0, 0]))); let keys = sphere_cell_list.get_keys(&[0, 0, 0, 0]); assert!(keys.len() == 1); @@ -923,7 +719,7 @@ mod tests { let key = key_distribution.sample(&mut rng); cell_list.insert(key, position); - reference.insert(key, cell_list.cell_index_from_position(&position)); + reference.insert(key, cell_list.0.cell_index_from_position(&position.point())); } else { let key = key_distribution.sample(&mut rng); cell_list.remove(&key); @@ -933,9 +729,9 @@ mod tests { // Validate that cell_index contains the expected keys and that // keys_map is consistent. - assert!(cell_list.cell_index.len() == reference.len()); + assert!(cell_list.0.cell_index.len() == reference.len()); for (reference_key, reference_value) in reference.drain() { - let value = cell_list.cell_index.get(&reference_key); + let value = cell_list.0.cell_index.get(&reference_key); check!(value == Some(&CellIndex(reference_value))); let keys = cell_list.get_keys(&reference_value); @@ -943,8 +739,8 @@ mod tests { } // Ensure that there are no extra values in keys_map. - let total = cell_list.keys_map.iter().map(Vec::len).sum(); - check!(cell_list.cell_index.len() == total); + let total = cell_list.0.keys_map.iter().map(Vec::len).sum(); + check!(cell_list.0.cell_index.len() == total); check!(total > 300); } @@ -976,7 +772,7 @@ mod tests { let key = key_distribution.sample(&mut rng); cell_list.insert(key, position); - reference.insert(key, cell_list.cell_index_from_position(&position)); + reference.insert(key, cell_list.0.cell_index_from_position(&position.point())); } else { let key = key_distribution.sample(&mut rng); cell_list.remove(&key); @@ -986,9 +782,9 @@ mod tests { // Validate that cell_index contains the expected keys and that // keys_map is consistent. - assert!(cell_list.cell_index.len() == reference.len()); + assert!(cell_list.0.cell_index.len() == reference.len()); for (reference_key, reference_value) in reference.drain() { - let value = cell_list.cell_index.get(&reference_key); + let value = cell_list.0.cell_index.get(&reference_key); check!(value == Some(&CellIndex(reference_value))); let keys = cell_list.get_keys(&reference_value); @@ -996,8 +792,8 @@ mod tests { } // Ensure that there are no extra values in keys_map. - let total = cell_list.keys_map.iter().map(Vec::len).sum(); - check!(cell_list.cell_index.len() == total); + let total = cell_list.0.keys_map.iter().map(Vec::len).sum(); + check!(cell_list.0.cell_index.len() == total); check!(total > 300); } } diff --git a/hoomd-spatial/src/vec_cell.rs b/hoomd-spatial/src/vec_cell.rs index 4877f24ae..7c8b8d42d 100644 --- a/hoomd-spatial/src/vec_cell.rs +++ b/hoomd-spatial/src/vec_cell.rs @@ -97,16 +97,16 @@ where K: Eq + Hash, { /// The width of each cell. - cell_width: PositiveReal, + pub(crate) cell_width: PositiveReal, /// A map from cell indices to cell contents. - keys_map: Vec>, + pub(crate) keys_map: Vec>, /// A map from particle indices to cell indices. - cell_index: FxHashMap>, + pub(crate) cell_index: FxHashMap>, /// The shape of `keys_map` is `(half_extent * 2 + 1).powi(D)`. - half_extent: u32, + pub(crate) half_extent: u32, /// Pre-computed stencils. #[serde_as(as = "Vec>")] @@ -206,13 +206,13 @@ pub(crate) fn generate_all_stencils(max_radius: u32) -> Vec { /// Most commonly used search radius. - nominal_search_radius: PositiveReal, + pub(crate) nominal_search_radius: PositiveReal, /// Largest possible search radius. - maximum_search_radius: f64, + pub(crate) maximum_search_radius: f64, /// Track the key type. - phantom_key: PhantomData, + pub(crate) phantom_key: PhantomData, } impl VecCellBuilder @@ -353,7 +353,7 @@ where { /// Compute the cell index given a position in space. #[inline] - fn cell_index_from_position(&self, position: &Cartesian) -> [i64; D] { + pub(crate) fn cell_index_from_position(&self, position: &Cartesian) -> [i64; D] { std::array::from_fn(|j| (position.coordinates[j] / self.cell_width.get()).floor() as i64) } @@ -361,7 +361,7 @@ where /// /// Returns `None` when the index is out of bounds. #[inline] - fn map_index_from_cell(half_extent: u32, cell_index: &[i64; D]) -> Option { + pub(crate) fn map_index_from_cell(half_extent: u32, cell_index: &[i64; D]) -> Option { assert!(D > 1); let mut vec_index: usize = 0; @@ -385,7 +385,7 @@ where /// Get the keys in a given cell index #[cfg(test)] #[inline] - fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { + pub(crate) fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { let index = Self::map_index_from_cell(self.half_extent, cell_index) .expect("cell_index should be in bounds"); &self.keys_map[index] @@ -625,7 +625,7 @@ where pub(crate) keys: Option<&'a Vec>, /// The cell list we are iterating in. - pub(crate) cell_list: &'a X, // VecCell, + pub(crate) cell_list: &'a X, /// Current location of the iteration in the cell. pub(crate) index_in_current_cell: usize, From a5de993fd9f72f5476bf5f33e19f311e669deca8 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 26 Jun 2026 15:01:24 -0400 Subject: [PATCH 15/16] renamed ClosedSpherical to OpenSpherical --- hoomd-mc/src/hypercuboid.rs | 4 +-- hoomd-microstate/src/boundary.rs | 4 +-- ...{closed_spherical.rs => open_spherical.rs} | 25 ++++++++----------- 3 files changed, 15 insertions(+), 18 deletions(-) rename hoomd-microstate/src/boundary/{closed_spherical.rs => open_spherical.rs} (56%) diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index dff841abd..970146a96 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -20,7 +20,7 @@ use std::array; use hoomd_geometry::shape::Hypercuboid; use hoomd_manifold::Spherical; -use hoomd_microstate::boundary::{Closed, ClosedSpherical, Periodic}; +use hoomd_microstate::boundary::{Closed, OpenSpherical, Periodic}; use hoomd_utility::valid::PositiveReal; use hoomd_vector::Cartesian; @@ -387,7 +387,7 @@ impl Cover> for Periodic> { } } -impl Cover> for ClosedSpherical { +impl Cover> for OpenSpherical { type Checkerboard = HypercuboidCheckerboard; #[inline] fn cover( diff --git a/hoomd-microstate/src/boundary.rs b/hoomd-microstate/src/boundary.rs index bcfe5589e..451e5e29d 100644 --- a/hoomd-microstate/src/boundary.rs +++ b/hoomd-microstate/src/boundary.rs @@ -24,13 +24,13 @@ use arrayvec::ArrayVec; use thiserror::Error; mod closed; -mod closed_spherical; mod open; +mod open_spherical; mod periodic; pub use closed::Closed; -pub use closed_spherical::ClosedSpherical; pub use open::Open; +pub use open_spherical::OpenSpherical; pub use periodic::Periodic; /// Enumerate possible sources of error in fallible boundary methods. diff --git a/hoomd-microstate/src/boundary/closed_spherical.rs b/hoomd-microstate/src/boundary/open_spherical.rs similarity index 56% rename from hoomd-microstate/src/boundary/closed_spherical.rs rename to hoomd-microstate/src/boundary/open_spherical.rs index 65730b43d..335d01540 100644 --- a/hoomd-microstate/src/boundary/closed_spherical.rs +++ b/hoomd-microstate/src/boundary/open_spherical.rs @@ -1,7 +1,7 @@ // Copyright (c) 2024-2026 The Regents of the University of Michigan. // Part of hoomd-rs, released under the BSD 3-Clause License. -//! Implement `ClosedSpherical` +//! Implement `OpenSpherical` use arrayvec::ArrayVec; use serde::{Deserialize, Serialize}; @@ -10,35 +10,32 @@ use super::{Error, GenerateGhosts, MAX_GHOSTS, Wrap}; use crate::property::Position; use hoomd_manifold::Spherical; -/// [`ClosedSpherical`] implements a hypercubic box enclosing a unit-radius -/// $`(N-1)`$-sphere. Use `ClosedSpherical` alongside [`SphericalVecCell`] to -/// implement [`ParallelSweep`] for [`Spherical`] bodies. `ClosedSpherical` is +/// [`OpenSpherical`] implements a hypercubic box enclosing a unit-radius +/// $`(N-1)`$-sphere. Use `OpenSpherical` alongside [`SphericalVecCell`] to +/// implement [`ParallelSweep`] for [`Spherical`] bodies. `OpenSpherical` is /// otherwise functionally identical to using [`Open`] for [`Spherical`] /// simulations. /// /// # Example /// ``` -/// use hoomd_microstate::boundary::ClosedSpherical; +/// use hoomd_microstate::boundary::OpenSpherical; /// /// # fn main() -> Result<(), Box> { -/// let closed_spherical: ClosedSpherical<3> = ClosedSpherical {}; +/// let closed_spherical: OpenSpherical<3> = OpenSpherical {}; /// # Ok(()) /// # } /// ``` /// -/// Similar to [`Closed`], `ClosedSpherical` does not wrap bodies and sites, -/// nor does it generate ghost sites. However, unlike [`Closed`], -/// `ClosedSpherical` does not implement `IsPointInside`, `MapPoint`, -/// `Scale`, or `Volume`, as these methods are unrelated to the boundary -/// conditions for `Spherical` simulations. +/// Similar to [`Open`], `OpenSpherical` does not wrap bodies and sites, +/// nor does it generate ghost sites. /// /// [`SphericalVecCell`]: hoomd_spatial::SphericalVecCell; /// [`ParallelSweep`]: hoomd_mc::ParallelSweep; /// [`Spherical`]: hoomd_manifold::Spherical; #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] -pub struct ClosedSpherical {} +pub struct OpenSpherical {} -impl Wrap for ClosedSpherical +impl Wrap for OpenSpherical where BS: Position>, { @@ -48,7 +45,7 @@ where } } -impl GenerateGhosts for ClosedSpherical +impl GenerateGhosts for OpenSpherical where S: Default, { From fb93ab4207cb332c7980d037be0f91f450f2b895 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 26 Jun 2026 15:18:55 -0400 Subject: [PATCH 16/16] clippy --- hoomd-spatial/src/spherical_vec_cell.rs | 40 +++++++++---------------- 1 file changed, 14 insertions(+), 26 deletions(-) diff --git a/hoomd-spatial/src/spherical_vec_cell.rs b/hoomd-spatial/src/spherical_vec_cell.rs index 211d262fa..734d0458d 100644 --- a/hoomd-spatial/src/spherical_vec_cell.rs +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -1,15 +1,6 @@ // Copyright (c) 2024-2026 The Regents of the University of Michigan. // Part of hoomd-rs, released under the BSD 3-Clause License. -#![expect( - clippy::cast_possible_truncation, - reason = "the necessary conversions are necessary and have been checked" -)] -#![expect( - clippy::cast_sign_loss, - reason = "the necessary conversions are necessary and have been checked" -)] - //! Implement `SphericalVecCell` use serde::{Deserialize, Serialize}; @@ -222,10 +213,6 @@ where /// # Ok(()) /// # } /// ``` - #[expect( - clippy::missing_panics_doc, - reason = "hard-coded constant will never panic" - )] #[inline] #[must_use] pub fn builder() -> SphericalVecCellBuilder { @@ -268,7 +255,7 @@ where /// ``` #[inline] fn insert(&mut self, key: K, position: Spherical) { - VecCell::::insert(&mut self.0, key, *position.point()) + VecCell::::insert(&mut self.0, key, *position.point()); } /// Remove the point with the given key. @@ -365,7 +352,7 @@ where /// ``` #[inline] fn clear(&mut self) { - self.0.clear() + self.0.clear(); } } @@ -519,18 +506,19 @@ mod tests { .build(); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0).point() + Spherical::<3>::from_polar_coordinates(PI / 2.0, 0.0).point() ) == [2, 0, 0] ); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<3>::from_polar_coordinates(PI / 2.0, PI / 2.0).point() + Spherical::<3>::from_polar_coordinates(PI / 2.0, PI / 2.0).point() ) == [0, 2, 0] ); check!( - spherical_cell_list.0.cell_index_from_position( - &Spherical::<3>::from_polar_coordinates(0.0, 0.0).point() - ) == [0, 0, 2] + spherical_cell_list + .0 + .cell_index_from_position(Spherical::<3>::from_polar_coordinates(0.0, 0.0).point()) + == [0, 0, 2] ); } @@ -541,22 +529,22 @@ mod tests { .build(); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<4>::from_polar_coordinates(PI / 2.0, 0.0, 0.0).point() + Spherical::<4>::from_polar_coordinates(PI / 2.0, 0.0, 0.0).point() ) == [2, 0, 0, 0] ); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, 0.0).point() + Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, 0.0).point() ) == [0, 2, 0, 0] ); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, PI / 2.0).point() + Spherical::<4>::from_polar_coordinates(PI / 2.0, PI / 2.0, PI / 2.0).point() ) == [0, 0, 2, 0] ); check!( spherical_cell_list.0.cell_index_from_position( - &Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0).point() + Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0).point() ) == [0, 0, 0, 2] ); } @@ -719,7 +707,7 @@ mod tests { let key = key_distribution.sample(&mut rng); cell_list.insert(key, position); - reference.insert(key, cell_list.0.cell_index_from_position(&position.point())); + reference.insert(key, cell_list.0.cell_index_from_position(position.point())); } else { let key = key_distribution.sample(&mut rng); cell_list.remove(&key); @@ -772,7 +760,7 @@ mod tests { let key = key_distribution.sample(&mut rng); cell_list.insert(key, position); - reference.insert(key, cell_list.0.cell_index_from_position(&position.point())); + reference.insert(key, cell_list.0.cell_index_from_position(position.point())); } else { let key = key_distribution.sample(&mut rng); cell_list.remove(&key);