qsym2/symmetry/symmetry_element/
mod.rs

1//! Geometrical symmetry elements.
2
3use std::convert::TryInto;
4use std::fmt;
5use std::hash::{Hash, Hasher};
6
7use approx;
8use derive_builder::Builder;
9use fraction;
10use log;
11use nalgebra::Vector3;
12use num::integer::gcd;
13use num_traits::{ToPrimitive, Zero};
14use serde::{Deserialize, Serialize};
15
16use crate::auxiliary::geometry;
17use crate::auxiliary::misc::{self, HashableFloat};
18use crate::symmetry::symmetry_element_order::ElementOrder;
19
20type F = fraction::GenericFraction<u32>;
21
22pub mod symmetry_operation;
23pub use symmetry_operation::*;
24
25#[cfg(test)]
26mod symmetry_element_tests;
27
28// ====================================
29// Enum definitions and implementations
30// ====================================
31
32/// Enumerated type to classify the type of the antiunitary term that contributes to a symmetry
33/// element.
34#[derive(Copy, Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
35pub enum AntiunitaryKind {
36    /// Variant for the antiunitary term being a complex-conjugation operation.
37    ComplexConjugation,
38
39    /// Variant for the antiunitary term being a time-reversal operation.
40    TimeReversal,
41}
42
43/// Enumerated type to classify the types of symmetry element.
44#[derive(Copy, Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
45pub enum SymmetryElementKind {
46    /// Proper symmetry element which consists of just a proper rotation axis.
47    ///
48    /// The associated option indicates the type of the associated antiunitary operation, if any.
49    Proper(Option<AntiunitaryKind>),
50
51    /// Improper symmetry element in the mirror-plane convention, which consists of a proper
52    /// rotation axis and an orthogonal mirror plane.
53    ///
54    /// The associated option indicates the type of the associated antiunitary operation, if any.
55    ImproperMirrorPlane(Option<AntiunitaryKind>),
56
57    /// Improper symmetry element in the inversion-centre convention, which consists of a proper
58    /// rotation axis and an inversion centre.
59    ///
60    /// The associated option indicates the type of the associated antiunitary operation, if any.
61    ImproperInversionCentre(Option<AntiunitaryKind>),
62}
63
64impl SymmetryElementKind {
65    /// Indicates if the antiunitary part of this element contains a time-reversal operation.
66    #[must_use]
67    pub fn contains_time_reversal(&self) -> bool {
68        match self {
69            Self::Proper(tr)
70            | Self::ImproperMirrorPlane(tr)
71            | Self::ImproperInversionCentre(tr) => *tr == Some(AntiunitaryKind::TimeReversal),
72        }
73    }
74
75    /// Indicates if an antiunitary operation is associated with this element.
76    #[must_use]
77    pub fn contains_antiunitary(&self) -> bool {
78        match self {
79            Self::Proper(au)
80            | Self::ImproperMirrorPlane(au)
81            | Self::ImproperInversionCentre(au) => au.is_some(),
82        }
83    }
84
85    /// Converts the current kind to the desired time-reversal form.
86    ///
87    /// # Arguments
88    ///
89    /// * `tr` - A boolean indicating whether time reversal is included (`true`) or removed
90    /// (`false`).
91    ///
92    /// # Returns
93    ///
94    /// A copy of the current kind with or without the time-reversal antiunitary operation.
95    #[must_use]
96    pub fn to_tr(&self, tr: bool) -> Self {
97        if tr {
98            self.to_antiunitary(Some(AntiunitaryKind::TimeReversal))
99        } else {
100            self.to_antiunitary(None)
101        }
102    }
103
104    /// Converts the associated antiunitary operation to the desired kind.
105    ///
106    /// # Arguments
107    ///
108    /// * `au` - An option containing the desired antiunitary kind.
109    ///
110    /// # Returns
111    ///
112    /// A new symmetry element kind with the desired antiunitary kind.
113    pub fn to_antiunitary(&self, au: Option<AntiunitaryKind>) -> Self {
114        match self {
115            Self::Proper(_) => Self::Proper(au),
116            Self::ImproperMirrorPlane(_) => Self::ImproperMirrorPlane(au),
117            Self::ImproperInversionCentre(_) => Self::ImproperInversionCentre(au),
118        }
119    }
120}
121
122/// Enumerated type to indicate the rotation group to which the unitary proper rotation part of a
123/// symmetry element belongs.
124#[derive(Clone, Hash, PartialEq, Eq, Debug, Serialize, Deserialize)]
125pub enum RotationGroup {
126    /// Variant indicating that the proper part of the symmetry element generates rotations in
127    /// $`\mathsf{SO}(3)`$.
128    SO3,
129
130    /// Variant indicating that the proper part of the symmetry element generates rotations in
131    /// $`\mathsf{SU}(2)`$. The associated boolean indicates whether the proper part of the element
132    /// itself is connected to the identity via a homotopy path of class 0 (`true`) or class 1
133    /// (`false`).
134    SU2(bool),
135}
136
137impl RotationGroup {
138    /// Indicates if the rotation is in $`\mathsf{SU}(2)`$.
139    pub fn is_su2(&self) -> bool {
140        matches!(self, RotationGroup::SU2(_))
141    }
142
143    /// Indicates if the rotation is in $`\mathsf{SU}(2)`$ and connected to the
144    /// identity via a homotopy path of class 1.
145    pub fn is_su2_class_1(&self) -> bool {
146        matches!(self, RotationGroup::SU2(false))
147    }
148}
149
150// ======================================
151// Struct definitions and implementations
152// ======================================
153
154pub struct SymmetryElementKindConversionError(String);
155
156impl fmt::Debug for SymmetryElementKindConversionError {
157    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
158        f.debug_struct("SymmetryElementKindConversionError")
159            .field("Message", &self.0)
160            .finish()
161    }
162}
163
164impl fmt::Display for SymmetryElementKindConversionError {
165    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
166        write!(
167            f,
168            "SymmetryElementKindConversionError with message: {}",
169            &self.0
170        )
171    }
172}
173
174impl std::error::Error for SymmetryElementKindConversionError {}
175
176impl TryInto<geometry::ImproperRotationKind> for SymmetryElementKind {
177    type Error = SymmetryElementKindConversionError;
178
179    fn try_into(self) -> Result<geometry::ImproperRotationKind, Self::Error> {
180        match self {
181            Self::Proper(_) => Err(SymmetryElementKindConversionError(
182                "Unable to convert a proper element to an `ImproperRotationKind`.".to_string(),
183            )),
184            Self::ImproperMirrorPlane(_) => Ok(geometry::ImproperRotationKind::MirrorPlane),
185            Self::ImproperInversionCentre(_) => Ok(geometry::ImproperRotationKind::InversionCentre),
186        }
187    }
188}
189
190impl fmt::Display for SymmetryElementKind {
191    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
192        match self {
193            Self::Proper(au) => match au {
194                None => write!(f, "Proper"),
195                Some(AntiunitaryKind::TimeReversal) => write!(f, "Time-reversed proper"),
196                Some(AntiunitaryKind::ComplexConjugation) => write!(f, "Complex-conjugated proper"),
197            },
198            Self::ImproperMirrorPlane(au) => match au {
199                None => write!(f, "Improper (mirror-plane convention)"),
200                Some(AntiunitaryKind::TimeReversal) => {
201                    write!(f, "Time-reversed improper (mirror-plane convention)")
202                }
203                Some(AntiunitaryKind::ComplexConjugation) => {
204                    write!(f, "Complex-conjugated improper (mirror-plane convention)")
205                }
206            },
207            Self::ImproperInversionCentre(au) => match au {
208                None => write!(f, "Improper (inversion-centre convention)"),
209                Some(AntiunitaryKind::TimeReversal) => {
210                    write!(f, "Time-reversed improper (inversion-centre convention)")
211                }
212                Some(AntiunitaryKind::ComplexConjugation) => {
213                    write!(
214                        f,
215                        "Complex-conjugated improper (inversion-centre convention)"
216                    )
217                }
218            },
219        }
220    }
221}
222
223/// Structure for storing and managing symmetry elements.
224///
225/// Each symmetry element is a geometrical object in $`\mathbb{R}^3`$ that encodes the following
226/// pieces of information:
227///
228/// * the axis of rotation $`\hat{\mathbf{n}}`$,
229/// * the angle of rotation $`\phi`$,
230/// * the associated improper operation, if any, and
231/// * the associated antiunitary operation, if any.
232///
233/// These pieces of information can be stored in the following representation of a symmetry element
234/// $`\hat{g}`$:
235///
236/// ```math
237/// \hat{g} = \hat{\alpha} \hat{\gamma} \hat{C}_n^k,
238/// ```
239///
240/// where
241/// * $`n \in \mathbb{N}_{+}`$, $`k \in \mathbb{Z}/n\mathbb{Z}`$ such that
242/// $`\lfloor -n/2 \rfloor < k \le \lfloor n/2 \rfloor`$,
243/// * $`\hat{\gamma}`$ is either the identity $`\hat{e}`$, the inversion operation $`\hat{i}`$, or
244/// a reflection operation $`\hat{\sigma}`$ perpendicular to the axis of rotation,
245/// * $`\hat{\alpha}`$ is either the identity $`\hat{e}`$, the complex conjugation $`\hat{K}`$, or
246/// the time reversal $`\hat{\theta}`$.
247///
248/// We shall refer to $`\hat{C}_n^k`$ as the *unitary proper rotation part*, $`\hat{\gamma}`$ the
249/// *improper rotation part*, and $`\hat{\alpha}`$ the *antiunitary part* of the symmetry element.
250///
251/// With this definition, the above pieces of information required to specify a geometrical symmetry
252/// element are given as follows:
253///
254/// * the axis of rotation $`\hat{\mathbf{n}}`$ is given by the axis of $`\hat{C}_n^k`$,
255/// * the angle of rotation $`\phi = 2\pi k/n \in (-\pi, \pi]`$,
256/// * the improper contribution $`\hat{\gamma}`$,
257/// * the antiunitary contribution $`\hat{\alpha}`$.
258///
259/// This definition also allows the unitary part of $`\hat{g}`$ (*i.e.* the composition
260/// $`\hat{\gamma} \hat{C}_n^k`$) to be interpreted as an element of either $`\mathsf{O}(3)`$ or
261/// $`\mathsf{SU}'(2)`$, which means that the unitary part of $`\hat{g}`$ is also a symmetry
262/// operation in the corresponding group, and a rather special one that can be used to generate
263/// other symmetry operations of the group. $`\hat{g}`$ thus serves as a bridge between molecular
264/// symmetry and abstract group theory.
265///
266/// There is one small caveat: for infinite-order elements, $`n`$ and $`k`$ can no longer be used
267/// to give the angle of rotation. There must thus be a mechanism to allow for infinite-order
268/// elements to be interpreted as an arbitrary finite-order one. An explicit specification of the
269/// angle of rotation $`\phi`$ seems to be the best way to do this. In other words, the angle of
270/// rotation of each element is specified by either a tuple of integers $`(k, n)`$ or a
271/// floating-point number $`\phi`$.
272#[derive(Builder, Clone, Serialize, Deserialize)]
273pub struct SymmetryElement {
274    /// The rotational order $`n`$ of the unitary proper rotation part of the symmetry element. This
275    /// can be finite or infinite, and will determine whether the proper power is `None` or
276    /// contains an integer value.
277    #[builder(setter(name = "proper_order"))]
278    raw_proper_order: ElementOrder,
279
280    /// The power $`k \in \mathbb{Z}/n\mathbb{Z}`$ of the unitary proper rotation part of the
281    /// symmetry element such that $`\lfloor -n/2 \rfloor < k <= \lfloor n/2 \rfloor`$. This is only
282    /// defined if [`Self::raw_proper_order`] is finite.
283    ///
284    /// The unitary proper rotation does not include any additional unitary rotations introduced by
285    /// the antiunitary part of the symmetry element (*e.g.* time reversal).
286    #[builder(setter(custom, name = "proper_power"), default = "None")]
287    raw_proper_power: Option<i32>,
288
289    /// The normalised axis of the unitary proper rotation part of the symmetry element whose
290    /// direction is as specified at construction.
291    #[builder(setter(custom))]
292    raw_axis: Vector3<f64>,
293
294    /// The spatial and antiunitary kind of the symmetry element.
295    #[builder(default = "SymmetryElementKind::Proper(None)")]
296    kind: SymmetryElementKind,
297
298    /// The rotation group in which the unitary proper rotation part of the symmetry element shall
299    /// be interpreted.
300    rotation_group: RotationGroup,
301
302    /// A boolean indicating whether the symmetry element is a generator of the group to which it
303    /// belongs.
304    #[builder(default = "false")]
305    generator: bool,
306
307    /// A threshold for approximate equality comparisons.
308    #[builder(setter(custom))]
309    threshold: f64,
310
311    /// An additional superscript for distinguishing symmetry elements.
312    #[builder(default = "String::new()")]
313    pub(crate) additional_superscript: String,
314
315    /// An additional subscript for distinguishing symmetry elements.
316    #[builder(default = "String::new()")]
317    pub(crate) additional_subscript: String,
318
319    /// The fraction $`k/n \in (-1/2, 1/2]`$ of the unitary proper rotation, represented exactly
320    /// for hashing and comparison purposes.
321    ///
322    /// The unitary proper rotation does not include any additional unitary rotations introduced by
323    /// the antiunitary part of the symmetry element (*e.g.* time reversal).
324    ///
325    /// This is not defined for infinite-order elements and cannot be set arbitrarily.
326    #[builder(setter(skip), default = "self.calc_proper_fraction()")]
327    proper_fraction: Option<F>,
328
329    /// The normalised proper angle $`\phi \in (-\pi, \pi]`$ corresponding to the unitary proper
330    /// rotation.
331    ///
332    /// The unitary proper rotation does not include any additional unitary rotations introduced by
333    /// the antiunitary part of the symmetry element (*e.g.* time reversal).
334    ///
335    /// This can be set arbitrarily only for infinite-order elements.
336    #[builder(setter(custom), default = "self.calc_proper_angle()")]
337    proper_angle: Option<f64>,
338}
339
340impl SymmetryElementBuilder {
341    /// Sets the power of the unitary proper rotation part of the element.
342    ///
343    /// # Arguments
344    ///
345    /// * `prop_pow` - A proper power to be set. This will be folded into the interval
346    /// $`(\lfloor -n/2 \rfloor, \lfloor n/2 \rfloor]`$.
347    pub fn proper_power(&mut self, prop_pow: i32) -> &mut Self {
348        let raw_proper_order = self
349            .raw_proper_order
350            .as_ref()
351            .expect("Proper order has not been set.");
352        self.raw_proper_power = match raw_proper_order {
353            ElementOrder::Int(io) => {
354                let io_i32 =
355                    i32::try_from(*io).expect("Unable to convert the integer order to `i32`.");
356                let residual = prop_pow.rem_euclid(io_i32);
357                if residual > io_i32.div_euclid(2) {
358                    Some(Some(residual - io_i32))
359                } else {
360                    Some(Some(residual))
361                }
362            }
363            ElementOrder::Inf => None,
364        };
365        self
366    }
367
368    /// Sets the rotation angle of the unitary proper rotation part of the infinite-order element.
369    ///
370    /// # Arguments
371    ///
372    /// * `ang` - A proper rotation angle to be set. This will be folded into the interval
373    /// $`(-\pi, \pi]`$.
374    ///
375    /// # Panics
376    ///
377    /// Panics when `self` is of finite order.
378    pub fn proper_angle(&mut self, ang: f64) -> &mut Self {
379        let proper_order = self
380            .raw_proper_order
381            .as_ref()
382            .expect("Proper order has not been set.");
383        self.proper_angle = match proper_order {
384            ElementOrder::Int(_) => panic!(
385                "Arbitrary proper rotation angles can only be set for infinite-order elements."
386            ),
387            ElementOrder::Inf => {
388                let (normalised_rotation_angle, _) = geometry::normalise_rotation_angle(
389                    ang,
390                    self.threshold.expect("Threshold value has not been set."),
391                );
392                Some(Some(normalised_rotation_angle))
393            }
394        };
395        self
396    }
397
398    /// Sets the raw axis of the element.
399    ///
400    /// # Arguments
401    ///
402    /// * `axs` - The raw axis which will be normalised.
403    pub fn raw_axis(&mut self, axs: Vector3<f64>) -> &mut Self {
404        let thresh = self.threshold.expect("Threshold value has not been set.");
405        if approx::abs_diff_eq!(axs.norm(), 1.0, epsilon = thresh) {
406            self.raw_axis = Some(axs);
407        } else {
408            log::warn!("Axis not normalised. Normalising...");
409            self.raw_axis = Some(axs.normalize());
410        }
411        self
412    }
413
414    /// Sets the comparison threshold of the element.
415    ///
416    /// # Arguments
417    ///
418    /// * `thresh` - The comparison threshold..
419    pub fn threshold(&mut self, thresh: f64) -> &mut Self {
420        if thresh >= 0.0 {
421            self.threshold = Some(thresh);
422        } else {
423            log::error!(
424                "Threshold value `{}` is invalid. Threshold must be non-negative.",
425                thresh
426            );
427            self.threshold = None;
428        }
429        self
430    }
431
432    /// Calculates the fraction $`k/n \in (-1/2, 1/2]`$ of the unitary proper rotation, represented
433    /// exactly for hashing and comparison purposes.
434    ///
435    /// The unitary proper rotation does not include any additional unitary rotations introduced by
436    /// the antiunitary part of the symmetry element (*e.g.* time reversal).
437    ///
438    /// This is not defined for infinite-order elements and cannot be set arbitrarily.
439    fn calc_proper_fraction(&self) -> Option<F> {
440        let raw_proper_order = self
441            .raw_proper_order
442            .as_ref()
443            .expect("Proper order has not been set.");
444        match raw_proper_order {
445            ElementOrder::Int(io) => {
446                // The generating element has a proper fraction, pp/n.
447                //
448                // If pp/n > 1/2, we seek a positive integer x such that
449                //  -1/2 < pp/n - x <= 1/2.
450                // It turns out that x ∈ [pp/n - 1/2, pp/n + 1/2).
451                //
452                // If pp/n <= -1/2, we seek a positive integer x such that
453                //  -1/2 < pp/n + x <= 1/2.
454                // It turns out that x ∈ (-pp/n - 1/2, -pp/n + 1/2].
455                //
456                // x is then used to bring pp/n back into the (-1/2, 1/2] interval.
457                //
458                // See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover
459                // Publications, Inc., New York, 2005) for further information.
460                let pp = self
461                    .raw_proper_power
462                    .expect("Proper power has not been set.")
463                    .expect("No proper powers found.");
464                let total_proper_fraction = if pp >= 0 {
465                    F::new(pp.unsigned_abs(), *io)
466                } else {
467                    F::new_neg(pp.unsigned_abs(), *io)
468                };
469                Some(total_proper_fraction)
470            }
471            ElementOrder::Inf => None,
472        }
473    }
474
475    /// Calculates the normalised proper angle $`\phi \in (-\pi, \pi]`$ corresponding to the unitary
476    /// proper rotation.
477    ///
478    /// The unitary proper rotation does not include any additional unitary rotations introduced by
479    /// the antiunitary part of the symmetry element (*e.g.* time reversal).
480    ///
481    /// This can be set arbitrarily only for infinite-order elements.
482    fn calc_proper_angle(&self) -> Option<f64> {
483        let raw_proper_order = self
484            .raw_proper_order
485            .as_ref()
486            .expect("Proper order has not been set.");
487        match raw_proper_order {
488            ElementOrder::Int(io) => {
489                let pp = self
490                    .raw_proper_power
491                    .expect("Proper power has not been set.")
492                    .expect("No proper powers found.");
493                let total_proper_fraction = if pp >= 0 {
494                    F::new(pp.unsigned_abs(), *io)
495                } else {
496                    F::new_neg(pp.unsigned_abs(), *io)
497                };
498                Some(
499                    total_proper_fraction
500                        .to_f64()
501                        .expect("Unable to convert the proper fraction to `f64`.")
502                        * 2.0
503                        * std::f64::consts::PI,
504                )
505            }
506            ElementOrder::Inf => self.proper_angle.unwrap_or(None),
507        }
508    }
509}
510
511impl SymmetryElement {
512    /// Returns a builder to construct a new symmetry element.
513    ///
514    /// # Returns
515    ///
516    /// A builder to construct a new symmetry element.
517    #[must_use]
518    pub fn builder() -> SymmetryElementBuilder {
519        SymmetryElementBuilder::default()
520    }
521
522    /// Returns the raw order of the unitary proper rotation. This might not be equal to the value of
523    /// $`n`$ if the fraction $`k/n`$ has been reduced.
524    pub fn raw_proper_order(&self) -> &ElementOrder {
525        &self.raw_proper_order
526    }
527
528    /// Returns the raw power of the unitary proper rotation. This might not be equal to the value of
529    /// $`k`$ if the fraction $`k/n`$ has been reduced.
530    pub fn raw_proper_power(&self) -> Option<&i32> {
531        self.raw_proper_power.as_ref()
532    }
533
534    /// Returns the raw axis of the unitary proper rotation.
535    pub fn raw_axis(&self) -> &Vector3<f64> {
536        &self.raw_axis
537    }
538
539    /// Returns the axis of the unitary proper rotation in the standard positive hemisphere.
540    pub fn standard_positive_axis(&self) -> Vector3<f64> {
541        geometry::get_standard_positive_pole(&self.raw_axis, self.threshold)
542    }
543
544    /// Returns the axis of the unitary proper rotation multiplied by the sign of the rotation angle.
545    /// If the proper rotation is a binary rotation, then the positive axis is always returned.
546    pub fn signed_axis(&self) -> Vector3<f64> {
547        let au = self.contains_antiunitary();
548        if self.is_o3_binary_rotation_axis(au) || self.is_o3_mirror_plane(au) {
549            self.standard_positive_axis()
550        } else {
551            self.proper_fraction
552                .map(|frac| {
553                    frac.signum()
554                        .to_f64()
555                        .expect("Unable to obtain the sign of the proper fraction.")
556                })
557                .or_else(|| self.proper_angle.map(|proper_angle| proper_angle.signum())).map(|signum| signum * self.raw_axis)
558                .unwrap_or_else(|| {
559                    log::warn!("No rotation signs could be obtained. The positive axis will be used for the signed axis.");
560                    self.standard_positive_axis()
561                })
562        }
563    }
564
565    /// Returns the pole of the unitary proper rotation part of the element while leaving any improper
566    /// and antiunitary contributions intact.
567    ///
568    /// If the unitary proper rotation part if a binary rotation, the pole is always in the standard
569    /// positive hemisphere.
570    ///
571    /// # Returns
572    ///
573    /// The position vector of the proper rotation pole.
574    pub fn proper_rotation_pole(&self) -> Vector3<f64> {
575        match *self.raw_proper_order() {
576            ElementOrder::Int(_) => {
577                let frac_1_2 = F::new(1u32, 2u32);
578                let proper_fraction = self.proper_fraction.expect("No proper fractions found.");
579                if proper_fraction == frac_1_2 {
580                    // Binary proper rotations
581                    geometry::get_standard_positive_pole(&self.raw_axis, self.threshold)
582                } else if proper_fraction > F::zero() {
583                    // Positive proper rotation angles
584                    self.raw_axis
585                } else if proper_fraction < F::zero() {
586                    // Negative proper rotation angles
587                    -self.raw_axis
588                } else {
589                    // Identity or inversion
590                    assert!(proper_fraction.is_zero());
591                    Vector3::zeros()
592                }
593            }
594            ElementOrder::Inf => {
595                if approx::abs_diff_eq!(
596                    self.proper_angle.expect("No proper angles found."),
597                    std::f64::consts::PI,
598                    epsilon = self.threshold
599                ) {
600                    // Binary proper rotations
601                    geometry::get_standard_positive_pole(&self.raw_axis, self.threshold)
602                } else if approx::abs_diff_ne!(
603                    self.proper_angle.expect("No proper angles found."),
604                    0.0,
605                    epsilon = self.threshold
606                ) {
607                    self.proper_angle.expect("No proper angles found.").signum() * self.raw_axis
608                } else {
609                    approx::assert_abs_diff_eq!(
610                        self.proper_angle.expect("No proper angles found."),
611                        0.0,
612                        epsilon = self.threshold
613                    );
614                    Vector3::zeros()
615                }
616            }
617        }
618    }
619
620    /// Returns the unitary proper fraction for this element, if any.
621    ///
622    /// The element lacks a unitary proper fraction if it is infinite-order.
623    pub fn proper_fraction(&self) -> Option<&F> {
624        self.proper_fraction.as_ref()
625    }
626
627    /// Returns the unitary proper angle for this element, if any.
628    ///
629    /// The element lacks a unitary proper angle if it is infinite-order and the rotation angle has
630    /// not been set.
631    pub fn proper_angle(&self) -> Option<f64> {
632        self.proper_angle
633    }
634
635    /// Returns the spatial and antiunitary kind of this element.
636    pub fn kind(&self) -> &SymmetryElementKind {
637        &self.kind
638    }
639
640    /// Returns the rotation group and possibly the identity-connected homotopy class in which the
641    /// proper rotation part of this element is to be interpreted.
642    pub fn rotation_group(&self) -> &RotationGroup {
643        &self.rotation_group
644    }
645
646    /// Returns a boolean indicating if the element is a generator of a group.
647    pub fn is_generator(&self) -> bool {
648        self.generator
649    }
650
651    /// Returns the threshold for approximate comparisons.
652    pub fn threshold(&self) -> f64 {
653        self.threshold
654    }
655
656    /// Checks if the symmetry element contains a time-reversal operator as the antiunitary part.
657    ///
658    /// # Returns
659    ///
660    /// A boolean indicating if the symmetry element contains a time-reversal operator as the
661    /// antiunitary part.
662    #[must_use]
663    pub fn contains_time_reversal(&self) -> bool {
664        self.kind.contains_time_reversal()
665    }
666
667    /// Checks if the symmetry element contains an antiunitary part.
668    ///
669    /// # Returns
670    ///
671    /// Returns `None` if the symmetry element has no antiunitary parts, or `Some` wrapping around
672    /// the antiunitary kind if the symmetry element contains an antiunitary part.
673    pub fn contains_antiunitary(&self) -> Option<AntiunitaryKind> {
674        match self.kind {
675            SymmetryElementKind::Proper(au)
676            | SymmetryElementKind::ImproperMirrorPlane(au)
677            | SymmetryElementKind::ImproperInversionCentre(au) => au,
678        }
679    }
680
681    /// Checks if the unitary proper rotation part of the element is in $`\mathsf{SU}(2)`$.
682    #[must_use]
683    pub fn is_su2(&self) -> bool {
684        self.rotation_group.is_su2()
685    }
686
687    /// Checks if the unitary proper rotation part of the element is in $`\mathsf{SU}(2)`$ and
688    /// connected to the identity via a homotopy path of class 1.
689    ///
690    /// See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover Publications, Inc., New
691    /// York, 2005) for further information.
692    #[must_use]
693    pub fn is_su2_class_1(&self) -> bool {
694        self.rotation_group.is_su2_class_1()
695    }
696
697    /// Checks if the spatial part of the symmetry element is proper and has the specified
698    /// antiunitary attribute.
699    ///
700    /// More specifically, this checks if the combination $`\hat{\gamma} \hat{C}_n^k`$ is proper,
701    /// and if $`\hat{\alpha}`$ is as specified by `au`.
702    ///
703    /// # Arguments
704    ///
705    /// * `au` - An `Option` for the desired antiunitary kind.
706    ///
707    /// # Returns
708    ///
709    /// A boolean indicating if the symmetry element is proper and has the specified antiunitary
710    /// attribute.
711    #[must_use]
712    pub fn is_o3_proper(&self, au: Option<AntiunitaryKind>) -> bool {
713        self.kind == SymmetryElementKind::Proper(au)
714    }
715
716    /// Checks if the symmetry element is spatially an identity element and has the specified
717    /// antiunitary attribute.
718    ///
719    /// More specifically, this checks if the combination $`\hat{\gamma} \hat{C}_n^k`$ is the
720    /// identity, and if $`\hat{\alpha}`$ is as specified by `au`.
721    ///
722    /// # Arguments
723    ///
724    /// * `au` - An `Option` for the desired antiunitary kind.
725    ///
726    /// # Returns
727    ///
728    /// A boolean indicating if this symmetry element is spatially an identity element and has the
729    /// specified antiunitary attribute.
730    #[must_use]
731    pub fn is_o3_identity(&self, au: Option<AntiunitaryKind>) -> bool {
732        self.kind == SymmetryElementKind::Proper(au)
733            && self
734                .proper_fraction
735                .map(|frac| frac.is_zero())
736                .or_else(|| {
737                    self.proper_angle.map(|proper_angle| {
738                        approx::abs_diff_eq!(proper_angle, 0.0, epsilon = self.threshold,)
739                    })
740                })
741                .unwrap_or(false)
742    }
743
744    /// Checks if the symmetry element is spatially an inversion centre and has the specified
745    /// antiunitary attribute.
746    ///
747    /// More specifically, this checks if the combination $`\hat{\gamma} \hat{C}_n^k`$ is the
748    /// spatial inversion, and if $`\hat{\alpha}`$ is as specified by `au`.
749    ///
750    /// # Arguments
751    ///
752    /// * `au` - An `Option` for the desired antiunitary kind.
753    ///
754    /// # Returns
755    ///
756    /// A boolean indicating if this symmetry element is spatially an inversion centre and has the
757    /// specified antiunitary attribute.
758    #[must_use]
759    pub fn is_o3_inversion_centre(&self, au: Option<AntiunitaryKind>) -> bool {
760        match self.kind {
761            SymmetryElementKind::ImproperMirrorPlane(sig_au) => {
762                sig_au == au
763                    && self
764                        .proper_fraction
765                        .map(|frac| frac == F::new(1u32, 2u32))
766                        .or_else(|| {
767                            self.proper_angle.map(|proper_angle| {
768                                approx::abs_diff_eq!(
769                                    proper_angle,
770                                    std::f64::consts::PI,
771                                    epsilon = self.threshold,
772                                )
773                            })
774                        })
775                        .unwrap_or(false)
776            }
777            SymmetryElementKind::ImproperInversionCentre(inv_au) => {
778                inv_au == au
779                    && self
780                        .proper_fraction
781                        .map(|frac| frac.is_zero())
782                        .or_else(|| {
783                            self.proper_angle.map(|proper_angle| {
784                                approx::abs_diff_eq!(proper_angle, 0.0, epsilon = self.threshold,)
785                            })
786                        })
787                        .unwrap_or(false)
788            }
789            _ => false,
790        }
791    }
792
793    /// Checks if the symmetry element is spatially a binary rotation axis and has the specified
794    /// antiunitary attribute.
795    ///
796    /// More specifically, this checks if the combination $`\hat{\gamma} \hat{C}_n^k`$ is a binary
797    /// rotation, and if $`\hat{\alpha}`$ is as specified by `au`.
798    ///
799    /// # Arguments
800    ///
801    /// * `au` - An `Option` for the desired antiunitary kind.
802    ///
803    /// # Returns
804    ///
805    /// A boolean indicating if this symmetry element is spatially a binary rotation axis and has
806    /// the specified antiunitary attribute.
807    #[must_use]
808    pub fn is_o3_binary_rotation_axis(&self, au: Option<AntiunitaryKind>) -> bool {
809        self.kind == SymmetryElementKind::Proper(au)
810            && self
811                .proper_fraction
812                .map(|frac| frac == F::new(1u32, 2u32))
813                .or_else(|| {
814                    self.proper_angle.map(|proper_angle| {
815                        approx::abs_diff_eq!(
816                            proper_angle,
817                            std::f64::consts::PI,
818                            epsilon = self.threshold,
819                        )
820                    })
821                })
822                .unwrap_or(false)
823    }
824
825    /// Checks if the symmetry element is spatially a mirror plane and has the specified
826    /// antiunitary attribute.
827    ///
828    /// More specifically, this checks if the combination $`\hat{\gamma} \hat{C}_n^k`$ is a
829    /// reflection, and if $`\hat{\alpha}`$ is as specified by `au`.
830    ///
831    /// # Arguments
832    ///
833    /// * `au` - An `Option` for the desired antiunitary kind.
834    ///
835    /// # Returns
836    ///
837    /// A boolean indicating if this symmetry element is spatially a mirror plane and has the
838    /// specified antiunitary attribute.
839    #[must_use]
840    pub fn is_o3_mirror_plane(&self, au: Option<AntiunitaryKind>) -> bool {
841        match self.kind {
842            SymmetryElementKind::ImproperMirrorPlane(sig_au) => {
843                sig_au == au
844                    && self
845                        .proper_fraction
846                        .map(|frac| frac.is_zero())
847                        .or_else(|| {
848                            self.proper_angle.map(|proper_angle| {
849                                approx::abs_diff_eq!(proper_angle, 0.0, epsilon = self.threshold,)
850                            })
851                        })
852                        .unwrap_or(false)
853            }
854            SymmetryElementKind::ImproperInversionCentre(inv_au) => {
855                inv_au == au
856                    && self
857                        .proper_fraction
858                        .map(|frac| frac == F::new(1u32, 2u32))
859                        .or_else(|| {
860                            self.proper_angle.map(|proper_angle| {
861                                approx::abs_diff_eq!(
862                                    proper_angle,
863                                    std::f64::consts::PI,
864                                    epsilon = self.threshold,
865                                )
866                            })
867                        })
868                        .unwrap_or(false)
869            }
870            _ => false,
871        }
872    }
873
874    /// Returns the full symbol for this symmetry element, which does not classify certain
875    /// improper rotation axes into inversion centres or mirror planes, but which does simplify
876    /// the power/order ratio, and which displays only the absolute value of the power since
877    /// symmetry elements do not distinguish senses of rotations since rotations of oposite
878    /// directions are inverses of each other, both of which must exist in the group.
879    ///
880    /// Some additional symbols that can be unconventional include:
881    ///
882    /// * `θ`: time reversal,
883    /// * `(Σ)`: the spatial part is in homotopy class 0 of $`\mathsf{SU}'(2)`$,
884    /// * `(QΣ)`: the spatial part is in homotopy class 1 of $`\mathsf{SU}'(2)`$.
885    ///
886    /// See [`RotationGroup`] for further information.
887    ///
888    /// # Returns
889    ///
890    /// The full symbol for this symmetry element.
891    #[must_use]
892    pub fn get_full_symbol(&self) -> String {
893        let tr_sym = if self.kind.contains_time_reversal() {
894            "θ·"
895        } else {
896            ""
897        };
898        let main_symbol: String = match self.kind {
899            SymmetryElementKind::Proper(_) => {
900                format!("{tr_sym}C")
901            }
902            SymmetryElementKind::ImproperMirrorPlane(_) => {
903                if *self.raw_proper_order() != ElementOrder::Inf
904                    && (*self
905                        .proper_fraction
906                        .expect("No proper fractions found for a finite-order element.")
907                        .numer()
908                        .expect("Unable to extract the numerator of the proper fraction.")
909                        == 1
910                        || self
911                            .proper_fraction
912                            .expect("No proper fractions found for a finite-order element.")
913                            .is_zero())
914                {
915                    format!("{tr_sym}S")
916                } else {
917                    format!("{tr_sym}σC")
918                }
919            }
920            SymmetryElementKind::ImproperInversionCentre(_) => {
921                if *self.raw_proper_order() != ElementOrder::Inf
922                    && (*self
923                        .proper_fraction
924                        .expect("No proper fractions found for a finite-order element.")
925                        .numer()
926                        .expect("Unable to extract the numerator of the proper fraction.")
927                        == 1
928                        || self
929                            .proper_fraction
930                            .expect("No proper fractions found for a finite-order element.")
931                            .is_zero())
932                {
933                    format!("{tr_sym}Ṡ")
934                } else {
935                    format!("{tr_sym}iC")
936                }
937            }
938        };
939
940        let su2_sym = if self.is_su2_class_1() {
941            "(QΣ)"
942        } else if self.is_su2() {
943            "(Σ)"
944        } else {
945            ""
946        };
947
948        if let Some(proper_fraction) = self.proper_fraction {
949            if proper_fraction.is_zero() {
950                format!("{main_symbol}1{su2_sym}")
951            } else {
952                let proper_order = proper_fraction
953                    .denom()
954                    .expect("Unable to extract the denominator of the proper fraction.")
955                    .to_string();
956                let proper_power = {
957                    let pow = *proper_fraction
958                        .numer()
959                        .expect("Unable to extract the numerator of the proper fraction.");
960                    if pow > 1 {
961                        format!("^{pow}")
962                    } else {
963                        String::new()
964                    }
965                };
966                format!("{main_symbol}{proper_order}{proper_power}{su2_sym}")
967            }
968        } else {
969            assert_eq!(*self.raw_proper_order(), ElementOrder::Inf);
970            let proper_angle = if let Some(proper_angle) = self.proper_angle {
971                format!("({:+.3})", proper_angle.abs())
972            } else {
973                String::new()
974            };
975            format!(
976                "{main_symbol}{}{proper_angle}{su2_sym}",
977                self.raw_proper_order()
978            )
979        }
980    }
981
982    /// Returns the simplified symbol for this symmetry element, which classifies special symmetry
983    /// elements (identity, inversion centre, time reversal, mirror planes), and which simplifies
984    /// the power/order ratio and displays only the absolute value of the power since symmetry
985    /// elements do not distinguish senses of rotations, as rotations of oposite directions are
986    /// inverses of each other, both of which must exist in the group.
987    ///
988    /// # Returns
989    ///
990    /// The simplified symbol for this symmetry element.
991    #[must_use]
992    pub fn get_simplified_symbol(&self) -> String {
993        let (main_symbol, needs_power) = match self.kind {
994            SymmetryElementKind::Proper(au) => {
995                if self.is_o3_identity(au) {
996                    match au {
997                        None => ("E".to_owned(), false),
998                        Some(AntiunitaryKind::TimeReversal) => ("θ".to_owned(), false),
999                        Some(AntiunitaryKind::ComplexConjugation) => ("K".to_owned(), false),
1000                    }
1001                } else {
1002                    let au_sym = match au {
1003                        None => "",
1004                        Some(AntiunitaryKind::TimeReversal) => "θ·",
1005                        Some(AntiunitaryKind::ComplexConjugation) => "K·",
1006                    };
1007                    (format!("{au_sym}C"), true)
1008                }
1009            }
1010            SymmetryElementKind::ImproperMirrorPlane(au) => {
1011                let au_sym = match au {
1012                    None => "",
1013                    Some(AntiunitaryKind::TimeReversal) => "θ·",
1014                    Some(AntiunitaryKind::ComplexConjugation) => "K·",
1015                };
1016                if self.is_o3_mirror_plane(au) {
1017                    (format!("{au_sym}σ"), false)
1018                } else if self.is_o3_inversion_centre(au) {
1019                    (format!("{au_sym}i"), false)
1020                } else if *self.raw_proper_order() == ElementOrder::Inf
1021                    || *self
1022                        .proper_fraction
1023                        .expect("No proper fractions found for a finite-order element.")
1024                        .numer()
1025                        .expect("Unable to extract the numerator of the proper fraction.")
1026                        == 1
1027                {
1028                    (format!("{au_sym}S"), false)
1029                } else {
1030                    (format!("{au_sym}σC"), true)
1031                }
1032            }
1033            SymmetryElementKind::ImproperInversionCentre(au) => {
1034                let au_sym = match au {
1035                    None => "",
1036                    Some(AntiunitaryKind::TimeReversal) => "θ·",
1037                    Some(AntiunitaryKind::ComplexConjugation) => "K·",
1038                };
1039                if self.is_o3_mirror_plane(au) {
1040                    (format!("{au_sym}σ"), false)
1041                } else if self.is_o3_inversion_centre(au) {
1042                    (format!("{au_sym}i"), false)
1043                } else if *self.raw_proper_order() == ElementOrder::Inf
1044                    || *self
1045                        .proper_fraction
1046                        .expect("No proper fractions found for a finite-order element.")
1047                        .numer()
1048                        .expect("Unable to extract the numerator of the proper fraction.")
1049                        == 1
1050                {
1051                    (format!("{au_sym}Ṡ"), false)
1052                } else {
1053                    (format!("{au_sym}iC"), true)
1054                }
1055            }
1056        };
1057
1058        let su2_sym = if self.is_su2_class_1() {
1059            "(QΣ)"
1060        } else if self.is_su2() {
1061            "(Σ)"
1062        } else {
1063            ""
1064        };
1065
1066        if let Some(proper_fraction) = self.proper_fraction {
1067            let au = self.contains_antiunitary();
1068            let proper_order = if self.is_o3_identity(au)
1069                || self.is_o3_inversion_centre(au)
1070                || self.is_o3_mirror_plane(au)
1071            {
1072                String::new()
1073            } else {
1074                proper_fraction
1075                    .denom()
1076                    .expect("Unable to extract the denominator of the proper fraction.")
1077                    .to_string()
1078            };
1079
1080            let proper_power = if needs_power {
1081                let pow = *proper_fraction
1082                    .numer()
1083                    .expect("Unable to extract the numerator of the proper fraction.");
1084                if pow > 1 {
1085                    format!("^{pow}")
1086                } else {
1087                    String::new()
1088                }
1089            } else {
1090                String::new()
1091            };
1092            format!(
1093                "{main_symbol}{}{proper_order}{proper_power}{}{su2_sym}",
1094                self.additional_superscript, self.additional_subscript
1095            )
1096        } else {
1097            assert_eq!(*self.raw_proper_order(), ElementOrder::Inf);
1098            let proper_angle = if let Some(proper_angle) = self.proper_angle {
1099                format!("({:+.3})", proper_angle.abs())
1100            } else {
1101                String::new()
1102            };
1103            format!(
1104                "{main_symbol}{}{}{proper_angle}{}{su2_sym}",
1105                self.additional_superscript,
1106                *self.raw_proper_order(),
1107                self.additional_subscript
1108            )
1109        }
1110    }
1111
1112    /// Returns the simplified symbol for this symmetry element, which classifies special symmetry
1113    /// elements (identity, inversion centre, mirror planes), and which simplifies the power/order
1114    /// ratio and displays only the absolute value of the power since symmetry elements do not
1115    /// distinguish senses of rotations. Rotations of oposite directions are inverses of each
1116    /// other, both of which must exist in the group.
1117    ///
1118    /// # Returns
1119    ///
1120    /// The simplified symbol for this symmetry element.
1121    #[must_use]
1122    pub fn get_simplified_symbol_signed_power(&self) -> String {
1123        let (main_symbol, needs_power) = match self.kind {
1124            SymmetryElementKind::Proper(au) => {
1125                if self.is_o3_identity(au) {
1126                    match au {
1127                        None => ("E".to_owned(), false),
1128                        Some(AntiunitaryKind::TimeReversal) => ("θ".to_owned(), false),
1129                        Some(AntiunitaryKind::ComplexConjugation) => ("K".to_owned(), false),
1130                    }
1131                } else {
1132                    let au_sym = match au {
1133                        None => "",
1134                        Some(AntiunitaryKind::TimeReversal) => "θ·",
1135                        Some(AntiunitaryKind::ComplexConjugation) => "K·",
1136                    };
1137                    (format!("{au_sym}C"), true)
1138                }
1139            }
1140            SymmetryElementKind::ImproperMirrorPlane(au) => {
1141                let au_sym = match au {
1142                    None => "",
1143                    Some(AntiunitaryKind::TimeReversal) => "θ·",
1144                    Some(AntiunitaryKind::ComplexConjugation) => "K·",
1145                };
1146                if self.is_o3_mirror_plane(au) {
1147                    (format!("{au_sym}σ"), false)
1148                } else if self.is_o3_inversion_centre(au) {
1149                    (format!("{au_sym}i"), false)
1150                } else if *self.raw_proper_order() == ElementOrder::Inf
1151                    || *self
1152                        .proper_fraction
1153                        .expect("No proper fractions found for a finite-order element.")
1154                        .numer()
1155                        .expect("Unable to extract the numerator of the proper fraction.")
1156                        == 1
1157                {
1158                    (format!("{au_sym}S"), true)
1159                } else {
1160                    (format!("{au_sym}σC"), true)
1161                }
1162            }
1163            SymmetryElementKind::ImproperInversionCentre(au) => {
1164                let au_sym = match au {
1165                    None => "",
1166                    Some(AntiunitaryKind::TimeReversal) => "θ·",
1167                    Some(AntiunitaryKind::ComplexConjugation) => "K·",
1168                };
1169                if self.is_o3_mirror_plane(au) {
1170                    (format!("{au_sym}σ"), false)
1171                } else if self.is_o3_inversion_centre(au) {
1172                    (format!("{au_sym}i"), false)
1173                } else if *self.raw_proper_order() == ElementOrder::Inf
1174                    || *self
1175                        .proper_fraction
1176                        .expect("No proper fractions found for a finite-order element.")
1177                        .numer()
1178                        .expect("Unable to extract the numerator of the proper fraction.")
1179                        == 1
1180                {
1181                    (format!("{au_sym}Ṡ"), true)
1182                } else {
1183                    (format!("{au_sym}iC"), true)
1184                }
1185            }
1186        };
1187
1188        let su2_sym = if self.is_su2_class_1() {
1189            "(QΣ)"
1190        } else if self.is_su2() {
1191            "(Σ)"
1192        } else {
1193            ""
1194        };
1195
1196        if let Some(proper_fraction) = self.proper_fraction {
1197            let au = self.contains_antiunitary();
1198            let proper_order = if self.is_o3_identity(au)
1199                || self.is_o3_inversion_centre(au)
1200                || self.is_o3_mirror_plane(au)
1201            {
1202                String::new()
1203            } else {
1204                proper_fraction
1205                    .denom()
1206                    .expect("Unable to extract the denominator of the proper fraction.")
1207                    .to_string()
1208            };
1209
1210            let proper_power = if needs_power {
1211                let pow = *proper_fraction
1212                    .numer()
1213                    .expect("Unable to extract the numerator of the proper fraction.");
1214                if !geometry::check_standard_positive_pole(
1215                    &self.proper_rotation_pole(),
1216                    self.threshold,
1217                ) {
1218                    format!("^(-{pow})")
1219                } else if pow > 1 {
1220                    format!("^{pow}")
1221                } else {
1222                    String::new()
1223                }
1224            } else {
1225                String::new()
1226            };
1227            format!(
1228                "{main_symbol}{}{proper_order}{proper_power}{}{su2_sym}",
1229                self.additional_superscript, self.additional_subscript
1230            )
1231        } else {
1232            assert_eq!(*self.raw_proper_order(), ElementOrder::Inf);
1233            let proper_angle = if let Some(proper_angle) = self.proper_angle {
1234                format!("({:+.3})", proper_angle.abs())
1235            } else {
1236                String::new()
1237            };
1238            format!(
1239                "{main_symbol}{}{}{proper_angle}{}{su2_sym}",
1240                self.additional_superscript,
1241                *self.raw_proper_order(),
1242                self.additional_subscript
1243            )
1244        }
1245    }
1246
1247    /// Returns a copy of the current improper symmetry element that has been converted to the
1248    /// required improper kind. For $`\mathsf{SU}'(2)`$ elements, the conversion will be carried
1249    /// out in the same homotopy class. The antiunitary part will be kept unchanged.
1250    ///
1251    /// To convert between the two improper kinds, we essentially seek integers
1252    /// $`n, n' \in \mathbb{N}_{+}`$ and $`k \in \mathbb{Z}/n\mathbb{Z}`$,
1253    /// $`k' \in \mathbb{Z}/n'\mathbb{Z}`$, such that
1254    ///
1255    /// ```math
1256    /// \sigma C_n^k = i C_{n'}^{k'},
1257    /// ```
1258    ///
1259    /// where the axes of all involved elements are parallel. By noting that
1260    /// $`\sigma = i C_2`$ and that $`k`$ and $`k'`$ must have opposite signs, we can easily show
1261    /// that, for $`k \ge 0, k' < 0`$,
1262    ///
1263    /// ```math
1264    /// \begin{aligned}
1265    ///     n' &= \frac{2n}{\operatorname{gcd}(2n, n - 2k)},\\
1266    ///     k' &= -\frac{n - 2k}{\operatorname{gcd}(2n, n - 2k)},
1267    /// \end{aligned}
1268    /// ```
1269    ///
1270    /// whereas for $`k < 0, k' \ge 0`$,
1271    ///
1272    /// ```math
1273    /// \begin{aligned}
1274    ///     n' &= \frac{2n}{\operatorname{gcd}(2n, n + 2k)},\\
1275    ///     k' &= \frac{n + 2k}{\operatorname{gcd}(2n, n + 2k)}.
1276    /// \end{aligned}
1277    /// ```
1278    ///
1279    /// The above relations are self-inverse. It can be further shown that
1280    /// $`\operatorname{gcd}(n', k') = 1`$. Hence, for symmetry *element* conversions, we can simply
1281    /// take $`k' = 1`$. This is because a symmetry element plays the role of a generator, and the
1282    /// coprimality of $`n'`$ and $`k'`$ means that $`i C_{n'}^{1}`$ is as valid a generator as
1283    /// $`i C_{n'}^{k'}`$.
1284    ///
1285    /// # Arguments
1286    ///
1287    /// * `improper_kind` - The improper kind to which `self` is to be converted. There is no need
1288    /// to make sure the time reversal specification in `improper_kind` matches that of `self` as
1289    /// the conversion will take care of this.
1290    /// * `preserves_power` - Boolean indicating if the proper rotation power $`k'`$ should be
1291    /// preserved or should be set to $`1`$.
1292    ///
1293    /// # Returns
1294    ///
1295    /// A copy of the current improper symmetry element that has been converted.
1296    ///
1297    /// # Panics
1298    ///
1299    /// Panics when `self` is not an improper element, or when `improper_kind` is not one of the
1300    /// improper variants.
1301    #[must_use]
1302    pub fn convert_to_improper_kind(
1303        &self,
1304        improper_kind: &SymmetryElementKind,
1305        preserves_power: bool,
1306    ) -> Self {
1307        let au = self.contains_antiunitary();
1308        assert!(
1309            !self.is_o3_proper(au),
1310            "Only improper elements can be converted."
1311        );
1312        let improper_kind = improper_kind.to_antiunitary(au);
1313        assert!(
1314            !matches!(improper_kind, SymmetryElementKind::Proper(_)),
1315            "`improper_kind` must be one of the improper variants."
1316        );
1317
1318        // self.kind and improper_kind must now have the same antiunitary part.
1319
1320        if self.kind == improper_kind {
1321            return self.clone();
1322        }
1323
1324        let (dest_order, dest_proper_power) = match *self.raw_proper_order() {
1325            ElementOrder::Int(_) => {
1326                let proper_fraction = self
1327                    .proper_fraction
1328                    .expect("No proper fractions found for an element with integer order.");
1329                let n = *proper_fraction.denom().unwrap();
1330                let k = if proper_fraction.is_sign_negative() {
1331                    -i32::try_from(*proper_fraction.numer().unwrap_or_else(|| {
1332                        panic!("Unable to retrieve the numerator of {proper_fraction:?}.")
1333                    }))
1334                    .expect("Unable to convert the numerator of the proper fraction to `i32`.")
1335                } else {
1336                    i32::try_from(*proper_fraction.numer().unwrap_or_else(|| {
1337                        panic!("Unable to retrieve the numerator of {proper_fraction:?}.")
1338                    }))
1339                    .expect("Unable to convert the numerator of the proper fraction to `i32`.")
1340                };
1341
1342                if k >= 0 {
1343                    // k >= 0, k2 < 0
1344                    let n_m_2k = n
1345                        .checked_sub(2 * k.unsigned_abs())
1346                        .expect("The value of `n - 2k` is negative.");
1347                    let n2 = ElementOrder::Int(2 * n / (gcd(2 * n, n_m_2k)));
1348                    let k2: i32 = if preserves_power {
1349                        -i32::try_from(n_m_2k / gcd(2 * n, n_m_2k))
1350                            .expect("Unable to convert `k'` to `i32`.")
1351                    } else {
1352                        1
1353                    };
1354                    (n2, k2)
1355                } else {
1356                    // k < 0, k2 >= 0
1357                    let n_p_2k = n
1358                        .checked_sub(2 * k.unsigned_abs())
1359                        .expect("The value of `n + 2k` is negative.");
1360                    let n2 = ElementOrder::Int(2 * n / (gcd(2 * n, n_p_2k)));
1361                    let k2: i32 = if preserves_power {
1362                        i32::try_from(n_p_2k / (gcd(2 * n, n_p_2k)))
1363                            .expect("Unable to convert `k'` to `i32`.")
1364                    } else {
1365                        1
1366                    };
1367                    (n2, k2)
1368                }
1369            }
1370            ElementOrder::Inf => (ElementOrder::Inf, 1),
1371        };
1372
1373        match dest_order {
1374            ElementOrder::Int(_) => Self::builder()
1375                .threshold(self.threshold)
1376                .proper_order(dest_order)
1377                .proper_power(dest_proper_power)
1378                .raw_axis(self.raw_axis)
1379                .kind(improper_kind)
1380                .rotation_group(self.rotation_group.clone())
1381                .generator(self.generator)
1382                .additional_superscript(self.additional_superscript.clone())
1383                .additional_subscript(self.additional_subscript.clone())
1384                .build()
1385                .expect("Unable to construct a symmetry element."),
1386
1387            ElementOrder::Inf => {
1388                if let Some(ang) = self.proper_angle {
1389                    Self::builder()
1390                        .threshold(self.threshold)
1391                        .proper_order(dest_order)
1392                        .proper_power(dest_proper_power)
1393                        .proper_angle(-std::f64::consts::PI + ang)
1394                        .raw_axis(self.raw_axis)
1395                        .kind(improper_kind)
1396                        .rotation_group(self.rotation_group.clone())
1397                        .generator(self.generator)
1398                        .additional_superscript(self.additional_superscript.clone())
1399                        .additional_subscript(self.additional_subscript.clone())
1400                        .build()
1401                        .expect("Unable to construct a symmetry element.")
1402                } else {
1403                    Self::builder()
1404                        .threshold(self.threshold)
1405                        .proper_order(dest_order)
1406                        .proper_power(dest_proper_power)
1407                        .raw_axis(self.raw_axis)
1408                        .kind(improper_kind)
1409                        .rotation_group(self.rotation_group.clone())
1410                        .generator(self.generator)
1411                        .additional_superscript(self.additional_superscript.clone())
1412                        .additional_subscript(self.additional_subscript.clone())
1413                        .build()
1414                        .expect("Unable to construct a symmetry element.")
1415                }
1416            }
1417        }
1418    }
1419
1420    /// Converts the symmetry element to one with the desired time-reversal property.
1421    ///
1422    /// # Arguments
1423    ///
1424    /// * `tr` - A boolean indicating if time reversal is to be included.
1425    ///
1426    /// # Returns
1427    ///
1428    /// A new symmetry element with or without time reversal as indicated by `tr`.
1429    pub fn to_tr(&self, tr: bool) -> Self {
1430        let mut c_self = self.clone();
1431        c_self.kind = c_self.kind.to_tr(tr);
1432        c_self
1433    }
1434
1435    /// Convert the proper rotation of the current element to one in $`\mathsf{SU}(2)`$.
1436    ///
1437    /// # Arguments
1438    ///
1439    /// * `normal` - A boolean indicating whether the resultant $`\mathsf{SU}(2)`$ proper rotation
1440    /// is of homotopy class 0 (`true`) or 1 (`false`) when connected to the identity.
1441    ///
1442    /// # Returns
1443    ///
1444    /// A symmetry element in $`\mathsf{SU}(2)`$, or `None` if the current symmetry element
1445    /// is already in $`\mathsf{SU}(2)`$.
1446    pub fn to_su2(&self, normal: bool) -> Option<Self> {
1447        if self.is_su2() {
1448            None
1449        } else {
1450            let mut element = self.clone();
1451            element.rotation_group = RotationGroup::SU2(normal);
1452            Some(element)
1453        }
1454    }
1455
1456    /// The closeness of the symmetry element's axis to one of the three space-fixed Cartesian axes.
1457    ///
1458    /// # Returns
1459    ///
1460    /// A tuple of two values:
1461    /// - A value $`\gamma \in [0, 1-1/\sqrt{3}]`$ indicating how close the axis is to one of the
1462    /// three Cartesian axes. The closer $`\gamma`$ is to $`0`$, the closer the alignment.
1463    /// - An index for the closest axis: `0` for $`z`$, `1` for $`y`$, `2` for $`x`$.
1464    ///
1465    /// # Panics
1466    ///
1467    /// Panics when $`\gamma`$ is outside the required closed interval $`[0, 1-1/\sqrt{3}]`$ by
1468    /// more than the threshold value in `self`.
1469    #[must_use]
1470    pub fn closeness_to_cartesian_axes(&self) -> (f64, usize) {
1471        let pos_axis = self.standard_positive_axis();
1472        let rev_pos_axis = Vector3::new(pos_axis[2], pos_axis[1], pos_axis[0]);
1473        let (amax_arg, amax_val) = rev_pos_axis.abs().argmax();
1474        let axis_closeness = 1.0 - amax_val;
1475        let thresh = self.threshold;
1476        assert!(
1477            -thresh <= axis_closeness && axis_closeness <= (1.0 - 1.0 / 3.0f64.sqrt() + thresh)
1478        );
1479
1480        // closest axis: 0 for z, 1 for y, 2 for x
1481        // This is so that z axes are preferred.
1482        let closest_axis = amax_arg;
1483        (axis_closeness, closest_axis)
1484    }
1485}
1486
1487impl fmt::Debug for SymmetryElement {
1488    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1489        let signed_axis = self.signed_axis();
1490        write!(
1491            f,
1492            "{}({:+.3}, {:+.3}, {:+.3})",
1493            self.get_full_symbol(),
1494            signed_axis[0] + 0.0,
1495            signed_axis[1] + 0.0,
1496            signed_axis[2] + 0.0
1497        )
1498    }
1499}
1500
1501impl fmt::Display for SymmetryElement {
1502    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1503        let au = self.contains_antiunitary();
1504        if self.is_o3_identity(au) || self.is_o3_inversion_centre(au) {
1505            write!(f, "{}", self.get_simplified_symbol())
1506        } else {
1507            let signed_axis = self.signed_axis();
1508            write!(
1509                f,
1510                "{}({:+.3}, {:+.3}, {:+.3})",
1511                self.get_simplified_symbol(),
1512                signed_axis[0] + 0.0,
1513                signed_axis[1] + 0.0,
1514                signed_axis[2] + 0.0
1515            )
1516        }
1517    }
1518}
1519
1520impl PartialEq for SymmetryElement {
1521    /// Two symmetry elements are equal if and only if the following conditions are all satisfied:
1522    ///
1523    /// * they are both in the same rotation group and belong to the same homotopy class;
1524    /// * they are both proper or improper;
1525    /// * they both have the same antiunitary properties;
1526    /// * their axes are either parallel or anti-parallel;
1527    /// * their proper rotation angles have equal absolute values.
1528    ///
1529    /// For improper elements, proper rotation angles are taken in the inversion centre convention.
1530    ///
1531    /// Thus, symmetry element equality is less strict than symmetry operation equality. This is so
1532    /// that parallel or anti-parallel symmetry elements with the same spatial and time-reversal
1533    /// parities and angle of rotation are deemed identical, thus facilitating symmetry detection
1534    /// where one does not yet care much about directions of rotations.
1535    #[allow(clippy::too_many_lines)]
1536    fn eq(&self, other: &Self) -> bool {
1537        if self.rotation_group != other.rotation_group {
1538            // Different rotation groups or homotopy classes.
1539            return false;
1540        }
1541
1542        if self.contains_antiunitary() != other.contains_antiunitary() {
1543            // Different anti-unitary kinds.
1544            return false;
1545        }
1546
1547        let au = self.contains_antiunitary();
1548
1549        if self.is_o3_proper(au) != other.is_o3_proper(au) {
1550            // Different spatial parities.
1551            return false;
1552        }
1553
1554        if self.is_o3_identity(au) && other.is_o3_identity(au) {
1555            // Both are spatial identity.
1556            // assert_eq!(
1557            //     misc::calculate_hash(self),
1558            //     misc::calculate_hash(other),
1559            //     "{self} and {other} have unequal hashes."
1560            // );
1561            // return true;
1562            return misc::calculate_hash(self) == misc::calculate_hash(other);
1563        }
1564
1565        if self.is_o3_inversion_centre(au) && other.is_o3_inversion_centre(au) {
1566            // Both are spatial inversion centre.
1567            // assert_eq!(
1568            //     misc::calculate_hash(self),
1569            //     misc::calculate_hash(other),
1570            //     "{self} and {other} have unequal hashes."
1571            // );
1572            // return true;
1573            return misc::calculate_hash(self) == misc::calculate_hash(other);
1574        }
1575
1576        let thresh = (self.threshold * other.threshold).sqrt();
1577
1578        let result = if self.is_o3_proper(au) {
1579            // Proper.
1580
1581            // Parallel or anti-parallel axes.
1582            let similar_poles = approx::relative_eq!(
1583                geometry::get_standard_positive_pole(&self.raw_axis, thresh),
1584                geometry::get_standard_positive_pole(&other.raw_axis, thresh),
1585                epsilon = thresh,
1586                max_relative = thresh
1587            );
1588
1589            // Same angle of rotation (irrespective of signs).
1590            let similar_angles = match (*self.raw_proper_order(), *other.raw_proper_order()) {
1591                (ElementOrder::Inf, ElementOrder::Inf) => {
1592                    match (self.proper_angle, other.proper_angle) {
1593                        (Some(s_angle), Some(o_angle)) => {
1594                            approx::relative_eq!(
1595                                s_angle.abs(),
1596                                o_angle.abs(),
1597                                epsilon = thresh,
1598                                max_relative = thresh
1599                            )
1600                        }
1601                        (None, None) => similar_poles,
1602                        _ => false,
1603                    }
1604                }
1605                (ElementOrder::Int(_), ElementOrder::Int(_)) => {
1606                    let c_proper_fraction = self
1607                        .proper_fraction
1608                        .expect("Proper fraction for `self` not found.");
1609                    let o_proper_fraction = other
1610                        .proper_fraction
1611                        .expect("Proper fraction for `other` not found.");
1612                    c_proper_fraction.abs() == o_proper_fraction.abs()
1613                }
1614                _ => false,
1615            };
1616
1617            similar_poles && similar_angles
1618        } else {
1619            // Improper => convert to inversion-centre convention.
1620            let inv_au = SymmetryElementKind::ImproperInversionCentre(au);
1621            let c_self = self.convert_to_improper_kind(&inv_au, false);
1622            let c_other = other.convert_to_improper_kind(&inv_au, false);
1623
1624            // Parallel or anti-parallel axes.
1625            let similar_poles = approx::relative_eq!(
1626                geometry::get_standard_positive_pole(&c_self.raw_axis, thresh),
1627                geometry::get_standard_positive_pole(&c_other.raw_axis, thresh),
1628                epsilon = thresh,
1629                max_relative = thresh
1630            );
1631
1632            // Same angle of rotation (irrespective of signs).
1633            let similar_angles = match (*c_self.raw_proper_order(), *c_other.raw_proper_order()) {
1634                (ElementOrder::Inf, ElementOrder::Inf) => {
1635                    match (c_self.proper_angle, c_other.proper_angle) {
1636                        (Some(s_angle), Some(o_angle)) => {
1637                            approx::relative_eq!(
1638                                s_angle.abs(),
1639                                o_angle.abs(),
1640                                epsilon = thresh,
1641                                max_relative = thresh
1642                            )
1643                        }
1644                        (None, None) => similar_poles,
1645                        _ => false,
1646                    }
1647                }
1648                (ElementOrder::Int(_), ElementOrder::Int(_)) => {
1649                    let c_proper_fraction = c_self
1650                        .proper_fraction
1651                        .expect("Proper fraction for `c_self` not found.");
1652                    let o_proper_fraction = c_other
1653                        .proper_fraction
1654                        .expect("Proper fraction for `c_other` not found.");
1655                    c_proper_fraction.abs() == o_proper_fraction.abs()
1656                }
1657                _ => false,
1658            };
1659
1660            similar_poles && similar_angles
1661        };
1662
1663        // if result {
1664        //     assert_eq!(
1665        //         misc::calculate_hash(self),
1666        //         misc::calculate_hash(other),
1667        //         "`{self}` and `{other}` have unequal hashes."
1668        //     );
1669        // }
1670        result && (misc::calculate_hash(self) == misc::calculate_hash(other))
1671    }
1672}
1673
1674impl Eq for SymmetryElement {}
1675
1676impl Hash for SymmetryElement {
1677    fn hash<H: Hasher>(&self, state: &mut H) {
1678        self.rotation_group.hash(state);
1679
1680        let au = self.contains_antiunitary();
1681        au.hash(state);
1682
1683        self.is_o3_proper(au).hash(state);
1684
1685        if self.is_o3_identity(au) || self.is_o3_inversion_centre(au) {
1686            true.hash(state);
1687        } else if self.kind == SymmetryElementKind::ImproperMirrorPlane(au) {
1688            let c_self = self
1689                .convert_to_improper_kind(&SymmetryElementKind::ImproperInversionCentre(au), false);
1690            let pole = geometry::get_standard_positive_pole(&c_self.raw_axis, c_self.threshold);
1691            pole[0]
1692                .round_factor(self.threshold)
1693                .integer_decode()
1694                .hash(state);
1695            pole[1]
1696                .round_factor(self.threshold)
1697                .integer_decode()
1698                .hash(state);
1699            pole[2]
1700                .round_factor(self.threshold)
1701                .integer_decode()
1702                .hash(state);
1703            if let ElementOrder::Inf = *c_self.raw_proper_order() {
1704                if let Some(angle) = c_self.proper_angle {
1705                    angle
1706                        .abs()
1707                        .round_factor(self.threshold)
1708                        .integer_decode()
1709                        .hash(state);
1710                } else {
1711                    0.hash(state);
1712                }
1713            } else {
1714                c_self
1715                    .proper_fraction
1716                    .expect("No proper fractions for `c_self` found.")
1717                    .abs()
1718                    .hash(state);
1719            };
1720        } else {
1721            let pole = geometry::get_standard_positive_pole(&self.raw_axis, self.threshold);
1722            pole[0]
1723                .round_factor(self.threshold)
1724                .integer_decode()
1725                .hash(state);
1726            pole[1]
1727                .round_factor(self.threshold)
1728                .integer_decode()
1729                .hash(state);
1730            pole[2]
1731                .round_factor(self.threshold)
1732                .integer_decode()
1733                .hash(state);
1734            if let ElementOrder::Inf = *self.raw_proper_order() {
1735                if let Some(angle) = self.proper_angle {
1736                    angle
1737                        .abs()
1738                        .round_factor(self.threshold)
1739                        .integer_decode()
1740                        .hash(state);
1741                } else {
1742                    0.hash(state);
1743                }
1744            } else {
1745                self.proper_fraction
1746                    .expect("No proper fractions for `self` found.")
1747                    .abs()
1748                    .hash(state);
1749            };
1750        };
1751    }
1752}
1753
1754/// Time-reversal antiunitary kind.
1755pub const TR: AntiunitaryKind = AntiunitaryKind::TimeReversal;
1756
1757/// Proper rotation symmetry element kind.
1758pub const ROT: SymmetryElementKind = SymmetryElementKind::Proper(None);
1759
1760/// Improper symmetry element kind in the mirror-plane convention.
1761pub const SIG: SymmetryElementKind = SymmetryElementKind::ImproperMirrorPlane(None);
1762
1763/// Improper symmetry element kind in the inversion-centre convention.
1764pub const INV: SymmetryElementKind = SymmetryElementKind::ImproperInversionCentre(None);
1765
1766/// Time-reversed proper rotation symmetry element kind.
1767pub const TRROT: SymmetryElementKind = SymmetryElementKind::Proper(Some(TR));
1768
1769/// Time-reversed improper symmetry element kind in the mirror-plane convention.
1770pub const TRSIG: SymmetryElementKind = SymmetryElementKind::ImproperMirrorPlane(Some(TR));
1771
1772/// Time-reversed improper symmetry element kind in the inversion-centre convention.
1773pub const TRINV: SymmetryElementKind = SymmetryElementKind::ImproperInversionCentre(Some(TR));
1774
1775/// Rotation group $`\mathsf{SO}(3)`$.
1776pub const SO3: RotationGroup = RotationGroup::SO3;
1777
1778/// Rotation group $`\mathsf{SU}(2)`$, homotopy path of class 0.
1779pub const SU2_0: RotationGroup = RotationGroup::SU2(true);
1780
1781/// Rotation group $`\mathsf{SU}(2)`$, homotopy path of class 1.
1782pub const SU2_1: RotationGroup = RotationGroup::SU2(false);