chematic_ewald/
lib.rs

1//! Smooth Particle Mesh Ewald (SPME) for long-range electrostatic interactions.
2//!
3//! Provides efficient computation of long-range Coulomb forces in periodic and non-periodic systems.
4//! Uses FFT-based reciprocal space summation combined with real-space cutoff.
5
6#![forbid(unsafe_code)]
7
8use std::fmt;
9
10pub mod pme;
11pub mod real;
12
13pub use pme::reciprocal_space_energy;
14pub use real::{K_COULOMB, direct_coulomb, direct_coulomb_cutoff, direct_coulomb_damped};
15
16/// Error type for SPME calculations.
17#[derive(Debug, Clone, PartialEq)]
18pub enum EwaldError {
19    /// Box matrix is singular or near-singular (determinant < 1e-10).
20    SingularBoxMatrix,
21}
22
23impl fmt::Display for EwaldError {
24    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
25        match self {
26            EwaldError::SingularBoxMatrix => {
27                write!(f, "singular or near-singular box matrix (determinant < 1e-10)")
28            }
29        }
30    }
31}
32
33impl std::error::Error for EwaldError {}
34
35/// Configuration for PME/SPME calculations.
36#[derive(Clone, Debug)]
37pub struct PmeConfig {
38    /// Ewald splitting parameter α (default: auto-computed from r_cut).
39    pub alpha: f64,
40    /// Reciprocal space cutoff: kmax per axis (default: [5, 5, 5]).
41    pub kmax: [usize; 3],
42    /// Real-space cutoff in Ångströms (default: 8.0).
43    pub r_cut: f64,
44    /// PME grid dimensions (default: [16, 16, 16]).
45    pub mesh: [usize; 3],
46    /// B-spline interpolation order (default: 4).
47    pub spline_order: u8,
48}
49
50impl Default for PmeConfig {
51    fn default() -> Self {
52        Self {
53            alpha: 0.0, // Will be auto-computed
54            kmax: [5, 5, 5],
55            r_cut: 8.0,
56            mesh: [16, 16, 16],
57            spline_order: 4,
58        }
59    }
60}
61
62/// Periodic box vectors for MD simulation (3x3 matrix, rows are box vectors).
63#[derive(Clone, Debug)]
64pub struct BoxVectors(pub [[f64; 3]; 3]);
65
66impl BoxVectors {
67    /// Create a cubic box with side length `a` Ångströms.
68    pub fn cubic(a: f64) -> Self {
69        BoxVectors([[a, 0.0, 0.0], [0.0, a, 0.0], [0.0, 0.0, a]])
70    }
71
72    /// Compute box volume (Ų).
73    pub fn volume(&self) -> f64 {
74        let a = self.0[0];
75        let b = self.0[1];
76        let c = self.0[2];
77        (a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0])
78            + a[2] * (b[0] * c[1] - b[1] * c[0]))
79            .abs()
80    }
81}
82
83/// Result of Ewald energy calculation.
84#[derive(Clone, Debug, PartialEq)]
85pub struct EwaldResult {
86    /// Real-space Coulomb energy (kcal/mol).
87    pub real_energy: f64,
88    /// Reciprocal-space (PME) contribution (kcal/mol).
89    pub reciprocal_energy: f64,
90    /// Self-energy correction for Gaussian clouds.
91    pub self_energy: f64,
92    /// Total Ewald energy = real + reciprocal + self.
93    pub total_energy: f64,
94}
95
96/// SPME Ewald energy for periodic system.
97///
98/// Computes long-range electrostatic energy using Smooth Particle Mesh Ewald (SPME).
99/// Splits calculation into real-space (short-range, direct sum with damping) and
100/// reciprocal-space (long-range, FFT-based).
101///
102/// # Arguments
103/// * `coords` - Atomic coordinates [[x0, y0, z0], [x1, y1, z1], ...]
104/// * `charges` - Partial charges (e)
105/// * `box_vecs` - Periodic box vectors
106/// * `config` - PME configuration (alpha, kmax, mesh, spline_order)
107///
108/// # Returns
109/// `Result<EwaldResult, EwaldError>` with breakdown:
110/// - real_energy: Short-range damped Coulomb (r < r_cut)
111/// - reciprocal_energy: Long-range FFT-based contribution
112/// - self_energy: Self-interaction correction
113/// - total_energy: Sum of all contributions
114///
115/// # Errors
116/// Returns `EwaldError::SingularBoxMatrix` if the box matrix is singular.
117///
118/// # Notes
119/// - Requires coordinates to be within the periodic box
120/// - Auto-computes alpha if not provided (α = 3.5 / r_cut)
121/// - Accurate for N up to several thousand atoms with standard mesh
122pub fn spme_energy(
123    coords: &[[f64; 3]],
124    charges: &[f64],
125    box_vecs: &BoxVectors,
126    config: &PmeConfig,
127) -> Result<EwaldResult, EwaldError> {
128    if box_vecs.volume().abs() < 1e-10 {
129        return Err(EwaldError::SingularBoxMatrix);
130    }
131
132    let alpha = if config.alpha > 0.0 {
133        config.alpha
134    } else {
135        3.5 / config.r_cut
136    };
137
138    // Real-space: short-range damped Coulomb
139    let real = real::direct_coulomb_damped(coords, charges, alpha);
140
141    // Reciprocal-space: FFT-based long-range
142    let reciprocal = pme::reciprocal_space_energy(coords, charges, box_vecs, config)?;
143
144    // Self-energy: correction for Gaussian charge distribution
145    let self_corr = compute_self_energy(charges, alpha);
146
147    Ok(EwaldResult {
148        real_energy: real,
149        reciprocal_energy: reciprocal,
150        self_energy: self_corr,
151        total_energy: real + reciprocal + self_corr,
152    })
153}
154
155/// Helper: Direct Coulomb with cutoff using flat coordinate array.
156#[allow(dead_code)]
157fn direct_coulomb_cutoff_coords(coords: &[[f64; 3]], charges: &[f64], r_cut: f64) -> f64 {
158    let mut energy = 0.0;
159    let n = coords.len();
160
161    for i in 0..n {
162        for j in (i + 1)..n {
163            let dx = coords[i][0] - coords[j][0];
164            let dy = coords[i][1] - coords[j][1];
165            let dz = coords[i][2] - coords[j][2];
166            let r = (dx * dx + dy * dy + dz * dz).sqrt();
167
168            if r > 1e-6 && r < r_cut {
169                energy += K_COULOMB * charges[i] * charges[j] / r;
170            }
171        }
172    }
173
174    energy
175}
176
177fn compute_self_energy(charges: &[f64], alpha: f64) -> f64 {
178    let alpha = if alpha > 0.0 { alpha } else { 3.5 / 8.0 };
179    let sqrt_pi = std::f64::consts::PI.sqrt();
180    let mut self_e = 0.0;
181    for &q in charges {
182        self_e -= K_COULOMB * alpha / sqrt_pi * q * q;
183    }
184    self_e
185}
186
187#[cfg(test)]
188mod tests {
189    use super::*;
190
191    #[test]
192    fn test_cubic_box_volume() {
193        let box_vecs = BoxVectors::cubic(10.0);
194        assert!((box_vecs.volume() - 1000.0).abs() < 1e-6);
195    }
196
197    #[test]
198    fn test_direct_coulomb_flat_coords() {
199        // Two unit charges 1 Å apart
200        let coords = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]];
201        let charges = vec![1.0, -1.0];
202        let energy = direct_coulomb(&coords, &charges);
203        // Expected: K * 1.0 * (-1.0) / 1.0 = -332.0637 kcal/mol (approx)
204        let expected = -K_COULOMB;
205        assert!(
206            (energy - expected).abs() < 1.0,
207            "energy={}, expected={}",
208            energy,
209            expected
210        );
211    }
212
213    #[test]
214    fn test_spme_singular_box_returns_error() {
215        let degenerate_box = BoxVectors([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]);
216        let coords = [[0.0, 0.0, 0.0]];
217        let charges = [1.0];
218        let result = spme_energy(&coords, &charges, &degenerate_box, &PmeConfig::default());
219        assert_eq!(result, Err(EwaldError::SingularBoxMatrix));
220    }
221
222    #[test]
223    fn test_spme_valid_box() {
224        let box_vecs = BoxVectors::cubic(10.0);
225        let coords = [[5.0, 5.0, 5.0]];
226        let charges = [1.0];
227        let result = spme_energy(&coords, &charges, &box_vecs, &PmeConfig::default());
228        assert!(result.is_ok());
229    }
230}
chematic_ewald/lib.rs

chematic_ewald/
lib.rs
