hoomd_interaction/univariate/
harmonic.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 [`Harmonic`]
5
6use serde::{Deserialize, Serialize};
7
8use super::{UnivariateEnergy, UnivariateForce};
9
10/// Quadratic potential well centered at a given separation distance.
11///
12/// ```math
13/// U = \frac{1}{2} k (r - r_0)^2
14/// ```
15///
16/// Compute the harmonic potential and force as a function of `r` with
17/// equilibrium spring length `r_0`.
18///
19/// # Examples
20///
21/// ```
22/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
23/// use hoomd_interaction::univariate::{
24///     Harmonic, UnivariateEnergy, UnivariateForce,
25/// };
26///
27/// let k = 2.0;
28/// let r_0 = 0.0;
29///
30/// let harmonic = Harmonic { k, r_0 };
31/// assert_abs_diff_eq!(harmonic.energy(0.0), 0.0);
32/// assert_relative_eq!(harmonic.energy(1.0), 1.0);
33/// assert_abs_diff_eq!(harmonic.force(1.0), -2.0, epsilon = 1e-12);
34/// ```
35///
36/// The parameters are public fields and may be accessed directly:
37///
38/// ```
39/// use hoomd_interaction::univariate::Harmonic;
40///
41/// let mut harmonic = Harmonic { k: 1.0, r_0: 0.0 };
42/// harmonic.k = 5.0;
43/// harmonic.r_0 = 1.0;
44/// ```
45#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
46pub struct Harmonic {
47    /// Spring constant $`[\mathrm{energy}] [\mathrm{length}]^{-2}`$.
48    pub k: f64,
49    /// Equilibrium spring length $`[\mathrm{length}]`$.
50    pub r_0: f64,
51}
52
53impl UnivariateEnergy for Harmonic {
54    #[inline]
55    fn energy(&self, r: f64) -> f64 {
56        0.5 * self.k * (r - self.r_0) * (r - self.r_0)
57    }
58}
59
60impl UnivariateForce for Harmonic {
61    #[inline]
62    fn force(&self, r: f64) -> f64 {
63        -self.k * (r - self.r_0)
64    }
65}
66
67#[cfg(test)]
68mod tests {
69    use super::*;
70    use approxim::assert_relative_eq;
71    use rstest::*;
72
73    #[rstest]
74    fn zero_energy_point(#[values(1.0, 2.0, 5.0, 10.0)] k: f64, #[values(0.0, 1.0, 2.0)] r_0: f64) {
75        let harmonic = Harmonic { k, r_0 };
76
77        assert_eq!(harmonic.k, k);
78        assert_eq!(harmonic.r_0, r_0);
79
80        assert_eq!(harmonic.energy(r_0), 0.0);
81        assert_eq!(harmonic.force(r_0), 0.0);
82    }
83
84    #[rstest]
85    fn general_case(#[values(1.0, 2.0, 5.0, 10.0)] k: f64, #[values(0.0, 1.0, 2.0)] r_0: f64) {
86        let r = 5.0;
87        let harmonic = Harmonic { k, r_0 };
88
89        assert_eq!(harmonic.k, k);
90        assert_eq!(harmonic.r_0, r_0);
91
92        let expected_energy = 0.5 * k * (r - r_0) * (r - r_0);
93        let expected_force = -k * (r - r_0);
94
95        assert_relative_eq!(harmonic.energy(r), expected_energy);
96        assert_relative_eq!(harmonic.force(r), expected_force);
97    }
98}