hoomd_mc/rotate/
versor.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement Rotate for Versor
5
6use rand::{Rng, distr::Distribution};
7use rand_distr::Normal;
8use std::f64::consts::PI;
9
10use hoomd_microstate::property::Orientation;
11use hoomd_utility::valid::PositiveReal;
12use hoomd_vector::{Cartesian, InnerProduct, Quaternion, Rotation, Versor};
13
14use crate::{Adjust, LocalTrial, Rotate};
15
16/// A normal distribution of random Versors, centered on a mean with
17/// some standard deviation.
18#[derive(Debug, Clone, Copy)]
19pub(crate) struct VersorDisplacement {
20    /// The standard deviation of the normal distribution of quaternions around the mean.
21    std_dev: f64,
22}
23
24impl From<f64> for VersorDisplacement {
25    #[inline]
26    fn from(value: f64) -> Self {
27        Self { std_dev: value }
28    }
29}
30impl Distribution<Versor> for VersorDisplacement {
31    /// Sample a random [`Versor`] displacement from a provided mean.
32    ///
33    /// Mathematically, we sample from a 3-dimensional Normal distribution
34    /// in the tangent space of SO(3), lift to the manifold, then rotate to center on
35    /// the mean of the [`VersorDisplacement`]. The result is a small displacement from
36    /// a quaternion input, with fast decay in the tails that make large displacements
37    /// unlikely. This is desirable for Monte Carlo, as large moves are very likely to
38    /// be rejected.
39    #[inline]
40    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Versor {
41        // Based on Karney 2007: doi.org/10.1016/j.jmgm.2006.04.002
42        const _SAFE_THETA: f64 = 1e-12;
43        loop {
44            // As in section 7, we select `s` from a 3D Gaussian distribution.
45            let normal =
46                Normal::new(0.0, self.std_dev).expect("Failed to create normal distribution.");
47
48            let s = Cartesian::from([(); 3].map(|()| normal.sample(rng)));
49            let theta = s.norm();
50
51            // Reject moves with |s| > pi to ensure detailed balance. Otherwise, there
52            // is a shorter path between the proposed move and the start (theta - pi),
53            // which has a different rejection probability.
54            if theta > PI {
55                continue;
56            }
57
58            // Lift the normally distributed values to SO(3) with the exponential map
59            let half_theta = 0.5 * theta;
60            let w = half_theta.cos();
61
62            let mut v_factor = half_theta.sin() / theta;
63            if !v_factor.is_finite() {
64                v_factor = 0.5;
65            }
66
67            let v = s * v_factor;
68
69            // We are normalized by construction, so call the tuple initializer
70            return (Quaternion::from([w, v[0], v[1], v[2]])).to_versor_unchecked();
71        }
72    }
73}
74
75impl<B> LocalTrial<B> for Rotate<Versor>
76where
77    B: Orientation<Rotation = Versor>,
78{
79    /// Perturb a body's orientation by a random amount.
80    ///
81    /// In three dimensions, we design this perturbation as a versor whose
82    /// distribution is centered on the existing orientation and whose distribution is
83    /// narrow. To do so, we sample from a 3-dimensional Normal distribution
84    /// in the tangent space of SO(3), lift to the manifold, then rotate to center on
85    /// the current orientation. The result is a small displacement from
86    /// a quaternion input, with fast decay in the tails that make large displacements
87    /// unlikely. This is desirable for Monte Carlo, as large moves are very likely to
88    /// be rejected. This sampling obeys detailed balance.
89    ///
90    /// # Example
91    ///
92    /// ```
93    /// use hoomd_mc::{LocalTrial, Rotate};
94    /// use hoomd_microstate::property::OrientedPoint;
95    /// use hoomd_vector::{Cartesian, Versor};
96    /// use rand::{Rng, SeedableRng, rngs::StdRng};
97    ///
98    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
99    /// let mut rng = StdRng::seed_from_u64(1);
100    /// let initial_orientation = Versor::default();
101    /// let body_properties = OrientedPoint {
102    ///     position: Cartesian::from([0.0, 0.0, 0.0]),
103    ///     orientation: initial_orientation,
104    /// };
105    /// let rotate = Rotate::with_maximum_rotation(0.1.try_into()?);
106    ///
107    /// let new_body_properties = rotate.propose(&mut rng, body_properties);
108    /// assert!(new_body_properties.orientation != initial_orientation);
109    /// # Ok(())
110    /// # }
111    /// ```
112    #[inline]
113    fn propose<R: Rng>(&self, rng: &mut R, body_properties: B) -> B {
114        let mut trial = body_properties;
115        let displacement = VersorDisplacement {
116            std_dev: self.maximum_rotation.get(),
117        };
118
119        let delta_quat = displacement.sample(rng);
120        *trial.orientation_mut() = delta_quat.combine(trial.orientation());
121
122        trial
123    }
124}
125
126impl Adjust for Rotate<Versor> {
127    /// Change the maximum trial move size by the given scale factor.
128    #[inline]
129    fn adjust(&mut self, factor: PositiveReal) {
130        self.maximum_rotation *= factor;
131
132        if self.maximum_rotation.get() > PI / 2.0 {
133            self.maximum_rotation = (PI / 2.0)
134                .try_into()
135                .expect("PI/2.0 should be a positive real");
136        }
137    }
138}
139
140#[cfg(test)]
141mod tests {
142    use super::*;
143    use assert2::check;
144    use hoomd_microstate::property::OrientedPoint;
145    use hoomd_vector::{Cartesian, Versor};
146    use rand::{SeedableRng, rngs::StdRng};
147    use rstest::*;
148
149    /// Number of trial moves to test.
150    const N: usize = 65536;
151
152    #[rstest]
153    fn rotate(#[values(0.1, 1.0, 1.53)] a: f64) {
154        // Ensure that `Rotate` proposes moves that rotate the body with a valid
155        // range of maximum distances.
156
157        let mut delta_thetas = Vec::with_capacity(N);
158
159        let mut rng = StdRng::seed_from_u64(1);
160        let body = OrientedPoint {
161            position: Cartesian::from([0.0, 0.0, 0.0]),
162            orientation: Versor::default(),
163        };
164        let rotate = Rotate::with_maximum_rotation(
165            a.try_into()
166                .expect("hard-coded constant should be a positive real"),
167        );
168
169        for _ in 0..N {
170            let trial = rotate.propose(&mut rng, body);
171            let delta_theta = trial.orientation.arc_distance(&body.orientation);
172            delta_thetas.push(delta_theta);
173        }
174
175        // Validate our distribution's mean is relatively close to 0.75 * a.
176        // Why this is the correct choice? I'm not sure. But most values of `a` result
177        // in a distribution whose mean is near 0.75 * a
178        let mean = (delta_thetas.iter().sum::<f64>() / N as f64).abs();
179        assert!(
180            (mean - (0.75 * a)).abs() < (0.1 * a),
181            "{} !< {}",
182            (mean - (0.75 * a)).abs(),
183            (0.1 * a)
184        );
185    }
186
187    #[test]
188    fn test_adjust() -> anyhow::Result<()> {
189        let mut rotate = Rotate::<Versor>::with_maximum_rotation(0.5.try_into()?);
190
191        rotate.adjust(2.0.try_into()?);
192        check!(rotate.maximum_rotation().get() == 1.0);
193
194        rotate.adjust(0.5.try_into()?);
195        check!(rotate.maximum_rotation().get() == 0.5);
196
197        rotate.adjust(10.0.try_into()?);
198        check!(rotate.maximum_rotation().get() == PI / 2.0);
199
200        Ok(())
201    }
202}
hoomd_mc/rotate/versor.rs

hoomd_mc/rotate/
versor.rs
