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}