1use 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 #[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 #[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 const N: usize = 1024;
82
83 #[rstest]
84 fn rotate(#[values(0.1, 1.0)] a: f64) {
85 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 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}