hoomd_interaction/univariate/
harmonic_repulsion.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 [`HarmonicRepulsion`]
5
6use serde::{Deserialize, Serialize};
7
8use super::{UnivariateEnergy, UnivariateForce};
9
10/// Repulsive half of a quadratic potential well.
11///
12/// [`HarmonicRepulsion`] is the conservative part of the dissipative particle
13/// dynamics soft potential. The parameter `a` controls the strength of the
14/// potential and `r_cut` sets the distance at which the energy goes to 0.
15///
16///
17/// The potential produce the radially repulsive force between particles as:
18/// ```math
19/// F(r) = \begin{cases}
20/// -\frac{A}{r_\mathrm{cut}} \left( r - r_\mathrm{cut}\right) & r < r_\mathrm{cut} \\
21///
22/// 0 & r \ge r_\mathrm{cut}
23/// \end{cases}
24/// ```
25///
26/// resulting from the potential energy:
27/// ```math
28/// U(r) = \begin{cases}
29/// \frac{1}{2} \frac{A}{r_\mathrm{cut}} \left(r - r_\mathrm{cut}\right)^2 & r \lt r_\mathrm{cut} \\
30///
31/// 0 & r \ge r_\mathrm{cut}
32/// \end{cases}
33/// ```
34///
35/// This potential forms the left half of a harmonic well centered at
36/// $`r = r_\mathrm{cut}`$ with a spring constant $`k = \frac{A}{r_\mathrm{cut}}`$.
37///
38/// Note that this potential has a maximum value of $`A`$ at $`r=0`$ and is
39/// therefore allows particles to completely overlap.
40///
41/// # Examples
42///
43///
44/// ```
45/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
46/// use hoomd_interaction::univariate::{
47///     HarmonicRepulsion, UnivariateEnergy, UnivariateForce,
48/// };
49///
50/// let a = 1.0;
51/// let r_cut = 1.0;
52///
53/// let h_repsulsion = HarmonicRepulsion { a, r_cut };
54/// assert_abs_diff_eq!(h_repsulsion.energy(1.5), 0.0);
55/// assert_relative_eq!(h_repsulsion.energy(0.5), 0.125);
56/// assert_abs_diff_eq!(h_repsulsion.force(0.5), 0.5, epsilon = 1e-12);
57/// ```
58///
59/// The parameters are public fields and may be accessed directly:
60///
61/// ```
62/// use hoomd_interaction::univariate::HarmonicRepulsion;
63///
64/// let mut h_repulsion = HarmonicRepulsion { a: 1.0, r_cut: 1.0 };
65/// h_repulsion.a = 5.0;
66/// h_repulsion.r_cut = 0.75;
67/// ```
68#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
69pub struct HarmonicRepulsion {
70    /// Potential strength $`[\mathrm{energy}] [\mathrm{length}]^{-1}`$.
71    pub a: f64,
72    /// Distance cut-off  $`[\mathrm{length}]`$.
73    pub r_cut: f64,
74}
75
76impl UnivariateEnergy for HarmonicRepulsion {
77    #[inline]
78    fn energy(&self, r: f64) -> f64 {
79        if r < self.r_cut {
80            0.5 * self.a / self.r_cut * (r - self.r_cut).powi(2)
81        } else {
82            0.0
83        }
84    }
85}
86
87impl UnivariateForce for HarmonicRepulsion {
88    #[inline]
89    fn force(&self, r: f64) -> f64 {
90        if r < self.r_cut {
91            self.a / self.r_cut * (self.r_cut - r)
92        } else {
93            0.0
94        }
95    }
96}
97
98#[cfg(test)]
99mod tests {
100    use super::*;
101    use approxim::assert_relative_eq;
102    use rstest::*;
103
104    #[rstest]
105    fn outside_rcut(#[values(1.0, 2.0, 5.0, 10.0)] a: f64, #[values(1.0, 2.0, 3.0)] r_cut: f64) {
106        let h_repulsion = HarmonicRepulsion { a, r_cut };
107
108        assert_eq!(h_repulsion.a, a);
109        assert_eq!(h_repulsion.r_cut, r_cut);
110
111        assert_eq!(h_repulsion.energy(r_cut + 1.0), 0.0);
112        assert_eq!(h_repulsion.force(r_cut + 1.0), 0.0);
113    }
114
115    #[rstest]
116    fn general_case(#[values(1.0, 2.0, 5.0, 10.0)] a: f64, #[values(1.0, 2.0, 3.0)] r_cut: f64) {
117        let r = 0.5;
118        let h_repulsion = HarmonicRepulsion { a, r_cut };
119
120        assert_eq!(h_repulsion.a, a);
121        assert_eq!(h_repulsion.r_cut, r_cut);
122
123        let expected_energy = a * (r_cut - r) - 0.5 * a / r_cut * (r_cut * r_cut - r * r);
124        let expected_force = a * (1.0 - r / r_cut);
125
126        assert_relative_eq!(h_repulsion.energy(r), expected_energy);
127        assert_relative_eq!(h_repulsion.force(r), expected_force);
128    }
129}
hoomd_interaction/univariate/harmonic_repulsion.rs

hoomd_interaction/univariate/
harmonic_repulsion.rs
