use crate::{math::Scalar, physics::molecular::potential::Potential};
#[derive(Clone, Debug)]
pub struct Morse {
pub rest_length: Scalar,
pub depth: Scalar,
pub parameter: Scalar,
}
impl Potential for Morse {
fn energy(&self, length: Scalar) -> Scalar {
self.depth * (1.0 - (self.parameter * (self.rest_length - length)).exp()).powi(2)
}
fn force(&self, length: Scalar) -> Scalar {
let exp = (self.parameter * (self.rest_length - length)).exp();
2.0 * self.parameter * self.depth * exp * (1.0 - exp)
}
fn stiffness(&self, length: Scalar) -> Scalar {
let exp = (self.parameter * (self.rest_length - length)).exp();
2.0 * self.parameter.powi(2) * self.depth * exp * (2.0 * exp - 1.0)
}
fn anharmonicity(&self, length: Scalar) -> Scalar {
let exp = (self.parameter * (self.rest_length - length)).exp();
2.0 * self.parameter.powi(3) * self.depth * exp * (1.0 - 4.0 * exp)
}
fn extension(&self, force: Scalar) -> Scalar {
let y = force / self.peak_force();
if (0.0..=1.0).contains(&y) {
(2.0 / (1.0 + (1.0 - y).sqrt())).ln() / self.parameter
} else {
Scalar::NAN
}
}
fn compliance(&self, force: Scalar) -> Scalar {
let y = force / self.peak_force();
if (0.0..1.0).contains(&y) {
let s = (1.0 - y).sqrt();
1.0 / (self.parameter.powi(2) * self.depth) / (s * (1.0 + s))
} else if y == 0.0 {
Scalar::INFINITY
} else {
Scalar::NAN
}
}
fn peak(&self) -> Scalar {
self.rest_length + 2.0_f64.ln() / self.parameter
}
fn peak_force(&self) -> Scalar {
0.5 * self.parameter * self.depth
}
fn rest_length(&self) -> Scalar {
self.rest_length
}
}