potentials/meta/
softcore.rs

1//! # Softcore Potential Wrapper
2//!
3//! Removes the singularity at r=0 for use in Free Energy Perturbation (FEP)
4//! and thermodynamic integration calculations.
5//!
6//! ## Formula
7//!
8//! ```text
9//! r_eff^2 = alpha * sigma^2 * lambda^n + r^2
10//!
11//! V_sc(r) = lambda^m * V(r_eff)
12//! ```
13//!
14//! Common choices:
15//! - alpha = 0.5
16//! - n = 1, m = 1 (linear coupling)
17//! - sigma = Lennard-Jones sigma or typical interaction range
18//!
19//! ## Properties
20//!
21//! - At lambda=0: r_eff = sqrt(alpha)*sigma, potential decoupled
22//! - At lambda=1: r_eff = r, original potential recovered
23//! - Finite energy at r=0 for all lambda > 0
24//!
25//! ## References
26//!
27//! Beutler et al., Chem. Phys. Lett. 222, 529 (1994)
28
29use crate::base::Potential2;
30use crate::math::Vector;
31
32/// Softcore potential for FEP calculations.
33///
34/// ## Type Parameters
35///
36/// - `P`: The underlying potential type
37/// - `T`: The vector type
38#[derive(Clone, Copy, Debug)]
39pub struct Softcore<P, T> {
40    inner: P,
41    alpha_sigma_sq: T, // alpha * sigma^2
42    lambda: T,
43    lambda_n: T, // lambda^n
44    lambda_m: T, // lambda^m
45}
46
47impl<P, T: Vector> Softcore<P, T> {
48    /// Creates a new softcore potential.
49    ///
50    /// ## Arguments
51    ///
52    /// - `inner`: The potential to modify
53    /// - `alpha`: Softcore parameter (typically 0.5)
54    /// - `sigma`: Length scale for softcore (typically LJ sigma) (length)
55    /// - `lambda`: Coupling parameter (0 = decoupled, 1 = coupled)
56    /// - `n`: Exponent for distance modification (typically 1)
57    /// - `m`: Exponent for energy scaling (typically 1)
58    ///
59    /// ## Example
60    ///
61    /// ```
62    /// use potentials::{pair::Lj, meta::Softcore};
63    ///
64    /// let lj = Lj::<f64>::new(1.0, 3.4);
65    /// let lj_soft: Softcore<_, f64> = Softcore::new(lj, 0.5, 3.4, 0.5, 1, 1);
66    /// ```
67    #[inline]
68    pub fn new(inner: P, alpha: f64, sigma: f64, lambda: f64, n: u32, m: u32) -> Self {
69        let lambda_n = lambda.powi(n as i32);
70        let lambda_m = lambda.powi(m as i32);
71
72        Self {
73            inner,
74            alpha_sigma_sq: T::splat(alpha * sigma * sigma),
75            lambda: T::splat(lambda),
76            lambda_n: T::splat(lambda_n),
77            lambda_m: T::splat(lambda_m),
78        }
79    }
80
81    /// Updates the lambda value.
82    ///
83    /// Use this during FEP simulations to change the coupling strength.
84    #[inline]
85    pub fn set_lambda(&mut self, lambda: f64, n: u32, m: u32) {
86        self.lambda = T::splat(lambda);
87        self.lambda_n = T::splat(lambda.powi(n as i32));
88        self.lambda_m = T::splat(lambda.powi(m as i32));
89    }
90
91    /// Computes the effective r^2 for softcore interaction.
92    #[inline(always)]
93    fn effective_r_sq(&self, r_sq: T) -> T {
94        self.alpha_sigma_sq * self.lambda_n + r_sq
95    }
96}
97
98impl<P: Potential2<T>, T: Vector> Potential2<T> for Softcore<P, T> {
99    #[inline(always)]
100    fn energy(&self, r_sq: T) -> T {
101        let r_eff_sq = self.effective_r_sq(r_sq);
102        self.lambda_m * self.inner.energy(r_eff_sq)
103    }
104
105    #[inline(always)]
106    fn force_factor(&self, r_sq: T) -> T {
107        let r_eff_sq = self.effective_r_sq(r_sq);
108        self.lambda_m * self.inner.force_factor(r_eff_sq)
109    }
110
111    #[inline(always)]
112    fn energy_force(&self, r_sq: T) -> (T, T) {
113        let r_eff_sq = self.effective_r_sq(r_sq);
114        let (e, f) = self.inner.energy_force(r_eff_sq);
115        (self.lambda_m * e, self.lambda_m * f)
116    }
117}
118
119#[cfg(test)]
120mod tests {
121    use super::*;
122    use crate::pair::Lj;
123    use approx::assert_relative_eq;
124
125    #[test]
126    fn test_softcore_lambda_one() {
127        let lj: Lj<f64> = Lj::new(1.0, 3.4);
128        let lj_soft = Softcore::new(lj, 0.5, 3.4, 1.0, 1, 1);
129
130        let r_sq = 16.0;
131
132        let (e, f) = lj_soft.energy_force(r_sq);
133        let r_eff_sq = 0.5 * 3.4 * 3.4 + r_sq;
134        let (e_expected, f_expected) = lj.energy_force(r_eff_sq);
135
136        assert_relative_eq!(e, e_expected, epsilon = 1e-10);
137        assert_relative_eq!(f, f_expected, epsilon = 1e-10);
138    }
139
140    #[test]
141    fn test_softcore_lambda_zero() {
142        let lj: Lj<f64> = Lj::new(1.0, 3.4);
143        let lj_soft = Softcore::new(lj, 0.5, 3.4, 0.0, 1, 1);
144
145        let r_sq = 16.0;
146        let e = lj_soft.energy(r_sq);
147
148        assert_relative_eq!(e, 0.0, epsilon = 1e-10);
149    }
150
151    #[test]
152    fn test_softcore_finite_at_zero() {
153        let lj: Lj<f64> = Lj::new(1.0, 3.4);
154        let lj_soft = Softcore::new(lj, 0.5, 3.4, 0.5, 1, 1);
155
156        let e = lj_soft.energy(0.0);
157
158        assert!(e.is_finite());
159        assert!(e.abs() < 1e10);
160    }
161
162    #[test]
163    fn test_softcore_numerical_derivative() {
164        let lj: Lj<f64> = Lj::new(1.0, 3.4);
165        let lj_soft = Softcore::new(lj, 0.5, 3.4, 0.5, 1, 1);
166
167        let r = 4.0;
168
169        let h = 1e-7;
170        let e_plus = lj_soft.energy((r + h) * (r + h));
171        let e_minus = lj_soft.energy((r - h) * (r - h));
172        let dv_dr_numerical = (e_plus - e_minus) / (2.0 * h);
173
174        let f = lj_soft.force_factor(r * r);
175        let dv_dr_analytical = -f * r;
176
177        assert_relative_eq!(dv_dr_analytical, dv_dr_numerical, epsilon = 1e-6);
178    }
179}
potentials/meta/softcore.rs

potentials/meta/
softcore.rs
