hoomd_interaction/
pairwise_cutoff.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 `PairwiseCutoff`
5
6use serde::{Deserialize, Serialize};
7
8use hoomd_microstate::{
9    Body, Microstate, Site, SiteKey, Transform, boundary::Wrap, property::Position,
10};
11use hoomd_spatial::PointsNearBall;
12
13use crate::{
14    DeltaEnergyInsert, DeltaEnergyOne, DeltaEnergyRemove, MaximumInteractionRange, SitePairEnergy,
15    TotalEnergy,
16};
17
18/// Short-ranged pairwise interactions between sites.
19///
20/// A [`PairwiseCutoff`] newtype wrapped around a type that implements
21/// [`SitePairEnergy`] represents:
22///
23/// ```math
24/// U_\mathrm{total} = \sum_{i=0}^{N-1}\sum_{j=i+1}^{N-1} U\left(s_i, s_j \right) \left[ \left|\vec{r}_j - \vec{r}_i\right| \lt r_\mathrm{cut} \right]\left[b_i \ne b_j\right]
25/// ```
26/// where $`U(s_i, s_j)`$ is the potential computed by [`SitePairEnergy`],
27/// $`s_i`$ is the full set of site properties for site i, $`\vec{r}_i`$ is
28/// the position of site i, $`b_i`$ is the body tag that holds site *i*, and
29/// $`\left[ \  \right]`$ denotes the Iverson bracket.
30///
31/// In other words, [`PairwiseCutoff`] sums the energy for all pairs that are
32/// separated by a distance less than the maximum interaction range `r_cut` and
33/// belong to different bodies.
34///
35/// Use [`PairwiseCutoff`] with [`Anisotropic`], [`Isotropic`], [`HardShape`], or
36/// your own custom type.
37///
38/// TODO: Reword this when [`PairwiseCutoff`] also implements `SitePairForce`.
39///
40/// [`Anisotropic`]: crate::pairwise::Anisotropic
41/// [`Isotropic`]: crate::pairwise::Isotropic
42/// [`HardShape`]: crate::pairwise::HardShape
43///
44/// # Example
45///
46/// Basic usage:
47/// ```
48/// use hoomd_interaction::{
49///     PairwiseCutoff, pairwise::Isotropic, univariate::LennardJones,
50/// };
51///
52/// let lennard_jones: LennardJones = LennardJones {
53///     epsilon: 1.5,
54///     sigma: 2.0,
55/// };
56/// let pairwise_cutoff = PairwiseCutoff(Isotropic {
57///     interaction: lennard_jones,
58///     r_cut: 5.0,
59/// });
60/// ```
61///
62/// Set a custom potential using a closure:
63/// ```
64/// use hoomd_interaction::{PairwiseCutoff, pairwise::Isotropic};
65///
66/// let pairwise_cutoff = PairwiseCutoff(Isotropic {
67///     interaction: |r: f64| 1.0 / (r.powi(12)),
68///     r_cut: 3.0,
69/// });
70/// ```
71///
72/// Implement a custom potential via a type:
73/// ```
74/// use hoomd_interaction::{
75///     PairwiseCutoff, pairwise::Isotropic, univariate::UnivariateEnergy,
76/// };
77///
78/// struct Custom {
79///     a: f64,
80/// }
81///
82/// impl UnivariateEnergy for Custom {
83///     fn energy(&self, r: f64) -> f64 {
84///         self.a / r.powi(12)
85///     }
86/// }
87///
88/// let custom = Custom { a: 2.0 };
89/// let pairwise_cutoff = PairwiseCutoff(Isotropic {
90///     interaction: custom,
91///     r_cut: 2.0,
92/// });
93/// ```
94///
95/// Hard sphere:
96/// ```
97/// use hoomd_interaction::{PairwiseCutoff, pairwise::HardSphere};
98/// use hoomd_microstate::property::Point;
99/// use hoomd_vector::Cartesian;
100///
101/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
102/// let hard_sphere = PairwiseCutoff(HardSphere { diameter: 1.0 });
103/// # Ok(())
104/// # }
105/// ```
106///
107/// Hard ellipse:
108///
109/// ```
110/// use hoomd_geometry::shape::Ellipse;
111/// use hoomd_interaction::{PairwiseCutoff, pairwise::HardShape};
112/// use hoomd_microstate::property::Point;
113/// use hoomd_vector::Cartesian;
114/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
115/// let ellipse = Ellipse::with_semi_axes([4.0.try_into()?, 1.0.try_into()?]);
116/// let hard_ellipse = PairwiseCutoff(HardShape(ellipse));
117/// # Ok(())
118/// # }
119/// ```
120#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
121pub struct PairwiseCutoff<E>(pub E);
122
123impl<E> PairwiseCutoff<E> {
124    /// Compute the pair energy between two sites.
125    ///
126    /// Use this method to compute an individual term in the total pair energy,
127    /// subject to the the maximum interaction range `r_cut` and inter-body checks:
128    ///
129    /// ```math
130    /// U\left(s_i, s_j \right) \left[ \left|\vec{r}_j - \vec{r}_i\right| \lt r_\mathrm{cut} \right]\left[b_i \ne b_j\right]
131    /// ```
132    ///
133    /// # Example
134    /// ```
135    /// use approxim::assert_relative_eq;
136    ///
137    /// use hoomd_interaction::{
138    ///     PairwiseCutoff, pairwise::Isotropic, univariate::LennardJones,
139    /// };
140    /// use hoomd_microstate::{Body, Microstate, Site};
141    /// use hoomd_vector::Cartesian;
142    ///
143    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
144    /// let lennard_jones: LennardJones = LennardJones {
145    ///     epsilon: 1.0,
146    ///     sigma: 1.0,
147    /// };
148    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
149    ///     interaction: lennard_jones,
150    ///     r_cut: 2.5,
151    /// });
152    ///
153    /// let body_a = Body::point(Cartesian::from([0.0, 0.0]));
154    /// let body_b = Body::point(Cartesian::from([0.0, 3.0]));
155    /// let body_c = Body::point(Cartesian::from([0.0, -2.0_f64.powf(1.0 / 6.0)]));
156    ///
157    /// let microstate = Microstate::builder()
158    ///     .bodies([body_a, body_b, body_c])
159    ///     .try_build()?;
160    ///
161    /// let sites = microstate.sites();
162    /// let energy_ab = pairwise_cutoff.site_pair_energy(&sites[0], &sites[1]);
163    /// let energy_ac = pairwise_cutoff.site_pair_energy(&sites[0], &sites[2]);
164    ///
165    /// assert_eq!(energy_ab, 0.0);
166    /// assert_relative_eq!(energy_ac, -1.0);
167    /// # Ok(())
168    /// # }
169    /// ```
170    #[inline]
171    pub fn site_pair_energy<S>(&self, site_i: &Site<S>, site_j: &Site<S>) -> f64
172    where
173        E: SitePairEnergy<S>,
174    {
175        if site_i.body_tag == site_j.body_tag {
176            return 0.0;
177        }
178
179        self.0
180            .site_pair_energy(&site_i.properties, &site_j.properties)
181    }
182
183    /// Compute the filtered energy contribution of a single site (`AllPairs` specialization)
184    #[inline(always)]
185    fn filtered_site_energy_all<B, S, X, C, F, F2>(
186        &self,
187        microstate: &Microstate<B, S, X, C>,
188        site_i_properties: &S,
189        filter: F,
190        site_pair_energy: F2,
191    ) -> f64
192    where
193        E: SitePairEnergy<S>,
194        F: Fn(&Site<S>) -> bool,
195        F2: Fn(&E, &S, &S) -> f64,
196    {
197        let mut energy = 0.0;
198
199        for site_j in microstate.sites().iter().chain(microstate.ghosts()) {
200            if filter(site_j) {
201                let one = site_pair_energy(&self.0, site_i_properties, &site_j.properties);
202                if one == f64::INFINITY {
203                    return one;
204                }
205
206                energy += one;
207            }
208        }
209
210        energy
211    }
212
213    /// Compute the filtered energy contribution of a single site (spatial data specialization)
214    #[inline(always)]
215    fn filtered_site_energy_spatial<P, B, S, X, C, F, F2>(
216        &self,
217        microstate: &Microstate<B, S, X, C>,
218        site_i_properties: &S,
219        filter: F,
220        site_pair_energy: F2,
221    ) -> f64
222    where
223        E: SitePairEnergy<S> + MaximumInteractionRange,
224        S: Position<Position = P>,
225        X: PointsNearBall<P, SiteKey>,
226        F: Fn(&Site<S>) -> bool,
227        F2: Fn(&E, &S, &S) -> f64,
228    {
229        let mut energy = 0.0;
230
231        for site_j in microstate.iter_sites_near(
232            site_i_properties.position(),
233            self.0.maximum_interaction_range(),
234        ) {
235            if filter(site_j) {
236                let one = site_pair_energy(&self.0, site_i_properties, &site_j.properties);
237                if one == f64::INFINITY {
238                    return one;
239                }
240
241                energy += one;
242            }
243        }
244
245        energy
246    }
247
248    /// Compute the filtered energy contribution of a single site.
249    #[inline(always)]
250    fn filtered_site_energy<P, B, S, X, C, F, F2>(
251        &self,
252        microstate: &Microstate<B, S, X, C>,
253        site_i_properties: &S,
254        filter: F,
255        site_pair_energy: F2,
256    ) -> f64
257    where
258        E: SitePairEnergy<S> + MaximumInteractionRange,
259        S: Position<Position = P>,
260        X: PointsNearBall<P, SiteKey>,
261        F: Fn(&Site<S>) -> bool,
262        F2: Fn(&E, &S, &S) -> f64,
263    {
264        if X::is_all_pairs() {
265            self.filtered_site_energy_all(microstate, site_i_properties, filter, site_pair_energy)
266        } else {
267            self.filtered_site_energy_spatial(
268                microstate,
269                site_i_properties,
270                filter,
271                site_pair_energy,
272            )
273        }
274    }
275
276    /// Compute the final energy of a body in the microstate.
277    #[inline(always)]
278    fn filtered_body_energy_final<P, B, S, X, C, F>(
279        &self,
280        microstate: &Microstate<B, S, X, C>,
281        body: &Body<B, S>,
282        filter: F,
283    ) -> f64
284    where
285        E: SitePairEnergy<S> + MaximumInteractionRange,
286        B: Transform<S>,
287        S: Position<Position = P>,
288        X: PointsNearBall<P, SiteKey>,
289        C: Wrap<B> + Wrap<S>,
290        F: Fn(&Site<S>) -> bool,
291    {
292        let mut energy_final = 0.0;
293        for s in &body.sites {
294            match microstate.boundary().wrap(body.properties.transform(s)) {
295                Err(_) => return f64::INFINITY,
296                Ok(site_i_properties) => {
297                    let one = self.filtered_site_energy(
298                        microstate,
299                        &site_i_properties,
300                        &filter,
301                        E::site_pair_energy,
302                    );
303                    if one == f64::INFINITY {
304                        return one;
305                    }
306
307                    energy_final += one;
308                }
309            }
310        }
311        energy_final
312    }
313
314    /// Compute the initial energy of a body in the microstate.
315    #[inline(always)]
316    fn filtered_body_energy_initial<P, B, S, X, C, F>(
317        &self,
318        microstate: &Microstate<B, S, X, C>,
319        body_index: usize,
320        filter: F,
321    ) -> f64
322    where
323        E: SitePairEnergy<S> + MaximumInteractionRange,
324        S: Position<Position = P>,
325        X: PointsNearBall<P, SiteKey>,
326        F: Fn(&Site<S>) -> bool,
327    {
328        let mut energy_initial = 0.0;
329        if !E::is_only_infinite_or_zero() {
330            for site_i in microstate.iter_body_sites(body_index) {
331                let one = self.filtered_site_energy(
332                    microstate,
333                    &site_i.properties,
334                    &filter,
335                    E::site_pair_energy_initial,
336                );
337                if one == f64::INFINITY {
338                    return one;
339                }
340
341                energy_initial += one;
342            }
343        }
344        energy_initial
345    }
346}
347
348impl<P, B, S, X, C, E> TotalEnergy<Microstate<B, S, X, C>> for PairwiseCutoff<E>
349where
350    E: SitePairEnergy<S> + MaximumInteractionRange,
351    S: Position<Position = P>,
352    X: PointsNearBall<P, SiteKey>,
353{
354    /// Compute the total energy of the microstate contributed by functions on pairs of sites.
355    ///
356    /// # Examples
357    ///
358    /// Lennard-Jones:
359    /// ```
360    /// use hoomd_interaction::{
361    ///     PairwiseCutoff, SitePairEnergy, TotalEnergy, pairwise::Isotropic,
362    ///     univariate::LennardJones,
363    /// };
364    /// use hoomd_microstate::{
365    ///     Body, Microstate,
366    ///     property::{Point, Position},
367    /// };
368    /// use hoomd_vector::{Cartesian, InnerProduct};
369    ///
370    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
371    /// let mut microstate = Microstate::new();
372    /// microstate.extend_bodies([
373    ///     Body::point(Cartesian::from([0.0, 0.0])),
374    ///     Body::point(Cartesian::from([1.0, 0.0])),
375    ///     Body::point(Cartesian::from([0.0, 5.0])),
376    ///     Body::point(Cartesian::from([-1.0, 5.0])),
377    /// ])?;
378    ///
379    /// let lennard_jones: LennardJones = LennardJones {
380    ///     epsilon: 1.5,
381    ///     sigma: 1.0 / 2.0_f64.powf(1.0 / 6.0),
382    /// };
383    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
384    ///     interaction: lennard_jones,
385    ///     r_cut: 2.5,
386    /// });
387    ///
388    /// let total_energy = pairwise_cutoff.total_energy(&microstate);
389    /// assert_eq!(total_energy, -3.0);
390    /// # Ok(())
391    /// # }
392    /// ```
393    ///
394    /// Hard circle:
395    /// ```
396    /// use hoomd_geometry::shape::Circle;
397    /// use hoomd_interaction::{
398    ///     PairwiseCutoff, TotalEnergy, pairwise::HardSphere,
399    /// };
400    /// use hoomd_microstate::{Body, Microstate, property::Point};
401    /// use hoomd_vector::{Angle, Cartesian};
402    ///
403    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
404    /// let mut microstate = Microstate::new();
405    /// microstate.extend_bodies([
406    ///     Body::point(Cartesian::from([0.0, 0.0])),
407    ///     Body::point(Cartesian::from([0.4, 0.0])),
408    /// ])?;
409    ///
410    /// let hard_circle = PairwiseCutoff(HardSphere { diameter: 1.0 });
411    ///
412    /// let total_energy = hard_circle.total_energy(&microstate);
413    /// assert_eq!(total_energy, f64::INFINITY);
414    ///
415    /// microstate.update_body_properties(0, Point::new([0.0, -2.0].into()));
416    /// let total_energy = hard_circle.total_energy(&microstate);
417    /// assert_eq!(total_energy, 0.0);
418    /// # Ok(())
419    /// # }
420    /// ```
421    #[inline]
422    fn total_energy(&self, microstate: &Microstate<B, S, X, C>) -> f64 {
423        let mut total = 0.0;
424
425        // If needed, total_energy could specialize further in the all-pairs
426        // code path. The current implementation performs many unneeded
427        // site_i.site_tag < site_j.site_tag checks. However, the solution is
428        // non-trivial. It would need to loop over the j sites *by tag*
429        // to avoid looping over unneeded sites. The ghost loop would need
430        // to be separate and include the tag filter.
431
432        for site_i in microstate.sites() {
433            let one = self.filtered_site_energy(
434                microstate,
435                &site_i.properties,
436                |site_j| site_i.site_tag < site_j.site_tag && site_i.body_tag != site_j.body_tag,
437                E::site_pair_energy,
438            );
439            if one == f64::INFINITY {
440                return one;
441            }
442
443            total += one;
444        }
445
446        total
447    }
448
449    /// Compute the difference in energy between two microstates.
450    ///
451    /// Returns `$ E_\mathrm{final} - E_\mathrm{initial} $`.
452    ///
453    /// # Example
454    ///
455    /// ```
456    /// use hoomd_interaction::{
457    ///     PairwiseCutoff, SitePairEnergy, TotalEnergy, pairwise::Isotropic,
458    ///     univariate::LennardJones,
459    /// };
460    /// use hoomd_microstate::{
461    ///     Body, Microstate,
462    ///     property::{Point, Position},
463    /// };
464    /// use hoomd_vector::{Cartesian, InnerProduct};
465    ///
466    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
467    /// let mut microstate_a = Microstate::new();
468    /// microstate_a.extend_bodies([
469    ///     Body::point(Cartesian::from([0.0, 0.0])),
470    ///     Body::point(Cartesian::from([1.0, 0.0])),
471    /// ])?;
472    ///
473    /// let mut microstate_b = Microstate::new();
474    /// microstate_b.extend_bodies([
475    ///     Body::point(Cartesian::from([0.0, 0.0])),
476    ///     Body::point(Cartesian::from([5.0, 0.0])),
477    /// ])?;
478    ///
479    /// let lennard_jones: LennardJones = LennardJones {
480    ///     epsilon: 1.5,
481    ///     sigma: 1.0 / 2.0_f64.powf(1.0 / 6.0),
482    /// };
483    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
484    ///     interaction: lennard_jones,
485    ///     r_cut: 2.5,
486    /// });
487    ///
488    /// let delta_energy_total =
489    ///     pairwise_cutoff.delta_energy_total(&microstate_a, &microstate_b);
490    /// assert_eq!(delta_energy_total, 1.5);
491    /// # Ok(())
492    /// # }
493    /// ```
494    #[inline]
495    fn delta_energy_total(
496        &self,
497        initial_microstate: &Microstate<B, S, X, C>,
498        final_microstate: &Microstate<B, S, X, C>,
499    ) -> f64 {
500        let mut energy_final = 0.0;
501
502        for site_i in final_microstate.sites() {
503            let one = self.filtered_site_energy(
504                final_microstate,
505                &site_i.properties,
506                |site_j| site_i.site_tag < site_j.site_tag && site_i.body_tag != site_j.body_tag,
507                E::site_pair_energy,
508            );
509            if one == f64::INFINITY {
510                return one;
511            }
512            energy_final += one;
513        }
514
515        let mut energy_initial = 0.0;
516        if !E::is_only_infinite_or_zero() {
517            for site_i in initial_microstate.sites() {
518                let one = self.filtered_site_energy(
519                    initial_microstate,
520                    &site_i.properties,
521                    |site_j| {
522                        site_i.site_tag < site_j.site_tag && site_i.body_tag != site_j.body_tag
523                    },
524                    E::site_pair_energy_initial,
525                );
526                if one == f64::INFINITY {
527                    return -one;
528                }
529                energy_initial += one;
530            }
531        }
532
533        energy_final - energy_initial
534    }
535}
536
537impl<P, B, S, X, C, E> DeltaEnergyOne<B, S, X, C> for PairwiseCutoff<E>
538where
539    E: SitePairEnergy<S> + MaximumInteractionRange,
540    B: Transform<S>,
541    S: Position<Position = P>,
542    X: PointsNearBall<P, SiteKey>,
543    C: Wrap<B> + Wrap<S>,
544{
545    /// Evaluate the change in energy contributed by `PairwiseCutoff` when one body is updated.
546    ///
547    /// # Examples
548    ///
549    /// Boxcar:
550    /// ```
551    /// use hoomd_interaction::{
552    ///     DeltaEnergyOne, PairwiseCutoff, pairwise::Isotropic, univariate::Boxcar,
553    /// };
554    /// use hoomd_microstate::{Body, Microstate, property::Point};
555    /// use hoomd_vector::Cartesian;
556    ///
557    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
558    /// let mut microstate = Microstate::new();
559    /// microstate.extend_bodies([
560    ///     Body::point(Cartesian::from([0.0, 0.0])),
561    ///     Body::point(Cartesian::from([1.0, 0.0])),
562    /// ])?;
563    ///
564    /// let epsilon = 2.0;
565    /// let (left, right) = (0.0, 1.5);
566    /// let boxcar = Boxcar {
567    ///     epsilon,
568    ///     left,
569    ///     right,
570    /// };
571    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
572    ///     interaction: boxcar,
573    ///     r_cut: 1.5,
574    /// });
575    ///
576    /// let delta_energy = pairwise_cutoff.delta_energy_one(
577    ///     &microstate,
578    ///     0,
579    ///     &Body::point([-1.0, 0.0].into()),
580    /// );
581    /// assert_eq!(delta_energy, -2.0);
582    /// # Ok(())
583    /// # }
584    /// ```
585    ///
586    /// Hard circle:
587    /// ```
588    /// use hoomd_geometry::shape::Circle;
589    /// use hoomd_interaction::{
590    ///     DeltaEnergyOne, PairwiseCutoff, pairwise::HardSphere,
591    /// };
592    /// use hoomd_microstate::{Body, Microstate, property::Point};
593    /// use hoomd_vector::{Angle, Cartesian};
594    ///
595    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
596    /// let mut microstate = Microstate::new();
597    /// microstate.extend_bodies([
598    ///     Body::point(Cartesian::from([0.0, 0.0])),
599    ///     Body::point(Cartesian::from([2.0, 0.0])),
600    /// ])?;
601    ///
602    /// let hard_circle = PairwiseCutoff(HardSphere { diameter: 1.0 });
603    ///
604    /// let delta_energy = hard_circle.delta_energy_one(
605    ///     &microstate,
606    ///     1,
607    ///     &Body::point([0.4, 0.0].into()),
608    /// );
609    /// assert_eq!(delta_energy, f64::INFINITY);
610    ///
611    /// let delta_energy = hard_circle.delta_energy_one(
612    ///     &microstate,
613    ///     1,
614    ///     &Body::point([1.5, 0.0].into()),
615    /// );
616    /// assert_eq!(delta_energy, 0.0);
617    /// # Ok(())
618    /// # }
619    /// ```
620    #[inline]
621    fn delta_energy_one(
622        &self,
623        initial_microstate: &Microstate<B, S, X, C>,
624        body_index: usize,
625        final_body: &Body<B, S>,
626    ) -> f64 {
627        let body_tag = initial_microstate.bodies()[body_index].tag;
628
629        let energy_final =
630            self.filtered_body_energy_final(initial_microstate, final_body, |site_j| {
631                body_tag != site_j.body_tag
632            });
633        if energy_final == f64::INFINITY {
634            return energy_final;
635        }
636
637        let energy_initial =
638            self.filtered_body_energy_initial(initial_microstate, body_index, |site_j| {
639                body_tag != site_j.body_tag
640            });
641
642        energy_final - energy_initial
643    }
644}
645
646impl<P, B, S, X, C, E> DeltaEnergyInsert<B, S, X, C> for PairwiseCutoff<E>
647where
648    E: SitePairEnergy<S> + MaximumInteractionRange,
649    B: Transform<S>,
650    S: Position<Position = P>,
651    X: PointsNearBall<P, SiteKey>,
652    C: Wrap<B> + Wrap<S>,
653{
654    /// Evaluate the change in energy contributed by `PairwiseCutoff` when one body is inserted.
655    ///
656    /// # Example
657    ///
658    /// Boxcar:
659    /// ```
660    /// use hoomd_interaction::{
661    ///     DeltaEnergyInsert, PairwiseCutoff, pairwise::Isotropic,
662    ///     univariate::Boxcar,
663    /// };
664    /// use hoomd_microstate::{Body, Microstate, property::Point};
665    /// use hoomd_vector::Cartesian;
666    ///
667    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
668    /// let mut microstate = Microstate::new();
669    /// microstate.extend_bodies([
670    ///     Body::point(Cartesian::from([0.0, 0.0])),
671    ///     Body::point(Cartesian::from([1.0, 0.0])),
672    /// ])?;
673    ///
674    /// let epsilon = 2.0;
675    /// let (left, right) = (0.0, 1.5);
676    /// let boxcar = Boxcar {
677    ///     epsilon,
678    ///     left,
679    ///     right,
680    /// };
681    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
682    ///     interaction: boxcar,
683    ///     r_cut: 1.5,
684    /// });
685    ///
686    /// let delta_energy = pairwise_cutoff
687    ///     .delta_energy_insert(&microstate, &Body::point([-1.0, 0.0].into()));
688    /// assert_eq!(delta_energy, 2.0);
689    /// # Ok(())
690    /// # }
691    /// ```
692    ///
693    /// Hard circle:
694    /// ```
695    /// use hoomd_geometry::shape::Circle;
696    /// use hoomd_interaction::{
697    ///     DeltaEnergyInsert, PairwiseCutoff, pairwise::HardSphere,
698    /// };
699    /// use hoomd_microstate::{Body, Microstate, property::Point};
700    /// use hoomd_vector::{Angle, Cartesian};
701    ///
702    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
703    /// let mut microstate = Microstate::new();
704    /// microstate.extend_bodies([Body::point(Cartesian::from([0.0, 0.0]))])?;
705    ///
706    /// let hard_circle = PairwiseCutoff(HardSphere { diameter: 1.0 });
707    ///
708    /// let delta_energy = hard_circle
709    ///     .delta_energy_insert(&microstate, &Body::point([0.4, 0.0].into()));
710    /// assert_eq!(delta_energy, f64::INFINITY);
711    ///
712    /// let delta_energy = hard_circle
713    ///     .delta_energy_insert(&microstate, &Body::point([1.5, 0.0].into()));
714    /// assert_eq!(delta_energy, 0.0);
715    /// # Ok(())
716    /// # }
717    /// ```
718    #[inline]
719    fn delta_energy_insert(
720        &self,
721        initial_microstate: &Microstate<B, S, X, C>,
722        new_body: &Body<B, S>,
723    ) -> f64 {
724        // The new body is not yet in the microstate, so there is no need to
725        // filter matching body tags. The new body does not yet have a tag.
726        self.filtered_body_energy_final(initial_microstate, new_body, |_| true)
727    }
728}
729
730impl<P, B, S, X, C, E> DeltaEnergyRemove<B, S, X, C> for PairwiseCutoff<E>
731where
732    E: SitePairEnergy<S> + MaximumInteractionRange,
733    S: Position<Position = P>,
734    X: PointsNearBall<P, SiteKey>,
735{
736    /// Evaluate the change in energy contributed by `PairwiseCutoff` when one body is removed.
737    ///
738    /// # Example
739    ///
740    /// Boxcar:
741    /// ```
742    /// use hoomd_interaction::{
743    ///     DeltaEnergyRemove, PairwiseCutoff, pairwise::Isotropic,
744    ///     univariate::Boxcar,
745    /// };
746    /// use hoomd_microstate::{Body, Microstate, property::Point};
747    /// use hoomd_vector::Cartesian;
748    ///
749    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
750    /// let mut microstate = Microstate::new();
751    /// microstate.extend_bodies([
752    ///     Body::point(Cartesian::from([0.0, 0.0])),
753    ///     Body::point(Cartesian::from([1.0, 0.0])),
754    /// ])?;
755    ///
756    /// let epsilon = 2.0;
757    /// let (left, right) = (0.0, 1.5);
758    /// let boxcar = Boxcar {
759    ///     epsilon,
760    ///     left,
761    ///     right,
762    /// };
763    /// let pairwise_cutoff = PairwiseCutoff(Isotropic {
764    ///     interaction: boxcar,
765    ///     r_cut: 1.5,
766    /// });
767    ///
768    /// let delta_energy = pairwise_cutoff.delta_energy_remove(&microstate, 0);
769    /// assert_eq!(delta_energy, -2.0);
770    /// # Ok(())
771    /// # }
772    /// ```
773    ///
774    /// Hard circle:
775    /// ```
776    /// use hoomd_geometry::shape::Circle;
777    /// use hoomd_interaction::{
778    ///     DeltaEnergyRemove, PairwiseCutoff, pairwise::HardSphere,
779    /// };
780    /// use hoomd_microstate::{Body, Microstate, property::Point};
781    /// use hoomd_vector::{Angle, Cartesian};
782    ///
783    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
784    /// let mut microstate = Microstate::new();
785    /// microstate.extend_bodies([
786    ///     Body::point(Cartesian::from([0.0, 0.0])),
787    ///     Body::point(Cartesian::from([2.0, 0.0])),
788    /// ])?;
789    ///
790    /// let hard_circle = PairwiseCutoff(HardSphere { diameter: 1.0 });
791    ///
792    /// let delta_energy = hard_circle.delta_energy_remove(&microstate, 1);
793    /// assert_eq!(delta_energy, 0.0);
794    /// # Ok(())
795    /// # }
796    /// ```
797    #[inline]
798    fn delta_energy_remove(
799        &self,
800        initial_microstate: &Microstate<B, S, X, C>,
801        body_index: usize,
802    ) -> f64 {
803        let body_tag = initial_microstate.bodies()[body_index].tag;
804        let energy_initial =
805            self.filtered_body_energy_initial(initial_microstate, body_index, |site_j| {
806                body_tag != site_j.body_tag
807            });
808
809        -energy_initial
810    }
811}
812
813#[cfg(test)]
814mod tests_finite {
815    use super::*;
816    use crate::{TotalEnergy, pairwise::Isotropic, univariate::HarmonicRepulsion};
817    use assert2::check;
818    use hoomd_geometry::shape::Hypercuboid;
819    use hoomd_microstate::{
820        boundary::{Closed, Open},
821        property::Point,
822    };
823    use hoomd_spatial::{AllPairs, VecCell};
824    use hoomd_vector::{Cartesian, distribution::Ball};
825
826    use approxim::assert_relative_eq;
827    use rand::{
828        RngExt, SeedableRng,
829        distr::{Distribution, Uniform},
830        rngs::StdRng,
831    };
832    use rstest::*;
833    use std::f64::consts::PI;
834
835    #[fixture]
836    fn square() -> Closed<Hypercuboid<2>> {
837        let cuboid = Hypercuboid {
838            edge_lengths: [
839                4.0.try_into()
840                    .expect("hard-coded constant should be positive"),
841                4.0.try_into()
842                    .expect("hard-coded constant should be positive"),
843            ],
844        };
845        Closed(cuboid)
846    }
847
848    mod pairwise_cutoff {
849        use super::*;
850        use crate::pairwise::Isotropic;
851
852        #[fixture]
853        fn microstate()
854        -> Microstate<Point<Cartesian<2>>, Point<Cartesian<2>>, AllPairs<SiteKey>, Open> {
855            let mut microstate = Microstate::new();
856            microstate
857                .extend_bodies([
858                    Body::point(Cartesian::from([0.0, 0.0])),
859                    Body::point(Cartesian::from([1.0, 0.0])),
860                    Body::point(Cartesian::from([0.0, 5.0])),
861                    Body::point(Cartesian::from([1.0, 5.0])),
862                ])
863                .expect("hard-coded bodies should be in the boundary");
864            microstate
865        }
866
867        #[rstest]
868        fn blanket_fn(
869            microstate: Microstate<
870                Point<Cartesian<2>>,
871                Point<Cartesian<2>>,
872                AllPairs<SiteKey>,
873                Open,
874            >,
875        ) {
876            // Ensure that closures can be used as UnivariateEnergy
877            let pairwise_cutoff = PairwiseCutoff(Isotropic {
878                interaction: |r| 1.0 / (r * 2.0),
879                r_cut: 2.0,
880            });
881
882            // Two pairs at a distance of 1.0 each with energy 1/2.
883            assert_eq!(pairwise_cutoff.total_energy(&microstate), 1.0);
884
885            let sites = microstate.sites();
886            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[0]) == 0.0);
887            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[1]) == 0.5);
888            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[2]) == 0.0);
889            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[3]) == 0.0);
890            check!(pairwise_cutoff.site_pair_energy(&sites[1], &sites[1]) == 0.0);
891            check!(pairwise_cutoff.site_pair_energy(&sites[1], &sites[2]) == 0.0);
892            check!(pairwise_cutoff.site_pair_energy(&sites[1], &sites[3]) == 0.0);
893            check!(pairwise_cutoff.site_pair_energy(&sites[2], &sites[2]) == 0.0);
894            check!(pairwise_cutoff.site_pair_energy(&sites[2], &sites[3]) == 0.5);
895            check!(pairwise_cutoff.site_pair_energy(&sites[3], &sites[3]) == 0.0);
896        }
897
898        #[rstest]
899        fn large_r_cut(
900            microstate: Microstate<
901                Point<Cartesian<2>>,
902                Point<Cartesian<2>>,
903                AllPairs<SiteKey>,
904                Open,
905            >,
906        ) {
907            // Ensure that PairwiseCutoff respects the r_cut value set.
908            let pairwise_cutoff = PairwiseCutoff(Isotropic {
909                interaction: |r| 1.0 / (r * 2.0),
910                r_cut: 5.0_f64.next_up(),
911            });
912
913            // Two pairs at a distance of 1.0 each with energy 1/2.
914            // Plus two pairs at a distance of 5.0 with energy 1/10
915            check!(pairwise_cutoff.total_energy(&microstate) == 1.2);
916        }
917
918        #[test]
919        fn body_exclusion() -> anyhow::Result<()> {
920            // Ensure that PairwiseCutoff excludes pairs in the same body.
921            let body_a = Body {
922                properties: Point::new(Cartesian::from([0.0, 0.0])),
923                sites: [
924                    Point::new(Cartesian::from([1.0, 1.0])),
925                    Point::new(Cartesian::from([1.0, -1.0])),
926                    Point::new(Cartesian::from([-1.0, 1.0])),
927                    Point::new(Cartesian::from([-1.0, -1.0])),
928                ]
929                .into(),
930            };
931            let body_b = Body {
932                properties: Point::new(Cartesian::from([3.0, 0.0])),
933                sites: body_a.sites.clone(),
934            };
935
936            let mut microstate = Microstate::new();
937            microstate.extend_bodies([body_a, body_b])?;
938
939            let pairwise_cutoff = PairwiseCutoff(Isotropic {
940                interaction: |_r| 1.0,
941                r_cut: 1.0_f64.next_up(),
942            });
943
944            // Of all the pairs a distance 1.0 apart, only 2 are interbody pairs.
945            check!(pairwise_cutoff.total_energy(&microstate) == 2.0);
946
947            let sites = microstate.sites();
948            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[0]) == 0.0);
949            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[1]) == 0.0);
950            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[2]) == 0.0);
951            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[3]) == 0.0);
952
953            check!(pairwise_cutoff.site_pair_energy(&sites[4], &sites[4]) == 0.0);
954            check!(pairwise_cutoff.site_pair_energy(&sites[4], &sites[5]) == 0.0);
955            check!(pairwise_cutoff.site_pair_energy(&sites[4], &sites[6]) == 0.0);
956            check!(pairwise_cutoff.site_pair_energy(&sites[4], &sites[7]) == 0.0);
957
958            check!(pairwise_cutoff.site_pair_energy(&sites[0], &sites[6]) == 1.0);
959            check!(pairwise_cutoff.site_pair_energy(&sites[1], &sites[7]) == 1.0);
960
961            Ok(())
962        }
963    }
964
965    mod delta_energy_one {
966        use super::*;
967
968        #[rstest]
969        fn site_outside(square: Closed<Hypercuboid<2>>) -> anyhow::Result<()> {
970            let body = Body {
971                properties: Point::new(Cartesian::from([0.0, 0.0])),
972                sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
973            };
974            let mut final_body = body.clone();
975            final_body.properties.position[0] = 1.0;
976
977            let microstate = Microstate::builder()
978                .boundary(square)
979                .bodies([body])
980                .try_build()?;
981
982            let energy = PairwiseCutoff(Isotropic {
983                interaction: |_r| 0.0,
984                r_cut: 0.0,
985            });
986
987            check!(energy.delta_energy_one(&microstate, 0, &final_body) == f64::INFINITY);
988
989            Ok(())
990        }
991
992        #[test]
993        fn body_exclusion() -> anyhow::Result<()> {
994            // Ensure that PairwiseCutoff.delta_energy_one excludes pairs in the same body.
995            let body_a = Body {
996                properties: Point::new(Cartesian::from([0.0, 0.0])),
997                sites: [
998                    Point::new(Cartesian::from([1.0, 1.0])),
999                    Point::new(Cartesian::from([1.0, -1.0])),
1000                    Point::new(Cartesian::from([-1.0, 1.0])),
1001                    Point::new(Cartesian::from([-1.0, -1.0])),
1002                ]
1003                .into(),
1004            };
1005            let body_b = Body {
1006                properties: Point::new(Cartesian::from([3.0, 0.0])),
1007                sites: body_a.sites.clone(),
1008            };
1009            let body_a_final = Body {
1010                properties: Point::new(Cartesian::from([-1.0, 0.0])),
1011                sites: body_a.sites.clone(),
1012            };
1013
1014            let mut microstate = Microstate::new();
1015            microstate.extend_bodies([body_a, body_b])?;
1016
1017            let pairwise_cutoff = PairwiseCutoff(Isotropic {
1018                interaction: |_r| 1.0,
1019                r_cut: 1.0_f64.next_up(),
1020            });
1021
1022            // Of all the pairs a distance 1.0 apart, only 2 are interbody pairs.
1023            // Moving body 0 to the left results in a -2.0 energy difference.
1024            check!(pairwise_cutoff.delta_energy_one(&microstate, 0, &body_a_final) == -2.0);
1025
1026            Ok(())
1027        }
1028
1029        #[test]
1030        fn random_moves() -> anyhow::Result<()> {
1031            // Ensure that PairwiseCutoff.delta_energy_one is consistent with TotalEnergy
1032            let body_template = Body {
1033                properties: Point::new(Cartesian::from([0.0, 0.0])),
1034                sites: [
1035                    Point::new(Cartesian::from([0.0, 1.0])),
1036                    Point::new(Cartesian::from([-1.0, 1.0])),
1037                    Point::new(Cartesian::from([-1.0, -1.0])),
1038                ]
1039                .into(),
1040            };
1041            let body_a = Body {
1042                properties: Point::new(Cartesian::from([3.0, 0.0])),
1043                sites: body_template.sites.clone(),
1044            };
1045            let body_b = body_template.clone();
1046
1047            let microstate_initial = Microstate::builder().bodies([body_a, body_b]).try_build()?;
1048
1049            let mut microstate_final = microstate_initial.clone();
1050            let harmonic_repulsion: HarmonicRepulsion = HarmonicRepulsion { a: 5.0, r_cut: 5.0 };
1051            let pairwise_cutoff = PairwiseCutoff(Isotropic {
1052                interaction: harmonic_repulsion,
1053                r_cut: 5.0,
1054            });
1055
1056            check!(pairwise_cutoff.total_energy(&microstate_initial) != 0.0);
1057
1058            // Use `HarmonicRepulsion` for validation because it is a varies
1059            // with r and will therefore show some changes for any moves (unlike
1060            // `BoxCar`). HarmonicRepulsion avoids numerical errors when two
1061            // sites get too close.
1062            let mut rng = StdRng::seed_from_u64(0);
1063            let r_distribution = Uniform::new(3.0, 6.0)?;
1064            let theta_distribution = Uniform::new(0.0, 2.0 * PI)?;
1065
1066            let mut new_body = body_template.clone();
1067            for _ in 0..1024 {
1068                let r = rng.sample(r_distribution);
1069                let theta = rng.sample(theta_distribution);
1070                new_body.properties.position = [r * theta.cos(), r * theta.sin()].into();
1071
1072                let delta_energy_one =
1073                    pairwise_cutoff.delta_energy_one(&microstate_initial, 0, &new_body);
1074                microstate_final.update_body_properties(0, new_body.properties)?;
1075                let delta_energy_total = pairwise_cutoff.total_energy(&microstate_final)
1076                    - pairwise_cutoff.total_energy(&microstate_initial);
1077
1078                assert_relative_eq!(delta_energy_one, delta_energy_total, epsilon = 1e-10);
1079                assert_relative_eq!(
1080                    pairwise_cutoff.delta_energy_total(&microstate_initial, &microstate_final),
1081                    delta_energy_total,
1082                    epsilon = 1e-10
1083                );
1084            }
1085
1086            Ok(())
1087        }
1088    }
1089
1090    mod delta_energy_insert_remove {
1091        use super::*;
1092
1093        #[rstest]
1094        fn site_outside(square: Closed<Hypercuboid<2>>) -> anyhow::Result<()> {
1095            let body = Body {
1096                properties: Point::new(Cartesian::from([0.0, 0.0])),
1097                sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
1098            };
1099            let mut new_body = body.clone();
1100            new_body.properties.position[0] = 1.0;
1101
1102            let microstate = Microstate::builder()
1103                .boundary(square)
1104                .bodies([body])
1105                .try_build()?;
1106
1107            let energy = PairwiseCutoff(Isotropic {
1108                interaction: |_r| 0.0,
1109                r_cut: 0.0,
1110            });
1111
1112            check!(energy.delta_energy_insert(&microstate, &new_body) == f64::INFINITY);
1113
1114            Ok(())
1115        }
1116
1117        #[test]
1118        fn body_exclusion() -> anyhow::Result<()> {
1119            // Ensure that PairwiseCutoff.delta_energy_insert excludes pairs in the same body.
1120            let body_a_new = Body {
1121                properties: Point::new(Cartesian::from([0.0, 0.0])),
1122                sites: [
1123                    Point::new(Cartesian::from([1.0, 1.0])),
1124                    Point::new(Cartesian::from([1.0, -1.0])),
1125                    Point::new(Cartesian::from([-1.0, 1.0])),
1126                    Point::new(Cartesian::from([-1.0, -1.0])),
1127                ]
1128                .into(),
1129            };
1130            let body_b = Body {
1131                properties: Point::new(Cartesian::from([3.0, 0.0])),
1132                sites: body_a_new.sites.clone(),
1133            };
1134
1135            let mut microstate = Microstate::new();
1136            microstate.extend_bodies([body_b])?;
1137
1138            let pairwise_cutoff = PairwiseCutoff(Isotropic {
1139                interaction: |_r| 1.0,
1140                r_cut: 1.0_f64.next_up(),
1141            });
1142
1143            // Of all the pairs a distance 1.0 apart, only 2 are interbody pairs.
1144            // Moving body 0 to the left results in a -2.0 energy difference.
1145            check!(pairwise_cutoff.delta_energy_insert(&microstate, &body_a_new) == 2.0);
1146
1147            microstate.add_body(body_a_new)?;
1148            check!(pairwise_cutoff.delta_energy_remove(&microstate, 1) == -2.0);
1149
1150            Ok(())
1151        }
1152
1153        #[test]
1154        fn random_moves() -> anyhow::Result<()> {
1155            // Ensure that PairwiseCutoff.delta_energy_insert is consistent with TotalEnergy
1156            let body_template = Body {
1157                properties: Point::new(Cartesian::from([0.0, 0.0])),
1158                sites: [
1159                    Point::new(Cartesian::from([0.0, 1.0])),
1160                    Point::new(Cartesian::from([-1.0, 1.0])),
1161                    Point::new(Cartesian::from([-1.0, -1.0])),
1162                ]
1163                .into(),
1164            };
1165            let body_a = Body {
1166                properties: Point::new(Cartesian::from([3.0, 0.0])),
1167                sites: body_template.sites.clone(),
1168            };
1169
1170            let microstate_initial = Microstate::builder().bodies([body_a]).try_build()?;
1171
1172            let mut microstate_final = microstate_initial.clone();
1173            let harmonic_repulsion: HarmonicRepulsion = HarmonicRepulsion { a: 5.0, r_cut: 5.0 };
1174            let pairwise_cutoff = PairwiseCutoff(Isotropic {
1175                interaction: harmonic_repulsion,
1176                r_cut: 5.0,
1177            });
1178
1179            // Use `HarmonicRepulsion` for validation because it is a varies
1180            // with r and will therefore show some changes for any moves (unlike
1181            // `BoxCar`). HarmonicRepulsion avoids numerical errors when two
1182            // sites get too close.
1183            let mut rng = StdRng::seed_from_u64(2);
1184            let r_distribution = Uniform::new(3.0, 6.0)?;
1185            let theta_distribution = Uniform::new(0.0, 2.0 * PI)?;
1186
1187            for _ in 0..1024 {
1188                let r = rng.sample(r_distribution);
1189                let theta = rng.sample(theta_distribution);
1190                let mut new_body = body_template.clone();
1191                new_body.properties.position = [r * theta.cos(), r * theta.sin()].into();
1192
1193                let delta_energy_insert =
1194                    pairwise_cutoff.delta_energy_insert(&microstate_initial, &new_body);
1195                let tag = microstate_final.add_body(new_body)?;
1196                let delta_energy_total = pairwise_cutoff.total_energy(&microstate_final)
1197                    - pairwise_cutoff.total_energy(&microstate_initial);
1198
1199                assert_relative_eq!(delta_energy_insert, delta_energy_total, epsilon = 1e-6);
1200
1201                let delta_energy_remove = pairwise_cutoff.delta_energy_remove(&microstate_final, 1);
1202                assert_relative_eq!(delta_energy_remove, -delta_energy_total, epsilon = 1e-6);
1203
1204                microstate_final.remove_body(
1205                    microstate_final.body_indices()[tag].expect("tag should be present"),
1206                );
1207            }
1208
1209            Ok(())
1210        }
1211    }
1212
1213    #[rstest]
1214    fn spatial_data_consistency(square: Closed<Hypercuboid<2>>) -> anyhow::Result<()> {
1215        const N_BODIES: usize = 2_000;
1216        let r_cut = 0.5;
1217        let mut rng = StdRng::seed_from_u64(0);
1218
1219        let mut microstate_all_pairs = Microstate::builder()
1220            .spatial_data(AllPairs::<SiteKey>::default())
1221            .boundary(square.clone())
1222            .try_build()?;
1223
1224        let cell_list = VecCell::builder()
1225            .nominal_search_radius(r_cut.try_into()?)
1226            .build();
1227        let mut microstate_vec_cell = Microstate::builder()
1228            .boundary(square.clone())
1229            .spatial_data(cell_list)
1230            .try_build()?;
1231
1232        let body_template = Body {
1233            properties: Point::new(Cartesian::from([0.0, 0.0])),
1234            sites: [Point::default()].into(),
1235        };
1236        for _ in 0..N_BODIES {
1237            let mut new_body = body_template.clone();
1238            new_body.properties.position = square.sample(&mut rng);
1239            microstate_all_pairs.add_body(new_body.clone())?;
1240            microstate_vec_cell.add_body(new_body)?;
1241        }
1242
1243        let harmonic_repulsion: HarmonicRepulsion = HarmonicRepulsion { a: 5.0, r_cut };
1244        let pairwise_cutoff = PairwiseCutoff(Isotropic {
1245            interaction: harmonic_repulsion,
1246            r_cut,
1247        });
1248
1249        assert_relative_eq!(
1250            pairwise_cutoff.total_energy(&microstate_all_pairs),
1251            pairwise_cutoff.total_energy(&microstate_vec_cell)
1252        );
1253
1254        let move_distribution = Ball {
1255            radius: 0.2.try_into()?,
1256        };
1257        for i in (0..N_BODIES).step_by(4) {
1258            assert_relative_eq!(
1259                pairwise_cutoff.delta_energy_remove(&microstate_all_pairs, i),
1260                pairwise_cutoff.delta_energy_remove(&microstate_vec_cell, i),
1261                epsilon = 1e-10
1262            );
1263
1264            let mut final_body = microstate_all_pairs.bodies()[i].clone();
1265            final_body.item.properties.position += move_distribution.sample(&mut rng);
1266
1267            assert_relative_eq!(
1268                pairwise_cutoff.delta_energy_one(&microstate_all_pairs, i, &final_body.item),
1269                pairwise_cutoff.delta_energy_one(&microstate_vec_cell, i, &final_body.item),
1270                epsilon = 1e-10
1271            );
1272        }
1273
1274        for _ in 0..N_BODIES / 4 {
1275            let mut new_body = body_template.clone();
1276            new_body.properties.position = square.sample(&mut rng);
1277
1278            assert_relative_eq!(
1279                pairwise_cutoff.delta_energy_insert(&microstate_all_pairs, &new_body),
1280                pairwise_cutoff.delta_energy_insert(&microstate_vec_cell, &new_body),
1281                epsilon = 1e-10
1282            );
1283        }
1284
1285        Ok(())
1286    }
1287}
1288
1289impl<E> MaximumInteractionRange for PairwiseCutoff<E>
1290where
1291    E: MaximumInteractionRange,
1292{
1293    #[inline]
1294    fn maximum_interaction_range(&self) -> f64 {
1295        self.0.maximum_interaction_range()
1296    }
1297}
1298
1299#[cfg(test)]
1300mod test_infinite {
1301    use super::*;
1302    use crate::{
1303        TotalEnergy,
1304        pairwise::{HardShape, HardSphere},
1305    };
1306    use assert2::check;
1307    use hoomd_geometry::shape::{Ellipse, Hypercuboid};
1308    use hoomd_microstate::{
1309        boundary::Closed,
1310        property::{OrientedPoint, Point},
1311    };
1312    use hoomd_vector::{Angle, Cartesian};
1313
1314    use rstest::*;
1315
1316    #[fixture]
1317    fn square() -> Closed<Hypercuboid<2>> {
1318        let cuboid = Hypercuboid {
1319            edge_lengths: [
1320                4.0.try_into()
1321                    .expect("hard-coded constant should be positive"),
1322                4.0.try_into()
1323                    .expect("hard-coded constant should be positive"),
1324            ],
1325        };
1326        Closed(cuboid)
1327    }
1328
1329    mod pairwise_cutoff {
1330        use super::*;
1331
1332        #[test]
1333        fn large_r_cut() -> anyhow::Result<()> {
1334            let mut microstate = Microstate::new();
1335            microstate.extend_bodies([
1336                Body::point(Cartesian::from([0.0, 0.0])),
1337                Body::point(Cartesian::from([0.0, 5.0])),
1338            ])?;
1339
1340            // Ensure that PairwiseCutoff respects the r_cut value set.
1341            let r_cut = 5.0_f64.next_up();
1342            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1343
1344            check!(pairwise_cutoff.total_energy(&microstate) == f64::INFINITY);
1345
1346            let r_cut = 5.0_f64;
1347            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1348
1349            check!(pairwise_cutoff.total_energy(&microstate) == 0.0);
1350
1351            Ok(())
1352        }
1353
1354        #[test]
1355        fn body_exclusion() -> anyhow::Result<()> {
1356            // Ensure that PairwiseCutoff excludes pairs in the same body.
1357            let body_a = Body {
1358                properties: Point::new(Cartesian::from([0.0, 0.0])),
1359                sites: [
1360                    Point::new(Cartesian::from([1.0, 1.0])),
1361                    Point::new(Cartesian::from([1.0, -1.0])),
1362                    Point::new(Cartesian::from([-1.0, 1.0])),
1363                    Point::new(Cartesian::from([-1.0, -1.0])),
1364                ]
1365                .into(),
1366            };
1367            let body_b = Body {
1368                properties: Point::new(Cartesian::from([4.0, 0.0])),
1369                sites: body_a.sites.clone(),
1370            };
1371
1372            let mut microstate = Microstate::new();
1373            microstate.extend_bodies([body_a, body_b])?;
1374
1375            let r_cut = 1.0_f64.next_up();
1376            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1377
1378            check!(pairwise_cutoff.total_energy(&microstate) == 0.0);
1379
1380            let r_cut = 2.0_f64.next_up();
1381            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1382
1383            check!(pairwise_cutoff.total_energy(&microstate) == f64::INFINITY);
1384
1385            Ok(())
1386        }
1387    }
1388
1389    mod delta_energy_one {
1390        use super::*;
1391
1392        #[rstest]
1393        fn site_outside(square: Closed<Hypercuboid<2>>) -> anyhow::Result<()> {
1394            let body = Body {
1395                properties: Point::new(Cartesian::from([0.0, 0.0])),
1396                sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
1397            };
1398            let mut final_body = body.clone();
1399            final_body.properties.position[0] = 1.0;
1400
1401            let microstate = Microstate::builder()
1402                .boundary(square)
1403                .bodies([body])
1404                .try_build()?;
1405
1406            let energy = PairwiseCutoff(HardSphere { diameter: 0.0 });
1407
1408            check!(energy.delta_energy_one(&microstate, 0, &final_body) == f64::INFINITY);
1409
1410            Ok(())
1411        }
1412
1413        #[test]
1414        fn body_exclusion() -> anyhow::Result<()> {
1415            // Ensure that PairwiseCutoff.delta_energy_one excludes pairs in the same body.
1416            let body_a = Body {
1417                properties: Point::new(Cartesian::from([-1.0, 0.0])),
1418                sites: [
1419                    Point::new(Cartesian::from([1.0, 1.0])),
1420                    Point::new(Cartesian::from([1.0, -1.0])),
1421                    Point::new(Cartesian::from([-1.0, 1.0])),
1422                    Point::new(Cartesian::from([-1.0, -1.0])),
1423                ]
1424                .into(),
1425            };
1426            let body_b = Body {
1427                properties: Point::new(Cartesian::from([3.0, 0.0])),
1428                sites: body_a.sites.clone(),
1429            };
1430            let body_a_overlap = Body {
1431                properties: Point::new(Cartesian::from([0.0, 0.0])),
1432                sites: body_a.sites.clone(),
1433            };
1434            let body_a_no_overlap = Body {
1435                properties: Point::new(Cartesian::from([-1.0, -1.0])),
1436                sites: body_a.sites.clone(),
1437            };
1438
1439            let mut microstate = Microstate::new();
1440            microstate.extend_bodies([body_a, body_b])?;
1441
1442            let r_cut = 1.0_f64.next_up();
1443            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1444
1445            // moving body a to the right generates overlaps
1446            check!(
1447                pairwise_cutoff.delta_energy_one(&microstate, 0, &body_a_overlap) == f64::INFINITY
1448            );
1449
1450            // moving body away results in no overlaps
1451            check!(pairwise_cutoff.delta_energy_one(&microstate, 0, &body_a_no_overlap) == 0.0);
1452
1453            Ok(())
1454        }
1455    }
1456
1457    mod delta_energy_insert_remove {
1458        use super::*;
1459
1460        #[rstest]
1461        fn site_outside(square: Closed<Hypercuboid<2>>) -> anyhow::Result<()> {
1462            let body = Body {
1463                properties: Point::new(Cartesian::from([0.0, 0.0])),
1464                sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
1465            };
1466            let mut new_body = body.clone();
1467            new_body.properties.position[0] = 1.0;
1468
1469            let microstate = Microstate::builder()
1470                .boundary(square)
1471                .bodies([body])
1472                .try_build()?;
1473
1474            let energy = PairwiseCutoff(HardSphere { diameter: 0.0 });
1475
1476            check!(energy.delta_energy_insert(&microstate, &new_body) == f64::INFINITY);
1477
1478            Ok(())
1479        }
1480
1481        #[test]
1482        fn body_exclusion() -> anyhow::Result<()> {
1483            // Ensure that PairwiseCutoff.delta_energy_insert excludes pairs in the same body.
1484            let body_a_new = Body {
1485                properties: Point::new(Cartesian::from([0.0, 0.0])),
1486                sites: [
1487                    Point::new(Cartesian::from([1.0, 1.0])),
1488                    Point::new(Cartesian::from([1.0, -1.0])),
1489                    Point::new(Cartesian::from([-1.0, 1.0])),
1490                    Point::new(Cartesian::from([-1.0, -1.0])),
1491                ]
1492                .into(),
1493            };
1494            let body_b = Body {
1495                properties: Point::new(Cartesian::from([3.0, 0.0])),
1496                sites: body_a_new.sites.clone(),
1497            };
1498
1499            let mut microstate = Microstate::new();
1500            microstate.extend_bodies([body_b])?;
1501
1502            let r_cut = 1.0_f64.next_up();
1503            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: r_cut });
1504
1505            check!(pairwise_cutoff.delta_energy_insert(&microstate, &body_a_new) == f64::INFINITY);
1506
1507            microstate.add_body(body_a_new)?;
1508            check!(pairwise_cutoff.delta_energy_remove(&microstate, 1) == 0.0);
1509
1510            Ok(())
1511        }
1512
1513        #[test]
1514        fn hard_shape_initial() -> anyhow::Result<()> {
1515            // Ensure that HardShape always evaluates 0 initial energies.
1516            let a = OrientedPoint {
1517                position: Cartesian::from([0.0, 0.0]),
1518                orientation: Angle::default(),
1519            };
1520            let b = OrientedPoint {
1521                position: Cartesian::from([1.5, 0.0]),
1522                orientation: Angle::default(),
1523            };
1524            let body_a = Body {
1525                properties: a,
1526                sites: [a].into(),
1527            };
1528            let body_b = Body {
1529                properties: b,
1530                sites: [a].into(),
1531            };
1532            let mut microstate = Microstate::new();
1533            microstate.extend_bodies([body_a, body_b.clone()])?;
1534
1535            let ellipse = Ellipse::with_semi_axes([1.0.try_into()?, 2.0.try_into()?]);
1536            let hard_ellipse = PairwiseCutoff(HardShape(ellipse));
1537
1538            // The initial configuration should have infinite energy.
1539            check!(hard_ellipse.total_energy(&microstate) == f64::INFINITY);
1540            check!(hard_ellipse.delta_energy_one(&microstate, 1, &body_b) == f64::INFINITY);
1541
1542            let mut new_body_b = body_b;
1543            new_body_b.properties.position.coordinates = [2.1, 0.0];
1544
1545            // That infinity should be ignored, resulting in a delta E of 0
1546            // when the body is moved into a non-overlapping state.
1547            check!(hard_ellipse.delta_energy_one(&microstate, 1, &new_body_b) == 0.0);
1548
1549            Ok(())
1550        }
1551
1552        #[test]
1553        fn delta_energy_total() -> anyhow::Result<()> {
1554            let mut microstate_0 = Microstate::new();
1555            microstate_0.extend_bodies([
1556                Body::point(Cartesian::from([0.0, 0.0])),
1557                Body::point(Cartesian::from([0.0, 1.125])),
1558            ])?;
1559
1560            let mut microstate_inf = Microstate::new();
1561            microstate_inf.extend_bodies([
1562                Body::point(Cartesian::from([0.0, 0.0])),
1563                Body::point(Cartesian::from([0.0, 0.875])),
1564            ])?;
1565
1566            let pairwise_cutoff = PairwiseCutoff(HardSphere { diameter: 1.0 });
1567
1568            check!(pairwise_cutoff.delta_energy_total(&microstate_0, &microstate_0) == 0.0);
1569            check!(
1570                pairwise_cutoff.delta_energy_total(&microstate_0, &microstate_inf) == f64::INFINITY
1571            );
1572            check!(pairwise_cutoff.delta_energy_total(&microstate_inf, &microstate_0) == 0.0);
1573            check!(
1574                pairwise_cutoff.delta_energy_total(&microstate_inf, &microstate_inf)
1575                    == f64::INFINITY
1576            );
1577
1578            Ok(())
1579        }
1580    }
1581}
hoomd_interaction/pairwise_cutoff.rs

hoomd_interaction/
pairwise_cutoff.rs
