qsym2/symmetry/symmetry_element/
symmetry_operation.rs

1//! Symmetry operations.
2
3use std::fmt;
4use std::hash::{Hash, Hasher};
5use std::ops::Mul;
6
7use anyhow::{self, format_err};
8use approx;
9use derive_builder::Builder;
10use fraction;
11use nalgebra::{Point3, Vector3};
12use ndarray::{Array2, Axis, ShapeBuilder};
13use num::Complex;
14use num_traits::{Inv, Pow, Zero};
15use ordered_float::OrderedFloat;
16use serde::{Deserialize, Serialize};
17
18use crate::angmom::spinor_rotation_3d::dmat_angleaxis_gen_single;
19use crate::auxiliary::geometry::{
20    self, improper_rotation_matrix, proper_rotation_matrix, PositiveHemisphere, Transform, IMINV,
21};
22use crate::auxiliary::misc::{self, HashableFloat};
23use crate::group::FiniteOrder;
24use crate::permutation::{IntoPermutation, PermutableCollection, Permutation};
25use crate::symmetry::symmetry_element::{
26    AntiunitaryKind, SymmetryElement, SymmetryElementKind, INV, ROT, SIG, SO3, SU2_0, SU2_1, TRINV,
27    TRROT, TRSIG,
28};
29use crate::symmetry::symmetry_element_order::ElementOrder;
30
31type F = fraction::GenericFraction<u32>;
32type Quaternion = (f64, Vector3<f64>);
33
34#[cfg(test)]
35#[path = "symmetry_operation_tests.rs"]
36mod symmetry_operation_tests;
37
38// =================
39// Trait definitions
40// =================
41
42/// Trait for special symmetry transformations.
43pub trait SpecialSymmetryTransformation {
44    // =================
45    // Group-theoretical
46    // =================
47
48    /// Checks if the proper rotation part of the symmetry operation is in $`\mathsf{SU}(2)`$.
49    ///
50    /// # Returns
51    ///
52    /// A boolean indicating if this symmetry operation contains an $`\mathsf{SU}(2)`$ proper
53    /// rotation.
54    fn is_su2(&self) -> bool;
55
56    /// Checks if the proper rotation part of the symmetry operation is in $`\mathsf{SU}(2)`$ and
57    /// connected to the identity via a homotopy path of class 1.
58    ///
59    /// # Returns
60    ///
61    /// A boolean indicating if this symmetry operation contains an $`\mathsf{SU}(2)`$ proper
62    /// rotation connected to the identity via a homotopy path of class 1.
63    fn is_su2_class_1(&self) -> bool;
64
65    // ============
66    // Spatial part
67    // ============
68
69    /// Checks if the spatial part of the symmetry operation is proper.
70    ///
71    /// # Returns
72    ///
73    /// A boolean indicating if the spatial part of the symmetry operation is proper.
74    fn is_proper(&self) -> bool;
75
76    /// Checks if the spatial part of the symmetry operation is the spatial identity.
77    ///
78    /// # Returns
79    ///
80    /// A boolean indicating if the spatial part of the symmetry operation is the spatial identity.
81    fn is_spatial_identity(&self) -> bool;
82
83    /// Checks if the spatial part of the symmetry operation is a spatial binary rotation.
84    ///
85    /// # Returns
86    ///
87    /// A boolean indicating if the spatial part of the symmetry operation is a spatial binary
88    /// rotation.
89    fn is_spatial_binary_rotation(&self) -> bool;
90
91    /// Checks if the spatial part of the symmetry operation is the spatial inversion.
92    ///
93    /// # Returns
94    ///
95    /// A boolean indicating if the spatial part of the symmetry operation is the spatial inversion.
96    fn is_spatial_inversion(&self) -> bool;
97
98    /// Checks if the spatial part of the symmetry operation is a spatial reflection.
99    ///
100    /// # Returns
101    ///
102    /// A boolean indicating if the spatial part of the symmetry operation is a spatial reflection.
103    fn is_spatial_reflection(&self) -> bool;
104
105    // ==================
106    // Time-reversal part
107    // ==================
108
109    /// Checks if the symmetry operation contains time reversal.
110    ///
111    /// # Returns
112    ///
113    /// A boolean indicating if the symmetry oppperation contains time reversal.
114    fn contains_time_reversal(&self) -> bool;
115
116    // ==========================
117    // Overall - provided methods
118    // ==========================
119
120    /// Checks if the symmetry operation is the identity in $`\mathsf{O}(3)`$, `E`, or
121    /// in $`\mathsf{SU}(2)`$, `E(Σ)`.
122    ///
123    /// # Returns
124    ///
125    /// A boolean indicating if this symmetry operation is the identity.
126    fn is_identity(&self) -> bool {
127        self.is_spatial_identity() && !self.contains_time_reversal() && !self.is_su2_class_1()
128    }
129
130    /// Checks if the symmetry operation is a pure time-reversal in $`\mathsf{O}(3)`$, `θ`, or
131    /// in $`\mathsf{SU}(2)`$, `θ(Σ)`.
132    ///
133    /// # Returns
134    ///
135    /// A boolean indicating if this symmetry operation is a pure time-reversal.
136    fn is_time_reversal(&self) -> bool {
137        self.is_spatial_identity() && self.contains_time_reversal() && !self.is_su2_class_1()
138    }
139
140    /// Checks if the symmetry operation is an inversion in $`\mathsf{O}(3)`$, `i`, but not in
141    /// $`\mathsf{SU}(2)`$, `i(Σ)`.
142    ///
143    /// # Returns
144    ///
145    /// A boolean indicating if this symmetry operation is an inversion in $`\mathsf{O}(3)`$.
146    fn is_inversion(&self) -> bool {
147        self.is_spatial_inversion() && !self.contains_time_reversal() && !self.is_su2()
148    }
149
150    /// Checks if the symmetry operation is a reflection in $`\mathsf{O}(3)`$, `σ`, but not in
151    /// $`\mathsf{SU}(2)`$, `σ(Σ)`.
152    ///
153    /// # Returns
154    ///
155    /// A boolean indicating if this symmetry operation is a reflection in $`\mathsf{O}(3)`$.
156    fn is_reflection(&self) -> bool {
157        self.is_spatial_reflection() && !self.contains_time_reversal() && !self.is_su2()
158    }
159}
160
161// ======================================
162// Struct definitions and implementations
163// ======================================
164
165/// Structure for managing symmetry operations generated from symmetry elements.
166///
167/// A symmetry element serves as a generator for symmetry operations. Thus, a symmetry element
168/// together with a signed integer indicating the number of times the symmetry element is applied
169/// (positively or negatively) specifies a symmetry operation.
170#[derive(Builder, Clone, Serialize, Deserialize)]
171pub struct SymmetryOperation {
172    /// The generating symmetry element for this symmetry operation.
173    pub generating_element: SymmetryElement,
174
175    /// The integral power indicating the number of times
176    /// [`Self::generating_element`] is applied to form the symmetry operation.
177    pub power: i32,
178
179    /// The total proper rotation angle associated with this operation (after taking into account
180    /// the power of the operation).
181    ///
182    /// This is simply the proper rotation angle of [`Self::generating_element`] multiplied by
183    /// [`Self::power`] and then folded onto the open interval $`(-\pi, \pi]`$.
184    ///
185    /// This angle lies in the open interval $`(-\pi, \pi]`$. For improper operations, this angle
186    /// depends on the convention used to describe the [`Self::generating_element`].
187    #[builder(setter(skip), default = "self.calc_total_proper_angle()")]
188    total_proper_angle: f64,
189
190    /// The fraction $`pk/n \in (-1/2, 1/2]`$ of the proper rotation, represented exactly for
191    /// hashing and comparison purposes.
192    ///
193    /// This is not defined for operations with infinite-order generating elements.
194    #[builder(setter(skip), default = "self.calc_total_proper_fraction()")]
195    pub(crate) total_proper_fraction: Option<F>,
196
197    /// The positive hemisphere used for distinguishing positive and negative rotation poles. If
198    /// `None`, the standard positive hemisphere as defined in S.L. Altmann, Rotations,
199    /// Quaternions, and Double Groups (Dover Publications, Inc., New York, 2005) is used.
200    #[builder(default = "None")]
201    pub positive_hemisphere: Option<PositiveHemisphere>,
202}
203
204impl SymmetryOperationBuilder {
205    fn calc_total_proper_angle(&self) -> f64 {
206        let (total_proper_angle, _) = geometry::normalise_rotation_angle(
207            self.generating_element
208                .as_ref()
209                .expect("Generating element has not been set.")
210                .proper_angle
211                .expect("Proper angle has not been set.")
212                * (f64::from(self.power.expect("Power has not been set."))),
213            self.generating_element
214                .as_ref()
215                .expect("Generating element has not been set.")
216                .threshold,
217        );
218        total_proper_angle
219    }
220
221    fn calc_total_proper_fraction(&self) -> Option<F> {
222        match self
223            .generating_element
224            .as_ref()
225            .expect("Generating element has not been set.")
226            .proper_fraction
227        {
228            Some(frac) => {
229                let pow = self.power.expect("Power has not been set.");
230                let (total_proper_fraction, _) =
231                    geometry::normalise_rotation_fraction(frac * F::from(pow));
232                Some(total_proper_fraction)
233            }
234            None => None,
235        }
236    }
237}
238
239impl SymmetryOperation {
240    /// Returns a builder to construct a new symmetry operation.
241    ///
242    /// # Returns
243    ///
244    /// A builder to construct a new symmetry operation.
245    #[must_use]
246    pub(crate) fn builder() -> SymmetryOperationBuilder {
247        SymmetryOperationBuilder::default()
248    }
249
250    /// Constructs a finite-order-element-generated symmetry operation from a quaternion.
251    ///
252    /// The rotation angle encoded in the quaternion is taken to be non-negative and assigned as
253    /// the proper rotation angle associated with the element generating the operation.
254    ///
255    /// If an improper operation is required, its generator will be constructed in the
256    /// inversion-centre convention.
257    ///
258    /// # Arguments
259    ///
260    /// * `qtn` - A quaternion encoding the proper rotation associated with the
261    /// generating element of the operation to be constructed.
262    /// * `proper` - A boolean indicating if the operation is proper or improper.
263    /// * `thresh` - Threshold for comparisons.
264    /// * `tr` - A boolean indicating if the resulting symmetry operation should be accompanied by
265    /// a time-reversal operation.
266    /// * `su2` - A boolean indicating if the resulting symmetry operation is to contain a proper
267    /// rotation in $`\mathsf{SU}(2)`$. The homotopy class of the operation will be deduced from
268    /// the specified quaternion.
269    /// * `poshem` - An option containing any custom positive hemisphere used to distinguish
270    /// positive and negative rotation poles.
271    ///
272    /// # Returns
273    ///
274    /// The constructed symmetry operation.
275    ///
276    /// # Panics
277    ///
278    /// Panics when the scalar part of the provided quaternion lies outside $`[0, 1]`$ by more than
279    /// the specified threshold `thresh`, or when the rotation angle associated with the quaternion
280    /// cannot be gracefully converted into an integer tuple of order and power.
281    #[must_use]
282    pub fn from_quaternion(
283        qtn: Quaternion,
284        proper: bool,
285        thresh: f64,
286        max_trial_power: u32,
287        tr: bool,
288        su2: bool,
289        poshem: Option<PositiveHemisphere>,
290    ) -> Self {
291        let (scalar_part, vector_part) = qtn;
292        let kind = if proper {
293            if tr {
294                TRROT
295            } else {
296                ROT
297            }
298        } else if tr {
299            TRINV
300        } else {
301            INV
302        };
303        let element = if su2 {
304            // SU(2)
305            assert!(
306                -1.0 - thresh <= scalar_part && scalar_part <= 1.0 + thresh,
307                "The scalar part of the quaternion must be in the interval [-1, +1]."
308            );
309            let (axis, order, power, su2_grp) =
310                if approx::abs_diff_eq!(scalar_part, 1.0, epsilon = thresh,) {
311                    // Zero-degree rotation, i.e. identity or inversion, class 0
312                    (Vector3::z(), 1u32, 1u32, SU2_0)
313                } else if approx::abs_diff_eq!(scalar_part, -1.0, epsilon = thresh,) {
314                    // 360-degree rotation, i.e. identity or inversion, class 1
315                    (Vector3::z(), 1u32, 1u32, SU2_1)
316                } else if approx::abs_diff_eq!(scalar_part, 0.0, epsilon = thresh,) {
317                    // 180-degree rotation, i.e. binary rotation or reflection. Whether the resultant
318                    // operation is in class 0 or class 1 depends on whether the vector part is in the
319                    // positive hemisphere or negative hemisphere.
320                    let positive_axis = poshem
321                        .as_ref()
322                        .cloned()
323                        .unwrap_or_default()
324                        .get_positive_pole(&vector_part, thresh);
325                    (
326                        positive_axis,
327                        2u32,
328                        1u32,
329                        if poshem
330                            .as_ref()
331                            .cloned()
332                            .unwrap_or_default()
333                            .check_positive_pole(&vector_part, thresh)
334                        {
335                            SU2_0
336                        } else {
337                            SU2_1
338                        },
339                    )
340                } else {
341                    // scalar_part != 0, 1, or -1
342                    let (standardised_scalar_part, standardised_vector_part, su2_grp) =
343                        if scalar_part > 0.0 {
344                            (scalar_part, vector_part, SU2_0)
345                        } else {
346                            (-scalar_part, -vector_part, SU2_1)
347                        };
348                    let half_proper_angle = standardised_scalar_part.acos();
349                    let proper_angle = 2.0 * half_proper_angle;
350                    let axis = standardised_vector_part / half_proper_angle.sin();
351                    let proper_fraction =
352                        geometry::get_proper_fraction(proper_angle, thresh, max_trial_power)
353                            .unwrap_or_else(|| {
354                                panic!(
355                                    "No proper fraction could be found for angle `{proper_angle}`."
356                                )
357                            });
358                    (
359                        axis,
360                        *proper_fraction.denom().unwrap_or_else(|| {
361                            panic!("Unable to extract the denominator of `{proper_fraction}`.")
362                        }),
363                        *proper_fraction.numer().unwrap_or_else(|| {
364                            panic!("Unable to extract the numerator of `{proper_fraction}`.")
365                        }),
366                        su2_grp,
367                    )
368                };
369            SymmetryElement::builder()
370                .threshold(thresh)
371                .proper_order(ElementOrder::Int(order))
372                .proper_power(
373                    power
374                        .try_into()
375                        .expect("Unable to convert the proper power to `i32`."),
376                )
377                .raw_axis(axis)
378                .kind(kind)
379                .rotation_group(su2_grp)
380                .build()
381                .unwrap_or_else(|_|
382                    panic!("Unable to construct a symmetry element of kind `{kind}` with the proper part in SU(2).")
383                )
384        } else {
385            // SO(3)
386            assert!(
387                -thresh <= scalar_part && scalar_part <= 1.0 + thresh,
388                "The scalar part of the quaternion must be in the interval [0, +1] when only SO(3) rotations are considered."
389            );
390            let (axis, order, power) = if approx::abs_diff_eq!(scalar_part, 1.0, epsilon = thresh,)
391            {
392                // Zero-degree rotation, i.e. identity or inversion
393                (Vector3::z(), 1u32, 1i32)
394            } else {
395                let half_proper_angle = scalar_part.acos(); // acos returns values in [0, π]
396                let proper_angle = 2.0 * half_proper_angle;
397                let axis = vector_part / half_proper_angle.sin();
398                let proper_fraction =
399                    geometry::get_proper_fraction(proper_angle, thresh, max_trial_power)
400                        .unwrap_or_else(|| {
401                            panic!("No proper fraction could be found for angle `{proper_angle}`.")
402                        });
403                let proper_power = if proper_fraction.is_sign_positive() {
404                    i32::try_from(*proper_fraction.numer().unwrap_or_else(|| {
405                        panic!("Unable to extract the numerator of `{proper_fraction}`.")
406                    }))
407                    .expect("Unable to convert the numerator of the proper fraction to `i32`.")
408                } else {
409                    -i32::try_from(*proper_fraction.numer().unwrap_or_else(|| {
410                        panic!("Unable to extract the numerator of `{proper_fraction}`.")
411                    }))
412                    .expect("Unable to convert the numerator of the proper fraction to `i32`.")
413                };
414                (
415                    axis,
416                    *proper_fraction.denom().unwrap_or_else(|| {
417                        panic!("Unable to extract the denominator of `{proper_fraction}`.")
418                    }),
419                    proper_power,
420                )
421            };
422
423            SymmetryElement::builder()
424                .threshold(thresh)
425                .proper_order(ElementOrder::Int(order))
426                .proper_power(power)
427                .raw_axis(axis)
428                .kind(kind)
429                .rotation_group(SO3)
430                .build()
431                .expect(
432                    "Unable to construct a symmetry element without an associated spin rotation.",
433                )
434        };
435
436        SymmetryOperation::builder()
437            .generating_element(element)
438            .power(1)
439            .positive_hemisphere(poshem)
440            .build()
441            .unwrap_or_else(|_|
442                panic!(
443                    "Unable to construct a symmetry operation of kind `{kind}` with {} rotation from the quaternion `{qtn:?}`.",
444                    if su2 { "SU(2)" } else { "SO(3)" }
445                )
446            )
447    }
448
449    /// Finds the quaternion associated with this operation.
450    ///
451    /// The rotation angle encoded in the quaternion is taken to be non-negative and assigned as
452    /// the proper rotation angle associated with the element generating the operation.
453    ///
454    /// If this is an operation generated from an improper element, the inversion-centre convention
455    /// will be used to determine the angle of proper rotation.
456    ///
457    /// Both $`\mathsf{SO}(3)`$ and $`\mathsf{SU}(2)`$ proper rotations are supported. For
458    /// $`\mathsf{SO}(3)`$ proper rotations, only quaternions in the standardised form are
459    /// returned.
460    ///
461    /// See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover Publications, Inc., New
462    /// York, 2005) (Chapter 9) for further information.
463    ///
464    /// # Returns
465    ///
466    /// The quaternion associated with this operation.
467    ///
468    /// # Panics
469    ///
470    /// Panics if the calculated scalar part of the quaternion lies outside the closed interval
471    /// $`[0, 1]`$ by more than the threshold value stored in the generating element in `self`.
472    #[must_use]
473    pub fn calc_quaternion(&self) -> Quaternion {
474        let c_self = if self.is_proper() {
475            self.clone()
476        } else {
477            // Time-reversal does not matter here.
478            self.convert_to_improper_kind(&INV)
479        };
480        debug_assert_eq!(
481            self.is_su2_class_1(),
482            c_self.is_su2_class_1(),
483            "`{self}` and `{c_self}` are in different homotopy classes."
484        );
485
486        // We only need the absolute value of the angle. Its sign information is
487        // encoded in the pole. `abs_angle` thus lies in [0, π], and so
488        //  cos(abs_angle/2) >= 0 and sin(abs_angle/2) >= 0.
489        // The scalar part is guaranteed to be in [0, 1].
490        // For binary rotations, the scalar part is zero, but the definition of pole ensures that
491        // the vector part still lies in the positive hemisphere.
492        let abs_angle = c_self.total_proper_angle.abs();
493        let scalar_part = (0.5 * abs_angle).cos();
494        let vector_part = (0.5 * abs_angle).sin() * c_self.calc_pole().coords;
495        debug_assert!(
496            -self.generating_element.threshold <= scalar_part
497                && scalar_part <= 1.0 + self.generating_element.threshold
498        );
499        debug_assert!(if approx::abs_diff_eq!(
500            scalar_part,
501            0.0,
502            epsilon = c_self.generating_element.threshold
503        ) {
504            c_self
505                .positive_hemisphere
506                .as_ref()
507                .cloned()
508                .unwrap_or_default()
509                .check_positive_pole(&vector_part, c_self.generating_element.threshold)
510        } else {
511            true
512        },);
513
514        if self.is_su2_class_1() {
515            (-scalar_part, -vector_part)
516        } else {
517            (scalar_part, vector_part)
518        }
519    }
520
521    /// Finds the pole associated with this operation with respect to the positive hemisphere
522    /// defined in [`Self::positive_hemisphere`].
523    ///
524    /// This is the point on the unit sphere that is left invariant by the operation.
525    ///
526    /// For improper operations, the inversion-centre convention is used to define
527    /// the pole. This allows a proper rotation and its improper partner to have the
528    /// same pole, thus facilitating the consistent specification of poles for the
529    /// identity / inversion and binary rotations / reflections.
530    ///
531    /// Note that binary rotations / reflections have unique poles on the positive
532    /// hemisphere (*i.e.*, $`C_2(\hat{\mathbf{n}}) = C_2^{-1}(\hat{\mathbf{n}})`$
533    /// and $`\sigma(\hat{\mathbf{n}}) = \sigma^{-1}(\hat{\mathbf{n}})`$).
534    ///
535    /// See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover
536    /// Publications, Inc., New York, 2005) (Chapter 9) for further information.
537    ///
538    /// # Returns
539    ///
540    /// The pole associated with this operation.
541    ///
542    /// # Panics
543    ///
544    /// Panics when no total proper fractions could be found for this operation.
545    #[must_use]
546    pub fn calc_pole(&self) -> Point3<f64> {
547        let op = if self.is_proper() {
548            self.clone()
549        } else {
550            // Time-reversal does not matter here.
551            self.convert_to_improper_kind(&INV)
552        };
553        match *op.generating_element.raw_proper_order() {
554            ElementOrder::Int(_) => {
555                let frac_1_2 = F::new(1u32, 2u32);
556                let total_proper_fraction = op
557                    .total_proper_fraction
558                    .expect("No total proper fractions found.");
559                if total_proper_fraction == frac_1_2 {
560                    // Binary rotations or reflections
561                    Point3::from(
562                        self.positive_hemisphere
563                            .as_ref()
564                            .cloned()
565                            .unwrap_or_default()
566                            .get_positive_pole(
567                                &op.generating_element.raw_axis,
568                                op.generating_element.threshold,
569                            ),
570                    )
571                } else if total_proper_fraction > F::zero() {
572                    // Positive rotation angles
573                    Point3::from(op.generating_element.raw_axis)
574                } else if total_proper_fraction < F::zero() {
575                    // Negative rotation angles
576                    Point3::from(-op.generating_element.raw_axis)
577                } else {
578                    // Identity or inversion
579                    assert!(total_proper_fraction.is_zero());
580                    Point3::from(Vector3::z())
581                }
582            }
583            ElementOrder::Inf => {
584                if approx::abs_diff_eq!(
585                    op.total_proper_angle,
586                    std::f64::consts::PI,
587                    epsilon = op.generating_element.threshold
588                ) {
589                    // Binary rotations or reflections
590                    Point3::from(
591                        self.positive_hemisphere
592                            .as_ref()
593                            .cloned()
594                            .unwrap_or_default()
595                            .get_positive_pole(
596                                &op.generating_element.raw_axis,
597                                op.generating_element.threshold,
598                            ),
599                    )
600                } else if approx::abs_diff_ne!(
601                    op.total_proper_angle,
602                    0.0,
603                    epsilon = op.generating_element.threshold
604                ) {
605                    Point3::from(op.total_proper_angle.signum() * op.generating_element.raw_axis)
606                } else {
607                    approx::assert_abs_diff_eq!(
608                        op.total_proper_angle,
609                        0.0,
610                        epsilon = op.generating_element.threshold
611                    );
612                    Point3::from(Vector3::z())
613                }
614            }
615        }
616    }
617
618    /// Finds the pole associated with the proper rotation of this operation.
619    ///
620    /// This is the point on the unit sphere that is left invariant by the proper rotation part of
621    /// the operation.
622    ///
623    /// For improper operations, no conversions will be performed, unlike in [`Self::calc_pole`].
624    ///
625    /// Note that binary rotations have unique poles on the positive hemisphere (*i.e.*,
626    /// $`C_2(\hat{\mathbf{n}}) = C_2^{-1}(\hat{\mathbf{n}})`$ and
627    /// $`\sigma(\hat{\mathbf{n}}) = \sigma^{-1}(\hat{\mathbf{n}})`$).
628    ///
629    /// See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover
630    /// Publications, Inc., New York, 2005) (Chapter 9) for further information.
631    ///
632    /// # Returns
633    ///
634    /// The pole associated with the proper rotation of this operation.
635    ///
636    /// # Panics
637    ///
638    /// Panics when no total proper fractions could be found for this operation.
639    #[must_use]
640    pub fn calc_proper_rotation_pole(&self) -> Point3<f64> {
641        match *self.generating_element.raw_proper_order() {
642            ElementOrder::Int(_) => {
643                let frac_1_2 = F::new(1u32, 2u32);
644                let total_proper_fraction = self
645                    .total_proper_fraction
646                    .expect("No total proper fractions found.");
647                if total_proper_fraction == frac_1_2 {
648                    // Binary rotations or reflections
649                    Point3::from(
650                        self.positive_hemisphere
651                            .as_ref()
652                            .cloned()
653                            .unwrap_or_default()
654                            .get_positive_pole(
655                                &self.generating_element.raw_axis,
656                                self.generating_element.threshold,
657                            ),
658                    )
659                } else if total_proper_fraction > F::zero() {
660                    // Positive rotation angles
661                    Point3::from(self.generating_element.raw_axis)
662                } else if total_proper_fraction < F::zero() {
663                    // Negative rotation angles
664                    Point3::from(-self.generating_element.raw_axis)
665                } else {
666                    // Identity or inversion
667                    assert!(total_proper_fraction.is_zero());
668                    Point3::from(Vector3::z())
669                }
670            }
671            ElementOrder::Inf => {
672                if approx::abs_diff_eq!(
673                    self.total_proper_angle,
674                    std::f64::consts::PI,
675                    epsilon = self.generating_element.threshold
676                ) {
677                    // Binary rotations or reflections
678                    Point3::from(
679                        self.positive_hemisphere
680                            .as_ref()
681                            .cloned()
682                            .unwrap_or_default()
683                            .get_positive_pole(
684                                &self.generating_element.raw_axis,
685                                self.generating_element.threshold,
686                            ),
687                    )
688                } else if approx::abs_diff_ne!(
689                    self.total_proper_angle,
690                    0.0,
691                    epsilon = self.generating_element.threshold
692                ) {
693                    Point3::from(
694                        self.total_proper_angle.signum() * self.generating_element.raw_axis,
695                    )
696                } else {
697                    approx::assert_abs_diff_eq!(
698                        self.total_proper_angle,
699                        0.0,
700                        epsilon = self.generating_element.threshold
701                    );
702                    Point3::from(Vector3::z())
703                }
704            }
705        }
706    }
707
708    /// Finds the pole angle associated with this operation.
709    ///
710    /// This is the non-negative angle that, together with the pole, uniquely determines the proper
711    /// part of this operation. This angle lies in the interval $`[0, \pi]`$.
712    ///
713    /// For improper operations, the inversion-centre convention is used to define the pole angle.
714    /// This allows a proper rotation and its improper partner to have the same pole angle, thus
715    /// facilitating the consistent specification of poles for the identity / inversion and binary
716    /// rotations / reflections.
717    ///
718    /// # Returns
719    ///
720    /// The pole angle associated with this operation.
721    ///
722    /// # Panics
723    ///
724    /// Panics when no total proper fractions could be found for this operation.
725    #[must_use]
726    pub fn calc_pole_angle(&self) -> f64 {
727        let c_self = if self.is_proper() {
728            self.clone()
729        } else {
730            // Time-reversal does not matter here.
731            self.convert_to_improper_kind(&INV)
732        };
733
734        c_self.total_proper_angle.abs()
735    }
736
737    /// Returns a copy of the current symmetry operation with the generating element
738    /// converted to the requested improper kind (power-preserving), provided that
739    /// it is an improper element.
740    ///
741    /// # Arguments
742    ///
743    /// * `improper_kind` - The improper kind to which `self` is to be converted. There is no need
744    /// to make sure the time reversal specification in `improper_kind` matches that of the
745    /// generating element of `self` as the conversion will take care of this.
746    ///
747    /// # Panics
748    ///
749    /// Panics if the converted symmetry operation cannot be constructed.
750    #[must_use]
751    pub fn convert_to_improper_kind(&self, improper_kind: &SymmetryElementKind) -> Self {
752        let c_element = self
753            .generating_element
754            .convert_to_improper_kind(improper_kind, true);
755        debug_assert_eq!(
756            self.generating_element.is_su2_class_1(),
757            c_element.is_su2_class_1()
758        );
759        Self::builder()
760            .generating_element(c_element)
761            .power(self.power)
762            .positive_hemisphere(self.positive_hemisphere.clone())
763            .build()
764            .expect("Unable to construct a symmetry operation.")
765    }
766
767    /// Converts the current symmetry operation $`O`$ to an equivalent symmetry element $`E`$ such
768    /// that $`O = E^1`$.
769    ///
770    /// The proper rotation axis of $`E`$ is the proper rotation pole (*not* the overall pole) of
771    /// $`O`$, and the proper rotation angle of $`E`$ is the total proper rotation angle of $`O`$,
772    /// either as an (order, power) integer tuple or an angle floating-point number.
773    ///
774    /// If $`O`$ is improper, then the improper generating element for $`E`$ is the same as that in
775    /// the generating element of $`O`$.
776    ///
777    /// # Returns
778    ///
779    /// The equivalent symmetry element $`E`$.
780    pub fn to_symmetry_element(&self) -> SymmetryElement {
781        let kind = if self.is_proper() {
782            let tr = self.contains_time_reversal();
783            if tr {
784                TRROT
785            } else {
786                ROT
787            }
788        } else {
789            self.generating_element.kind
790        };
791        let additional_superscript = if self.is_proper() {
792            String::new()
793        } else {
794            self.generating_element.additional_superscript.clone()
795        };
796        let additional_subscript = if self.is_proper() {
797            String::new()
798        } else {
799            self.generating_element.additional_subscript.clone()
800        };
801        let rotation_group = if self.is_su2_class_1() {
802            SU2_1
803        } else if self.is_su2() {
804            SU2_0
805        } else {
806            SO3
807        };
808        let axis = if self.is_spatial_reflection() {
809            self.positive_hemisphere
810                .as_ref()
811                .cloned()
812                .unwrap_or_default()
813                .get_positive_pole(
814                    self.generating_element.raw_axis(),
815                    self.generating_element.threshold,
816                )
817        } else {
818            self.calc_proper_rotation_pole().coords
819        };
820
821        if let Some(total_proper_fraction) = self.total_proper_fraction {
822            let proper_order = *total_proper_fraction
823                .denom()
824                .expect("Unable to extract the denominator of the total proper fraction.");
825            let numer = *total_proper_fraction
826                .numer()
827                .expect("Unable to extract the numerator of the total proper fraction.");
828            let proper_power =
829                i32::try_from(numer).expect("Unable to convert the numerator to `i32`.");
830            SymmetryElement::builder()
831                .threshold(self.generating_element.threshold())
832                .proper_order(ElementOrder::Int(proper_order))
833                .proper_power(proper_power)
834                .raw_axis(axis)
835                .kind(kind)
836                .rotation_group(rotation_group)
837                .additional_superscript(additional_superscript)
838                .additional_subscript(additional_subscript)
839                .build()
840                .unwrap()
841        } else {
842            let proper_angle = self.total_proper_angle;
843            SymmetryElement::builder()
844                .threshold(self.generating_element.threshold())
845                .proper_order(ElementOrder::Inf)
846                .proper_angle(proper_angle)
847                .raw_axis(axis)
848                .kind(kind)
849                .rotation_group(rotation_group)
850                .additional_superscript(additional_superscript)
851                .additional_subscript(additional_subscript)
852                .build()
853                .unwrap()
854        }
855    }
856
857    /// Generates the abbreviated symbol for this symmetry operation.
858    #[must_use]
859    pub fn get_abbreviated_symbol(&self) -> String {
860        self.to_symmetry_element()
861            .get_simplified_symbol_signed_power()
862    }
863
864    /// Returns the representation matrix for the spatial part of this symmetry operation.
865    ///
866    /// This representation matrix is in the basis of coordinate *functions* $`(y, z, x)`$.
867    #[must_use]
868    pub fn get_3d_spatial_matrix(&self) -> Array2<f64> {
869        if self.is_proper() {
870            if self.is_identity() || self.is_time_reversal() {
871                Array2::<f64>::eye(3)
872            } else {
873                let angle = self.calc_pole_angle();
874                let axis = self.calc_pole().coords;
875                let mat = proper_rotation_matrix(angle, &axis, 1);
876
877                // nalgebra matrix iter is column-major.
878                Array2::<f64>::from_shape_vec(
879                    (3, 3).f(),
880                    mat.iter().copied().collect::<Vec<_>>(),
881                )
882                .unwrap_or_else(
883                    |_| panic!(
884                        "Unable to construct a three-dimensional rotation matrix for angle {angle} and axis {axis}."
885                    )
886                )
887                .select(Axis(0), &[1, 2, 0])
888                .select(Axis(1), &[1, 2, 0])
889            }
890        } else if self.is_spatial_inversion() {
891            -Array2::<f64>::eye(3)
892        } else {
893            // Pole and pole angle are obtained in the inversion-centre convention.
894            let angle = self.calc_pole_angle();
895            let axis = self.calc_pole().coords;
896            let mat = improper_rotation_matrix(angle, &axis, 1, &IMINV);
897
898            // nalgebra matrix iter is column-major.
899            Array2::<f64>::from_shape_vec(
900                (3, 3).f(),
901                mat.iter().copied().collect::<Vec<_>>(),
902            )
903            .unwrap_or_else(
904                |_| panic!(
905                    "Unable to construct a three-dimensional improper rotation matrix for angle {angle} and axis {axis}."
906                )
907            )
908            .select(Axis(0), &[1, 2, 0])
909            .select(Axis(1), &[1, 2, 0])
910        }
911    }
912
913    /// Returns the representation matrix for this symmetry operation in a $`j`$-adapted basis.
914    ///
915    /// This representation matrix is in the basis of $`\ket{j, m_j}`$. For improper rotations in
916    /// half-odd-integer $`j`$ bases, the Pauli gauge is used where the inversion operation is
917    /// represented by a trivial identity matrix.
918    ///
919    /// # Arguments
920    ///
921    /// * `two_j` - The value of $`2j`$ for the basis.
922    /// * `increasingm` - A boolean indicating if the $`\ket{j, m_j}`$ kets are arranged in
923    /// increasing $`m_j`$ values.
924    ///
925    /// # Errors
926    ///
927    /// This function returns an error if `two_j` indicates a half-odd-integer spinor basis but the
928    /// symmetry operation is one in $`\mathsf{SO}(3)`$.
929    #[must_use]
930    pub fn get_wigner_matrix(
931        &self,
932        two_j: u32,
933        increasingm: bool,
934    ) -> Result<Array2<Complex<f64>>, anyhow::Error> {
935        let spinor_basis = two_j.rem_euclid(2) == 1;
936        if spinor_basis && !self.is_su2() {
937            Err(format_err!(
938                "Unable to construct a Wigner matrix for an SO(3) rotation in a spinor basis."
939            ))
940        } else {
941            let odd_sh_basis = (!spinor_basis) && two_j.rem_euclid(4) == 2;
942            let dmat_rotation = {
943                let angle = self.calc_pole_angle();
944                let axis = self.calc_pole().coords;
945                if (spinor_basis && self.is_su2_class_1()) || (odd_sh_basis && !self.is_proper()) {
946                    -dmat_angleaxis_gen_single(two_j, angle, axis, increasingm)
947                } else {
948                    dmat_angleaxis_gen_single(two_j, angle, axis, increasingm)
949                }
950            };
951            if self.contains_time_reversal() {
952                // The complex conjugation is required by corepresentation theory.
953                Ok(dmat_angleaxis_gen_single(
954                    two_j,
955                    std::f64::consts::PI,
956                    Vector3::y(),
957                    increasingm,
958                )
959                .dot(&dmat_rotation.map(|x| x.conj())))
960            } else {
961                Ok(dmat_rotation)
962            }
963        }
964    }
965
966    /// Convert the proper rotation of the current operation to one in hopotopy class 0 of
967    /// $`\mathsf{SU}(2)`$.
968    ///
969    /// # Returns
970    ///
971    /// A symmetry element in $`\mathsf{SU}(2)`$.
972    pub fn to_su2_class_0(&self) -> Self {
973        let q_identity = Self::from_quaternion(
974            (-1.0, -Vector3::z()),
975            true,
976            self.generating_element.threshold(),
977            1,
978            false,
979            true,
980            None,
981        );
982        if self.is_su2() {
983            if self.is_su2_class_1() {
984                self * q_identity
985            } else {
986                self.clone()
987            }
988        } else {
989            let mut op = self.clone();
990            op.generating_element.rotation_group = SU2_0;
991            if op.is_su2_class_1() {
992                let mut q_op = op * q_identity;
993                if !q_op.is_proper() {
994                    q_op = q_op.convert_to_improper_kind(&SIG);
995                }
996                q_op
997            } else {
998                op
999            }
1000        }
1001    }
1002
1003    /// Sets the positive hemisphere governing this symmetry operation.
1004    ///
1005    /// # Arguments
1006    ///
1007    /// * `poshem` - An `Option` containing a custom positive hemisphere, if any.
1008    pub fn set_positive_hemisphere(&mut self, poshem: Option<&PositiveHemisphere>) {
1009        self.positive_hemisphere = poshem.cloned();
1010    }
1011}
1012
1013// =====================
1014// Trait implementations
1015// =====================
1016
1017impl FiniteOrder for SymmetryOperation {
1018    type Int = u32;
1019
1020    /// Calculates the order of this symmetry operation.
1021    fn order(&self) -> Self::Int {
1022        let denom = *self
1023            .total_proper_fraction
1024            .expect("No total proper fractions found.")
1025            .denom()
1026            .expect("Unable to extract the denominator.");
1027        let spatial_order =
1028            if (self.is_proper() && !self.contains_time_reversal()) || denom.rem_euclid(2) == 0 {
1029                denom
1030            } else {
1031                2 * denom
1032            };
1033        if self.is_su2() {
1034            2 * spatial_order
1035        } else {
1036            spatial_order
1037        }
1038    }
1039}
1040
1041impl SpecialSymmetryTransformation for SymmetryOperation {
1042    // ============
1043    // Spatial part
1044    // ============
1045
1046    /// Checks if the spatial part of the symmetry operation is proper.
1047    ///
1048    /// # Returns
1049    ///
1050    /// A boolean indicating if the spatial part of the symmetry operation is proper.
1051    fn is_proper(&self) -> bool {
1052        let au = self.generating_element.contains_antiunitary();
1053        self.generating_element.is_o3_proper(au) || self.power.rem_euclid(2) == 0
1054    }
1055
1056    /// Checks if the spatial part of the symmetry operation is the spatial identity.
1057    ///
1058    /// # Returns
1059    ///
1060    /// A boolean indicating if the spatial part of the symmetry operation is the spatial identity.
1061    fn is_spatial_identity(&self) -> bool {
1062        self.is_proper()
1063            && match *self.generating_element.raw_proper_order() {
1064                ElementOrder::Int(_) => self
1065                    .total_proper_fraction
1066                    .expect("Total proper fraction not found for a finite-order operation.")
1067                    .is_zero(),
1068                ElementOrder::Inf => approx::abs_diff_eq!(
1069                    self.total_proper_angle,
1070                    0.0,
1071                    epsilon = self.generating_element.threshold
1072                ),
1073            }
1074    }
1075
1076    /// Checks if the spatial part of the symmetry operation is a spatial binary rotation.
1077    ///
1078    /// # Returns
1079    ///
1080    /// A boolean indicating if the spatial part of the symmetry operation is a spatial binary
1081    /// rotation.
1082    fn is_spatial_binary_rotation(&self) -> bool {
1083        self.is_proper()
1084            && match *self.generating_element.raw_proper_order() {
1085                ElementOrder::Int(_) => {
1086                    self.total_proper_fraction
1087                        .expect("Total proper fraction not found for a finite-order operation.")
1088                        == F::new(1u32, 2u32)
1089                }
1090                ElementOrder::Inf => {
1091                    approx::abs_diff_eq!(
1092                        self.total_proper_angle,
1093                        std::f64::consts::PI,
1094                        epsilon = self.generating_element.threshold
1095                    )
1096                }
1097            }
1098    }
1099
1100    /// Checks if the spatial part of the symmetry operation is the spatial inversion.
1101    ///
1102    /// # Returns
1103    ///
1104    /// A boolean indicating if the spatial part of the symmetry operation is the spatial inversion.
1105    fn is_spatial_inversion(&self) -> bool {
1106        !self.is_proper()
1107            && match self.generating_element.kind {
1108                SymmetryElementKind::ImproperMirrorPlane(_) => {
1109                    if let ElementOrder::Int(_) = *self.generating_element.raw_proper_order() {
1110                        self.total_proper_fraction
1111                            .expect("Total proper fraction not found for a finite-order operation.")
1112                            == F::new(1u32, 2u32)
1113                    } else {
1114                        approx::abs_diff_eq!(
1115                            self.total_proper_angle,
1116                            std::f64::consts::PI,
1117                            epsilon = self.generating_element.threshold
1118                        )
1119                    }
1120                }
1121                SymmetryElementKind::ImproperInversionCentre(_) => {
1122                    if let ElementOrder::Int(_) = *self.generating_element.raw_proper_order() {
1123                        self.total_proper_fraction
1124                            .expect("Total proper fraction not found for a finite-order operation.")
1125                            .is_zero()
1126                    } else {
1127                        approx::abs_diff_eq!(
1128                            self.total_proper_angle,
1129                            0.0,
1130                            epsilon = self.generating_element.threshold
1131                        )
1132                    }
1133                }
1134                _ => false,
1135            }
1136    }
1137
1138    /// Checks if the spatial part of the symmetry operation is a spatial reflection.
1139    ///
1140    /// # Returns
1141    ///
1142    /// A boolean indicating if the spatial part of the symmetry operation is a spatial reflection.
1143    fn is_spatial_reflection(&self) -> bool {
1144        !self.is_proper()
1145            && match self.generating_element.kind {
1146                SymmetryElementKind::ImproperMirrorPlane(_) => {
1147                    if let ElementOrder::Int(_) = *self.generating_element.raw_proper_order() {
1148                        self.total_proper_fraction
1149                            .expect("Total proper fraction not found for a finite-order operation.")
1150                            .is_zero()
1151                    } else {
1152                        approx::abs_diff_eq!(
1153                            self.total_proper_angle,
1154                            0.0,
1155                            epsilon = self.generating_element.threshold
1156                        )
1157                    }
1158                }
1159                SymmetryElementKind::ImproperInversionCentre(_) => {
1160                    if let ElementOrder::Int(_) = self.generating_element.raw_proper_order() {
1161                        self.total_proper_fraction
1162                            .expect("Total proper fraction not found for a finite-order operation.")
1163                            == F::new(1u32, 2u32)
1164                    } else {
1165                        approx::abs_diff_eq!(
1166                            self.total_proper_angle,
1167                            std::f64::consts::PI,
1168                            epsilon = self.generating_element.threshold
1169                        )
1170                    }
1171                }
1172                _ => false,
1173            }
1174    }
1175
1176    // ==================
1177    // Time-reversal part
1178    // ==================
1179
1180    /// Checks if the symmetry operation is antiunitary or not.
1181    ///
1182    /// # Returns
1183    ///
1184    /// A boolean indicating if the symmetry oppperation is antiunitary.
1185    fn contains_time_reversal(&self) -> bool {
1186        self.generating_element.contains_time_reversal() && self.power.rem_euclid(2) == 1
1187    }
1188
1189    // ==================
1190    // Spin rotation part
1191    // ==================
1192
1193    /// Checks if the proper rotation part of the symmetry operation is in $`\mathsf{SU}(2)`$.
1194    ///
1195    /// # Returns
1196    ///
1197    /// A boolean indicating if this symmetry operation contains an $`\mathsf{SU}(2)`$ proper
1198    /// rotation.
1199    fn is_su2(&self) -> bool {
1200        self.generating_element.rotation_group.is_su2()
1201    }
1202
1203    /// Checks if the proper rotation part of the symmetry operation (in the inversion convention)
1204    /// is in $`\mathsf{SU}(2)`$ and connected to the identity via a homotopy path of class 1.
1205    ///
1206    /// # Returns
1207    ///
1208    /// A boolean indicating if this symmetry operation contains an $`\mathsf{SU}(2)`$ proper
1209    /// rotation connected to the identity via a homotopy path of class 1.
1210    fn is_su2_class_1(&self) -> bool {
1211        if self.is_su2() {
1212            // The following is wrong, because `self.is_proper()` takes into account the power applied
1213            // to the spatial part, but not yet to the spin rotation part. Then, for example,
1214            // [QΣ·S3(+0.816, -0.408, +0.408)]^2 would become Σ'·[C3(+0.816, -0.408, +0.408)]^2 where
1215            // Σ' is the associated spin rotation of [C3(+0.816, -0.408, +0.408)]^2, which is not the
1216            // same as Σ^2.
1217            // let c_self = if self.is_proper() {
1218            //     self.clone()
1219            // } else {
1220            //     self.convert_to_improper_kind(&INV)
1221            // };
1222            //
1223            // The following is correct.
1224            let c_self = match self.generating_element.kind {
1225                SymmetryElementKind::Proper(_)
1226                | SymmetryElementKind::ImproperInversionCentre(_) => self.clone(),
1227                SymmetryElementKind::ImproperMirrorPlane(au) => {
1228                    self.convert_to_improper_kind(&SymmetryElementKind::ImproperInversionCentre(au))
1229                }
1230            };
1231            let generating_element_au = c_self.generating_element.contains_antiunitary();
1232            let spatial_proper_identity = c_self
1233                .generating_element
1234                .is_o3_identity(generating_element_au)
1235                || c_self
1236                    .generating_element
1237                    .is_o3_inversion_centre(generating_element_au);
1238
1239            let inverse_from_time_reversal =
1240                if self.is_su2() && generating_element_au == Some(AntiunitaryKind::TimeReversal) {
1241                    self.power.rem_euclid(4) == 2 || self.power.rem_euclid(4) == 3
1242                } else {
1243                    false
1244                };
1245
1246            let inverse_from_rotation_group = if spatial_proper_identity {
1247                // The proper part of the generating element is the identity. In this case, no
1248                // matter the value of proper power, the result is always the identity.
1249                false
1250            } else {
1251                let thresh = c_self.generating_element.threshold;
1252                let odd_jumps_from_angle = c_self
1253                    .generating_element
1254                    .proper_fraction
1255                    .map(|frac| {
1256                        let pow = c_self.power;
1257                        let total_unormalised_proper_fraction = frac * F::from(pow);
1258                        let (_, x) = geometry::normalise_rotation_fraction(
1259                            total_unormalised_proper_fraction,
1260                        );
1261                        x.rem_euclid(2) == 1
1262                    })
1263                    .unwrap_or_else(|| {
1264                        let total_unormalised_proper_angle = c_self
1265                            .generating_element
1266                            .proper_angle
1267                            .expect("Proper angle of generating element not found.")
1268                            * f64::from(c_self.power);
1269                        let (_, x) = geometry::normalise_rotation_angle(
1270                            total_unormalised_proper_angle,
1271                            thresh,
1272                        );
1273                        x.rem_euclid(2) == 1
1274                    });
1275                let single_jump_from_c2 = (c_self.is_spatial_binary_rotation()
1276                    || c_self.is_spatial_reflection())
1277                    && !self
1278                        .positive_hemisphere
1279                        .as_ref()
1280                        .cloned()
1281                        .unwrap_or_default()
1282                        .check_positive_pole(c_self.generating_element.raw_axis(), thresh);
1283                odd_jumps_from_angle != single_jump_from_c2
1284            };
1285
1286            let intrinsic_inverse = c_self.generating_element.rotation_group().is_su2_class_1()
1287                && c_self.power.rem_euclid(2) == 1;
1288
1289            let inverse_count = [
1290                inverse_from_time_reversal,
1291                inverse_from_rotation_group,
1292                intrinsic_inverse,
1293            ]
1294            .into_iter()
1295            .filter(|&inverse| inverse)
1296            .count();
1297
1298            inverse_count.rem_euclid(2) == 1
1299        } else {
1300            false
1301        }
1302    }
1303}
1304
1305impl fmt::Debug for SymmetryOperation {
1306    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1307        if self.power == 1 {
1308            write!(f, "{:?}", self.generating_element)
1309        } else if self.power >= 0 {
1310            write!(f, "[{:?}]^{}", self.generating_element, self.power)
1311        } else {
1312            write!(f, "[{:?}]^({})", self.generating_element, self.power)
1313        }
1314    }
1315}
1316
1317impl fmt::Display for SymmetryOperation {
1318    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1319        if self.power == 1 {
1320            write!(f, "{}", self.generating_element)
1321        } else if self.power >= 0 {
1322            write!(f, "[{}]^{}", self.generating_element, self.power)
1323        } else {
1324            write!(f, "[{}]^({})", self.generating_element, self.power)
1325        }
1326    }
1327}
1328
1329impl PartialEq for SymmetryOperation {
1330    fn eq(&self, other: &Self) -> bool {
1331        if (*self.generating_element.raw_proper_order() == ElementOrder::Inf)
1332            != (*other.generating_element.raw_proper_order() == ElementOrder::Inf)
1333        {
1334            // We disable comparisons between operations with infinite-order and
1335            // finite-order generating elements, because they cannot be made to
1336            // have the same hashes without losing the fidelity of exact-fraction
1337            // representations for operations with finite-order generating elements.
1338            return false;
1339        }
1340
1341        // =================
1342        // Group-theoretical
1343        // =================
1344
1345        if self.is_su2() != other.is_su2() {
1346            return false;
1347        }
1348
1349        if self.is_su2_class_1() != other.is_su2_class_1() {
1350            return false;
1351        }
1352
1353        // ==========================
1354        // Special general operations
1355        // ==========================
1356
1357        if self.is_proper() != other.is_proper() {
1358            return false;
1359        }
1360
1361        if self.contains_time_reversal() != other.contains_time_reversal() {
1362            return false;
1363        }
1364
1365        // ===========================
1366        // Special specific operations
1367        // ===========================
1368
1369        // At this stage, `self` and `other` must have the same spatial parity, unitarity, and
1370        // SO3/SU2 properties.
1371        if self.is_spatial_identity() && other.is_spatial_identity() {
1372            // assert_eq!(misc::calculate_hash(self), misc::calculate_hash(other));
1373            // return true;
1374            return misc::calculate_hash(self) == misc::calculate_hash(other);
1375        }
1376
1377        // ======
1378        // Others
1379        // ======
1380
1381        let thresh =
1382            (self.generating_element.threshold * other.generating_element.threshold).sqrt();
1383
1384        let result = if (self.is_spatial_binary_rotation() && other.is_spatial_binary_rotation())
1385            || (self.is_spatial_reflection() && other.is_spatial_reflection())
1386        {
1387            approx::relative_eq!(
1388                self.calc_pole(),
1389                other.calc_pole(),
1390                epsilon = thresh,
1391                max_relative = thresh
1392            )
1393        } else {
1394            let c_self = if self.is_proper() {
1395                self.clone()
1396            } else {
1397                // Time-reversal does not matter here.
1398                self.convert_to_improper_kind(&INV)
1399            };
1400            let c_other = if other.is_proper() {
1401                other.clone()
1402            } else {
1403                // Time-reversal does not matter here.
1404                other.convert_to_improper_kind(&INV)
1405            };
1406
1407            let angle_comparison = if let (Some(s_frac), Some(o_frac)) =
1408                (c_self.total_proper_fraction, c_other.total_proper_fraction)
1409            {
1410                s_frac.abs() == o_frac.abs()
1411            } else {
1412                approx::relative_eq!(
1413                    c_self.total_proper_angle.abs(),
1414                    c_other.total_proper_angle.abs(),
1415                    epsilon = thresh,
1416                    max_relative = thresh
1417                )
1418            };
1419
1420            angle_comparison
1421                && approx::relative_eq!(
1422                    self.calc_pole(),
1423                    other.calc_pole(),
1424                    epsilon = thresh,
1425                    max_relative = thresh
1426                )
1427        };
1428
1429        // if result {
1430        //     assert_eq!(
1431        //         misc::calculate_hash(self),
1432        //         misc::calculate_hash(other),
1433        //         "`{self}` and `{other}` have unequal hashes.",
1434        //     );
1435        // }
1436        result && (misc::calculate_hash(self) == misc::calculate_hash(other))
1437    }
1438}
1439
1440impl Eq for SymmetryOperation {}
1441
1442impl Hash for SymmetryOperation {
1443    fn hash<H: Hasher>(&self, state: &mut H) {
1444        let c_self = match self.generating_element.kind {
1445            SymmetryElementKind::Proper(_) | SymmetryElementKind::ImproperInversionCentre(_) => {
1446                self.clone()
1447            }
1448            SymmetryElementKind::ImproperMirrorPlane(_) => self.convert_to_improper_kind(&INV),
1449        };
1450        // ==========================
1451        // Special general operations
1452        // ==========================
1453        c_self.is_proper().hash(state);
1454        c_self.contains_time_reversal().hash(state);
1455        c_self.is_su2().hash(state);
1456        c_self.is_su2_class_1().hash(state);
1457
1458        // ===========================
1459        // Special specific operations
1460        // ===========================
1461        if c_self.is_spatial_identity() {
1462            true.hash(state);
1463        } else {
1464            let pole = c_self.calc_pole();
1465            pole[0]
1466                .round_factor(c_self.generating_element.threshold)
1467                .integer_decode()
1468                .hash(state);
1469            pole[1]
1470                .round_factor(c_self.generating_element.threshold)
1471                .integer_decode()
1472                .hash(state);
1473            pole[2]
1474                .round_factor(c_self.generating_element.threshold)
1475                .integer_decode()
1476                .hash(state);
1477
1478            if !c_self.is_spatial_binary_rotation() && !c_self.is_spatial_reflection() {
1479                if let Some(frac) = c_self.total_proper_fraction {
1480                    // self.total_proper_fraction lies in (-1/2, 0) ∪ (0, 1/2).
1481                    // 0 and 1/2 are excluded because this is not an identity,
1482                    // inversion, binary rotation, or reflection.
1483                    frac.abs().hash(state);
1484                } else {
1485                    // self.total_proper_angle lies in (-π, 0) ∪ (0, π).
1486                    // 0 and π are excluded because this is not an identity,
1487                    // inversion, binary rotation, or reflection.
1488                    let abs_ang = c_self.total_proper_angle.abs();
1489                    abs_ang
1490                        .round_factor(c_self.generating_element.threshold)
1491                        .integer_decode()
1492                        .hash(state);
1493                };
1494            }
1495        };
1496    }
1497}
1498
1499impl Mul<&'_ SymmetryOperation> for &SymmetryOperation {
1500    type Output = SymmetryOperation;
1501
1502    fn mul(self, rhs: &SymmetryOperation) -> Self::Output {
1503        assert_eq!(
1504            self.is_su2(),
1505            rhs.is_su2(),
1506            "`self` and `rhs` must both have or not have associated spin rotations."
1507        );
1508        assert_eq!(self.positive_hemisphere, rhs.positive_hemisphere);
1509        let su2 = self.is_su2();
1510        let (q1_s, q1_v) = self.calc_quaternion();
1511        let (q2_s, q2_v) = rhs.calc_quaternion();
1512
1513        let q3_s = q1_s * q2_s - q1_v.dot(&q2_v);
1514        let q3_v = q1_s * q2_v + q2_s * q1_v + q1_v.cross(&q2_v);
1515
1516        // Is the resulting operation proper?
1517        let proper = self.is_proper() == rhs.is_proper();
1518
1519        // Does the resulting operation contain a time reversal?
1520        let tr = self.contains_time_reversal() != rhs.contains_time_reversal();
1521
1522        // Does the resulting operation pick up a quaternion sign change due to θ^2?
1523        let tr2 = self.contains_time_reversal() && rhs.contains_time_reversal();
1524
1525        let thresh = (self.generating_element.threshold * rhs.generating_element.threshold).sqrt();
1526        let max_trial_power = u32::MAX;
1527
1528        let q3 = if su2 {
1529            if tr2 {
1530                (-q3_s, -q3_v)
1531            } else {
1532                (q3_s, q3_v)
1533            }
1534        } else if q3_s >= 0.0 {
1535            (q3_s, q3_v)
1536        } else {
1537            (-q3_s, -q3_v)
1538        };
1539
1540        SymmetryOperation::from_quaternion(
1541            q3,
1542            proper,
1543            thresh,
1544            max_trial_power,
1545            tr,
1546            su2,
1547            self.positive_hemisphere.clone(),
1548        )
1549    }
1550}
1551
1552impl Mul<&'_ SymmetryOperation> for SymmetryOperation {
1553    type Output = SymmetryOperation;
1554
1555    fn mul(self, rhs: &SymmetryOperation) -> Self::Output {
1556        &self * rhs
1557    }
1558}
1559
1560impl Mul<SymmetryOperation> for SymmetryOperation {
1561    type Output = SymmetryOperation;
1562
1563    fn mul(self, rhs: SymmetryOperation) -> Self::Output {
1564        &self * &rhs
1565    }
1566}
1567
1568impl Mul<SymmetryOperation> for &SymmetryOperation {
1569    type Output = SymmetryOperation;
1570
1571    fn mul(self, rhs: SymmetryOperation) -> Self::Output {
1572        self * &rhs
1573    }
1574}
1575
1576impl Pow<i32> for &SymmetryOperation {
1577    type Output = SymmetryOperation;
1578
1579    fn pow(self, rhs: i32) -> Self::Output {
1580        SymmetryOperation::builder()
1581            .generating_element(self.generating_element.clone())
1582            .power(self.power * rhs)
1583            .positive_hemisphere(self.positive_hemisphere.clone())
1584            .build()
1585            .expect("Unable to construct a symmetry operation.")
1586    }
1587}
1588
1589impl Pow<i32> for SymmetryOperation {
1590    type Output = SymmetryOperation;
1591
1592    fn pow(self, rhs: i32) -> Self::Output {
1593        (&self).pow(rhs)
1594    }
1595}
1596
1597impl Inv for &SymmetryOperation {
1598    type Output = SymmetryOperation;
1599
1600    fn inv(self) -> Self::Output {
1601        SymmetryOperation::builder()
1602            .generating_element(self.generating_element.clone())
1603            .power(-self.power)
1604            .positive_hemisphere(self.positive_hemisphere.clone())
1605            .build()
1606            .expect("Unable to construct an inverse symmetry operation.")
1607    }
1608}
1609
1610impl Inv for SymmetryOperation {
1611    type Output = SymmetryOperation;
1612
1613    fn inv(self) -> Self::Output {
1614        (&self).inv()
1615    }
1616}
1617
1618impl<M> IntoPermutation<M> for SymmetryOperation
1619where
1620    M: Transform + PermutableCollection<Rank = usize>,
1621{
1622    fn act_permute(&self, rhs: &M) -> Option<Permutation<usize>> {
1623        let angle = self.calc_pole_angle();
1624        let axis = self.calc_pole().coords;
1625        let mut t_mol = if self.is_proper() {
1626            rhs.rotate(angle, &axis)
1627        } else {
1628            rhs.improper_rotate(angle, &axis, &IMINV)
1629        };
1630        if self.contains_time_reversal() {
1631            t_mol.reverse_time_mut();
1632        }
1633        rhs.get_perm_of(&t_mol)
1634    }
1635}
1636
1637// =================
1638// Utility functions
1639// =================
1640/// Sorts symmetry operations in-place based on:
1641///
1642/// * whether they are unitary or antiunitary
1643/// * whether they are proper or improper
1644/// * whether they are the identity or inversion
1645/// * whether they are a spatial binary rotation or spatial reflection
1646/// * their orders
1647/// * their powers
1648/// * their closeness to Cartesian axes
1649/// * the axes of closest inclination
1650/// * whether they are of homotopy class 1 in $`\mathsf{SU}'(2)`$.
1651///
1652/// # Arguments
1653///
1654/// * `operations` - A mutable reference to a vector of symmetry operations.
1655pub(crate) fn sort_operations(operations: &mut [SymmetryOperation]) {
1656    operations.sort_by_key(|op| {
1657        let (axis_closeness, closest_axis) = op.generating_element.closeness_to_cartesian_axes();
1658        let c_op = if op.is_proper()
1659            || op.generating_element.kind == SIG
1660            || op.generating_element.kind == TRSIG
1661        {
1662            op.clone()
1663        } else if op.contains_time_reversal() {
1664            op.convert_to_improper_kind(&TRSIG)
1665        } else {
1666            op.convert_to_improper_kind(&SIG)
1667        };
1668
1669        let total_proper_fraction = c_op
1670            .total_proper_fraction
1671            .expect("No total proper fractions found.");
1672        let denom = i64::try_from(
1673            *total_proper_fraction
1674                .denom()
1675                .expect("The denominator of the total proper fraction cannot be extracted."),
1676        )
1677        .unwrap_or_else(|_| {
1678            panic!("Unable to convert the denominator of `{total_proper_fraction:?}` to `i64`.")
1679        });
1680        let numer = i64::try_from(
1681            *total_proper_fraction
1682                .numer()
1683                .expect("The numerator of the total proper fraction cannot be extracted."),
1684        )
1685        .unwrap_or_else(|_| {
1686            panic!("Unable to convert the numerator of `{total_proper_fraction:?}` to `i64`.")
1687        });
1688
1689        let negative_rotation = !c_op
1690            .positive_hemisphere
1691            .as_ref()
1692            .cloned()
1693            .unwrap_or_default()
1694            .check_positive_pole(
1695                &c_op.calc_proper_rotation_pole().coords,
1696                c_op.generating_element.threshold(),
1697            );
1698
1699        (
1700            c_op.contains_time_reversal(),
1701            !c_op.is_proper(),
1702            !(c_op.is_spatial_identity() || c_op.is_spatial_inversion()),
1703            c_op.is_spatial_binary_rotation() || c_op.is_spatial_reflection(),
1704            -denom,
1705            negative_rotation,
1706            if negative_rotation { -numer } else { numer },
1707            OrderedFloat(axis_closeness),
1708            closest_axis,
1709            c_op.is_su2_class_1(),
1710        )
1711    });
1712}