hoomd_interaction/univariate/
shifted.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 [`Shifted`]
5
6use serde::{Deserialize, Serialize};
7
8use super::{UnivariateEnergy, UnivariateForce};
9
10/// Shift another potential to 0 at a given `r`.
11///
12/// ```math
13/// U(r) = f(r) - f(r_\mathrm{shift})
14/// ```
15///
16/// # Example
17///
18/// Shifted Lennard-Jones:
19/// ```
20/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
21/// use hoomd_interaction::univariate::{
22///     LennardJones, Shifted, UnivariateEnergy,
23/// };
24///
25/// let epsilon = 1.5;
26/// let sigma = 1.0;
27/// let r_shift = 2.5;
28/// let lj: LennardJones = LennardJones { epsilon, sigma };
29/// let shifted_lj = Shifted {
30///     f: lj.clone(),
31///     r_shift,
32/// };
33///
34/// assert_abs_diff_eq!(shifted_lj.energy(r_shift), 0.0);
35/// assert_relative_eq!(
36///     shifted_lj.energy(2.0_f64.powf(1.0 / 6.0) * sigma),
37///     -epsilon - lj.energy(r_shift)
38/// );
39/// ```
40///
41/// Fields can be accessed directly, including those of the original potential `f`:
42/// ```
43/// use hoomd_interaction::univariate::{LennardJones, Shifted};
44///
45/// let epsilon = 1.5;
46/// let sigma = 1.0;
47/// let r_shift = 2.5;
48/// let mut shifted_lj = Shifted {
49///     f: LennardJones::<12, 6> { epsilon, sigma },
50///     r_shift,
51/// };
52///
53/// shifted_lj.r_shift = 3.0;
54/// shifted_lj.f.sigma = 1.2;
55/// ```
56#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
57pub struct Shifted<F> {
58    /// The original potential.
59    pub f: F,
60    /// `r` value *(\[length\])* where the shifted potential will be 0.
61    pub r_shift: f64,
62}
63
64impl<F> Default for Shifted<F>
65where
66    F: Default,
67{
68    /// Construct a shifted potential with default parameters
69    ///
70    /// The defaults are:
71    /// `f = F::default()`
72    /// `r_shift = 0.0`
73    ///
74    /// # Example
75    ///
76    /// ```
77    /// use hoomd_interaction::univariate::{LennardJones, Shifted};
78    ///
79    /// let shifted_lj = Shifted::<LennardJones>::default();
80    /// ```
81    #[inline]
82    fn default() -> Self {
83        Self {
84            f: F::default(),
85            r_shift: 0.0,
86        }
87    }
88}
89
90impl<F: UnivariateEnergy> UnivariateEnergy for Shifted<F> {
91    #[inline]
92    fn energy(&self, r: f64) -> f64 {
93        self.f.energy(r) - self.f.energy(self.r_shift)
94    }
95}
96
97impl<F: UnivariateForce> UnivariateForce for Shifted<F> {
98    #[inline]
99    fn force(&self, r: f64) -> f64 {
100        self.f.force(r)
101    }
102}
103
104#[cfg(test)]
105mod tests {
106    use super::*;
107    use approxim::{assert_abs_diff_eq, assert_relative_eq};
108    use rstest::*;
109
110    use crate::univariate::LennardJones;
111
112    #[rstest]
113    fn special_points_12_6(
114        #[values(1.0, 2.0, 12.125, 0.25)] epsilon: f64,
115        #[values(1.0, 2.0, 0.5)] sigma: f64,
116    ) {
117        let lj: LennardJones = LennardJones { epsilon, sigma };
118        let r_shift = 2.5 * sigma;
119        let shifted_lj = Shifted {
120            f: lj.clone(),
121            r_shift,
122        };
123
124        assert_eq!(shifted_lj.f.epsilon, epsilon);
125        assert_eq!(shifted_lj.f.sigma, sigma);
126        assert_eq!(shifted_lj.r_shift, r_shift);
127
128        let delta_e = lj.energy(r_shift);
129
130        // Zero crossing
131        assert_abs_diff_eq!(shifted_lj.energy(sigma), 0.0 - delta_e);
132        assert_relative_eq!(shifted_lj.force(sigma), 24.0 * epsilon / sigma);
133
134        // Bottom of the well
135        assert_relative_eq!(
136            shifted_lj.energy(2.0_f64.powf(1.0 / 6.0) * sigma),
137            -epsilon - delta_e
138        );
139        assert_abs_diff_eq!(
140            shifted_lj.force(2.0_f64.powf(1.0 / 6.0) * sigma),
141            0.0,
142            epsilon = 1e-12
143        );
144
145        // r = 2 sigma
146        assert_relative_eq!(
147            shifted_lj.energy(2.0 * sigma),
148            -63.0 / 1024.0 * epsilon - delta_e
149        );
150        assert_relative_eq!(
151            shifted_lj.force(2.0 * sigma),
152            -93.0 / 512.0 * epsilon / sigma
153        );
154    }
155}