hoomd_mc/rotate/
angle.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 Angle
5
6use rand::{
7    Rng,
8    distr::{Distribution, Uniform},
9};
10use std::f64::consts::PI;
11
12use hoomd_microstate::property::Orientation;
13use hoomd_utility::valid::PositiveReal;
14use hoomd_vector::Angle;
15
16use crate::{Adjust, LocalTrial, Rotate};
17
18impl<B> LocalTrial<B> for Rotate<Angle>
19where
20    B: Orientation<Rotation = Angle>,
21{
22    /// Perturb a body's orientation by a random amount.
23    ///
24    /// # Example
25    ///
26    /// ```
27    /// use hoomd_mc::{LocalTrial, Rotate};
28    /// use hoomd_microstate::property::OrientedPoint;
29    /// use hoomd_vector::{Angle, Cartesian};
30    /// use rand::{Rng, SeedableRng, rngs::StdRng};
31    ///
32    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
33    /// let mut rng = StdRng::seed_from_u64(1);
34    /// let body_properties = OrientedPoint {
35    ///     position: Cartesian::from([0.0, 0.0]),
36    ///     orientation: Angle::from(0.0),
37    /// };
38    /// let rotate = Rotate::with_maximum_rotation(0.1.try_into()?);
39    ///
40    /// let new_body_properties = rotate.propose(&mut rng, body_properties);
41    /// assert!(new_body_properties.orientation.theta.abs() < 0.1);
42    /// # Ok(())
43    /// # }
44    /// ```
45    #[inline]
46    fn propose<R: Rng>(&self, rng: &mut R, body_properties: B) -> B {
47        let mut trial = body_properties;
48
49        let a = self.maximum_rotation.get();
50        let uniform = Uniform::new(-a, a).expect("a should be positive");
51
52        let delta_theta = uniform.sample(rng);
53        trial.orientation_mut().theta += delta_theta;
54
55        trial
56    }
57}
58
59impl Adjust for Rotate<Angle> {
60    /// Change the maximum trial move size by the given scale factor.
61    #[inline]
62    fn adjust(&mut self, factor: PositiveReal) {
63        self.maximum_rotation *= factor;
64
65        if self.maximum_rotation.get() > PI {
66            self.maximum_rotation = (PI).try_into().expect("PI should be a positive real");
67        }
68    }
69}
70
71#[cfg(test)]
72mod tests {
73    use super::*;
74    use assert2::check;
75    use hoomd_microstate::property::OrientedPoint;
76    use hoomd_vector::{Angle, Cartesian};
77    use rand::{SeedableRng, rngs::StdRng};
78    use rstest::*;
79
80    /// Number of trial moves to test.
81    const N: usize = 1024;
82
83    #[rstest]
84    fn rotate(#[values(0.1, 1.0)] a: f64) {
85        // Ensure that `Rotate` proposes moves that rotate the body with a valid
86        // range of maximum distances.
87
88        let mut total = 0.0;
89        let mut min_norm = f64::INFINITY;
90        let mut max_norm = 0.0_f64;
91
92        let mut rng = StdRng::seed_from_u64(1);
93        let body = OrientedPoint {
94            position: Cartesian::from([0.0, 0.0]),
95            orientation: Angle::default(),
96        };
97        let rotate = Rotate::with_maximum_rotation(
98            a.try_into()
99                .expect("hard-coded constant should be a positive real"),
100        );
101
102        for _ in 0..N {
103            let trial = rotate.propose(&mut rng, body);
104
105            let delta_theta = trial.orientation.theta - body.orientation.theta;
106            total += delta_theta;
107            min_norm = min_norm.min(delta_theta.abs());
108            max_norm = max_norm.max(delta_theta.abs());
109        }
110
111        assert!(max_norm <= a);
112
113        let average = total / N as f64;
114
115        // Validate with appropriately loose tolerances to account for the small sample size.
116        assert!(min_norm < a * 0.1);
117        assert!(max_norm > a * 0.9);
118        assert!(average.abs() < a * 0.1);
119    }
120
121    #[test]
122    fn test_adjust() -> anyhow::Result<()> {
123        let mut rotate = Rotate::<Angle>::with_maximum_rotation(0.5.try_into()?);
124
125        rotate.adjust(2.0.try_into()?);
126        check!(rotate.maximum_rotation().get() == 1.0);
127
128        rotate.adjust(0.5.try_into()?);
129        check!(rotate.maximum_rotation().get() == 0.5);
130
131        rotate.adjust(10.0.try_into()?);
132        check!(rotate.maximum_rotation().get() == PI);
133
134        Ok(())
135    }
136}