molframe_spatial/geometry/
periodic.rs

1//! Orthorhombic and triclinic minimum-image displacement.
2
3use crate::SpatialError;
4use crate::numeric::{f64_f32, i64_f64, rounded_i64};
5use molframe_core::structure::UnitCell;
6
7const ORTHOGONAL_ANGLE_TOLERANCE: f64 = 1.0e-10;
8
9const IMAGE_OFFSETS: [[f64; 3]; 27] = [
10    [-1.0, -1.0, -1.0],
11    [-1.0, -1.0, 0.0],
12    [-1.0, -1.0, 1.0],
13    [-1.0, 0.0, -1.0],
14    [-1.0, 0.0, 0.0],
15    [-1.0, 0.0, 1.0],
16    [-1.0, 1.0, -1.0],
17    [-1.0, 1.0, 0.0],
18    [-1.0, 1.0, 1.0],
19    [0.0, -1.0, -1.0],
20    [0.0, -1.0, 0.0],
21    [0.0, -1.0, 1.0],
22    [0.0, 0.0, -1.0],
23    [0.0, 0.0, 0.0],
24    [0.0, 0.0, 1.0],
25    [0.0, 1.0, -1.0],
26    [0.0, 1.0, 0.0],
27    [0.0, 1.0, 1.0],
28    [1.0, -1.0, -1.0],
29    [1.0, -1.0, 0.0],
30    [1.0, -1.0, 1.0],
31    [1.0, 0.0, -1.0],
32    [1.0, 0.0, 0.0],
33    [1.0, 0.0, 1.0],
34    [1.0, 1.0, -1.0],
35    [1.0, 1.0, 0.0],
36    [1.0, 1.0, 1.0],
37];
38
39/// An invertible unit-cell basis and its inverse.
40#[derive(Clone, Copy, PartialEq, Debug)]
41pub struct PeriodicBox {
42    basis: [[f64; 3]; 3],
43    inverse: [[f64; 3]; 3],
44    orthogonal: bool,
45}
46
47/// A shortest displacement on a periodic flat torus and the chosen image.
48#[derive(Clone, Copy, PartialEq, Debug)]
49pub struct PeriodicImage {
50    /// Cartesian displacement from the left point to the chosen right image.
51    pub displacement: [f32; 3],
52    /// Lattice translation subtracted from the un-wrapped right point.
53    pub lattice_shift: [i64; 3],
54}
55
56impl PeriodicBox {
57    /// Builds a periodic box from crystallographic lengths and angles.
58    ///
59    /// # Errors
60    ///
61    /// Returns [`SpatialError::InvalidCell`] when any parameter is non-finite,
62    /// an edge is not positive, or the vectors are coplanar.
63    pub fn from_cell(cell: UnitCell) -> Result<Self, SpatialError> {
64        validate_cell(&cell)?;
65
66        let orthogonal = cell
67            .angles
68            .iter()
69            .all(|angle| (*angle - 90.0).abs() < ORTHOGONAL_ANGLE_TOLERANCE);
70
71        let basis = basis_from_cell(&cell)?;
72        let Some(inverse) = inverse3(basis) else {
73            return Err(SpatialError::InvalidCell);
74        };
75
76        Ok(Self {
77            basis,
78            inverse,
79            orthogonal,
80        })
81    }
82
83    /// The shortest periodic displacement from `left` to `right`.
84    ///
85    /// The nearest rounded fractional image is sufficient for an orthogonal
86    /// cell. In a skewed cell a neighbouring lattice translation can be
87    /// shorter, so all twenty-seven candidates around it are checked in a
88    /// fixed order.
89    #[must_use]
90    pub fn displacement(&self, left: [f32; 3], right: [f32; 3]) -> [f32; 3] {
91        self.minimum_image(left, right).displacement
92    }
93
94    /// The shortest displacement and exact lattice image selected for it.
95    ///
96    /// Exact distance ties retain the deposited image. Triclinic candidates
97    /// are visited in a fixed order, so the selected image is deterministic.
98    #[must_use]
99    pub fn minimum_image(&self, left: [f32; 3], right: [f32; 3]) -> PeriodicImage {
100        let (best, lattice_shift) = self.image_f64(left.map(f64::from), right.map(f64::from));
101
102        PeriodicImage {
103            displacement: best.map(f64_f32),
104            lattice_shift,
105        }
106    }
107
108    /// Minimum-image displacement without rounding generated sample positions
109    /// to stored-coordinate precision. Image selection matches `minimum_image`.
110    #[must_use]
111    pub fn displacement_f64(&self, left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
112        self.image_f64(left, right).0
113    }
114
115    fn image_f64(&self, left: [f64; 3], right: [f64; 3]) -> ([f64; 3], [i64; 3]) {
116        let cartesian = cartesian_delta(left, right);
117        let fractional = multiply(self.inverse, cartesian);
118        if self.orthogonal {
119            orthogonal_image(self.basis, cartesian, fractional)
120        } else {
121            triclinic_image(self.basis, cartesian, fractional)
122        }
123    }
124
125    /// Squared minimum-image separation.
126    #[must_use]
127    pub fn distance_squared(&self, left: [f32; 3], right: [f32; 3]) -> f32 {
128        let displacement = self.displacement(left, right);
129
130        displacement[0] * displacement[0]
131            + displacement[1] * displacement[1]
132            + displacement[2] * displacement[2]
133    }
134
135    /// Converts a Cartesian position to fractional cell coordinates.
136    #[must_use]
137    pub fn fractional(&self, position: [f32; 3]) -> [f64; 3] {
138        multiply(self.inverse, position.map(f64::from))
139    }
140
141    /// Converts fractional cell coordinates to Cartesian ångström.
142    #[must_use]
143    pub fn cartesian(&self, fractional: [f64; 3]) -> [f32; 3] {
144        multiply(self.basis, fractional).map(f64_f32)
145    }
146
147    /// Wraps one position into the primary fractional interval `[0, 1)`.
148    #[must_use]
149    pub fn wrap(&self, position: [f32; 3]) -> [f32; 3] {
150        let fractional = self.fractional(position).map(|value| value.rem_euclid(1.0));
151
152        self.cartesian(fractional)
153    }
154
155    /// Interpolates along the shortest periodic geodesic and wraps the result.
156    ///
157    /// Returns `None` for a non-finite interpolation parameter. Values outside
158    /// `[0, 1]` extrapolate along the same selected image.
159    #[must_use]
160    pub fn interpolate(&self, left: [f32; 3], right: [f32; 3], amount: f64) -> Option<[f32; 3]> {
161        if !amount.is_finite() {
162            return None;
163        }
164        let displacement = self.minimum_image(left, right).displacement;
165        let point = std::array::from_fn(|axis| {
166            f64_f32(f64::from(left[axis]) + amount * f64::from(displacement[axis]))
167        });
168        Some(self.wrap(point))
169    }
170
171    /// Maximum fractional-coordinate changes produced by a Cartesian radius.
172    ///
173    /// The bounds are the Euclidean norms of the reciprocal-basis rows times
174    /// the radius. They conservatively bound cell and tree image searches.
175    #[must_use]
176    pub(crate) fn fractional_cutoff_bounds(&self, cutoff: f32) -> [f64; 3] {
177        self.inverse
178            .map(|row| squared(row).sqrt() * f64::from(cutoff))
179    }
180
181    /// Translates a Cartesian point by an exact lattice shift.
182    pub(crate) fn translated(
183        &self,
184        point: [f32; 3],
185        shift: [i64; 3],
186    ) -> Result<[f32; 3], SpatialError> {
187        let convert = |value| i64_f64(value).ok_or(SpatialError::NumericRangeExceeded);
188        let translation = multiply(
189            self.basis,
190            [convert(shift[0])?, convert(shift[1])?, convert(shift[2])?],
191        );
192        Ok(std::array::from_fn(|axis| {
193            f64_f32(f64::from(point[axis]) + translation[axis])
194        }))
195    }
196}
197
198/// Validates finite cell lengths/angles and positive edge lengths.
199///
200/// Runtime and auxiliary space are `O(1)`.
201fn validate_cell(cell: &UnitCell) -> Result<(), SpatialError> {
202    let finite = cell
203        .lengths
204        .iter()
205        .chain(&cell.angles)
206        .all(|value| value.is_finite());
207
208    let positive = cell.lengths.iter().all(|length| *length > 0.0);
209
210    if finite && positive {
211        Ok(())
212    } else {
213        Err(SpatialError::InvalidCell)
214    }
215}
216
217/// Builds the crystallographic Cartesian basis matrix.
218///
219/// # Errors
220///
221/// Returns [`SpatialError::InvalidCell`] for singular or numerically invalid
222/// angle combinations.
223fn basis_from_cell(cell: &UnitCell) -> Result<[[f64; 3]; 3], SpatialError> {
224    let [a, b, c] = cell.lengths;
225    let [alpha, beta, gamma] = cell.angles.map(f64::to_radians);
226
227    let sin_gamma = gamma.sin();
228
229    if !sin_gamma.is_finite() || sin_gamma.abs() <= f64::EPSILON {
230        return Err(SpatialError::InvalidCell);
231    }
232
233    let cos_gamma = gamma.cos();
234    let cos_beta = beta.cos();
235    let cos_alpha = alpha.cos();
236
237    let cx = c * cos_beta;
238    let cy = c * (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
239    let cz_squared = c.mul_add(c, -cx.mul_add(cx, cy * cy));
240
241    if !cz_squared.is_finite() || cz_squared <= f64::EPSILON {
242        return Err(SpatialError::InvalidCell);
243    }
244
245    Ok([
246        [a, b * cos_gamma, cx],
247        [0.0, b * sin_gamma, cy],
248        [0.0, 0.0, cz_squared.sqrt()],
249    ])
250}
251
252/// Computes Cartesian `right - left` in `f64` arithmetic.
253///
254/// Runtime and auxiliary space are `O(1)`.
255#[inline]
256fn cartesian_delta(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
257    [right[0] - left[0], right[1] - left[1], right[2] - left[2]]
258}
259
260/// Computes the minimum image for an orthogonal cell.
261///
262/// The original Cartesian displacement is retained on exact distance ties,
263/// preserving the displacement sign convention at half-cell boundaries.
264fn orthogonal_image(
265    basis: [[f64; 3]; 3],
266    cartesian: [f64; 3],
267    fractional: [f64; 3],
268) -> ([f64; 3], [i64; 3]) {
269    let rounded = fractional.map(f64::round);
270    let wrapped = [
271        fractional[0] - rounded[0],
272        fractional[1] - rounded[1],
273        fractional[2] - rounded[2],
274    ];
275
276    let candidate = multiply(basis, wrapped);
277
278    if squared(candidate) < squared(cartesian) {
279        (candidate, integer_shift(rounded))
280    } else {
281        (cartesian, [0; 3])
282    }
283}
284
285/// Computes the minimum image for a skewed triclinic cell.
286///
287/// All 27 translations around the rounded fractional image are evaluated in a
288/// deterministic order. Runtime and auxiliary space are `O(1)`.
289fn triclinic_image(
290    basis: [[f64; 3]; 3],
291    cartesian: [f64; 3],
292    fractional: [f64; 3],
293) -> ([f64; 3], [i64; 3]) {
294    let centre = fractional.map(f64::round);
295    let base = [
296        fractional[0] - centre[0],
297        fractional[1] - centre[1],
298        fractional[2] - centre[2],
299    ];
300
301    let mut best = cartesian;
302    let mut best_squared = squared(best);
303    let mut best_shift = [0; 3];
304
305    for offset in IMAGE_OFFSETS {
306        let image = [
307            base[0] - offset[0],
308            base[1] - offset[1],
309            base[2] - offset[2],
310        ];
311
312        let candidate = multiply(basis, image);
313        let candidate_squared = squared(candidate);
314
315        if candidate_squared < best_squared {
316            best = candidate;
317            best_squared = candidate_squared;
318            best_shift = integer_shift([
319                centre[0] + offset[0],
320                centre[1] + offset[1],
321                centre[2] + offset[2],
322            ]);
323        }
324    }
325
326    (best, best_shift)
327}
328
329fn integer_shift(shift: [f64; 3]) -> [i64; 3] {
330    shift.map(rounded_i64)
331}
332
333/// Multiplies a three-by-three matrix by a three-component vector.
334///
335/// Runtime and auxiliary space are `O(1)`.
336#[inline]
337fn multiply(matrix: [[f64; 3]; 3], vector: [f64; 3]) -> [f64; 3] {
338    [
339        matrix[0][0] * vector[0] + matrix[0][1] * vector[1] + matrix[0][2] * vector[2],
340        matrix[1][0] * vector[0] + matrix[1][1] * vector[1] + matrix[1][2] * vector[2],
341        matrix[2][0] * vector[0] + matrix[2][1] * vector[1] + matrix[2][2] * vector[2],
342    ]
343}
344
345/// Returns a three-dimensional squared vector length.
346///
347/// Runtime and auxiliary space are `O(1)`.
348#[inline]
349fn squared(vector: [f64; 3]) -> f64 {
350    vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]
351}
352
353/// Computes all cofactors of a three-by-three matrix.
354///
355/// Runtime and auxiliary space are `O(1)`.
356fn cofactors3(matrix: [[f64; 3]; 3]) -> [[f64; 3]; 3] {
357    [
358        [
359            matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1],
360            matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2],
361            matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0],
362        ],
363        [
364            matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2],
365            matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0],
366            matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1],
367        ],
368        [
369            matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1],
370            matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2],
371            matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0],
372        ],
373    ]
374}
375
376/// Computes a three-by-three inverse using the adjugate formula.
377///
378/// Returns `None` for non-finite or effectively singular matrices. Runtime and
379/// auxiliary space are `O(1)`.
380fn inverse3(matrix: [[f64; 3]; 3]) -> Option<[[f64; 3]; 3]> {
381    let cofactors = cofactors3(matrix);
382
383    let determinant = matrix[0][0] * cofactors[0][0]
384        + matrix[0][1] * cofactors[0][1]
385        + matrix[0][2] * cofactors[0][2];
386
387    if !determinant.is_finite() || determinant.abs() <= f64::EPSILON {
388        return None;
389    }
390
391    let reciprocal = determinant.recip();
392
393    Some([
394        [
395            cofactors[0][0] * reciprocal,
396            cofactors[1][0] * reciprocal,
397            cofactors[2][0] * reciprocal,
398        ],
399        [
400            cofactors[0][1] * reciprocal,
401            cofactors[1][1] * reciprocal,
402            cofactors[2][1] * reciprocal,
403        ],
404        [
405            cofactors[0][2] * reciprocal,
406            cofactors[1][2] * reciprocal,
407            cofactors[2][2] * reciprocal,
408        ],
409    ])
410}
411
412#[cfg(test)]
413#[path = "periodic_tests.rs"]
414mod tests;
molframe_spatial/geometry/periodic.rs

molframe_spatial/geometry/
periodic.rs
