From 9935ce72ea7b93803b2a3d4f574beb96eb94b10b Mon Sep 17 00:00:00 2001 From: mthran Date: Wed, 17 Dec 2025 16:08:35 -0500 Subject: [PATCH 01/13] translate for 3-sphere --- hoomd-interaction/src/univariate.rs | 2 + .../src/univariate/zetterling.rs | 102 +++++++++++++++ hoomd-manifold/src/sphere.rs | 120 ++++++++++++++---- hoomd-mc/src/translate/sphere.rs | 11 +- hoomd-microstate/src/property/point.rs | 43 ++++--- hoomd-vector/src/quaternion.rs | 9 ++ 6 files changed, 236 insertions(+), 51 deletions(-) create mode 100644 hoomd-interaction/src/univariate/zetterling.rs diff --git a/hoomd-interaction/src/univariate.rs b/hoomd-interaction/src/univariate.rs index bcc88870e..ea6ee6d75 100644 --- a/hoomd-interaction/src/univariate.rs +++ b/hoomd-interaction/src/univariate.rs @@ -13,6 +13,7 @@ mod overlap_penalty; mod shifted; mod weeks_chandler_anderson; mod xplor; +mod zetterling; pub use boxcar::Boxcar; pub use expanded::Expanded; @@ -24,6 +25,7 @@ pub use overlap_penalty::OverlapPenalty; pub use shifted::Shifted; pub use weeks_chandler_anderson::WeeksChandlerAnderson; pub use xplor::Xplor; +pub use zetterling::Zetterling; /// Computes energy as a function of one variable. /// diff --git a/hoomd-interaction/src/univariate/zetterling.rs b/hoomd-interaction/src/univariate/zetterling.rs new file mode 100644 index 000000000..7c70f39e2 --- /dev/null +++ b/hoomd-interaction/src/univariate/zetterling.rs @@ -0,0 +1,102 @@ +// Copyright (c) 2024-2025 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +//! Implement [`Zetterling`] + +use super::{UnivariateEnergy, UnivariateForce}; + +/// `Zetterling` computes the oscillating pair potential between every pair of +/// particles in the simulation state. +/// ```math +/// U(r) = \epsilon \frac{\exp(\alpha r/\ell)\cos(2k_Fr/\ell)}{(r/\ell)^3} + \beta\left(\frac{\sigma \ell}{r}\right)^n +/// ``` +#[derive(Clone, Debug, PartialEq)] +pub struct Zetterling { + /// Energy scale of the first term *(\[energy\])*. + pub epsilon: f64, + /// Screening factor *(\[unitless\])*. + pub alpha: f64, + /// Wave number to mimic the Friedel oscillations effect *(\[unitless\])*. + pub kf: f64, + /// Energy scale of the second term *(\[energy\])*. + pub beta: f64, + /// Repulsive core size *(\[unitless\])*. + pub sigma: f64, + /// The power to take sigma/r in the second term *(\[unitless\])*. + pub n: f64, + /// The length scale of the distances *(\[length\])* + pub ell: f64, +} + +impl UnivariateEnergy for Zetterling { + #[inline] + fn energy(&self, r: f64) -> f64 { + let sigma_ell_r = self.sigma * self.ell / r; + let r_ell = r / self.ell; + + self.epsilon * ((self.alpha * r_ell).exp()) * ((2.0 * self.kf * r_ell).cos()) / (r_ell.powi(3)) + self.beta * (sigma_ell_r.powf(self.n)) + } +} + +impl UnivariateForce for Zetterling { + #[inline] + fn force(&self, r: f64) -> f64 { + let r_ell = r / self.ell; + let r_inv = r.recip(); + let cos = (2.0 * self.kf * r_ell).cos(); + let sin = (2.0 * self.kf * r_ell).sin(); + let exp = (self.alpha * r_ell).exp(); + + let first = -self.epsilon * self.alpha * (self.ell.powi(2)) * exp * cos * (r_inv.powi(3)); + let second = self.epsilon * 2.0 * self.kf * (self.ell.powi(2)) * exp * sin * (r_inv.powi(3)); + let third = 3.0 * self.epsilon * (self.ell.powi(3)) * exp * cos * (r_inv.powi(4)); + let fourth = self.beta * self.n * ((self.sigma * self.ell).powf(self.n)) * (r_inv.powf(self.n+1.0)); + + first + second + third + fourth + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approxim::assert_relative_eq; + + #[test] + fn select_points() { + let zetterling1 = Zetterling { + epsilon: 1.58, + alpha: -0.22, + kf: 4.12, + beta: 0.95533, + sigma: 1.0, + n: 18.0, + ell: 1.0, + }; + + let (u_1, r_1) = (-0.742_644_124_392_870, 1.130_547_632_166_212); + assert_relative_eq!(u_1, zetterling1.energy(r_1), epsilon = 1e-12); + assert_relative_eq!(0.0, zetterling1.force(r_1), epsilon = 1e-6); + + let (u_2, r_2) = (-0.153_547_971_965_573, 1.879_994_132_512_336); + assert_relative_eq!(u_2, zetterling1.energy(r_2), epsilon = 1e-12); + assert_relative_eq!(0.0, zetterling1.force(r_2), epsilon = 1e-10); + + let zetterling2 = Zetterling { + epsilon: 1.04, + alpha: 0.33, + kf: 4.139, + beta: 0.94656, + sigma: 1.0, + n: 14.5, + ell: 1.0, + }; + + let (u_1, r_1) = (-0.883_354_387_732_971, 1.132_662_993_677_647); + assert_relative_eq!(u_1, zetterling2.energy(r_1), epsilon = 1e-12); + assert_relative_eq!(0.0, zetterling2.force(r_1), epsilon = 1e-6); + + let (u_2, r_2) = (-0.287_906_300_605_616, 1.879_255_603_675_727); + assert_relative_eq!(u_2, zetterling2.energy(r_2), epsilon = 1e-12); + assert_relative_eq!(0.0, zetterling2.force(r_2), epsilon = 1e-6); + } +} \ No newline at end of file diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index ea43fc3d0..78a40e4b5 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -12,7 +12,7 @@ use rand::{ }; use hoomd_utility::valid::PositiveReal; -use hoomd_vector::{Cartesian, InnerProduct, Metric}; +use hoomd_vector::{Cartesian, InnerProduct, Metric, Quaternion, Versor}; /// Point on the surface of a sphere. /// @@ -61,8 +61,24 @@ impl Spherical { assert_relative_eq!(rad, radius); Spherical { point, radius } } +} - /// Implements a stereographic projection from the N-sphere to an N-dimensional plane. +impl Spherical<3> { + /// Create a 2-sphere from spherical coordinates + #[inline] + #[must_use] + pub fn from_polar_coordinates(r: f64, theta: f64, phi: f64) -> Spherical<3> { + let theta_mod = theta.rem_euclid(PI); + let phi_mod = phi.rem_euclid(2.0 * PI); + let point = Cartesian::from([ + r * (theta_mod.sin()) * (phi_mod.cos()), + r * (theta_mod.sin()) * (phi_mod.sin()), + r * (theta_mod.cos()), + ]); + Spherical::from_cartesian_coordinates(point, r) + } + /// Implements a stereographic projection from the 2-sphere to an 2-dimensional + /// plane by projecting through the $`(0,0,1)`$ axis. /// /// # Example /// ``` @@ -84,47 +100,99 @@ impl Spherical { #[inline] #[must_use] pub fn stereographic_projection(&self) -> Vec { - (0..N - 1) + (0..2) .collect::>() .iter() .map(|i| { - self.point.coordinates[*i] / (1.0 - self.point.coordinates[N - 1] / self.radius) + self.point.coordinates[*i] / (1.0 - self.point.coordinates[2] / self.radius) }) .collect::>() } } -impl Spherical<3> { - /// Create a 2-sphere from spherical coordinates +impl Spherical<4> { + /// Create a point on a 3-sphere from spherical coordinates. Note that this uses + /// the convention + /// ```math + /// \begin{pmatrix}r\cos\psi + /// \\ r\sin\psi\cos\theta + /// \\ r\sin\psi\sin\theta\cos\phi + /// \\ r\sin\psi\sin\theta\sin\phi + /// \end{pmatrix} + /// ``` + /// where $`\psi`$ and $`theta`$ both run over the range $`0`$ to $`\pi`$ and $`\phi`$ + /// runs from $`0`$ to $`2\pi`$. #[inline] #[must_use] - pub fn from_polar_coordinates(r: f64, theta: f64, phi: f64) -> Spherical<3> { + pub fn from_polar_coordinates(r: f64, psi: f64, theta: f64, phi: f64) -> Spherical<4> { + let psi_mod = psi.rem_euclid(PI); let theta_mod = theta.rem_euclid(PI); let phi_mod = phi.rem_euclid(2.0 * PI); let point = Cartesian::from([ - r * (theta_mod.sin()) * (phi_mod.cos()), - r * (theta_mod.sin()) * (phi_mod.sin()), - r * (theta_mod.cos()), + r * (psi_mod.cos()), + r * (psi_mod.sin()) * (theta_mod.cos()), + r * (psi_mod.sin()) * (theta_mod.sin()) * (phi_mod.cos()), + r * (psi_mod.sin()) * (theta_mod.sin()) * (phi_mod.sin()), ]); Spherical::from_cartesian_coordinates(point, r) } -} - -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(r: f64, theta: f64, phi_1: f64, phi_2: f64) -> Spherical<4> { - let theta_mod = theta.rem_euclid(PI); - let phi_1_mod = phi_1.rem_euclid(PI); - let phi_2_mod = phi_2.rem_euclid(2.0 * PI); - let point = Cartesian::from([ - r * (theta_mod.sin()) * (phi_1_mod.cos()), - r * (theta_mod.sin()) * (phi_1_mod.sin()) * (phi_2_mod.cos()), - r * (theta_mod.sin()) * (phi_1_mod.sin()) * (phi_2_mod.sin()), - r * (theta_mod.cos()), - ]); - Spherical::from_cartesian_coordinates(point, r) + 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]), + 1.0) + } + /// Create a versor which maps $`(1,0,0,0)`$ to the target `Spherical<4>` point. + /// # Example + /// ``` + /// use hoomd_manifold::Spherical; + /// use hoomd_vector::{Cartesian, Quaternion, Versor}; + /// use approxim::assert_relative_eq; + /// use std::f64::consts::PI; + /// + /// # fn main() -> Result<(), Box> { + /// let radius = 1.0; + /// let x = Spherical::<4>::from_polar_coordinates(1.0,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]); + println!("phi {} theta {} psi {}", phi, theta, psi); + 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) + } + /// Implements a stereographic projection from the 2-sphere to an 2-dimensional + /// plane by projecting through the $`(1,0,0,0)`$ axis. + #[inline] + #[must_use] + pub fn stereographic_projection(&self) -> Vec { + (1..4) + .collect::>() + .iter() + .map(|i| { + self.point.coordinates[*i] / (1.0 - self.point.coordinates[0] / self.radius) + }) + .collect::>() } } @@ -240,7 +308,7 @@ impl Default for Spherical { #[inline] fn default() -> Self { let mut zero = Cartesian::::default(); - zero.coordinates[N - 1] = 1.0; + zero.coordinates[0] = 1.0; Spherical { point: zero, radius: 1.0_f64, diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index 123542a6d..d38b053b4 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -7,13 +7,10 @@ use rand::{Rng, distr::Distribution}; use crate::{LocalTrial, Translate}; use hoomd_manifold::{Spherical, SphericalDisk}; -use hoomd_microstate::property::Position; -use hoomd_vector::InnerProduct; +use hoomd_microstate::property::{Point, Position}; +use hoomd_vector::{InnerProduct, Quaternion, Versor}; -impl LocalTrial for Translate> -where - B: Position>, -{ +impl LocalTrial>> for Translate>> { /// Propose local trial moves for a body on the surface of a sphere /// /// # Example @@ -50,7 +47,7 @@ where /// # } /// ``` #[inline] - fn propose(&self, rng: &mut R, body_properties: B) -> B { + fn propose(&self, rng: &mut R, body_properties: Point>) -> Point> { let mut trial = body_properties; let disk = SphericalDisk { disk_radius: *self.maximum_distance(), diff --git a/hoomd-microstate/src/property/point.rs b/hoomd-microstate/src/property/point.rs index 913046011..9572fe094 100644 --- a/hoomd-microstate/src/property/point.rs +++ b/hoomd-microstate/src/property/point.rs @@ -6,7 +6,8 @@ use super::Position; use crate::Transform; use hoomd_manifold::{Hyperbolic, Minkowski, Spherical}; -use hoomd_vector::Cartesian; +use hoomd_vector::{Cartesian, Quaternion, Versor}; +use std::ops::Mul; /// A position in space and nothing more. /// @@ -153,6 +154,7 @@ impl Transform>> for Point> { /// local body frame. fn transform(&self, site_properties: &Point>) -> Point> { let radius = self.position.radius(); + // NOTE: this does not correctly conjugate the site let body_point = self.position.coordinates(); let body_phi = body_point[1].atan2(body_point[0]); let body_theta = (body_point[2] / radius).acos(); @@ -175,7 +177,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 @@ -184,6 +186,17 @@ impl Transform>> for Point> { #[inline] fn transform(&self, site_properties: &Point>) -> Point> { let radius = self.position.radius(); + let body_versor = self.position.to_versor(); + let site_versor = Quaternion::from(*site_properties.position.coordinates()).to_versor().expect("not a null vector"); + let transformation = ((*body_versor.get()) * (*site_versor.get()) * (*body_versor.get())).to_versor().expect("spherical points cannot be null"); + let sphere_point = Spherical::<4>::from_versor(transformation); + let rescaled = Spherical::<4>::from_cartesian_coordinates( + sphere_point.point().mul(radius), radius); + Point { + position: rescaled, + } + +/* let radius = self.position.radius(); let body_point = self.position.coordinates(); let body_phi_1 = (body_point[2].powi(2) + body_point[1].powi(2)) .sqrt() @@ -208,7 +221,7 @@ impl Transform>> for Point> { -trial_coords[0] * (body_theta.sin()) + trial_coords[3] * (body_theta.cos()), ]); let new_sphere = Spherical::from_cartesian_coordinates(transformed_point, radius); - Point::new(new_sphere) + Point::new(new_sphere) */ } } @@ -309,25 +322,19 @@ mod tests { #[test] fn transform_s3_point() { - let theta = PI / 5.0; - let blip = PI / 10.0; - let body = Point::new(Spherical::<4>::from_polar_coordinates(1.0, theta, 0.0, 0.0)); - let site = Point::new(Spherical::<4>::from_polar_coordinates( - 1.0, - blip, - PI / 2.0, - 0.0, - )); + let psi = 3.0*PI/5.0; + let blip = 0.13; + let body = Point::new(Spherical::<4>::from_polar_coordinates(1.0, psi, PI/4.0, 0.0)); + let site = Point::new(Spherical::<4>::from_polar_coordinates(1.0, blip, PI/4.0, 0.0)); let transformed_site = body.transform(&site); assert_relative_eq!( *transformed_site.position().point(), [ - (theta.sin()) * (blip.cos()), - blip.sin(), - 0.0, - (theta.cos()) * (blip.cos()) - ] - .into() + (psi+blip).cos(), + ((psi+blip).sin())*((PI/4.0).cos()), + ((psi+blip).sin())*((PI/4.0).sin()), + 0.0 + ].into() ); } } diff --git a/hoomd-vector/src/quaternion.rs b/hoomd-vector/src/quaternion.rs index 5c63e2646..5d4ce3732 100644 --- a/hoomd-vector/src/quaternion.rs +++ b/hoomd-vector/src/quaternion.rs @@ -536,6 +536,15 @@ impl Versor { }) } + /// Get the components of a [`Versor`]. + #[inline] + #[must_use] + pub fn get_components(&self) -> (f64, f64, f64, f64) { + let quaternion = self.get(); + let vec_components = quaternion.vector; + (quaternion.scalar, vec_components[0], vec_components[1], vec_components[2]) + } + /// Normalize the versor. /// /// Nominally, all [`Versor`] instances retain a unit norm. Due to limited From c4b818d24f318bc357aa7e970832ffc576714e1f Mon Sep 17 00:00:00 2001 From: mthran Date: Wed, 17 Dec 2025 17:13:55 -0500 Subject: [PATCH 02/13] trial moves for 3-sphere --- hoomd-manifold/src/sphere.rs | 56 +++++++++++++++++++++++++++++--- hoomd-mc/src/translate/sphere.rs | 55 ++++++++++++++++++++++++++++++- 2 files changed, 105 insertions(+), 6 deletions(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 78a40e4b5..96ac6078d 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -4,6 +4,7 @@ //! Implement vector and curved manifold types on a sphere. use std::f64::consts::PI; +use std::ops::Mul; use approxim::{approx_derive::RelativeEq, assert_relative_eq}; use rand::{ @@ -12,7 +13,7 @@ use rand::{ }; use hoomd_utility::valid::PositiveReal; -use hoomd_vector::{Cartesian, InnerProduct, Metric, Quaternion, Versor}; +use hoomd_vector::{Cartesian, InnerProduct, Metric, Quaternion, Rotate, Versor}; /// Point on the surface of a sphere. /// @@ -297,11 +298,11 @@ impl Metric for Spherical<4> { /// # Ok(()) /// # } /// ``` -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 { @@ -316,7 +317,7 @@ 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] @@ -348,6 +349,29 @@ 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 radius = self.point.radius; + let max_trans = (self.disk_radius.get()) / radius; + 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); + let transformed_point = Spherical::<4>::from_cartesian_coordinates( + sphere_point.point().mul(radius), radius); + transformed_point + } +} + #[cfg(test)] mod tests { use super::*; @@ -413,7 +437,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; @@ -434,4 +458,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]), 1.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); + } + } } diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index d38b053b4..76e4c4d55 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -63,4 +63,57 @@ impl LocalTrial>> for Translate>> { } } -// 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, 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 radius: f64 = 1.0; + /// let initial_point = Point::new( + /// Spherical::<4>::from_polar_coordinates(1.0, 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_eq!(new_body_properties.position().radius(), radius,); + /// + /// // 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 disk = SphericalDisk { + disk_radius: *self.maximum_distance(), + point: *trial.position_mut(), + }; + *trial.position_mut() = disk.sample(rng); + let rescale = trial.position().radius() / trial.position().point().norm(); + *trial.position_mut() = Spherical::from_cartesian_coordinates( + *trial.position().point() * rescale, + trial.position().radius(), + ); + trial + } +} + +// TODO: test code \ No newline at end of file From 1ac905cedf471e563be78857a49feae721a517f1 Mon Sep 17 00:00:00 2001 From: mthran Date: Thu, 18 Dec 2025 13:13:31 -0500 Subject: [PATCH 03/13] Update sphere.rs --- hoomd-manifold/src/sphere.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 96ac6078d..e8e9ed703 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -178,7 +178,6 @@ impl Spherical<4> { 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]); - println!("phi {} theta {} psi {}", phi, theta, psi); 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) } From 3e2b0a0f0a646b3ad71edcf3b95a9e3de295cc38 Mon Sep 17 00:00:00 2001 From: mthran Date: Mon, 16 Feb 2026 12:46:29 -0500 Subject: [PATCH 04/13] Merge with trunk --- hoomd-mc/src/translate/sphere.rs | 2 +- hoomd-microstate/src/property/point.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index 76e4c4d55..b868e576c 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -8,7 +8,7 @@ use rand::{Rng, distr::Distribution}; use crate::{LocalTrial, Translate}; use hoomd_manifold::{Spherical, SphericalDisk}; use hoomd_microstate::property::{Point, Position}; -use hoomd_vector::{InnerProduct, Quaternion, Versor}; +use hoomd_vector::InnerProduct; impl LocalTrial>> for Translate>> { /// Propose local trial moves for a body on the surface of a sphere diff --git a/hoomd-microstate/src/property/point.rs b/hoomd-microstate/src/property/point.rs index 9572fe094..ead0ea380 100644 --- a/hoomd-microstate/src/property/point.rs +++ b/hoomd-microstate/src/property/point.rs @@ -6,7 +6,7 @@ use super::Position; use crate::Transform; use hoomd_manifold::{Hyperbolic, Minkowski, Spherical}; -use hoomd_vector::{Cartesian, Quaternion, Versor}; +use hoomd_vector::{Cartesian, Quaternion}; use std::ops::Mul; /// A position in space and nothing more. From 6ae23c9311392575a82f984be694fa8dd1e2fe25 Mon Sep 17 00:00:00 2001 From: mthran Date: Mon, 16 Feb 2026 14:43:18 -0500 Subject: [PATCH 05/13] Change trial moves on sphere to exponential map --- Cargo.lock | 1 + hoomd-mc/Cargo.toml | 1 + hoomd-mc/src/translate/sphere.rs | 69 +++++++++++++++++++++----------- 3 files changed, 47 insertions(+), 24 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index f14947023..409bdaa1d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3196,6 +3196,7 @@ dependencies = [ "hoomd-vector", "itertools 0.14.0", "rand 0.9.2", + "rand_distr 0.5.1", "rayon", "rstest", ] diff --git a/hoomd-mc/Cargo.toml b/hoomd-mc/Cargo.toml index c6c4697be..cfbf8f272 100644 --- a/hoomd-mc/Cargo.toml +++ b/hoomd-mc/Cargo.toml @@ -16,6 +16,7 @@ workspace = true [dependencies] itertools.workspace = true rand.workspace = true +rand_distr.workspace = true rayon.workspace = true hoomd-geometry.workspace = true diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index b868e576c..f9dbf41f5 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; 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}; +use rand_distr::StandardNormal; impl LocalTrial>> for Translate>> { /// Propose local trial moves for a body on the surface of a sphere @@ -24,10 +25,10 @@ impl LocalTrial>> for Translate>> { /// /// # fn main() -> Result<(), Box> { /// let mut rng = StdRng::seed_from_u64(14); - /// let radius: f64 = 2.0; + /// let radius: f64 = 1.0; /// let initial_point = Point::new(Spherical::from_cartesian_coordinates( - /// [2.0_f64.sqrt(), 2.0_f64.sqrt(), 0.0].into(), - /// 2.0_f64, + /// [0.5_f64.sqrt(), 0.5_f64.sqrt(), 0.0].into(), + /// 1.0_f64, /// )); /// let d = 0.1; /// let translate = Translate::with_maximum_distance(d.try_into()?); @@ -49,16 +50,25 @@ impl LocalTrial>> for Translate>> { #[inline] fn propose(&self, rng: &mut R, 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 = trial.position().radius() / trial.position().point().norm(); + 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()); + println!("proj: {:}", proj); + 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"); + println!("new unit vector: {:?}", unit); + 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( - *trial.position().point() * rescale, - trial.position().radius(), - ); + new, 1.0); trial } } @@ -102,16 +112,27 @@ impl LocalTrial>> for Translate>> { #[inline] fn propose(&self, rng: &mut R, 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 = trial.position().radius() / trial.position().point().norm(); + 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()); + println!("proj: {:}", proj); + 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"); + println!("new unit vector: {:?}", unit); + 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( - *trial.position().point() * rescale, - trial.position().radius(), - ); + new, 1.0); trial } } From 9a1c925080dfb3efc29b61f459dc1f1d511fd860 Mon Sep 17 00:00:00 2001 From: mthran Date: Mon, 16 Feb 2026 16:03:18 -0500 Subject: [PATCH 06/13] Reduce precision on radius check --- hoomd-mc/src/translate/sphere.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index f9dbf41f5..a5d2b619f 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -116,7 +116,6 @@ impl LocalTrial>> for Translate>> { 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()); - println!("proj: {:}", proj); let tangent = Cartesian::from([ vec[0] - proj * trial.position.coordinates()[0], vec[1] - proj * trial.position.coordinates()[1], @@ -124,7 +123,6 @@ impl LocalTrial>> for Translate>> { vec[3] - proj * trial.position.coordinates()[3], ]); let (unit, _norm) = tangent.to_unit().expect("cannot be null"); - println!("new unit vector: {:?}", unit); let new = Cartesian::from([ trial.position.coordinates()[0] * cs + unit.get().coordinates[0]*sn, trial.position.coordinates()[1] * cs + unit.get().coordinates[1]*sn, From da75b054b7a68b368351a4e4cdbfb7d5a28cd63f Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 20 Feb 2026 13:03:51 -0500 Subject: [PATCH 07/13] Change precision on radius check --- hoomd-manifold/src/sphere.rs | 26 +++++++++++++------------- hoomd-mc/src/translate/sphere.rs | 10 ++++------ 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index e8e9ed703..3924d32d0 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -59,7 +59,7 @@ impl Spherical { #[must_use] pub fn from_cartesian_coordinates(point: Cartesian, radius: f64) -> Spherical { let rad = point.norm(); - assert_relative_eq!(rad, radius); + assert_relative_eq!(rad, radius, epsilon=1e-8); Spherical { point, radius } } } @@ -78,7 +78,7 @@ impl Spherical<3> { ]); Spherical::from_cartesian_coordinates(point, r) } - /// Implements a stereographic projection from the 2-sphere to an 2-dimensional + /// Implements a stereographic projection from the 2-sphere to an 2-dimensional /// plane by projecting through the $`(0,0,1)`$ axis. /// /// # Example @@ -112,16 +112,16 @@ impl Spherical<3> { } impl Spherical<4> { - /// Create a point on a 3-sphere from spherical coordinates. Note that this uses - /// the convention - /// ```math - /// \begin{pmatrix}r\cos\psi - /// \\ r\sin\psi\cos\theta - /// \\ r\sin\psi\sin\theta\cos\phi + /// Create a point on a 3-sphere from spherical coordinates. Note that this uses + /// the convention + /// ```math + /// \begin{pmatrix}r\cos\psi + /// \\ r\sin\psi\cos\theta + /// \\ r\sin\psi\sin\theta\cos\phi /// \\ r\sin\psi\sin\theta\sin\phi /// \end{pmatrix} /// ``` - /// where $`\psi`$ and $`theta`$ both run over the range $`0`$ to $`\pi`$ and $`\phi`$ + /// where $`\psi`$ and $`theta`$ both run over the range $`0`$ to $`\pi`$ and $`\phi`$ /// runs from $`0`$ to $`2\pi`$. #[inline] #[must_use] @@ -157,14 +157,14 @@ impl Spherical<4> { /// # fn main() -> Result<(), Box> { /// let radius = 1.0; /// let x = Spherical::<4>::from_polar_coordinates(1.0,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); @@ -181,7 +181,7 @@ impl Spherical<4> { 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) } - /// Implements a stereographic projection from the 2-sphere to an 2-dimensional + /// Implements a stereographic projection from the 2-sphere to an 2-dimensional /// plane by projecting through the $`(1,0,0,0)`$ axis. #[inline] #[must_use] @@ -361,7 +361,7 @@ impl Distribution> for SphericalDisk<4> { 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); diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index a5d2b619f..f81ea2c24 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -51,17 +51,15 @@ impl LocalTrial>> for Translate>> { 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 (sn, cs) = ((displacement.abs()).sin(), (displacement.abs()).cos()); let vec = Cartesian::<3>::from(std::array::from_fn(|_| rng.sample(StandardNormal))); let proj = vec.dot(trial.position.point()); - println!("proj: {:}", proj); 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"); - println!("new unit vector: {:?}", unit); let new = Cartesian::from([ trial.position.coordinates()[0] * cs + unit.get().coordinates[0]*sn, trial.position.coordinates()[1] * cs + unit.get().coordinates[1]*sn, @@ -75,7 +73,7 @@ impl LocalTrial>> for Translate>> { impl LocalTrial>> for Translate>> { /// Propose local trial moves for a body on the surface of a 3-sphere - /// + /// /// # Example /// ``` /// use approxim::assert_relative_eq; @@ -113,7 +111,7 @@ impl LocalTrial>> for Translate>> { 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 (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([ @@ -135,4 +133,4 @@ impl LocalTrial>> for Translate>> { } } -// TODO: test code \ No newline at end of file +// TODO: test code From ee80e136897ef4bb5ac45f594e92acb0b8aaf47a Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Tue, 24 Feb 2026 15:33:06 -0500 Subject: [PATCH 08/13] fixed force on zetterling potential --- .../src/univariate/zetterling.rs | 19 ++++++++----- hoomd-mc/src/translate/sphere.rs | 27 ++++++++++++++++--- 2 files changed, 35 insertions(+), 11 deletions(-) diff --git a/hoomd-interaction/src/univariate/zetterling.rs b/hoomd-interaction/src/univariate/zetterling.rs index d0f5930f4..8dbc8e68d 100644 --- a/hoomd-interaction/src/univariate/zetterling.rs +++ b/hoomd-interaction/src/univariate/zetterling.rs @@ -8,7 +8,7 @@ use super::{UnivariateEnergy, UnivariateForce}; /// `Zetterling` computes the oscillating pair potential between every pair of /// particles in the simulation state. /// ```math -/// U(r) = \epsilon \frac{\exp(\alpha r/\ell)\cos(2k_Fr/\ell)}{(r/\ell)^3} + \beta\left(\frac{\sigma \ell}{r}\right)^n +/// U(r) = \epsilon \frac{\exp(\alpha r/\ell)\cos(2k_Fr/\ell)}{(r/\ell)^3} + \beta\left(\frac{\sigma \ell}{r}\right)^n + v_0 /// ``` #[derive(Clone, Debug, PartialEq)] pub struct Zetterling { @@ -26,6 +26,8 @@ pub struct Zetterling { pub n: f64, /// The length scale of the distances *(\[length\])* pub ell: f64, + /// The energy offset, chosen so the truncated potential is continuous *(\[energy]\)* + pub v_0: f64, } impl UnivariateEnergy for Zetterling { @@ -37,6 +39,7 @@ impl UnivariateEnergy for Zetterling { self.epsilon * ((self.alpha * r_ell).exp()) * ((2.0 * self.kf * r_ell).cos()) / (r_ell.powi(3)) + self.beta * (sigma_ell_r.powf(self.n)) + + self.v_0 } } @@ -44,19 +47,19 @@ impl UnivariateForce for Zetterling { #[inline] fn force(&self, r: f64) -> f64 { let r_ell = r / self.ell; - let r_inv = r.recip(); + let r_ell_inv = self.ell / r; let cos = (2.0 * self.kf * r_ell).cos(); let sin = (2.0 * self.kf * r_ell).sin(); let exp = (self.alpha * r_ell).exp(); - let first = -self.epsilon * self.alpha * (self.ell.powi(2)) * exp * cos * (r_inv.powi(3)); + let first = -self.epsilon * self.alpha * exp * cos * (r_ell_inv.powi(3)); let second = - self.epsilon * 2.0 * self.kf * (self.ell.powi(2)) * exp * sin * (r_inv.powi(3)); - let third = 3.0 * self.epsilon * (self.ell.powi(3)) * exp * cos * (r_inv.powi(4)); + self.epsilon * 2.0 * self.kf * exp * sin * (r_ell_inv.powi(3)); + let third = 3.0 * self.epsilon * exp * cos * (r_ell_inv.powi(4)); let fourth = self.beta * self.n - * ((self.sigma * self.ell).powf(self.n)) - * (r_inv.powf(self.n + 1.0)); + * (self.sigma.powf(self.n)) + * (r_ell_inv.powf(self.n + 1.0)); first + second + third + fourth } @@ -77,6 +80,7 @@ mod tests { sigma: 1.0, n: 18.0, ell: 1.0, + v_0: 0.0, }; let (u_1, r_1) = (-0.742_644_124_392_870, 1.130_547_632_166_212); @@ -95,6 +99,7 @@ mod tests { sigma: 1.0, n: 14.5, ell: 1.0, + v_0: 0.0, }; let (u_1, r_1) = (-0.883_354_387_732_971, 1.132_662_993_677_647); diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index 47a0bb428..997abc6ff 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -4,10 +4,10 @@ //! Implement Translation moves on curved surfaces use rand::{Rng, RngExt}; -use rand_distr::{Distribution, StandardNormal}; +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::{Cartesian, InnerProduct}; @@ -26,7 +26,7 @@ impl LocalTrial>> for Translate>> { /// # fn main() -> Result<(), Box> { /// 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()?); @@ -53,6 +53,25 @@ impl LocalTrial>> for Translate>> { 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::<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(), @@ -61,7 +80,7 @@ impl LocalTrial>> for Translate>> { let rescale = 1.0 / trial.position().point().norm(); *trial.position_mut() = Spherical::from_cartesian_coordinates(*trial.position().point() * rescale); - trial + trial */ } } From 6bd654a337db59b856ca9fef6ab71f0245bf5e11 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Wed, 18 Mar 2026 10:48:26 -0400 Subject: [PATCH 09/13] Fix Windows doc error in biquaternion --- hoomd-bevy/src/representation/hyperbolic_polygon.wgsl | 8 ++++++-- hoomd-interaction/src/univariate/lennard_jones.rs | 5 +++-- hoomd-manifold/src/lib.rs | 11 ++++++++--- 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/hoomd-bevy/src/representation/hyperbolic_polygon.wgsl b/hoomd-bevy/src/representation/hyperbolic_polygon.wgsl index 42c192a4d..c5d9bb694 100644 --- a/hoomd-bevy/src/representation/hyperbolic_polygon.wgsl +++ b/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; @group(2) @binding(5) var image_color_sampler: sampler; -@group(2) @binding(6) var n_sides: f32; +struct PolygonParams { + n_sides: f32, +}; +@group(2) @binding(6) +var 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(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 diff --git a/hoomd-interaction/src/univariate/lennard_jones.rs b/hoomd-interaction/src/univariate/lennard_jones.rs index e862f4b69..ebc24c0c7 100644 --- a/hoomd-interaction/src/univariate/lennard_jones.rs +++ b/hoomd-interaction/src/univariate/lennard_jones.rs @@ -56,7 +56,8 @@ use super::{UnivariateEnergy, UnivariateForce}; /// assert_abs_diff_eq!(lennard_jones.energy(sigma), 0.0); /// assert_relative_eq!( /// lennard_jones.energy(2.0_f64.powf(1.0 / 4.0) * sigma), -/// -epsilon +/// -epsilon, +/// epsilon=1e-12, /// ); /// assert_abs_diff_eq!( /// lennard_jones.force(2.0_f64.powf(1.0 / 4.0) * sigma), @@ -170,7 +171,7 @@ mod tests { assert_relative_eq!(lj.force(sigma), 16.0 * epsilon / sigma); // Bottom of the well - assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma), -epsilon); + assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma), -epsilon, epsilon=1e-8); assert_abs_diff_eq!( lj.force(2.0_f64.powf(1.0 / 4.0) * sigma), 0.0, diff --git a/hoomd-manifold/src/lib.rs b/hoomd-manifold/src/lib.rs index 44faa98f2..2e34e9e97 100644 --- a/hoomd-manifold/src/lib.rs +++ b/hoomd-manifold/src/lib.rs @@ -136,18 +136,23 @@ //! }; //! use num::complex::Complex; //! use std::f64::consts::PI; +//! use approxim::assert_relative_eq; //! //! # fn main() -> Result<(), Box> { //! let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]); //! let q = Biquaternion::from([ -//! Complex::new(0.0, PI / 4.0).sinh(), +//! Complex::new(0.0, (0.2_f64).sinh()), //! Complex::new(0.0, 0.0), //! Complex::new(0.0, 0.0), -//! Complex::new(0.0, PI / 4.0).cosh(), +//! Complex::new((0.2_f64).cosh(), 0.0), //! ]); //! let v = q.to_unit()?; //! let boosted = v.hyperbolic_rotate(&x); -//! // boosted is approximately [(PI/2.0).sinh(), 0.0, 0.0, (PI/2.0).cosh()] +//! assert_relative_eq!( +//! boosted, +//! [(0.4_f64).sinh(), 0.0, 0.0, (0.4_f64).cosh()].into(), +//! epsilon = 1e-12 +//! ); //! # Ok(()) //! # } //! ``` From 656b6b88c5fde5049a82753256d50af1b6a18219 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Tue, 31 Mar 2026 13:23:16 -0400 Subject: [PATCH 10/13] gsd file for spherical points --- hoomd-manifold/src/sphere.rs | 3 +- hoomd-microstate/src/append.rs | 69 +++++++++++++++++++++++++++++++++- 2 files changed, 70 insertions(+), 2 deletions(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 8a8b2b05a..e6839e69a 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -210,7 +210,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 = if arg >= 1.0 {1.0} else if arg <= -1.0 {-1.0} else {arg}; + arg_clipped.acos() } #[inline] fn distance_squared(&self, other: &Self) -> f64 { diff --git a/hoomd-microstate/src/append.rs b/hoomd-microstate/src/append.rs index f4ba09d2a..5715a7b3d 100644 --- a/hoomd-microstate/src/append.rs +++ b/hoomd-microstate/src/append.rs @@ -5,11 +5,12 @@ use hoomd_geometry::shape::Hypercuboid; use hoomd_gsd::hoomd::{AppendError, Dimensions, Frame, HoomdGsdFile}; +use hoomd_manifold::Spherical; use hoomd_vector::{Angle, Cartesian, Versor}; use crate::{ AppendMicrostate, Microstate, - boundary::{Closed, Periodic}, + boundary::{Closed, Open, Periodic}, property::{OrientedPoint, Point}, }; @@ -202,8 +203,29 @@ impl AppendMicrostate, Versor>, X, Periodic< } } +impl AppendMicrostate>, X, Open> for HoomdGsdFile { + #[inline] + fn append_microstate( + &mut self, + microstate: &Microstate>, X, Open>, + ) -> Result, AppendError> { + self.append_frame(microstate.step())? + .configuration_box([5.0, 5.0, 5.0, 0.0, 0.0, 0.0])? + .configuration_dimensions(Dimensions::Three)? + .particles_position( + microstate + .iter_sites_tag_order() + .map(|s| { + let proj: Vec = s.properties.position.stereographic_projection(); + [proj[0], proj[1], proj[2]].into() + }), + ) + } +} + #[cfg(test)] mod test { + use approxim::assert_relative_eq; use assert2::assert; use std::f64::consts::PI; use tempfile::tempdir; @@ -298,6 +320,51 @@ mod test { Ok(()) } + #[test] + fn point_hypersphere() -> anyhow::Result<()> { + let microstate = Microstate::builder() + .boundary(Open) + .bodies([ + Body::point(Spherical::from_cartesian_coordinates( + Cartesian::from([1.0, 0.0, 0.0, 0.0])) + ), + Body::point(Spherical::from_cartesian_coordinates( + Cartesian::from([0.0, 1.0, 0.0, 0.0])) + ), + Body::point(Spherical::from_cartesian_coordinates( + Cartesian::from([0.0, 0.0, 1.0, 0.0])) + ), + ]) + .step(1234) + .try_build()?; + + let tmp_dir = tempdir()?; + let path = tmp_dir.path().join("test.gsd"); + let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?; + hoomd_gsd_file.append_microstate(µstate)?; + + drop(hoomd_gsd_file); + + let gsd_file = GsdFile::open(path, Mode::Read)?; + + assert!(gsd_file.n_frames() == 1); + + let step = gsd_file.iter_scalars::(0, "configuration/step")?; + itertools::assert_equal(step, [1234]); + + let positions: Vec<[f32;3]> = gsd_file.iter_arrays::(0, "particles/position")?.collect(); + assert_relative_eq!(positions[0][0], 1.0); + assert_relative_eq!(positions[0][1], 0.0); + assert_relative_eq!(positions[0][2], 0.0); + assert_relative_eq!(positions[1][0], 0.0); + assert_relative_eq!(positions[1][1], 1.0); + assert_relative_eq!(positions[1][2], 0.0); + assert_relative_eq!(positions[2][0], 0.0); + assert_relative_eq!(positions[2][1], 0.0); + assert_relative_eq!(positions[2][2], 1.0); + Ok(()) + } + #[test] fn oriented_point_closed_rectangle_2d() -> anyhow::Result<()> { let boundary = Closed(Rectangle { From 9ef4e24d17f889e807c80b49464f57c97d3bdbe2 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 3 Apr 2026 12:35:43 -0400 Subject: [PATCH 11/13] append for hyperbolic points --- hoomd-bevy/katex.html | 2 +- hoomd-derive/katex.html | 2 +- hoomd-geometry/katex.html | 2 +- hoomd-gsd/katex.html | 2 +- hoomd-interaction/katex.html | 2 +- .../src/univariate/lennard_jones.rs | 8 +- .../src/univariate/zetterling.rs | 9 +- hoomd-linear-algebra/katex.html | 2 +- hoomd-manifold/katex.html | 2 +- hoomd-manifold/src/lib.rs | 2 +- hoomd-manifold/src/sphere.rs | 8 +- hoomd-mc/katex.html | 2 +- hoomd-mc/src/translate/sphere.rs | 33 +- hoomd-microstate/src/append.rs | 301 +++++++++++++++++- hoomd-rand/katex.html | 2 +- hoomd-simulation/katex.html | 2 +- hoomd-spatial/katex.html | 2 +- hoomd-utility/katex.html | 2 +- hoomd-vector/katex.html | 2 +- 19 files changed, 330 insertions(+), 57 deletions(-) diff --git a/hoomd-bevy/katex.html b/hoomd-bevy/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-bevy/katex.html +++ b/hoomd-bevy/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-derive/katex.html b/hoomd-derive/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-derive/katex.html +++ b/hoomd-derive/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-geometry/katex.html b/hoomd-geometry/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-geometry/katex.html +++ b/hoomd-geometry/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-gsd/katex.html b/hoomd-gsd/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-gsd/katex.html +++ b/hoomd-gsd/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-interaction/katex.html b/hoomd-interaction/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-interaction/katex.html +++ b/hoomd-interaction/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-interaction/src/univariate/lennard_jones.rs b/hoomd-interaction/src/univariate/lennard_jones.rs index ebc24c0c7..79cdc1451 100644 --- a/hoomd-interaction/src/univariate/lennard_jones.rs +++ b/hoomd-interaction/src/univariate/lennard_jones.rs @@ -57,7 +57,7 @@ use super::{UnivariateEnergy, UnivariateForce}; /// assert_relative_eq!( /// lennard_jones.energy(2.0_f64.powf(1.0 / 4.0) * sigma), /// -epsilon, -/// epsilon=1e-12, +/// epsilon = 1e-12, /// ); /// assert_abs_diff_eq!( /// lennard_jones.force(2.0_f64.powf(1.0 / 4.0) * sigma), @@ -171,7 +171,11 @@ mod tests { assert_relative_eq!(lj.force(sigma), 16.0 * epsilon / sigma); // Bottom of the well - assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma), -epsilon, epsilon=1e-8); + assert_relative_eq!( + lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma), + -epsilon, + epsilon = 1e-8 + ); assert_abs_diff_eq!( lj.force(2.0_f64.powf(1.0 / 4.0) * sigma), 0.0, diff --git a/hoomd-interaction/src/univariate/zetterling.rs b/hoomd-interaction/src/univariate/zetterling.rs index 8dbc8e68d..7e50c8be9 100644 --- a/hoomd-interaction/src/univariate/zetterling.rs +++ b/hoomd-interaction/src/univariate/zetterling.rs @@ -53,13 +53,10 @@ impl UnivariateForce for Zetterling { let exp = (self.alpha * r_ell).exp(); let first = -self.epsilon * self.alpha * exp * cos * (r_ell_inv.powi(3)); - let second = - self.epsilon * 2.0 * self.kf * exp * sin * (r_ell_inv.powi(3)); + let second = self.epsilon * 2.0 * self.kf * exp * sin * (r_ell_inv.powi(3)); let third = 3.0 * self.epsilon * exp * cos * (r_ell_inv.powi(4)); - let fourth = self.beta - * self.n - * (self.sigma.powf(self.n)) - * (r_ell_inv.powf(self.n + 1.0)); + let fourth = + self.beta * self.n * (self.sigma.powf(self.n)) * (r_ell_inv.powf(self.n + 1.0)); first + second + third + fourth } diff --git a/hoomd-linear-algebra/katex.html b/hoomd-linear-algebra/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-linear-algebra/katex.html +++ b/hoomd-linear-algebra/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-manifold/katex.html b/hoomd-manifold/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-manifold/katex.html +++ b/hoomd-manifold/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-manifold/src/lib.rs b/hoomd-manifold/src/lib.rs index 672709ee1..a95000838 100644 --- a/hoomd-manifold/src/lib.rs +++ b/hoomd-manifold/src/lib.rs @@ -131,12 +131,12 @@ //! Boost point in 3D hyperbolic space in x direction using biquaternion //! algebra: //! ``` +//! use approxim::assert_relative_eq; //! use hoomd_manifold::{ //! Biquaternion, HyperbolicRotate, Minkowski, UnitBiquaternion, //! }; //! use num::complex::Complex; //! use std::f64::consts::PI; -//! use approxim::assert_relative_eq; //! //! # fn main() -> Result<(), Box> { //! let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]); diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 0401f6a75..605336e8a 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -210,7 +210,13 @@ impl Metric for Spherical<4> { #[inline] fn distance(&self, other: &Self) -> f64 { let arg = Cartesian::dot(&self.point, &other.point); - let arg_clipped = if arg >= 1.0 {1.0} else if arg <= -1.0 {-1.0} else {arg}; + let arg_clipped = if arg >= 1.0 { + 1.0 + } else if arg <= -1.0 { + -1.0 + } else { + arg + }; arg_clipped.acos() } #[inline] diff --git a/hoomd-mc/katex.html b/hoomd-mc/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-mc/katex.html +++ b/hoomd-mc/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-mc/src/translate/sphere.rs b/hoomd-mc/src/translate/sphere.rs index 997abc6ff..edf6930bb 100644 --- a/hoomd-mc/src/translate/sphere.rs +++ b/hoomd-mc/src/translate/sphere.rs @@ -26,7 +26,7 @@ impl LocalTrial>> for Translate>> { /// # fn main() -> Result<(), Box> { /// let mut rng = StdRng::seed_from_u64(14); /// let initial_point = Point::new(Spherical::from_cartesian_coordinates( - /// [0.5_f64.sqrt(),0.0, 0.5_f64.sqrt()].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()?); @@ -53,7 +53,7 @@ impl LocalTrial>> for Translate>> { body_properties: Point>, ) -> Point> { let mut trial = body_properties; - let displacement = (self.maximum_distance().get())*rng.sample::(StandardNormal); + 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()); @@ -64,23 +64,22 @@ impl LocalTrial>> for Translate>> { ]); 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()[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.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 */ + // 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 } } diff --git a/hoomd-microstate/src/append.rs b/hoomd-microstate/src/append.rs index 5715a7b3d..be50ef7dc 100644 --- a/hoomd-microstate/src/append.rs +++ b/hoomd-microstate/src/append.rs @@ -5,13 +5,13 @@ use hoomd_geometry::shape::Hypercuboid; use hoomd_gsd::hoomd::{AppendError, Dimensions, Frame, HoomdGsdFile}; -use hoomd_manifold::Spherical; +use hoomd_manifold::{Hyperbolic, Spherical}; use hoomd_vector::{Angle, Cartesian, Versor}; use crate::{ AppendMicrostate, Microstate, boundary::{Closed, Open, Periodic}, - property::{OrientedPoint, Point}, + property::{OrientedHyperbolicPoint, OrientedPoint, Point}, }; impl AppendMicrostate>, X, Closed>> for HoomdGsdFile { @@ -212,17 +212,62 @@ impl AppendMicrostate>, X, Open> for HoomdGsdFile { self.append_frame(microstate.step())? .configuration_box([5.0, 5.0, 5.0, 0.0, 0.0, 0.0])? .configuration_dimensions(Dimensions::Three)? + .particles_position(microstate.iter_sites_tag_order().map(|s| { + let proj: Vec = s.properties.position.stereographic_projection(); + [proj[0], proj[1], proj[2]].into() + })) + } +} + +impl AppendMicrostate, X, C> for HoomdGsdFile { + #[inline] + fn append_microstate( + &mut self, + microstate: &Microstate, X, C>, + ) -> Result, AppendError> { + self.append_frame(microstate.step())? + .configuration_box([2.0, 2.0, 2.0, 0.0, 0.0, 0.0])? + .configuration_dimensions(Dimensions::Two)? .particles_position( microstate .iter_sites_tag_order() - .map(|s| { - let proj: Vec = s.properties.position.stereographic_projection(); - [proj[0], proj[1], proj[2]].into() + .map(|s| s.properties.position.to_poincare()) + .map(|p| [p[0], p[1], 0.0].into()), + )? + .particles_orientation( + microstate + .iter_sites_tag_order() + .map(|s| s.properties.orientation.theta) + .map(|a| { + Versor::from_axis_angle( + [0.0, 0.0, 1.0] + .try_into() + .expect("hard-coded vector can be normalized"), + a, + ) }), ) } } +impl AppendMicrostate>, X, C> for HoomdGsdFile { + #[inline] + fn append_microstate( + &mut self, + microstate: &Microstate>, X, C>, + ) -> Result, AppendError> { + self.append_frame(microstate.step())? + .configuration_box([2.0, 2.0, 2.0, 0.0, 0.0, 0.0])? + .configuration_dimensions(Dimensions::Two)? + .particles_position( + microstate + .iter_sites_tag_order() + .map(|s| s.properties.position.to_poincare()) + .map(|p| [p[0], p[1], 0.0].into()), + ) + } +} + #[cfg(test)] mod test { use approxim::assert_relative_eq; @@ -232,7 +277,7 @@ mod test { use super::*; use crate::Body; - use hoomd_geometry::shape::Rectangle; + use hoomd_geometry::shape::{EightEight, Rectangle}; use hoomd_gsd::file_layer::{GsdFile, Mode}; #[test] @@ -321,19 +366,19 @@ mod test { } #[test] - fn point_hypersphere() -> anyhow::Result<()> { + fn point_spherical_3d() -> anyhow::Result<()> { let microstate = Microstate::builder() .boundary(Open) .bodies([ - Body::point(Spherical::from_cartesian_coordinates( - Cartesian::from([1.0, 0.0, 0.0, 0.0])) - ), - Body::point(Spherical::from_cartesian_coordinates( - Cartesian::from([0.0, 1.0, 0.0, 0.0])) - ), - Body::point(Spherical::from_cartesian_coordinates( - Cartesian::from([0.0, 0.0, 1.0, 0.0])) - ), + Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([ + 1.0, 0.0, 0.0, 0.0, + ]))), + Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([ + 0.0, 1.0, 0.0, 0.0, + ]))), + Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([ + 0.0, 0.0, 1.0, 0.0, + ]))), ]) .step(1234) .try_build()?; @@ -352,7 +397,9 @@ mod test { let step = gsd_file.iter_scalars::(0, "configuration/step")?; itertools::assert_equal(step, [1234]); - let positions: Vec<[f32;3]> = gsd_file.iter_arrays::(0, "particles/position")?.collect(); + let positions: Vec<[f32; 3]> = gsd_file + .iter_arrays::(0, "particles/position")? + .collect(); assert_relative_eq!(positions[0][0], 1.0); assert_relative_eq!(positions[0][1], 0.0); assert_relative_eq!(positions[0][2], 0.0); @@ -365,6 +412,226 @@ mod test { Ok(()) } + #[test] + fn point_hyperbolic_2d_open() -> anyhow::Result<()> { + let microstate = Microstate::builder() + .boundary(Open) + .bodies([ + Body::point(Hyperbolic::<3>::from_polar_coordinates(1.2, 0.0)), + Body::point(Hyperbolic::<3>::from_polar_coordinates(0.6, 1.5)), + ]) + .step(1234) + .try_build()?; + + let tmp_dir = tempdir()?; + let path = tmp_dir.path().join("test.gsd"); + let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?; + hoomd_gsd_file.append_microstate(µstate)?; + + drop(hoomd_gsd_file); + + let gsd_file = GsdFile::open(path, Mode::Read)?; + + assert!(gsd_file.n_frames() == 1); + + let step = gsd_file.iter_scalars::(0, "configuration/step")?; + itertools::assert_equal(step, [1234]); + + let positions: Vec<[f32; 3]> = gsd_file + .iter_arrays::(0, "particles/position")? + .collect(); + assert_relative_eq!(positions[0][0], (f32::sinh(1.2)) / (1.0 + f32::cosh(1.2))); + assert_relative_eq!(positions[0][1], 0.0); + assert_relative_eq!( + positions[1][0], + (f32::sinh(0.6) * f32::cos(1.5)) / (1.0 + f32::cosh(0.6)) + ); + assert_relative_eq!( + positions[1][1], + (f32::sinh(0.6) * f32::sin(1.5)) / (1.0 + f32::cosh(0.6)) + ); + Ok(()) + } + + #[test] + fn point_hyperbolic_2d_eighteight_periodic() -> anyhow::Result<()> { + let boundary = Periodic::new(0.6, EightEight {})?; + let microstate = Microstate::builder() + .boundary(boundary) + .bodies([ + Body::point(Hyperbolic::<3>::from_polar_coordinates(0.2, 0.3)), + Body::point(Hyperbolic::<3>::from_polar_coordinates(0.6, 0.7)), + ]) + .step(1234) + .try_build()?; + + let tmp_dir = tempdir()?; + let path = tmp_dir.path().join("test.gsd"); + let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?; + hoomd_gsd_file.append_microstate(µstate)?; + + drop(hoomd_gsd_file); + + let gsd_file = GsdFile::open(path, Mode::Read)?; + + assert!(gsd_file.n_frames() == 1); + + let step = gsd_file.iter_scalars::(0, "configuration/step")?; + itertools::assert_equal(step, [1234]); + + let positions: Vec<[f32; 3]> = gsd_file + .iter_arrays::(0, "particles/position")? + .collect(); + assert_relative_eq!( + positions[0][0], + (f32::sinh(0.2) * f32::cos(0.3)) / (1.0 + f32::cosh(0.2)) + ); + assert_relative_eq!( + positions[0][1], + (f32::sinh(0.2) * f32::sin(0.3)) / (1.0 + f32::cosh(0.2)) + ); + assert_relative_eq!( + positions[1][0], + (f32::sinh(0.6) * f32::cos(0.7)) / (1.0 + f32::cosh(0.6)) + ); + assert_relative_eq!( + positions[1][1], + (f32::sinh(0.6) * f32::sin(0.7)) / (1.0 + f32::cosh(0.6)) + ); + Ok(()) + } + + #[test] + fn oriented_point_hyperbolic_2d_open() -> anyhow::Result<()> { + let microstate = Microstate::builder() + .boundary(Open) + .bodies([ + Body { + properties: OrientedHyperbolicPoint { + position: Hyperbolic::<3>::from_polar_coordinates(0.5, 0.6), + orientation: Angle::from(0.3), + }, + sites: vec![OrientedHyperbolicPoint { + position: Hyperbolic::<3>::default(), + orientation: Angle::default(), + }], + }, + Body { + properties: OrientedHyperbolicPoint { + position: Hyperbolic::<3>::from_polar_coordinates(0.9, 0.3), + orientation: Angle::from(1.2), + }, + sites: vec![OrientedHyperbolicPoint { + position: Hyperbolic::<3>::default(), + orientation: Angle::default(), + }], + }, + ]) + .step(1234) + .try_build()?; + + let tmp_dir = tempdir()?; + let path = tmp_dir.path().join("test.gsd"); + let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?; + hoomd_gsd_file.append_microstate(µstate)?; + + drop(hoomd_gsd_file); + + let gsd_file = GsdFile::open(path, Mode::Read)?; + + assert!(gsd_file.n_frames() == 1); + + let step = gsd_file.iter_scalars::(0, "configuration/step")?; + itertools::assert_equal(step, [1234]); + + let positions: Vec<[f32; 3]> = gsd_file + .iter_arrays::(0, "particles/position")? + .collect(); + assert_relative_eq!( + positions[0][0], + (f32::sinh(0.5) * f32::cos(0.6)) / (1.0 + f32::cosh(0.5)) + ); + assert_relative_eq!( + positions[0][1], + (f32::sinh(0.5) * f32::sin(0.6)) / (1.0 + f32::cosh(0.5)) + ); + assert_relative_eq!( + positions[1][0], + (f32::sinh(0.9) * f32::cos(0.3)) / (1.0 + f32::cosh(0.9)) + ); + assert_relative_eq!( + positions[1][1], + (f32::sinh(0.9) * f32::sin(0.3)) / (1.0 + f32::cosh(0.9)) + ); + Ok(()) + } + + #[test] + fn oriented_point_hyperbolic_2d_periodic_eighteight() -> anyhow::Result<()> { + let boundary = Periodic::new(0.6, EightEight {})?; + let microstate = Microstate::builder() + .boundary(boundary) + .bodies([ + Body { + properties: OrientedHyperbolicPoint { + position: Hyperbolic::<3>::from_polar_coordinates(0.5, 0.6), + orientation: Angle::from(0.3), + }, + sites: vec![OrientedHyperbolicPoint { + position: Hyperbolic::<3>::default(), + orientation: Angle::default(), + }], + }, + Body { + properties: OrientedHyperbolicPoint { + position: Hyperbolic::<3>::from_polar_coordinates(0.9, 0.3), + orientation: Angle::from(1.2), + }, + sites: vec![OrientedHyperbolicPoint { + position: Hyperbolic::<3>::default(), + orientation: Angle::default(), + }], + }, + ]) + .step(1234) + .try_build()?; + + let tmp_dir = tempdir()?; + let path = tmp_dir.path().join("test.gsd"); + let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?; + hoomd_gsd_file.append_microstate(µstate)?; + + drop(hoomd_gsd_file); + + let gsd_file = GsdFile::open(path, Mode::Read)?; + + assert!(gsd_file.n_frames() == 1); + + let step = gsd_file.iter_scalars::(0, "configuration/step")?; + itertools::assert_equal(step, [1234]); + + let positions: Vec<[f32; 3]> = gsd_file + .iter_arrays::(0, "particles/position")? + .collect(); + assert_relative_eq!( + positions[0][0], + (f32::sinh(0.5) * f32::cos(0.6)) / (1.0 + f32::cosh(0.5)) + ); + assert_relative_eq!( + positions[0][1], + (f32::sinh(0.5) * f32::sin(0.6)) / (1.0 + f32::cosh(0.5)) + ); + assert_relative_eq!( + positions[1][0], + (f32::sinh(0.9) * f32::cos(0.3)) / (1.0 + f32::cosh(0.9)) + ); + assert_relative_eq!( + positions[1][1], + (f32::sinh(0.9) * f32::sin(0.3)) / (1.0 + f32::cosh(0.9)) + ); + Ok(()) + } + #[test] fn oriented_point_closed_rectangle_2d() -> anyhow::Result<()> { let boundary = Closed(Rectangle { diff --git a/hoomd-rand/katex.html b/hoomd-rand/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-rand/katex.html +++ b/hoomd-rand/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-simulation/katex.html b/hoomd-simulation/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-simulation/katex.html +++ b/hoomd-simulation/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-spatial/katex.html b/hoomd-spatial/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-spatial/katex.html +++ b/hoomd-spatial/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-utility/katex.html b/hoomd-utility/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-utility/katex.html +++ b/hoomd-utility/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html diff --git a/hoomd-vector/katex.html b/hoomd-vector/katex.html index 9917b05c1..f06a1b81a 120000 --- a/hoomd-vector/katex.html +++ b/hoomd-vector/katex.html @@ -1 +1 @@ -../katex.html \ No newline at end of file +../katex.html From 57fcf9c9df65b32e090c2995f779ddd40f4e9df3 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 3 Apr 2026 14:07:16 -0400 Subject: [PATCH 12/13] remove zetterling files --- hoomd-interaction/src/univariate.rs | 2 - .../src/univariate/zetterling.rs | 110 ------------------ hoomd-manifold/src/sphere.rs | 2 +- 3 files changed, 1 insertion(+), 113 deletions(-) delete mode 100644 hoomd-interaction/src/univariate/zetterling.rs diff --git a/hoomd-interaction/src/univariate.rs b/hoomd-interaction/src/univariate.rs index 66980638c..2bb4d3d5f 100644 --- a/hoomd-interaction/src/univariate.rs +++ b/hoomd-interaction/src/univariate.rs @@ -13,7 +13,6 @@ mod overlap_penalty; mod shifted; mod weeks_chandler_anderson; mod xplor; -mod zetterling; pub use boxcar::Boxcar; pub use expanded::Expanded; @@ -25,7 +24,6 @@ pub use overlap_penalty::OverlapPenalty; pub use shifted::Shifted; pub use weeks_chandler_anderson::WeeksChandlerAnderson; pub use xplor::Xplor; -pub use zetterling::Zetterling; /// Computes energy as a function of one variable. /// diff --git a/hoomd-interaction/src/univariate/zetterling.rs b/hoomd-interaction/src/univariate/zetterling.rs deleted file mode 100644 index 7e50c8be9..000000000 --- a/hoomd-interaction/src/univariate/zetterling.rs +++ /dev/null @@ -1,110 +0,0 @@ -// Copyright (c) 2024-2026 The Regents of the University of Michigan. -// Part of hoomd-rs, released under the BSD 3-Clause License. - -//! Implement [`Zetterling`] - -use super::{UnivariateEnergy, UnivariateForce}; - -/// `Zetterling` computes the oscillating pair potential between every pair of -/// particles in the simulation state. -/// ```math -/// U(r) = \epsilon \frac{\exp(\alpha r/\ell)\cos(2k_Fr/\ell)}{(r/\ell)^3} + \beta\left(\frac{\sigma \ell}{r}\right)^n + v_0 -/// ``` -#[derive(Clone, Debug, PartialEq)] -pub struct Zetterling { - /// Energy scale of the first term *(\[energy\])*. - pub epsilon: f64, - /// Screening factor *(\[unitless\])*. - pub alpha: f64, - /// Wave number to mimic the Friedel oscillations effect *(\[unitless\])*. - pub kf: f64, - /// Energy scale of the second term *(\[energy\])*. - pub beta: f64, - /// Repulsive core size *(\[unitless\])*. - pub sigma: f64, - /// The power to take sigma/r in the second term *(\[unitless\])*. - pub n: f64, - /// The length scale of the distances *(\[length\])* - pub ell: f64, - /// The energy offset, chosen so the truncated potential is continuous *(\[energy]\)* - pub v_0: f64, -} - -impl UnivariateEnergy for Zetterling { - #[inline] - fn energy(&self, r: f64) -> f64 { - let sigma_ell_r = self.sigma * self.ell / r; - let r_ell = r / self.ell; - - self.epsilon * ((self.alpha * r_ell).exp()) * ((2.0 * self.kf * r_ell).cos()) - / (r_ell.powi(3)) - + self.beta * (sigma_ell_r.powf(self.n)) - + self.v_0 - } -} - -impl UnivariateForce for Zetterling { - #[inline] - fn force(&self, r: f64) -> f64 { - let r_ell = r / self.ell; - let r_ell_inv = self.ell / r; - let cos = (2.0 * self.kf * r_ell).cos(); - let sin = (2.0 * self.kf * r_ell).sin(); - let exp = (self.alpha * r_ell).exp(); - - let first = -self.epsilon * self.alpha * exp * cos * (r_ell_inv.powi(3)); - let second = self.epsilon * 2.0 * self.kf * exp * sin * (r_ell_inv.powi(3)); - let third = 3.0 * self.epsilon * exp * cos * (r_ell_inv.powi(4)); - let fourth = - self.beta * self.n * (self.sigma.powf(self.n)) * (r_ell_inv.powf(self.n + 1.0)); - - first + second + third + fourth - } -} - -#[cfg(test)] -mod tests { - use super::*; - use approxim::assert_relative_eq; - - #[test] - fn select_points() { - let zetterling1 = Zetterling { - epsilon: 1.58, - alpha: -0.22, - kf: 4.12, - beta: 0.95533, - sigma: 1.0, - n: 18.0, - ell: 1.0, - v_0: 0.0, - }; - - let (u_1, r_1) = (-0.742_644_124_392_870, 1.130_547_632_166_212); - assert_relative_eq!(u_1, zetterling1.energy(r_1), epsilon = 1e-12); - assert_relative_eq!(0.0, zetterling1.force(r_1), epsilon = 1e-6); - - let (u_2, r_2) = (-0.153_547_971_965_573, 1.879_994_132_512_336); - assert_relative_eq!(u_2, zetterling1.energy(r_2), epsilon = 1e-12); - assert_relative_eq!(0.0, zetterling1.force(r_2), epsilon = 1e-10); - - let zetterling2 = Zetterling { - epsilon: 1.04, - alpha: 0.33, - kf: 4.139, - beta: 0.94656, - sigma: 1.0, - n: 14.5, - ell: 1.0, - v_0: 0.0, - }; - - let (u_1, r_1) = (-0.883_354_387_732_971, 1.132_662_993_677_647); - assert_relative_eq!(u_1, zetterling2.energy(r_1), epsilon = 1e-12); - assert_relative_eq!(0.0, zetterling2.force(r_1), epsilon = 1e-6); - - let (u_2, r_2) = (-0.287_906_300_605_616, 1.879_255_603_675_727); - assert_relative_eq!(u_2, zetterling2.energy(r_2), epsilon = 1e-12); - assert_relative_eq!(0.0, zetterling2.force(r_2), epsilon = 1e-6); - } -} diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 605336e8a..5f0c0f5dc 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -269,7 +269,7 @@ 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 { From 9b6173ba003dae8a3c2e904e2a73a85ced453204 Mon Sep 17 00:00:00 2001 From: Michelle Thran Date: Fri, 3 Apr 2026 15:45:31 -0400 Subject: [PATCH 13/13] Fix clippy warnings --- hoomd-manifold/src/sphere.rs | 53 +++++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/hoomd-manifold/src/sphere.rs b/hoomd-manifold/src/sphere.rs index 5f0c0f5dc..eaa91a6c5 100644 --- a/hoomd-manifold/src/sphere.rs +++ b/hoomd-manifold/src/sphere.rs @@ -125,26 +125,48 @@ impl Spherical<4> { /// 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 approxim::assert_relative_eq; /// 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 = 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()) + /// let pole_versor = Quaternion::from([1.0, 0.0, 0.0, 0.0]) /// .to_versor() - /// .expect("Hard-coded example is valid"); + /// .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); + /// 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(()) /// # } /// ``` @@ -210,13 +232,7 @@ impl Metric for Spherical<4> { #[inline] fn distance(&self, other: &Self) -> f64 { let arg = Cartesian::dot(&self.point, &other.point); - let arg_clipped = if arg >= 1.0 { - 1.0 - } else if arg <= -1.0 { - -1.0 - } else { - arg - }; + let arg_clipped = arg.clamp(-1.0, 1.0); arg_clipped.acos() } #[inline] @@ -357,8 +373,7 @@ impl Distribution> for SphericalDisk<4> { .to_versor() .expect("sperical points cannot be null"); let sphere_point = Spherical::<4>::from_versor(transformation); - let transformed_point = Spherical::<4>::from_cartesian_coordinates(*sphere_point.point()); - transformed_point + Spherical::<4>::from_cartesian_coordinates(*sphere_point.point()) } }