hoomd_interaction/pairwise/
hard_shape.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 `HardShape`
5
6use serde::{Deserialize, Serialize};
7
8use crate::{MaximumInteractionRange, SitePairEnergy};
9use hoomd_geometry::{BoundingSphereRadius, IntersectsAtGlobal};
10use hoomd_microstate::property::{Orientation, Position};
11use hoomd_vector::Metric;
12
13/// Infinite energy when sites overlap, 0 when they don't (*not differentiable*).
14///
15/// [`HardShape`] represents each site with a hard orientable shape.
16///
17/// The generic type names are:
18/// * `G`: The [`shape`](hoomd_geometry::shape) type.
19///
20/// # Example
21///
22/// ```
23/// use hoomd_geometry::{Convex, shape::Rectangle};
24/// use hoomd_interaction::pairwise::HardShape;
25///
26/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
27/// let square = Rectangle::with_equal_edges(1.0.try_into()?);
28/// let hard_shape = HardShape(Convex(square));
29/// # Ok(())
30/// # }
31/// ```
32#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
33pub struct HardShape<G>(pub G);
34
35impl<G, S, R, P> SitePairEnergy<S> for HardShape<G>
36where
37    S: Position<Position = P> + Orientation<Rotation = R>,
38    G: IntersectsAtGlobal<G, P, R>,
39{
40    /// Compute the energy contribution from a pair of sites.
41    ///
42    /// A pair of hard shapes contributes an infinite energy when they overlap,
43    /// and zero when they do not.
44    ///
45    /// # Example
46    ///
47    /// ```
48    /// use hoomd_geometry::{Convex, shape::Rectangle};
49    /// use hoomd_interaction::{SitePairEnergy, pairwise::HardShape};
50    /// use hoomd_microstate::property::OrientedPoint;
51    /// use hoomd_vector::{Angle, Cartesian};
52    /// use std::f64::consts::PI;
53    ///
54    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
55    /// let square = Rectangle::with_equal_edges(1.0.try_into()?);
56    /// let hard_shape = HardShape(Convex(square));
57    ///
58    /// let a = OrientedPoint {
59    ///     position: Cartesian::from([1.0, -1.0]),
60    ///     orientation: Angle::from(PI / 2.0),
61    /// };
62    /// let b = OrientedPoint {
63    ///     position: Cartesian::from([2.0, 0.0]),
64    ///     orientation: Angle::from(PI / 4.0),
65    /// };
66    ///
67    /// assert_eq!(hard_shape.site_pair_energy(&a, &b), 0.0);
68    ///
69    /// let c = OrientedPoint {
70    ///     position: Cartesian::from([1.5, -0.5]),
71    ///     orientation: Angle::from(PI / 4.0),
72    /// };
73    ///
74    /// assert_eq!(hard_shape.site_pair_energy(&a, &c), f64::INFINITY);
75    /// # Ok(())
76    /// # }
77    /// ```
78    #[inline]
79    fn site_pair_energy(&self, site_properties_i: &S, site_properties_j: &S) -> f64 {
80        // Use the global form of `intersects_at` so that the circumsphere early
81        // rejection test can be performed before the expensive transformation
82        // into the local coordinate system.
83        if self.0.intersects_at_global(
84            &self.0,
85            site_properties_i.position(),
86            site_properties_i.orientation(),
87            site_properties_j.position(),
88            site_properties_j.orientation(),
89        ) {
90            f64::INFINITY
91        } else {
92            0.0
93        }
94    }
95
96    /// Evaluate the energy contribution from a pair of sites *in the initial state*.
97    ///
98    /// Hard shapes are assumed to be non-overlapping in the initial state.
99    /// This method always returns zero.
100    #[inline]
101    fn site_pair_energy_initial(&self, _site_properties_i: &S, _site_properties_j: &S) -> f64 {
102        0.0
103    }
104
105    #[inline]
106    fn is_only_infinite_or_zero() -> bool {
107        true
108    }
109}
110
111impl<G> MaximumInteractionRange for HardShape<G>
112where
113    G: BoundingSphereRadius,
114{
115    #[inline]
116    fn maximum_interaction_range(&self) -> f64 {
117        self.0.bounding_sphere_radius().get() * 2.0
118    }
119}
120
121/// Model infinitely hard spheres when used with [`PairwiseCutoff`] (*not differentiable*).
122///
123/// [`HardSphere`] represents each site as a hard sphere with a diameter given
124/// by the `diameter` field.
125///
126/// [`PairwiseCutoff`]: crate::PairwiseCutoff
127#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
128pub struct HardSphere {
129    /// The sphere's diameter.
130    pub diameter: f64,
131}
132
133impl<S, P> SitePairEnergy<S> for HardSphere
134where
135    S: Position<Position = P>,
136    P: Metric,
137{
138    /// Compute the energy contribution from a pair of sites.
139    ///
140    /// The interaction energy is infinite when two spheres overlap and zero
141    /// when they do not.
142    #[inline]
143    fn site_pair_energy(&self, site_properties_i: &S, site_properties_j: &S) -> f64 {
144        if site_properties_i
145            .position()
146            .distance_squared(site_properties_j.position())
147            < self.diameter.powi(2)
148        {
149            f64::INFINITY
150        } else {
151            0.0
152        }
153    }
154
155    /// Evaluate the energy contribution from a pair of sites *in the initial state*.
156    ///
157    /// Hard shapes are assumed to be non-overlapping in the initial state.
158    /// This implementation always returns zero.
159    #[inline]
160    fn site_pair_energy_initial(&self, _site_properties_i: &S, _site_properties_j: &S) -> f64 {
161        0.0
162    }
163
164    #[inline]
165    fn is_only_infinite_or_zero() -> bool {
166        true
167    }
168}
169
170impl MaximumInteractionRange for HardSphere {
171    #[inline]
172    fn maximum_interaction_range(&self) -> f64 {
173        self.diameter
174    }
175}
hoomd_interaction/pairwise/hard_shape.rs

hoomd_interaction/pairwise/
hard_shape.rs
