interatomic/twobody/
mie.rs

1// Copyright 2023 Mikael Lund
2//
3// Licensed under the Apache license, version 2.0 (the "license");
4// you may not use this file except in compliance with the license.
5// You may obtain a copy of the license at
6//
7//     http://www.apache.org/licenses/license-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the license is distributed on an "as is" basis,
11// without warranties or conditions of any kind, either express or implied.
12// See the license for the specific language governing permissions and
13// limitations under the license.
14
15use crate::twobody::IsotropicTwobodyEnergy;
16use crate::Cutoff;
17
18#[cfg(feature = "serde")]
19use serde::{Deserialize, Serialize};
20
21/// Mie potential
22///
23/// This is a generalization of the Lennard-Jones potential due to G. Mie,
24/// ["Zur kinetischen Theorie der einatomigen Körper"](https://doi.org/10.1002/andp.19033160802).
25/// The energy is
26/// $$ u(r) = ε C \left [\left (\frac{σ}{r}\right )^n - \left (\frac{σ}{r}\right )^m \right ]$$
27/// where $C = \frac{n}{n-m} \cdot \left (\frac{n}{m}\right )^{\frac{m}{n-m}}$ and $n > m$.
28/// The Lennard-Jones potential is recovered for $n = 12$ and $m = 6$.
29///
30/// # Examples:
31/// ~~~
32/// use interatomic::twobody::*;
33/// let (epsilon, sigma, r2) = (1.5, 2.0, 2.5);
34/// let mie = Mie::<12, 6>::new(epsilon, sigma);
35/// let lj = LennardJones::new(epsilon, sigma);
36/// assert_eq!(mie.isotropic_twobody_energy(r2), lj.isotropic_twobody_energy(r2));
37/// ~~~
38
39#[derive(Clone, Debug, PartialEq, Copy)]
40#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
41pub struct Mie<const N: u32, const M: u32> {
42    /// Interaction strength, ε
43    #[cfg_attr(feature = "serde", serde(alias = "eps", alias = "ε"))]
44    epsilon: f64,
45    /// Diameter, σ
46    #[cfg_attr(feature = "serde", serde(alias = "σ"))]
47    sigma: f64,
48}
49
50impl<const N: u32, const M: u32> Mie<N, M> {
51    const C: f64 = (N / (N - M) * (N / M).pow(M / (N - M))) as f64;
52
53    /// Compile-time optimization if N and M are divisible by 2
54    const OPTIMIZE: bool = N.is_multiple_of(2) && M.is_multiple_of(2);
55    const N_OVER_M: i32 = (N / M) as i32;
56    const M_HALF: i32 = (M / 2) as i32;
57
58    /// Create a new Mie potential with epsilon (well depth) and sigma (diameter).
59    pub const fn new(epsilon: f64, sigma: f64) -> Self {
60        assert!(M > 0);
61        assert!(N > M);
62        Self { epsilon, sigma }
63    }
64}
65
66impl<const N: u32, const M: u32> IsotropicTwobodyEnergy for Mie<N, M> {
67    #[inline]
68    fn isotropic_twobody_energy(&self, distance_squared: f64) -> f64 {
69        if Self::OPTIMIZE {
70            let mth_power = (self.sigma * self.sigma / distance_squared).powi(Self::M_HALF); // (σ/r)^m
71            return Self::C * self.epsilon * (mth_power.powi(Self::N_OVER_M) - mth_power);
72        }
73        let s_over_r = self.sigma / distance_squared.sqrt(); // (σ/r)
74        Self::C * self.epsilon * (s_over_r.powi(N as i32) - s_over_r.powi(M as i32))
75    }
76
77    #[inline]
78    fn isotropic_twobody_force(&self, distance_squared: f64) -> f64 {
79        if Self::OPTIMIZE {
80            let x = (self.sigma * self.sigma / distance_squared).powi(Self::M_HALF);
81            let xn = x.powi(Self::N_OVER_M);
82            return Self::C * self.epsilon * (N as f64 * xn - M as f64 * x)
83                / (2.0 * distance_squared);
84        }
85        let s_over_r = self.sigma / distance_squared.sqrt();
86        Self::C
87            * self.epsilon
88            * (N as f64 * s_over_r.powi(N as i32) - M as f64 * s_over_r.powi(M as i32))
89            / (2.0 * distance_squared)
90    }
91}
92
93impl<const N: u32, const M: u32> Cutoff for Mie<N, M> {
94    fn cutoff(&self) -> f64 {
95        f64::INFINITY
96    }
97    fn cutoff_squared(&self) -> f64 {
98        f64::INFINITY
99    }
100    fn lower_cutoff(&self) -> f64 {
101        self.sigma * 0.6
102    }
103}
104
105#[cfg(test)]
106mod tests {
107    use super::*;
108    use crate::twobody::LennardJones;
109    use approx::assert_relative_eq;
110
111    #[test]
112    fn test_mie_force() {
113        let (eps, sigma) = (1.5, 2.0);
114        let mie = Mie::<12, 6>::new(eps, sigma);
115        let lj = LennardJones::new(eps, sigma);
116        for r in [1.5, 2.0, 2.5, 3.0, 5.0] {
117            let r2 = r * r;
118            assert_relative_eq!(
119                mie.isotropic_twobody_force(r2),
120                lj.isotropic_twobody_force(r2),
121                epsilon = 1e-10
122            );
123        }
124    }
125}
interatomic/twobody/mie.rs

interatomic/twobody/
mie.rs
