diff --git a/doc/src/release-notes.md b/doc/src/release-notes.md index d456d6d90..7886cf270 100644 --- a/doc/src/release-notes.md +++ b/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>` (#287). * `[hoomd-utility]`: Implement `Eq`, `PartialOrd`, and `Ord` for `PositiveReal` (#287). *Changed:* +* `[hoomd-mc]`: Improve the numerical stability of translation moves for `Point>` (#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:* diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 9145b7ced..522376e4c 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/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 Spherical { 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> { + /// 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 { /// Max distance away from point. pub disk_radius: PositiveReal, /// The center of the disk. - pub point: Spherical<3>, + pub point: Spherical, } impl Default for Spherical { @@ -222,7 +309,9 @@ impl Default for Spherical { } } -impl Distribution> for SphericalDisk { +impl Distribution> for SphericalDisk<3> { + /// Translates 3-dimensional cartesian vector named "point" along the + /// surface of a sphere by maximum distance of r. #[inline] fn sample(&self, rng: &mut R) -> Spherical<3> { let max_angle = self.disk_radius.get(); @@ -274,6 +363,34 @@ impl Distribution> for SphericalDisk { } } +impl Distribution> for SphericalDisk<4> { + /// Translates 3-dimensional cartesian vector named "point" along the + /// surface of a sphere by maximum distance of r. + #[inline] + fn sample(&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 + ); + } } diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index 9c44dfa7e..d6ff3b0f6 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -3,12 +3,13 @@ //! Implement Translation moves on curved surfaces -use rand::{Rng, distr::Distribution}; +use rand::{Rng, RngExt}; +use rand_distr::StandardNormal; use crate::{LocalTrial, Translate}; -use hoomd_manifold::{Spherical, SphericalDisk}; +use hoomd_manifold::Spherical; use hoomd_microstate::property::{Point, Position}; -use hoomd_vector::InnerProduct; +use hoomd_vector::{Cartesian, InnerProduct}; impl LocalTrial>> for Translate>> { /// Propose local trial moves for a body on the surface of a sphere @@ -52,16 +53,92 @@ impl LocalTrial>> for Translate>> { body_properties: Point>, ) -> Point> { let mut trial = body_properties; - let disk = SphericalDisk { - disk_radius: *self.maximum_distance(), - point: *trial.position_mut(), - }; - *trial.position_mut() = disk.sample(rng); - let rescale = 1.0 / trial.position().point().norm(); - *trial.position_mut() = - Spherical::from_cartesian_coordinates(*trial.position().point() * rescale); + let displacement = (self.maximum_distance().get()) * rng.sample::(StandardNormal); + let (sn, cs) = (displacement.sin(), displacement.cos()); + let vec = Cartesian::<3>::from(std::array::from_fn(|_| rng.sample(StandardNormal))); + let proj = vec.dot(trial.position.point()); + let tangent = Cartesian::from([ + vec[0] - proj * trial.position.coordinates()[0], + vec[1] - proj * trial.position.coordinates()[1], + vec[2] - proj * trial.position.coordinates()[2], + ]); + let (unit, _norm) = tangent.to_unit().expect("cannot be null"); + let new = Cartesian::from([ + trial.position.coordinates()[0] * cs + unit.get().coordinates[0] * sn, + trial.position.coordinates()[1] * cs + unit.get().coordinates[1] * sn, + trial.position.coordinates()[2] * cs + unit.get().coordinates[2] * sn, + ]); + *trial.position_mut() = Spherical::from_cartesian_coordinates(new); trial } } -// TODO: test +impl LocalTrial>> for Translate>> { + /// Propose local trial moves for a body on the surface of a 3-sphere + /// + /// # Example + /// ``` + /// use approxim::assert_relative_eq; + /// use hoomd_manifold::{Spherical, SphericalDisk}; + /// use hoomd_mc::{LocalTrial, Translate}; + /// use hoomd_microstate::property::{Point, Position}; + /// use hoomd_vector::{Cartesian, InnerProduct, Metric, Vector}; + /// use rand::{Rng, SeedableRng, rngs::StdRng}; + /// use std::f64::consts::PI; + /// + /// # fn main() -> Result<(), Box> { + /// let mut rng = StdRng::seed_from_u64(14); + /// let initial_point = Point::new(Spherical::<4>::from_polar_coordinates( + /// PI / 4.0, + /// PI / 10.0, + /// 5.0 * PI / 4.0, + /// )); + /// let d = 0.1; + /// let translate = Translate::with_maximum_distance(d.try_into()?); + /// + /// let new_body_properties = translate.propose(&mut rng, initial_point); + /// + /// // Translation move keeps point on the surface of the sphere + /// assert_relative_eq!( + /// new_body_properties.position().point().norm(), + /// 1.0_f64, + /// epsilon = 1e-8 + /// ); + /// + /// // Translation move does not translate the point more than a distance d away + /// assert!( + /// d > new_body_properties + /// .position() + /// .distance(&initial_point.position()) + /// ); + /// # Ok(()) + /// # } + /// ``` + #[inline] + fn propose( + &self, + rng: &mut R, + body_properties: Point>, + ) -> Point> { + let mut trial = body_properties; + let displacement = (self.maximum_distance().get()) * rng.sample::(StandardNormal); + let (sn, cs) = (displacement.sin(), displacement.cos()); + let vec = Cartesian::<4>::from(std::array::from_fn(|_| rng.sample(StandardNormal))); + let proj = vec.dot(trial.position.point()); + let tangent = Cartesian::from([ + vec[0] - proj * trial.position.coordinates()[0], + vec[1] - proj * trial.position.coordinates()[1], + vec[2] - proj * trial.position.coordinates()[2], + vec[3] - proj * trial.position.coordinates()[3], + ]); + let (unit, _norm) = tangent.to_unit().expect("cannot be null"); + let new = Cartesian::from([ + trial.position.coordinates()[0] * cs + unit.get().coordinates[0] * sn, + trial.position.coordinates()[1] * cs + unit.get().coordinates[1] * sn, + trial.position.coordinates()[2] * cs + unit.get().coordinates[2] * sn, + trial.position.coordinates()[3] * cs + unit.get().coordinates[3] * sn, + ]); + *trial.position_mut() = Spherical::from_cartesian_coordinates(new); + trial + } +} diff --git a/hoomd-microstate/src/property/point.rs b/hoomd-microstate/src/property/point.rs index df19713c4..133856f20 100644 --- a/hoomd-microstate/src/property/point.rs +++ b/hoomd-microstate/src/property/point.rs @@ -222,7 +222,7 @@ impl Transform>> for Point> { /// Move `Point>` properties from the local body frame to the /// system frame. /// - /// All positions on the 3-sphere are associated with some $`SO(4)`$ + /// All positions on the 3-sphere are associated with some $`SU(2)`$ /// transformation which translates the origin to that position. The local /// body frame is the frame in which the body position is the origin. The /// position of the sites in the system frame is obtained by applying the