potentials/meta/
cutoff.rs

1//! # Hard Cutoff Wrapper
2//!
3//! Sets potential to zero beyond a cutoff distance.
4//!
5//! ## Formula
6//!
7//! ```text
8//! V_cut(r) = V(r) if r < rc
9//!          = 0    if r >= rc
10//! ```
11//!
12//! ## Notes
13//!
14//! - Simple truncation (discontinuous at rc)
15//! - For smooth cutoffs, use [`Switch`] or [`Shift`]
16//! - Implemented branchlessly using mask operations
17
18use crate::base::Potential2;
19use crate::math::{Mask, Vector};
20
21/// Hard cutoff wrapper.
22///
23/// ## Type Parameters
24///
25/// - `P`: The underlying potential type
26/// - `T`: The vector type
27#[derive(Clone, Copy, Debug)]
28pub struct Cutoff<P, T> {
29    inner: P,
30    rc_sq: T,
31}
32
33impl<P, T: Vector> Cutoff<P, T> {
34    /// Creates a new cutoff wrapper.
35    ///
36    /// ## Arguments
37    ///
38    /// - `inner`: The potential to wrap
39    /// - `rc`: Cutoff distance (length)
40    ///
41    /// ## Example
42    ///
43    /// ```
44    /// use potentials::{pair::Lj, meta::Cutoff};
45    ///
46    /// let lj = Lj::<f64>::new(1.0, 3.4);
47    /// let lj_cut: Cutoff<_, f64> = Cutoff::new(lj, 12.0);
48    /// ```
49    #[inline]
50    pub fn new(inner: P, rc: f64) -> Self {
51        Self {
52            inner,
53            rc_sq: T::splat(rc * rc),
54        }
55    }
56
57    /// Returns a reference to the inner potential.
58    #[inline]
59    pub fn inner(&self) -> &P {
60        &self.inner
61    }
62}
63
64impl<P: Potential2<T>, T: Vector> Potential2<T> for Cutoff<P, T> {
65    #[inline(always)]
66    fn energy(&self, r_sq: T) -> T {
67        let e = self.inner.energy(r_sq);
68        let mask = r_sq.lt(self.rc_sq);
69        mask.select(e, T::zero())
70    }
71
72    #[inline(always)]
73    fn force_factor(&self, r_sq: T) -> T {
74        let s = self.inner.force_factor(r_sq);
75        let mask = r_sq.lt(self.rc_sq);
76        mask.select(s, T::zero())
77    }
78
79    #[inline(always)]
80    fn energy_force(&self, r_sq: T) -> (T, T) {
81        let (e, s) = self.inner.energy_force(r_sq);
82        let mask = r_sq.lt(self.rc_sq);
83        (mask.select(e, T::zero()), mask.select(s, T::zero()))
84    }
85}
86
87#[cfg(test)]
88mod tests {
89    use super::*;
90    use crate::pair::Lj;
91    use approx::assert_relative_eq;
92
93    #[test]
94    fn test_cutoff_inside() {
95        let lj: Lj<f64> = Lj::new(1.0, 3.4);
96        let lj_cut = Cutoff::new(lj, 12.0);
97
98        let r_sq = 16.0;
99
100        let e_base = lj.energy(r_sq);
101        let e_cut = lj_cut.energy(r_sq);
102
103        assert_relative_eq!(e_base, e_cut, epsilon = 1e-10);
104    }
105
106    #[test]
107    fn test_cutoff_outside() {
108        let lj: Lj<f64> = Lj::new(1.0, 3.4);
109        let lj_cut = Cutoff::new(lj, 10.0);
110
111        let r_sq = 121.0;
112
113        let e = lj_cut.energy(r_sq);
114        assert_relative_eq!(e, 0.0, epsilon = 1e-10);
115
116        let s = lj_cut.force_factor(r_sq);
117        assert_relative_eq!(s, 0.0, epsilon = 1e-10);
118    }
119
120    #[test]
121    fn test_cutoff_at_boundary() {
122        let lj: Lj<f64> = Lj::new(1.0, 3.4);
123        let rc = 10.0;
124        let lj_cut = Cutoff::new(lj, rc);
125
126        let r_in = rc - 0.001;
127        let e_in = lj_cut.energy(r_in * r_in);
128        assert!(e_in != 0.0);
129
130        let r_out = rc + 0.001;
131        let e_out = lj_cut.energy(r_out * r_out);
132        assert_relative_eq!(e_out, 0.0, epsilon = 1e-10);
133    }
134}
potentials/meta/cutoff.rs

potentials/meta/
cutoff.rs
