Skip to main content

molframe_spatial/backends/
neighbor.rs

1//! Reusable Verlet neighbour candidates.
2
3use crate::brute::distance_squared;
4use crate::{CellList, NeighborListOptions, NeighborPair, PeriodicBox, SpatialError};
5use molframe_core::coords::CoordinateGeneration;
6
7const SAFE_DISPLACEMENT_FRACTION: f32 = 0.5;
8
9/// Candidate pairs built with a displacement skin for reuse across frames.
10#[derive(Debug)]
11pub struct NeighborList {
12    reference: Vec<[f32; 3]>,
13    candidates: Vec<(u32, u32)>,
14    cutoff: f32,
15    skin: f32,
16    periodic: Option<PeriodicBox>,
17    generation: CoordinateGeneration,
18}
19
20impl NeighborList {
21    /// Builds candidates out to `cutoff + skin`.
22    ///
23    /// # Errors
24    ///
25    /// Returns an error for invalid distances or atom indices.
26    pub fn build(
27        positions: &[[f32; 3]],
28        left: &[u32],
29        right: &[u32],
30        cutoff: f32,
31        skin: f32,
32        periodic: Option<&PeriodicBox>,
33        generation: CoordinateGeneration,
34    ) -> Result<Self, SpatialError> {
35        Self::build_with_options(
36            positions,
37            left,
38            right,
39            cutoff,
40            periodic,
41            generation,
42            NeighborListOptions::with_skin(skin),
43        )
44    }
45
46    /// Builds candidates under an explicit cell-grid memory policy.
47    ///
48    /// # Errors
49    ///
50    /// Returns an error for invalid options, distances or atom indices.
51    pub fn build_with_options(
52        positions: &[[f32; 3]],
53        left: &[u32],
54        right: &[u32],
55        cutoff: f32,
56        periodic: Option<&PeriodicBox>,
57        generation: CoordinateGeneration,
58        options: NeighborListOptions,
59    ) -> Result<Self, SpatialError> {
60        validate_distance(cutoff)?;
61        options.validate()?;
62
63        let outer = cutoff + options.skin;
64
65        if !outer.is_finite() {
66            return Err(SpatialError::InvalidCutoff);
67        }
68
69        let index = CellList::build_with_options(
70            positions,
71            right,
72            outer,
73            periodic.copied(),
74            options.cell_grid,
75        )?;
76        let pairs = if left == right {
77            crate::backends::cell::pairs_same_selection(&index, left, outer)?
78        } else {
79            index.pairs(left, outer)?
80        };
81
82        Ok(Self {
83            reference: positions.to_vec(),
84            candidates: pairs
85                .into_iter()
86                .map(|pair| (pair.first, pair.second))
87                .collect(),
88            cutoff,
89            skin: options.skin,
90            periodic: periodic.copied(),
91            generation,
92        })
93    }
94
95    /// The coordinate generation used to build this candidate set.
96    #[must_use]
97    pub const fn generation(&self) -> CoordinateGeneration {
98        self.generation
99    }
100
101    /// Whether a structure-local cache entry belongs to `generation`.
102    #[must_use]
103    pub const fn is_current(&self, generation: CoordinateGeneration) -> bool {
104        self.generation.get() == generation.get()
105    }
106
107    /// Whether no atom has moved more than half the skin.
108    ///
109    /// Reuse additionally requires the same coordinate count and the same
110    /// periodic box used to construct the candidate set.
111    #[must_use]
112    pub fn can_reuse(&self, positions: &[[f32; 3]], periodic: Option<&PeriodicBox>) -> bool {
113        if positions.len() != self.reference.len()
114            || !same_periodic_box(self.periodic.as_ref(), periodic)
115        {
116            return false;
117        }
118
119        let limit = self.skin * SAFE_DISPLACEMENT_FRACTION;
120        let limit_squared = limit * limit;
121
122        self.reference.iter().zip(positions).all(|(before, after)| {
123            displacement_within_limit(*before, *after, limit_squared, periodic)
124        })
125    }
126
127    /// Filters the candidates at the requested inner cutoff.
128    ///
129    /// # Errors
130    ///
131    /// Returns an error when positions moved far enough that an unseen pair
132    /// could have entered the cutoff, or when the cutoff is invalid.
133    pub fn pairs(
134        &self,
135        positions: &[[f32; 3]],
136        cutoff: f32,
137        periodic: Option<&PeriodicBox>,
138    ) -> Result<Vec<NeighborPair>, SpatialError> {
139        validate_query_cutoff(cutoff, self.cutoff)?;
140
141        if !self.can_reuse(positions, periodic) {
142            return Err(SpatialError::StaleNeighborList);
143        }
144
145        let cutoff_squared = cutoff * cutoff;
146        let mut found = Vec::new();
147
148        for &(left, right) in &self.candidates {
149            append_candidate_if_within(
150                positions,
151                left,
152                right,
153                cutoff_squared,
154                periodic,
155                &mut found,
156            )?;
157        }
158
159        // `candidates` originates from canonicalised CellList output.
160        // Filtering preserves both its ordering and uniqueness.
161        Ok(found)
162    }
163}
164
165/// Adds one cached candidate when its current distance remains within cutoff.
166///
167/// # Errors
168///
169/// Returns an out-of-bounds error if a cached candidate cannot be resolved in
170/// the supplied coordinate frame.
171fn append_candidate_if_within(
172    positions: &[[f32; 3]],
173    left: u32,
174    right: u32,
175    cutoff_squared: f32,
176    periodic: Option<&PeriodicBox>,
177    found: &mut Vec<NeighborPair>,
178) -> Result<(), SpatialError> {
179    let left_position = candidate_position(positions, left)?;
180    let right_position = candidate_position(positions, right)?;
181
182    let squared = distance_squared(left_position, right_position, periodic);
183
184    if squared <= cutoff_squared {
185        found.push(NeighborPair::new(left, right, squared));
186    }
187
188    Ok(())
189}
190
191/// Returns a cached candidate position or an out-of-bounds error.
192///
193/// Runtime and auxiliary space are `O(1)`.
194fn candidate_position(positions: &[[f32; 3]], atom: u32) -> Result<[f32; 3], SpatialError> {
195    let index = usize::try_from(atom).map_err(|_| SpatialError::NumericRangeExceeded)?;
196    positions
197        .get(index)
198        .copied()
199        .ok_or(SpatialError::AtomOutOfBounds(atom))
200}
201
202/// Tests whether one atom moved within the permitted squared displacement.
203///
204/// Runtime and auxiliary space are `O(1)`.
205#[inline]
206fn displacement_within_limit(
207    before: [f32; 3],
208    after: [f32; 3],
209    limit_squared: f32,
210    periodic: Option<&PeriodicBox>,
211) -> bool {
212    distance_squared(before, after, periodic) <= limit_squared
213}
214
215/// Compares the periodic cell used to build and query the neighbour list.
216///
217/// Value equality is used rather than pointer identity.
218#[inline]
219fn same_periodic_box(built: Option<&PeriodicBox>, current: Option<&PeriodicBox>) -> bool {
220    built == current
221}
222
223/// Validates a non-negative finite distance.
224///
225/// Runtime and auxiliary space are `O(1)`.
226fn validate_distance(value: f32) -> Result<(), SpatialError> {
227    if value.is_finite() && value >= 0.0 {
228        Ok(())
229    } else {
230        Err(SpatialError::InvalidCutoff)
231    }
232}
233
234/// Validates that a requested cutoff is covered by the candidate build cutoff.
235///
236/// Restricting queries to the original inner cutoff is required for the
237/// half-skin Verlet reuse guarantee.
238fn validate_query_cutoff(cutoff: f32, build_cutoff: f32) -> Result<(), SpatialError> {
239    if cutoff.is_finite() && cutoff >= 0.0 && cutoff <= build_cutoff {
240        Ok(())
241    } else {
242        Err(SpatialError::InvalidCutoff)
243    }
244}
245
246#[cfg(test)]
247#[path = "neighbor_tests.rs"]
248mod tests;