hoomd_interaction/pairwise/approximate_shape_overlap.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 [`ApproximateShapeOverlap`].
5
6use serde::{Deserialize, Serialize};
7
8use super::AnisotropicEnergy;
9use crate::univariate::UnivariateEnergy;
10use hoomd_geometry::IntersectsAt;
11use hoomd_utility::valid::PositiveReal;
12use hoomd_vector::{InnerProduct, Rotate, Rotation};
13
14/// Apply an isotropic potential evaluated at the approximate signed overlap distance.
15///
16/// Use [`ApproximateShapeOverlap`] combined with [`OverlapPenalty`] and
17/// `QuickInsert` to quickly insert new bodies into the microstate or
18/// `QuickCompress` to quickly compress the system to a target density.
19/// In both cases, [`ApproximateShapeOverlap`] will allow partial overlaps
20/// that can be resolved by later trial moves.
21///
22/// [`ApproximateShapeOverlap`] approximates the distance `r` that shape j
23/// must be moved to remove any overlaps. If the shapes are already separate,
24/// `r` is 0. Then it returns `isotropic.energy(-r)`. `resolution` sets the
25/// steps between `r` values.
26///
27/// The overlap distance is *approximate* and expensive to evaluate. Do not
28/// use this potential for equilibration or sampling. It is suitable *only*
29/// for use during a brief initialization phase when `QuickInsert` is adding
30/// bodies or `QuickCompress` is compressing the system.
31///
32/// [`OverlapPenalty`]: crate::univariate::OverlapPenalty
33///
34/// # Example
35///
36/// ```
37/// use hoomd_geometry::{Convex, shape::ConvexPolygon};
38/// use hoomd_interaction::{
39/// pairwise::ApproximateShapeOverlap, univariate::OverlapPenalty,
40/// };
41///
42/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
43/// let approximate_shape_overlap = ApproximateShapeOverlap::new(
44/// Convex(ConvexPolygon::regular(6)),
45/// OverlapPenalty::default(),
46/// 0.01.try_into()?,
47/// );
48/// # Ok(())
49/// # }
50/// ```
51#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
52pub struct ApproximateShapeOverlap<E, A, B = A> {
53 /// The site i's shape.
54 pub shape_i: A,
55 /// The site j's shape.
56 pub shape_j: B,
57 /// Potential to apply as a function of the signed distance between shapes.
58 pub isotropic: E,
59 /// Resolve the separation distance in multiples of this resolution.
60 pub resolution: PositiveReal,
61}
62
63impl<E, A, B, V, R> AnisotropicEnergy<V, R> for ApproximateShapeOverlap<E, A, B>
64where
65 E: UnivariateEnergy,
66 V: InnerProduct,
67 R: Rotation + Rotate<V>,
68 A: IntersectsAt<B, V, R>,
69{
70 /// Evaluate the energy contribution from a pair of sites.
71 ///
72 /// ```
73 /// use hoomd_geometry::{Convex, shape::ConvexPolygon};
74 /// use hoomd_interaction::{
75 /// SitePairEnergy,
76 /// pairwise::{Anisotropic, ApproximateShapeOverlap},
77 /// univariate::OverlapPenalty,
78 /// };
79 /// use hoomd_microstate::property::OrientedPoint;
80 /// use hoomd_vector::{Angle, Cartesian};
81 ///
82 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
83 /// let approximate_shape_overlap = Anisotropic {
84 /// interaction: ApproximateShapeOverlap::new(
85 /// Convex(ConvexPolygon::regular(6)),
86 /// OverlapPenalty::default(),
87 /// 0.01.try_into()?,
88 /// ),
89 /// r_cut: 1.0,
90 /// };
91 ///
92 /// let a = OrientedPoint {
93 /// position: Cartesian::from([0.0, 0.0]),
94 /// orientation: Angle::default(),
95 /// };
96 /// let b = OrientedPoint {
97 /// position: Cartesian::from([0.6, 0.0]),
98 /// orientation: Angle::default(),
99 /// };
100 ///
101 /// let energy = approximate_shape_overlap.site_pair_energy(&a, &b);
102 /// assert!(energy >= 100.0);
103 ///
104 /// let c = OrientedPoint {
105 /// position: Cartesian::from([0.9, 0.0]),
106 /// orientation: Angle::default(),
107 /// };
108 ///
109 /// let energy = approximate_shape_overlap.site_pair_energy(&a, &c);
110 /// assert!(energy >= 0.0);
111 /// assert!(energy < f64::INFINITY);
112 ///
113 /// let d = OrientedPoint {
114 /// position: Cartesian::from([1.0, 0.0]),
115 /// orientation: Angle::default(),
116 /// };
117 ///
118 /// let energy = approximate_shape_overlap.site_pair_energy(&a, &d);
119 /// assert_eq!(energy, 0.0);
120 /// # Ok(())
121 /// # }
122 /// ```
123 #[inline]
124 fn energy(&self, r_ij: &V, o_ij: &R) -> f64 {
125 let r = self.shape_i.approximate_separation_distance(
126 &self.shape_j,
127 r_ij,
128 o_ij,
129 self.resolution,
130 );
131
132 self.isotropic.energy(-r)
133 }
134}
135
136impl<E, G> ApproximateShapeOverlap<E, G> {
137 /// Construct a new approximate shape overlap potential.
138 ///
139 /// `new()` sets both `shape_i` and `shape_j` to `shape`. Use struct
140 /// initialization syntax when you need to set these separately.
141 ///
142 /// # Example
143 ///
144 /// ```
145 /// use hoomd_geometry::{Convex, shape::ConvexPolygon};
146 /// use hoomd_interaction::{
147 /// pairwise::ApproximateShapeOverlap, univariate::OverlapPenalty,
148 /// };
149 ///
150 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
151 /// let approximate_shape_overlap = ApproximateShapeOverlap::new(
152 /// Convex(ConvexPolygon::regular(6)),
153 /// OverlapPenalty::default(),
154 /// 0.01.try_into()?,
155 /// );
156 /// # Ok(())
157 /// # }
158 /// ```
159 #[inline]
160 pub fn new(shape: G, isotropic: E, resolution: PositiveReal) -> Self
161 where
162 G: Clone,
163 {
164 Self {
165 shape_i: shape.clone(),
166 shape_j: shape,
167 isotropic,
168 resolution,
169 }
170 }
171}