dreid_kernel/potentials/nonbonded/
vdw.rs

1use crate::math::Real;
2use crate::traits::PairKernel;
3use crate::types::EnergyDiff;
4
5/// Lennard-Jones 12-6 potential for Van der Waals interactions.
6///
7/// # Physics
8///
9/// Models the short-range repulsion and long-range attraction between non-bonded atoms.
10///
11/// - **Formula**: $$ E = D_0 \left[ \left( \frac{R_0}{r} \right)^{12} - 2 \left( \frac{R_0}{r} \right)^6 \right] $$
12/// - **Derivative Factor (`diff`)**: $$ D = -\frac{1}{r} \frac{dE}{dr} = \frac{12 D_0}{r^2} \left[ \left( \frac{R_0}{r} \right)^{12} - \left( \frac{R_0}{r} \right)^6 \right] $$
13///
14/// # Parameters
15///
16/// - `d0`: The energy well depth $D_0$.
17/// - `r0_sq`: The squared equilibrium distance $R_0^2$.
18///
19/// # Inputs
20///
21/// - `r_sq`: Squared distance $r^2$ between two atoms.
22///
23/// # Implementation Notes
24///
25/// - This implementation avoids `sqrt` entirely by operating on squared distances.
26/// - The power chain `s -> s3 -> s6` is used for efficient calculation of $r^{-6}$ and $r^{-12}$ terms.
27/// - All intermediate calculations are shared between energy and force computations.
28/// - Branchless and panic-free.
29#[derive(Clone, Copy, Debug, Default)]
30pub struct LennardJones;
31
32impl<T: Real> PairKernel<T> for LennardJones {
33    type Params = (T, T);
34
35    /// Computes only the potential energy.
36    ///
37    /// # Formula
38    ///
39    /// $$ E = D_0 (s^6 - 2s^3), \quad \text{where } s = (R_0/r)^2 $$
40    #[inline(always)]
41    fn energy(r_sq: T, (d0, r0_sq): Self::Params) -> T {
42        let s = r0_sq * r_sq.recip();
43        let s3 = s * s * s;
44        let s6 = s3 * s3;
45
46        d0 * (s6 - T::from(2.0) * s3)
47    }
48
49    /// Computes only the force pre-factor $D$.
50    ///
51    /// # Formula
52    ///
53    /// $$ D = \frac{12 D_0}{r^2} (s^6 - s^3), \quad \text{where } s = (R_0/r)^2 $$
54    ///
55    /// This factor is defined such that the force vector can be computed
56    /// by a single vector multiplication: $\vec{F} = -D \cdot \vec{r}$.
57    #[inline(always)]
58    fn diff(r_sq: T, (d0, r0_sq): Self::Params) -> T {
59        let inv_r2 = r_sq.recip();
60        let s = r0_sq * inv_r2;
61        let s3 = s * s * s;
62        let s6 = s3 * s3;
63
64        T::from(12.0) * d0 * inv_r2 * (s6 - s3)
65    }
66
67    /// Computes both energy and force pre-factor efficiently.
68    ///
69    /// This method reuses intermediate calculations to minimize operations.
70    #[inline(always)]
71    fn compute(r_sq: T, (d0, r0_sq): Self::Params) -> EnergyDiff<T> {
72        let inv_r2 = r_sq.recip();
73        let s = r0_sq * inv_r2;
74        let s3 = s * s * s;
75        let s6 = s3 * s3;
76
77        let energy = d0 * (s6 - T::from(2.0) * s3);
78
79        let diff = T::from(12.0) * d0 * inv_r2 * (s6 - s3);
80
81        EnergyDiff { energy, diff }
82    }
83}
84
85/// Buckingham (Exponential-6) potential for Van der Waals interactions.
86///
87/// # Physics
88///
89/// Models short-range repulsion exponentially and long-range attraction with an $r^{-6}$ term,
90/// providing a more physically accurate repulsion wall than Lennard-Jones.
91///
92/// - **Formula**: $$ E = D_0 \left[ \frac{6}{\zeta-6} \exp\left(\zeta(1 - \frac{r}{R_0})\right) - \frac{\zeta}{\zeta-6} \left(\frac{R_0}{r}\right)^6 \right] $$
93/// - **Derivative Factor (`diff`)**: $$ D = -\frac{1}{r} \frac{dE}{dr} = \frac{6\zeta D_0}{r(\zeta-6)R_0} \left[ \exp\left(\zeta(1 - \frac{r}{R_0})\right) - \left(\frac{R_0}{r}\right)^7 \right] $$
94///
95/// # Parameters
96///
97/// For computational efficiency, the physical parameters ($D_0, R_0, \zeta$) are pre-computed
98/// into the standard Buckingham form ($A, B, C$):
99/// - `a`: The repulsion pre-factor $A = \frac{6 D_0}{\zeta-6} e^{\zeta}$.
100/// - `b`: The repulsion decay constant $B = \zeta / R_0$.
101/// - `c`: The attraction pre-factor $C = \frac{\zeta D_0 R_0^6}{\zeta-6}$.
102/// - `r_fusion_sq`: The squared distance threshold for regularization.
103///
104/// # Inputs
105///
106/// - `r_sq`: Squared distance $r^2$ between two atoms.
107///
108/// # Implementation Notes
109///
110/// - The kernel operates on the computationally efficient $A, B, C$ form.
111/// - A branchless regularization is applied for $r^2 < r_{fusion}^2$ using a mathematical mask
112///   to prevent energy collapse at very short distances, ensuring numerical stability.
113/// - Requires one `sqrt` and one `exp` call, making it computationally more demanding than LJ.
114/// - Power chain `inv_r2 -> inv_r6 -> inv_r8` is used for the attractive term.
115#[derive(Clone, Copy, Debug, Default)]
116pub struct Buckingham;
117
118impl<T: Real> PairKernel<T> for Buckingham {
119    type Params = (T, T, T, T);
120
121    /// Computes only the potential energy.
122    ///
123    /// # Formula
124    ///
125    /// $$ E = A e^{-B r} - C / r^6 $$
126    #[inline(always)]
127    fn energy(r_sq: T, (a, b, c, r_fusion_sq): Self::Params) -> T {
128        let is_safe = T::from((r_sq >= r_fusion_sq) as u8 as f32);
129        let effective_r_sq = r_sq.max(r_fusion_sq);
130
131        let r = effective_r_sq.sqrt();
132        let inv_r2 = effective_r_sq.recip();
133        let inv_r6 = inv_r2 * inv_r2 * inv_r2;
134
135        let repulsion = a * T::exp(-b * r);
136        let attraction = c * inv_r6;
137        let energy_unsafe = repulsion - attraction;
138
139        const FUSION_ENERGY_PENALTY: f32 = 1.0e6;
140        let penalty = T::from(FUSION_ENERGY_PENALTY);
141
142        energy_unsafe * is_safe + penalty * (T::from(1.0) - is_safe)
143    }
144
145    /// Computes only the force pre-factor $D$.
146    ///
147    /// # Formula
148    ///
149    /// $$ D = \frac{A B e^{-B r}}{r} - \frac{6 C}{r^8} $$
150    ///
151    /// This factor is defined such that the force vector can be computed
152    /// by a single vector multiplication: $\vec{F} = -D \cdot \vec{r}$.
153    #[inline(always)]
154    fn diff(r_sq: T, (a, b, c, r_fusion_sq): Self::Params) -> T {
155        let is_safe = T::from((r_sq >= r_fusion_sq) as u8 as f32);
156        let effective_r_sq = r_sq.max(r_fusion_sq);
157
158        let inv_r = effective_r_sq.rsqrt();
159        let r = effective_r_sq * inv_r;
160        let inv_r2 = inv_r * inv_r;
161        let inv_r4 = inv_r2 * inv_r2;
162        let inv_r8 = inv_r4 * inv_r4;
163
164        let exp_term = T::exp(-b * r);
165
166        let repulsion_factor = a * b * exp_term * inv_r;
167        let attraction_factor = T::from(6.0) * c * inv_r8;
168        let diff_unsafe = repulsion_factor - attraction_factor;
169
170        const FUSION_FORCE_PENALTY: f32 = 1.0e6;
171        let penalty = T::from(FUSION_FORCE_PENALTY);
172
173        diff_unsafe * is_safe + penalty * (T::from(1.0) - is_safe)
174    }
175
176    /// Computes both energy and force pre-factor efficiently.
177    ///
178    /// This method reuses intermediate calculations to minimize operations.
179    #[inline(always)]
180    fn compute(r_sq: T, (a, b, c, r_fusion_sq): Self::Params) -> EnergyDiff<T> {
181        let is_safe = T::from((r_sq >= r_fusion_sq) as u8 as f32);
182        let effective_r_sq = r_sq.max(r_fusion_sq);
183
184        let inv_r = effective_r_sq.rsqrt();
185        let r = effective_r_sq * inv_r;
186        let inv_r2 = inv_r * inv_r;
187        let inv_r4 = inv_r2 * inv_r2;
188        let inv_r6 = inv_r4 * inv_r2;
189        let inv_r8 = inv_r6 * inv_r2;
190
191        let exp_term = T::exp(-b * r);
192
193        let repulsion_energy = a * exp_term;
194        let attraction_energy = c * inv_r6;
195        let energy_unsafe = repulsion_energy - attraction_energy;
196
197        let repulsion_force = repulsion_energy * b * inv_r;
198        let attraction_force = T::from(6.0) * c * inv_r8;
199        let diff_unsafe = repulsion_force - attraction_force;
200
201        const FUSION_ENERGY_PENALTY: f32 = 1.0e6;
202        const FUSION_FORCE_PENALTY: f32 = 1.0e6;
203
204        let energy_penalty = T::from(FUSION_ENERGY_PENALTY);
205        let force_penalty = T::from(FUSION_FORCE_PENALTY);
206
207        EnergyDiff {
208            energy: energy_unsafe * is_safe + energy_penalty * (T::from(1.0) - is_safe),
209            diff: diff_unsafe * is_safe + force_penalty * (T::from(1.0) - is_safe),
210        }
211    }
212}
213
214/// Splined Buckingham (Exp-6) potential with $C^2$ continuous short-range regularization.
215///
216/// # Physics
217///
218/// This kernel enhances the standard Buckingham potential by replacing the problematic
219/// short-range region ($r < r_{spline}$) with a quintic (5th degree) polynomial, $P_5(r)$.
220/// This guarantees that the energy, force ($C^1$), and force derivative ($C^2$) are
221/// continuous everywhere, which is critical for stable molecular dynamics simulations.
222///
223/// - **Formula**: $$ E(r) = \begin{cases} D_0 \left[ \frac{6}{\zeta-6} \exp\left(\zeta(1 - \frac{r}{R_0})\right) - \frac{\zeta}{\zeta-6} \left(\frac{R_0}{r}\right)^6 \right] & r \ge r_{spline} \\ P_5(r) = \sum_{i=0}^{5} p_i r^i & r < r_{spline} \end{cases} $$
224/// - **Derivative Factor (`diff`)**: $$ D = -\frac{1}{r} \frac{dE}{dr} = \begin{cases} \frac{6\zeta D_0}{r(\zeta-6)R_0} \left[ \exp\left(\zeta(1 - \frac{r}{R_0})\right) - \left(\frac{R_0}{r}\right)^7 \right] & r \ge r_{spline} \\ -\frac{1}{r} \sum_{i=1}^{5} i \cdot p_i r^{i-1} & r < r_{spline} \end{cases} $$
225///
226/// # Parameters
227///
228/// For computational efficiency, the physical parameters ($D_0, R_0, \zeta$) are pre-computed
229/// into two sets for the long-range and short-range parts:
230///
231/// - **Long-Range Part ($r \ge r_{spline}$)**:
232///   - `a`: The repulsion pre-factor $A = \frac{6 D_0}{\zeta-6} e^{\zeta}$.
233///   - `b`: The repulsion decay constant $B = \zeta / R_0$.
234///   - `c`: The attraction pre-factor $C = \frac{\zeta D_0 R_0^6}{\zeta-6}$.
235///
236/// - **Short-Range Part ($r < r_{spline}$)**:
237///   - `r_spline_sq`: The squared distance threshold for switching to the polynomial.
238///   - `p0..p5`: The six coefficients of the quintic polynomial $P_5(r)$, pre-computed
239///     to satisfy $C^2$ continuity at $r_{spline}$ and boundary conditions at $r=0$.
240///
241/// # Inputs
242///
243/// - `r_sq`: Squared distance $r^2$ between two atoms.
244///
245/// # Implementation Notes
246///
247/// - A branchless selection mechanism is used to switch between the Buckingham and
248///   polynomial forms, making it SIMD-friendly.
249/// - The polynomial is evaluated using Horner's method for improved numerical stability
250///   and reduced floating-point operations.
251/// - The entire computation, including both paths, is executed to avoid pipeline stalls,
252///   making the runtime performance constant and predictable.
253#[derive(Clone, Copy, Debug, Default)]
254pub struct SplinedBuckingham;
255
256impl<T: Real> PairKernel<T> for SplinedBuckingham {
257    type Params = (T, T, T, T, T, T, T, T, T, T);
258
259    /// Computes only the potential energy, selecting between Exp-6 and polynomial forms.
260    ///
261    /// # Formula
262    ///
263    /// $$ E(r) = \begin{cases} A e^{-Br} - C/r^6 & r \ge r_{spline} \\ P_5(r) & r < r_{spline} \end{cases} $$
264    #[inline(always)]
265    fn energy(r_sq: T, params: Self::Params) -> T {
266        let (a, b, c, r_spline_sq, p0, p1, p2, p3, p4, p5) = params;
267
268        let mask = T::from((r_sq >= r_spline_sq) as u8 as f32);
269
270        let r = r_sq.sqrt();
271
272        let r_sq_long = r_sq.max(r_spline_sq);
273        let inv_r2_long = r_sq_long.recip();
274        let inv_r6_long = inv_r2_long * inv_r2_long * inv_r2_long;
275
276        let energy_long = a * T::exp(-b * r) - c * inv_r6_long;
277
278        let energy_short = ((((p5 * r + p4) * r + p3) * r + p2) * r + p1) * r + p0;
279
280        energy_long * mask + energy_short * (T::from(1.0) - mask)
281    }
282
283    /// Computes only the force pre-factor $D$, selecting between Exp-6 and polynomial forms.
284    ///
285    /// # Formula
286    ///
287    /// $$ D(r) = \begin{cases} \frac{A B e^{-B r}}{r} - \frac{6 C}{r^8} & r \ge r_{spline} \\ -\frac{P'_5(r)}{r} & r < r_{spline} \end{cases} $$
288    ///
289    /// This factor is defined such that the force vector can be computed
290    /// by a single vector multiplication: $\vec{F} = -D \cdot \vec{r}$.
291    #[inline(always)]
292    fn diff(r_sq: T, params: Self::Params) -> T {
293        let (a, b, c, r_spline_sq, _, p1, p2, p3, p4, p5) = params;
294
295        let mask = T::from((r_sq >= r_spline_sq) as u8 as f32);
296
297        let r = r_sq.sqrt();
298
299        let r_sq_long = r_sq.max(r_spline_sq);
300        let inv_r_long = r_sq_long.rsqrt();
301        let inv_r2_long = inv_r_long * inv_r_long;
302        let inv_r4_long = inv_r2_long * inv_r2_long;
303        let inv_r8_long = inv_r4_long * inv_r4_long;
304
305        let diff_long = a * b * T::exp(-b * r) * inv_r_long - T::from(6.0) * c * inv_r8_long;
306
307        let r2 = r_sq;
308        let r3 = r2 * r;
309        let r4 = r2 * r2;
310        let poly_deriv = p1
311            + T::from(2.0) * p2 * r
312            + T::from(3.0) * p3 * r2
313            + T::from(4.0) * p4 * r3
314            + T::from(5.0) * p5 * r4;
315        let diff_short = -(poly_deriv * r.recip());
316
317        diff_long * mask + diff_short * (T::from(1.0) - mask)
318    }
319
320    /// Computes both energy and force pre-factor efficiently.
321    ///
322    /// This method reuses intermediate calculations to minimize operations.
323    #[inline(always)]
324    fn compute(r_sq: T, params: Self::Params) -> EnergyDiff<T> {
325        let (a, b, c, r_spline_sq, p0, p1, p2, p3, p4, p5) = params;
326
327        let mask = T::from((r_sq >= r_spline_sq) as u8 as f32);
328
329        let r = r_sq.sqrt();
330
331        let r_sq_long = r_sq.max(r_spline_sq);
332        let inv_r_long = r_sq_long.rsqrt();
333        let inv_r2_long = inv_r_long * inv_r_long;
334        let inv_r4_long = inv_r2_long * inv_r2_long;
335        let inv_r6_long = inv_r4_long * inv_r2_long;
336
337        let exp_term_long = T::exp(-b * r);
338        let energy_long = a * exp_term_long - c * inv_r6_long;
339
340        let inv_r8_long = inv_r4_long * inv_r4_long;
341        let diff_long = a * b * exp_term_long * inv_r_long - T::from(6.0) * c * inv_r8_long;
342
343        let energy_short = ((((p5 * r + p4) * r + p3) * r + p2) * r + p1) * r + p0;
344        let r2 = r_sq;
345        let poly_deriv = p1
346            + T::from(2.0) * p2 * r
347            + T::from(3.0) * p3 * r2
348            + T::from(4.0) * p4 * (r2 * r)
349            + T::from(5.0) * p5 * (r2 * r2);
350        let diff_short = -(poly_deriv * r.recip());
351
352        EnergyDiff {
353            energy: energy_long * mask + energy_short * (T::from(1.0) - mask),
354            diff: diff_long * mask + diff_short * (T::from(1.0) - mask),
355        }
356    }
357}
358
359#[cfg(test)]
360mod tests {
361    use super::*;
362    use approx::assert_relative_eq;
363
364    // ------------------------------------------------------------------------
365    // Test Constants
366    // ------------------------------------------------------------------------
367
368    const H: f64 = 1e-6;
369    const TOL_DIFF: f64 = 1e-4;
370
371    // Typical LJ parameters: D0 = 0.1 kcal/mol, R0 = 3.5 Å
372    const D0: f64 = 0.1;
373    const R0: f64 = 3.5;
374    const R0_SQ: f64 = R0 * R0;
375
376    // Buckingham parameters: A, B, C, r_fusion_sq
377    const BUCK_A: f64 = 1000.0;
378    const BUCK_B: f64 = 3.0;
379    const BUCK_C: f64 = 50.0;
380    const BUCK_FUSION_SQ: f64 = 0.5;
381
382    // ========================================================================
383    // Lennard-Jones Tests
384    // ========================================================================
385
386    mod lennard_jones {
387        use super::*;
388
389        // --------------------------------------------------------------------
390        // 1. Sanity Checks
391        // --------------------------------------------------------------------
392
393        #[test]
394        fn sanity_compute_equals_separate() {
395            let r_sq = 9.0_f64;
396            let params = (D0, R0_SQ);
397
398            let result = LennardJones::compute(r_sq, params);
399            let energy_only = LennardJones::energy(r_sq, params);
400            let diff_only = LennardJones::diff(r_sq, params);
401
402            assert_relative_eq!(result.energy, energy_only, epsilon = 1e-14);
403            assert_relative_eq!(result.diff, diff_only, epsilon = 1e-14);
404        }
405
406        #[test]
407        fn sanity_f32_f64_consistency() {
408            let r_sq_64 = 12.25_f64;
409            let r_sq_32 = 12.25_f32;
410            let params_64 = (D0, R0_SQ);
411            let params_32 = (D0 as f32, R0_SQ as f32);
412
413            let e64 = LennardJones::energy(r_sq_64, params_64);
414            let e32 = LennardJones::energy(r_sq_32, params_32);
415
416            assert_relative_eq!(e64, e32 as f64, epsilon = 1e-5);
417        }
418
419        #[test]
420        fn sanity_equilibrium_energy_minimum() {
421            let e = LennardJones::energy(R0_SQ, (D0, R0_SQ));
422            assert_relative_eq!(e, -D0, epsilon = 1e-10);
423        }
424
425        #[test]
426        fn sanity_equilibrium_zero_force() {
427            let d = LennardJones::diff(R0_SQ, (D0, R0_SQ));
428            assert_relative_eq!(d, 0.0, epsilon = 1e-10);
429        }
430
431        // --------------------------------------------------------------------
432        // 2. Numerical Stability
433        // --------------------------------------------------------------------
434
435        #[test]
436        fn stability_large_distance() {
437            let r_sq = 1e6_f64;
438            let result = LennardJones::compute(r_sq, (D0, R0_SQ));
439
440            assert!(result.energy.is_finite());
441            assert!(result.diff.is_finite());
442            assert!(result.energy.abs() < 1e-10);
443        }
444
445        #[test]
446        fn stability_small_distance() {
447            let r_sq = 1.0_f64;
448            let result = LennardJones::compute(r_sq, (D0, R0_SQ));
449
450            assert!(result.energy.is_finite());
451            assert!(result.diff.is_finite());
452            assert!(result.energy > 0.0);
453        }
454
455        // --------------------------------------------------------------------
456        // 3. Finite Difference Verification
457        // --------------------------------------------------------------------
458
459        fn finite_diff_check(r: f64, params: (f64, f64)) {
460            let r_sq = r * r;
461
462            let r_plus = r + H;
463            let r_minus = r - H;
464            let e_plus = LennardJones::energy(r_plus * r_plus, params);
465            let e_minus = LennardJones::energy(r_minus * r_minus, params);
466            let de_dr_numerical = (e_plus - e_minus) / (2.0 * H);
467
468            let d_analytic = LennardJones::diff(r_sq, params);
469            let de_dr_analytic = -d_analytic * r;
470
471            assert_relative_eq!(de_dr_numerical, de_dr_analytic, epsilon = TOL_DIFF);
472        }
473
474        #[test]
475        fn finite_diff_repulsion_region() {
476            finite_diff_check(2.5, (D0, R0_SQ));
477        }
478
479        #[test]
480        fn finite_diff_equilibrium_region() {
481            finite_diff_check(R0, (D0, R0_SQ));
482        }
483
484        #[test]
485        fn finite_diff_attraction_region() {
486            finite_diff_check(5.0, (D0, R0_SQ));
487        }
488
489        #[test]
490        fn finite_diff_long_range() {
491            finite_diff_check(10.0, (D0, R0_SQ));
492        }
493
494        // --------------------------------------------------------------------
495        // 4. LJ-Specific: Physical Behavior
496        // --------------------------------------------------------------------
497
498        #[test]
499        fn specific_repulsion_positive_energy() {
500            let e = LennardJones::energy(4.0, (D0, R0_SQ));
501            assert!(e > 0.0);
502        }
503
504        #[test]
505        fn specific_attraction_negative_energy() {
506            let e = LennardJones::energy(25.0, (D0, R0_SQ));
507            assert!(e < 0.0);
508        }
509
510        #[test]
511        fn specific_diff_sign_repulsion() {
512            let d = LennardJones::diff(4.0, (D0, R0_SQ));
513            assert!(d > 0.0);
514        }
515
516        #[test]
517        fn specific_diff_sign_attraction() {
518            let d = LennardJones::diff(25.0, (D0, R0_SQ));
519            assert!(d < 0.0);
520        }
521    }
522
523    // ========================================================================
524    // Buckingham Tests
525    // ========================================================================
526
527    mod buckingham {
528        use super::*;
529
530        fn params() -> (f64, f64, f64, f64) {
531            (BUCK_A, BUCK_B, BUCK_C, BUCK_FUSION_SQ)
532        }
533
534        // --------------------------------------------------------------------
535        // 1. Sanity Checks
536        // --------------------------------------------------------------------
537
538        #[test]
539        fn sanity_compute_equals_separate() {
540            let r_sq = 4.0_f64;
541            let p = params();
542
543            let result = Buckingham::compute(r_sq, p);
544            let energy_only = Buckingham::energy(r_sq, p);
545            let diff_only = Buckingham::diff(r_sq, p);
546
547            assert_relative_eq!(result.energy, energy_only, epsilon = 1e-12);
548            assert_relative_eq!(result.diff, diff_only, epsilon = 1e-12);
549        }
550
551        #[test]
552        fn sanity_f32_f64_consistency() {
553            let r_sq = 4.0;
554            let p64 = params();
555            let p32 = (
556                BUCK_A as f32,
557                BUCK_B as f32,
558                BUCK_C as f32,
559                BUCK_FUSION_SQ as f32,
560            );
561
562            let e64 = Buckingham::energy(r_sq, p64);
563            let e32 = Buckingham::energy(r_sq as f32, p32);
564
565            assert_relative_eq!(e64, e32 as f64, epsilon = 1e-3);
566        }
567
568        // --------------------------------------------------------------------
569        // 2. Numerical Stability
570        // --------------------------------------------------------------------
571
572        #[test]
573        fn stability_fusion_region() {
574            let r_sq = 0.1_f64;
575            let result = Buckingham::compute(r_sq, params());
576
577            assert!(result.energy.is_finite());
578            assert!(result.diff.is_finite());
579            assert!(result.energy > 1e5);
580        }
581
582        #[test]
583        fn stability_large_distance() {
584            let r_sq = 1e4_f64;
585            let result = Buckingham::compute(r_sq, params());
586
587            assert!(result.energy.is_finite());
588            assert!(result.diff.is_finite());
589        }
590
591        // --------------------------------------------------------------------
592        // 3. Finite Difference Verification
593        // --------------------------------------------------------------------
594
595        fn finite_diff_check(r: f64) {
596            let p = params();
597            let r_sq = r * r;
598
599            let r_plus = r + H;
600            let r_minus = r - H;
601            let e_plus = Buckingham::energy(r_plus * r_plus, p);
602            let e_minus = Buckingham::energy(r_minus * r_minus, p);
603            let de_dr_numerical = (e_plus - e_minus) / (2.0 * H);
604
605            let d_analytic = Buckingham::diff(r_sq, p);
606            let de_dr_analytic = -d_analytic * r;
607
608            assert_relative_eq!(de_dr_numerical, de_dr_analytic, epsilon = TOL_DIFF);
609        }
610
611        #[test]
612        fn finite_diff_short_range() {
613            finite_diff_check(1.5);
614        }
615
616        #[test]
617        fn finite_diff_medium_range() {
618            finite_diff_check(3.0);
619        }
620
621        #[test]
622        fn finite_diff_long_range() {
623            finite_diff_check(8.0);
624        }
625
626        // --------------------------------------------------------------------
627        // 4. Buckingham-Specific
628        // --------------------------------------------------------------------
629
630        #[test]
631        fn specific_exponential_dominates_short_range() {
632            let e1 = Buckingham::energy(0.81, params());
633            let e2 = Buckingham::energy(1.0, params());
634            assert!(e1.is_finite());
635            assert!(e2.is_finite());
636        }
637    }
638
639    // ========================================================================
640    // SplinedBuckingham Tests
641    // ========================================================================
642
643    mod splined_buckingham {
644        use super::*;
645
646        fn params() -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) {
647            let a = 1000.0;
648            let b = 3.0;
649            let c = 50.0;
650            let r_spline_sq = 1.0;
651
652            let p0 = 100.0;
653            let p1 = -50.0;
654            let p2 = 10.0;
655            let p3 = 0.0;
656            let p4 = 0.0;
657            let p5 = 0.0;
658
659            (a, b, c, r_spline_sq, p0, p1, p2, p3, p4, p5)
660        }
661
662        // --------------------------------------------------------------------
663        // 1. Sanity Checks
664        // --------------------------------------------------------------------
665
666        #[test]
667        fn sanity_compute_equals_separate() {
668            let r_sq = 4.0_f64;
669            let p = params();
670
671            let result = SplinedBuckingham::compute(r_sq, p);
672            let energy_only = SplinedBuckingham::energy(r_sq, p);
673            let diff_only = SplinedBuckingham::diff(r_sq, p);
674
675            assert_relative_eq!(result.energy, energy_only, epsilon = 1e-12);
676            assert_relative_eq!(result.diff, diff_only, epsilon = 1e-12);
677        }
678
679        // --------------------------------------------------------------------
680        // 2. Numerical Stability
681        // --------------------------------------------------------------------
682
683        #[test]
684        fn stability_inside_spline_region() {
685            let r_sq = 0.25_f64;
686            let result = SplinedBuckingham::compute(r_sq, params());
687
688            assert!(result.energy.is_finite());
689            assert!(result.diff.is_finite());
690        }
691
692        #[test]
693        fn stability_outside_spline_region() {
694            let r_sq = 4.0_f64;
695            let result = SplinedBuckingham::compute(r_sq, params());
696
697            assert!(result.energy.is_finite());
698            assert!(result.diff.is_finite());
699        }
700
701        // --------------------------------------------------------------------
702        // 3. Finite Difference Verification
703        // --------------------------------------------------------------------
704
705        fn finite_diff_check(r: f64) {
706            let p = params();
707            let r_sq = r * r;
708
709            let r_plus = r + H;
710            let r_minus = r - H;
711            let e_plus = SplinedBuckingham::energy(r_plus * r_plus, p);
712            let e_minus = SplinedBuckingham::energy(r_minus * r_minus, p);
713            let de_dr_numerical = (e_plus - e_minus) / (2.0 * H);
714
715            let d_analytic = SplinedBuckingham::diff(r_sq, p);
716            let de_dr_analytic = -d_analytic * r;
717
718            assert_relative_eq!(de_dr_numerical, de_dr_analytic, epsilon = TOL_DIFF);
719        }
720
721        #[test]
722        fn finite_diff_inside_spline() {
723            finite_diff_check(0.5);
724        }
725
726        #[test]
727        fn finite_diff_outside_spline() {
728            finite_diff_check(2.0);
729        }
730
731        #[test]
732        fn finite_diff_long_range() {
733            finite_diff_check(5.0);
734        }
735
736        // --------------------------------------------------------------------
737        // 4. Spline-Specific: Region Selection
738        // --------------------------------------------------------------------
739
740        #[test]
741        fn specific_uses_polynomial_inside() {
742            let p = params();
743            let r_sq = 0.25_f64;
744
745            let r = r_sq.sqrt();
746            let expected = p.4 + p.5 * r + p.6 * r * r;
747            let actual = SplinedBuckingham::energy(r_sq, p);
748
749            assert_relative_eq!(actual, expected, epsilon = 1e-10);
750        }
751
752        fn c2_continuous_params() -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) {
753            let a = 1000.0_f64;
754            let b = 3.0_f64;
755            let c = 50.0_f64;
756            let r_s = 1.2_f64;
757            let r_spline_sq = r_s * r_s;
758
759            let exp_term = (-b * r_s).exp();
760            let inv_r = 1.0 / r_s;
761            let inv_r2 = inv_r * inv_r;
762            let inv_r6 = inv_r2 * inv_r2 * inv_r2;
763            let inv_r7 = inv_r6 * inv_r;
764            let inv_r8 = inv_r6 * inv_r2;
765
766            let e_s = a * exp_term - c * inv_r6;
767            let de_s = -a * b * exp_term + 6.0 * c * inv_r7;
768            let d2e_s = a * b * b * exp_term - 42.0 * c * inv_r8;
769
770            let e_max = 1e4_f64;
771            let p0 = e_max;
772            let p1 = 0.0;
773            let p5 = 0.0;
774
775            let r2 = r_s * r_s;
776            let r3 = r2 * r_s;
777            let r4 = r2 * r2;
778
779            let rhs1 = e_s - p0;
780            let rhs2 = de_s;
781            let rhs3 = d2e_s;
782
783            let det = r2 * (3.0 * r2 * 12.0 * r2 - 4.0 * r3 * 6.0 * r_s)
784                - r3 * (2.0 * r_s * 12.0 * r2 - 4.0 * r3 * 2.0)
785                + r4 * (2.0 * r_s * 6.0 * r_s - 3.0 * r2 * 2.0);
786
787            let det_p2 = rhs1 * (3.0 * r2 * 12.0 * r2 - 4.0 * r3 * 6.0 * r_s)
788                - r3 * (rhs2 * 12.0 * r2 - 4.0 * r3 * rhs3)
789                + r4 * (rhs2 * 6.0 * r_s - 3.0 * r2 * rhs3);
790
791            let det_p3 = r2 * (rhs2 * 12.0 * r2 - 4.0 * r3 * rhs3)
792                - rhs1 * (2.0 * r_s * 12.0 * r2 - 4.0 * r3 * 2.0)
793                + r4 * (2.0 * r_s * rhs3 - rhs2 * 2.0);
794
795            let det_p4 = r2 * (3.0 * r2 * rhs3 - rhs2 * 6.0 * r_s)
796                - r3 * (2.0 * r_s * rhs3 - rhs2 * 2.0)
797                + rhs1 * (2.0 * r_s * 6.0 * r_s - 3.0 * r2 * 2.0);
798
799            let p2 = det_p2 / det;
800            let p3 = det_p3 / det;
801            let p4 = det_p4 / det;
802
803            (a, b, c, r_spline_sq, p0, p1, p2, p3, p4, p5)
804        }
805
806        #[test]
807        fn verify_spline_coefficients() {
808            let p = c2_continuous_params();
809            let (a, b, c, r_spline_sq, p0, p1, p2, p3, p4, p5) = p;
810            let r_s = r_spline_sq.sqrt();
811
812            let exp_term = (-b * r_s).exp();
813            let inv_r = 1.0 / r_s;
814            let inv_r2 = inv_r * inv_r;
815            let inv_r6 = inv_r2 * inv_r2 * inv_r2;
816            let inv_r7 = inv_r6 * inv_r;
817            let inv_r8 = inv_r6 * inv_r2;
818
819            let e_buck = a * exp_term - c * inv_r6;
820            let de_buck = -a * b * exp_term + 6.0 * c * inv_r7;
821            let d2e_buck = a * b * b * exp_term - 42.0 * c * inv_r8;
822
823            let r2 = r_s * r_s;
824            let r3 = r2 * r_s;
825            let r4 = r2 * r2;
826            let r5 = r4 * r_s;
827
828            let e_poly = p0 + p1 * r_s + p2 * r2 + p3 * r3 + p4 * r4 + p5 * r5;
829            let de_poly = p1 + 2.0 * p2 * r_s + 3.0 * p3 * r2 + 4.0 * p4 * r3 + 5.0 * p5 * r4;
830            let d2e_poly = 2.0 * p2 + 6.0 * p3 * r_s + 12.0 * p4 * r2 + 20.0 * p5 * r3;
831
832            assert_relative_eq!(e_poly, e_buck, epsilon = 1e-8);
833            assert_relative_eq!(de_poly, de_buck, epsilon = 1e-8);
834            assert_relative_eq!(d2e_poly, d2e_buck, epsilon = 1e-8);
835        }
836
837        #[test]
838        fn continuity_c0_energy_at_boundary() {
839            let p = c2_continuous_params();
840            let r_spline_sq = p.3;
841            let r_s = r_spline_sq.sqrt();
842            let eps = 1e-8;
843
844            let r_inside = r_s - eps;
845            let r_outside = r_s + eps;
846
847            let e_inside = SplinedBuckingham::energy(r_inside * r_inside, p);
848            let e_outside = SplinedBuckingham::energy(r_outside * r_outside, p);
849
850            assert_relative_eq!(e_inside, e_outside, epsilon = 1e-4);
851        }
852
853        #[test]
854        fn continuity_c1_force_at_boundary() {
855            let p = c2_continuous_params();
856            let r_spline_sq = p.3;
857            let r_s = r_spline_sq.sqrt();
858            let eps = 1e-8;
859
860            let r_inside = r_s - eps;
861            let r_outside = r_s + eps;
862
863            let d_inside = SplinedBuckingham::diff(r_inside * r_inside, p);
864            let d_outside = SplinedBuckingham::diff(r_outside * r_outside, p);
865
866            let de_dr_inside = -d_inside * r_inside;
867            let de_dr_outside = -d_outside * r_outside;
868
869            assert_relative_eq!(de_dr_inside, de_dr_outside, epsilon = 1e-3);
870        }
871
872        #[test]
873        fn continuity_c2_second_derivative_at_boundary() {
874            let p = c2_continuous_params();
875            let (a, b, c, r_spline_sq, _p0, _p1, p2, p3, p4, p5) = p;
876            let r_s = r_spline_sq.sqrt();
877
878            let exp_term = (-b * r_s).exp();
879            let inv_r8 = 1.0 / r_s.powi(8);
880            let d2e_buck_analytical = a * b * b * exp_term - 42.0 * c * inv_r8;
881
882            let d2e_poly_analytical =
883                2.0 * p2 + 6.0 * p3 * r_s + 12.0 * p4 * r_s * r_s + 20.0 * p5 * r_s.powi(3);
884
885            assert_relative_eq!(d2e_poly_analytical, d2e_buck_analytical, epsilon = 1e-6);
886
887            let h = 1e-6;
888
889            let e_m = SplinedBuckingham::energy((r_s - h).powi(2), p);
890            let e_0 = SplinedBuckingham::energy(r_s.powi(2), p);
891            let e_p = SplinedBuckingham::energy((r_s + h).powi(2), p);
892            let d2e_numerical_straddling = (e_p - 2.0 * e_0 + e_m) / (h * h);
893
894            let relative_error =
895                (d2e_numerical_straddling - d2e_buck_analytical).abs() / d2e_buck_analytical.abs();
896            assert!(
897                relative_error < 0.1,
898                "C² continuity check: relative error {} > 0.1",
899                relative_error
900            );
901        }
902
903        #[test]
904        fn continuity_finite_diff_across_boundary() {
905            let p = c2_continuous_params();
906            let r_s = p.3.sqrt();
907
908            let r = r_s;
909            let r_sq = r * r;
910
911            let r_plus = r + H;
912            let r_minus = r - H;
913            let e_plus = SplinedBuckingham::energy(r_plus * r_plus, p);
914            let e_minus = SplinedBuckingham::energy(r_minus * r_minus, p);
915            let de_dr_numerical = (e_plus - e_minus) / (2.0 * H);
916
917            let d_analytic = SplinedBuckingham::diff(r_sq, p);
918            let de_dr_analytic = -d_analytic * r;
919
920            assert_relative_eq!(de_dr_numerical, de_dr_analytic, epsilon = TOL_DIFF);
921        }
922    }
923}
dreid_kernel/potentials/nonbonded/vdw.rs

dreid_kernel/potentials/nonbonded/
vdw.rs
