hoomd_interaction/univariate/
overlap_penalty.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 `OverlapPenalty`.
5
6use serde::{Deserialize, Serialize};
7
8use super::UnivariateEnergy;
9
10/// Monotonically non-decreasing potential to push sites apart (*not differentiable*).
11///
12/// [`OverlapPenalty`] is specifically designed to work with the `QuickInsert`
13/// and `QuickCompress` algorithms to quickly prepare states with non-overlapping
14/// particles. Combine it with [`ApproximateShapeOverlap`] to compute an energy that
15/// penalizes hard particle overlaps.
16///
17/// The potential has three regions:
18/// ```math
19/// U(r) = \begin{cases}
20/// \infty & r < -d_\mathrm{max} \\
21/// \frac{1}{2} kr^2 + \varepsilon_\mathrm{shoulder} & r < 0 \\
22/// 0 & r \ge 0
23/// \end{cases}
24/// ```
25/// The first region describes a completely hard interaction when sites overlap
26/// too far. This prevents `QuickInsert` from creating too much strain with an
27/// insertion. The second part applies a harmonic potential that allows trial moves
28/// to gradually resolve overlaps. In the third region, sites are allowed to move
29/// freely when not overlapping. The shoulder potential prevents trial moves from
30/// creating new overlaps.
31///
32/// [`ApproximateShapeOverlap`]: crate::pairwise::ApproximateShapeOverlap
33///
34/// # Example
35///
36/// ```
37/// use hoomd_interaction::univariate::OverlapPenalty;
38///
39/// let overlap_penalty = OverlapPenalty::default();
40/// ```
41#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
42pub struct OverlapPenalty {
43    /// Spring stiffness $`[\mathrm{energy}] [\mathrm{length}]^{-2}`$.
44    pub k: f64,
45
46    /// The largest overlap distance to allow  $`[\mathrm{length}]`$.
47    pub maximum_allowed_overlap: f64,
48
49    /// Height of the potential as $`r`$ approaches 0 from the left $`[\mathrm{energy}]`$
50    pub epsilon_shoulder: f64,
51}
52
53impl Default for OverlapPenalty {
54    /// Default overlap penalty parameters.
55    ///
56    /// The default values are tuned for use with `QuickInsert` and `QuickCompress`
57    /// applied to systems of spherical particles with diameter approximately 1.
58    ///
59    /// * $`k = 10,000`$
60    /// * $`d_\mathrm{max} = 0.2`$
61    /// * $`\varepsilon_\mathrm{shoulder} = 100`$
62    ///
63    /// Call [`OverlapPenalty::scaled_default`] to initialize with values scaled
64    /// for use with non-unit diameter sites.
65    ///
66    /// # Example
67    ///
68    /// ```
69    /// use hoomd_interaction::univariate::OverlapPenalty;
70    ///
71    /// let overlap_penalty = OverlapPenalty::default();
72    /// ```
73    #[inline]
74    fn default() -> Self {
75        Self {
76            k: 10_000.0,
77            maximum_allowed_overlap: 0.2,
78            epsilon_shoulder: 100.0,
79        }
80    }
81}
82
83impl OverlapPenalty {
84    /// Default overlap penalty parameters for a given diameter.
85    ///
86    /// Construct an [`OverlapPenalty`] with default parameters scaled to suit
87    /// a site with the given diameter.
88    ///
89    /// # Example
90    ///
91    /// ```
92    /// use hoomd_interaction::univariate::OverlapPenalty;
93    ///
94    /// let overlap_penalty = OverlapPenalty::scaled_default(2.0);
95    ///
96    /// assert_eq!(overlap_penalty.maximum_allowed_overlap, 0.4);
97    /// ```
98    #[must_use]
99    #[inline]
100    pub fn scaled_default(diameter: f64) -> Self {
101        let mut overlap_penalty = Self::default();
102        overlap_penalty.maximum_allowed_overlap *= diameter;
103        overlap_penalty
104    }
105}
106
107impl UnivariateEnergy for OverlapPenalty {
108    #[inline]
109    fn energy(&self, r: f64) -> f64 {
110        match r {
111            x if x < -self.maximum_allowed_overlap => f64::INFINITY,
112            x if x < 0.0 => 0.5 * self.k * r.powi(2) + self.epsilon_shoulder,
113            _ => 0.0,
114        }
115    }
116}
117
118#[cfg(test)]
119mod tests {
120    use super::*;
121
122    #[test]
123    fn shoulder() {
124        let k = 0.0;
125        let epsilon_shoulder = 15.0;
126        let maximum_allowed_overlap = 1.0;
127
128        let overlap_penalty = OverlapPenalty {
129            k,
130            maximum_allowed_overlap,
131            epsilon_shoulder,
132        };
133
134        assert_eq!(overlap_penalty.energy(1.0), 0.0);
135        assert_eq!(overlap_penalty.energy(0.0), 0.0);
136        assert_eq!(
137            overlap_penalty.energy(0.0_f64.next_down()),
138            epsilon_shoulder
139        );
140        assert_eq!(overlap_penalty.energy(-0.5), epsilon_shoulder);
141    }
142
143    #[test]
144    fn core() {
145        let k = 0.0;
146        let epsilon_shoulder = 0.0;
147        let maximum_allowed_overlap = 0.5;
148
149        let overlap_penalty = OverlapPenalty {
150            k,
151            maximum_allowed_overlap,
152            epsilon_shoulder,
153        };
154
155        assert_eq!(overlap_penalty.energy(1.0), 0.0);
156        assert_eq!(overlap_penalty.energy(0.0), 0.0);
157        assert_eq!(overlap_penalty.energy(0.0_f64.next_down()), 0.0);
158        assert_eq!(overlap_penalty.energy(-0.5), 0.0);
159        assert_eq!(
160            overlap_penalty.energy((-0.5_f64).next_down()),
161            f64::INFINITY
162        );
163        assert_eq!(overlap_penalty.energy(-1.0), f64::INFINITY);
164    }
165
166    #[test]
167    fn harmonic() {
168        let k = 6.0;
169        let epsilon_shoulder = 0.0;
170        let maximum_allowed_overlap = 10.0;
171
172        let overlap_penalty = OverlapPenalty {
173            k,
174            maximum_allowed_overlap,
175            epsilon_shoulder,
176        };
177
178        assert_eq!(overlap_penalty.energy(1.0), 0.0);
179        assert_eq!(overlap_penalty.energy(0.0), 0.0);
180        assert_eq!(overlap_penalty.energy(-0.125), 0.5 * k * 0.125_f64.powi(2));
181    }
182}
hoomd_interaction/univariate/overlap_penalty.rs

hoomd_interaction/univariate/
overlap_penalty.rs
