conspire/physics/molecular/potential/morse/
mod.rs

1use crate::{
2    math::{Quantity, Scalar},
3    physics::molecular::potential::Potential,
4    units::{
5        Energy, Force, ForcePerLength, Length, ReciprocalForcePerLength, ReciprocalLength, Stress,
6    },
7};
8
9/// The Morse potential.[^1]
10/// [^1]: P.M. Morse, [Physical Review **34**, 57 (1929)](https://doi.org/10.1103/PhysRev.34.57).
11#[derive(Clone, Debug)]
12pub struct Morse {
13    /// The rest length $`x_0`$.
14    pub rest_length: Scalar,
15    /// The potential depth $`u_0`$.
16    pub depth: Scalar,
17    /// The Morse parameter $`a`$.
18    pub parameter: Scalar,
19}
20
21impl Morse {
22    /// Returns the potential depth.
23    fn depth(&self) -> Quantity<Energy> {
24        Quantity::new(self.depth)
25    }
26    /// Returns the Morse parameter.
27    fn parameter(&self) -> Quantity<ReciprocalLength> {
28        Quantity::new(self.parameter)
29    }
30}
31
32impl Potential for Morse {
33    /// ```math
34    /// u(x) = u_0\left[1 - e^{-a(x - x_0)}\right]^2
35    /// ```
36    fn energy(&self, length: Quantity<Length>) -> Quantity<Energy> {
37        let exp = (self.parameter() * (self.rest_length() - length)).exp();
38        self.depth() * (1.0 - exp).powi(2)
39    }
40    /// ```math
41    /// f(x) = 2au_0e^{-a(x - x_0)}\left[1 - e^{-a(x - x_0)}\right]
42    /// ```
43    fn force(&self, length: Quantity<Length>) -> Quantity<Force> {
44        let exp = (self.parameter() * (self.rest_length() - length)).exp();
45        2.0 * self.parameter() * self.depth() * exp * (1.0 - exp)
46    }
47    /// ```math
48    /// f(u) = \pm 2a u_0\sqrt{u/u_0}\left(1 \mp \sqrt{u/u_0}\right)
49    /// ```
50    fn forces_at_energy(&self, energy: Quantity<Energy>) -> [Quantity<Force>; 2] {
51        let y = energy / self.depth();
52        let f = 2.0 * self.parameter() * self.depth() * y.sqrt();
53        let tensile = if (0.0..=1.0).contains(&y) {
54            f * (1.0 - y.sqrt())
55        } else {
56            Quantity::new(Scalar::NAN)
57        };
58        [tensile, -f * (1.0 + y.sqrt())]
59    }
60    /// ```math
61    /// k(x) = 2a^2u_0e^{-a(x - x_0)}\left[2e^{-a(x - x_0)} - 1\right]
62    /// ```
63    fn stiffness(&self, length: Quantity<Length>) -> Quantity<ForcePerLength> {
64        let exp = (self.parameter() * (self.rest_length() - length)).exp();
65        2.0 * (self.parameter() * self.parameter()) * self.depth() * exp * (2.0 * exp - 1.0)
66    }
67    /// ```math
68    /// h(x) = 2a^3u_0e^{-a(x - x_0)}\left[1 - 4e^{-a(x - x_0)}\right]
69    /// ```
70    fn anharmonicity(&self, length: Quantity<Length>) -> Quantity<Stress> {
71        let exp = (self.parameter() * (self.rest_length() - length)).exp();
72        2.0 * (self.parameter() * self.parameter())
73            * self.depth()
74            * self.parameter()
75            * exp
76            * (1.0 - 4.0 * exp)
77    }
78    /// ```math
79    /// \Delta x(f) = \frac{1}{a}\,\ln\left(\frac{2}{1 + \sqrt{1 - f/f_\mathrm{max}}}\right)
80    /// ```
81    fn extension(&self, force: Quantity<Force>) -> Quantity<Length> {
82        let y = force / self.peak_force();
83        if y <= 1.0 {
84            (2.0 / (1.0 + (1.0 - y).sqrt())).ln() / self.parameter()
85        } else {
86            Quantity::new(Scalar::NAN)
87        }
88    }
89    /// ```math
90    /// \Delta x(u) = \frac{1}{a}\,\ln\left(\frac{1}{1\mp\sqrt{u/u_0}}\right)
91    /// ```
92    fn extensions_at_energy(&self, energy: Quantity<Energy>) -> [Quantity<Length>; 2] {
93        let y = energy / self.depth();
94        let tensile = if (0.0..=1.0).contains(&y) {
95            (1.0 / (1.0 - y.sqrt())).ln() / self.parameter()
96        } else {
97            Quantity::new(Scalar::NAN)
98        };
99        [tensile, (1.0 / (1.0 + y.sqrt())).ln() / self.parameter()]
100    }
101    /// ```math
102    /// c(f) = \frac{1}{a^2u_0}\,\frac{\left(1-f/f_\mathrm{max}\right)^{-1/2}}{1+\sqrt{1-f/f_\mathrm{max}}}
103    /// ```
104    fn compliance(&self, force: Quantity<Force>) -> Quantity<ReciprocalForcePerLength> {
105        let y = force / self.peak_force();
106        if (0.0..1.0).contains(&y) {
107            let s = (1.0 - y).sqrt();
108            1.0 / (self.parameter() * self.parameter() * self.depth()) / (s * (1.0 + s))
109        } else if y == 0.0 {
110            Quantity::new(Scalar::INFINITY)
111        } else {
112            Quantity::new(Scalar::NAN)
113        }
114    }
115    /// ```math
116    /// \text{arg max }u(x) = x_0 + \frac{1}{a}\,\ln(2)
117    /// ```
118    fn peak(&self) -> Quantity<Length> {
119        self.rest_length() + 2.0_f64.ln() / self.parameter()
120    }
121    /// ```math
122    /// f(x_\mathrm{peak}) = \frac{au_0}{2}
123    /// ```
124    fn peak_force(&self) -> Quantity<Force> {
125        0.5 * self.parameter() * self.depth()
126    }
127    /// ```math
128    /// \text{arg min }u(x) = x_0
129    /// ```
130    fn rest_length(&self) -> Quantity<Length> {
131        Quantity::new(self.rest_length)
132    }
133}
conspire/physics/molecular/potential/morse/mod.rs

conspire/physics/molecular/potential/morse/
mod.rs
