Merge pull request #285 from glotzerlab/spherical

Spherical

Author: joaander

doc/src/release-notes.md

@@ -4,16 +4,22 @@
 
 *Added:*
 
+* `[hoomd-manifold]`: Add `Spherical<4>::from_versor` and the corresponding `::to_versor` (#285).
+* `[hoomd-mc]`: Implement translation moves for `Point<Spherical<4>>` (#287).
 * `[hoomd-utility]`: Implement `Eq`, `PartialOrd`, and `Ord` for `PositiveReal` (#287).
 
 *Changed:*
 
+* `[hoomd-mc]`: Improve the numerical stability of translation moves for `Point<Spherical<3>>` (#287).
+
 *Deprecated:*
 
 *Removed:*
 
 *Fixed:*
 
+* `[hoomd-manifold]`: Fixed numerical stability issue in `Spherical<3>::distance` where the dot product could result in an out of bounds value.
+
 ## 1.1.0 (2026-04-17)
 
 *Highlights:*

hoomd-manifold/src/sphere.rs

@@ -5,14 +5,14 @@
 
 use approxim::{approx_derive::RelativeEq, assert_relative_eq};
 use rand::{
-    Rng,
+    Rng, RngExt,
     distr::{Distribution, Uniform},
 };
 use serde::{Deserialize, Serialize};
 use std::f64::consts::PI;
 
 use hoomd_utility::valid::PositiveReal;
-use hoomd_vector::{Cartesian, InnerProduct, Metric};
+use hoomd_vector::{Cartesian, InnerProduct, Metric, Quaternion, Rotate, Versor};
 
 /// Point on the surface of a sphere.
 ///
@@ -53,8 +53,8 @@ impl<const N: usize> Spherical<N> {
         assert_relative_eq!(rad, 1.0_f64, epsilon = 1e-6);
         Spherical { point }
     }
-
-    /// Implements a stereographic projection from the N-sphere to an N-dimensional plane.
+    /// Implements a stereographic projection from the N-sphere to an n-dimensional
+    /// plane by projecting through the $`(0,\cdots, 0,1)`$ axis.
     ///
     /// # Example
     /// ```
@@ -100,7 +100,18 @@ impl Spherical<3> {
 }
 
 impl Spherical<4> {
-    /// Create a 3-sphere from spherical coordinates
+    /// Create a point on a 3-sphere from spherical coordinates. Note that this uses
+    /// the convention
+    /// ```math
+    /// \begin{pmatrix}
+    /// \sin\theta \cos\phi_1
+    /// \\ \sin\theta \sin\phi_1 \cos\phi_2
+    /// \\ \sin\theta \sin\phi_1 \sin\phi_2
+    /// \\ \cos\theta
+    /// \end{pmatrix}
+    /// ```
+    /// where $`\theta`$ and $`phi_1`$ both run over the range $`0`$ to $`\pi`$ and $`\phi_2`$
+    /// runs from $`0`$ to $`2\pi`$.
     #[inline]
     #[must_use]
     pub fn from_polar_coordinates(theta: f64, phi_1: f64, phi_2: f64) -> Spherical<4> {
@@ -115,6 +126,80 @@ impl Spherical<4> {
         ]);
         Spherical::from_cartesian_coordinates(point)
     }
+    /// Create a point on a unit-radius 3-sphere from a unit quaternion.
+    #[inline]
+    #[must_use]
+    pub fn from_versor(versor: Versor) -> Spherical<4> {
+        let quat = versor.get();
+        let (a, b, c, d) = (quat.scalar, quat.vector[0], quat.vector[1], quat.vector[2]);
+        Spherical::<4>::from_cartesian_coordinates(Cartesian::from([a, b, c, d]))
+    }
+    /// Create a versor which maps $`(1,0,0,0)`$ to the target `Spherical<4>` point.
+    /// # Example
+    /// ```
+    /// use approxim::assert_relative_eq;
+    /// use hoomd_manifold::Spherical;
+    /// use hoomd_vector::{Cartesian, Quaternion, Versor};
+    /// use std::f64::consts::PI;
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// let radius = 1.0;
+    /// let x = Spherical::<4>::from_polar_coordinates(
+    ///     PI / 4.0,
+    ///     PI / 8.0,
+    ///     5.0 * PI / 4.0,
+    /// );
+    /// let x_versor = x.to_versor();
+    /// let pole_versor = Quaternion::from([1.0, 0.0, 0.0, 0.0])
+    ///     .to_versor()
+    ///     .expect("not a null vector");
+    /// let transformation =
+    ///     (*x_versor.get() * *pole_versor.get() * *x_versor.get())
+    ///         .to_versor()
+    ///         .expect("Hard-coded example is valid");
+    /// let mapped_pole = Spherical::<4>::from_versor(transformation);
+    /// assert_relative_eq!(
+    ///     mapped_pole.coordinates()[0],
+    ///     x.coordinates()[0],
+    ///     epsilon = 1e-12
+    /// );
+    /// assert_relative_eq!(
+    ///     mapped_pole.coordinates()[1],
+    ///     x.coordinates()[1],
+    ///     epsilon = 1e-12
+    /// );
+    /// assert_relative_eq!(
+    ///     mapped_pole.coordinates()[2],
+    ///     x.coordinates()[2],
+    ///     epsilon = 1e-12
+    /// );
+    /// assert_relative_eq!(
+    ///     mapped_pole.coordinates()[3],
+    ///     x.coordinates()[3],
+    ///     epsilon = 1e-12
+    /// );
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    #[must_use]
+    pub fn to_versor(&self) -> Versor {
+        let phi = self.coordinates()[3].atan2(self.coordinates()[2]);
+        let theta = ((self.coordinates()[3].powi(2) + self.coordinates()[2].powi(2)).sqrt())
+            .atan2(self.coordinates()[1]);
+        let psi = ((self.coordinates()[3].powi(2)
+            + self.coordinates()[2].powi(2)
+            + self.coordinates()[1].powi(2))
+        .sqrt())
+        .atan2(self.coordinates()[0]);
+        let n_hat = Cartesian::from([
+            theta.cos(),
+            (theta.sin()) * (phi.cos()),
+            (theta.sin()) * (phi.sin()),
+        ])
+        .to_unit_unchecked();
+        Versor::from_axis_angle(n_hat.0, psi)
+    }
 }
 
 impl Metric for Spherical<3> {
@@ -131,7 +216,8 @@ impl Metric for Spherical<3> {
     #[inline]
     fn distance(&self, other: &Self) -> f64 {
         let arg = Cartesian::dot(&self.point, &other.point);
-        arg.acos()
+        let arg_clipped = arg.clamp(-1.0, 1.0);
+        arg_clipped.acos()
     }
     #[inline]
     fn distance_squared(&self, other: &Self) -> f64 {
@@ -158,7 +244,8 @@ impl Metric for Spherical<4> {
     #[inline]
     fn distance(&self, other: &Self) -> f64 {
         let arg = Cartesian::dot(&self.point, &other.point);
-        arg.acos()
+        let arg_clipped = arg.clamp(-1.0, 1.0);
+        arg_clipped.acos()
     }
     #[inline]
     fn distance_squared(&self, other: &Self) -> f64 {
@@ -206,11 +293,11 @@ impl Metric for Spherical<4> {
 /// # }
 /// ```
 #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
-pub struct SphericalDisk {
+pub struct SphericalDisk<const N: usize> {
     /// Max distance away from point.
     pub disk_radius: PositiveReal,
     /// The center of the disk.
-    pub point: Spherical<3>,
+    pub point: Spherical<N>,
 }
 
 impl<const N: usize> Default for Spherical<N> {
@@ -222,7 +309,9 @@ impl<const N: usize> Default for Spherical<N> {
     }
 }
 
-impl Distribution<Spherical<3>> for SphericalDisk {
+impl Distribution<Spherical<3>> for SphericalDisk<3> {
+    /// Translates 3-dimensional cartesian vector named "point" along the
+    /// surface of a sphere by maximum distance of r.
     #[inline]
     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Spherical<3> {
         let max_angle = self.disk_radius.get();
@@ -274,6 +363,34 @@ impl Distribution<Spherical<3>> for SphericalDisk {
     }
 }
 
+impl Distribution<Spherical<4>> for SphericalDisk<4> {
+    /// Translates 3-dimensional cartesian vector named "point" along the
+    /// surface of a sphere by maximum distance of r.
+    #[inline]
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Spherical<4> {
+        let max_trans = self.disk_radius.get();
+        let point = self.point;
+        // generate random unit cartesian vector
+        let v: Versor = rng.random();
+        let b_hat = v
+            .rotate(&Cartesian::from([1.0, 0.0, 0.0]))
+            .to_unit()
+            .expect("hard coded non-null vector");
+        let eta = Uniform::new(0.0, max_trans).expect("hard coded non-negative");
+        let translation_versor = Versor::from_axis_angle(b_hat.0, eta.sample(rng));
+
+        let position_versor = Quaternion::from(*point.coordinates())
+            .to_versor()
+            .expect("spherical points cannot be null");
+        let transformation =
+            ((*translation_versor.get()) * (*position_versor.get()) * (*translation_versor.get()))
+                .to_versor()
+                .expect("spherical points cannot be null");
+        let sphere_point = Spherical::<4>::from_versor(transformation);
+        Spherical::<4>::from_cartesian_coordinates(*sphere_point.point())
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -329,7 +446,7 @@ mod tests {
     }
 
     #[test]
-    fn random_sphere() {
+    fn random_two_sphere() {
         // Generate ten random points on the Hyperbolic
         let mut rng = StdRng::seed_from_u64(42);
         let d = 0.1;
@@ -350,4 +467,70 @@ mod tests {
             assert!(d > distance);
         }
     }
+    #[test]
+    fn random_three_sphere() {
+        // Generate ten random points on the Hyperbolic
+        let mut rng = StdRng::seed_from_u64(42);
+        let d = 0.1;
+        let n_pole = Spherical::from_cartesian_coordinates(Cartesian::from([1.0, 0.0, 0.0, 0.0]));
+        for _n in 0..10 {
+            let disk = SphericalDisk {
+                disk_radius: d.try_into().expect("hard-coded positive number"),
+                point: n_pole,
+            };
+            let random_point: Spherical<4> = disk.sample(&mut rng);
+
+            // check that points remain on Sphere
+            let rho = random_point.point.norm_squared();
+            assert_relative_eq!(rho, 1.0, epsilon = 1e-12);
+
+            // check that points are within distance d of north pole
+            let distance = random_point.point().distance(n_pole.point());
+            assert!(d > distance);
+        }
+    }
+    #[test]
+    fn from_versor() {
+        // generate a 3-sphere point from a versor
+        let mut rng = StdRng::seed_from_u64(358);
+        let v: Versor = rng.random();
+        let sphere_pt = Spherical::<4>::from_versor(v);
+        assert_eq!(sphere_pt.coordinates()[0], v.get().scalar);
+        assert_eq!(sphere_pt.coordinates()[1], v.get().vector.coordinates[0]);
+        assert_eq!(sphere_pt.coordinates()[2], v.get().vector.coordinates[1]);
+        assert_eq!(sphere_pt.coordinates()[3], v.get().vector.coordinates[2]);
+    }
+    #[test]
+    fn to_versor() {
+        // generate a 3-sphere point from a versor
+        let mut rng = StdRng::seed_from_u64(5121);
+        let v: Versor = rng.random();
+        let sphere_pt = Spherical::<4>::from_versor(v);
+        // get versor which maps identity versor to 3-sphere point
+        let versor_from_spherical = sphere_pt.to_versor();
+        let pole_versor = Quaternion::from([1.0, 0.0, 0.0, 0.0])
+            .to_versor()
+            .expect("not a null vector");
+        let transformation =
+            (*versor_from_spherical.get() * *pole_versor.get() * *versor_from_spherical.get())
+                .to_versor()
+                .expect("Hard-coded example is valid");
+        let q = transformation.get();
+        assert_relative_eq!(q.scalar, v.get().scalar, epsilon = 1e-12);
+        assert_relative_eq!(
+            q.vector.coordinates[0],
+            v.get().vector.coordinates[0],
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            q.vector.coordinates[1],
+            v.get().vector.coordinates[1],
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            q.vector.coordinates[2],
+            v.get().vector.coordinates[2],
+            epsilon = 1e-12
+        );
+    }
 }