Curved gsd

hoomd-bevy/src/representation/hyperbolic_polygon.wgsl

@@ -19,7 +19,11 @@ const PI: f32 = 3.141592653589793238462643;
 
 @group(2) @binding(4) var image_color_texture: texture_2d<f32>;
 @group(2) @binding(5) var image_color_sampler: sampler;
-@group(2) @binding(6) var<storage, read> n_sides: f32;
+struct PolygonParams {
+    n_sides: f32,
+};
+@group(2) @binding(6)
+var<uniform> polygon: PolygonParams;
 
 struct VertexOutput {
     // this is `clip position` when the struct is used as a vertex stage output
@@ -75,7 +79,7 @@ fn vertex(vertex: Vertex) -> VertexOutput {
         vec4<f32>(vertex.position, 1.0)
     );
     out.position = mesh_functions::mesh2d_position_world_to_clip(out.world_position);
-    out.n_sides = n_sides;
+    out.n_sides = polygon.n_sides;
 #endif
 
 #ifdef VERTEX_NORMALS

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 2-sphere to an 2-dimensional
+    /// plane by projecting through the $`(0,0,1)`$ axis.
     ///
     /// # Example
     /// ```
@@ -100,7 +100,7 @@ impl Spherical<3> {
 }
 
 impl Spherical<4> {
-    /// Create a 3-sphere from spherical coordinates
+    /// Create a point on a unit-radius 3-sphere from a unit quaternion.
     #[inline]
     #[must_use]
     pub fn from_polar_coordinates(theta: f64, phi_1: f64, phi_2: f64) -> Spherical<4> {
@@ -115,6 +115,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 (a, b, c, d) = versor.get_components();
+        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> {
@@ -158,7 +232,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 +281,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 +297,7 @@ impl<const N: usize> Default for Spherical<N> {
     }
 }
 
-impl Distribution<Spherical<3>> for SphericalDisk {
+impl Distribution<Spherical<3>> for SphericalDisk<3> {
     #[inline]
     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Spherical<3> {
         let max_angle = self.disk_radius.get();
@@ -274,6 +349,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("sperical 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 +432,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 +453,26 @@ 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);
+        }
+    }
 }

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<Point<Spherical<3>>> for Translate<Point<Spherical<3>>> {
     /// Propose local trial moves for a body on the surface of a sphere
@@ -25,7 +26,7 @@ impl LocalTrial<Point<Spherical<3>>> for Translate<Point<Spherical<3>>> {
     /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
     /// let mut rng = StdRng::seed_from_u64(14);
     /// let initial_point = Point::new(Spherical::from_cartesian_coordinates(
-    ///     [0.5_f64.sqrt(), 0.5_f64.sqrt(), 0.0].into(),
+    ///     [0.5_f64.sqrt(), 0.0, 0.5_f64.sqrt()].into(),
     /// ));
     /// let d = 0.1;
     /// let translate = Translate::with_maximum_distance(d.try_into()?);
@@ -52,16 +53,104 @@ impl LocalTrial<Point<Spherical<3>>> for Translate<Point<Spherical<3>>> {
         body_properties: Point<Spherical<3>>,
     ) -> Point<Spherical<3>> {
         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::<f64, _>(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
+        //         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);
+        // trial
     }
 }
 
-// TODO: test
+impl LocalTrial<Point<Spherical<4>>> for Translate<Point<Spherical<4>>> {
+    /// 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<dyn std::error::Error>> {
+    /// 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<R: Rng>(
+        &self,
+        rng: &mut R,
+        body_properties: Point<Spherical<4>>,
+    ) -> Point<Spherical<4>> {
+        let mut trial = body_properties;
+        let displacement = (self.maximum_distance().get()) * rng.sample::<f64, _>(StandardNormal);
+        let (sn, cs) = ((displacement.abs()).sin(), (displacement.abs()).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
+    }
+}
+
+// TODO: test code