Skip to main content

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}