hoomd_interaction/external_type.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 `External`
5
6use crate::{
7 DeltaEnergyInsert, DeltaEnergyOne, DeltaEnergyRemove, MaximumInteractionRange,
8 NetSiteForceAndVirial, NetSiteForceVirialAndTorque, SiteEnergy, SiteForceAndVirial,
9 SiteForceVirialAndTorque, TotalEnergy,
10};
11use hoomd_microstate::{Body, Microstate, Transform, boundary::Wrap, property::Position};
12use hoomd_vector::{Outer, Wedge};
13use serde::{Deserialize, Serialize};
14
15/// Interactions between sites and external fields.
16///
17/// An [`External`] newtype wrapping a type that implements [`SiteEnergy`] represents:
18///
19/// ```math
20/// U_\mathrm{total} = \sum_{i=0}^{N-1} U\left( s_i \right)
21/// ```
22/// where $`s_i`$ is the full set of site properties for site i.
23///
24/// An [`External`] newtype wrapping a type that implements [`SiteForceAndVirial`] and/or
25/// [`SiteForceVirialAndTorque`] represents:
26/// ```math
27/// \vec{F}_i = \vec{F}\left(s_i\right)
28/// ```
29/// ```math
30/// \vec{\tau}_i = \vec{\tau}\left(s_i\right)
31/// ```
32/// where $`\vec{F}(s_i)`$ is the force computed by [`SiteForceAndVirial`]
33/// (or [`SiteForceVirialAndTorque`]) and $`\vec{\tau}(s_i)`$ is the torque computed by
34/// [`SiteForceVirialAndTorque`].
35///
36/// A type that implements *both* [`SiteEnergy`] and [`SiteForceAndVirial`]
37/// (or [`SiteForceVirialAndTorque`]) *must* compute forces and torques that are
38/// derivatives of the energy.
39///
40/// Use [`External`] with [`ConstantForce`] or your own custom type that
41/// implements [`SiteEnergy`], [`SiteForceAndVirial`] and/or
42/// [`SiteForceVirialAndTorque`].
43///
44/// [`ConstantForce`]: crate::external::ConstantForce
45///
46/// # Examples
47///
48/// A linear external potential given by a constant force:
49/// ```
50/// use hoomd_interaction::{External, TotalEnergy, external::ConstantForce};
51/// use hoomd_microstate::{Body, Microstate, property::Point};
52/// use hoomd_vector::Cartesian;
53///
54/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
55/// let mut microstate = Microstate::new();
56/// microstate.extend_bodies([
57/// Body::point(Cartesian::from([1.0, 0.0])),
58/// Body::point(Cartesian::from([-1.0, 2.0])),
59/// ])?;
60///
61/// let constant_force = External(ConstantForce {
62/// force: Cartesian::from([0.0, -1.0]),
63/// r_0: Cartesian::default(),
64/// });
65///
66/// let total_energy = constant_force.total_energy(µstate);
67/// assert_eq!(total_energy, 2.0);
68/// # Ok(())
69/// # }
70/// ```
71///
72/// Infinite interaction with a wall:
73/// ```
74/// use hoomd_interaction::{External, SiteEnergy, TotalEnergy};
75/// use hoomd_microstate::{Body, Microstate, property::Point};
76/// use hoomd_vector::Cartesian;
77///
78/// struct Wall;
79///
80/// impl SiteEnergy<Point<Cartesian<2>>> for Wall {
81/// fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
82/// if site_properties.position[1].abs() < 1.0 {
83/// f64::INFINITY
84/// } else {
85/// 0.0
86/// }
87/// }
88///
89/// fn is_only_infinite_or_zero() -> bool {
90/// true
91/// }
92/// }
93///
94/// fn main() -> Result<(), Box<dyn std::error::Error>> {
95/// let mut microstate = Microstate::new();
96/// microstate.extend_bodies([
97/// Body::point(Cartesian::from([1.0, 1.25])),
98/// Body::point(Cartesian::from([-1.0, 2.0])),
99/// ])?;
100///
101/// let wall = External(Wall);
102///
103/// let total_energy = wall.total_energy(µstate);
104/// assert_eq!(total_energy, 0.0);
105/// Ok(())
106/// }
107/// ```
108#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
109pub struct External<E>(pub E);
110
111impl<B, S, X, C, E> TotalEnergy<Microstate<B, S, X, C>> for External<E>
112where
113 E: SiteEnergy<S>,
114{
115 /// Compute the total energy of the microstate contributed by functions of a single site.
116 ///
117 /// The sum over sites differs from HOOMD-blue where external energies are
118 /// evaluated only at the body centers. In general, hoomd-rs interactions apply
119 /// to sites. Use a custom implementation to compute energies over body centers.
120 ///
121 /// # Examples
122 ///
123 /// A linear external potential given by a constant force:
124 /// ```
125 /// use hoomd_interaction::{External, TotalEnergy, external::ConstantForce};
126 /// use hoomd_microstate::{Body, Microstate, property::Point};
127 /// use hoomd_vector::Cartesian;
128 ///
129 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
130 /// let mut microstate = Microstate::new();
131 /// microstate.extend_bodies([
132 /// Body::point(Cartesian::from([1.0, 0.0])),
133 /// Body::point(Cartesian::from([-1.0, 2.0])),
134 /// ])?;
135 ///
136 /// let constant_force = External(ConstantForce {
137 /// force: Cartesian::from([0.0, -1.0]),
138 /// r_0: Cartesian::default(),
139 /// });
140 ///
141 /// let total_energy = constant_force.total_energy(µstate);
142 /// assert_eq!(total_energy, 2.0);
143 /// # Ok(())
144 /// # }
145 /// ```
146 ///
147 /// Infinite interaction with a wall:
148 /// ```
149 /// use hoomd_interaction::{External, SiteEnergy, TotalEnergy};
150 /// use hoomd_microstate::{Body, Microstate, property::Point};
151 /// use hoomd_vector::Cartesian;
152 ///
153 /// struct Wall;
154 ///
155 /// impl SiteEnergy<Point<Cartesian<2>>> for Wall {
156 /// fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
157 /// if site_properties.position[1].abs() < 1.0 {
158 /// f64::INFINITY
159 /// } else {
160 /// 0.0
161 /// }
162 /// }
163 ///
164 /// fn is_only_infinite_or_zero() -> bool {
165 /// true
166 /// }
167 /// }
168 ///
169 /// fn main() -> Result<(), Box<dyn std::error::Error>> {
170 /// let mut microstate = Microstate::new();
171 /// microstate.extend_bodies([
172 /// Body::point(Cartesian::from([1.0, 1.25])),
173 /// Body::point(Cartesian::from([-1.0, 2.0])),
174 /// ])?;
175 ///
176 /// let wall = External(Wall);
177 ///
178 /// let total_energy = wall.total_energy(µstate);
179 /// assert_eq!(total_energy, 0.0);
180 /// Ok(())
181 /// }
182 /// ```
183 #[inline]
184 fn total_energy(&self, microstate: &Microstate<B, S, X, C>) -> f64 {
185 let mut total = 0.0;
186 for site in microstate.sites() {
187 let one = self.0.site_energy(&site.properties);
188 if one == f64::INFINITY {
189 return one;
190 }
191 total += one;
192 }
193
194 total
195 }
196
197 /// Compute the difference in energy between two microstates.
198 ///
199 /// Returns $` E_\mathrm{final} - E_\mathrm{initial} `$.
200 ///
201 /// # Example
202 ///
203 /// ```
204 /// use hoomd_interaction::{External, TotalEnergy, external::ConstantForce};
205 /// use hoomd_microstate::{Body, Microstate, property::Point};
206 /// use hoomd_vector::Cartesian;
207 ///
208 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
209 /// let mut microstate_a = Microstate::new();
210 /// microstate_a.extend_bodies([
211 /// Body::point(Cartesian::from([1.0, 0.0])),
212 /// Body::point(Cartesian::from([-1.0, 2.0])),
213 /// ])?;
214 ///
215 /// let mut microstate_b = Microstate::new();
216 /// microstate_b.extend_bodies([
217 /// Body::point(Cartesian::from([1.0, 1.0])),
218 /// Body::point(Cartesian::from([-1.0, 2.0])),
219 /// ])?;
220 ///
221 /// let constant_force = External(ConstantForce {
222 /// force: Cartesian::from([0.0, -1.0]),
223 /// r_0: Cartesian::default(),
224 /// });
225 ///
226 /// let delta_energy_total =
227 /// constant_force.delta_energy_total(µstate_a, µstate_b);
228 /// assert_eq!(delta_energy_total, 1.0);
229 /// # Ok(())
230 /// # }
231 /// ```
232 #[inline]
233 fn delta_energy_total(
234 &self,
235 initial_microstate: &Microstate<B, S, X, C>,
236 final_microstate: &Microstate<B, S, X, C>,
237 ) -> f64 {
238 let mut energy_final = 0.0;
239 for site in final_microstate.sites() {
240 let one = self.0.site_energy(&site.properties);
241 if one == f64::INFINITY {
242 return one;
243 }
244 energy_final += one;
245 }
246
247 let mut energy_initial = 0.0;
248 if !E::is_only_infinite_or_zero() {
249 for site in initial_microstate.sites() {
250 let one = self.0.site_energy_initial(&site.properties);
251 if one == f64::INFINITY {
252 return -one;
253 }
254 energy_initial += one;
255 }
256 }
257
258 energy_final - energy_initial
259 }
260}
261
262impl<P, B, S, X, C, E> DeltaEnergyOne<B, S, X, C> for External<E>
263where
264 E: SiteEnergy<S>,
265 B: Transform<S>,
266 S: Position<Position = P>,
267 C: Wrap<B> + Wrap<S>,
268{
269 /// Evaluate the change in energy contributed by `External` when a single body is updated.
270 ///
271 /// # Examples
272 ///
273 /// A linear external potential given by a constant force:
274 /// ```
275 /// use hoomd_interaction::{
276 /// DeltaEnergyOne, External, external::ConstantForce,
277 /// };
278 /// use hoomd_microstate::{Body, Microstate, property::Point};
279 /// use hoomd_vector::Cartesian;
280 ///
281 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
282 /// let mut microstate = Microstate::new();
283 /// microstate.add_body(Body::point(Cartesian::from([0.0, 0.0])))?;
284 ///
285 /// let constant_force = External(ConstantForce {
286 /// force: Cartesian::from([0.0, -1.0]),
287 /// r_0: Cartesian::default(),
288 /// });
289 ///
290 /// let delta_energy = constant_force.delta_energy_one(
291 /// µstate,
292 /// 0,
293 /// &Body::point([0.0, -1.0].into()),
294 /// );
295 /// assert_eq!(delta_energy, -1.0);
296 /// # Ok(())
297 /// # }
298 /// ```
299 ///
300 /// Infinite interaction with a wall:
301 /// ```
302 /// use hoomd_interaction::{DeltaEnergyOne, External, SiteEnergy};
303 /// use hoomd_microstate::{Body, Microstate, property::Point};
304 /// use hoomd_vector::Cartesian;
305 ///
306 /// struct Wall;
307 ///
308 /// impl SiteEnergy<Point<Cartesian<2>>> for Wall {
309 /// fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
310 /// if site_properties.position[1].abs() < 1.0 {
311 /// f64::INFINITY
312 /// } else {
313 /// 0.0
314 /// }
315 /// }
316 ///
317 /// fn is_only_infinite_or_zero() -> bool {
318 /// true
319 /// }
320 /// }
321 ///
322 /// fn main() -> Result<(), Box<dyn std::error::Error>> {
323 /// let mut microstate = Microstate::new();
324 /// microstate.extend_bodies([
325 /// Body::point(Cartesian::from([1.0, 1.25])),
326 /// Body::point(Cartesian::from([-1.0, 2.0])),
327 /// ])?;
328 ///
329 /// let wall = External(Wall);
330 ///
331 /// let delta_energy = wall.delta_energy_one(
332 /// µstate,
333 /// 0,
334 /// &Body::point([0.0, -0.5].into()),
335 /// );
336 /// assert_eq!(delta_energy, f64::INFINITY);
337 /// Ok(())
338 /// }
339 /// ```
340 #[inline]
341 fn delta_energy_one(
342 &self,
343 initial_microstate: &Microstate<B, S, X, C>,
344 body_index: usize,
345 final_body: &Body<B, S>,
346 ) -> f64 {
347 let mut energy_final = 0.0;
348 for s in &final_body.sites {
349 match initial_microstate
350 .boundary()
351 .wrap(final_body.properties.transform(s))
352 {
353 Ok(wrapped_site) => {
354 let one = self.0.site_energy(&wrapped_site);
355 if one == f64::INFINITY {
356 return one;
357 }
358 energy_final += one;
359 }
360 Err(_) => return f64::INFINITY,
361 }
362 }
363
364 let energy_initial = if E::is_only_infinite_or_zero() {
365 0.0
366 } else {
367 initial_microstate
368 .iter_body_sites(body_index)
369 .fold(0.0, |total, s| {
370 total + self.0.site_energy_initial(&s.properties)
371 })
372 };
373
374 energy_final - energy_initial
375 }
376}
377
378impl<P, B, S, X, C, E> DeltaEnergyInsert<B, S, X, C> for External<E>
379where
380 E: SiteEnergy<S>,
381 B: Transform<S>,
382 S: Position<Position = P>,
383 C: Wrap<B> + Wrap<S>,
384{
385 /// Evaluate the change in energy contributed by `External` when a single body is inserted.
386 ///
387 /// # Examples
388 ///
389 /// A linear external potential given by a constant force:
390 /// ```
391 /// use hoomd_interaction::{
392 /// DeltaEnergyInsert, External, external::ConstantForce,
393 /// };
394 /// use hoomd_microstate::{Body, Microstate, property::Point};
395 /// use hoomd_vector::Cartesian;
396 ///
397 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
398 /// let mut microstate = Microstate::new();
399 /// microstate.add_body(Body::point(Cartesian::from([0.0, 0.0])))?;
400 ///
401 /// let constant_force = External(ConstantForce {
402 /// force: Cartesian::from([0.0, -1.0]),
403 /// r_0: Cartesian::default(),
404 /// });
405 ///
406 /// let delta_energy = constant_force
407 /// .delta_energy_insert(µstate, &Body::point([0.0, -1.0].into()));
408 /// assert_eq!(delta_energy, -1.0);
409 /// # Ok(())
410 /// # }
411 /// ```
412 ///
413 /// Infinite interaction with a wall:
414 /// ```
415 /// use hoomd_interaction::{DeltaEnergyInsert, External, SiteEnergy};
416 /// use hoomd_microstate::{Body, Microstate, property::Point};
417 /// use hoomd_vector::Cartesian;
418 ///
419 /// struct Wall;
420 ///
421 /// impl SiteEnergy<Point<Cartesian<2>>> for Wall {
422 /// fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
423 /// if site_properties.position[1].abs() < 1.0 {
424 /// f64::INFINITY
425 /// } else {
426 /// 0.0
427 /// }
428 /// }
429 ///
430 /// fn is_only_infinite_or_zero() -> bool {
431 /// true
432 /// }
433 /// }
434 ///
435 /// fn main() -> Result<(), Box<dyn std::error::Error>> {
436 /// let mut microstate = Microstate::new();
437 /// microstate.extend_bodies([
438 /// Body::point(Cartesian::from([1.0, 1.25])),
439 /// Body::point(Cartesian::from([-1.0, 2.0])),
440 /// ])?;
441 ///
442 /// let wall = External(Wall);
443 ///
444 /// let delta_energy = wall
445 /// .delta_energy_insert(µstate, &Body::point([0.0, -0.5].into()));
446 /// assert_eq!(delta_energy, f64::INFINITY);
447 /// Ok(())
448 /// }
449 /// ```
450 #[inline]
451 fn delta_energy_insert(
452 &self,
453 initial_microstate: &Microstate<B, S, X, C>,
454 new_body: &Body<B, S>,
455 ) -> f64 {
456 let mut energy_final = 0.0;
457 for s in &new_body.sites {
458 match initial_microstate
459 .boundary()
460 .wrap(new_body.properties.transform(s))
461 {
462 Ok(wrapped_site) => {
463 let one = self.0.site_energy(&wrapped_site);
464 if one == f64::INFINITY {
465 return one;
466 }
467 energy_final += one;
468 }
469 Err(_) => return f64::INFINITY,
470 }
471 }
472
473 energy_final
474 }
475}
476
477impl<B, S, X, C, E> DeltaEnergyRemove<B, S, X, C> for External<E>
478where
479 E: SiteEnergy<S>,
480{
481 /// Evaluate the change in energy contributed by `External` when a single body is removed.
482 ///
483 /// # Examples
484 ///
485 /// A linear external potential given by a constant force:
486 /// ```
487 /// use hoomd_interaction::{
488 /// DeltaEnergyRemove, External, external::ConstantForce,
489 /// };
490 /// use hoomd_microstate::{Body, Microstate, property::Point};
491 /// use hoomd_vector::Cartesian;
492 ///
493 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
494 /// let mut microstate = Microstate::new();
495 /// microstate.add_body(Body::point(Cartesian::from([0.0, 1.0])))?;
496 ///
497 /// let constant_force = External(ConstantForce {
498 /// force: Cartesian::from([0.0, -1.0]),
499 /// r_0: Cartesian::default(),
500 /// });
501 ///
502 /// let delta_energy = constant_force.delta_energy_remove(µstate, 0);
503 /// assert_eq!(delta_energy, -1.0);
504 /// # Ok(())
505 /// # }
506 /// ```
507 ///
508 /// Infinite interaction with a wall:
509 /// ```
510 /// use hoomd_interaction::{DeltaEnergyRemove, External, SiteEnergy};
511 /// use hoomd_microstate::{Body, Microstate, property::Point};
512 /// use hoomd_vector::Cartesian;
513 ///
514 /// struct Wall;
515 ///
516 /// impl SiteEnergy<Point<Cartesian<2>>> for Wall {
517 /// fn site_energy(&self, site_properties: &Point<Cartesian<2>>) -> f64 {
518 /// if site_properties.position[1].abs() < 1.0 {
519 /// f64::INFINITY
520 /// } else {
521 /// 0.0
522 /// }
523 /// }
524 ///
525 /// fn is_only_infinite_or_zero() -> bool {
526 /// true
527 /// }
528 /// }
529 ///
530 /// fn main() -> Result<(), Box<dyn std::error::Error>> {
531 /// let mut microstate = Microstate::new();
532 /// microstate.extend_bodies([
533 /// Body::point(Cartesian::from([1.0, 1.25])),
534 /// Body::point(Cartesian::from([-1.0, 2.0])),
535 /// ])?;
536 ///
537 /// let wall = External(Wall);
538 ///
539 /// let delta_energy = wall.delta_energy_remove(µstate, 0);
540 /// assert_eq!(delta_energy, 0.0);
541 /// Ok(())
542 /// }
543 /// ```
544 #[inline]
545 fn delta_energy_remove(
546 &self,
547 initial_microstate: &Microstate<B, S, X, C>,
548 body_index: usize,
549 ) -> f64 {
550 if E::is_only_infinite_or_zero() {
551 return 0.0;
552 }
553
554 let energy_initial = initial_microstate
555 .iter_body_sites(body_index)
556 .fold(0.0, |total, s| {
557 total + self.0.site_energy_initial(&s.properties)
558 });
559
560 -energy_initial
561 }
562}
563
564impl<V, B, S, X, C, E> NetSiteForceAndVirial<B, S, X, C> for External<E>
565where
566 V: Outer,
567 S: Position<Position = V>,
568 E: SiteForceAndVirial<S, Force = V>,
569{
570 type Force = V;
571
572 /// Compute the net force and virial on a given site.
573 ///
574 /// # Example
575 ///
576 /// ```
577 /// use hoomd_interaction::{
578 /// External, NetSiteForceAndVirial, external::ConstantForce,
579 /// };
580 /// use hoomd_linear_algebra::matrix::Matrix;
581 /// use hoomd_microstate::{Body, Microstate, property::Point};
582 /// use hoomd_vector::Cartesian;
583 ///
584 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
585 /// let mut microstate = Microstate::new();
586 /// microstate.add_body(Body::point(Cartesian::from([2.0, 0.0])))?;
587 ///
588 /// let constant_force = External(ConstantForce {
589 /// force: Cartesian::from([0.0, -1.0]),
590 /// r_0: Cartesian::default(),
591 /// });
592 ///
593 /// let (force, virial) =
594 /// constant_force.net_site_force_and_virial(µstate, 0);
595 /// assert_eq!(force, [0.0, -1.0].into());
596 /// assert_eq!(
597 /// virial,
598 /// Matrix {
599 /// rows: [[0.0, 0.0], [-2.0, 0.0]]
600 /// }
601 /// );
602 /// # Ok(())
603 /// # }
604 /// ```
605 #[inline]
606 fn net_site_force_and_virial(
607 &self,
608 microstate: &Microstate<B, S, X, C>,
609 site_index: usize,
610 ) -> (V, V::Tensor) {
611 let site = µstate.sites()[site_index];
612 self.0.site_force_and_virial(&site.properties)
613 }
614}
615
616impl<V, B, S, X, C, E> NetSiteForceVirialAndTorque<B, S, X, C> for External<E>
617where
618 V: Wedge + Outer,
619 S: Position<Position = V>,
620 E: SiteForceVirialAndTorque<S, Force = V>,
621{
622 type Force = V;
623
624 /// Compute the net force, virial, and torque on a given site.
625 ///
626 /// # Example
627 ///
628 /// ```
629 /// use hoomd_interaction::{
630 /// External, NetSiteForceVirialAndTorque, external::ConstantForce,
631 /// };
632 /// use hoomd_linear_algebra::matrix::Matrix;
633 /// use hoomd_microstate::{Body, Microstate, property::Point};
634 /// use hoomd_vector::Cartesian;
635 ///
636 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
637 /// let mut microstate = Microstate::new();
638 /// microstate.add_body(Body::point(Cartesian::from([2.0, 0.0])))?;
639 ///
640 /// let constant_force = External(ConstantForce {
641 /// force: Cartesian::from([0.0, -1.0]),
642 /// r_0: Cartesian::default(),
643 /// });
644 ///
645 /// let (force, virial, torque) =
646 /// constant_force.net_site_force_virial_and_torque(µstate, 0);
647 /// assert_eq!(force, [0.0, -1.0].into());
648 /// assert_eq!(
649 /// virial,
650 /// Matrix {
651 /// rows: [[0.0, 0.0], [-2.0, 0.0]]
652 /// }
653 /// );
654 /// assert_eq!(torque, 0.0);
655 /// # Ok(())
656 /// # }
657 /// ```
658 #[inline]
659 fn net_site_force_virial_and_torque(
660 &self,
661 microstate: &Microstate<B, S, X, C>,
662 site_index: usize,
663 ) -> (V, V::Tensor, V::Bivector) {
664 let site = µstate.sites()[site_index];
665 let (force, virial, torque) = self.0.site_force_virial_and_torque(&site.properties);
666 (force, virial, torque)
667 }
668}
669
670impl<E> MaximumInteractionRange for External<E> {
671 #[inline]
672 fn maximum_interaction_range(&self) -> f64 {
673 // External interactions are not applied between pairs of particles.
674 0.0
675 }
676}
677
678#[cfg(test)]
679mod test_finite {
680 use super::*;
681
682 use crate::external::ConstantForce;
683 use assert2::check;
684 use hoomd_geometry::shape::Rectangle;
685 use hoomd_microstate::{
686 Body, Microstate,
687 boundary::{Closed, Open},
688 property::{Point, Position},
689 };
690 use hoomd_vector::Cartesian;
691 use rstest::*;
692
693 struct TestSE;
694
695 impl<S> SiteEnergy<S> for TestSE
696 where
697 S: Position<Position = Cartesian<2>>,
698 {
699 fn site_energy(&self, site_properties: &S) -> f64 {
700 site_properties.position()[0] + site_properties.position()[1]
701 }
702 }
703
704 mod site_energy {
705 use super::*;
706 use hoomd_microstate::SiteKey;
707 use hoomd_spatial::AllPairs;
708
709 #[fixture]
710 fn microstate()
711 -> Microstate<Point<Cartesian<2>>, Point<Cartesian<2>>, AllPairs<SiteKey>, Open> {
712 let mut microstate = Microstate::new();
713 microstate
714 .extend_bodies([
715 Body::point(Cartesian::from([1.0, 0.0])),
716 Body::point(Cartesian::from([-1.0, 3.0])),
717 ])
718 .expect("hard-coded bodies should be in the boundary");
719 microstate
720 }
721
722 #[rstest]
723 fn single_total(
724 microstate: Microstate<
725 Point<Cartesian<2>>,
726 Point<Cartesian<2>>,
727 AllPairs<SiteKey>,
728 Open,
729 >,
730 ) {
731 let test_se = TestSE;
732 let single = External(test_se);
733
734 check!(single.total_energy(µstate) == 3.0);
735 }
736
737 #[rstest]
738 fn single_site(
739 microstate: Microstate<
740 Point<Cartesian<2>>,
741 Point<Cartesian<2>>,
742 AllPairs<SiteKey>,
743 Open,
744 >,
745 ) {
746 let test_se = TestSE;
747 let single = External(test_se);
748
749 check!(single.0.site_energy(µstate.sites()[0].properties) == 1.0);
750 check!(single.0.site_energy(µstate.sites()[1].properties) == 2.0);
751 }
752 }
753 mod delta_energy {
754 use super::*;
755
756 struct Zero;
757
758 impl SiteEnergy<Point<Cartesian<2>>> for Zero {
759 fn site_energy(&self, _site_properties: &Point<Cartesian<2>>) -> f64 {
760 0.0
761 }
762 }
763
764 #[test]
765 fn site_outside() {
766 let cuboid = Rectangle::with_equal_edges(
767 4.0.try_into()
768 .expect("hard-coded constant should be positive"),
769 );
770 let square = Closed(cuboid);
771
772 let body = Body {
773 properties: Point::new(Cartesian::from([0.0, 0.0])),
774 sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
775 };
776 let mut final_body = body.clone();
777 final_body.properties.position[0] = 1.0;
778
779 let microstate = Microstate::builder()
780 .boundary(square)
781 .bodies([body])
782 .try_build()
783 .expect("the hard-coded bodies should be in the boundary");
784
785 let energy = External(Zero);
786
787 check!(energy.delta_energy_one(µstate, 0, &final_body) == f64::INFINITY);
788 check!(energy.delta_energy_insert(µstate, &final_body) == f64::INFINITY);
789 }
790
791 #[test]
792 fn delta_energy() -> anyhow::Result<()> {
793 let cuboid = Rectangle::with_equal_edges(
794 4.0.try_into()
795 .expect("hard-coded constant should be positive"),
796 );
797 let square = Closed(cuboid);
798
799 let body = Body {
800 properties: Point::new(Cartesian::from([0.0, 0.0])),
801 sites: [Point::new(Cartesian::from([0.0, 0.0]))].into(),
802 };
803 let mut final_body = body.clone();
804 final_body.properties.position[1] = 0.5;
805
806 let microstate = Microstate::builder()
807 .boundary(square)
808 .bodies([body])
809 .try_build()
810 .expect("the hard-coded bodies should be in the boundary");
811
812 let energy = External(ConstantForce {
813 r_0: [0.0, -1.0].into(),
814 force: [0.0, -4.0].into(),
815 });
816
817 check!(energy.delta_energy_one(µstate, 0, &final_body) == 2.0);
818 check!(energy.delta_energy_insert(µstate, &final_body) == 6.0);
819 check!(energy.delta_energy_remove(µstate, 0) == -4.0);
820
821 let mut microstate_final = microstate.clone();
822 microstate_final.update_body_properties(0, final_body.properties)?;
823
824 check!(energy.delta_energy_total(µstate, µstate_final) == 2.0);
825
826 Ok(())
827 }
828 }
829}
830
831#[cfg(test)]
832mod test_infinite {
833 use super::*;
834 use assert2::check;
835 use hoomd_geometry::shape::Rectangle;
836 use hoomd_microstate::{
837 Body, Microstate,
838 boundary::{Closed, Open},
839 property::{Point, Position},
840 };
841 use hoomd_vector::Cartesian;
842 use rstest::*;
843
844 struct TestSO;
845
846 impl<S> SiteEnergy<S> for TestSO
847 where
848 S: Position<Position = Cartesian<2>>,
849 {
850 fn site_energy(&self, site_properties: &S) -> f64 {
851 if site_properties.position()[1].abs() < 1.0 {
852 f64::INFINITY
853 } else {
854 0.0
855 }
856 }
857
858 fn is_only_infinite_or_zero() -> bool {
859 true
860 }
861 }
862
863 mod site_energy {
864 use super::*;
865 use hoomd_microstate::SiteKey;
866 use hoomd_spatial::AllPairs;
867
868 #[fixture]
869 fn microstate()
870 -> Microstate<Point<Cartesian<2>>, Point<Cartesian<2>>, AllPairs<SiteKey>, Open> {
871 let mut microstate = Microstate::new();
872 microstate
873 .extend_bodies([
874 Body::point(Cartesian::from([1.0, -2.0])),
875 Body::point(Cartesian::from([-1.0, 3.0])),
876 ])
877 .expect("hard-coded bodies should be in the boundary");
878 microstate
879 }
880
881 #[fixture]
882 fn overlapping_microstate()
883 -> Microstate<Point<Cartesian<2>>, Point<Cartesian<2>>, AllPairs<SiteKey>, Open> {
884 let mut microstate = Microstate::new();
885 microstate
886 .extend_bodies([
887 Body::point(Cartesian::from([1.0, 0.75])),
888 Body::point(Cartesian::from([-1.0, 3.0])),
889 ])
890 .expect("hard-coded bodies should be in the boundary");
891 microstate
892 }
893
894 #[rstest]
895 fn single_total_0(
896 microstate: Microstate<
897 Point<Cartesian<2>>,
898 Point<Cartesian<2>>,
899 AllPairs<SiteKey>,
900 Open,
901 >,
902 ) {
903 let single = External(TestSO);
904
905 check!(single.total_energy(µstate) == 0.0);
906 }
907
908 #[rstest]
909 fn single_total_inf(
910 overlapping_microstate: Microstate<
911 Point<Cartesian<2>>,
912 Point<Cartesian<2>>,
913 AllPairs<SiteKey>,
914 Open,
915 >,
916 ) {
917 let single = External(TestSO);
918
919 check!(single.total_energy(&overlapping_microstate) == f64::INFINITY);
920 }
921
922 #[rstest]
923 fn single_site_0(
924 microstate: Microstate<
925 Point<Cartesian<2>>,
926 Point<Cartesian<2>>,
927 AllPairs<SiteKey>,
928 Open,
929 >,
930 ) {
931 let single = External(TestSO);
932
933 check!(single.0.site_energy(µstate.sites()[0].properties) == 0.0);
934 check!(single.0.site_energy(µstate.sites()[1].properties) == 0.0);
935 }
936
937 #[rstest]
938 fn single_site_inf(
939 overlapping_microstate: Microstate<
940 Point<Cartesian<2>>,
941 Point<Cartesian<2>>,
942 AllPairs<SiteKey>,
943 Open,
944 >,
945 ) {
946 let single = External(TestSO);
947
948 check!(
949 single
950 .0
951 .site_energy(&overlapping_microstate.sites()[0].properties)
952 == f64::INFINITY
953 );
954 check!(
955 single
956 .0
957 .site_energy(&overlapping_microstate.sites()[1].properties)
958 == 0.0
959 );
960 }
961 }
962 mod delta_energy {
963 use super::*;
964
965 struct Zero;
966
967 impl SiteEnergy<Point<Cartesian<2>>> for Zero {
968 fn site_energy(&self, _site_properties: &Point<Cartesian<2>>) -> f64 {
969 0.0
970 }
971
972 fn is_only_infinite_or_zero() -> bool {
973 true
974 }
975 }
976
977 #[test]
978 fn site_outside() {
979 let cuboid = Rectangle::with_equal_edges(
980 4.0.try_into()
981 .expect("hard-coded constant should be positive"),
982 );
983 let square = Closed(cuboid);
984
985 let body = Body {
986 properties: Point::new(Cartesian::from([0.0, 0.0])),
987 sites: [Point::new(Cartesian::from([1.0, 0.0]))].into(),
988 };
989 let mut final_body = body.clone();
990 final_body.properties.position[0] = 1.0;
991
992 let microstate = Microstate::builder()
993 .boundary(square)
994 .bodies([body])
995 .try_build()
996 .expect("the hard-coded bodies should be in the boundary");
997
998 let energy = External(Zero);
999
1000 check!(energy.delta_energy_one(µstate, 0, &final_body) == f64::INFINITY);
1001 check!(energy.delta_energy_insert(µstate, &final_body) == f64::INFINITY);
1002 }
1003
1004 #[test]
1005 fn delta_energy() -> anyhow::Result<()> {
1006 let cuboid = Rectangle::with_equal_edges(
1007 4.0.try_into()
1008 .expect("hard-coded constant should be positive"),
1009 );
1010 let square = Closed(cuboid);
1011
1012 let body = Body {
1013 properties: Point::new(Cartesian::from([1.5, 1.5])),
1014 sites: [Point::new(Cartesian::from([0.0, 0.0]))].into(),
1015 };
1016 let mut final_body_inf = body.clone();
1017 final_body_inf.properties.position[1] = 0.5;
1018
1019 let mut final_body_0 = body.clone();
1020 final_body_0.properties.position[1] = -1.5;
1021
1022 let microstate = Microstate::builder()
1023 .boundary(square)
1024 .bodies([body])
1025 .try_build()
1026 .expect("the hard-coded bodies should be in the boundary");
1027
1028 let energy = External(TestSO);
1029
1030 check!(energy.delta_energy_one(µstate, 0, &final_body_0) == 0.0);
1031 check!(energy.delta_energy_one(µstate, 0, &final_body_inf) == f64::INFINITY);
1032 check!(energy.delta_energy_insert(µstate, &final_body_0) == 0.0);
1033 check!(energy.delta_energy_insert(µstate, &final_body_inf) == f64::INFINITY);
1034 check!(energy.delta_energy_remove(µstate, 0) == 0.0);
1035
1036 let mut microstate_inf = microstate.clone();
1037 microstate_inf.update_body_properties(0, final_body_inf.properties)?;
1038
1039 let mut microstate_0 = microstate.clone();
1040 microstate_0.update_body_properties(0, final_body_0.properties)?;
1041
1042 check!(energy.delta_energy_total(µstate_0, µstate_0) == 0.0);
1043 check!(energy.delta_energy_total(µstate_0, µstate_inf) == f64::INFINITY);
1044 check!(energy.delta_energy_total(µstate_inf, µstate_0) == 0.0);
1045 check!(energy.delta_energy_total(µstate_inf, µstate_inf) == f64::INFINITY);
1046
1047 Ok(())
1048 }
1049 }
1050}