Skip to main content

molrs/builder/
occupancy.rs

1//! Lattice/grid occupancy tracking for self-avoidance.
2//!
3//! Overlap is decided **purely by cell occupancy** — never by pairwise distance
4//! or a neighbour list. Two occupancy models are supported, selected per
5//! [`GrowthStrategy`](super::GrowthStrategy):
6//!
7//! - [`OccupancyMode::SameCell`] — reject only if the candidate's own cell is
8//!   already taken. Used by lattice strategies whose step geometry already
9//!   guarantees a minimum separation (e.g. FCC: any two distinct sites are
10//!   `>= bond_length` apart), so the grid only has to forbid re-occupying a
11//!   site. The cell edge is chosen so each lattice site maps to a unique cell.
12//! - [`OccupancyMode::BlockClear`] — reject if the candidate's cell or any of
13//!   its 26 neighbours (excluding the bonding tip's cell) is occupied. With a
14//!   cell edge equal to the excluded radius this guarantees every pair of
15//!   non-bonded monomers is at least one cell — i.e. `excluded_radius` — apart.
16
17use std::collections::HashSet;
18
19use crate::spatial::simbox::SimBox;
20use crate::types::{F, Pbc3};
21
22/// How cell occupancy decides whether a candidate position overlaps.
23#[derive(Clone, Copy, Debug, PartialEq)]
24pub enum OccupancyMode {
25    /// Reject iff the candidate's own cell is occupied. `cell` is the grid edge.
26    SameCell {
27        /// Grid cell edge length.
28        cell: F,
29    },
30    /// Reject iff the candidate's cell or any 26-neighbour (minus the tip's
31    /// cell) is occupied. `cell` equals the guaranteed minimum separation.
32    BlockClear {
33        /// Grid cell edge length, equal to the enforced minimum separation.
34        cell: F,
35    },
36}
37
38impl OccupancyMode {
39    fn cell(self) -> F {
40        match self {
41            OccupancyMode::SameCell { cell } | OccupancyMode::BlockClear { cell } => cell,
42        }
43    }
44}
45
46/// A sparse periodic/reflective occupancy grid over the simulation box.
47pub(crate) struct OccupancyGrid {
48    mode: OccupancyMode,
49    ncells: [i64; 3],
50    pbc: Pbc3,
51    occupied: HashSet<[i64; 3]>,
52}
53
54impl OccupancyGrid {
55    /// Build an empty grid sized to `simbox` with the given boundary flags.
56    pub fn new(mode: OccupancyMode, simbox: &SimBox, pbc: Pbc3) -> Self {
57        let l = simbox.lengths();
58        let cell = mode.cell();
59        let ncells = [
60            (l[0] / cell).floor().max(1.0) as i64,
61            (l[1] / cell).floor().max(1.0) as i64,
62            (l[2] / cell).floor().max(1.0) as i64,
63        ];
64        Self {
65            mode,
66            ncells,
67            pbc,
68            occupied: HashSet::new(),
69        }
70    }
71
72    fn normalize_axis(&self, mut i: i64, ax: usize) -> i64 {
73        let n = self.ncells[ax];
74        if self.pbc[ax] {
75            i = i.rem_euclid(n);
76        } else {
77            i = i.clamp(0, n - 1);
78        }
79        i
80    }
81
82    fn cell_of(&self, p: [F; 3]) -> [i64; 3] {
83        let cell = self.mode.cell();
84        [
85            self.normalize_axis((p[0] / cell).floor() as i64, 0),
86            self.normalize_axis((p[1] / cell).floor() as i64, 1),
87            self.normalize_axis((p[2] / cell).floor() as i64, 2),
88        ]
89    }
90
91    fn normalize_cell(&self, c: [i64; 3]) -> [i64; 3] {
92        [
93            self.normalize_axis(c[0], 0),
94            self.normalize_axis(c[1], 1),
95            self.normalize_axis(c[2], 2),
96        ]
97    }
98
99    /// Is `p` free to occupy? `tip` (the bonding partner) is exempt from the
100    /// neighbour scan so a bonded step is never self-blocked.
101    pub fn is_free(&self, p: [F; 3], tip: Option<[F; 3]>) -> bool {
102        let c = self.cell_of(p);
103        match self.mode {
104            OccupancyMode::SameCell { .. } => !self.occupied.contains(&c),
105            OccupancyMode::BlockClear { .. } => {
106                let tip_cell = tip.map(|t| self.cell_of(t));
107                for dx in -1..=1 {
108                    for dy in -1..=1 {
109                        for dz in -1..=1 {
110                            let nc = self.normalize_cell([c[0] + dx, c[1] + dy, c[2] + dz]);
111                            if Some(nc) == tip_cell {
112                                continue;
113                            }
114                            if self.occupied.contains(&nc) {
115                                return false;
116                            }
117                        }
118                    }
119                }
120                true
121            }
122        }
123    }
124
125    /// Mark `p`'s cell occupied.
126    pub fn insert(&mut self, p: [F; 3]) {
127        let c = self.cell_of(p);
128        self.occupied.insert(c);
129    }
130
131    /// Free `p`'s cell (used when backtracking).
132    pub fn remove(&mut self, p: [F; 3]) {
133        let c = self.cell_of(p);
134        self.occupied.remove(&c);
135    }
136}