qsym2/analysis/
mod.rs

1//! Symmetry analysis via representation and corepresentation theories.
2
3use std::cmp::Ordering;
4use std::fmt;
5use std::ops::Mul;
6
7use anyhow::{self, format_err, Context};
8use duplicate::duplicate_item;
9use itertools::Itertools;
10use log;
11use ndarray::{s, Array, Array1, Array2, Axis, Dimension, Ix0, Ix2};
12use ndarray_einsum::*;
13use ndarray_linalg::{solve::Inverse, types::Lapack};
14use num_complex::{Complex, ComplexFloat};
15use num_traits::{ToPrimitive, Zero};
16#[cfg(feature = "python")]
17use pyo3::prelude::*;
18use rayon::prelude::*;
19use serde::{Deserialize, Serialize};
20
21use crate::chartab::chartab_group::CharacterProperties;
22use crate::chartab::{CharacterTable, DecompositionError, SubspaceDecomposable};
23use crate::group::{class::ClassProperties, GroupProperties};
24use crate::io::format::{log_subtitle, qsym2_output};
25use crate::symmetry::symmetry_group::UnitaryRepresentedSymmetryGroup;
26
27// ======
28// Metric
29// ======
30
31/// Structure to handle various choices for the specification of the metric matrices.
32#[derive(Clone)]
33pub(crate) struct Metric<'a, T, D: Dimension> {
34    /// The  Hermitian atomic-orbital overlap matrix of the underlying basis set used to describe
35    /// the determinant. For restricted, unrestricted, and generalised spin constraints, only the
36    /// spatial overlap matrix is required for the full overlap matrix to be deduced. For
37    /// relativistic SOC spin constraint, the full overlap matrix must be specified.
38    hermitian: &'a Array<T, D>,
39
40    /// The full (w.r.t. the pertinent spin constraint) complex-symmetric atomic-orbital overlap
41    /// matrix of the underlying basis set used to describe the determinant. For restricted,
42    /// unrestricted, and generalised spin constraints, only the spatial overlap matrix is required
43    /// for the full overlap matrix to be deduced. For relativistic SOC spin constraint, the full
44    /// overlap matrix must be specified.
45    complex_symmetric: Option<&'a Array<T, D>>,
46}
47
48// =======
49// Overlap
50// =======
51
52// ----------------
53// Trait definition
54// ----------------
55
56/// Trait for computing the inner product
57/// $`\langle \hat{\iota} \mathbf{v}_i, \mathbf{v}_j \rangle`$ between two linear-space quantities
58/// $`\mathbf{v}_i`$ and $`\mathbf{v}_j`$. The involutory operator $`\hat{\iota}`$ determines
59/// whether the inner product is a sesquilinear form or a bilinear form.
60pub trait Overlap<T, D>
61where
62    T: ComplexFloat + fmt::Debug + Lapack,
63    D: Dimension,
64{
65    /// If `true`, the inner product is bilinear and $`\hat{\iota} = \hat{\kappa}`$. If `false`,
66    /// the inner product is sesquilinear and $`\hat{\iota} = \mathrm{id}`$.
67    fn complex_symmetric(&self) -> bool;
68
69    /// Returns the overlap between `self` and `other`, with respect to a metric `metric` (and
70    /// possibly its complex-symmetric version `metric_h`) of the underlying basis in which `self`
71    /// and `other` are expressed.
72    fn overlap(
73        &self,
74        other: &Self,
75        metric: Option<&Array<T, D>>,
76        metric_h: Option<&Array<T, D>>,
77    ) -> Result<T, anyhow::Error>;
78
79    /// Returns a string describing the mathematical definition of the overlap (*i.e.* the inner
80    /// product) between two quantities of this type.
81    fn overlap_definition(&self) -> String;
82}
83
84// =====
85// Orbit
86// =====
87
88// --------------------------------------
89// Struct definitions and implementations
90// --------------------------------------
91
92/// Lazy iterator for orbits generated by the action of a group on an origin.
93pub struct OrbitIterator<'a, G, I>
94where
95    G: GroupProperties,
96{
97    /// A mutable iterator over the elements of the group. Each element will be applied on the
98    /// origin to yield a corresponding item in the orbit.
99    group_iter: <<G as GroupProperties>::ElementCollection as IntoIterator>::IntoIter,
100
101    /// The origin of the orbit.
102    origin: &'a I,
103
104    /// A function defining the action of each group element on the origin.
105    action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
106}
107
108impl<'a, G, I> OrbitIterator<'a, G, I>
109where
110    G: GroupProperties,
111{
112    /// Creates and returns a new orbit iterator.
113    ///
114    /// # Arguments
115    ///
116    /// * `group` - A group.
117    /// * `origin` - An origin.
118    /// * `action` - A function or closure defining the action of each group element on the origin.
119    ///
120    /// # Returns
121    ///
122    /// An orbit iterator.
123    pub fn new(
124        group: &G,
125        origin: &'a I,
126        action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
127    ) -> Self {
128        Self {
129            group_iter: group.elements().clone().into_iter(),
130            origin,
131            action,
132        }
133    }
134}
135
136impl<'a, G, I> Iterator for OrbitIterator<'a, G, I>
137where
138    G: GroupProperties,
139{
140    type Item = Result<I, anyhow::Error>;
141
142    fn next(&mut self) -> Option<Self::Item> {
143        self.group_iter
144            .next()
145            .map(|op| (self.action)(&op, self.origin))
146    }
147}
148
149// ----------------
150// Trait definition
151// ----------------
152
153/// Trait for orbits arising from group actions.
154pub trait Orbit<G, I>
155where
156    G: GroupProperties,
157{
158    /// Type of the iterator over items in the orbit.
159    type OrbitIter: Iterator<Item = Result<I, anyhow::Error>>;
160
161    /// The group generating the orbit.
162    fn group(&self) -> &G;
163
164    /// The origin of the orbit.
165    fn origin(&self) -> &I;
166
167    /// An iterator over items in the orbit arising from the action of the group on the origin.
168    fn iter(&self) -> Self::OrbitIter;
169}
170
171// ========
172// Analysis
173// ========
174
175// ---------------
176// Enum definition
177// ---------------
178
179/// Enumerated type specifying the comparison mode for filtering out orbit overlap
180/// eigenvalues.
181#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
182#[cfg_attr(feature = "python", pyclass(eq, eq_int))]
183pub enum EigenvalueComparisonMode {
184    /// Compares the eigenvalues using only their real parts.
185    Real,
186
187    /// Compares the eigenvalues using their moduli.
188    Modulus,
189}
190
191impl Default for EigenvalueComparisonMode {
192    fn default() -> Self {
193        Self::Modulus
194    }
195}
196
197impl fmt::Display for EigenvalueComparisonMode {
198    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
199        match self {
200            EigenvalueComparisonMode::Real => write!(f, "Real part"),
201            EigenvalueComparisonMode::Modulus => write!(f, "Modulus"),
202        }
203    }
204}
205
206// ----------------
207// Trait definition
208// ----------------
209
210// ~~~~~~~~~~~
211// RepAnalysis
212// ~~~~~~~~~~~
213
214/// Trait for representation or corepresentation analysis on an orbit of items spanning a
215/// linear space.
216pub trait RepAnalysis<G, I, T, D>: Orbit<G, I>
217where
218    T: ComplexFloat + Lapack + fmt::Debug,
219    <T as ComplexFloat>::Real: ToPrimitive,
220    G: GroupProperties + ClassProperties + CharacterProperties,
221    G::GroupElement: fmt::Display,
222    G::CharTab: SubspaceDecomposable<T>,
223    D: Dimension,
224    I: Overlap<T, D> + Clone,
225    Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
226{
227    // ----------------
228    // Required methods
229    // ----------------
230
231    /// Sets the overlap matrix between the items in the orbit.
232    ///
233    /// # Arguments
234    ///
235    /// * `smat` - The overlap matrix between the items in the orbit.
236    fn set_smat(&mut self, smat: Array2<T>);
237
238    /// Returns the overlap matrix between the items in the orbit.
239    #[must_use]
240    fn smat(&self) -> Option<&Array2<T>>;
241
242    /// Returns the transformation matrix $`\mathbf{X}`$ for the overlap matrix $`\mathbf{S}`$
243    /// between the items in the orbit.
244    ///
245    /// The matrix $`\mathbf{X}`$ serves to bring $`\mathbf{S}`$ to full rank, *i.e.*, the matrix
246    /// $`\tilde{\mathbf{S}}`$ defined by
247    ///
248    /// ```math
249    ///     \tilde{\mathbf{S}} = \mathbf{X}^{\dagger\lozenge} \mathbf{S} \mathbf{X}
250    /// ```
251    ///
252    /// is a full-rank matrix.
253    ///
254    /// If the overlap between items is complex-symmetric (see [`Overlap::complex_symmetric`]), then
255    /// $`\lozenge = *`$ is the complex-conjugation operation, otherwise, $`\lozenge`$ is the
256    /// identity.
257    ///
258    /// Depending on how $`\mathbf{X}`$ has been computed, $`\tilde{\mathbf{S}}`$ might also be
259    /// orthogonal. Either way, $`\tilde{\mathbf{S}}`$ is always guaranteed to be of full-rank.
260    #[must_use]
261    fn xmat(&self) -> &Array2<T>;
262
263    /// Returns the norm-preserving scalar map $`f`$ for every element of the generating group
264    /// defined by
265    ///
266    /// ```math
267    ///     \langle \hat{\iota} \mathbf{v}_w, \hat{g}_i \mathbf{v}_x \rangle
268    ///     = f \left( \langle \hat{\iota} \hat{g}_i^{-1} \mathbf{v}_w, \mathbf{v}_x \rangle \right).
269    /// ```
270    ///
271    /// Typically, for $`\hat{\iota} = \mathrm{id}`$, if $`\hat{g}_i`$ is unitary, then $`f`$ is the
272    /// identity, and if $`\hat{g}_i`$ is antiunitary, then $`f`$ is the complex-conjugation
273    /// operation. Either way, the norm of the inner product is preserved.
274    ///
275    /// If the overlap between items is complex-symmetric (*i.e.* $`\hat{\iota} = \hat{\kappa}`$,
276    /// see [`Overlap::complex_symmetric`]), then this map is currently unsupported because it is
277    /// currently unclear if the unitary and antiunitary symmetry operators in QSym² commute with
278    /// $`\hat{\kappa}`$. This thus precludes the use of the Cayley table to speed up the computation
279    /// of the orbit overlap matrix.
280    #[must_use]
281    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error>;
282
283    /// Returns the threshold for integrality checks of irreducible representation or
284    /// corepresentation multiplicities.
285    #[must_use]
286    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real;
287
288    /// Returns the enumerated type specifying the comparison mode for filtering out orbit overlap
289    /// eigenvalues.
290    #[must_use]
291    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode;
292
293    // ----------------
294    // Provided methods
295    // ----------------
296
297    /// Calculates and stores the overlap matrix between items in the orbit, with respect to a
298    /// metric of the basis in which these items are expressed.
299    ///
300    /// # Arguments
301    ///
302    /// * `metric` - The metric of the basis in which the orbit items are expressed.
303    /// * `metric_h` - The complex-symmetric metric of the basis in which the orbit items are
304    /// expressed. This is required if antiunitary operations are involved.
305    /// * `use_cayley_table` - A boolean indicating if the Cayley table of the group should be used
306    /// to speed up the computation of the overlap matrix.
307    fn calc_smat(
308        &mut self,
309        metric: Option<&Array<T, D>>,
310        metric_h: Option<&Array<T, D>>,
311        use_cayley_table: bool,
312    ) -> Result<&mut Self, anyhow::Error> {
313        let order = self.group().order();
314        let mut smat = Array2::<T>::zeros((order, order));
315        let item_0 = self.origin();
316        if let (Some(ctb), true) = (self.group().cayley_table(), use_cayley_table) {
317            log::debug!("Cayley table available. Group closure will be used to speed up overlap matrix computation.");
318            let ovs = self
319                .iter()
320                .map(|item_res| {
321                    let item = item_res?;
322                    item.overlap(item_0, metric, metric_h)
323                })
324                .collect::<Result<Vec<_>, _>>()?;
325            for (i, j) in (0..order).cartesian_product(0..order) {
326                let jinv = ctb
327                    .slice(s![.., j])
328                    .iter()
329                    .position(|&x| x == 0)
330                    .ok_or(format_err!(
331                        "Unable to find the inverse of group element `{j}`."
332                    ))?;
333                let jinv_i = ctb[(jinv, i)];
334                smat[(i, j)] = self.norm_preserving_scalar_map(jinv)?(ovs[jinv_i]);
335            }
336        } else {
337            log::debug!("Cayley table not available or the use of Cayley table not requested. Overlap matrix will be constructed without group-closure speed-up.");
338            for pair in self
339                .iter()
340                .map(|item_res| item_res.map_err(|err| err.to_string()))
341                .enumerate()
342                .combinations_with_replacement(2)
343            {
344                let (w, item_w_res) = &pair[0];
345                let (x, item_x_res) = &pair[1];
346                let item_w = item_w_res
347                    .as_ref()
348                    .map_err(|err| format_err!(err.clone()))
349                    .with_context(|| "One of the items in the orbit is not available")?;
350                let item_x = item_x_res
351                    .as_ref()
352                    .map_err(|err| format_err!(err.clone()))
353                    .with_context(|| "One of the items in the orbit is not available")?;
354                smat[(*w, *x)] = item_w.overlap(item_x, metric, metric_h).map_err(|err| {
355                    log::warn!("{err}");
356                    log::warn!(
357                        "Unable to calculate the overlap between items `{w}` and `{x}` in the orbit."
358                    );
359                    err
360                })?;
361                if *w != *x {
362                    smat[(*x, *w)] = item_x.overlap(item_w, metric, metric_h).map_err(|err| {
363                        log::warn!("{err}");
364                        log::warn!(
365                            "Unable to calculate the overlap between items `{x}` and `{w}` in the orbit."
366                        );
367                        err
368                    })?;
369                }
370            }
371        }
372        if self.origin().complex_symmetric() {
373            self.set_smat((smat.clone() + smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())))
374        } else {
375            self.set_smat(
376                (smat.clone() + smat.t().to_owned().mapv(|x| x.conj()))
377                    .mapv(|x| x / (T::one() + T::one())),
378            )
379        }
380        Ok(self)
381    }
382
383    /// Normalises overlap matrix between items in the orbit such that its diagonal entries are
384    /// unity.
385    ///
386    /// # Errors
387    ///
388    /// Errors if no orbit overlap matrix can be found, of if linear-algebraic errors are
389    /// encountered.
390    fn normalise_smat(&mut self) -> Result<&mut Self, anyhow::Error> {
391        let smat = self
392            .smat()
393            .ok_or(format_err!("No orbit overlap matrix to normalise."))?;
394        let norm = smat.diag().mapv(|x| <T as ComplexFloat>::sqrt(x));
395        let nspatial = norm.len();
396        let norm_col = norm
397            .clone()
398            .into_shape_with_order([nspatial, 1])
399            .map_err(|err| format_err!(err))?;
400        let norm_row = norm
401            .into_shape_with_order([1, nspatial])
402            .map_err(|err| format_err!(err))?;
403        let norm_mat = norm_col.dot(&norm_row);
404        let normalised_smat = smat / norm_mat;
405        self.set_smat(normalised_smat);
406        Ok(self)
407    }
408
409    /// Computes the $`\mathbf{T}(g)`$ matrix for a particular element $`g`$ of the generating
410    /// group.
411    ///
412    /// The elements of this matrix are given by
413    ///
414    /// ```math
415    ///     T_{wx}(g)
416    ///         = \langle \hat{\iota} \hat{g}_w \mathbf{v}_0, \hat{g} \hat{g}_x \mathbf{v}_0 \rangle.
417    /// ```
418    ///
419    /// This means that $`\mathbf{T}(g)`$ is just the orbit overlap matrix $`\mathbf{S}`$ with its
420    /// columns permuted according to the way $`g`$ composites on the elements in the group from
421    /// the left.
422    ///
423    /// # Arguments
424    ///
425    /// * `op` - The element $`g`$ in the generating group.
426    ///
427    /// # Returns
428    ///
429    /// The matrix $`\mathbf{T}(g)`$.
430    #[must_use]
431    fn calc_tmat(&self, op: &G::GroupElement) -> Result<Array2<T>, anyhow::Error> {
432        let ctb = self
433            .group()
434            .cayley_table()
435            .expect("The Cayley table for the group cannot be found.");
436        let i = self.group().get_index_of(op).unwrap_or_else(|| {
437            panic!("Unable to retrieve the index of element `{op}` in the group.")
438        });
439        let ix = ctb.slice(s![i, ..]).iter().cloned().collect::<Vec<_>>();
440        let twx = self
441            .smat()
442            .ok_or(format_err!("No orbit overlap matrix found."))?
443            .select(Axis(1), &ix);
444        Ok(twx)
445    }
446
447    /// Computes the representation or corepresentation matrix $`\mathbf{D}(g)`$ for a particular
448    /// element $`g`$ in the generating group in the basis of the orbit.
449    ///
450    /// The matrix $`\mathbf{D}(g)`$ is defined by
451    ///
452    /// ```math
453    ///     \hat{g} \mathcal{G} \cdot \mathbf{v}_0 = \mathcal{G} \cdot \mathbf{v}_0 \mathbf{D}(g),
454    /// ```
455    ///
456    /// where $`\mathcal{G} \cdot \mathbf{v}_0`$ is the orbit generated by the action of the group
457    /// $`\mathcal{G}`$ on the origin $`\mathbf{v}_0`$.
458    ///
459    /// # Arguments
460    ///
461    /// * `op` - The element $`g`$ of the generating group.
462    ///
463    /// # Returns
464    ///
465    /// The matrix $`\mathbf{D}(g)`$.
466    #[must_use]
467    fn calc_dmat(&self, op: &G::GroupElement) -> Result<Array2<T>, anyhow::Error> {
468        let complex_symmetric = self.origin().complex_symmetric();
469        let xmath = if complex_symmetric {
470            self.xmat().t().to_owned()
471        } else {
472            self.xmat().t().mapv(|x| x.conj())
473        };
474        let smattilde = xmath
475            .dot(
476                self.smat()
477                    .ok_or(format_err!("No orbit overlap matrix found."))?,
478            )
479            .dot(self.xmat());
480        let smattilde_inv = smattilde
481            .inv()
482            .expect("The inverse of S~ could not be found.");
483        let dmat = einsum(
484            "ij,jk,kl,lm->im",
485            &[&smattilde_inv, &xmath, &self.calc_tmat(op)?, self.xmat()],
486        )
487        .map_err(|err| format_err!(err))
488        .with_context(|| "Unable to compute the matrix product [(S~)^(-1) X† T X].")?
489        .into_dimensionality::<Ix2>()
490        .map_err(|err| format_err!(err))
491        .with_context(|| {
492            "Unable to convert the matrix product [(S~)^(-1) X† T X] to two dimensions."
493        });
494        dmat
495    }
496
497    /// Computes the character of a particular element $`g`$ in the generating group in the basis
498    /// of the orbit.
499    ///
500    /// See [`Self::calc_dmat`] for more information.
501    ///
502    /// # Arguments
503    ///
504    /// * `op` - The element $`g`$ of the generating group.
505    ///
506    /// # Returns
507    ///
508    /// The character $`\chi(g)`$.
509    #[must_use]
510    fn calc_character(&self, op: &G::GroupElement) -> Result<T, anyhow::Error> {
511        let complex_symmetric = self.origin().complex_symmetric();
512        let xmath = if complex_symmetric {
513            self.xmat().t().to_owned()
514        } else {
515            self.xmat().t().mapv(|x| x.conj())
516        };
517        let smattilde = xmath
518            .dot(
519                self.smat()
520                    .ok_or(format_err!("No orbit overlap matrix found."))?,
521            )
522            .dot(self.xmat());
523        let smattilde_inv = smattilde
524            .inv()
525            .expect("The inverse of S~ could not be found.");
526        let chi = einsum(
527            "ij,jk,kl,li",
528            &[&smattilde_inv, &xmath, &self.calc_tmat(op)?, self.xmat()],
529        )
530        .map_err(|err| format_err!(err))
531        .with_context(|| "Unable to compute the trace of the matrix product [(S~)^(-1) X† T X].")?
532        .into_dimensionality::<Ix0>()
533        .map_err(|err| format_err!(err))
534        .with_context(|| "Unable to convert the trace of the matrix product [(S~)^(-1) X† T X] to zero dimensions.")?;
535        chi.into_iter().next().ok_or(format_err!(
536            "Unable to extract the character from the representation matrix."
537        ))
538    }
539
540    /// Computes the characters of the elements in a conjugacy-class transversal of the generating
541    /// group in the basis of the orbit.
542    ///
543    /// See [`Self::calc_dmat`]  and [`Self::calc_character`] for more information.
544    ///
545    /// # Returns
546    ///
547    /// The conjugacy class symbols and the corresponding characters.
548    #[must_use]
549    fn calc_characters(
550        &self,
551    ) -> Result<Vec<(<G as ClassProperties>::ClassSymbol, T)>, anyhow::Error> {
552        let complex_symmetric = self.origin().complex_symmetric();
553        let xmath = if complex_symmetric {
554            self.xmat().t().to_owned()
555        } else {
556            self.xmat().t().mapv(|x| x.conj())
557        };
558        let smattilde = xmath
559            .dot(
560                self.smat()
561                    .ok_or(format_err!("No orbit overlap matrix found."))?,
562            )
563            .dot(self.xmat());
564        let smattilde_inv = smattilde
565            .inv()
566            .expect("The inverse of S~ could not be found.");
567        let chis = (0..self.group().class_number()).map(|cc_i| {
568            let cc = self.group().get_cc_symbol_of_index(cc_i).unwrap();
569            let op = self.group().get_cc_transversal(cc_i).unwrap();
570            let chi = einsum(
571                "ij,jk,kl,li",
572                &[&smattilde_inv, &xmath, &self.calc_tmat(&op)?, self.xmat()],
573            )
574            .map_err(|err| format_err!(err))
575            .with_context(|| "Unable to compute the trace of the matrix product [(S~)^(-1) X† T X].")?
576            .into_dimensionality::<Ix0>()
577            .map_err(|err| format_err!(err))
578            .with_context(|| "Unable to convert the trace of the matrix product [(S~)^(-1) X† T X] to zero dimensions.")?;
579            let chi_val = chi.into_iter().next().ok_or(format_err!(
580                "Unable to extract the character from the representation matrix."
581            ))?;
582            Ok((cc, chi_val))
583        }).collect::<Result<Vec<_>, _>>();
584        chis
585    }
586
587    /// Reduces the representation or corepresentation spanned by the items in the orbit to a
588    /// direct sum of the irreducible representations or corepresentations of the generating group.
589    ///
590    /// # Returns
591    ///
592    /// The decomposed result.
593    ///
594    /// # Errors
595    ///
596    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
597    /// are non-integral.
598    fn analyse_rep(
599        &self,
600    ) -> Result<
601        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
602        DecompositionError,
603    > {
604        let chis = self
605            .calc_characters()
606            .map_err(|err| DecompositionError(err.to_string()))?;
607        let res = self.group().character_table().reduce_characters(
608            &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
609            self.integrality_threshold(),
610        );
611        res
612    }
613
614    /// Converts a slice of tuples of class symbols and characters to a nicely formatted table.
615    ///
616    /// # Arguments
617    ///
618    /// * `chis` - A slice of tuples of class symbols and characters.
619    /// * `integrality_threshold` - Threshold for ascertaining the integrality of characters.
620    ///
621    /// # Returns
622    ///
623    /// A string of a nicely formatted table.
624    fn characters_to_string(
625        &self,
626        chis: &[(<G as ClassProperties>::ClassSymbol, T)],
627        integrality_threshold: <T as ComplexFloat>::Real,
628    ) -> String
629    where
630        T: ComplexFloat + Lapack + fmt::Debug,
631        <T as ComplexFloat>::Real: ToPrimitive,
632    {
633        let ndigits = (-integrality_threshold.log10()).to_usize().unwrap_or(6) + 1;
634        let (ccs, characters): (Vec<_>, Vec<_>) = chis
635            .iter()
636            .map(|(cc, chi)| (cc.to_string(), format!("{chi:+.ndigits$}")))
637            .unzip();
638        let cc_width = ccs
639            .iter()
640            .map(|cc| cc.chars().count())
641            .max()
642            .unwrap_or(5)
643            .max(5);
644        let characters_width = characters
645            .iter()
646            .map(|chi| chi.chars().count())
647            .max()
648            .unwrap_or(9)
649            .max(9);
650
651        let divider = "┈".repeat(cc_width + characters_width + 4);
652        let header = format!(" {:<cc_width$}  {:<}", "Class", "Character");
653        let body = Itertools::intersperse(
654            chis.iter()
655                .map(|(cc, chi)| format!(" {:<cc_width$}  {:<+.ndigits$}", cc.to_string(), chi)),
656            "\n".to_string(),
657        )
658        .collect::<String>();
659
660        Itertools::intersperse(
661            [divider.clone(), header, divider.clone(), body, divider].into_iter(),
662            "\n".to_string(),
663        )
664        .collect::<String>()
665    }
666}
667
668// ~~~~~~~~~~~~~~~~~~~~~~~
669// ProjectionDecomposition
670// ~~~~~~~~~~~~~~~~~~~~~~~
671
672/// Trait for decomposing the origin of an orbit into the irreducible components of the generating
673/// group using the projection operator.
674pub trait ProjectionDecomposition<G, I, T, D>: RepAnalysis<G, I, T, D>
675where
676    T: ComplexFloat + Lapack + fmt::Debug,
677    <T as ComplexFloat>::Real: ToPrimitive,
678    G: GroupProperties + ClassProperties + CharacterProperties + Sync + Send,
679    <G as CharacterProperties>::RowSymbol: Sync + Send,
680    G::GroupElement: fmt::Display,
681    G::CharTab: SubspaceDecomposable<T>,
682    D: Dimension,
683    I: Overlap<T, D> + Clone,
684    Self: Sync + Send,
685    Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
686    for<'a> Complex<f64>: Mul<&'a T, Output = Complex<f64>>,
687{
688    /// Calculates the irreducible compositions of the origin $`\mathbf{v}`$ of the orbit using the
689    /// projection operator:
690    ///
691    /// ```math
692    /// p_{\Gamma} =
693    ///     \frac{\dim \Gamma}{\lvert \mathcal{G} \rvert}
694    ///     \sum_{i=1}^{\lvert \mathcal{G} \rvert}
695    ///         \chi_{\Gamma}^{*}(\hat{g}_i)
696    ///         \braket{ \mathbf{v} | \hat{g}_i \mathbf{v} }
697    /// ```
698    ///
699    /// At the moment, this only works for unitary irreducible representations because of the
700    /// presence of the characters $`\chi_{\Gamma}(\hat{g}_i)`$ which can only be tabulated for
701    /// unitary operations.
702    ///
703    /// # Returns
704    ///
705    /// A vector of tuples, each of which contains the irreducible row label and the corresponding
706    /// projection value.
707    fn calc_projection_compositions(
708        &self,
709    ) -> Result<
710        Vec<(
711            <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
712            Complex<f64>,
713        )>,
714        anyhow::Error,
715    > {
716        self.group()
717            .character_table()
718            .get_all_rows()
719            .iter()
720            .map(|row_label| {
721                let group = self.group();
722                let chartab = group.character_table();
723                let id_class = group
724                    .get_cc_symbol_of_index(0)
725                    .ok_or(format_err!("Unable to retrieve the identity class."))?;
726                let dim = chartab
727                    .get_character(&row_label, &id_class).complex_value();
728                let group_order = group.order().to_f64().ok_or(
729                    DecompositionError("The group order cannot be converted to `f64`.".to_string())
730                )?;
731                let projection: Complex<f64> = (0..group.order())
732                    .into_par_iter()
733                    .try_fold(|| Complex::<f64>::zero(), |acc, i| {
734                        let cc_i = group
735                            .get_cc_of_element_index(i)
736                            .ok_or(format_err!("Unable to retrieve the conjugacy class index of element index {i}."))?;
737                        let cc = group
738                            .get_cc_symbol_of_index(cc_i)
739                            .ok_or(format_err!("Unable to retrieve the conjugacy class symbol of conjugacy class index {cc_i}."))?;
740                        let chi_i_star = group.character_table().get_character(&row_label, &cc).complex_conjugate();
741                        let s_0i = self
742                            .smat()
743                            .ok_or(format_err!("No orbit overlap matrix found."))?
744                            .get((0, i))
745                            .ok_or(format_err!("Unable to retrieve the overlap matrix element with index `(0, {i})`."))?;
746                        Ok::<_, anyhow::Error>(acc + chi_i_star.complex_value() * s_0i)
747                    })
748                    .try_reduce(|| Complex::<f64>::zero(), |a, s| Ok(a + s))?;
749                Ok::<_, anyhow::Error>((row_label.clone(), &dim * projection / group_order))
750            }).collect::<Result<Vec<_>, anyhow::Error>>()
751    }
752
753    /// Converts a slice of tuples of irreducible row symbols and projection values to a nicely
754    /// formatted table.
755    ///
756    /// # Arguments
757    ///
758    /// * `projections` - A slice of tuples of irreducible row symbols and projection values.
759    /// * `integrality_threshold` - Threshold for ascertaining the integrality of projection values.
760    ///
761    /// # Returns
762    ///
763    /// A string of a nicely formatted table.
764    fn projections_to_string(
765        &self,
766        projections: &[(
767            <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
768            Complex<f64>,
769        )],
770        integrality_threshold: <T as ComplexFloat>::Real,
771    ) -> String {
772        let ndigits = (-integrality_threshold.log10()).to_usize().unwrap_or(6) + 1;
773        let (row_labels, projection_values): (Vec<_>, Vec<_>) = projections
774            .iter()
775            .map(|(row, proj)| (row.to_string(), format!("{proj:+.ndigits$}")))
776            .unzip();
777        let row_width = row_labels
778            .iter()
779            .map(|row| row.chars().count())
780            .max()
781            .unwrap_or(9)
782            .max(9);
783        let proj_width = projection_values
784            .iter()
785            .map(|proj| proj.chars().count())
786            .max()
787            .unwrap_or(10)
788            .max(10);
789
790        let divider = "┈".repeat(row_width + proj_width + 4);
791        let header = format!(" {:<row_width$}  {:<}", "Component", "Projection");
792        let body = Itertools::intersperse(
793            projections.iter().map(|(row, proj)| {
794                format!(" {:<row_width$}  {:<+.ndigits$}", row.to_string(), proj)
795            }),
796            "\n".to_string(),
797        )
798        .collect::<String>();
799
800        Itertools::intersperse(
801            [divider.clone(), header, divider.clone(), body, divider].into_iter(),
802            "\n".to_string(),
803        )
804        .collect::<String>()
805    }
806}
807
808// Partial blanket implementation for orbits generated by unitary-represented symmetry groups over
809// real or complex numbers.
810#[duplicate_item(
811    duplicate!{
812        [ dtype_nested; [f64]; [Complex<f64>] ]
813        [
814            gtype_ [ UnitaryRepresentedSymmetryGroup ]
815            dtype_ [ dtype_nested ]
816        ]
817    }
818)]
819impl<O, I, D> ProjectionDecomposition<gtype_, I, dtype_, D> for O
820where
821    O: RepAnalysis<gtype_, I, dtype_, D>,
822    D: Dimension,
823    I: Overlap<dtype_, D> + Clone,
824    Self: Sync + Send,
825    Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
826{
827}
828
829// =================
830// Macro definitions
831// =================
832
833macro_rules! fn_calc_xmat_real {
834    ( $(#[$meta:meta])* $vis:vis $func:ident ) => {
835        $(#[$meta])*
836        $vis fn $func(&mut self, preserves_full_rank: bool) -> Result<&mut Self, anyhow::Error> {
837            // Real, symmetric S
838            let thresh = self.linear_independence_threshold;
839            let smat = self
840                .smat
841                .as_ref()
842                .ok_or(format_err!("No overlap matrix found for this orbit."))?;
843            use ndarray_linalg::norm::Norm;
844            if (smat.to_owned() - smat.t()).norm_l2() > thresh {
845                Err(format_err!("Overlap matrix is not symmetric."))
846            } else {
847                let (s_eig, umat) = smat.eigh(UPLO::Lower).map_err(|err| format_err!(err))?;
848                let nonzero_s_indices = match self.eigenvalue_comparison_mode {
849                    EigenvalueComparisonMode::Modulus => {
850                        s_eig.iter().positions(|x| x.abs() > thresh).collect_vec()
851                    }
852                    EigenvalueComparisonMode::Real => {
853                        s_eig.iter().positions(|x| *x > thresh).collect_vec()
854                    }
855                };
856                let nonzero_s_eig = s_eig.select(Axis(0), &nonzero_s_indices);
857                let nonzero_umat = umat.select(Axis(1), &nonzero_s_indices);
858                let nullity = smat.shape()[0] - nonzero_s_indices.len();
859                let xmat = if nullity == 0 && preserves_full_rank {
860                    Array2::eye(smat.shape()[0])
861                } else {
862                    let s_s = Array2::<f64>::from_diag(&nonzero_s_eig.mapv(|x| 1.0 / x.sqrt()));
863                    nonzero_umat.dot(&s_s)
864                };
865                self.smat_eigvals = Some(s_eig);
866                self.xmat = Some(xmat);
867                Ok(self)
868            }
869        }
870    }
871}
872
873macro_rules! fn_calc_xmat_complex {
874    ( $(#[$meta:meta])* $vis:vis $func:ident ) => {
875        $(#[$meta])*
876        $vis fn $func(&mut self, preserves_full_rank: bool) -> Result<&mut Self, anyhow::Error> {
877            // Complex S, symmetric or Hermitian
878            let thresh = self.linear_independence_threshold;
879            let smat = self
880                .smat
881                .as_ref()
882                .ok_or(format_err!("No overlap matrix found for this orbit."))?;
883            let (s_eig, umat_nonortho) = smat.eig().map_err(|err| format_err!(err))?;
884
885            let nonzero_s_indices = match self.eigenvalue_comparison_mode {
886                EigenvalueComparisonMode::Modulus => s_eig
887                    .iter()
888                    .positions(|x| ComplexFloat::abs(*x) > thresh)
889                    .collect_vec(),
890                EigenvalueComparisonMode::Real => {
891                    if s_eig
892                        .iter()
893                        .any(|x| Float::abs(ComplexFloat::im(*x)) > thresh)
894                    {
895                        log::warn!("Comparing eigenvalues using the real parts, but not all eigenvalues are real.");
896                    }
897                    s_eig
898                        .iter()
899                        .positions(|x| ComplexFloat::re(*x) > thresh)
900                        .collect_vec()
901                }
902            };
903            let nonzero_s_eig = s_eig.select(Axis(0), &nonzero_s_indices);
904            let nonzero_umat_nonortho = umat_nonortho.select(Axis(1), &nonzero_s_indices);
905
906            // `eig` does not guarantee orthogonality of `nonzero_umat_nonortho`.
907            // Gram--Schmidt is therefore required.
908            let nonzero_umat = complex_modified_gram_schmidt(
909                &nonzero_umat_nonortho,
910                self.origin.complex_symmetric(),
911                thresh,
912            )
913            .map_err(
914                |_| format_err!("Unable to orthonormalise the linearly-independent eigenvectors of the overlap matrix.")
915            )?;
916
917            let nullity = smat.shape()[0] - nonzero_s_indices.len();
918            let xmat = if nullity == 0 && preserves_full_rank {
919                Array2::<Complex<T>>::eye(smat.shape()[0])
920            } else {
921                let s_s = Array2::<Complex<T>>::from_diag(
922                    &nonzero_s_eig.mapv(|x| Complex::<T>::from(T::one()) / x.sqrt()),
923                );
924                nonzero_umat.dot(&s_s)
925            };
926            self.smat_eigvals = Some(s_eig);
927            self.xmat = Some(xmat);
928            Ok(self)
929        }
930    }
931}
932
933pub(crate) use fn_calc_xmat_complex;
934pub(crate) use fn_calc_xmat_real;
935
936// =================
937// Utility functions
938// =================
939
940/// Logs overlap eigenvalues nicely and indicates where the threshold has been crossed.
941///
942/// # Arguments
943///
944/// * `eigvals` - The eigenvalues.
945/// * `thresh` - The cut-off threshold to be marked out.
946pub(crate) fn log_overlap_eigenvalues<T>(
947    title: &str,
948    eigvals: &Array1<T>,
949    thresh: <T as ComplexFloat>::Real,
950    eigenvalue_comparison_mode: &EigenvalueComparisonMode,
951) where
952    T: std::fmt::LowerExp + ComplexFloat,
953    <T as ComplexFloat>::Real: std::fmt::LowerExp,
954{
955    let mut eigvals_sorted = eigvals.iter().collect::<Vec<_>>();
956    match eigenvalue_comparison_mode {
957        EigenvalueComparisonMode::Modulus => {
958            eigvals_sorted.sort_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap());
959        }
960        EigenvalueComparisonMode::Real => {
961            eigvals_sorted.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
962        }
963    }
964    eigvals_sorted.reverse();
965    let eigvals_str = eigvals_sorted
966        .iter()
967        .map(|v| format!("{v:+.3e}"))
968        .collect::<Vec<_>>();
969    log_subtitle(title);
970    qsym2_output!("");
971
972    match eigenvalue_comparison_mode {
973        EigenvalueComparisonMode::Modulus => {
974            qsym2_output!("Eigenvalues are sorted in decreasing order of their moduli.");
975        }
976        EigenvalueComparisonMode::Real => {
977            qsym2_output!("Eigenvalues are sorted in decreasing order of their real parts.");
978        }
979    }
980    let count_length = usize::try_from(eigvals.len().ilog10() + 2).unwrap_or(2);
981    let eigval_length = eigvals_str
982        .iter()
983        .map(|v| v.chars().count())
984        .max()
985        .unwrap_or(20);
986    qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
987    qsym2_output!("{:>count_length$}  Eigenvalue", "#");
988    qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
989    let mut write_thresh = false;
990    for (i, eigval) in eigvals_str.iter().enumerate() {
991        let cmp = match eigenvalue_comparison_mode {
992            EigenvalueComparisonMode::Modulus => {
993                eigvals_sorted[i].abs().partial_cmp(&thresh).expect(
994                    "Unable to compare the modulus of an eigenvalue with the specified threshold.",
995                )
996            }
997            EigenvalueComparisonMode::Real => eigvals_sorted[i].re().partial_cmp(&thresh).expect(
998                "Unable to compare the real part of an eigenvalue with the specified threshold.",
999            ),
1000        };
1001        if cmp == Ordering::Less && !write_thresh {
1002            let comparison_mode_str = match eigenvalue_comparison_mode {
1003                EigenvalueComparisonMode::Modulus => "modulus-based",
1004                EigenvalueComparisonMode::Real => "real-part-based",
1005            };
1006            qsym2_output!(
1007                "{} <-- linear independence threshold ({comparison_mode_str}): {:+.3e}",
1008                "-".repeat(count_length + 3 + eigval_length),
1009                thresh
1010            );
1011            write_thresh = true;
1012        }
1013        qsym2_output!("{i:>count_length$}  {eigval}",);
1014    }
1015    qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
1016}