hoomd_mc/
sweep.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 Sweep
5
6use rand::RngExt;
7use serde::{Deserialize, Serialize};
8use std::fmt::Display;
9
10use hoomd_interaction::DeltaEnergyOne;
11use hoomd_microstate::{
12    Body, Microstate, SiteKey, Transform,
13    boundary::{GenerateGhosts, Wrap},
14    property::Position,
15};
16use hoomd_simulation::macrostate::Temperature;
17use hoomd_spatial::PointUpdate;
18use hoomd_utility::valid::OpenUnitIntervalNumber;
19
20use super::{Adjust, Count, LocalTrial, Trial, Tune, tune_local::tune_local_trial};
21
22/// Apply a local trial move to every body in the microstate.
23///
24/// The first field in the tuple struct determines what trial moves `Sweep`
25/// attempts.
26///
27/// # Example
28///
29/// ```
30/// use hoomd_mc::{Sweep, Translate};
31/// use hoomd_vector::Cartesian;
32///
33/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
34/// let d = 0.1;
35/// let translate =
36///     Translate::<Cartesian<2>>::with_maximum_distance(d.try_into()?);
37/// let translate_sweep = Sweep(translate);
38/// # Ok(())
39/// # }
40/// ```
41#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
42pub struct Sweep<L>(pub L);
43
44impl<P, B, S, X, C, L, H, MA> Trial<Microstate<B, S, X, C>, H, MA> for Sweep<L>
45where
46    P: Copy,
47    B: Copy + Default + Transform<S> + Position<Position = P>,
48    S: Copy + Default + Position<Position = P>,
49    X: PointUpdate<P, SiteKey>,
50    L: LocalTrial<B>,
51    H: DeltaEnergyOne<B, S, X, C>,
52    C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
53    MA: Temperature,
54{
55    type Count = Count;
56
57    /// Apply a local trial move to each body in the microstate.
58    ///
59    /// Each trial move is accepted when:
60    /// ```math
61    /// r < \exp\left(\frac{-\Delta H}{kT}\right)
62    /// ```
63    /// where `r` is a random value uniformly distributed in `[0,1)`, $`\Delta H`$ is
64    /// the change in energy computed by the given `hamiltonian` and $`kT`$ is the
65    /// `temperature` given in `macrostate`.
66    ///
67    /// # Example
68    ///
69    /// ```
70    /// use hoomd_interaction::Zero;
71    /// use hoomd_mc::{Sweep, Translate, Trial};
72    /// use hoomd_microstate::{Body, Microstate, property::Position};
73    /// use hoomd_simulation::macrostate::Isothermal;
74    /// use hoomd_vector::Cartesian;
75    ///
76    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
77    /// let mut microstate = Microstate::new();
78    /// microstate.add_body(Body::point(Cartesian::from([0.0, 0.0])));
79    /// let d = 0.1;
80    /// let translate = Translate::with_maximum_distance(d.try_into()?);
81    /// let mut translate_sweep = Sweep(translate);
82    ///
83    /// let hamiltonian = Zero;
84    /// let macrostate = Isothermal { temperature: 1.0 };
85    ///
86    /// for _ in 0..1_000 {
87    ///     translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
88    ///     microstate.increment_step();
89    /// }
90    /// # Ok(())
91    /// # }
92    /// ```
93    #[inline]
94    fn apply(
95        &mut self,
96        microstate: &mut Microstate<B, S, X, C>,
97        hamiltonian: &H,
98        macrostate: &MA,
99    ) -> Self::Count {
100        let kt = macrostate.temperature();
101        let mut rng = microstate.counter().make_rng();
102        let mut count = Self::Count::default();
103        let mut trial = Body::<B, S>::default();
104
105        // For loop over a range instead of bodies().iter() as the latter holds an immutable borrow.
106        // The call to `update_body_properties` makes a mutable borrow of microstate.
107        for body_index in 0..microstate.bodies().len() {
108            trial.clone_from(&microstate.bodies()[body_index].item);
109
110            // Wrap the body position here. The site positions will be wrapped
111            // by the delta_energy methods and again by update_body_properties.
112            // We could early reject if we checked the site properties first. If
113            // performance becomes an issue, we could also wrap site properties
114            // once here and pass those to unchecked variants of the delta
115            // energy and update methods.
116            match microstate
117                .boundary()
118                .wrap(self.0.propose(&mut rng, trial.properties))
119            {
120                Ok(new_properties) => {
121                    trial.properties = new_properties;
122
123                    let delta_h = hamiltonian.delta_energy_one(microstate, body_index, &trial);
124                    if delta_h != f64::INFINITY
125                        && (delta_h <= 0.0 || rng.random::<f64>() < (-delta_h / kt).exp())
126                        && microstate
127                            .update_body_properties(body_index, trial.properties)
128                            .is_ok()
129                    {
130                        count.accepted += 1;
131                    } else {
132                        count.rejected += 1;
133                    }
134                }
135                Err(_) => count.rejected += 1,
136            }
137        }
138
139        microstate.increment_substep();
140        count
141    }
142}
143
144impl<P, B, S, X, C, L, H, MA> Tune<P, B, S, X, C, L, H, MA> for Sweep<L>
145where
146    P: Copy,
147    B: Copy + Default + Transform<S> + Position<Position = P>,
148    S: Copy + Default + Position<Position = P>,
149    X: PointUpdate<P, SiteKey>,
150    L: LocalTrial<B> + Adjust + Display,
151    H: DeltaEnergyOne<B, S, X, C>,
152    C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
153    MA: Temperature,
154{
155    /// Tune the trial move maximum size to achieve a given acceptance ratio.
156    ///
157    /// Use [`tune_default`] unless you have a specific need to adjust the
158    /// tuning parameters.
159    ///
160    /// [`tune_default`]: Self::tune_default
161    ///
162    /// # Example
163    ///
164    /// ```
165    /// use hoomd_geometry::shape::Rectangle;
166    /// use hoomd_interaction::{
167    ///     MaximumInteractionRange, PairwiseCutoff, pairwise::HardSphere,
168    /// };
169    /// use hoomd_mc::{Sweep, Translate, Trial, Tune};
170    /// use hoomd_microstate::{
171    ///     Body, Microstate, boundary::Periodic, property::Position,
172    /// };
173    /// use hoomd_simulation::macrostate::Isothermal;
174    /// use hoomd_vector::Cartesian;
175    ///
176    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
177    /// let square = Rectangle::with_equal_edges(2.2.try_into()?);
178    /// let mut microstate = Microstate::builder()
179    ///     .boundary(Periodic::new(1.0, square)?)
180    ///     .try_build()?;
181    /// microstate.add_body(Body::point(Cartesian::from([-0.6, -0.6])))?;
182    /// microstate.add_body(Body::point(Cartesian::from([-0.6, 0.6])))?;
183    /// microstate.add_body(Body::point(Cartesian::from([0.6, -0.6])))?;
184    /// microstate.add_body(Body::point(Cartesian::from([0.6, 0.6])))?;
185    /// let d = 0.1;
186    /// let translate = Translate::with_maximum_distance(d.try_into()?);
187    /// let mut translate_sweep = Sweep(translate);
188    ///
189    /// let hamiltonian = PairwiseCutoff(HardSphere { diameter: 1.0 });
190    /// let macrostate = Isothermal { temperature: 1.0 };
191    ///
192    /// translate_sweep.tune_default(&microstate, &hamiltonian, &macrostate);
193    ///
194    /// # Ok(())
195    /// # }
196    /// ```
197    #[inline]
198    fn tune(
199        &mut self,
200        microstate: &Microstate<B, S, X, C>,
201        hamiltonian: &H,
202        macrostate: &MA,
203        target_acceptance: OpenUnitIntervalNumber,
204        samples: usize,
205        steps: usize,
206    ) {
207        tune_local_trial(
208            &mut self.0,
209            microstate,
210            hamiltonian,
211            macrostate,
212            target_acceptance,
213            samples,
214            steps,
215        );
216    }
217}
218
219#[cfg(test)]
220mod tests {
221    use super::*;
222    use crate::Translate;
223    use approxim::assert_relative_eq;
224    use hoomd_geometry::shape::Hypercuboid;
225    use hoomd_interaction::{External, SiteEnergy, TotalEnergy, Zero};
226    use hoomd_microstate::{boundary::Closed, property::Point};
227    use hoomd_simulation::macrostate::Isothermal;
228    use hoomd_vector::{Cartesian, InnerProduct};
229    use rand::Rng;
230    use rstest::*;
231
232    const K: f64 = 2.0;
233    const N_STEPS: u64 = 200_000;
234
235    struct Harmonic(Cartesian<2>);
236
237    impl SiteEnergy<Point<Cartesian<2>>> for Harmonic {
238        fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
239            1.0 / 2.0 * K * (site_properties.position - self.0).norm_squared()
240        }
241    }
242
243    struct Right;
244    impl LocalTrial<Point<Cartesian<2>>> for Right {
245        fn propose<R: Rng>(
246            &self,
247            _rng: &mut R,
248            body_properties: Point<Cartesian<2>>,
249        ) -> Point<Cartesian<2>> {
250            let mut trial = body_properties;
251            trial.position_mut()[0] += 1.0;
252            trial
253        }
254    }
255
256    #[rstest]
257    fn harmonic_oscillators(#[values(1.0, 2.5)] kt: f64) {
258        // Model a harmonic oscillator and validate the average position and energy distribution.
259        // Check with a relatively large tolerance because N_STEPS is relatively
260        // small to keep the test run short.
261        const EPSILON: f64 = 0.3;
262
263        let origin = Cartesian::from([1.0, -2.0]);
264
265        let mut microstate = Microstate::new();
266        microstate
267            .add_body(Body::point(origin))
268            .expect("the hard-coded body should be inside the boundary");
269        let hamiltonian = External(Harmonic(origin));
270        let macrostate = Isothermal { temperature: kt };
271
272        let d = 0.1;
273        let translate = Translate::with_maximum_distance(
274            d.try_into()
275                .expect("hard-coded constant should be positive"),
276        );
277        let mut translate_sweep = Sweep(translate);
278
279        let mut position_accumulator = Cartesian::default();
280        let mut energy_accumulator = 0.0;
281
282        for _ in 0..N_STEPS {
283            translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
284
285            position_accumulator += microstate.bodies()[0].item.properties.position;
286            energy_accumulator += hamiltonian.total_energy(&microstate);
287
288            microstate.increment_step();
289        }
290
291        let position_average = position_accumulator / (N_STEPS as f64);
292        assert_relative_eq!(position_average[0], origin[0], epsilon = EPSILON);
293        assert_relative_eq!(position_average[1], origin[1], epsilon = EPSILON);
294
295        let energy_average = energy_accumulator / (N_STEPS as f64);
296        assert_relative_eq!(energy_average, kt, epsilon = EPSILON);
297    }
298
299    #[test]
300    fn reject_boundary_body() {
301        let cuboid = Hypercuboid {
302            edge_lengths: [
303                4.0.try_into()
304                    .expect("hard-coded constant should be positive"),
305                4.0.try_into()
306                    .expect("hard-coded constant should be positive"),
307            ],
308        };
309        let square = Closed(cuboid);
310
311        let mut microstate = Microstate::builder()
312            .boundary(square)
313            .bodies([Body::point([0.0, 0.0].into())])
314            .try_build()
315            .expect("the hard-coded bodies should be in the boundary");
316        let hamiltonian = Zero;
317        let translate = Right;
318        let mut translate_sweep = Sweep(translate);
319        let macrostate = Isothermal { temperature: 1.0 };
320
321        // The first move to the right ends in the boundary and should be accepted.
322        let counter = translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
323        assert_eq!(counter.accepted, 1);
324        assert_eq!(counter.rejected, 0);
325
326        // The second move to the right places the body just on the boundary and should be
327        // rejected.
328        let counter = translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
329        assert_eq!(counter.accepted, 0);
330        assert_eq!(counter.rejected, 1);
331    }
332
333    #[test]
334    fn reject_boundary_site() {
335        let body = Body {
336            properties: Point::new(Cartesian::from([0.0, 0.0])),
337            sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
338        };
339
340        let cuboid = Hypercuboid {
341            edge_lengths: [
342                6.0.try_into()
343                    .expect("hard-coded constant should be positive"),
344                6.0.try_into()
345                    .expect("hard-coded constant should be positive"),
346            ],
347        };
348        let square = Closed(cuboid);
349        let mut microstate = Microstate::builder()
350            .boundary(square)
351            .bodies([body])
352            .try_build()
353            .expect("the hard-coded bodies should be in the boundary");
354        let hamiltonian = Zero;
355        let translate = Right;
356        let mut translate_sweep = Sweep(translate);
357        let macrostate = Isothermal { temperature: 1.0 };
358
359        // The first move to the right ends in the boundary and should be accepted.
360        let counter = translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
361        assert_eq!(counter.accepted, 1);
362        assert_eq!(counter.rejected, 0);
363
364        // The second move to the right places the body just on the boundary and should be
365        // rejected.
366        let counter = translate_sweep.apply(&mut microstate, &hamiltonian, &macrostate);
367        assert_eq!(counter.accepted, 0);
368        assert_eq!(counter.rejected, 1);
369    }
370}
hoomd_mc/sweep.rs

hoomd_mc/
sweep.rs
