hoomd_interaction/univariate/
lennard_jones_gauss.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 [`LennardJonesGauss`]
5
6use serde::{Deserialize, Serialize};
7
8use super::{UnivariateEnergy, UnivariateForce};
9
10/// Double-well potential with a steep repulsive core
11///
12/// ```math
13/// U(r) = 1[\mathrm{energy}]\cdot\left[ \left(\frac{\ell}{r}\right)^{12} - 2\left(\frac{\ell}{r}\right)^{6}\right] - \varepsilon \exp \left(-\frac{(r/\ell-r_0)^2}{2\sigma^2}\right)
14/// ```
15///
16/// Compute the Lennard-Jones-Gauss (LJG) potential and force as a function of `r`.
17///
18/// # Examples
19///
20/// ```
21/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
22/// use hoomd_interaction::univariate::{
23///     LennardJonesGauss, UnivariateEnergy, UnivariateForce,
24/// };
25///
26/// let epsilon = 0.5;
27/// let sigma_squared = 0.5;
28/// let r_0 = 0.5_f64.powf(1.0 / 6.0);
29/// let scale = 1.0_f64;
30///
31/// let lennard_jones_gauss: LennardJonesGauss = LennardJonesGauss {
32///     epsilon,
33///     sigma_squared,
34///     r_0,
35///     scale,
36/// };
37/// assert_relative_eq!(
38///     lennard_jones_gauss.energy(0.5_f64.powf(1.0 / 6.0)),
39///     -epsilon,
40///     epsilon = 1e-12
41/// );
42/// ```
43///
44/// The parameters are public fields and may be accessed directly:
45///
46/// ```
47/// use hoomd_interaction::univariate::LennardJonesGauss;
48///
49/// let mut lennard_jones_gauss: LennardJonesGauss = LennardJonesGauss {
50///     epsilon: 1.5,
51///     sigma_squared: 0.02,
52///     r_0: 3.2,
53///     scale: 1.0,
54/// };
55/// lennard_jones_gauss.epsilon = 1.5;
56/// lennard_jones_gauss.sigma_squared = 0.02;
57/// lennard_jones_gauss.r_0 = 3.2;
58/// ```
59#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
60pub struct LennardJonesGauss {
61    /// Scale of Gaussian, in units of energy
62    pub epsilon: f64,
63    /// Width of Gaussian, sigma^2, unitless
64    pub sigma_squared: f64,
65    /// Gaussian center, unitless
66    pub r_0: f64,
67    /// unit of the length scale
68    pub scale: f64,
69}
70
71impl UnivariateEnergy for LennardJonesGauss {
72    #[inline]
73    fn energy(&self, r: f64) -> f64 {
74        let r_inv = self.scale / r;
75        let arg = -((r / self.scale) - self.r_0).powi(2) / (2.0 * self.sigma_squared);
76        r_inv.powi(12) - 2.0 * r_inv.powi(6) - self.epsilon * arg.exp()
77    }
78}
79
80impl UnivariateForce for LennardJonesGauss {
81    #[inline]
82    fn force(&self, r: f64) -> f64 {
83        let r_inv = self.scale / r;
84        let arg = -((r / self.scale) - self.r_0).powi(2) / (2.0 * self.sigma_squared);
85        12.0 * (r_inv.powi(13) - r_inv.powi(7))
86            - (self.epsilon * ((r / self.scale) - self.r_0) / self.sigma_squared) * arg.exp()
87    }
88}
89
90#[cfg(test)]
91mod tests {
92    use super::*;
93    use approxim::{assert_abs_diff_eq, assert_relative_eq};
94    use rstest::*;
95
96    #[rstest]
97    fn select_points_first_set() {
98        let lj_gauss: LennardJonesGauss = LennardJonesGauss {
99            epsilon: 2.0,
100            sigma_squared: 0.5,
101            r_0: 3.0,
102            scale: 1.0,
103        };
104
105        assert_eq!(lj_gauss.epsilon, 2.0);
106        assert_eq!(lj_gauss.sigma_squared, 0.5);
107        assert_eq!(lj_gauss.r_0, 3.0);
108
109        // numeric tests
110        assert_relative_eq!(
111            lj_gauss.energy(1.5_f64),
112            -0.378_674_092_892_274,
113            epsilon = 1e-12
114        );
115        assert_abs_diff_eq!(
116            lj_gauss.force(1.5_f64),
117            -0.008_277_841_185_963,
118            epsilon = 1e-12
119        );
120        assert_relative_eq!(
121            lj_gauss.energy(3.2_f64),
122            -1.923_440_656_092_139,
123            epsilon = 1e-12
124        );
125        assert_abs_diff_eq!(
126            lj_gauss.force(3.2_f64),
127            -0.772_120_758_370_149,
128            epsilon = 1e-12
129        );
130    }
131    #[rstest]
132    fn select_points_second_set() {
133        let lj_gauss: LennardJonesGauss = LennardJonesGauss {
134            epsilon: 10.0,
135            sigma_squared: 0.1,
136            r_0: 5.0,
137            scale: 1.0,
138        };
139
140        assert_eq!(lj_gauss.epsilon, 10.0);
141        assert_eq!(lj_gauss.sigma_squared, 0.1);
142        assert_eq!(lj_gauss.r_0, 5.0);
143
144        // numeric tests
145        assert_relative_eq!(
146            lj_gauss.energy(1.5_f64),
147            -0.167_875_643_768_546,
148            epsilon = 1e-12
149        );
150        assert_abs_diff_eq!(
151            lj_gauss.force(1.5_f64),
152            -0.640_673_188_557_149,
153            epsilon = 1e-12
154        );
155        assert_relative_eq!(
156            lj_gauss.energy(3.2_f64),
157            -0.001_862_699_147_576_425,
158            epsilon = 1e-12
159        );
160        assert_abs_diff_eq!(
161            lj_gauss.force(3.2_f64),
162            -0.003_472_622_566_788_369,
163            epsilon = 1e-12
164        );
165    }
166}
hoomd_interaction/univariate/lennard_jones_gauss.rs

hoomd_interaction/univariate/
lennard_jones_gauss.rs
