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, Tagged, 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, TuneOptions, tune_local::tune_local_trial};
21
22/// Apply a local trial move to bodies in the microstate.
23///
24/// The first field in the tuple struct determines what trial moves `Sweep`
25/// attempts. [`Sweep::apply`] applies the trial move to each body in the
26/// microstate once. [`Sweep::apply_with_filter`] applies the trial move
27/// to select bodies.
28///
29/// # Example
30///
31/// ```
32/// use hoomd_mc::{Sweep, Translate};
33/// use hoomd_vector::Cartesian;
34///
35/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
36/// let d = 0.1;
37/// let translate =
38/// Translate::<Cartesian<2>>::with_maximum_distance(d.try_into()?);
39/// let translate_sweep = Sweep(translate);
40/// # Ok(())
41/// # }
42/// ```
43#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
44pub struct Sweep<L>(pub L);
45
46impl<P, B, S, X, C, L, H, MA> Trial<Microstate<B, S, X, C>, H, MA> for Sweep<L>
47where
48 P: Copy,
49 B: Copy + Default + Transform<S> + Position<Position = P>,
50 S: Copy + Default + Position<Position = P>,
51 X: PointUpdate<P, SiteKey>,
52 L: LocalTrial<B>,
53 H: DeltaEnergyOne<B, S, X, C>,
54 C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
55 MA: Temperature,
56{
57 type Count = Count;
58
59 /// Apply a local trial move to each body in the microstate.
60 ///
61 /// Each trial move is accepted when:
62 /// ```math
63 /// r < \exp\left(\frac{-\Delta H}{kT}\right)
64 /// ```
65 /// where `r` is a random value uniformly distributed in `[0,1)`, $`\Delta H`$ is
66 /// the change in energy computed by the given `hamiltonian` and $`kT`$ is the
67 /// `temperature` given in `macrostate`.
68 ///
69 /// # Example
70 ///
71 /// ```
72 /// use hoomd_interaction::Zero;
73 /// use hoomd_mc::{Sweep, Translate, Trial};
74 /// use hoomd_microstate::{Body, Microstate, property::Position};
75 /// use hoomd_simulation::macrostate::Isothermal;
76 /// use hoomd_vector::Cartesian;
77 ///
78 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
79 /// let mut microstate = Microstate::new();
80 /// microstate.add_body(Body::point(Cartesian::from([0.0, 0.0])));
81 /// let d = 0.1;
82 /// let translate = Translate::with_maximum_distance(d.try_into()?);
83 /// let mut translate_sweep = Sweep(translate);
84 ///
85 /// let hamiltonian = Zero;
86 /// let macrostate = Isothermal { temperature: 1.0 };
87 ///
88 /// for _ in 0..1_000 {
89 /// translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
90 /// microstate.increment_step();
91 /// }
92 /// # Ok(())
93 /// # }
94 /// ```
95 #[inline]
96 fn apply(
97 &mut self,
98 microstate: &mut Microstate<B, S, X, C>,
99 hamiltonian: &H,
100 macrostate: &MA,
101 ) -> Self::Count {
102 self.apply_with_filter(microstate, hamiltonian, macrostate, |_| true)
103 }
104}
105
106impl<L> Sweep<L> {
107 /// Apply a local trial move to select bodies in the microstate.
108 ///
109 /// `apply_with_filter` applies trial moves to bodies where `should_move_body`
110 /// returns `true`.
111 ///
112 /// Each trial move is accepted when:
113 /// ```math
114 /// r < \exp\left(\frac{-\Delta H}{kT}\right)
115 /// ```
116 /// where `r` is a random value uniformly distributed in `[0,1)`, $`\Delta H`$ is
117 /// the change in energy computed by the given `hamiltonian` and $`kT`$ is the
118 /// `temperature` given in `macrostate`.
119 ///
120 /// # Example
121 ///
122 /// ```
123 /// use hoomd_interaction::Zero;
124 /// use hoomd_mc::{Sweep, Translate, Trial};
125 /// use hoomd_microstate::{Body, Microstate, property::Position};
126 /// use hoomd_simulation::macrostate::Isothermal;
127 /// use hoomd_vector::Cartesian;
128 ///
129 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
130 /// let mut microstate = Microstate::new();
131 /// microstate.add_body(Body::point(Cartesian::from([0.0, 0.0])));
132 /// microstate.add_body(Body::point(Cartesian::from([1.0, 0.0])));
133 /// microstate.add_body(Body::point(Cartesian::from([0.0, 1.0])));
134 /// let d = 0.1;
135 /// let translate = Translate::with_maximum_distance(d.try_into()?);
136 /// let mut translate_sweep = Sweep(translate);
137 ///
138 /// let hamiltonian = Zero;
139 /// let macrostate = Isothermal { temperature: 1.0 };
140 ///
141 /// for _ in 0..1_000 {
142 /// translate_sweep.apply_with_filter(
143 /// &mut microstate,
144 /// &hamiltonian,
145 /// ¯ostate,
146 /// |b| b.tag != 0,
147 /// );
148 /// microstate.increment_step();
149 /// }
150 ///
151 /// assert_eq!(
152 /// microstate.bodies()[0].item.properties.position,
153 /// [0.0, 0.0].into()
154 /// );
155 /// # Ok(())
156 /// # }
157 /// ```
158 #[inline]
159 pub fn apply_with_filter<P, B, S, X, C, H, MA, F>(
160 &self,
161 microstate: &mut Microstate<B, S, X, C>,
162 hamiltonian: &H,
163 macrostate: &MA,
164 should_move_body: F,
165 ) -> Count
166 where
167 P: Copy,
168 B: Copy + Default + Transform<S> + Position<Position = P>,
169 S: Copy + Default + Position<Position = P>,
170 X: PointUpdate<P, SiteKey>,
171 L: LocalTrial<B>,
172 H: DeltaEnergyOne<B, S, X, C>,
173 C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
174 MA: Temperature,
175 F: Fn(&Tagged<Body<B, S>>) -> bool,
176 {
177 let kt = macrostate.temperature();
178 let mut rng = microstate.counter().make_rng();
179 let mut count = Count::default();
180 let mut trial = Body::<B, S>::default();
181
182 // For loop over a range instead of bodies().iter() as the latter holds an immutable borrow.
183 // The call to `update_body_properties` makes a mutable borrow of microstate.
184 for body_index in 0..microstate.bodies().len() {
185 let body = µstate.bodies()[body_index];
186 if !should_move_body(body) {
187 continue;
188 }
189
190 trial.clone_from(&body.item);
191
192 // Wrap the body position here. The site positions will be wrapped
193 // by the delta_energy methods and again by update_body_properties.
194 // We could early reject if we checked the site properties first. If
195 // performance becomes an issue, we could also wrap site properties
196 // once here and pass those to unchecked variants of the delta
197 // energy and update methods.
198 match microstate
199 .boundary()
200 .wrap(self.0.propose(&mut rng, trial.properties))
201 {
202 Ok(new_properties) => {
203 trial.properties = new_properties;
204
205 let delta_h = hamiltonian.delta_energy_one(microstate, body_index, &trial);
206 if delta_h != f64::INFINITY
207 && (delta_h <= 0.0 || rng.random::<f64>() < (-delta_h / kt).exp())
208 && microstate
209 .update_body_properties(body_index, trial.properties)
210 .is_ok()
211 {
212 count.accepted += 1;
213 } else {
214 count.rejected += 1;
215 }
216 }
217 Err(_) => count.rejected += 1,
218 }
219 }
220
221 microstate.increment_substep();
222 count
223 }
224
225 /// Tune the trial move maximum size to achieve a given acceptance ratio.
226 ///
227 /// `tune_with_options_and_filter` applies trial moves to bodies where
228 /// `should_move_body` returns `true`.
229 ///
230 /// # Example
231 ///
232 /// ```
233 /// use hoomd_geometry::shape::Rectangle;
234 /// use hoomd_interaction::{
235 /// MaximumInteractionRange, PairwiseCutoff, pairwise::HardSphere,
236 /// };
237 /// use hoomd_mc::{Sweep, Translate, Trial, Tune, TuneOptions};
238 /// use hoomd_microstate::{
239 /// Body, Microstate, boundary::Periodic, property::Position,
240 /// };
241 /// use hoomd_simulation::macrostate::Isothermal;
242 /// use hoomd_vector::Cartesian;
243 ///
244 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
245 /// let square = Rectangle::with_equal_edges(2.2.try_into()?);
246 /// let mut microstate = Microstate::builder()
247 /// .boundary(Periodic::new(1.0, square)?)
248 /// .try_build()?;
249 /// microstate.add_body(Body::point(Cartesian::from([-0.6, -0.6])))?;
250 /// microstate.add_body(Body::point(Cartesian::from([-0.6, 0.6])))?;
251 /// microstate.add_body(Body::point(Cartesian::from([0.6, -0.6])))?;
252 /// microstate.add_body(Body::point(Cartesian::from([0.6, 0.6])))?;
253 /// let d = 0.1;
254 /// let translate = Translate::with_maximum_distance(d.try_into()?);
255 /// let mut translate_sweep = Sweep(translate);
256 ///
257 /// let hamiltonian = PairwiseCutoff(HardSphere { diameter: 1.0 });
258 /// let macrostate = Isothermal { temperature: 1.0 };
259 ///
260 /// translate_sweep.tune_with_options_and_filter(
261 /// µstate,
262 /// &hamiltonian,
263 /// ¯ostate,
264 /// &TuneOptions::default(),
265 /// |b| b.tag >= 2,
266 /// );
267 ///
268 /// # Ok(())
269 /// # }
270 /// ```
271 #[inline]
272 pub fn tune_with_options_and_filter<P, B, S, X, C, H, MA, F>(
273 &mut self,
274 microstate: &Microstate<B, S, X, C>,
275 hamiltonian: &H,
276 macrostate: &MA,
277 options: &TuneOptions,
278 should_move_body: F,
279 ) where
280 P: Copy,
281 B: Copy + Default + Transform<S> + Position<Position = P>,
282 S: Copy + Default + Position<Position = P>,
283 X: PointUpdate<P, SiteKey>,
284 L: LocalTrial<B> + Adjust + Display,
285 H: DeltaEnergyOne<B, S, X, C>,
286 C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
287 MA: Temperature,
288 F: Fn(&Tagged<Body<B, S>>) -> bool,
289 {
290 tune_local_trial(
291 &mut self.0,
292 microstate,
293 hamiltonian,
294 macrostate,
295 options,
296 should_move_body,
297 );
298 }
299}
300
301impl<P, B, S, X, C, L, H, MA> Tune<P, B, S, X, C, L, H, MA> for Sweep<L>
302where
303 P: Copy,
304 B: Copy + Default + Transform<S> + Position<Position = P>,
305 S: Copy + Default + Position<Position = P>,
306 X: PointUpdate<P, SiteKey>,
307 L: LocalTrial<B> + Adjust + Display,
308 H: DeltaEnergyOne<B, S, X, C>,
309 C: Wrap<B> + Wrap<S> + GenerateGhosts<S>,
310 MA: Temperature,
311{
312 /// Tune the trial move maximum size to achieve a given acceptance ratio.
313 ///
314 /// # Example
315 ///
316 /// ```
317 /// use hoomd_geometry::shape::Rectangle;
318 /// use hoomd_interaction::{
319 /// MaximumInteractionRange, PairwiseCutoff, pairwise::HardSphere,
320 /// };
321 /// use hoomd_mc::{Sweep, Translate, Trial, Tune, TuneOptions};
322 /// use hoomd_microstate::{
323 /// Body, Microstate, boundary::Periodic, property::Position,
324 /// };
325 /// use hoomd_simulation::macrostate::Isothermal;
326 /// use hoomd_vector::Cartesian;
327 ///
328 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
329 /// let square = Rectangle::with_equal_edges(2.2.try_into()?);
330 /// let mut microstate = Microstate::builder()
331 /// .boundary(Periodic::new(1.0, square)?)
332 /// .try_build()?;
333 /// microstate.add_body(Body::point(Cartesian::from([-0.6, -0.6])))?;
334 /// microstate.add_body(Body::point(Cartesian::from([-0.6, 0.6])))?;
335 /// microstate.add_body(Body::point(Cartesian::from([0.6, -0.6])))?;
336 /// microstate.add_body(Body::point(Cartesian::from([0.6, 0.6])))?;
337 /// let d = 0.1;
338 /// let translate = Translate::with_maximum_distance(d.try_into()?);
339 /// let mut translate_sweep = Sweep(translate);
340 ///
341 /// let hamiltonian = PairwiseCutoff(HardSphere { diameter: 1.0 });
342 /// let macrostate = Isothermal { temperature: 1.0 };
343 ///
344 /// translate_sweep.tune_with_options(
345 /// µstate,
346 /// &hamiltonian,
347 /// ¯ostate,
348 /// &TuneOptions::default(),
349 /// );
350 ///
351 /// # Ok(())
352 /// # }
353 /// ```
354 #[inline]
355 fn tune(
356 &mut self,
357 microstate: &Microstate<B, S, X, C>,
358 hamiltonian: &H,
359 macrostate: &MA,
360 target_acceptance: OpenUnitIntervalNumber,
361 samples: usize,
362 steps: usize,
363 ) {
364 tune_local_trial(
365 &mut self.0,
366 microstate,
367 hamiltonian,
368 macrostate,
369 &TuneOptions {
370 target_acceptance,
371 samples,
372 steps,
373 },
374 |_| true,
375 );
376 }
377}
378
379#[cfg(test)]
380mod tests {
381 use super::*;
382 use crate::Translate;
383 use approxim::assert_relative_eq;
384 use assert2::check;
385 use hoomd_geometry::shape::Hypercuboid;
386 use hoomd_interaction::{External, SiteEnergy, TotalEnergy, Zero};
387 use hoomd_microstate::{boundary::Closed, property::Point};
388 use hoomd_simulation::macrostate::Isothermal;
389 use hoomd_vector::{Cartesian, InnerProduct};
390 use rand::Rng;
391 use rstest::*;
392
393 const K: f64 = 2.0;
394 const N_STEPS: u64 = 200_000;
395
396 struct Harmonic(Cartesian<2>);
397
398 impl SiteEnergy<Point<Cartesian<2>>> for Harmonic {
399 fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
400 1.0 / 2.0 * K * (site_properties.position - self.0).norm_squared()
401 }
402 }
403
404 struct Right;
405 impl LocalTrial<Point<Cartesian<2>>> for Right {
406 fn propose<R: Rng>(
407 &self,
408 _rng: &mut R,
409 body_properties: Point<Cartesian<2>>,
410 ) -> Point<Cartesian<2>> {
411 let mut trial = body_properties;
412 trial.position_mut()[0] += 1.0;
413 trial
414 }
415 }
416
417 #[rstest]
418 fn harmonic_oscillators(#[values(1.0, 2.5)] kt: f64) {
419 // Model a harmonic oscillator and validate the average position and energy distribution.
420 // Check with a relatively large tolerance because N_STEPS is relatively
421 // small to keep the test run short.
422 const EPSILON: f64 = 0.3;
423
424 let origin = Cartesian::from([1.0, -2.0]);
425
426 let mut microstate = Microstate::new();
427 microstate
428 .add_body(Body::point(origin))
429 .expect("the hard-coded body should be inside the boundary");
430 let hamiltonian = External(Harmonic(origin));
431 let macrostate = Isothermal { temperature: kt };
432
433 let d = 0.1;
434 let translate = Translate::with_maximum_distance(
435 d.try_into()
436 .expect("hard-coded constant should be positive"),
437 );
438 let mut translate_sweep = Sweep(translate);
439
440 let mut position_accumulator = Cartesian::default();
441 let mut energy_accumulator = 0.0;
442
443 for _ in 0..N_STEPS {
444 translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
445
446 position_accumulator += microstate.bodies()[0].item.properties.position;
447 energy_accumulator += hamiltonian.total_energy(µstate);
448
449 microstate.increment_step();
450 }
451
452 let position_average = position_accumulator / (N_STEPS as f64);
453 assert_relative_eq!(position_average[0], origin[0], epsilon = EPSILON);
454 assert_relative_eq!(position_average[1], origin[1], epsilon = EPSILON);
455
456 let energy_average = energy_accumulator / (N_STEPS as f64);
457 assert_relative_eq!(energy_average, kt, epsilon = EPSILON);
458 }
459
460 #[test]
461 fn reject_boundary_body() {
462 let cuboid = Hypercuboid {
463 edge_lengths: [
464 4.0.try_into()
465 .expect("hard-coded constant should be positive"),
466 4.0.try_into()
467 .expect("hard-coded constant should be positive"),
468 ],
469 };
470 let square = Closed(cuboid);
471
472 let mut microstate = Microstate::builder()
473 .boundary(square)
474 .bodies([Body::point([0.0, 0.0].into())])
475 .try_build()
476 .expect("the hard-coded bodies should be in the boundary");
477 let hamiltonian = Zero;
478 let translate = Right;
479 let mut translate_sweep = Sweep(translate);
480 let macrostate = Isothermal { temperature: 1.0 };
481
482 // The first move to the right ends in the boundary and should be accepted.
483 let counter = translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
484 assert_eq!(counter.accepted, 1);
485 assert_eq!(counter.rejected, 0);
486
487 // The second move to the right places the body just on the boundary and should be
488 // rejected.
489 let counter = translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
490 assert_eq!(counter.accepted, 0);
491 assert_eq!(counter.rejected, 1);
492 }
493
494 #[test]
495 fn reject_boundary_site() {
496 let body = Body {
497 properties: Point::new(Cartesian::from([0.0, 0.0])),
498 sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
499 };
500
501 let cuboid = Hypercuboid {
502 edge_lengths: [
503 6.0.try_into()
504 .expect("hard-coded constant should be positive"),
505 6.0.try_into()
506 .expect("hard-coded constant should be positive"),
507 ],
508 };
509 let square = Closed(cuboid);
510 let mut microstate = Microstate::builder()
511 .boundary(square)
512 .bodies([body])
513 .try_build()
514 .expect("the hard-coded bodies should be in the boundary");
515 let hamiltonian = Zero;
516 let translate = Right;
517 let mut translate_sweep = Sweep(translate);
518 let macrostate = Isothermal { temperature: 1.0 };
519
520 // The first move to the right ends in the boundary and should be accepted.
521 let counter = translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
522 assert_eq!(counter.accepted, 1);
523 assert_eq!(counter.rejected, 0);
524
525 // The second move to the right places the body just on the boundary and should be
526 // rejected.
527 let counter = translate_sweep.apply(&mut microstate, &hamiltonian, ¯ostate);
528 assert_eq!(counter.accepted, 0);
529 assert_eq!(counter.rejected, 1);
530 }
531
532 #[test]
533 fn filter() -> anyhow::Result<()> {
534 let cuboid = Hypercuboid::with_equal_edges(4.0.try_into()?);
535 let square = Closed(cuboid);
536
537 let mut microstate = Microstate::builder()
538 .boundary(square)
539 .bodies([
540 Body::point([0.0, 0.0].into()),
541 Body::point([0.0, 0.0].into()),
542 Body::point([0.0, 0.0].into()),
543 ])
544 .try_build()?;
545 let hamiltonian = Zero;
546 let translate = Translate::with_maximum_distance(0.1.try_into()?);
547 let translate_sweep = Sweep(translate);
548 let macrostate = Isothermal { temperature: 1.0 };
549
550 let counter =
551 translate_sweep.apply_with_filter(&mut microstate, &hamiltonian, ¯ostate, |body| {
552 body.tag != 0
553 });
554 check!(counter.accepted == 2);
555 check!(counter.rejected == 0);
556 check!(microstate.bodies()[0].item.properties.position == [0.0, 0.0].into());
557
558 Ok(())
559 }
560}