potentials/pair/
mie.rs

1//! # Generalized Mie Potential (n-m)
2//!
3//! The Mie potential generalizes Lennard-Jones with arbitrary integer exponents.
4//!
5//! ## Formula
6//!
7//! ```text
8//! V(r) = C * epsilon * [(sigma/r)^n - (sigma/r)^m]
9//!
10//! where:
11//!     C = (n / (n - m)) * (n / m)^(m / (n - m))
12//! ```
13//!
14//! The prefactor C ensures that:
15//! - The minimum energy is exactly -epsilon
16//! - r_min = (n/m)^(1/(n-m)) * sigma
17//!
18//! ## Special Cases
19//!
20//! | N  | M | Name |
21//! |----|---|------|
22//! | 12 | 6 | Lennard-Jones |
23//! | 9  | 6 | Soft-core LJ |
24//! | 14 | 7 | Harder repulsion |
25//!
26//! ## Implementation Notes
27//!
28//! - Uses const generics for compile-time optimization of exponents
29//! - Exponent computations are fully unrolled by the compiler
30//! - Zero runtime overhead for exponent handling
31
32use crate::base::Potential2;
33use crate::math::Vector;
34
35/// Generalized Mie (n-m) potential.
36///
37/// ## Type Parameters
38///
39/// - `T`: Numeric type (f32, f64, or SIMD vector)
40/// - `N`: Repulsive exponent (must be > M)
41/// - `M`: Attractive exponent (must be > 0)
42///
43/// ## Parameters
44///
45/// - `epsilon`: Well depth (energy)
46/// - `sigma`: Characteristic length (length)
47///
48/// ## Precomputed Values
49///
50/// Stores `cn = C * epsilon * sigma^n` and `cm = C * epsilon * sigma^m`
51/// for efficient evaluation.
52#[derive(Clone, Copy, Debug, PartialEq)]
53pub struct Mie<T, const N: u32, const M: u32> {
54    /// Coefficient for r^-n term: C * epsilon * sigma^n
55    cn: T,
56    /// Coefficient for r^-m term: C * epsilon * sigma^m
57    cm: T,
58}
59
60impl<T: Vector, const N: u32, const M: u32> Mie<T, N, M> {
61    /// Creates a new Mie potential.
62    ///
63    /// ## Arguments
64    ///
65    /// - `epsilon`: Well depth (energy units)
66    /// - `sigma`: Characteristic length (length units)
67    ///
68    /// ## Panics
69    ///
70    /// Panics if `N <= M` or `M == 0`.
71    ///
72    /// ## Example
73    ///
74    /// ```
75    /// use potentials::pair::Mie;
76    ///
77    /// // Standard LJ 12-6
78    /// let lj = Mie::<f64, 12, 6>::new(1.0, 1.0);
79    ///
80    /// // Softer 9-6 potential
81    /// let soft = Mie::<f64, 9, 6>::new(1.0, 1.0);
82    /// ```
83    #[inline]
84    pub fn new(epsilon: f64, sigma: f64) -> Self {
85        assert!(N > M, "Mie potential requires N > M, got N={}, M={}", N, M);
86        assert!(M > 0, "Mie potential requires M > 0, got M={}", M);
87
88        let n = N as f64;
89        let m = M as f64;
90
91        // Prefactor: C = (n/(n-m)) * (n/m)^(m/(n-m))
92        let n_minus_m = n - m;
93        let prefactor = (n / n_minus_m) * (n / m).powf(m / n_minus_m);
94        let c_eps = prefactor * epsilon;
95
96        // Precompute sigma^n and sigma^m using integer powers
97        let sigma_n = int_pow(sigma, N);
98        let sigma_m = int_pow(sigma, M);
99
100        Self {
101            cn: T::splat(c_eps * sigma_n),
102            cm: T::splat(c_eps * sigma_m),
103        }
104    }
105
106    /// Creates a new Mie potential from precomputed coefficients.
107    ///
108    /// ## Arguments
109    ///
110    /// - `cn`: Coefficient for r^-n term (C * epsilon * sigma^n)
111    /// - `cm`: Coefficient for r^-m term (C * epsilon * sigma^m)
112    ///
113    /// This is useful when mixing rules produce cn/cm directly.
114    #[inline]
115    pub fn from_coefficients(cn: f64, cm: f64) -> Self {
116        Self {
117            cn: T::splat(cn),
118            cm: T::splat(cm),
119        }
120    }
121
122    /// Returns the cn coefficient.
123    #[inline]
124    pub fn cn(&self) -> T {
125        self.cn
126    }
127
128    /// Returns the cm coefficient.
129    #[inline]
130    pub fn cm(&self) -> T {
131        self.cm
132    }
133
134    /// Returns the repulsive exponent N.
135    #[inline]
136    pub const fn n(&self) -> u32 {
137        N
138    }
139
140    /// Returns the attractive exponent M.
141    #[inline]
142    pub const fn m(&self) -> u32 {
143        M
144    }
145}
146
147impl<T: Vector, const N: u32, const M: u32> Potential2<T> for Mie<T, N, M> {
148    /// Computes the potential energy.
149    ///
150    /// ```text
151    /// V(r) = cn/r^n - cm/r^m
152    /// ```
153    #[inline(always)]
154    fn energy(&self, r_sq: T) -> T {
155        let r_inv = r_sq.rsqrt();
156
157        // Compute r^-n and r^-m using optimized integer powers
158        let r_neg_n = int_pow_vec::<T, N>(r_inv);
159        let r_neg_m = int_pow_vec::<T, M>(r_inv);
160
161        self.cn * r_neg_n - self.cm * r_neg_m
162    }
163
164    /// Computes the force factor.
165    ///
166    /// ```text
167    /// dV/dr = -n*cn/r^(n+1) + m*cm/r^(m+1)
168    /// S = -(dV/dr)/r = n*cn/r^(n+2) - m*cm/r^(m+2)
169    /// ```
170    #[inline(always)]
171    fn force_factor(&self, r_sq: T) -> T {
172        let r_inv = r_sq.rsqrt();
173        let r_sq_inv = r_sq.recip();
174
175        let r_neg_n = int_pow_vec::<T, N>(r_inv);
176        let r_neg_m = int_pow_vec::<T, M>(r_inv);
177
178        let n = T::splat(N as f64);
179        let m = T::splat(M as f64);
180
181        (n * self.cn * r_neg_n - m * self.cm * r_neg_m) * r_sq_inv
182    }
183
184    /// Computes energy and force factor together (optimized).
185    ///
186    /// Shares the computation of r^-n and r^-m.
187    #[inline(always)]
188    fn energy_force(&self, r_sq: T) -> (T, T) {
189        let r_inv = r_sq.rsqrt();
190        let r_sq_inv = r_sq.recip();
191
192        let r_neg_n = int_pow_vec::<T, N>(r_inv);
193        let r_neg_m = int_pow_vec::<T, M>(r_inv);
194
195        let cn_rn = self.cn * r_neg_n;
196        let cm_rm = self.cm * r_neg_m;
197
198        let energy = cn_rn - cm_rm;
199
200        let n = T::splat(N as f64);
201        let m = T::splat(M as f64);
202        let force = (n * cn_rn - m * cm_rm) * r_sq_inv;
203
204        (energy, force)
205    }
206}
207
208/// Computes x^n using binary exponentiation for f64.
209///
210/// Used for precomputation of sigma^n in constructor.
211#[inline]
212fn int_pow(x: f64, n: u32) -> f64 {
213    match n {
214        0 => 1.0,
215        1 => x,
216        2 => x * x,
217        3 => x * x * x,
218        4 => {
219            let x2 = x * x;
220            x2 * x2
221        }
222        5 => {
223            let x2 = x * x;
224            x2 * x2 * x
225        }
226        6 => {
227            let x2 = x * x;
228            x2 * x2 * x2
229        }
230        7 => {
231            let x2 = x * x;
232            let x4 = x2 * x2;
233            x4 * x2 * x
234        }
235        8 => {
236            let x2 = x * x;
237            let x4 = x2 * x2;
238            x4 * x4
239        }
240        9 => {
241            let x2 = x * x;
242            let x4 = x2 * x2;
243            let x8 = x4 * x4;
244            x8 * x
245        }
246        10 => {
247            let x2 = x * x;
248            let x4 = x2 * x2;
249            let x8 = x4 * x4;
250            x8 * x2
251        }
252        11 => {
253            let x2 = x * x;
254            let x4 = x2 * x2;
255            let x8 = x4 * x4;
256            x8 * x2 * x
257        }
258        12 => {
259            let x2 = x * x;
260            let x4 = x2 * x2;
261            x4 * x4 * x4
262        }
263        13 => {
264            let x2 = x * x;
265            let x4 = x2 * x2;
266            x4 * x4 * x4 * x
267        }
268        14 => {
269            let x2 = x * x;
270            let x4 = x2 * x2;
271            let x8 = x4 * x4;
272            x8 * x4 * x2
273        }
274        _ => {
275            // Binary exponentiation for large n
276            let mut result = 1.0;
277            let mut base = x;
278            let mut exp = n;
279            while exp > 0 {
280                if exp & 1 == 1 {
281                    result *= base;
282                }
283                base *= base;
284                exp >>= 1;
285            }
286            result
287        }
288    }
289}
290
291/// Computes x^N for Vector types at compile time using const generics.
292///
293/// The compiler completely eliminates the match at compile time,
294/// generating optimal code for each specific power.
295#[inline(always)]
296fn int_pow_vec<T: Vector, const N: u32>(x: T) -> T {
297    match N {
298        0 => T::splat(1.0),
299        1 => x,
300        2 => x * x,
301        3 => x * x * x,
302        4 => {
303            let x2 = x * x;
304            x2 * x2
305        }
306        5 => {
307            let x2 = x * x;
308            x2 * x2 * x
309        }
310        6 => {
311            let x2 = x * x;
312            x2 * x2 * x2
313        }
314        7 => {
315            let x2 = x * x;
316            let x4 = x2 * x2;
317            x4 * x2 * x
318        }
319        8 => {
320            let x2 = x * x;
321            let x4 = x2 * x2;
322            x4 * x4
323        }
324        9 => {
325            let x2 = x * x;
326            let x4 = x2 * x2;
327            let x8 = x4 * x4;
328            x8 * x
329        }
330        10 => {
331            let x2 = x * x;
332            let x4 = x2 * x2;
333            let x8 = x4 * x4;
334            x8 * x2
335        }
336        11 => {
337            let x2 = x * x;
338            let x4 = x2 * x2;
339            let x8 = x4 * x4;
340            x8 * x2 * x
341        }
342        12 => {
343            let x2 = x * x;
344            let x4 = x2 * x2;
345            x4 * x4 * x4
346        }
347        13 => {
348            let x2 = x * x;
349            let x4 = x2 * x2;
350            x4 * x4 * x4 * x
351        }
352        14 => {
353            let x2 = x * x;
354            let x4 = x2 * x2;
355            let x8 = x4 * x4;
356            x8 * x4 * x2
357        }
358        _ => {
359            // Binary exponentiation for large N
360            let mut result = T::splat(1.0);
361            let mut base = x;
362            let mut exp = N;
363            while exp > 0 {
364                if exp & 1 == 1 {
365                    result = result * base;
366                }
367                base = base * base;
368                exp >>= 1;
369            }
370            result
371        }
372    }
373}
374
375#[cfg(test)]
376mod tests {
377    use super::*;
378    use approx::assert_relative_eq;
379
380    #[test]
381    fn test_mie_reduces_to_lj() {
382        let mie: Mie<f64, 12, 6> = Mie::new(1.0, 1.0);
383        let lj = crate::pair::Lj::<f64>::new(1.0, 1.0);
384
385        let r_sq = 1.5 * 1.5;
386        let e_mie = mie.energy(r_sq);
387        let e_lj = lj.energy(r_sq);
388
389        assert_relative_eq!(e_mie, e_lj, epsilon = 1e-10);
390    }
391
392    #[test]
393    fn test_mie_at_minimum() {
394        let epsilon = 1.5;
395        let sigma = 2.0;
396        let mie: Mie<f64, 10, 5> = Mie::new(epsilon, sigma);
397
398        let r_min = 2.0_f64.powf(1.0 / 5.0) * sigma;
399        let r_sq = r_min * r_min;
400        let energy = mie.energy(r_sq);
401
402        assert_relative_eq!(energy, -epsilon, epsilon = 1e-10);
403    }
404
405    #[test]
406    fn test_mie_force_at_minimum() {
407        let mie: Mie<f64, 12, 6> = Mie::new(1.0, 1.0);
408
409        let r_min = 2.0_f64.powf(1.0 / 6.0);
410        let r_sq = r_min * r_min;
411        let force = mie.force_factor(r_sq);
412
413        assert_relative_eq!(force, 0.0, epsilon = 1e-10);
414    }
415
416    #[test]
417    fn test_mie_energy_force_consistency() {
418        let mie: Mie<f64, 9, 6> = Mie::new(0.5, 2.0);
419        let r_sq = 6.25;
420
421        let (e1, f1) = mie.energy_force(r_sq);
422        let e2 = mie.energy(r_sq);
423        let f2 = mie.force_factor(r_sq);
424
425        assert_relative_eq!(e1, e2, epsilon = 1e-12);
426        assert_relative_eq!(f1, f2, epsilon = 1e-12);
427    }
428
429    #[test]
430    fn test_mie_numerical_derivative() {
431        let mie: Mie<f64, 9, 6> = Mie::new(1.0, 1.0);
432        let r = 1.2;
433        let r_sq = r * r;
434
435        let h = 1e-6;
436        let v_plus = mie.energy((r + h) * (r + h));
437        let v_minus = mie.energy((r - h) * (r - h));
438        let dv_dr_numerical = (v_plus - v_minus) / (2.0 * h);
439
440        let s_numerical = -dv_dr_numerical / r;
441        let s_analytical = mie.force_factor(r_sq);
442
443        assert_relative_eq!(s_analytical, s_numerical, epsilon = 1e-6);
444    }
445}
potentials/pair/mie.rs

potentials/pair/
mie.rs
