qsym2/symmetry/
symmetry_symbols.rs

1//! Symbols for enumerating rows and columns of character tables of symmetry groups.
2
3use std::cmp::Ordering;
4use std::collections::HashSet;
5use std::fmt;
6use std::hash::{Hash, Hasher};
7use std::str::FromStr;
8
9use approx;
10use derive_builder::Builder;
11use indexmap::{IndexMap, IndexSet};
12use itertools::Itertools;
13use log;
14use nalgebra::Vector3;
15use ndarray::{Array2, ArrayView2, Axis};
16use num_traits::ToPrimitive;
17use phf::{phf_map, phf_set};
18use serde::{Deserialize, Serialize};
19
20use crate::chartab::CharacterTable;
21use crate::chartab::character::Character;
22use crate::chartab::chartab_group::CharacterProperties;
23use crate::chartab::chartab_symbols::{
24    CollectionSymbol, DecomposedSymbol, DecomposedSymbolBuilderError, GenericSymbol,
25    GenericSymbolParsingError, LinearSpaceSymbol, MathematicalSymbol, ReducibleLinearSpaceSymbol,
26    disambiguate_linspace_symbols,
27};
28use crate::chartab::unityroot::UnityRoot;
29use crate::group::FiniteOrder;
30use crate::symmetry::symmetry_element::SymmetryElement;
31use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
32use crate::symmetry::symmetry_element_order::ORDER_1;
33use crate::symmetry::symmetry_group::SymmetryGroupProperties;
34
35#[cfg(test)]
36#[path = "symmetry_symbols_tests.rs"]
37mod symmetry_symbols_tests;
38
39// =========
40// Constants
41// =========
42
43/// Map from Mulliken degeneracy labels to degeneracy degrees.
44static MULLIKEN_IRREP_DEGENERACIES: phf::Map<&'static str, u64> = phf_map! {
45    "A" => 1u64,
46    "B" => 1u64,
47    "Σ" => 1u64,
48    "Γ" => 1u64,
49    "E" => 2u64,
50    "Π" => 2u64,
51    "Δ" => 2u64,
52    "Φ" => 2u64,
53    "T" => 3u64,
54    "F" => 4u64,
55    "H" => 5u64,
56    "I" => 6u64,
57    "J" => 7u64,
58    "K" => 8u64,
59    "L" => 9u64,
60    "M" => 10u64,
61    "A~" => 1u64,
62    "B~" => 1u64,
63    "Σ~" => 1u64,
64    "Γ~" => 1u64,
65    "E~" => 2u64,
66    "Π~" => 2u64,
67    "Δ~" => 2u64,
68    "Φ~" => 2u64,
69    "T~" => 3u64,
70    "F~" => 4u64,
71    "H~" => 5u64,
72    "I~" => 6u64,
73    "J~" => 7u64,
74    "K~" => 8u64,
75    "L~" => 9u64,
76    "M~" => 10u64,
77};
78
79/// Map from degeneracy degrees to Mulliken degeneracy labels.
80static INV_MULLIKEN_IRREP_DEGENERACIES: phf::Map<u64, &'static str> = phf_map! {
81     2u64 => "E",
82     3u64 => "T",
83     4u64 => "F",
84     5u64 => "H",
85     6u64 => "I",
86     7u64 => "J",
87     8u64 => "K",
88     9u64 => "L",
89     10u64 => "M",
90};
91
92/// Set of groups for which the principal axes are forced to be $`C_3`$ to satisfy conventions.
93pub static FORCED_C3_PRINCIPAL_GROUPS: phf::Set<&'static str> = phf_set! {
94    "O",
95    "Oh",
96    "Td",
97    "O + θ·O",
98    "Oh + θ·Oh",
99    "Td + θ·Td",
100    "u[O + θ·O]",
101    "u[Oh + θ·Oh]",
102    "u[Td + θ·Td]",
103    "O*",
104    "Oh*",
105    "Td*",
106    "(O + θ·O)*",
107    "(Oh + θ·Oh)*",
108    "(Td + θ·Td)*",
109    "u[(O + θ·O)*]",
110    "u[(Oh + θ·Oh)*]",
111    "u[(Td + θ·Td)*]",
112};
113
114// ==================
115// Struct definitions
116// ==================
117
118// -------------------
119// MullikenIrrepSymbol
120// -------------------
121
122/// Structure to handle Mulliken irreducible representation symbols.
123#[derive(Builder, Debug, Clone, PartialEq, Eq, Hash, PartialOrd, Serialize, Deserialize)]
124pub struct MullikenIrrepSymbol {
125    /// The generic part of the symbol.
126    generic_symbol: GenericSymbol,
127}
128
129impl MullikenIrrepSymbol {
130    pub fn builder() -> MullikenIrrepSymbolBuilder {
131        MullikenIrrepSymbolBuilder::default()
132    }
133
134    /// Parses a string representing a Mulliken irrep symbol.
135    ///
136    /// Some permissible Mulliken irrep symbols:
137    ///
138    /// ```text
139    /// "T"
140    /// "||T|_(2g)|"
141    /// "|^(3)|T|_(2g)|"
142    /// ```
143    ///
144    /// # Arguments
145    ///
146    /// * `symstr` - A string to be parsed to give a Mulliken symbol.
147    ///
148    /// # Errors
149    ///
150    /// Errors when the string cannot be parsed as a generic symbol.
151    pub fn new(symstr: &str) -> Result<Self, MullikenIrrepSymbolBuilderError> {
152        Self::from_str(symstr)
153    }
154}
155
156// ---------------------
157// MullikenIrcorepSymbol
158// ---------------------
159
160/// Structure to handle Mulliken irreducible corepresentation symbols.
161#[derive(Builder, Debug, Clone, PartialEq, Eq, Hash, PartialOrd, Serialize, Deserialize)]
162pub struct MullikenIrcorepSymbol {
163    /// The decomposed symbol containing the irreps inducing this ircorep.
164    inducing_irreps: DecomposedSymbol<MullikenIrrepSymbol>,
165}
166
167impl MullikenIrcorepSymbol {
168    /// Returns a builder to construct a new [`MullikenIrcorepSymbol`] structure.
169    pub fn builder() -> MullikenIrcorepSymbolBuilder {
170        MullikenIrcorepSymbolBuilder::default()
171    }
172
173    /// Parses a string representing a Mulliken ircorep symbol, which consists of one or more
174    /// Mulliken irrep symbols joined together by `"+"`.
175    ///
176    /// See [`Self::from_str`] for more information.
177    ///
178    /// # Arguments
179    ///
180    /// * `symstr` - A string to be parsed to give a Mulliken ircorep symbol.
181    ///
182    /// # Errors
183    ///
184    /// Errors when the string cannot be parsed as a Mulliken ircorep symbol.
185    pub fn new(symstr: &str) -> Result<Self, MullikenIrcorepSymbolBuilderError> {
186        Self::from_str(symstr)
187    }
188}
189
190// -------------------
191// SymmetryClassSymbol
192// -------------------
193
194/// Structure to handle conjugacy class symbols.
195#[derive(Builder, Debug, Clone, Serialize, Deserialize)]
196pub struct SymmetryClassSymbol<R: Clone + Serialize> {
197    /// The generic part of the symbol.
198    generic_symbol: GenericSymbol,
199
200    /// One or more representative elements in the class.
201    representatives: Option<Vec<R>>,
202}
203
204impl<R: Clone + Serialize> SymmetryClassSymbol<R> {
205    pub fn builder() -> SymmetryClassSymbolBuilder<R> {
206        SymmetryClassSymbolBuilder::default()
207    }
208
209    /// Creates a class symbol from a string and some representative elements.
210    ///
211    /// Some permissible conjugacy class symbols:
212    ///
213    /// ```text
214    /// "1||C3||"
215    /// "1||C3|^(2)|"
216    /// "12||C2|^(5)|"
217    /// "2||S|^(z)|(α)"
218    /// ```
219    ///
220    /// Note that the prefactor is required.
221    ///
222    /// # Arguments
223    ///
224    /// * `symstr` - A string to be parsed to give a class symbol.
225    /// * `rep` - Optional representative elements for this class.
226    ///
227    /// # Returns
228    ///
229    /// A [`Result`] wrapping the constructed class symbol.
230    ///
231    /// # Panics
232    ///
233    /// Panics when unable to construct a class symbol from the specified string.
234    ///
235    /// # Errors
236    ///
237    /// Errors when the string contains no parsable class size prefactor, or when the string cannot
238    /// be parsed as a generic symbol.
239    pub fn new(symstr: &str, rep: Option<Vec<R>>) -> Result<Self, GenericSymbolParsingError> {
240        let generic_symbol = GenericSymbol::from_str(symstr)?;
241        if generic_symbol.multiplicity().is_none() {
242            Err(GenericSymbolParsingError(format!(
243                "{symstr} contains no class multiplicity prefactor."
244            )))
245        } else {
246            Ok(Self::builder()
247                .generic_symbol(generic_symbol)
248                .representatives(rep)
249                .build()
250                .unwrap_or_else(|_| panic!("Unable to construct a class symbol from `{symstr}`.")))
251        }
252    }
253}
254
255// =====================
256// Trait implementations
257// =====================
258
259// -------------------
260// MullikenIrrepSymbol
261// -------------------
262
263impl MathematicalSymbol for MullikenIrrepSymbol {
264    /// The main part of the symbol, which primarily denotes the dimensionality of the irrep space.
265    fn main(&self) -> String {
266        self.generic_symbol.main()
267    }
268
269    /// The pre-superscript part of the symbol, which can be used to denote antiunitary symmetries
270    /// or spin multiplicities.
271    fn presuper(&self) -> String {
272        self.generic_symbol.presuper()
273    }
274
275    fn presub(&self) -> String {
276        self.generic_symbol.presub()
277    }
278
279    /// The post-superscript part of the symbol, which denotes reflection parity.
280    fn postsuper(&self) -> String {
281        self.generic_symbol.postsuper()
282    }
283
284    /// The post-subscript part of the symbol, which denotes inversion parity when available and
285    /// which disambiguates similar irreps.
286    fn postsub(&self) -> String {
287        self.generic_symbol.postsub()
288    }
289
290    /// The prefactor part of the symbol, which is always `"1"` implicitly because of irreducibility.
291    fn prefactor(&self) -> String {
292        String::new()
293    }
294
295    /// The postfactor part of the symbol, which is always empty.
296    fn postfactor(&self) -> String {
297        String::new()
298    }
299
300    /// The dimensionality of the irreducible representation.
301    fn multiplicity(&self) -> Option<usize> {
302        if let Some(&mult) = MULLIKEN_IRREP_DEGENERACIES.get(&self.main()) {
303            Some(
304                mult.try_into()
305                    .unwrap_or_else(|_| panic!("Unable to convert {mult} to `usize`.")),
306            )
307        } else {
308            None
309        }
310    }
311}
312
313impl FromStr for MullikenIrrepSymbol {
314    type Err = MullikenIrrepSymbolBuilderError;
315
316    /// Parses a string representing a Mulliken irrep symbol.
317    ///
318    /// Some permissible Mulliken irrep symbols:
319    ///
320    /// ```text
321    /// "T"
322    /// "||T|_(2g)|"
323    /// "|^(3)|T|_(2g)|"
324    /// ```
325    /// # Arguments
326    ///
327    /// * `symstr` - A string to be parsed to give a Mulliken symbol.
328    ///
329    /// # Returns
330    ///
331    /// A [`Result`] wrapping the constructed Mulliken symbol.
332    ///
333    /// # Panics
334    ///
335    /// Panics when unable to construct a Mulliken symbol from the specified string.
336    ///
337    /// # Errors
338    ///
339    /// Errors when the string cannot be parsed as a generic symbol.
340    fn from_str(symstr: &str) -> Result<Self, Self::Err> {
341        let generic_symbol = GenericSymbol::from_str(symstr)
342            .unwrap_or_else(|_| panic!("Unable to parse {symstr} as a generic symbol."));
343        Self::builder().generic_symbol(generic_symbol).build()
344    }
345}
346
347impl LinearSpaceSymbol for MullikenIrrepSymbol {
348    fn dimensionality(&self) -> usize {
349        usize::try_from(
350            *MULLIKEN_IRREP_DEGENERACIES
351                .get(&self.main())
352                .unwrap_or_else(|| {
353                    panic!(
354                        "Unknown dimensionality for Mulliken symbol {}.",
355                        self.main()
356                    )
357                }),
358        )
359        .expect("Unable to convert the dimensionality of this irrep to `usize`.")
360    }
361
362    fn set_dimensionality(&mut self, dim: usize) -> bool {
363        let dim_u64 = u64::try_from(dim).unwrap_or_else(|err| {
364            log::error!("{err}");
365            panic!("Unable to convert `{dim}` to `u64`.")
366        });
367        let dim_from_current_main = *MULLIKEN_IRREP_DEGENERACIES.get(&self.main()).unwrap_or(&0);
368        if dim_from_current_main == dim_u64 {
369            log::debug!(
370                "The current main symbol `{}` already has the right dimension of `{dim}`. No new Mulliken symbols will be set.",
371                self.main()
372            );
373            false
374        } else {
375            let main_opt = INV_MULLIKEN_IRREP_DEGENERACIES.get(&dim_u64);
376            if let Some(main) = main_opt {
377                self.generic_symbol.set_main(main);
378                true
379            } else {
380                log::warn!(
381                    "Unable to retrieve an unambiguous Mulliken symbol for dimensionality `{dim_u64}`. Main symbol of {self} will be kept unchanged."
382                );
383                false
384            }
385        }
386    }
387}
388
389impl fmt::Display for MullikenIrrepSymbol {
390    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
391        write!(f, "{}", self.generic_symbol)
392    }
393}
394
395// ---------------------
396// MullikenIrcorepSymbol
397// ---------------------
398
399impl From<DecomposedSymbolBuilderError> for MullikenIrcorepSymbolBuilderError {
400    fn from(value: DecomposedSymbolBuilderError) -> Self {
401        match value {
402            DecomposedSymbolBuilderError::UninitializedField(msg) => {
403                MullikenIrcorepSymbolBuilderError::UninitializedField(msg)
404            }
405            DecomposedSymbolBuilderError::ValidationError(msg) => {
406                MullikenIrcorepSymbolBuilderError::ValidationError(msg)
407            }
408        }
409    }
410}
411
412impl FromStr for MullikenIrcorepSymbol {
413    type Err = MullikenIrcorepSymbolBuilderError;
414
415    /// Parses a string representing a Mulliken ircorep symbol. A valid string representing a
416    /// Mulliken ircorep symbol is one consisting of one or more Mulliken irrep symbol strings,
417    /// separated by a `+` character.
418    ///
419    /// # Arguments
420    ///
421    /// * `symstr` - A string to be parsed to give a Mulliken ircorep symbol.
422    ///
423    /// # Returns
424    ///
425    /// A [`Result`] wrapping the constructed Mulliken ircorep symbol.
426    ///
427    /// # Panics
428    ///
429    /// Panics when unable to construct a Mulliken ircorep symbol from the specified string.
430    ///
431    /// # Errors
432    ///
433    /// Errors when the string cannot be parsed.
434    fn from_str(symstr: &str) -> Result<Self, Self::Err> {
435        MullikenIrcorepSymbol::builder()
436            .inducing_irreps(DecomposedSymbol::from_str(symstr)?)
437            .build()
438    }
439}
440
441impl ReducibleLinearSpaceSymbol for MullikenIrcorepSymbol {
442    type Subspace = MullikenIrrepSymbol;
443
444    fn from_subspaces(irreps: &[(Self::Subspace, usize)]) -> Self {
445        Self::builder()
446            .inducing_irreps(DecomposedSymbol::from_subspaces(irreps))
447            .build()
448            .expect("Unable to construct a Mulliken ircorep symbol from a slice of irrep symbols.")
449    }
450
451    fn subspaces(&self) -> Vec<(&Self::Subspace, &usize)> {
452        self.inducing_irreps.subspaces()
453    }
454}
455
456impl fmt::Display for MullikenIrcorepSymbol {
457    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
458        write!(f, "D[{}]", self.main())
459    }
460}
461
462// -------------------
463// SymmetryClassSymbol
464// -------------------
465
466impl<R: Clone + Serialize> PartialEq for SymmetryClassSymbol<R> {
467    fn eq(&self, other: &Self) -> bool {
468        self.generic_symbol == other.generic_symbol
469    }
470}
471
472impl<R: Clone + Serialize> Eq for SymmetryClassSymbol<R> {}
473
474impl<R: Clone + Serialize> Hash for SymmetryClassSymbol<R> {
475    fn hash<H: Hasher>(&self, state: &mut H) {
476        self.generic_symbol.hash(state);
477    }
478}
479
480impl<R: Clone + Serialize> MathematicalSymbol for SymmetryClassSymbol<R> {
481    /// The main part of the symbol, which denotes the representative symmetry operation.
482    fn main(&self) -> String {
483        self.generic_symbol.main()
484    }
485
486    /// The pre-superscript part of the symbol, which is empty.
487    fn presuper(&self) -> String {
488        String::new()
489    }
490
491    /// The pre-subscript part of the symbol, which is empty.
492    fn presub(&self) -> String {
493        String::new()
494    }
495
496    /// The post-superscript part of the symbol, which is empty.
497    fn postsuper(&self) -> String {
498        String::new()
499    }
500
501    /// The post-subscript part of the symbol, which is empty.
502    fn postsub(&self) -> String {
503        String::new()
504    }
505
506    /// The prefactor part of the symbol, which denotes the size of the class.
507    fn prefactor(&self) -> String {
508        self.generic_symbol.prefactor()
509    }
510
511    /// The postfactor part of the symbol, which is empty.
512    fn postfactor(&self) -> String {
513        String::new()
514    }
515
516    /// The number of times the representative elements are 'duplicated' to give the size of the
517    /// class.
518    fn multiplicity(&self) -> Option<usize> {
519        self.generic_symbol.multiplicity()
520    }
521}
522
523impl<R: Clone + Serialize> CollectionSymbol for SymmetryClassSymbol<R> {
524    type CollectionElement = R;
525
526    fn from_reps(
527        symstr: &str,
528        reps: Option<Vec<Self::CollectionElement>>,
529    ) -> Result<Self, GenericSymbolParsingError> {
530        Self::new(symstr, reps)
531    }
532
533    fn representative(&self) -> Option<&Self::CollectionElement> {
534        self.representatives.as_ref().map(|reps| &reps[0])
535    }
536
537    fn representatives(&self) -> Option<&Vec<Self::CollectionElement>> {
538        self.representatives.as_ref()
539    }
540}
541
542impl<R: SpecialSymmetryTransformation + Clone + Serialize> SpecialSymmetryTransformation
543    for SymmetryClassSymbol<R>
544{
545    // ============
546    // Spatial part
547    // ============
548
549    /// Checks if this class is proper.
550    ///
551    /// # Returns
552    ///
553    /// A flag indicating if this class is proper.
554    fn is_proper(&self) -> bool {
555        self.representative()
556            .as_ref()
557            .expect("No representative element found for this class.")
558            .is_proper()
559    }
560
561    fn is_spatial_identity(&self) -> bool {
562        self.representative()
563            .as_ref()
564            .expect("No representative element found for this class.")
565            .is_spatial_identity()
566    }
567
568    fn is_spatial_binary_rotation(&self) -> bool {
569        self.representative()
570            .as_ref()
571            .expect("No representative element found for this class.")
572            .is_spatial_binary_rotation()
573    }
574
575    fn is_spatial_inversion(&self) -> bool {
576        self.representative()
577            .as_ref()
578            .expect("No representative element found for this class.")
579            .is_spatial_inversion()
580    }
581
582    fn is_spatial_reflection(&self) -> bool {
583        self.representative()
584            .as_ref()
585            .expect("No representative element found for this class.")
586            .is_spatial_reflection()
587    }
588
589    // ==================
590    // Time-reversal part
591    // ==================
592
593    /// Checks if this class contains a time reversal operation.
594    ///
595    /// # Returns
596    ///
597    /// A flag indicating if this class contains a time reversal operation.
598    fn contains_time_reversal(&self) -> bool {
599        self.representative()
600            .as_ref()
601            .expect("No representative element found for this class.")
602            .contains_time_reversal()
603    }
604
605    // ==================
606    // Spin rotation part
607    // ==================
608
609    /// Checks if this class contains an active associated spin rotation (normal or inverse).
610    ///
611    /// # Returns
612    ///
613    /// A flag indicating if this class contains an active associated spin rotation.
614    fn is_su2(&self) -> bool {
615        self.representative()
616            .as_ref()
617            .expect("No representative element found for this class.")
618            .is_su2()
619    }
620
621    /// Checks if this class contains an active and inverse associated spin rotation.
622    ///
623    /// # Returns
624    ///
625    /// A flag indicating if this class contains an active and inverse associated spin rotation.
626    fn is_su2_class_1(&self) -> bool {
627        self.representatives()
628            .expect("No representative element found for this class.")
629            .iter()
630            .all(|rep| rep.is_su2_class_1())
631    }
632}
633
634impl<R: FiniteOrder + Clone + Serialize> FiniteOrder for SymmetryClassSymbol<R> {
635    type Int = R::Int;
636
637    fn order(&self) -> Self::Int {
638        self.representative()
639            .as_ref()
640            .expect("No representative element found for this class.")
641            .order()
642    }
643}
644
645impl<R: Clone + Serialize> fmt::Display for SymmetryClassSymbol<R> {
646    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
647        write!(f, "{}", self.generic_symbol)
648    }
649}
650
651// =======
652// Methods
653// =======
654
655/// Reorder the rows so that the characters are in increasing order as we
656/// go down the table, with the first column being used as the primary sort
657/// order, then the second column, and so on, with the exceptions of some
658/// special columns, if available. These are:
659///
660/// * time-reversal class, $`\theta`$,
661/// * inversion class, $`i`$, or horizonal mirror plane $`\sigma_h`$ if $`i`$ not available.
662///
663/// # Arguments
664///
665/// * `char_arr` - A view of a two-dimensional square array containing the characters where
666///   each column is for one conjugacy class and each row one irrep.
667/// * `class_symbols` - An index map containing the conjugacy class symbols for the columns of
668///   `char_arr`. The keys are the symbols and the values are the column indices.
669///
670/// # Returns
671///
672/// A new array with the rows correctly sorted.
673///
674/// # Panics
675///
676/// Panics when expected classes cannot be found.
677pub(super) fn sort_irreps<R: Clone + Serialize + SpecialSymmetryTransformation>(
678    char_arr: &ArrayView2<Character>,
679    frobenius_schur_indicators: &[i8],
680    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
681    principal_classes: &[SymmetryClassSymbol<R>],
682) -> (Array2<Character>, Vec<i8>) {
683    log::debug!("Sorting irreducible representations...");
684    let su2_0 = if class_symbols.keys().any(|cc_sym| cc_sym.is_su2()) {
685        "(Σ)"
686    } else {
687        ""
688    };
689    let class_e = class_symbols
690        .first()
691        .expect("No class symbols found.")
692        .0
693        .clone();
694    let class_e1: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||E(QΣ)||", None)
695        .expect("Unable to construct class symbol `1||E(QΣ)||`.");
696    let class_i: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||i{su2_0}||"), None)
697        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||i{su2_0}||`."));
698    let class_ti: SymmetryClassSymbol<R> =
699        SymmetryClassSymbol::new(&format!("1||θ·i{su2_0}||"), None)
700            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·i{su2_0}||`."));
701    let class_s: SymmetryClassSymbol<R> =
702        SymmetryClassSymbol::new(&format!("1||σh{su2_0}||"), None)
703            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh{su2_0}||`."));
704    let class_s2: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||σh(Σ), σh(QΣ)||", None)
705        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh(Σ), σh(QΣ)||`."));
706    let class_ts: SymmetryClassSymbol<R> =
707        SymmetryClassSymbol::new(&format!("1||θ·σh{su2_0}||"), None)
708            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·σh{su2_0}||`."));
709    let class_t: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||θ{su2_0}||"), None)
710        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ{su2_0}||`."));
711
712    let mut leading_classes: IndexSet<SymmetryClassSymbol<R>> = IndexSet::new();
713    let mut sign_only_classes: HashSet<SymmetryClassSymbol<R>> = HashSet::new();
714
715    // Highest priority: SU(2) class 1
716    if class_symbols.contains_key(&class_e1) {
717        leading_classes.insert(class_e1.clone());
718        sign_only_classes.insert(class_e1);
719    }
720
721    // Second highest priority: time-reversal
722    if class_symbols.contains_key(&class_t) {
723        leading_classes.insert(class_t.clone());
724        sign_only_classes.insert(class_t);
725    }
726
727    // Third highest priority: inversion, or horizontal mirror plane if inversion not available,
728    // or time-reversed horizontal mirror plane if non-time-reversed version not available.
729    if class_symbols.contains_key(&class_i) {
730        leading_classes.insert(class_i.clone());
731        sign_only_classes.insert(class_i);
732    } else if class_symbols.contains_key(&class_ti) {
733        leading_classes.insert(class_ti.clone());
734        sign_only_classes.insert(class_ti);
735    } else if class_symbols.contains_key(&class_s) {
736        leading_classes.insert(class_s.clone());
737        sign_only_classes.insert(class_s);
738    } else if class_symbols.contains_key(&class_s2) {
739        leading_classes.insert(class_s2.clone());
740        sign_only_classes.insert(class_s2);
741    } else if class_symbols.contains_key(&class_ts) {
742        leading_classes.insert(class_ts.clone());
743        sign_only_classes.insert(class_ts);
744    };
745
746    // Forth highest priority: identity
747    leading_classes.insert(class_e);
748
749    // Fifth highest priority: principal classes, if not yet encountered
750    leading_classes.extend(principal_classes.iter().cloned());
751
752    log::debug!("Irreducible representation sort order:");
753    for leading_cc in leading_classes.iter() {
754        log::debug!(
755            "  {}{}",
756            leading_cc,
757            if sign_only_classes.contains(leading_cc) {
758                " (sign only)"
759            } else {
760                ""
761            }
762        );
763    }
764
765    let leading_idxs: IndexSet<usize> = leading_classes
766        .iter()
767        .map(|cc| {
768            *class_symbols
769                .get(cc)
770                .unwrap_or_else(|| panic!("Class `{cc}` cannot be found."))
771        })
772        .collect();
773
774    let n_rows = char_arr.nrows();
775    let mut col_idxs: Vec<usize> = Vec::with_capacity(n_rows);
776    col_idxs.extend(leading_idxs.iter());
777    col_idxs.extend((1..n_rows).filter(|i| !leading_idxs.contains(i)));
778    let mut sort_arr = char_arr.select(Axis(1), &col_idxs);
779    let one = Character::new(&[(UnityRoot::new(0, 2), 1)]);
780    let m_one = Character::new(&[(UnityRoot::new(1, 2), 1)]);
781    sign_only_classes.iter().for_each(|class| {
782        let col_idx = leading_classes
783            .get_index_of(class)
784            .unwrap_or_else(|| panic!("Unable to obtain the column index of class `{class}`."));
785        sort_arr.column_mut(col_idx).mapv_inplace(|character| {
786            let character_c = character.complex_value();
787            if approx::relative_eq!(
788                character_c.im,
789                0.0,
790                max_relative = character.threshold(),
791                epsilon = character.threshold()
792            ) {
793                if character_c.re > 0.0 {
794                    one.clone()
795                } else {
796                    m_one.clone()
797                }
798            } else if approx::relative_eq!(
799                character_c.re,
800                0.0,
801                max_relative = character.threshold(),
802                epsilon = character.threshold()
803            ) {
804                if character_c.im > 0.0 {
805                    one.clone()
806                } else {
807                    m_one.clone()
808                }
809            } else {
810                panic!("Character {character} is neither purely real nor purely imaginary for sign-only sorting.")
811            }
812        });
813    });
814    // let sort_arr = if has_e1 {
815    //     // Maps E(QΣ) characters to ±1, otherwise all linear/projective irreps of the same
816    //     // degeneracy would be grouped together.
817    //     let mut sort_arr = char_arr.select(Axis(1), &col_idxs);
818    //     let one = Character::new(&[(UnityRoot::new(0, 2), 1)]);
819    //     let m_one = Character::new(&[(UnityRoot::new(1, 2), 1)]);
820
821    //     sort_arr.column_mut(0).mapv_inplace(|e1_character| {
822    //         if e1_character.complex_value().re > 0.0 {
823    //             one.clone()
824    //         } else {
825    //             m_one.clone()
826    //         }
827    //     });
828    //     sort_arr
829    // } else {
830    //     char_arr.select(Axis(1), &col_idxs)
831    // };
832
833    let sort_row_indices: Vec<_> = (0..n_rows)
834        .sorted_by(|&i, &j| {
835            let keys_i = sort_arr.row(i).iter().cloned().collect_vec();
836            let keys_j = sort_arr.row(j).iter().cloned().collect_vec();
837            keys_i
838                .partial_cmp(&keys_j)
839                .unwrap_or_else(|| panic!("`{keys_i:?}` and `{keys_j:?}` cannot be compared."))
840        })
841        .collect();
842    let char_arr = char_arr.select(Axis(0), &sort_row_indices);
843    let old_fs = frobenius_schur_indicators.iter().collect::<Vec<_>>();
844    let sorted_fs = sort_row_indices.iter().map(|&i| *old_fs[i]).collect_vec();
845    log::debug!("Sorting irreducible representations... Done.");
846    (char_arr, sorted_fs)
847}
848
849/// Determines the principal classes given a list of class symbols and any forcing conditions.
850///
851/// By default, the principal classes are those with the highest order, regardless of whether they
852/// are proper or improper, unitary or antiunitary.
853///
854/// # Arguments
855///
856/// * `class_symbols` - An indexmap of class symbols and their corresponding indices.
857/// * `force_proper_principal` - A flag indicating if the principal classes are forced to be
858///   proper.
859/// * `force_principal` - An option containing specific classes that are forced to be principal.
860///
861/// # Returns
862///
863/// A vector of symbols of principal classes.
864///
865/// # Panics
866///
867/// Panics when:
868///
869/// * both `force_proper_principal` and `force_principal` are specified;
870/// * classes specified in `force_principal` cannot be found in `class_symbols`;
871/// * no principal classes can be found.
872pub(super) fn deduce_principal_classes<R, P>(
873    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
874    force_principal_predicate: Option<P>,
875    force_principal: Option<SymmetryClassSymbol<R>>,
876) -> Vec<SymmetryClassSymbol<R>>
877where
878    R: fmt::Debug + SpecialSymmetryTransformation + FiniteOrder + Clone + Serialize,
879    P: Copy + Fn(&SymmetryClassSymbol<R>) -> bool,
880{
881    log::debug!("Determining principal classes...");
882    let principal_classes = force_principal.map_or_else(|| {
883        // sorted_class_symbols contains class symbols sorted in the order:
884        // - decreasing operation order
885        // - unitary operations, then antiunitary operations
886        // - proper operations, then improper operations in each operation order
887        // - if SU(2), then homotopy class 0, then homotopy class 1
888        let mut sorted_class_symbols = class_symbols
889            .clone()
890            .sorted_by(|cc1, _, cc2, _| {
891                PartialOrd::partial_cmp(
892                    &(cc1.order(), !cc1.contains_time_reversal(), cc1.is_proper(), !cc1.is_su2_class_1()),
893                    &(cc2.order(), !cc2.contains_time_reversal(), cc2.is_proper(), !cc2.is_su2_class_1()),
894                )
895                .expect("Unable to sort class symbols.")
896            })
897            .rev();
898        let (principal_rep, _) = if let Some(predicate) = force_principal_predicate {
899            // Larger order always prioritised, regardless of proper or improper,
900            // unless force_proper_principal is set, in which case
901            sorted_class_symbols
902                .find(|(cc, _)| predicate(cc))
903                .expect("`Unable to find classes fulfilling the specified predicate.`")
904        } else {
905            // No force_proper_principal, so proper prioritised, which has been ensured by the sort
906            // in the construction of `sorted_class_symbols`.
907            sorted_class_symbols
908                .next()
909                .expect("Unexpected empty `sorted_class_symbols`.")
910        };
911
912        // Now find all principal classes with the same order, parity, unitarity, and homotopy
913        // class as `principal_rep`.
914        class_symbols
915            .keys()
916            .filter(|cc| {
917                cc.contains_time_reversal() == principal_rep.contains_time_reversal()
918                && cc.is_proper() == principal_rep.is_proper()
919                && cc.order() == principal_rep.order()
920                && cc.is_su2_class_1() == principal_rep.is_su2_class_1()
921            })
922            .cloned()
923            .collect::<Vec<_>>()
924    }, |principal_cc| {
925        assert!(
926            force_principal_predicate.is_none(),
927            "`force_principal_predicate` and `force_principal` cannot be both provided."
928        );
929        assert!(
930            class_symbols.contains_key(&principal_cc),
931            "Forcing principal-axis class to be `{principal_cc}`, but `{principal_cc}` is not a valid class in the group."
932        );
933
934        log::warn!(
935            "Principal-axis class forced to be `{principal_cc}`. Auto-detection of principal-axis classes will be skipped.",
936        );
937        vec![principal_cc]
938    });
939
940    assert!(!principal_classes.is_empty());
941    if principal_classes.len() == 1 {
942        log::debug!("Principal-axis class found: {}", principal_classes[0]);
943    } else {
944        log::debug!("Principal-axis classes found:");
945        for princc in &principal_classes {
946            log::debug!("  {}", princc);
947        }
948    }
949    log::debug!("Determining principal classes... Done.");
950    principal_classes
951}
952
953/// Deduces irreducible representation symbols based on Mulliken's convention.
954///
955/// # Arguments
956///
957/// * `char_arr` - A view of a two-dimensional square array containing the characters where
958///   each column is for one conjugacy class and each row one irrep.
959/// * `class_symbols` - An index map containing the conjugacy class symbols for the columns of
960///   `char_arr`. The keys are the symbols and the values are the column indices.
961/// * `principal_classes` - The principal classes to be used for irrep symbol deduction.
962///
963/// # Returns
964///
965/// A vector of Mulliken symbols corresponding to the rows of `char_arr`.
966///
967/// # Panics
968///
969/// Panics when expected classes cannot be found in `class_symbols`.
970#[allow(clippy::too_many_lines)]
971pub(super) fn deduce_mulliken_irrep_symbols<R>(
972    char_arr: &ArrayView2<Character>,
973    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
974    principal_classes: &[SymmetryClassSymbol<R>],
975) -> Vec<MullikenIrrepSymbol>
976where
977    R: fmt::Debug + SpecialSymmetryTransformation + FiniteOrder + Clone + Serialize,
978{
979    log::debug!("Generating Mulliken irreducible representation symbols...");
980
981    let su2_0 = if class_symbols.keys().any(|cc_sym| cc_sym.is_su2()) {
982        "(Σ)"
983    } else {
984        ""
985    };
986
987    let e_cc = class_symbols
988        .first()
989        .expect("No class symbols found.")
990        .0
991        .clone();
992    let e1_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||E(QΣ)||", None)
993        .expect("Unable to construct class symbol `1||E(QΣ)||`.");
994    let i_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||i{su2_0}||"), None)
995        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||i{su2_0}||`."));
996    let ti_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||θ·i{su2_0}||"), None)
997        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·i{su2_0}||`."));
998    let s_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||σh{su2_0}||"), None)
999        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh{su2_0}||`."));
1000    let s2_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||σh(Σ), σh(QΣ)||", None)
1001        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh(Σ), σh(QΣ)||`."));
1002    let ts_cc: SymmetryClassSymbol<R> =
1003        SymmetryClassSymbol::new(&format!("1||θ·σh{su2_0}||"), None)
1004            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·σh{su2_0}||`."));
1005    let t_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||θ{su2_0}||"), None)
1006        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ{su2_0}||`."));
1007
1008    // SU(2) class 1 parity?
1009    let su2_1_parity = class_symbols.contains_key(&e1_cc);
1010
1011    // Inversion parity?
1012    let i_parity = class_symbols.contains_key(&i_cc);
1013
1014    // Time-reversed inversion parity?
1015    let ti_parity = class_symbols.contains_key(&ti_cc);
1016
1017    // Reflection parity?
1018    let s_parity = class_symbols.contains_key(&s_cc) || class_symbols.contains_key(&s2_cc);
1019
1020    // Time-reversed reflection parity?
1021    let ts_parity = class_symbols.contains_key(&ts_cc);
1022
1023    // SU(2) class 1?
1024    if su2_1_parity {
1025        log::debug!("E(QΣ) found. This will be used to determine projective representations.");
1026    }
1027
1028    // i_parity takes priority.
1029    if i_parity {
1030        log::debug!("Inversion centre found. This will be used for g/u ordering.");
1031    } else if ti_parity {
1032        log::debug!(
1033            "Time-reversed inversion centre found (but no inversion centre). This will be used for g/u ordering."
1034        );
1035    } else if s_parity {
1036        log::debug!(
1037            "Horizontal mirror plane found (but no inversion centre). This will be used for '/'' ordering."
1038        );
1039    } else if ts_parity {
1040        log::debug!(
1041            "Time-reversed horizontal mirror plane found (but no inversion centre). This will be used for '/'' ordering."
1042        );
1043    }
1044
1045    // Time-reversal?
1046    let t_parity = class_symbols.contains_key(&t_cc);
1047    if t_parity {
1048        log::debug!("Time reversal found. This will be used for magnetic ordering.");
1049    }
1050
1051    let e2p1 = UnityRoot::new(1u32, 2u32);
1052    let e2p2 = UnityRoot::new(2u32, 2u32);
1053    let char_p1 = Character::new(&[(e2p2, 1usize)]);
1054    let char_m1 = Character::new(&[(e2p1, 1usize)]);
1055
1056    // First pass: assign irrep symbols based on Mulliken's convention as much as possible.
1057    log::debug!("First pass: assign symbols from rules");
1058
1059    let mut raw_irrep_symbols = char_arr.rows().into_iter().map(|irrep| {
1060        // Determine the main symmetry
1061        let dim = irrep[
1062            *class_symbols.get(&e_cc).unwrap_or_else(|| panic!("Class `{}` not found.", &e_cc))
1063        ].complex_value();
1064        assert!(
1065            approx::relative_eq!(dim.im, 0.0)
1066                && approx::relative_eq!(dim.re.round(), dim.re)
1067                && dim.re.round() > 0.0
1068        );
1069
1070        assert!(dim.re.round() > 0.0);
1071        #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
1072        let dim = dim.re.round() as u64;
1073        assert!(dim > 0);
1074        let main = if dim >= 2 {
1075            // Degenerate irreps
1076            INV_MULLIKEN_IRREP_DEGENERACIES.get(&dim).unwrap_or_else(|| {
1077                log::warn!("{} cannot be assigned a standard dimensionality symbol. A generic 'Λ' will be used instead.", dim);
1078                &"Λ"
1079            })
1080        } else {
1081            // Non-degenerate irreps
1082            let complex = irrep.map(|character| character.complex_conjugate()) != irrep;
1083            if complex {
1084                "Γ"
1085            } else {
1086                let char_rots: HashSet<_> = principal_classes
1087                    .iter()
1088                    .map(|cc| {
1089                        irrep[
1090                            *class_symbols.get(cc).unwrap_or_else(|| panic!("Class `{cc}` not found."))
1091                        ].clone()
1092                    })
1093                    .collect();
1094                if char_rots.len() == 1 && *char_rots.iter().next().expect("No rotation classes found.") == char_p1 {
1095                    "A"
1096                } else if char_rots
1097                    .iter()
1098                    .all(|char_rot| char_rot.clone() == char_p1 || char_rot.clone() == char_m1)
1099                {
1100                    "B"
1101                } else {
1102                    panic!("There are principal rotations but with unexpected non-(±1) characters.");
1103                    // // There are principal rotations but with non-(±1) characters. These must be
1104                    // // complex.
1105                    // "Γ"
1106                }
1107            }
1108        };
1109
1110        let projective = if su2_1_parity {
1111            let char_e1 = irrep[
1112                *class_symbols
1113                    .get(&e1_cc)
1114                    .unwrap_or_else(|| {
1115                        panic!("The character for `{e1_cc}` could not be found.")
1116                    })
1117            ].clone();
1118            let char_e1_c = char_e1.complex_value();
1119            assert!(
1120                approx::relative_eq!(
1121                    char_e1_c.im,
1122                    0.0,
1123                    epsilon = char_e1.threshold(),
1124                    max_relative = char_e1.threshold()
1125                ) && approx::relative_eq!(
1126                    char_e1_c.re.round(),
1127                    char_e1_c.re,
1128                    epsilon = char_e1.threshold(),
1129                    max_relative = char_e1.threshold()
1130                ),
1131            );
1132
1133            #[allow(clippy::cast_possible_truncation)]
1134            let char_e1_c = char_e1_c
1135                .re
1136                .round()
1137                .to_i32()
1138                .unwrap_or_else(|| panic!("Unable to convert the real part of the character for `{e1_cc}` to `i32`."));
1139            let dim_i32 = i32::try_from(dim).expect("Unable to convert the irrep dimensionality to `i32`.");
1140            if char_e1_c == dim_i32 {
1141                false
1142            } else {
1143                assert_eq!(char_e1_c, -dim_i32);
1144                true
1145            }
1146        } else {
1147            false
1148        };
1149        let projective_str = if projective { "~" } else { "" };
1150
1151        let (inv, mir) = if i_parity || ti_parity {
1152            // Determine inversion symmetry
1153            // Inversion symmetry trumps reflection symmetry.
1154            let char_inv = irrep[
1155                *class_symbols
1156                    .get(&i_cc)
1157                    .unwrap_or_else(|| {
1158                        class_symbols
1159                            .get(&ti_cc)
1160                            .unwrap_or_else(|| {
1161                                panic!("Neither `{}` nor `{}` found.", &i_cc, &ti_cc)
1162                            })
1163                    })
1164            ].clone();
1165            let char_inv_c = char_inv.complex_value();
1166            assert!(
1167                approx::relative_eq!(
1168                    char_inv_c.im,
1169                    0.0,
1170                    epsilon = char_inv.threshold(),
1171                    max_relative = char_inv.threshold()
1172                ) && approx::relative_eq!(
1173                    char_inv_c.re.round(),
1174                    char_inv_c.re,
1175                    epsilon = char_inv.threshold(),
1176                    max_relative = char_inv.threshold()
1177                ),
1178            );
1179
1180            #[allow(clippy::cast_possible_truncation)]
1181            let char_inv_c = char_inv_c.re.round() as i32;
1182            match char_inv_c.cmp(&0) {
1183                Ordering::Greater => ("g", ""),
1184                Ordering::Less => ("u", ""),
1185                Ordering::Equal => panic!("Inversion character must not be zero."),
1186            }
1187        } else if !projective && (s_parity || ts_parity) {
1188            // Determine reflection symmetry
1189            // Note that for projective irreps, characters under σh are either imaginary or zero,
1190            // thus no '/'' classifications.
1191            let char_ref = irrep[
1192                *class_symbols
1193                    .get(&s_cc)
1194                    .unwrap_or_else(|| {
1195                        class_symbols
1196                            .get(&s2_cc)
1197                            .unwrap_or_else(|| {
1198                                class_symbols
1199                                    .get(&ts_cc)
1200                                    .unwrap_or_else(|| {
1201                                        panic!("Neither `{}`, `{}`, nor `{}` found.", &s_cc, &s2_cc, &ts_cc)
1202                                        })
1203                        })
1204                    })
1205            ].clone();
1206            let char_ref_c = char_ref.complex_value();
1207            assert!(
1208                approx::relative_eq!(
1209                    char_ref_c.im,
1210                    0.0,
1211                    epsilon = char_ref.threshold(),
1212                    max_relative = char_ref.threshold()
1213                ) && approx::relative_eq!(
1214                    char_ref_c.re.round(),
1215                    char_ref_c.re,
1216                    epsilon = char_ref.threshold(),
1217                    max_relative = char_ref.threshold()
1218                ),
1219            );
1220
1221            #[allow(clippy::cast_possible_truncation)]
1222            let char_ref_c = char_ref_c.re.round() as i32;
1223            match char_ref_c.cmp(&0) {
1224                Ordering::Greater => ("", "'"),
1225                Ordering::Less => ("", "''"),
1226                Ordering::Equal => panic!("Reflection character must not be zero for linear irreducible representations."),
1227            }
1228        } else {
1229            ("", "")
1230        };
1231
1232
1233        let trev = if t_parity {
1234            // Determine time-reversal symmetry
1235            let char_trev = irrep[
1236                *class_symbols.get(&t_cc).unwrap_or_else(|| {
1237                    panic!("Class `{}` not found.", &t_cc)
1238                })
1239            ].clone();
1240            let char_trev_c = char_trev.complex_value();
1241            if approx::relative_eq!(
1242                char_trev_c.im,
1243                0.0,
1244                epsilon = char_trev.threshold(),
1245                max_relative = char_trev.threshold()
1246            ) && approx::relative_eq!(
1247                char_trev_c.re.round(),
1248                char_trev_c.re,
1249                epsilon = char_trev.threshold(),
1250                max_relative = char_trev.threshold()
1251            ) {
1252                // Real, integral time-reversal character
1253                #[allow(clippy::cast_possible_truncation)]
1254                let char_trev_c = char_trev_c.re.round() as i32;
1255                match char_trev_c.cmp(&0) {
1256                    Ordering::Greater => "+",
1257                    Ordering::Less => "-",
1258                    Ordering::Equal => panic!("Real time-reversal character must not be zero."),
1259                }
1260            } else {
1261                // Non-real or non-integral time-reversal character
1262                ""
1263            }
1264        } else {
1265            ""
1266        };
1267
1268        MullikenIrrepSymbol::new(format!("|^({trev})|{main}{projective_str}|^({mir})_({inv})|").as_str())
1269            .unwrap_or_else(|_| {
1270                panic!(
1271                    "Unable to construct symmetry symbol `|^({trev})|{main}{projective_str}|^({mir})_({inv})|`."
1272                )
1273            })
1274    }).collect_vec();
1275
1276    let mut complex_irrep_indices = char_arr
1277        .rows()
1278        .into_iter()
1279        .enumerate()
1280        .filter_map(|(i, irrep)| {
1281            let complex = irrep.map(|character| character.complex_conjugate()) != irrep;
1282            if complex { Some(i) } else { None }
1283        })
1284        .collect::<IndexSet<_>>();
1285    complex_irrep_indices.reverse();
1286
1287    let cc_pairs = if !complex_irrep_indices.is_empty() {
1288        log::debug!("Grouping pairs of complex-conjugate irreps...");
1289        let mut cc_pairs: Vec<(usize, usize)> = vec![];
1290        while !complex_irrep_indices.is_empty() {
1291            let complex_irrep_index = complex_irrep_indices.pop().unwrap();
1292            let complex_irrep = char_arr.row(complex_irrep_index);
1293            let complex_conj_irrep = complex_irrep.map(|character| character.complex_conjugate());
1294            let complex_conj_irrep_index = char_arr
1295                .rows()
1296                .into_iter()
1297                .enumerate()
1298                .find_map(|(i, irrep)| {
1299                    if irrep == complex_conj_irrep {
1300                        Some(i)
1301                    } else {
1302                        None
1303                    }
1304                })
1305                .expect("Unable to find the complex-conjugate irrep.");
1306            complex_irrep_indices.shift_remove(&complex_conj_irrep_index);
1307            cc_pairs.push((complex_irrep_index, complex_conj_irrep_index));
1308            raw_irrep_symbols[complex_irrep_index]
1309                .generic_symbol
1310                .set_presub("a");
1311            raw_irrep_symbols[complex_conj_irrep_index]
1312                .generic_symbol
1313                .set_presub("b");
1314        }
1315        log::debug!(
1316            "There {} {} {} of complex-conjugate irreps.",
1317            if cc_pairs.len() == 1 { "is" } else { "are" },
1318            cc_pairs.len(),
1319            if cc_pairs.len() == 1 { "pair" } else { "pairs" },
1320        );
1321        log::debug!("Grouping pairs of complex-conjugate irreps... Done.");
1322        Some(cc_pairs)
1323    } else {
1324        None
1325    };
1326
1327    log::debug!("Second pass: disambiguate identical cases not distinguishable by rules");
1328    let mut irrep_symbols = disambiguate_linspace_symbols(raw_irrep_symbols.into_iter());
1329
1330    if let Some(cc_pairs_vec) = cc_pairs {
1331        log::debug!("Equalising post-subscripts in complex-conjugate pairs...");
1332        for (complex_irrep_index, complex_conj_irrep_index) in &cc_pairs_vec {
1333            let complex_irrep_postsub = irrep_symbols[*complex_irrep_index].postsub();
1334            irrep_symbols[*complex_conj_irrep_index]
1335                .generic_symbol
1336                .set_postsub(&complex_irrep_postsub);
1337        }
1338        log::debug!("Equalising post-subscripts in complex-conjugate pairs... Done.");
1339    }
1340
1341    log::debug!("Generating Mulliken irreducible representation symbols... Done.");
1342
1343    irrep_symbols
1344}
1345
1346/// Determines the mirror-plane symbol given a principal axis.
1347///
1348/// # Arguments
1349///
1350/// * `sigma_axis` - The normalised normal vector of a mirror plane.
1351/// * `principal_axis` - The normalised principal rotation axis.
1352/// * `thresh` - Threshold for comparisons.
1353/// * `force_d` - Flag indicating if vertical mirror planes should be given the $`d`$ symbol
1354///   instead of $`v`$.
1355///
1356/// # Returns
1357///
1358/// The mirror-plane symbol.
1359#[must_use]
1360pub(super) fn deduce_sigma_symbol(
1361    sigma_axis: &Vector3<f64>,
1362    principal_element: &SymmetryElement,
1363    thresh: f64,
1364    force_d: bool,
1365) -> Option<String> {
1366    if approx::relative_eq!(
1367        principal_element.raw_axis().dot(sigma_axis).abs(),
1368        0.0,
1369        epsilon = thresh,
1370        max_relative = thresh
1371    ) && *principal_element.raw_proper_order() != ORDER_1
1372    {
1373        // Vertical plane containing principal axis
1374        if force_d {
1375            Some("d".to_owned())
1376        } else {
1377            Some("v".to_owned())
1378        }
1379    } else if approx::relative_eq!(
1380        principal_element.raw_axis().cross(sigma_axis).norm(),
1381        0.0,
1382        epsilon = thresh,
1383        max_relative = thresh
1384    ) && *principal_element.raw_proper_order() != ORDER_1
1385    {
1386        // Horizontal plane perpendicular to principal axis
1387        Some("h".to_owned())
1388    } else {
1389        None
1390    }
1391}
1392
1393/// Enumerated type specifying the parity under a mirror plane.
1394#[derive(Clone, PartialEq, Eq, Hash, Debug)]
1395pub enum MirrorParity {
1396    /// Variant for even parity.
1397    Even,
1398
1399    /// Variant for odd parity.
1400    Odd,
1401
1402    /// Variant for no parity.
1403    Neither,
1404}
1405
1406impl fmt::Display for MirrorParity {
1407    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1408        match self {
1409            MirrorParity::Even => write!(f, "(+)"),
1410            MirrorParity::Odd => write!(f, "(-)"),
1411            MirrorParity::Neither => write!(f, "( )"),
1412        }
1413    }
1414}
1415
1416/// Given a group with a character table, deduces the parities of a specified representation under
1417/// all classes in the group that contain a spatial reflection.
1418///
1419/// # Arguments
1420///
1421/// * `group` - A group with a character table.
1422/// * `rep` - A representation (that has been decomposed as irreps or ircoreps of `group`) for
1423///   which mirror parities are to be deduced.
1424///
1425/// # Returns
1426///
1427/// An indexmap whose keys are reflection classes and whose values are the corresponding parities of
1428/// `rep`.
1429///
1430/// # Panics
1431///
1432/// Panics on unexpected errors.
1433pub(crate) fn deduce_mirror_parities<G, R>(
1434    group: &G,
1435    rep: &R,
1436) -> IndexMap<<<G as CharacterProperties>::CharTab as CharacterTable>::ColSymbol, MirrorParity>
1437where
1438    G: SymmetryGroupProperties,
1439    R: ReducibleLinearSpaceSymbol<
1440        Subspace = <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
1441    >,
1442{
1443    let id_cc = group
1444        .get_cc_symbol_of_index(0)
1445        .expect("Unable to obtain the conjugacy class symbol of the identity.");
1446    let mirrors = group.filter_cc_symbols(|cc| cc.is_spatial_reflection());
1447    mirrors
1448        .iter()
1449        .map(|sigma_cc| {
1450            let sigma_cc_chars = rep
1451                .subspaces()
1452                .iter()
1453                .map(|(irrep, _)| {
1454                    let id_char = group.character_table().get_character(irrep, &id_cc);
1455                    let char = group.character_table().get_character(irrep, sigma_cc);
1456                    if *char == *id_char {
1457                        MirrorParity::Even
1458                    } else if *char == -id_char {
1459                        MirrorParity::Odd
1460                    } else {
1461                        MirrorParity::Neither
1462                    }
1463                })
1464                .collect::<HashSet<_>>();
1465            let sigma_cc_parity = if sigma_cc_chars.len() == 1 {
1466                sigma_cc_chars
1467                    .into_iter()
1468                    .next()
1469                    .expect("Unable to extract the mirror parity.")
1470            } else {
1471                MirrorParity::Neither
1472            };
1473            (sigma_cc.clone(), sigma_cc_parity)
1474        })
1475        .collect::<IndexMap<_, _>>()
1476}