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::character::Character;
21use crate::chartab::chartab_group::CharacterProperties;
22use crate::chartab::chartab_symbols::{
23    disambiguate_linspace_symbols, CollectionSymbol, DecomposedSymbol,
24    DecomposedSymbolBuilderError, GenericSymbol, GenericSymbolParsingError, LinearSpaceSymbol,
25    MathematicalSymbol, ReducibleLinearSpaceSymbol,
26};
27use crate::chartab::unityroot::UnityRoot;
28use crate::chartab::CharacterTable;
29use crate::group::FiniteOrder;
30use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
31use crate::symmetry::symmetry_element::SymmetryElement;
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    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    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    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!("Unable to retrieve an unambiguous Mulliken symbol for dimensionality `{dim_u64}`. Main symbol of {self} will be kept unchanged.");
381                false
382            }
383        }
384    }
385}
386
387impl fmt::Display for MullikenIrrepSymbol {
388    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
389        write!(f, "{}", self.generic_symbol)
390    }
391}
392
393// ---------------------
394// MullikenIrcorepSymbol
395// ---------------------
396
397impl From<DecomposedSymbolBuilderError> for MullikenIrcorepSymbolBuilderError {
398    fn from(value: DecomposedSymbolBuilderError) -> Self {
399        match value {
400            DecomposedSymbolBuilderError::UninitializedField(msg) => {
401                MullikenIrcorepSymbolBuilderError::UninitializedField(msg)
402            }
403            DecomposedSymbolBuilderError::ValidationError(msg) => {
404                MullikenIrcorepSymbolBuilderError::ValidationError(msg)
405            }
406        }
407    }
408}
409
410impl FromStr for MullikenIrcorepSymbol {
411    type Err = MullikenIrcorepSymbolBuilderError;
412
413    /// Parses a string representing a Mulliken ircorep symbol. A valid string representing a
414    /// Mulliken ircorep symbol is one consisting of one or more Mulliken irrep symbol strings,
415    /// separated by a `+` character.
416    ///
417    /// # Arguments
418    ///
419    /// * `symstr` - A string to be parsed to give a Mulliken ircorep symbol.
420    ///
421    /// # Returns
422    ///
423    /// A [`Result`] wrapping the constructed Mulliken ircorep symbol.
424    ///
425    /// # Panics
426    ///
427    /// Panics when unable to construct a Mulliken ircorep symbol from the specified string.
428    ///
429    /// # Errors
430    ///
431    /// Errors when the string cannot be parsed.
432    fn from_str(symstr: &str) -> Result<Self, Self::Err> {
433        MullikenIrcorepSymbol::builder()
434            .inducing_irreps(DecomposedSymbol::from_str(symstr)?)
435            .build()
436    }
437}
438
439impl ReducibleLinearSpaceSymbol for MullikenIrcorepSymbol {
440    type Subspace = MullikenIrrepSymbol;
441
442    fn from_subspaces(irreps: &[(Self::Subspace, usize)]) -> Self {
443        Self::builder()
444            .inducing_irreps(DecomposedSymbol::from_subspaces(irreps))
445            .build()
446            .expect("Unable to construct a Mulliken ircorep symbol from a slice of irrep symbols.")
447    }
448
449    fn subspaces(&self) -> Vec<(&Self::Subspace, &usize)> {
450        self.inducing_irreps.subspaces()
451    }
452}
453
454impl fmt::Display for MullikenIrcorepSymbol {
455    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
456        write!(f, "D[{}]", self.main())
457    }
458}
459
460// -------------------
461// SymmetryClassSymbol
462// -------------------
463
464impl<R: Clone + Serialize> PartialEq for SymmetryClassSymbol<R> {
465    fn eq(&self, other: &Self) -> bool {
466        self.generic_symbol == other.generic_symbol
467    }
468}
469
470impl<R: Clone + Serialize> Eq for SymmetryClassSymbol<R> {}
471
472impl<R: Clone + Serialize> Hash for SymmetryClassSymbol<R> {
473    fn hash<H: Hasher>(&self, state: &mut H) {
474        self.generic_symbol.hash(state);
475    }
476}
477
478impl<R: Clone + Serialize> MathematicalSymbol for SymmetryClassSymbol<R> {
479    /// The main part of the symbol, which denotes the representative symmetry operation.
480    fn main(&self) -> String {
481        self.generic_symbol.main()
482    }
483
484    /// The pre-superscript part of the symbol, which is empty.
485    fn presuper(&self) -> String {
486        String::new()
487    }
488
489    /// The pre-subscript part of the symbol, which is empty.
490    fn presub(&self) -> String {
491        String::new()
492    }
493
494    /// The post-superscript part of the symbol, which is empty.
495    fn postsuper(&self) -> String {
496        String::new()
497    }
498
499    /// The post-subscript part of the symbol, which is empty.
500    fn postsub(&self) -> String {
501        String::new()
502    }
503
504    /// The prefactor part of the symbol, which denotes the size of the class.
505    fn prefactor(&self) -> String {
506        self.generic_symbol.prefactor()
507    }
508
509    /// The postfactor part of the symbol, which is empty.
510    fn postfactor(&self) -> String {
511        String::new()
512    }
513
514    /// The number of times the representative elements are 'duplicated' to give the size of the
515    /// class.
516    fn multiplicity(&self) -> Option<usize> {
517        self.generic_symbol.multiplicity()
518    }
519}
520
521impl<R: Clone + Serialize> CollectionSymbol for SymmetryClassSymbol<R> {
522    type CollectionElement = R;
523
524    fn from_reps(
525        symstr: &str,
526        reps: Option<Vec<Self::CollectionElement>>,
527    ) -> Result<Self, GenericSymbolParsingError> {
528        Self::new(symstr, reps)
529    }
530
531    fn representative(&self) -> Option<&Self::CollectionElement> {
532        self.representatives.as_ref().map(|reps| &reps[0])
533    }
534
535    fn representatives(&self) -> Option<&Vec<Self::CollectionElement>> {
536        self.representatives.as_ref()
537    }
538}
539
540impl<R: SpecialSymmetryTransformation + Clone + Serialize> SpecialSymmetryTransformation
541    for SymmetryClassSymbol<R>
542{
543    // ============
544    // Spatial part
545    // ============
546
547    /// Checks if this class is proper.
548    ///
549    /// # Returns
550    ///
551    /// A flag indicating if this class is proper.
552    fn is_proper(&self) -> bool {
553        self.representative()
554            .as_ref()
555            .expect("No representative element found for this class.")
556            .is_proper()
557    }
558
559    fn is_spatial_identity(&self) -> bool {
560        self.representative()
561            .as_ref()
562            .expect("No representative element found for this class.")
563            .is_spatial_identity()
564    }
565
566    fn is_spatial_binary_rotation(&self) -> bool {
567        self.representative()
568            .as_ref()
569            .expect("No representative element found for this class.")
570            .is_spatial_binary_rotation()
571    }
572
573    fn is_spatial_inversion(&self) -> bool {
574        self.representative()
575            .as_ref()
576            .expect("No representative element found for this class.")
577            .is_spatial_inversion()
578    }
579
580    fn is_spatial_reflection(&self) -> bool {
581        self.representative()
582            .as_ref()
583            .expect("No representative element found for this class.")
584            .is_spatial_reflection()
585    }
586
587    // ==================
588    // Time-reversal part
589    // ==================
590
591    /// Checks if this class contains a time reversal operation.
592    ///
593    /// # Returns
594    ///
595    /// A flag indicating if this class contains a time reversal operation.
596    fn contains_time_reversal(&self) -> bool {
597        self.representative()
598            .as_ref()
599            .expect("No representative element found for this class.")
600            .contains_time_reversal()
601    }
602
603    // ==================
604    // Spin rotation part
605    // ==================
606
607    /// Checks if this class contains an active associated spin rotation (normal or inverse).
608    ///
609    /// # Returns
610    ///
611    /// A flag indicating if this class contains an active associated spin rotation.
612    fn is_su2(&self) -> bool {
613        self.representative()
614            .as_ref()
615            .expect("No representative element found for this class.")
616            .is_su2()
617    }
618
619    /// Checks if this class contains an active and inverse associated spin rotation.
620    ///
621    /// # Returns
622    ///
623    /// A flag indicating if this class contains an active and inverse associated spin rotation.
624    fn is_su2_class_1(&self) -> bool {
625        self.representatives()
626            .expect("No representative element found for this class.")
627            .iter()
628            .all(|rep| rep.is_su2_class_1())
629    }
630}
631
632impl<R: FiniteOrder + Clone + Serialize> FiniteOrder for SymmetryClassSymbol<R> {
633    type Int = R::Int;
634
635    fn order(&self) -> Self::Int {
636        self.representative()
637            .as_ref()
638            .expect("No representative element found for this class.")
639            .order()
640    }
641}
642
643impl<R: Clone + Serialize> fmt::Display for SymmetryClassSymbol<R> {
644    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
645        write!(f, "{}", self.generic_symbol)
646    }
647}
648
649// =======
650// Methods
651// =======
652
653/// Reorder the rows so that the characters are in increasing order as we
654/// go down the table, with the first column being used as the primary sort
655/// order, then the second column, and so on, with the exceptions of some
656/// special columns, if available. These are:
657///
658/// * time-reversal class, $`\theta`$,
659/// * inversion class, $`i`$, or horizonal mirror plane $`\sigma_h`$ if $`i`$ not available.
660///
661/// # Arguments
662///
663/// * `char_arr` - A view of a two-dimensional square array containing the characters where
664/// each column is for one conjugacy class and each row one irrep.
665/// * `class_symbols` - An index map containing the conjugacy class symbols for the columns of
666/// `char_arr`. The keys are the symbols and the values are the column indices.
667///
668/// # Returns
669///
670/// A new array with the rows correctly sorted.
671///
672/// # Panics
673///
674/// Panics when expected classes cannot be found.
675pub(super) fn sort_irreps<R: Clone + Serialize + SpecialSymmetryTransformation>(
676    char_arr: &ArrayView2<Character>,
677    frobenius_schur_indicators: &[i8],
678    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
679    principal_classes: &[SymmetryClassSymbol<R>],
680) -> (Array2<Character>, Vec<i8>) {
681    log::debug!("Sorting irreducible representations...");
682    let su2_0 = if class_symbols.keys().any(|cc_sym| cc_sym.is_su2()) {
683        "(Σ)"
684    } else {
685        ""
686    };
687    let class_e = class_symbols
688        .first()
689        .expect("No class symbols found.")
690        .0
691        .clone();
692    let class_e1: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||E(QΣ)||", None)
693        .expect("Unable to construct class symbol `1||E(QΣ)||`.");
694    let class_i: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||i{su2_0}||"), None)
695        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||i{su2_0}||`."));
696    let class_ti: SymmetryClassSymbol<R> =
697        SymmetryClassSymbol::new(&format!("1||θ·i{su2_0}||"), None)
698            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·i{su2_0}||`."));
699    let class_s: SymmetryClassSymbol<R> =
700        SymmetryClassSymbol::new(&format!("1||σh{su2_0}||"), None)
701            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh{su2_0}||`."));
702    let class_s2: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||σh(Σ), σh(QΣ)||", None)
703        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh(Σ), σh(QΣ)||`."));
704    let class_ts: SymmetryClassSymbol<R> =
705        SymmetryClassSymbol::new(&format!("1||θ·σh{su2_0}||"), None)
706            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·σh{su2_0}||`."));
707    let class_t: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||θ{su2_0}||"), None)
708        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ{su2_0}||`."));
709
710    let mut leading_classes: IndexSet<SymmetryClassSymbol<R>> = IndexSet::new();
711    let mut sign_only_classes: HashSet<SymmetryClassSymbol<R>> = HashSet::new();
712
713    // Highest priority: SU(2) class 1
714    if class_symbols.contains_key(&class_e1) {
715        leading_classes.insert(class_e1.clone());
716        sign_only_classes.insert(class_e1);
717    }
718
719    // Second highest priority: time-reversal
720    if class_symbols.contains_key(&class_t) {
721        leading_classes.insert(class_t.clone());
722        sign_only_classes.insert(class_t);
723    }
724
725    // Third highest priority: inversion, or horizontal mirror plane if inversion not available,
726    // or time-reversed horizontal mirror plane if non-time-reversed version not available.
727    if class_symbols.contains_key(&class_i) {
728        leading_classes.insert(class_i.clone());
729        sign_only_classes.insert(class_i);
730    } else if class_symbols.contains_key(&class_ti) {
731        leading_classes.insert(class_ti.clone());
732        sign_only_classes.insert(class_ti);
733    } else if class_symbols.contains_key(&class_s) {
734        leading_classes.insert(class_s.clone());
735        sign_only_classes.insert(class_s);
736    } else if class_symbols.contains_key(&class_s2) {
737        leading_classes.insert(class_s2.clone());
738        sign_only_classes.insert(class_s2);
739    } else if class_symbols.contains_key(&class_ts) {
740        leading_classes.insert(class_ts.clone());
741        sign_only_classes.insert(class_ts);
742    };
743
744    // Forth highest priority: identity
745    leading_classes.insert(class_e);
746
747    // Fifth highest priority: principal classes, if not yet encountered
748    leading_classes.extend(principal_classes.iter().cloned());
749
750    log::debug!("Irreducible representation sort order:");
751    for leading_cc in leading_classes.iter() {
752        log::debug!(
753            "  {}{}",
754            leading_cc,
755            if sign_only_classes.contains(leading_cc) {
756                " (sign only)"
757            } else {
758                ""
759            }
760        );
761    }
762
763    let leading_idxs: IndexSet<usize> = leading_classes
764        .iter()
765        .map(|cc| {
766            *class_symbols
767                .get(cc)
768                .unwrap_or_else(|| panic!("Class `{cc}` cannot be found."))
769        })
770        .collect();
771
772    let n_rows = char_arr.nrows();
773    let mut col_idxs: Vec<usize> = Vec::with_capacity(n_rows);
774    col_idxs.extend(leading_idxs.iter());
775    col_idxs.extend((1..n_rows).filter(|i| !leading_idxs.contains(i)));
776    let mut sort_arr = char_arr.select(Axis(1), &col_idxs);
777    let one = Character::new(&[(UnityRoot::new(0, 2), 1)]);
778    let m_one = Character::new(&[(UnityRoot::new(1, 2), 1)]);
779    sign_only_classes.iter().for_each(|class| {
780        let col_idx = leading_classes
781            .get_index_of(class)
782            .unwrap_or_else(|| panic!("Unable to obtain the column index of class `{class}`."));
783        sort_arr.column_mut(col_idx).mapv_inplace(|character| {
784            let character_c = character.complex_value();
785            if approx::relative_eq!(
786                character_c.im,
787                0.0,
788                max_relative = character.threshold(),
789                epsilon = character.threshold()
790            ) {
791                if character_c.re > 0.0 {
792                    one.clone()
793                } else {
794                    m_one.clone()
795                }
796            } else if approx::relative_eq!(
797                character_c.re,
798                0.0,
799                max_relative = character.threshold(),
800                epsilon = character.threshold()
801            ) {
802                if character_c.im > 0.0 {
803                    one.clone()
804                } else {
805                    m_one.clone()
806                }
807            } else {
808                panic!("Character {character} is neither purely real nor purely imaginary for sign-only sorting.")
809            }
810        });
811    });
812    // let sort_arr = if has_e1 {
813    //     // Maps E(QΣ) characters to ±1, otherwise all linear/projective irreps of the same
814    //     // degeneracy would be grouped together.
815    //     let mut sort_arr = char_arr.select(Axis(1), &col_idxs);
816    //     let one = Character::new(&[(UnityRoot::new(0, 2), 1)]);
817    //     let m_one = Character::new(&[(UnityRoot::new(1, 2), 1)]);
818
819    //     sort_arr.column_mut(0).mapv_inplace(|e1_character| {
820    //         if e1_character.complex_value().re > 0.0 {
821    //             one.clone()
822    //         } else {
823    //             m_one.clone()
824    //         }
825    //     });
826    //     sort_arr
827    // } else {
828    //     char_arr.select(Axis(1), &col_idxs)
829    // };
830
831    let sort_row_indices: Vec<_> = (0..n_rows)
832        .sorted_by(|&i, &j| {
833            let keys_i = sort_arr.row(i).iter().cloned().collect_vec();
834            let keys_j = sort_arr.row(j).iter().cloned().collect_vec();
835            keys_i
836                .partial_cmp(&keys_j)
837                .unwrap_or_else(|| panic!("`{keys_i:?}` and `{keys_j:?}` cannot be compared."))
838        })
839        .collect();
840    let char_arr = char_arr.select(Axis(0), &sort_row_indices);
841    let old_fs = frobenius_schur_indicators.iter().collect::<Vec<_>>();
842    let sorted_fs = sort_row_indices.iter().map(|&i| *old_fs[i]).collect_vec();
843    log::debug!("Sorting irreducible representations... Done.");
844    (char_arr, sorted_fs)
845}
846
847/// Determines the principal classes given a list of class symbols and any forcing conditions.
848///
849/// By default, the principal classes are those with the highest order, regardless of whether they
850/// are proper or improper, unitary or antiunitary.
851///
852/// # Arguments
853///
854/// * `class_symbols` - An indexmap of class symbols and their corresponding indices.
855/// * `force_proper_principal` - A flag indicating if the principal classes are forced to be
856/// proper.
857/// * `force_principal` - An option containing specific classes that are forced to be principal.
858///
859/// # Returns
860///
861/// A vector of symbols of principal classes.
862///
863/// # Panics
864///
865/// Panics when:
866///
867/// * both `force_proper_principal` and `force_principal` are specified;
868/// * classes specified in `force_principal` cannot be found in `class_symbols`;
869/// * no principal classes can be found.
870pub(super) fn deduce_principal_classes<R, P>(
871    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
872    force_principal_predicate: Option<P>,
873    force_principal: Option<SymmetryClassSymbol<R>>,
874) -> Vec<SymmetryClassSymbol<R>>
875where
876    R: fmt::Debug + SpecialSymmetryTransformation + FiniteOrder + Clone + Serialize,
877    P: Copy + Fn(&SymmetryClassSymbol<R>) -> bool,
878{
879    log::debug!("Determining principal classes...");
880    let principal_classes = force_principal.map_or_else(|| {
881        // sorted_class_symbols contains class symbols sorted in the order:
882        // - decreasing operation order
883        // - unitary operations, then antiunitary operations
884        // - proper operations, then improper operations in each operation order
885        // - if SU(2), then homotopy class 0, then homotopy class 1
886        let mut sorted_class_symbols = class_symbols
887            .clone()
888            .sorted_by(|cc1, _, cc2, _| {
889                PartialOrd::partial_cmp(
890                    &(cc1.order(), !cc1.contains_time_reversal(), cc1.is_proper(), !cc1.is_su2_class_1()),
891                    &(cc2.order(), !cc2.contains_time_reversal(), cc2.is_proper(), !cc2.is_su2_class_1()),
892                )
893                .expect("Unable to sort class symbols.")
894            })
895            .rev();
896        let (principal_rep, _) = if let Some(predicate) = force_principal_predicate {
897            // Larger order always prioritised, regardless of proper or improper,
898            // unless force_proper_principal is set, in which case
899            sorted_class_symbols
900                .find(|(cc, _)| predicate(cc))
901                .expect("`Unable to find classes fulfilling the specified predicate.`")
902        } else {
903            // No force_proper_principal, so proper prioritised, which has been ensured by the sort
904            // in the construction of `sorted_class_symbols`.
905            sorted_class_symbols
906                .next()
907                .expect("Unexpected empty `sorted_class_symbols`.")
908        };
909
910        // Now find all principal classes with the same order, parity, unitarity, and homotopy
911        // class as `principal_rep`.
912        class_symbols
913            .keys()
914            .filter(|cc| {
915                cc.contains_time_reversal() == principal_rep.contains_time_reversal()
916                && cc.is_proper() == principal_rep.is_proper()
917                && cc.order() == principal_rep.order()
918                && cc.is_su2_class_1() == principal_rep.is_su2_class_1()
919            })
920            .cloned()
921            .collect::<Vec<_>>()
922    }, |principal_cc| {
923        assert!(
924            force_principal_predicate.is_none(),
925            "`force_principal_predicate` and `force_principal` cannot be both provided."
926        );
927        assert!(
928            class_symbols.contains_key(&principal_cc),
929            "Forcing principal-axis class to be `{principal_cc}`, but `{principal_cc}` is not a valid class in the group."
930        );
931
932        log::warn!(
933            "Principal-axis class forced to be `{principal_cc}`. Auto-detection of principal-axis classes will be skipped.",
934        );
935        vec![principal_cc]
936    });
937
938    assert!(!principal_classes.is_empty());
939    if principal_classes.len() == 1 {
940        log::debug!("Principal-axis class found: {}", principal_classes[0]);
941    } else {
942        log::debug!("Principal-axis classes found:");
943        for princc in &principal_classes {
944            log::debug!("  {}", princc);
945        }
946    }
947    log::debug!("Determining principal classes... Done.");
948    principal_classes
949}
950
951/// Deduces irreducible representation symbols based on Mulliken's convention.
952///
953/// # Arguments
954///
955/// * `char_arr` - A view of a two-dimensional square array containing the characters where
956/// each column is for one conjugacy class and each row one irrep.
957/// * `class_symbols` - An index map containing the conjugacy class symbols for the columns of
958/// `char_arr`. The keys are the symbols and the values are the column indices.
959/// * `principal_classes` - The principal classes to be used for irrep symbol deduction.
960///
961/// # Returns
962///
963/// A vector of Mulliken symbols corresponding to the rows of `char_arr`.
964///
965/// # Panics
966///
967/// Panics when expected classes cannot be found in `class_symbols`.
968#[allow(clippy::too_many_lines)]
969pub(super) fn deduce_mulliken_irrep_symbols<R>(
970    char_arr: &ArrayView2<Character>,
971    class_symbols: &IndexMap<SymmetryClassSymbol<R>, usize>,
972    principal_classes: &[SymmetryClassSymbol<R>],
973) -> Vec<MullikenIrrepSymbol>
974where
975    R: fmt::Debug + SpecialSymmetryTransformation + FiniteOrder + Clone + Serialize,
976{
977    log::debug!("Generating Mulliken irreducible representation symbols...");
978
979    let su2_0 = if class_symbols.keys().any(|cc_sym| cc_sym.is_su2()) {
980        "(Σ)"
981    } else {
982        ""
983    };
984
985    let e_cc = class_symbols
986        .first()
987        .expect("No class symbols found.")
988        .0
989        .clone();
990    let e1_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||E(QΣ)||", None)
991        .expect("Unable to construct class symbol `1||E(QΣ)||`.");
992    let i_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||i{su2_0}||"), None)
993        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||i{su2_0}||`."));
994    let ti_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 s_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||σh{su2_0}||"), None)
997        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh{su2_0}||`."));
998    let s2_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new("1||σh(Σ), σh(QΣ)||", None)
999        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||σh(Σ), σh(QΣ)||`."));
1000    let ts_cc: SymmetryClassSymbol<R> =
1001        SymmetryClassSymbol::new(&format!("1||θ·σh{su2_0}||"), None)
1002            .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ·σh{su2_0}||`."));
1003    let t_cc: SymmetryClassSymbol<R> = SymmetryClassSymbol::new(&format!("1||θ{su2_0}||"), None)
1004        .unwrap_or_else(|_| panic!("Unable to construct class symbol `1||θ{su2_0}||`."));
1005
1006    // SU(2) class 1 parity?
1007    let su2_1_parity = class_symbols.contains_key(&e1_cc);
1008
1009    // Inversion parity?
1010    let i_parity = class_symbols.contains_key(&i_cc);
1011
1012    // Time-reversed inversion parity?
1013    let ti_parity = class_symbols.contains_key(&ti_cc);
1014
1015    // Reflection parity?
1016    let s_parity = class_symbols.contains_key(&s_cc) || class_symbols.contains_key(&s2_cc);
1017
1018    // Time-reversed reflection parity?
1019    let ts_parity = class_symbols.contains_key(&ts_cc);
1020
1021    // SU(2) class 1?
1022    if su2_1_parity {
1023        log::debug!("E(QΣ) found. This will be used to determine projective representations.");
1024    }
1025
1026    // i_parity takes priority.
1027    if i_parity {
1028        log::debug!("Inversion centre found. This will be used for g/u ordering.");
1029    } else if ti_parity {
1030        log::debug!(
1031            "Time-reversed inversion centre found (but no inversion centre). This will be used for g/u ordering."
1032        );
1033    } else if s_parity {
1034        log::debug!(
1035            "Horizontal mirror plane found (but no inversion centre). This will be used for '/'' ordering."
1036        );
1037    } else if ts_parity {
1038        log::debug!(
1039            "Time-reversed horizontal mirror plane found (but no inversion centre). This will be used for '/'' ordering."
1040        );
1041    }
1042
1043    // Time-reversal?
1044    let t_parity = class_symbols.contains_key(&t_cc);
1045    if t_parity {
1046        log::debug!("Time reversal found. This will be used for magnetic ordering.");
1047    }
1048
1049    let e2p1 = UnityRoot::new(1u32, 2u32);
1050    let e2p2 = UnityRoot::new(2u32, 2u32);
1051    let char_p1 = Character::new(&[(e2p2, 1usize)]);
1052    let char_m1 = Character::new(&[(e2p1, 1usize)]);
1053
1054    // First pass: assign irrep symbols based on Mulliken's convention as much as possible.
1055    log::debug!("First pass: assign symbols from rules");
1056
1057    let mut raw_irrep_symbols = char_arr.rows().into_iter().map(|irrep| {
1058        // Determine the main symmetry
1059        let dim = irrep[
1060            *class_symbols.get(&e_cc).unwrap_or_else(|| panic!("Class `{}` not found.", &e_cc))
1061        ].complex_value();
1062        assert!(
1063            approx::relative_eq!(dim.im, 0.0)
1064                && approx::relative_eq!(dim.re.round(), dim.re)
1065                && dim.re.round() > 0.0
1066        );
1067
1068        assert!(dim.re.round() > 0.0);
1069        #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
1070        let dim = dim.re.round() as u64;
1071        assert!(dim > 0);
1072        let main = if dim >= 2 {
1073            // Degenerate irreps
1074            INV_MULLIKEN_IRREP_DEGENERACIES.get(&dim).map_or_else(|| {
1075                log::warn!("{} cannot be assigned a standard dimensionality symbol. A generic 'Λ' will be used instead.", dim);
1076                "Λ"
1077            }, |sym| sym)
1078        } else {
1079            // Non-degenerate irreps
1080            let complex = irrep.map(|character| character.complex_conjugate()) != irrep;
1081            if complex {
1082                "Γ"
1083            } else {
1084                let char_rots: HashSet<_> = principal_classes
1085                    .iter()
1086                    .map(|cc| {
1087                        irrep[
1088                            *class_symbols.get(cc).unwrap_or_else(|| panic!("Class `{cc}` not found."))
1089                        ].clone()
1090                    })
1091                    .collect();
1092                if char_rots.len() == 1 && *char_rots.iter().next().expect("No rotation classes found.") == char_p1 {
1093                    "A"
1094                } else if char_rots
1095                    .iter()
1096                    .all(|char_rot| char_rot.clone() == char_p1 || char_rot.clone() == char_m1)
1097                {
1098                    "B"
1099                } else {
1100                    panic!("");
1101                    // // There are principal rotations but with non-(±1) characters. These must be
1102                    // // complex.
1103                    // "Γ"
1104                }
1105            }
1106        };
1107
1108        let projective = if su2_1_parity {
1109            let char_e1 = irrep[
1110                *class_symbols
1111                    .get(&e1_cc)
1112                    .unwrap_or_else(|| {
1113                        panic!("The character for `{e1_cc}` could not be found.")
1114                    })
1115            ].clone();
1116            let char_e1_c = char_e1.complex_value();
1117            assert!(
1118                approx::relative_eq!(
1119                    char_e1_c.im,
1120                    0.0,
1121                    epsilon = char_e1.threshold(),
1122                    max_relative = char_e1.threshold()
1123                ) && approx::relative_eq!(
1124                    char_e1_c.re.round(),
1125                    char_e1_c.re,
1126                    epsilon = char_e1.threshold(),
1127                    max_relative = char_e1.threshold()
1128                ),
1129            );
1130
1131            #[allow(clippy::cast_possible_truncation)]
1132            let char_e1_c = char_e1_c
1133                .re
1134                .round()
1135                .to_i32()
1136                .unwrap_or_else(|| panic!("Unable to convert the real part of the character for `{e1_cc}` to `i32`."));
1137            let dim_i32 = i32::try_from(dim).expect("Unable to convert the irrep dimensionality to `i32`.");
1138            if char_e1_c == dim_i32 {
1139                false
1140            } else {
1141                assert_eq!(char_e1_c, -dim_i32);
1142                true
1143            }
1144        } else {
1145            false
1146        };
1147        let projective_str = if projective { "~" } else { "" };
1148
1149        let (inv, mir) = if i_parity || ti_parity {
1150            // Determine inversion symmetry
1151            // Inversion symmetry trumps reflection symmetry.
1152            let char_inv = irrep[
1153                *class_symbols
1154                    .get(&i_cc)
1155                    .unwrap_or_else(|| {
1156                        class_symbols
1157                            .get(&ti_cc)
1158                            .unwrap_or_else(|| {
1159                                panic!("Neither `{}` nor `{}` found.", &i_cc, &ti_cc)
1160                            })
1161                    })
1162            ].clone();
1163            let char_inv_c = char_inv.complex_value();
1164            assert!(
1165                approx::relative_eq!(
1166                    char_inv_c.im,
1167                    0.0,
1168                    epsilon = char_inv.threshold(),
1169                    max_relative = char_inv.threshold()
1170                ) && approx::relative_eq!(
1171                    char_inv_c.re.round(),
1172                    char_inv_c.re,
1173                    epsilon = char_inv.threshold(),
1174                    max_relative = char_inv.threshold()
1175                ),
1176            );
1177
1178            #[allow(clippy::cast_possible_truncation)]
1179            let char_inv_c = char_inv_c.re.round() as i32;
1180            match char_inv_c.cmp(&0) {
1181                Ordering::Greater => ("g", ""),
1182                Ordering::Less => ("u", ""),
1183                Ordering::Equal => panic!("Inversion character must not be zero."),
1184            }
1185        } else if !projective && (s_parity || ts_parity) {
1186            // Determine reflection symmetry
1187            // Note that for projective irreps, characters under σh are either imaginary or zero,
1188            // thus no '/'' classifications.
1189            let char_ref = irrep[
1190                *class_symbols
1191                    .get(&s_cc)
1192                    .unwrap_or_else(|| {
1193                        class_symbols
1194                            .get(&s2_cc)
1195                            .unwrap_or_else(|| {
1196                                class_symbols
1197                                    .get(&ts_cc)
1198                                    .unwrap_or_else(|| {
1199                                        panic!("Neither `{}`, `{}`, nor `{}` found.", &s_cc, &s2_cc, &ts_cc)
1200                                        })
1201                        })
1202                    })
1203            ].clone();
1204            let char_ref_c = char_ref.complex_value();
1205            assert!(
1206                approx::relative_eq!(
1207                    char_ref_c.im,
1208                    0.0,
1209                    epsilon = char_ref.threshold(),
1210                    max_relative = char_ref.threshold()
1211                ) && approx::relative_eq!(
1212                    char_ref_c.re.round(),
1213                    char_ref_c.re,
1214                    epsilon = char_ref.threshold(),
1215                    max_relative = char_ref.threshold()
1216                ),
1217            );
1218
1219            #[allow(clippy::cast_possible_truncation)]
1220            let char_ref_c = char_ref_c.re.round() as i32;
1221            match char_ref_c.cmp(&0) {
1222                Ordering::Greater => ("", "'"),
1223                Ordering::Less => ("", "''"),
1224                Ordering::Equal => panic!("Reflection character must not be zero for linear irreducible representations."),
1225            }
1226        } else {
1227            ("", "")
1228        };
1229
1230
1231        let trev = if t_parity {
1232            // Determine time-reversal symmetry
1233            let char_trev = irrep[
1234                *class_symbols.get(&t_cc).unwrap_or_else(|| {
1235                    panic!("Class `{}` not found.", &t_cc)
1236                })
1237            ].clone();
1238            let char_trev_c = char_trev.complex_value();
1239            if approx::relative_eq!(
1240                char_trev_c.im,
1241                0.0,
1242                epsilon = char_trev.threshold(),
1243                max_relative = char_trev.threshold()
1244            ) && approx::relative_eq!(
1245                char_trev_c.re.round(),
1246                char_trev_c.re,
1247                epsilon = char_trev.threshold(),
1248                max_relative = char_trev.threshold()
1249            ) {
1250                // Real, integral time-reversal character
1251                #[allow(clippy::cast_possible_truncation)]
1252                let char_trev_c = char_trev_c.re.round() as i32;
1253                match char_trev_c.cmp(&0) {
1254                    Ordering::Greater => "+",
1255                    Ordering::Less => "-",
1256                    Ordering::Equal => panic!("Real time-reversal character must not be zero."),
1257                }
1258            } else {
1259                // Non-real or non-integral time-reversal character
1260                ""
1261            }
1262        } else {
1263            ""
1264        };
1265
1266        MullikenIrrepSymbol::new(format!("|^({trev})|{main}{projective_str}|^({mir})_({inv})|").as_str())
1267            .unwrap_or_else(|_| {
1268                panic!(
1269                    "Unable to construct symmetry symbol `|^({trev})|{main}{projective_str}|^({mir})_({inv})|`."
1270                )
1271            })
1272    }).collect_vec();
1273
1274    let mut complex_irrep_indices = char_arr
1275        .rows()
1276        .into_iter()
1277        .enumerate()
1278        .filter_map(|(i, irrep)| {
1279            let complex = irrep.map(|character| character.complex_conjugate()) != irrep;
1280            if complex {
1281                Some(i)
1282            } else {
1283                None
1284            }
1285        })
1286        .collect::<IndexSet<_>>();
1287    complex_irrep_indices.reverse();
1288
1289    let cc_pairs = if !complex_irrep_indices.is_empty() {
1290        log::debug!("Grouping pairs of complex-conjugate irreps...");
1291        let mut cc_pairs: Vec<(usize, usize)> = vec![];
1292        while !complex_irrep_indices.is_empty() {
1293            let complex_irrep_index = complex_irrep_indices.pop().unwrap();
1294            let complex_irrep = char_arr.row(complex_irrep_index);
1295            let complex_conj_irrep = complex_irrep.map(|character| character.complex_conjugate());
1296            let complex_conj_irrep_index = char_arr
1297                .rows()
1298                .into_iter()
1299                .enumerate()
1300                .find_map(|(i, irrep)| {
1301                    if irrep == complex_conj_irrep {
1302                        Some(i)
1303                    } else {
1304                        None
1305                    }
1306                })
1307                .expect("Unable to find the complex-conjugate irrep.");
1308            complex_irrep_indices.shift_remove(&complex_conj_irrep_index);
1309            cc_pairs.push((complex_irrep_index, complex_conj_irrep_index));
1310            raw_irrep_symbols[complex_irrep_index]
1311                .generic_symbol
1312                .set_presub("a");
1313            raw_irrep_symbols[complex_conj_irrep_index]
1314                .generic_symbol
1315                .set_presub("b");
1316        }
1317        log::debug!(
1318            "There {} {} {} of complex-conjugate irreps.",
1319            if cc_pairs.len() == 1 { "is" } else { "are" },
1320            cc_pairs.len(),
1321            if cc_pairs.len() == 1 { "pair" } else { "pairs" },
1322        );
1323        log::debug!("Grouping pairs of complex-conjugate irreps... Done.");
1324        Some(cc_pairs)
1325    } else {
1326        None
1327    };
1328
1329    log::debug!("Second pass: disambiguate identical cases not distinguishable by rules");
1330    let mut irrep_symbols = disambiguate_linspace_symbols(raw_irrep_symbols.into_iter());
1331
1332    if let Some(cc_pairs_vec) = cc_pairs {
1333        log::debug!("Equalising post-subscripts in complex-conjugate pairs...");
1334        for (complex_irrep_index, complex_conj_irrep_index) in &cc_pairs_vec {
1335            let complex_irrep_postsub = irrep_symbols[*complex_irrep_index].postsub();
1336            irrep_symbols[*complex_conj_irrep_index]
1337                .generic_symbol
1338                .set_postsub(&complex_irrep_postsub);
1339        }
1340        log::debug!("Equalising post-subscripts in complex-conjugate pairs... Done.");
1341    }
1342
1343    log::debug!("Generating Mulliken irreducible representation symbols... Done.");
1344
1345    irrep_symbols
1346}
1347
1348/// Determines the mirror-plane symbol given a principal axis.
1349///
1350/// # Arguments
1351///
1352/// * `sigma_axis` - The normalised normal vector of a mirror plane.
1353/// * `principal_axis` - The normalised principal rotation axis.
1354/// * `thresh` - Threshold for comparisons.
1355/// * `force_d` - Flag indicating if vertical mirror planes should be given the $`d`$ symbol
1356/// instead of $`v`$.
1357///
1358/// # Returns
1359///
1360/// The mirror-plane symbol.
1361#[must_use]
1362pub(super) fn deduce_sigma_symbol(
1363    sigma_axis: &Vector3<f64>,
1364    principal_element: &SymmetryElement,
1365    thresh: f64,
1366    force_d: bool,
1367) -> Option<String> {
1368    if approx::relative_eq!(
1369        principal_element.raw_axis().dot(sigma_axis).abs(),
1370        0.0,
1371        epsilon = thresh,
1372        max_relative = thresh
1373    ) && *principal_element.raw_proper_order() != ORDER_1
1374    {
1375        // Vertical plane containing principal axis
1376        if force_d {
1377            Some("d".to_owned())
1378        } else {
1379            Some("v".to_owned())
1380        }
1381    } else if approx::relative_eq!(
1382        principal_element.raw_axis().cross(sigma_axis).norm(),
1383        0.0,
1384        epsilon = thresh,
1385        max_relative = thresh
1386    ) && *principal_element.raw_proper_order() != ORDER_1
1387    {
1388        // Horizontal plane perpendicular to principal axis
1389        Some("h".to_owned())
1390    } else {
1391        None
1392    }
1393}
1394
1395/// Enumerated type specifying the parity under a mirror plane.
1396#[derive(Clone, PartialEq, Eq, Hash, Debug)]
1397pub enum MirrorParity {
1398    /// Variant for even parity.
1399    Even,
1400
1401    /// Variant for odd parity.
1402    Odd,
1403
1404    /// Variant for no parity.
1405    Neither,
1406}
1407
1408impl fmt::Display for MirrorParity {
1409    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1410        match self {
1411            MirrorParity::Even => write!(f, "(+)"),
1412            MirrorParity::Odd => write!(f, "(-)"),
1413            MirrorParity::Neither => write!(f, "( )"),
1414        }
1415    }
1416}
1417
1418/// Given a group with a character table, deduces the parities of a specified representation under
1419/// all classes in the group that contain a spatial reflection.
1420///
1421/// # Arguments
1422///
1423/// * `group` - A group with a character table.
1424/// * `rep` - A representation (that has been decomposed as irreps or ircoreps of `group`) for
1425/// which mirror parities are to be deduced.
1426///
1427/// # Returns
1428///
1429/// An indexmap whose keys are reflection classes and whose values are the corresponding parities of
1430/// `rep`.
1431///
1432/// # Panics
1433///
1434/// Panics on unexpected errors.
1435pub(crate) fn deduce_mirror_parities<G, R>(
1436    group: &G,
1437    rep: &R,
1438) -> IndexMap<<<G as CharacterProperties>::CharTab as CharacterTable>::ColSymbol, MirrorParity>
1439where
1440    G: SymmetryGroupProperties,
1441    R: ReducibleLinearSpaceSymbol<
1442        Subspace = <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
1443    >,
1444{
1445    let id_cc = group
1446        .get_cc_symbol_of_index(0)
1447        .expect("Unable to obtain the conjugacy class symbol of the identity.");
1448    let mirrors = group.filter_cc_symbols(|cc| cc.is_spatial_reflection());
1449    mirrors
1450        .iter()
1451        .map(|sigma_cc| {
1452            let sigma_cc_chars = rep
1453                .subspaces()
1454                .iter()
1455                .map(|(irrep, _)| {
1456                    let id_char = group.character_table().get_character(irrep, &id_cc);
1457                    let char = group.character_table().get_character(irrep, sigma_cc);
1458                    if *char == *id_char {
1459                        MirrorParity::Even
1460                    } else if *char == -id_char {
1461                        MirrorParity::Odd
1462                    } else {
1463                        MirrorParity::Neither
1464                    }
1465                })
1466                .collect::<HashSet<_>>();
1467            let sigma_cc_parity = if sigma_cc_chars.len() == 1 {
1468                sigma_cc_chars
1469                    .into_iter()
1470                    .next()
1471                    .expect("Unable to extract the mirror parity.")
1472            } else {
1473                MirrorParity::Neither
1474            };
1475            (sigma_cc.clone(), sigma_cc_parity)
1476        })
1477        .collect::<IndexMap<_, _>>()
1478}