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