qsym2/rotsym/
mod.rs

1//! Rotational symmetry based on moments of inertia.
2
3use std::fmt;
4
5use approx;
6use log;
7use nalgebra as na;
8use serde::{Deserialize, Serialize};
9
10#[cfg(test)]
11#[path = "rotsym_tests.rs"]
12mod rotsym_tests;
13
14/// Enumerated type to classify the types of rotational symmetry of a molecular system based on its
15/// principal moments of inertia.
16#[derive(Clone, Debug, Serialize, Deserialize)]
17pub enum RotationalSymmetry {
18    /// All three principal moments of inertia are identical.
19    Spherical,
20    /// The unique principal moment of inertia is the largest, the other two
21    /// are equal and sum to the unique one.
22    OblatePlanar,
23    /// The unique principal moment of inertia is the largest, the other two
24    /// are equal but do not sum to the unique one.
25    OblateNonPlanar,
26    /// The unique principal moment of inertia is zero, the other two are equal.
27    ProlateLinear,
28    /// The unique principal moment of inertia is the smallest but non-zero,
29    /// the other two are equal.
30    ProlateNonLinear,
31    /// The largest principal moment of inertia is the sum of the other two,
32    /// but they are all distinct.
33    AsymmetricPlanar,
34    /// All principal moments of inertia are distinct and do not have any
35    /// special relations between them.
36    AsymmetricNonPlanar,
37}
38
39impl fmt::Display for RotationalSymmetry {
40    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
41        match self {
42            RotationalSymmetry::Spherical => write!(f, "Spherical"),
43            RotationalSymmetry::OblatePlanar => write!(f, "Oblate, planar"),
44            RotationalSymmetry::OblateNonPlanar => write!(f, "Oblate, non-planar"),
45            RotationalSymmetry::ProlateLinear => write!(f, "Prolate, linear"),
46            RotationalSymmetry::ProlateNonLinear => write!(f, "Prolate, non-linear"),
47            RotationalSymmetry::AsymmetricPlanar => write!(f, "Asymmetric, planar"),
48            RotationalSymmetry::AsymmetricNonPlanar => write!(f, "Asymmetric, non-planar"),
49        }
50    }
51}
52
53/// Determines the rotational symmetry given an inertia tensor.
54///
55/// # Arguments
56///
57/// * `inertia_tensor` - An inertia tensor which is a $`3 \times 3`$ matrix.
58/// * `thresh` - A threshold for comparing moments of inertia.
59///
60/// # Returns
61///
62/// The rotational symmetry as one of the [`RotationalSymmetry`] variants.
63///
64/// # Panics
65///
66/// Panics when the moments of inertia cannot be compared.
67#[must_use]
68pub fn calc_rotational_symmetry(
69    inertia_tensor: &na::Matrix3<f64>,
70    thresh: f64,
71) -> RotationalSymmetry {
72    let moi_mat = inertia_tensor.symmetric_eigenvalues();
73    let mut moi: Vec<&f64> = moi_mat.iter().collect();
74    moi.sort_by(|a, b| {
75        a.partial_cmp(b)
76            .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
77    });
78    log::debug!("Moments of inertia:");
79    for component in &moi {
80        log::debug!("  {component:+.14}");
81    }
82    if approx::relative_eq!(*moi[0], *moi[1], epsilon = thresh, max_relative = thresh) {
83        if approx::relative_eq!(*moi[1], *moi[2], epsilon = thresh, max_relative = thresh) {
84            return RotationalSymmetry::Spherical;
85        }
86        if approx::relative_eq!(
87            *moi[2],
88            *moi[0] + *moi[1],
89            epsilon = thresh,
90            max_relative = thresh
91        ) {
92            return RotationalSymmetry::OblatePlanar;
93        }
94        return RotationalSymmetry::OblateNonPlanar;
95    }
96    if approx::relative_eq!(*moi[1], *moi[2], epsilon = thresh, max_relative = thresh) {
97        if moi[0].abs() < thresh {
98            return RotationalSymmetry::ProlateLinear;
99        }
100        return RotationalSymmetry::ProlateNonLinear;
101    }
102    if approx::relative_eq!(
103        *moi[2],
104        *moi[0] + *moi[1],
105        epsilon = thresh,
106        max_relative = thresh
107    ) {
108        return RotationalSymmetry::AsymmetricPlanar;
109    }
110    RotationalSymmetry::AsymmetricNonPlanar
111}