potentials/bond/
g96.rs

1//! # GROMOS96 Bond Potential
2//!
3//! The GROMOS96 quartic bond potential, which avoids expensive
4//! square root operations by working directly with r^2.
5//!
6//! ## Formula
7//!
8//! ```text
9//! V(r) = (k/4) * (r^2 - r0^2)^2
10//! ```
11//!
12//! where:
13//! - `k`: Force constant (energy/length^4 units)
14//! - `r0`: Equilibrium bond length (length units)
15//!
16//! ## Force Factor
17//!
18//! ```text
19//! S = -(dV/dr) / r = -k * (r^2 - r0^2)
20//! ```
21//!
22//! ## Implementation Notes
23//!
24//! - Operates entirely on r^2, avoiding sqrt and division by r
25//! - Near equilibrium: V_G96 ≈ k * r0^2 * (r - r0)^2
26//! - Relation to harmonic: k_G96 = k_harm / r0^2
27
28use crate::base::Potential2;
29use crate::math::Vector;
30
31/// GROMOS96 bond potential.
32///
33/// ## Parameters
34///
35/// - `k`: Force constant (energy/length^4 units)
36/// - `r0`: Equilibrium distance (length units)
37///
38/// ## Precomputed Values
39///
40/// - `k_quarter`: Stores `k/4` for efficient computation
41/// - `r0_sq`: Stores `r0^2` to avoid sqrt
42#[derive(Clone, Copy, Debug, PartialEq)]
43pub struct G96<T> {
44    /// Force constant divided by 4
45    k_quarter: T,
46    /// Equilibrium distance squared
47    r0_sq: T,
48}
49
50impl<T: Vector> G96<T> {
51    /// Creates a new GROMOS96 bond potential.
52    ///
53    /// ## Arguments
54    ///
55    /// - `k`: Force constant (energy/length^4 units)
56    /// - `r0`: Equilibrium bond length (length units)
57    ///
58    /// ## Example
59    ///
60    /// ```
61    /// use potentials::bond::G96;
62    ///
63    /// // Convert from harmonic: k_g96 = k_harm / r0^2
64    /// let k_harm = 300.0;
65    /// let r0 = 1.54;
66    /// let bond = G96::<f64>::new(k_harm / (r0 * r0), r0);
67    /// ```
68    #[inline]
69    pub fn new(k: f64, r0: f64) -> Self {
70        Self {
71            k_quarter: T::splat(k / 4.0),
72            r0_sq: T::splat(r0 * r0),
73        }
74    }
75
76    /// Creates from r0 squared directly (avoids a multiplication).
77    #[inline]
78    pub fn from_r0_sq(k: f64, r0_sq: f64) -> Self {
79        Self {
80            k_quarter: T::splat(k / 4.0),
81            r0_sq: T::splat(r0_sq),
82        }
83    }
84
85    /// Returns the equilibrium distance squared.
86    #[inline]
87    pub fn r0_sq(&self) -> T {
88        self.r0_sq
89    }
90}
91
92impl<T: Vector> Potential2<T> for G96<T> {
93    /// Computes the potential energy.
94    ///
95    /// ```text
96    /// V(r) = (k/4) * (r^2 - r0^2)^2
97    /// ```
98    ///
99    /// Note: No sqrt required!
100    #[inline(always)]
101    fn energy(&self, r_sq: T) -> T {
102        let delta = r_sq - self.r0_sq;
103        self.k_quarter * delta * delta
104    }
105
106    /// Computes the force factor.
107    ///
108    /// ```text
109    /// dV/dr = (k/4) * 2 * (r^2 - r0^2) * 2r = k * r * (r^2 - r0^2)
110    /// S = -(dV/dr)/r = -k * (r^2 - r0^2)
111    /// ```
112    ///
113    /// Note: No sqrt or division required!
114    #[inline(always)]
115    fn force_factor(&self, r_sq: T) -> T {
116        let delta = r_sq - self.r0_sq;
117        let four = T::splat(4.0);
118        T::zero() - four * self.k_quarter * delta
119    }
120
121    /// Computes energy and force factor together (optimized).
122    ///
123    /// Shares the computation of `delta`.
124    #[inline(always)]
125    fn energy_force(&self, r_sq: T) -> (T, T) {
126        let delta = r_sq - self.r0_sq;
127        let four = T::splat(4.0);
128
129        let energy = self.k_quarter * delta * delta;
130        let force = T::zero() - four * self.k_quarter * delta;
131
132        (energy, force)
133    }
134}
135
136#[cfg(test)]
137mod tests {
138    use super::*;
139    use approx::assert_relative_eq;
140
141    #[test]
142    fn test_g96_at_equilibrium() {
143        let g96: G96<f64> = G96::new(100.0, 1.5);
144
145        let r0 = 1.5;
146        let energy = g96.energy(r0 * r0);
147        assert_relative_eq!(energy, 0.0, epsilon = 1e-10);
148    }
149
150    #[test]
151    fn test_g96_force_at_equilibrium() {
152        let g96: G96<f64> = G96::new(100.0, 1.5);
153        let r0 = 1.5;
154
155        let force = g96.force_factor(r0 * r0);
156        assert_relative_eq!(force, 0.0, epsilon = 1e-10);
157    }
158
159    #[test]
160    fn test_g96_displaced() {
161        let k = 4.0;
162        let r0 = 1.0;
163        let g96: G96<f64> = G96::new(k, r0);
164
165        let r_sq = 2.0;
166        let energy = g96.energy(r_sq);
167
168        assert_relative_eq!(energy, 1.0, epsilon = 1e-10);
169    }
170
171    #[test]
172    fn test_g96_numerical_derivative() {
173        let g96: G96<f64> = G96::new(100.0, 1.5);
174        let r = 1.6;
175        let r_sq = r * r;
176
177        let h = 1e-6;
178        let v_plus = g96.energy((r + h) * (r + h));
179        let v_minus = g96.energy((r - h) * (r - h));
180        let dv_dr_numerical = (v_plus - v_minus) / (2.0 * h);
181
182        let s_numerical = -dv_dr_numerical / r;
183        let s_analytical = g96.force_factor(r_sq);
184
185        assert_relative_eq!(s_analytical, s_numerical, epsilon = 1e-6);
186    }
187
188    #[test]
189    fn test_g96_from_r0_sq() {
190        let g96_a: G96<f64> = G96::new(100.0, 1.5);
191        let g96_b: G96<f64> = G96::from_r0_sq(100.0, 1.5 * 1.5);
192
193        let r_sq = 2.0;
194        assert_relative_eq!(g96_a.energy(r_sq), g96_b.energy(r_sq), epsilon = 1e-10);
195    }
196}
potentials/bond/g96.rs

potentials/bond/
g96.rs
