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

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

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