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///     epsilon = 1e-12
61/// );
62/// assert_abs_diff_eq!(
63///     lennard_jones.force(2.0_f64.powf(1.0 / 4.0) * sigma),
64///     0.0,
65///     epsilon = 1e-12
66/// );
67/// ```
68///
69/// The parameters are public fields and may be accessed directly:
70///
71/// ```
72/// use hoomd_interaction::univariate::LennardJones;
73///
74/// let mut lennard_jones: LennardJones = LennardJones::default();
75/// lennard_jones.epsilon = 1.5;
76/// lennard_jones.sigma = 3.0;
77/// ```
78#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
79pub struct LennardJones<const N: i32 = 12, const M: i32 = 6> {
80    /// Energy scale *(\[energy\])*.
81    pub epsilon: f64,
82    /// Interaction width *(\[length\])*.
83    pub sigma: f64,
84}
85
86impl<const N: i32, const M: i32> Default for LennardJones<N, M> {
87    /// Construct a [`LennardJones`] with default parameters (epsilon=1.0, sigma=1.0)
88    ///
89    /// # Example
90    ///
91    /// ```
92    /// use hoomd_interaction::univariate::LennardJones;
93    ///
94    /// let lennard_jones: LennardJones = LennardJones::default();
95    /// ```
96    #[inline]
97    fn default() -> Self {
98        Self {
99            epsilon: 1.0,
100            sigma: 1.0,
101        }
102    }
103}
104
105impl<const N: i32, const M: i32> UnivariateEnergy for LennardJones<N, M> {
106    #[inline]
107    fn energy(&self, r: f64) -> f64 {
108        let sigma_r = self.sigma / r;
109
110        4.0 * self.epsilon * (sigma_r.powi(N) - sigma_r.powi(M))
111    }
112}
113
114impl<const N: i32, const M: i32> UnivariateForce for LennardJones<N, M> {
115    #[inline]
116    fn force(&self, r: f64) -> f64 {
117        let r_inv = r.recip();
118        let sigma_r = self.sigma * r_inv;
119
120        self.epsilon
121            * r_inv
122            * (4.0 * f64::from(N) * sigma_r.powi(N) - 4.0 * f64::from(M) * sigma_r.powi(M))
123    }
124}
125
126#[cfg(test)]
127mod tests {
128    use super::*;
129    use approxim::{assert_abs_diff_eq, assert_relative_eq};
130    use rstest::*;
131
132    #[rstest]
133    fn special_points_12_6(
134        #[values(1.0, 2.0, 12.125, 0.25)] epsilon: f64,
135        #[values(1.0, 2.0, 0.5)] sigma: f64,
136    ) {
137        let lj: LennardJones = LennardJones { epsilon, sigma };
138
139        assert_eq!(lj.epsilon, epsilon);
140        assert_eq!(lj.sigma, sigma);
141
142        // Zero crossing
143        assert_abs_diff_eq!(lj.energy(sigma), 0.0);
144        assert_relative_eq!(lj.force(sigma), 24.0 * epsilon / sigma);
145
146        // Bottom of the well
147        assert_relative_eq!(lj.energy(2.0_f64.powf(1.0 / 6.0) * sigma), -epsilon);
148        assert_abs_diff_eq!(
149            lj.force(2.0_f64.powf(1.0 / 6.0) * sigma),
150            0.0,
151            epsilon = 1e-12
152        );
153
154        // r = 2 sigma
155        assert_relative_eq!(lj.energy(2.0 * sigma), -63.0 / 1024.0 * epsilon);
156        assert_relative_eq!(lj.force(2.0 * sigma), -93.0 / 512.0 * epsilon / sigma);
157    }
158
159    #[rstest]
160    fn special_points_8_4(
161        #[values(1.0, 2.0, 12.125, 0.25)] epsilon: f64,
162        #[values(1.0, 2.0, 0.5)] sigma: f64,
163    ) {
164        let lj: LennardJones<8, 4> = LennardJones { epsilon, sigma };
165
166        assert_eq!(lj.epsilon, epsilon);
167        assert_eq!(lj.sigma, sigma);
168
169        // Zero crossing
170        assert_abs_diff_eq!(lj.energy(sigma), 0.0);
171        assert_relative_eq!(lj.force(sigma), 16.0 * epsilon / sigma);
172
173        // Bottom of the well
174        assert_relative_eq!(
175            lj.energy(2.0_f64.powf(1.0 / 4.0) * sigma),
176            -epsilon,
177            epsilon = 1e-12
178        );
179        assert_abs_diff_eq!(
180            lj.force(2.0_f64.powf(1.0 / 4.0) * sigma),
181            0.0,
182            epsilon = 1e-12
183        );
184
185        // r = 2 sigma
186        assert_relative_eq!(lj.energy(2.0 * sigma), -15.0 / 64.0 * epsilon);
187        assert_relative_eq!(lj.force(2.0 * sigma), -7.0 / 16.0 * epsilon / sigma);
188    }
189}
hoomd_interaction/univariate/lennard_jones.rs

hoomd_interaction/univariate/
lennard_jones.rs
