diff --git a/Cargo.lock b/Cargo.lock index 5d7220e09..77f1cd115 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3894,7 +3894,9 @@ dependencies = [ name = "hoomd-spatial" version = "1.1.0" dependencies = [ + "approxim", "assert2", + "hoomd-manifold", "hoomd-utility", "hoomd-vector", "log", 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 } } diff --git a/hoomd-mc/src/hypercuboid.rs b/hoomd-mc/src/hypercuboid.rs index 370bdb1db..970146a96 100644 --- a/hoomd-mc/src/hypercuboid.rs +++ b/hoomd-mc/src/hypercuboid.rs @@ -19,7 +19,8 @@ use serde_with::serde_as; use std::array; use hoomd_geometry::shape::Hypercuboid; -use hoomd_microstate::boundary::{Closed, Periodic}; +use hoomd_manifold::Spherical; +use hoomd_microstate::boundary::{Closed, OpenSpherical, Periodic}; use hoomd_utility::valid::PositiveReal; use hoomd_vector::Cartesian; @@ -136,6 +137,26 @@ 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()), + ) + } + + #[inline] + fn space_indices_by_color(&self) -> &[Vec] { + >>::space_indices_by_color(self) + } + + #[inline] + fn num_spaces(&self) -> usize { + >>::num_spaces(self) + } +} + impl HypercuboidCheckerboard { /// Collapse a multi-dimensional index to a single value in `[0, num_spaces]`. #[inline] @@ -366,6 +387,38 @@ impl Cover> for Periodic> { } } +impl Cover> for OpenSpherical { + 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..451e5e29d 100644 --- a/hoomd-microstate/src/boundary.rs +++ b/hoomd-microstate/src/boundary.rs @@ -25,10 +25,12 @@ use thiserror::Error; mod closed; mod open; +mod open_spherical; mod periodic; pub use closed::Closed; 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/open_spherical.rs b/hoomd-microstate/src/boundary/open_spherical.rs new file mode 100644 index 000000000..335d01540 --- /dev/null +++ b/hoomd-microstate/src/boundary/open_spherical.rs @@ -0,0 +1,60 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! Implement `OpenSpherical` + +use arrayvec::ArrayVec; +use serde::{Deserialize, Serialize}; + +use super::{Error, GenerateGhosts, MAX_GHOSTS, Wrap}; +use crate::property::Position; +use hoomd_manifold::Spherical; + +/// [`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::OpenSpherical; +/// +/// # fn main() -> Result<(), Box> { +/// let closed_spherical: OpenSpherical<3> = OpenSpherical {}; +/// # Ok(()) +/// # } +/// ``` +/// +/// 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 OpenSpherical {} + +impl Wrap for OpenSpherical +where + BS: Position>, +{ + #[inline] + fn wrap(&self, properties: BS) -> Result { + Ok(properties) + } +} + +impl GenerateGhosts for OpenSpherical +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/Cargo.toml b/hoomd-spatial/Cargo.toml index d3d6bec1d..2a3c61eb7 100644 --- a/hoomd-spatial/Cargo.toml +++ b/hoomd-spatial/Cargo.toml @@ -23,10 +23,12 @@ rustc-hash.workspace = true serde.workspace = true serde_with.workspace = true +hoomd-manifold.workspace = true 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/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..734d0458d --- /dev/null +++ b/hoomd-spatial/src/spherical_vec_cell.rs @@ -0,0 +1,787 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! Implement `SphericalVecCell` + +use serde::{Deserialize, Serialize}; +use serde_with::serde_as; +use std::{array, cmp::Eq, f64::consts::PI, fmt, hash::Hash}; + +use hoomd_manifold::Spherical; +use hoomd_utility::valid::PositiveReal; + +use super::{PointUpdate, PointsNearBall, WithSearchRadius}; + +use crate::{ + IndexFromPosition, + vec_cell::{PointsIterator, VecCell, VecCellBuilder}, +}; + +/// 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(pub VecCell); + +/// Construct a [`SphericalVecCell`] with given parameters. +/// +/// # Example +/// +/// ``` +/// use hoomd_manifold::Spherical; +/// use hoomd_spatial::SphericalVecCell; +/// use std::f64::consts::PI; +/// +/// # fn main() -> 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(pub VecCellBuilder); + +impl SphericalVecCellBuilder +where + K: Copy + Eq + Hash, +{ + /// 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 spherical_nominal_search_radius( + mut self, + spherical_nominal_search_radius: PositiveReal, + ) -> Self { + 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"); + self + } + + /// 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 spherical_maximum_search_radius(mut self, spherical_maximum_search_radius: f64) -> Self { + self.0.maximum_search_radius = (spherical_maximum_search_radius.clamp(0.0, PI / 2.0)).sin(); + self + } + + /// 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.0.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.0.maximum_search_radius = euclidean_maximum_search_radius; + self + } + + /// Construct the [`SphericalVecCell`] with the chosen parameters. + #[inline] + #[must_use] + pub fn build(self) -> SphericalVecCell { + SphericalVecCell(VecCellBuilder::::build(self.0)) + } +} + +impl Default for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// 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() + } +} + +impl WithSearchRadius for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// 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(radius: PositiveReal) -> Self { + let euclidean_radius = (radius.get().clamp(0.0, PI / 2.0)).sin(); + Self::builder() + .euclidean_nominal_search_radius( + PositiveReal::try_from(euclidean_radius) + .expect("positive number given previous check"), + ) + .euclidean_maximum_search_radius(euclidean_radius) + .build() + } +} + +impl SphericalVecCell +where + K: Eq + Hash, +{ + /// Get the keys in a given cell index + #[cfg(test)] + #[inline] + fn get_keys(&self, cell_index: &[i64; D]) -> &[K] { + VecCell::::get_keys(&self.0, cell_index) + } +} + +impl SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Construct a `SphericalVecCell` builder. + /// + /// Use the builder to set any or all parameters and construct a [`SphericalVecCell`]. + /// + /// # Example + /// + /// ``` + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::SphericalVecCell; + /// + /// # fn main() -> Result<(), Box> { + /// let two_sphere_vec_cell = SphericalVecCell::::builder() + /// .euclidean_nominal_search_radius(0.2.try_into()?) + /// .euclidean_maximum_search_radius(0.4) + /// .build(); + /// # Ok(()) + /// # } + /// ``` + #[inline] + #[must_use] + pub fn builder() -> SphericalVecCellBuilder { + SphericalVecCellBuilder(VecCell::::builder()) + } + + /// 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) { + VecCell::::shrink_to_fit(&mut self.0); + } +} + +impl PointUpdate, K> for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Insert or update a point identified by a key. + /// + /// # 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) { + VecCell::::insert(&mut self.0, key, *position.point()); + } + + /// Remove the point with the given key. + /// + /// # 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 / 2.0, PI / 4.0, PI), + /// ); + /// + /// spherical_vec_cell.remove(&0) + /// ``` + #[inline] + fn remove(&mut self, key: &K) { + VecCell::::remove(&mut self.0, key); + } + + /// Get the number of points in the spatial data structure. + /// + /// # 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::<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!(spherical_vec_cell.len(), 2) + /// ``` + #[inline] + fn len(&self) -> usize { + self.0.cell_index.len() + } + + /// Test if the spatial data structure is empty. + /// + /// # Example + /// ``` + /// use hoomd_spatial::{PointUpdate, SphericalVecCell}; + /// + /// let mut spherical_vec_cell = SphericalVecCell::::default(); + /// + /// assert!(spherical_vec_cell.is_empty()); + /// ``` + #[inline] + fn is_empty(&self) -> bool { + self.0.cell_index.is_empty() + } + + /// Test if the spatial data structure contains a key. + /// ``` + /// 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::<3>::from_polar_coordinates(3.0 * PI / 5.0, 0.0), + /// ); + /// + /// assert!(spherical_vec_cell.contains_key(&0)); + /// ``` + #[inline] + fn contains_key(&self, key: &K) -> bool { + self.0.cell_index.contains_key(key) + } + + /// Remove all points. + /// + /// # 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, 0.0), + /// ); + /// + /// spherical_vec_cell.clear(); + /// ``` + #[inline] + fn clear(&mut self) { + self.0.clear(); + } +} + +impl Iterator for PointsIterator<'_, K, D, SphericalVecCell> +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 = + 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]); + } + } +} + +impl PointsNearBall, K> for SphericalVecCell +where + K: Copy + Eq + Hash, +{ + /// Find all the points that *might* be in the given ball with a specified + /// radius in the `Spherical` distance metric. + /// + /// `points_near_ball` will iterate over all points in the given ball *and + /// possibly others as well*. [`SphericalVecCell`] may iterate over the points in + /// any order. + /// + /// # Example + /// ``` + /// use hoomd_manifold::Spherical; + /// use hoomd_spatial::{PointUpdate, PointsNearBall, SphericalVecCell}; + /// use std::f64::consts::PI; + /// + /// let mut spherical_vec_cell = 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 spherical_vec_cell.points_near_ball( + /// &Spherical::<3>::from_cartesian_coordinates([0.0, 0.0, 1.0].into()), + /// PI / 4.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 { + // convert spherical distance to Euclidean distance + let euclidean_radius = radius.asin(); + VecCell::::points_near_ball(&self.0, position.point(), euclidean_radius) + } +} + +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::SphericalVecCell; + /// use log::info; + /// + /// let spherical_vec_cell = SphericalVecCell::::default(); + /// + /// info!("{spherical_vec_cell}"); + /// ``` + #[allow( + clippy::missing_inline_in_public_items, + reason = "no need to inline display" + )] + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + VecCell::::fmt(&self.0, f) + } +} + +impl IndexFromPosition> for SphericalVecCell +where + K: Eq + Hash, +{ + type Location = [i64; D]; + + #[inline] + fn location_from_position(&self, position: &Spherical) -> Self::Location { + 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, + distr::{Distribution, Uniform}, + rngs::StdRng, + }; + use rstest::*; + + use approxim::assert_relative_eq; + use assert2::{assert, check}; + use rustc_hash::FxHashMap; + 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.0.cell_index_from_position( + 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() + ) == [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] + ); + } + + #[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.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.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.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.0.cell_index_from_position( + Spherical::<4>::from_polar_coordinates(0.0, 0.0, 0.0).point() + ) == [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.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.0.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.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.0.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.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); + 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.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); + 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.0.cell_index_from_position(position.point())); + } 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.0.cell_index.len() == reference.len()); + for (reference_key, reference_value) in reference.drain() { + let value = cell_list.0.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.0.keys_map.iter().map(Vec::len).sum(); + check!(cell_list.0.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.0.cell_index_from_position(position.point())); + } 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.0.cell_index.len() == reference.len()); + for (reference_key, reference_value) in reference.drain() { + let value = cell_list.0.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.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 72f4a8ea7..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>")] @@ -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::>(); @@ -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] @@ -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, /// 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, {