phasesmith_crystallography/
cell.rs

1//! Unit-cell and reciprocal-metric calculations.
2
3use std::error::Error;
4use std::fmt::{Display, Formatter};
5
6/// Number and order of direct unit-cell parameters: `a`, `b`, `c`, `alpha`,
7/// `beta`, `gamma`.
8pub const CELL_PARAMETER_COUNT: usize = 6;
9
10/// Three-by-three row-major matrix.
11pub type Matrix3 = [[f64; 3]; 3];
12
13/// Direct unit cell in ångströms and degrees.
14#[derive(Clone, Copy, Debug, PartialEq)]
15pub struct UnitCell {
16    /// Direct `a` length in ångströms.
17    pub a_angstrom: f64,
18    /// Direct `b` length in ångströms.
19    pub b_angstrom: f64,
20    /// Direct `c` length in ångströms.
21    pub c_angstrom: f64,
22    /// Angle between `b` and `c`, in degrees.
23    pub alpha_deg: f64,
24    /// Angle between `a` and `c`, in degrees.
25    pub beta_deg: f64,
26    /// Angle between `a` and `b`, in degrees.
27    pub gamma_deg: f64,
28}
29
30/// Validated direct/reciprocal geometry derived from a [`UnitCell`].
31#[derive(Clone, Copy, Debug, PartialEq)]
32pub struct CellGeometry {
33    /// Cartesian direct basis vectors as matrix columns, in ångströms.
34    pub direct_basis: Matrix3,
35    /// Cartesian reciprocal basis vectors as matrix columns, in inverse
36    /// ångströms without a `2 pi` factor.
37    pub reciprocal_basis: Matrix3,
38    /// Direct metric tensor in square ångströms.
39    pub direct_metric: Matrix3,
40    /// Reciprocal metric tensor in inverse square ångströms.
41    pub reciprocal_metric: Matrix3,
42    /// Unit-cell volume in cubic ångströms.
43    pub volume_angstrom3: f64,
44    direct_metric_derivatives: [Matrix3; CELL_PARAMETER_COUNT],
45}
46
47/// Unit-cell validation or reciprocal-space evaluation failure.
48#[derive(Clone, Copy, Debug, PartialEq, Eq)]
49pub enum CellError {
50    /// A length or angle is not finite.
51    NonFiniteParameter,
52    /// A direct length is not positive.
53    NonPositiveLength,
54    /// An angle does not lie strictly within 0 and 180 degrees.
55    InvalidAngle,
56    /// The six parameters do not define a positive-volume cell.
57    DegenerateCell,
58    /// The zero Miller index does not define a finite d-spacing.
59    ZeroReflection,
60}
61
62impl Display for CellError {
63    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
64        formatter.write_str(match self {
65            Self::NonFiniteParameter => "unit-cell parameters must be finite",
66            Self::NonPositiveLength => "unit-cell lengths must be positive",
67            Self::InvalidAngle => "unit-cell angles must lie strictly within (0, 180) degrees",
68            Self::DegenerateCell => "unit-cell parameters must define a positive finite volume",
69            Self::ZeroReflection => "hkl = (0, 0, 0) has no finite d-spacing",
70        })
71    }
72}
73
74impl Error for CellError {}
75
76impl UnitCell {
77    /// Validate the six parameters and derive direct/reciprocal geometry.
78    ///
79    /// # Errors
80    ///
81    /// Returns [`CellError`] when values are non-finite, lengths/angles are
82    /// outside their domains, or the resulting metric is degenerate.
83    pub fn geometry(self) -> Result<CellGeometry, CellError> {
84        let parameters = [
85            self.a_angstrom,
86            self.b_angstrom,
87            self.c_angstrom,
88            self.alpha_deg,
89            self.beta_deg,
90            self.gamma_deg,
91        ];
92        if parameters.iter().any(|value| !value.is_finite()) {
93            return Err(CellError::NonFiniteParameter);
94        }
95        if parameters[..3].iter().any(|value| *value <= 0.0) {
96            return Err(CellError::NonPositiveLength);
97        }
98        if parameters[3..]
99            .iter()
100            .any(|value| !(*value > 0.0 && *value < 180.0))
101        {
102            return Err(CellError::InvalidAngle);
103        }
104
105        let radians_per_degree = std::f64::consts::PI / 180.0;
106        let alpha = self.alpha_deg * radians_per_degree;
107        let beta = self.beta_deg * radians_per_degree;
108        let gamma = self.gamma_deg * radians_per_degree;
109        let (cos_alpha, cos_beta, cos_gamma) = (alpha.cos(), beta.cos(), gamma.cos());
110        let sin_gamma = gamma.sin();
111        let direct_metric = [
112            [
113                self.a_angstrom * self.a_angstrom,
114                self.a_angstrom * self.b_angstrom * cos_gamma,
115                self.a_angstrom * self.c_angstrom * cos_beta,
116            ],
117            [
118                self.a_angstrom * self.b_angstrom * cos_gamma,
119                self.b_angstrom * self.b_angstrom,
120                self.b_angstrom * self.c_angstrom * cos_alpha,
121            ],
122            [
123                self.a_angstrom * self.c_angstrom * cos_beta,
124                self.b_angstrom * self.c_angstrom * cos_alpha,
125                self.c_angstrom * self.c_angstrom,
126            ],
127        ];
128        let determinant = determinant(direct_metric);
129        if !determinant.is_finite() || determinant <= 0.0 || sin_gamma <= 0.0 {
130            return Err(CellError::DegenerateCell);
131        }
132        let volume = determinant.sqrt();
133        let c_x = self.c_angstrom * cos_beta;
134        let c_y = self.c_angstrom * (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
135        let c_z = volume / (self.a_angstrom * self.b_angstrom * sin_gamma);
136        if !c_y.is_finite() || !c_z.is_finite() || c_z <= 0.0 {
137            return Err(CellError::DegenerateCell);
138        }
139        let direct_basis = [
140            [self.a_angstrom, self.b_angstrom * cos_gamma, c_x],
141            [0.0, self.b_angstrom * sin_gamma, c_y],
142            [0.0, 0.0, c_z],
143        ];
144        let reciprocal_metric = inverse(direct_metric).ok_or(CellError::DegenerateCell)?;
145        let reciprocal_basis = transpose(inverse(direct_basis).ok_or(CellError::DegenerateCell)?);
146        let direct_metric_derivatives = direct_metric_derivatives(self, alpha, beta, gamma);
147        Ok(CellGeometry {
148            direct_basis,
149            reciprocal_basis,
150            direct_metric,
151            reciprocal_metric,
152            volume_angstrom3: volume,
153            direct_metric_derivatives,
154        })
155    }
156}
157
158impl CellGeometry {
159    /// Return reciprocal-axis lengths and their direct-cell derivatives.
160    ///
161    /// The lengths are `a*`, `b*`, and `c*` without a `2 pi` factor. The
162    /// derivative rows follow reciprocal-axis order and columns follow direct
163    /// cell parameter order.
164    #[must_use]
165    pub fn reciprocal_axis_lengths_and_derivatives(
166        &self,
167    ) -> ([f64; 3], [[f64; CELL_PARAMETER_COUNT]; 3]) {
168        let lengths = [
169            self.reciprocal_metric[0][0].sqrt(),
170            self.reciprocal_metric[1][1].sqrt(),
171            self.reciprocal_metric[2][2].sqrt(),
172        ];
173        let mut derivatives = [[0.0; CELL_PARAMETER_COUNT]; 3];
174        for axis in 0..3 {
175            let reciprocal_column = [
176                self.reciprocal_metric[0][axis],
177                self.reciprocal_metric[1][axis],
178                self.reciprocal_metric[2][axis],
179            ];
180            for (parameter, derivative) in derivatives[axis].iter_mut().enumerate() {
181                let product =
182                    matrix_vector(self.direct_metric_derivatives[parameter], reciprocal_column);
183                let d_reciprocal_diagonal = -dot(reciprocal_column, product);
184                *derivative = 0.5 * d_reciprocal_diagonal / lengths[axis];
185            }
186        }
187        (lengths, derivatives)
188    }
189
190    /// Return `|g|² = hᵀ G* h` without constructing cell derivatives.
191    #[must_use]
192    pub fn q_squared(&self, hkl: [i32; 3]) -> f64 {
193        let h = [f64::from(hkl[0]), f64::from(hkl[1]), f64::from(hkl[2])];
194        dot(h, matrix_vector(self.reciprocal_metric, h))
195    }
196
197    /// Return `|g|^2 = h^T G* h` and its derivatives in direct-cell parameter
198    /// order. Unlike d-spacing evaluation, the zero reflection is permitted.
199    #[must_use]
200    pub fn q_squared_and_derivatives(&self, hkl: [i32; 3]) -> (f64, [f64; CELL_PARAMETER_COUNT]) {
201        let h = [f64::from(hkl[0]), f64::from(hkl[1]), f64::from(hkl[2])];
202        let reciprocal_h = matrix_vector(self.reciprocal_metric, h);
203        let q_squared = dot(h, reciprocal_h);
204        let mut derivatives = [0.0; CELL_PARAMETER_COUNT];
205        for (parameter, derivative) in derivatives.iter_mut().enumerate() {
206            let d_metric_h = matrix_vector(self.direct_metric_derivatives[parameter], reciprocal_h);
207            *derivative = -dot(reciprocal_h, d_metric_h);
208        }
209        (q_squared, derivatives)
210    }
211
212    /// Return d-spacing and derivatives in direct-cell parameter order.
213    ///
214    /// # Errors
215    ///
216    /// Returns [`CellError::ZeroReflection`] for `hkl = (0, 0, 0)`.
217    pub fn d_spacing_and_derivatives(
218        &self,
219        hkl: [i32; 3],
220    ) -> Result<(f64, [f64; CELL_PARAMETER_COUNT]), CellError> {
221        let (q_squared, q_derivatives) = self.q_squared_and_derivatives(hkl);
222        if q_squared <= 0.0 || !q_squared.is_finite() {
223            return Err(CellError::ZeroReflection);
224        }
225        let d_spacing = q_squared.sqrt().recip();
226        let factor = -0.5 * d_spacing.powi(3);
227        Ok((d_spacing, q_derivatives.map(|value| factor * value)))
228    }
229
230    /// Return volume derivatives in direct-cell parameter order.
231    #[must_use]
232    pub fn volume_derivatives(&self) -> [f64; CELL_PARAMETER_COUNT] {
233        let mut derivatives = [0.0; CELL_PARAMETER_COUNT];
234        for (parameter, derivative) in derivatives.iter_mut().enumerate() {
235            let product = multiply(
236                self.reciprocal_metric,
237                self.direct_metric_derivatives[parameter],
238            );
239            *derivative =
240                0.5 * self.volume_angstrom3 * (product[0][0] + product[1][1] + product[2][2]);
241        }
242        derivatives
243    }
244}
245
246fn direct_metric_derivatives(
247    cell: UnitCell,
248    alpha: f64,
249    beta: f64,
250    gamma: f64,
251) -> [Matrix3; CELL_PARAMETER_COUNT] {
252    let zero = [[0.0; 3]; 3];
253    let mut derivatives = [zero; CELL_PARAMETER_COUNT];
254    let (ca, cb, cg) = (alpha.cos(), beta.cos(), gamma.cos());
255    derivatives[0] = [
256        [
257            2.0 * cell.a_angstrom,
258            cell.b_angstrom * cg,
259            cell.c_angstrom * cb,
260        ],
261        [cell.b_angstrom * cg, 0.0, 0.0],
262        [cell.c_angstrom * cb, 0.0, 0.0],
263    ];
264    derivatives[1] = [
265        [0.0, cell.a_angstrom * cg, 0.0],
266        [
267            cell.a_angstrom * cg,
268            2.0 * cell.b_angstrom,
269            cell.c_angstrom * ca,
270        ],
271        [0.0, cell.c_angstrom * ca, 0.0],
272    ];
273    derivatives[2] = [
274        [0.0, 0.0, cell.a_angstrom * cb],
275        [0.0, 0.0, cell.b_angstrom * ca],
276        [
277            cell.a_angstrom * cb,
278            cell.b_angstrom * ca,
279            2.0 * cell.c_angstrom,
280        ],
281    ];
282    let radians_per_degree = std::f64::consts::PI / 180.0;
283    let d_alpha = -cell.b_angstrom * cell.c_angstrom * alpha.sin() * radians_per_degree;
284    derivatives[3][1][2] = d_alpha;
285    derivatives[3][2][1] = d_alpha;
286    let d_beta = -cell.a_angstrom * cell.c_angstrom * beta.sin() * radians_per_degree;
287    derivatives[4][0][2] = d_beta;
288    derivatives[4][2][0] = d_beta;
289    let d_gamma = -cell.a_angstrom * cell.b_angstrom * gamma.sin() * radians_per_degree;
290    derivatives[5][0][1] = d_gamma;
291    derivatives[5][1][0] = d_gamma;
292    derivatives
293}
294
295fn determinant(matrix: Matrix3) -> f64 {
296    matrix[0][0] * (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])
297        - matrix[0][1] * (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])
298        + matrix[0][2] * (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])
299}
300
301fn inverse(matrix: Matrix3) -> Option<Matrix3> {
302    let det = determinant(matrix);
303    if !det.is_finite() || det == 0.0 {
304        return None;
305    }
306    let inverse_det = det.recip();
307    Some([
308        [
309            (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) * inverse_det,
310            (matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2]) * inverse_det,
311            (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1]) * inverse_det,
312        ],
313        [
314            (matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2]) * inverse_det,
315            (matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0]) * inverse_det,
316            (matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2]) * inverse_det,
317        ],
318        [
319            (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * inverse_det,
320            (matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1]) * inverse_det,
321            (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]) * inverse_det,
322        ],
323    ])
324}
325
326fn transpose(matrix: Matrix3) -> Matrix3 {
327    [
328        [matrix[0][0], matrix[1][0], matrix[2][0]],
329        [matrix[0][1], matrix[1][1], matrix[2][1]],
330        [matrix[0][2], matrix[1][2], matrix[2][2]],
331    ]
332}
333
334fn multiply(left: Matrix3, right: Matrix3) -> Matrix3 {
335    let mut result = [[0.0; 3]; 3];
336    for (row, result_row) in result.iter_mut().enumerate() {
337        for (column, value) in result_row.iter_mut().enumerate() {
338            *value = (0..3)
339                .map(|inner| left[row][inner] * right[inner][column])
340                .sum();
341        }
342    }
343    result
344}
345
346fn matrix_vector(matrix: Matrix3, vector: [f64; 3]) -> [f64; 3] {
347    matrix.map(|row| dot(row, vector))
348}
349
350fn dot(left: [f64; 3], right: [f64; 3]) -> f64 {
351    left[0] * right[0] + left[1] * right[1] + left[2] * right[2]
352}
353
354#[cfg(test)]
355mod tests {
356    use super::*;
357
358    fn triclinic() -> UnitCell {
359        UnitCell {
360            a_angstrom: 4.3,
361            b_angstrom: 5.1,
362            c_angstrom: 6.2,
363            alpha_deg: 78.0,
364            beta_deg: 83.0,
365            gamma_deg: 71.0,
366        }
367    }
368
369    fn changed(mut cell: UnitCell, parameter: usize, delta: f64) -> UnitCell {
370        match parameter {
371            0 => cell.a_angstrom += delta,
372            1 => cell.b_angstrom += delta,
373            2 => cell.c_angstrom += delta,
374            3 => cell.alpha_deg += delta,
375            4 => cell.beta_deg += delta,
376            5 => cell.gamma_deg += delta,
377            _ => unreachable!(),
378        }
379        cell
380    }
381
382    #[test]
383    fn cubic_geometry_is_exact() {
384        let geometry = UnitCell {
385            a_angstrom: 4.0,
386            b_angstrom: 4.0,
387            c_angstrom: 4.0,
388            alpha_deg: 90.0,
389            beta_deg: 90.0,
390            gamma_deg: 90.0,
391        }
392        .geometry()
393        .expect("cubic cell");
394        assert!((geometry.volume_angstrom3 - 64.0).abs() < 1.0e-12);
395        assert!((geometry.reciprocal_metric[0][0] - 1.0 / 16.0).abs() < 1.0e-15);
396        let (d, _) = geometry
397            .d_spacing_and_derivatives([1, 1, 0])
398            .expect("nonzero reflection");
399        assert!((d - 4.0 / 2.0_f64.sqrt()).abs() < 1.0e-14);
400    }
401
402    #[test]
403    fn cell_derivatives_match_centered_differences() {
404        let cell = triclinic();
405        let geometry = cell.geometry().expect("triclinic cell");
406        let (d, derivatives) = geometry
407            .d_spacing_and_derivatives([2, -1, 3])
408            .expect("reflection");
409        let volume_derivatives = geometry.volume_derivatives();
410        let (reciprocal_lengths, reciprocal_derivatives) =
411            geometry.reciprocal_axis_lengths_and_derivatives();
412        for parameter in 0..CELL_PARAMETER_COUNT {
413            let step = if parameter < 3 { 1.0e-6 } else { 1.0e-5 };
414            let plus = changed(cell, parameter, step)
415                .geometry()
416                .expect("plus cell");
417            let minus = changed(cell, parameter, -step)
418                .geometry()
419                .expect("minus cell");
420            let plus_d = plus
421                .d_spacing_and_derivatives([2, -1, 3])
422                .expect("plus reflection")
423                .0;
424            let minus_d = minus
425                .d_spacing_and_derivatives([2, -1, 3])
426                .expect("minus reflection")
427                .0;
428            let finite_d = (plus_d - minus_d) / (2.0 * step);
429            let finite_volume = (plus.volume_angstrom3 - minus.volume_angstrom3) / (2.0 * step);
430            assert!((derivatives[parameter] - finite_d).abs() < 2.0e-8 * d.max(1.0));
431            assert!(
432                (volume_derivatives[parameter] - finite_volume).abs()
433                    < 2.0e-8 * geometry.volume_angstrom3
434            );
435            let plus_lengths = plus.reciprocal_axis_lengths_and_derivatives().0;
436            let minus_lengths = minus.reciprocal_axis_lengths_and_derivatives().0;
437            for axis in 0..3 {
438                let finite_reciprocal = (plus_lengths[axis] - minus_lengths[axis]) / (2.0 * step);
439                assert!(
440                    (reciprocal_derivatives[axis][parameter] - finite_reciprocal).abs()
441                        < 2.0e-8 * reciprocal_lengths[axis].max(1.0)
442                );
443            }
444        }
445    }
446
447    #[test]
448    fn invalid_cells_and_zero_reflection_are_rejected() {
449        let mut cell = triclinic();
450        cell.a_angstrom = 0.0;
451        assert_eq!(cell.geometry(), Err(CellError::NonPositiveLength));
452        let geometry = triclinic().geometry().expect("cell");
453        assert_eq!(
454            geometry.d_spacing_and_derivatives([0, 0, 0]),
455            Err(CellError::ZeroReflection)
456        );
457    }
458}
phasesmith_crystallography/cell.rs

phasesmith_crystallography/
cell.rs
