chemrust_core/data/lattice/cell_param/
unit_cell.rs

1use std::{cmp::Ordering, fmt::Display};
2
3use crystallographic_group::database::CrystalSystem;
4use nalgebra::Matrix3;
5#[derive(Debug, Clone, Copy)]
6/// Lattice constants.
7/// Consider using newtype to wrap the angles f64 to ensure in radian form.
8pub struct CellConstants {
9    pub(crate) a: f64,
10    pub(crate) b: f64,
11    pub(crate) c: f64,
12    pub(crate) alpha: f64,
13    pub(crate) beta: f64,
14    pub(crate) gamma: f64,
15}
16
17#[derive(Debug, Clone, Copy)]
18pub struct LatticeVectors {
19    pub(crate) tensor: Matrix3<f64>,
20}
21
22pub trait UnitCellParameters {
23    fn cell_volume(&self) -> f64 {
24        self.lattice_bases().determinant()
25    }
26    fn lattice_bases(&self) -> Matrix3<f64>;
27    fn metric_tensor(&self) -> Matrix3<f64> {
28        let mat = self.lattice_bases();
29        let mat_transpose = mat.transpose();
30        mat_transpose * mat
31    }
32    fn length_a(&self) -> f64 {
33        CellConstants::from(self.lattice_bases()).a
34    }
35    fn length_b(&self) -> f64 {
36        CellConstants::from(self.lattice_bases()).b
37    }
38    fn length_c(&self) -> f64 {
39        CellConstants::from(self.lattice_bases()).c
40    }
41    /// Should return radians!
42    fn angle_alpha(&self) -> f64 {
43        CellConstants::from(self.lattice_bases()).alpha
44    }
45    /// Should return radians!
46    fn angle_beta(&self) -> f64 {
47        CellConstants::from(self.lattice_bases()).beta
48    }
49
50    /// Should return radians!
51    fn angle_gamma(&self) -> f64 {
52        CellConstants::from(self.lattice_bases()).gamma
53    }
54    fn get_crystal_system(&self) -> CrystalSystem {
55        let (length_a, length_b, length_c) = (self.length_a(), self.length_b(), self.length_c());
56        let (alpha, beta, gamma) = (self.angle_alpha(), self.angle_beta(), self.angle_gamma());
57        let axis_length_equal_count = [
58            compare_f64(length_a, length_b),
59            compare_f64(length_b, length_c),
60            compare_f64(length_a, length_c),
61        ]
62        .iter()
63        .filter(|ord| matches!(ord, Ordering::Equal))
64        .count();
65
66        let angle_eq_90_count = [
67            compare_f64(90.0, alpha),
68            compare_f64(90.0, beta),
69            compare_f64(90.0, gamma),
70        ]
71        .iter()
72        .filter(|ord| matches!(ord, Ordering::Equal))
73        .count();
74
75        let angle_eq_120_count = [
76            compare_f64(120.0, alpha),
77            compare_f64(120.0, beta),
78            compare_f64(120.0, gamma),
79        ]
80        .iter()
81        .filter(|ord| matches!(ord, Ordering::Equal))
82        .count();
83        match axis_length_equal_count {
84            3 => {
85                if angle_eq_90_count == 3 {
86                    CrystalSystem::Cubic
87                } else {
88                    CrystalSystem::Trigonal // Rhombohedral belongs to Trigonal
89                }
90            }
91            1 => {
92                if angle_eq_90_count == 3 {
93                    CrystalSystem::Tetragonal
94                } else if angle_eq_90_count == 2 && angle_eq_120_count == 1 {
95                    CrystalSystem::Hexagonal
96                } else {
97                    CrystalSystem::Triclinic
98                }
99            }
100            2 => {
101                // Floating point accuracy issue : a = b && b = c in tol but a != c
102                CrystalSystem::Triclinic
103            }
104            _ => match angle_eq_90_count {
105                3 => CrystalSystem::Orthorhombic,
106                2 => CrystalSystem::Monoclinic,
107                _ => CrystalSystem::Triclinic,
108            },
109        }
110    }
111}
112
113impl CellConstants {
114    /// angles should be specified in degrees
115    pub fn new(a: f64, b: f64, c: f64, alpha: f64, beta: f64, gamma: f64) -> Self {
116        Self {
117            a,
118            b,
119            c,
120            alpha,
121            beta,
122            gamma,
123        }
124    }
125}
126
127impl UnitCellParameters for CellConstants {
128    fn cell_volume(&self) -> f64 {
129        let CellConstants {
130            a,
131            b,
132            c,
133            alpha,
134            beta,
135            gamma,
136        } = self;
137        let cos_a = alpha.cos();
138        let cos_b = beta.cos();
139        let cos_y = gamma.cos();
140        a * b
141            * c
142            * (1.0 - cos_a * cos_a - cos_b * cos_b - cos_y * cos_y + 2.0 * cos_a * cos_b * cos_y)
143                .sqrt()
144    }
145
146    fn lattice_bases(&self) -> Matrix3<f64> {
147        let CellConstants {
148            a,
149            b,
150            c,
151            alpha,
152            beta,
153            gamma,
154        } = self;
155        let volume = self.cell_volume();
156        let cos_a = alpha.cos();
157        let cos_b = beta.cos();
158        let cos_y = gamma.cos();
159        let sin_y = gamma.sin();
160        //     [a         bcosy                     ccosB]
161        // A = [0         bsiny   c(cosa - cosbcosy)/siny]
162        //     [0             0                v/(absiny)]
163        // The columns are `a`, `b` and `c` vectors;
164        Matrix3::new(
165            *a,
166            b * cos_y,
167            c * cos_b,
168            0.0,
169            b * sin_y,
170            c * (cos_a - cos_b * cos_y) / sin_y,
171            0.0,
172            0.0,
173            volume / (a * b * sin_y),
174        )
175    }
176
177    fn length_a(&self) -> f64 {
178        self.a
179    }
180
181    fn length_b(&self) -> f64 {
182        self.b
183    }
184
185    fn length_c(&self) -> f64 {
186        self.c
187    }
188
189    fn angle_alpha(&self) -> f64 {
190        self.alpha
191    }
192
193    fn angle_beta(&self) -> f64 {
194        self.beta
195    }
196
197    fn angle_gamma(&self) -> f64 {
198        self.gamma
199    }
200}
201
202impl From<Matrix3<f64>> for CellConstants {
203    fn from(mat: Matrix3<f64>) -> Self {
204        let (v_a, v_b, v_c) = (mat.column(0), mat.column(1), mat.column(2));
205        let (a, b, c) = (v_a.norm(), v_b.norm(), v_c.norm());
206        let (alpha, beta, gamma) = (v_b.angle(&v_c), v_a.angle(&v_c), v_a.angle(&v_b));
207        Self {
208            a,
209            b,
210            c,
211            alpha,
212            beta,
213            gamma,
214        }
215    }
216}
217
218impl From<&Matrix3<f64>> for CellConstants {
219    fn from(value: &Matrix3<f64>) -> Self {
220        Self::from(*value)
221    }
222}
223
224impl Display for CellConstants {
225    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
226        write!(f, "a_length: {:>20.18}; b_length: {:>20.18}; c_length: {:>20.18}; alpha: {} beta: {} gamma: {}", self.a, self.b, self.c, self.alpha, self.beta, self.gamma)
227    }
228}
229
230impl LatticeVectors {
231    pub fn new(tensor: Matrix3<f64>) -> Self {
232        Self { tensor }
233    }
234
235    pub fn tensor(&self) -> Matrix3<f64> {
236        self.tensor
237    }
238}
239
240impl UnitCellParameters for LatticeVectors {
241    fn cell_volume(&self) -> f64 {
242        self.tensor().determinant()
243    }
244
245    fn lattice_bases(&self) -> Matrix3<f64> {
246        self.tensor()
247    }
248}
249
250impl From<CellConstants> for LatticeVectors {
251    fn from(constants: CellConstants) -> Self {
252        Self::new(constants.lattice_bases())
253    }
254}
255
256fn compare_f64(v1: f64, v2: f64) -> Ordering {
257    if (v1 - v2).abs() < 1e-6 {
258        Ordering::Equal
259    } else if v1 - v2 < -1e-6 {
260        Ordering::Less
261    } else {
262        Ordering::Greater
263    }
264}
265
266#[cfg(test)]
267mod test {
268    use nalgebra::{Matrix3, Point3, Rotation3, Vector3};
269
270    use crate::{
271        data::lattice::cell_param::unit_cell::UnitCellParameters,
272        systems::crystal_model::rotated_lattice_tensor,
273    };
274
275    use super::CellConstants;
276
277    #[test]
278    fn cell_repr() {
279        let lattice_cart = Matrix3::new(
280            // a
281            18.931530020488704480,
282            -0.000000000000003553,
283            0.000000000000000000,
284            // b
285            -9.465765010246645517,
286            16.395185930251127360,
287            0.000000000000000000,
288            // c
289            0.000000000000000000,
290            0.000000000000000000,
291            9.999213039981000861,
292        );
293        let cell = CellConstants::from(lattice_cart);
294        println!("{}", cell);
295        println!("{:#>20.18}", cell.lattice_bases());
296        let p = Point3::new(0.07560343470042601, 0.0756034355668187, 0.5000000004346841);
297        let o: Point3<f64> = Point3::origin();
298        let b = cell.lattice_bases().column(1).xyz();
299        let j = Vector3::y_axis();
300        let rot = Rotation3::rotation_between(&b, &j).unwrap();
301        let mat = rotated_lattice_tensor(&cell, rot);
302        let cart_p = mat * p;
303        println!("{:#.5}", mat);
304        println!("{:?}", j);
305        println!("{:#}", cart_p);
306        let frac_p =
307            cell.lattice_bases().try_inverse().unwrap() * rot.matrix() * cell.lattice_bases() * p;
308        println!("{:#}", frac_p);
309        println!("{:#}", cell.lattice_bases() * frac_p);
310        let po_cart = (cart_p - o).norm_squared();
311        let metric_tensor = cell.metric_tensor();
312        let po = frac_p - o;
313        let po_norm_squared = po.transpose() * metric_tensor * po;
314        println!("{:#}", metric_tensor);
315        println!(
316            "V^2: {}, det(G) : {}",
317            cell.cell_volume().powi(2),
318            metric_tensor.determinant()
319        );
320        println!(
321            "cart_length: {}, frac_length by metric tensor: {}",
322            po_cart, po_norm_squared.x
323        );
324    }
325}
chemrust_core/data/lattice/cell_param/unit_cell.rs

chemrust_core/data/lattice/cell_param/
unit_cell.rs
