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}