hoomd_interaction/univariate/
lennard_jones.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 [`LennardJones`]
5
6use serde::{Deserialize, Serialize};
7
8use super::{UnivariateEnergy, UnivariateForce};
9
10/// Potential with a steep repulsive core and an attractive well.
11///
12/// ```math
13/// U(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{N} - \left( \frac{\sigma}{r} \right)^{M} \right]
14/// ```
15///
16/// Compute the Lennard-Jones (LJ) potential and force as a function of `r`.
17///
18/// # Examples
19///
20/// In basic usage, the exponents `N` and `M` default to 12 and 6, respectively:
21///
22/// ```
23/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
24/// use hoomd_interaction::univariate::{
25///     LennardJones, UnivariateEnergy, UnivariateForce,
26/// };
27///
28/// let epsilon = 1.5;
29/// let sigma = 2.5;
30///
31/// let lennard_jones: LennardJones = LennardJones { epsilon, sigma };
32/// assert_abs_diff_eq!(lennard_jones.energy(sigma), 0.0);
33/// assert_relative_eq!(
34///     lennard_jones.energy(2.0_f64.powf(1.0 / 6.0) * sigma),
35///     -epsilon
36/// );
37/// assert_abs_diff_eq!(
38///     lennard_jones.force(2.0_f64.powf(1.0 / 6.0) * sigma),
39///     0.0,
40///     epsilon = 1e-12
41/// );
42/// ```
43///
44/// You can choose any values for `N` and `M` _at compile time_:
45///
46/// ```
47/// use approxim::{assert_abs_diff_eq, assert_relative_eq};
48/// use hoomd_interaction::univariate::{
49///     LennardJones, UnivariateEnergy, UnivariateForce,
50/// };
51///
52/// let epsilon = 1.5;
53/// let sigma = 2.5;
54///
55/// let lennard_jones: LennardJones<8, 4> = LennardJones { epsilon, sigma };
56/// assert_abs_diff_eq!(lennard_jones.energy(sigma), 0.0);
57/// assert_relative_eq!(
58///     lennard_jones.energy(2.0_f64.powf(1.0 / 4.0) * sigma),
59///     -epsilon
60/// );
61/// assert_abs_diff_eq!(
62///     lennard_jones.force(2.0_f64.powf(1.0 / 4.0) * sigma),
63///     0.0,
64///     epsilon = 1e-12
65/// );
66/// ```
67///
68/// The parameters are public fields and may be accessed directly:
69///
70/// ```
71/// use hoomd_interaction::univariate::LennardJones;
72///
73/// let mut lennard_jones: LennardJones = LennardJones::default();
74/// lennard_jones.epsilon = 1.5;
75/// lennard_jones.sigma = 3.0;
76/// ```
77#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
78pub struct LennardJones<const N: i32 = 12, const M: i32 = 6> {
79    /// Energy scale *(\[energy\])*.
80    pub epsilon: f64,
81    /// Interaction width *(\[length\])*.
82    pub sigma: f64,
83}
84
85impl<const N: i32, const M: i32> Default for LennardJones<N, M> {
86    /// Construct a [`LennardJones`] with default parameters (epsilon=1.0, sigma=1.0)
87    ///
88    /// # Example
89    ///
90    /// ```
91    /// use hoomd_interaction::univariate::LennardJones;
92    ///
93    /// let lennard_jones: LennardJones = LennardJones::default();
94    /// ```
95    #[inline]
96    fn default() -> Self {
97        Self {
98            epsilon: 1.0,
99            sigma: 1.0,
100        }
101    }
102}
103
104impl<const N: i32, const M: i32> UnivariateEnergy for LennardJones<N, M> {
105    #[inline]
106    fn energy(&self, r: f64) -> f64 {
107        let sigma_r = self.sigma / r;
108
109        4.0 * self.epsilon * (sigma_r.powi(N) - sigma_r.powi(M))
110    }
111}
112
113impl<const N: i32, const M: i32> UnivariateForce for LennardJones<N, M> {
114    #[inline]
115    fn force(&self, r: f64) -> f64 {
116        let r_inv = r.recip();
117        let sigma_r = self.sigma * r_inv;
118
119        self.epsilon
120            * r_inv
121            * (4.0 * f64::from(N) * sigma_r.powi(N) - 4.0 * f64::from(M) * sigma_r.powi(M))
122    }
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128    use approxim::{assert_abs_diff_eq, assert_relative_eq};
129    use rstest::*;
130
131    #[rstest]
132    fn special_points_12_6(
133        #[values(1.0, 2.0, 12.125, 0.25)] epsilon: f64,
134        #[values(1.0, 2.0, 0.5)] sigma: f64,
135    ) {
136        let lj: LennardJones = LennardJones { epsilon, sigma };
137
138        assert_eq!(lj.epsilon, epsilon);
139        assert_eq!(lj.sigma, sigma);
140
141        // Zero crossing
142        assert_abs_diff_eq!(lj.energy(sigma), 0.0);
143        assert_relative_eq!(lj.force(sigma), 24.0 * epsilon / sigma);
144
145        // Bottom of the well
146        assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 6.0) * sigma), -epsilon);
147        assert_abs_diff_eq!(
148            lj.force(2.0_f64.powf(1.0 / 6.0) * sigma),
149            0.0,
150            epsilon = 1e-12
151        );
152
153        // r = 2 sigma
154        assert_relative_eq!(lj.energy(2.0 * sigma), -63.0 / 1024.0 * epsilon);
155        assert_relative_eq!(lj.force(2.0 * sigma), -93.0 / 512.0 * epsilon / sigma);
156    }
157
158    #[rstest]
159    fn special_points_8_4(
160        #[values(1.0, 2.0, 12.125, 0.25)] epsilon: f64,
161        #[values(1.0, 2.0, 0.5)] sigma: f64,
162    ) {
163        let lj: LennardJones<8, 4> = LennardJones { epsilon, sigma };
164
165        assert_eq!(lj.epsilon, epsilon);
166        assert_eq!(lj.sigma, sigma);
167
168        // Zero crossing
169        assert_abs_diff_eq!(lj.energy(sigma), 0.0);
170        assert_relative_eq!(lj.force(sigma), 16.0 * epsilon / sigma);
171
172        // Bottom of the well
173        assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma), -epsilon);
174        assert_abs_diff_eq!(
175            lj.force(2.0_f64.powf(1.0 / 4.0) * sigma),
176            0.0,
177            epsilon = 1e-12
178        );
179
180        // r = 2 sigma
181        assert_relative_eq!(lj.energy(2.0 * sigma), -15.0 / 64.0 * epsilon);
182        assert_relative_eq!(lj.force(2.0 * sigma), -7.0 / 16.0 * epsilon / sigma);
183    }
184}
hoomd_interaction/univariate/lennard_jones.rs

hoomd_interaction/univariate/
lennard_jones.rs
