qsym2/drivers/representation_analysis/
slater_determinant.rs

1//! Driver for symmetry analysis of Slater determinants.
2
3use std::collections::HashMap;
4use std::fmt;
5use std::ops::Mul;
6
7use anyhow::{self, bail, format_err};
8use approx::abs_diff_eq;
9use derive_builder::Builder;
10use duplicate::duplicate_item;
11use indexmap::IndexMap;
12use itertools::Itertools;
13use ndarray::{s, Array2, Array4};
14use ndarray_linalg::types::Lapack;
15use num::ToPrimitive;
16use num_complex::{Complex, ComplexFloat};
17use num_traits::real::Real;
18use num_traits::Float;
19use rayon::prelude::*;
20use serde::{Deserialize, Serialize};
21
22use crate::analysis::{
23    log_overlap_eigenvalues, EigenvalueComparisonMode, Orbit, Overlap, ProjectionDecomposition,
24    RepAnalysis,
25};
26use crate::angmom::spinor_rotation_3d::{SpinConstraint, SpinOrbitCoupled, StructureConstraint};
27use crate::chartab::chartab_group::CharacterProperties;
28use crate::chartab::chartab_symbols::ReducibleLinearSpaceSymbol;
29use crate::chartab::{CharacterTable, SubspaceDecomposable};
30use crate::drivers::representation_analysis::angular_function::{
31    find_angular_function_representation, find_spinor_function_representation,
32    AngularFunctionRepAnalysisParams,
33};
34use crate::drivers::representation_analysis::{
35    fn_construct_magnetic_group, fn_construct_unitary_group, log_bao, log_cc_transversal,
36    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
37};
38use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
39use crate::drivers::QSym2Driver;
40use crate::group::{GroupProperties, MagneticRepresentedGroup, UnitaryRepresentedGroup};
41use crate::io::format::{
42    log_subtitle, nice_bool, qsym2_output, write_subtitle, write_title, QSym2Output,
43};
44use crate::symmetry::symmetry_element::symmetry_operation::{
45    SpecialSymmetryTransformation, SymmetryOperation,
46};
47use crate::symmetry::symmetry_group::{
48    MagneticRepresentedSymmetryGroup, SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup,
49};
50use crate::symmetry::symmetry_symbols::{
51    deduce_mirror_parities, MirrorParity, MullikenIrcorepSymbol, SymmetryClassSymbol,
52};
53use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
54use crate::target::density::density_analysis::DensitySymmetryOrbit;
55use crate::target::determinant::determinant_analysis::SlaterDeterminantSymmetryOrbit;
56use crate::target::determinant::SlaterDeterminant;
57use crate::target::orbital::orbital_analysis::generate_det_mo_orbits;
58
59#[cfg(test)]
60#[path = "slater_determinant_tests.rs"]
61mod slater_determinant_tests;
62
63// ==================
64// Struct definitions
65// ==================
66
67// ----------
68// Parameters
69// ----------
70
71const fn default_true() -> bool {
72    true
73}
74const fn default_symbolic() -> Option<CharacterTableDisplay> {
75    Some(CharacterTableDisplay::Symbolic)
76}
77
78/// Structure containing control parameters for Slater determinant representation analysis.
79#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
80pub struct SlaterDeterminantRepAnalysisParams<T: From<f64>> {
81    /// Threshold for checking if subspace multiplicities are integral.
82    pub integrality_threshold: T,
83
84    /// Threshold for determining zero eigenvalues in the orbit overlap matrix.
85    pub linear_independence_threshold: T,
86
87    /// Boolean indicating if molecular orbital symmetries are to be analysed alongside the overall
88    /// determinantal symmetry.
89    #[builder(default = "true")]
90    #[serde(default = "default_true")]
91    pub analyse_mo_symmetries: bool,
92
93    /// Boolean indicating if molecular orbital mirror parities are to be analysed alongside
94    /// molecular orbital symmetries.
95    #[builder(default = "false")]
96    #[serde(default)]
97    pub analyse_mo_mirror_parities: bool,
98
99    /// Boolean indicating if symmetry projection should be performed for molecular orbitals.
100    #[builder(default = "false")]
101    #[serde(default)]
102    pub analyse_mo_symmetry_projections: bool,
103
104    /// Boolean indicating if density symmetries are to be analysed alongside wavefunction symmetries
105    /// for this determinant.
106    #[builder(default = "false")]
107    #[serde(default)]
108    pub analyse_density_symmetries: bool,
109
110    /// Option indicating if the magnetic group is to be used for symmetry analysis, and if so,
111    /// whether unitary representations or unitary-antiunitary corepresentations should be used.
112    #[builder(default = "None")]
113    #[serde(default)]
114    pub use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
115
116    /// Boolean indicating if the double group is to be used for symmetry analysis.
117    #[builder(default = "false")]
118    #[serde(default)]
119    pub use_double_group: bool,
120
121    /// Boolean indicating if the Cayley table of the group, if available, should be used to speed
122    /// up the computation of orbit overlap matrices.
123    #[builder(default = "true")]
124    #[serde(default = "default_true")]
125    pub use_cayley_table: bool,
126
127    /// The kind of symmetry transformation to be applied on the reference determinant to generate
128    /// the orbit for symmetry analysis.
129    #[builder(default = "SymmetryTransformationKind::Spatial")]
130    #[serde(default)]
131    pub symmetry_transformation_kind: SymmetryTransformationKind,
132
133    /// Option indicating if the character table of the group used for symmetry analysis is to be
134    /// printed out.
135    #[builder(default = "Some(CharacterTableDisplay::Symbolic)")]
136    #[serde(default = "default_symbolic")]
137    pub write_character_table: Option<CharacterTableDisplay>,
138
139    /// Boolean indicating if the eigenvalues of the orbit overlap matrix are to be printed out.
140    #[builder(default = "true")]
141    #[serde(default = "default_true")]
142    pub write_overlap_eigenvalues: bool,
143
144    /// The comparison mode for filtering out orbit overlap eigenvalues.
145    #[builder(default = "EigenvalueComparisonMode::Modulus")]
146    #[serde(default)]
147    pub eigenvalue_comparison_mode: EigenvalueComparisonMode,
148
149    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
150    /// subgroup of an infinite group can be used for the symmetry analysis.
151    #[builder(default = "None")]
152    #[serde(default)]
153    pub infinite_order_to_finite: Option<u32>,
154}
155
156impl<T> SlaterDeterminantRepAnalysisParams<T>
157where
158    T: Float + From<f64>,
159{
160    /// Returns a builder to construct a [`SlaterDeterminantRepAnalysisParams`] structure.
161    pub fn builder() -> SlaterDeterminantRepAnalysisParamsBuilder<T> {
162        SlaterDeterminantRepAnalysisParamsBuilder::default()
163    }
164}
165
166impl Default for SlaterDeterminantRepAnalysisParams<f64> {
167    fn default() -> Self {
168        Self::builder()
169            .integrality_threshold(1e-7)
170            .linear_independence_threshold(1e-7)
171            .build()
172            .expect("Unable to construct a default `SlaterDeterminantRepAnalysisParams<f64>`.")
173    }
174}
175
176impl<T> fmt::Display for SlaterDeterminantRepAnalysisParams<T>
177where
178    T: From<f64> + fmt::LowerExp + fmt::Debug,
179{
180    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
181        writeln!(
182            f,
183            "Integrality threshold: {:.3e}",
184            self.integrality_threshold
185        )?;
186        writeln!(
187            f,
188            "Linear independence threshold: {:.3e}",
189            self.linear_independence_threshold
190        )?;
191        writeln!(
192            f,
193            "Orbit eigenvalue comparison mode: {}",
194            self.eigenvalue_comparison_mode
195        )?;
196        writeln!(
197            f,
198            "Write overlap eigenvalues: {}",
199            nice_bool(self.write_overlap_eigenvalues)
200        )?;
201        writeln!(f)?;
202        writeln!(
203            f,
204            "Analyse molecular orbital symmetry: {}",
205            nice_bool(self.analyse_mo_symmetries)
206        )?;
207        writeln!(
208            f,
209            "Analyse molecular orbital symmetry projection: {}",
210            nice_bool(self.analyse_mo_symmetry_projections)
211        )?;
212        writeln!(
213            f,
214            "Analyse molecular orbital mirror parity: {}",
215            nice_bool(self.analyse_mo_mirror_parities)
216        )?;
217        writeln!(
218            f,
219            "Analyse density symmetry: {}",
220            nice_bool(self.analyse_density_symmetries)
221        )?;
222        writeln!(f)?;
223        writeln!(
224            f,
225            "Use magnetic group for analysis: {}",
226            match self.use_magnetic_group {
227                None => "no",
228                Some(MagneticSymmetryAnalysisKind::Representation) =>
229                    "yes, using unitary representations",
230                Some(MagneticSymmetryAnalysisKind::Corepresentation) =>
231                    "yes, using magnetic corepresentations",
232            }
233        )?;
234        writeln!(
235            f,
236            "Use double group for analysis: {}",
237            nice_bool(self.use_double_group)
238        )?;
239        writeln!(
240            f,
241            "Use Cayley table for orbit overlap matrices: {}",
242            nice_bool(self.use_cayley_table)
243        )?;
244        if let Some(finite_order) = self.infinite_order_to_finite {
245            writeln!(f, "Infinite order to finite: {finite_order}")?;
246        }
247        writeln!(
248            f,
249            "Symmetry transformation kind: {}",
250            self.symmetry_transformation_kind
251        )?;
252        writeln!(f)?;
253        writeln!(
254            f,
255            "Write character table: {}",
256            if let Some(chartab_display) = self.write_character_table.as_ref() {
257                format!("yes, {}", chartab_display.to_string().to_lowercase())
258            } else {
259                "no".to_string()
260            }
261        )?;
262
263        Ok(())
264    }
265}
266
267// ------
268// Result
269// ------
270
271/// Structure to contain Slater determinant representation analysis results.
272#[derive(Clone, Builder)]
273pub struct SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
274where
275    G: SymmetryGroupProperties + Clone,
276    G::CharTab: SubspaceDecomposable<T>,
277    T: ComplexFloat + Lapack,
278    SC: StructureConstraint + fmt::Display,
279    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
280{
281    /// The control parameters used to obtain this set of Slater determinant representation
282    /// analysis results.
283    parameters: &'a SlaterDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
284
285    /// The Slater determinant being analysed.
286    determinant: &'a SlaterDeterminant<'a, T, SC>,
287
288    /// The group used for the representation analysis.
289    group: G,
290
291    /// The deduced overall symmetry of the determinant.
292    determinant_symmetry: Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
293
294    /// The deduced symmetries of the molecular orbitals constituting the determinant, if required.
295    mo_symmetries: Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>>,
296
297    /// The deduced symmetry projections of the molecular orbitals constituting the determinant, if required.
298    mo_symmetry_projections: Option<
299        Vec<
300            Vec<
301                Option<
302                    Vec<(
303                        <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
304                        Complex<f64>,
305                    )>,
306                >,
307            >,
308        >,
309    >,
310
311    /// The deduced mirror parities of the molecular orbitals constituting the determinant, if required.
312    mo_mirror_parities:
313        Option<Vec<Vec<Option<IndexMap<SymmetryClassSymbol<SymmetryOperation>, MirrorParity>>>>>,
314
315    /// The overlap eigenvalues above and below the linear independence threshold for each
316    /// molecular orbital symmetry deduction.
317    mo_symmetries_thresholds: Option<Vec<Vec<(Option<T>, Option<T>)>>>,
318
319    /// The deduced symmetries of the various densities constructible from the determinant, if
320    /// required. In each tuple, the first element gives a description of the density corresponding
321    /// to the symmetry result.
322    determinant_density_symmetries: Option<
323        Vec<(
324            String,
325            Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
326        )>,
327    >,
328
329    /// The deduced symmetries of the total densities of the molecular orbitals constituting the
330    /// determinant, if required.
331    mo_density_symmetries:
332        Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>>,
333}
334
335impl<'a, G, T, SC> SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
336where
337    G: SymmetryGroupProperties + Clone,
338    G::CharTab: SubspaceDecomposable<T>,
339    T: ComplexFloat + Lapack,
340    SC: StructureConstraint + Clone + fmt::Display,
341    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
342{
343    /// Returns a builder to construct a new [`SlaterDeterminantRepAnalysisResultBuilder`]
344    /// structure.
345    fn builder() -> SlaterDeterminantRepAnalysisResultBuilder<'a, G, T, SC> {
346        SlaterDeterminantRepAnalysisResultBuilder::default()
347    }
348}
349
350impl<'a, G, T, SC> SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
351where
352    G: SymmetryGroupProperties + Clone,
353    G::CharTab: SubspaceDecomposable<T>,
354    T: ComplexFloat + Lapack,
355    SC: StructureConstraint + fmt::Display,
356    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
357{
358    /// Returns the group used for the representation analysis.
359    pub fn group(&self) -> &G {
360        &self.group
361    }
362
363    /// Returns the determinant symmetry obtained from the analysis result.
364    pub fn determinant_symmetry(
365        &self,
366    ) -> &Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String> {
367        &self.determinant_symmetry
368    }
369
370    /// Returns the deduced symmetries of the molecular orbitals constituting the determinant, if required.
371    pub fn mo_symmetries(
372        &self,
373    ) -> &Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>> {
374        &self.mo_symmetries
375    }
376
377    /// Returns the deduced symmetry projections of the molecular orbitals constituting the determinant, if required.
378    pub fn mo_symmetry_projections(
379        &self,
380    ) -> &Option<
381        Vec<
382            Vec<
383                Option<
384                    Vec<(
385                        <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
386                        Complex<f64>,
387                    )>,
388                >,
389            >,
390        >,
391    > {
392        &self.mo_symmetry_projections
393    }
394
395    /// Returns the deduced symmetries of the various densities constructible from the determinant,
396    /// if required. In each tuple, the first element gives a description of the density corresponding
397    /// to the symmetry result.
398    pub fn determinant_density_symmetries(
399        &self,
400    ) -> &Option<
401        Vec<(
402            String,
403            Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
404        )>,
405    > {
406        &self.determinant_density_symmetries
407    }
408
409    /// Returns the deduced symmetries of the total densities of the molecular orbitals constituting
410    /// the determinant, if required.
411    pub fn mo_density_symmetries(
412        &self,
413    ) -> &Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>> {
414        &self.mo_density_symmetries
415    }
416}
417
418impl<'a, G, T, SC> fmt::Display for SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
419where
420    G: SymmetryGroupProperties + Clone,
421    G::CharTab: SubspaceDecomposable<T>,
422    T: ComplexFloat + Lapack,
423    SC: StructureConstraint + fmt::Display,
424    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
425{
426    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
427        write_subtitle(f, "Orbit-based symmetry analysis results")?;
428        writeln!(f)?;
429        writeln!(
430            f,
431            "> Group: {} ({})",
432            self.group
433                .finite_subgroup_name()
434                .map(|subgroup_name| format!("{} > {}", self.group.name(), subgroup_name))
435                .unwrap_or(self.group.name()),
436            self.group.group_type().to_string().to_lowercase()
437        )?;
438        writeln!(f)?;
439        writeln!(f, "> Overall determinantal result")?;
440        writeln!(
441            f,
442            "  Energy  : {}",
443            self.determinant
444                .energy()
445                .map(|e| e.to_string())
446                .unwrap_or_else(|err| format!("-- ({err})"))
447        )?;
448        writeln!(
449            f,
450            "  Symmetry: {}",
451            self.determinant_symmetry
452                .as_ref()
453                .map(|s| s.to_string())
454                .unwrap_or_else(|err| format!("-- ({err})"))
455        )?;
456        writeln!(f)?;
457
458        if let Some(den_syms) = self.determinant_density_symmetries.as_ref() {
459            writeln!(f, "> Overall determinantal density result")?;
460            let den_type_width = den_syms
461                .iter()
462                .map(|(den_type, _)| den_type.chars().count())
463                .max()
464                .unwrap_or(7)
465                .max(7);
466            for (den_type, den_sym_res) in den_syms.iter() {
467                writeln!(
468                    f,
469                    "  {den_type:<den_type_width$}: {}",
470                    den_sym_res
471                        .as_ref()
472                        .map(|e| e.to_string())
473                        .unwrap_or_else(|err| format!("-- ({err})"))
474                )?;
475            }
476            writeln!(f)?;
477        }
478
479        if let Some(mo_symmetries) = self.mo_symmetries.as_ref() {
480            let mo_spin_index_length = 4;
481            let mo_index_length = mo_symmetries
482                .iter()
483                .map(|spin_mo_symmetries| spin_mo_symmetries.len())
484                .max()
485                .and_then(|max_mo_length| usize::try_from(max_mo_length.ilog10() + 2).ok())
486                .unwrap_or(4);
487            let mo_occ_length = 5;
488            let mo_symmetry_length = mo_symmetries
489                .iter()
490                .flat_map(|spin_mo_symmetries| {
491                    spin_mo_symmetries.iter().map(|mo_sym| {
492                        mo_sym
493                            .as_ref()
494                            .map(|sym| sym.to_string())
495                            .unwrap_or("--".to_string())
496                            .chars()
497                            .count()
498                    })
499                })
500                .max()
501                .unwrap_or(0)
502                .max(8);
503            let mo_energies_opt = self.determinant.mo_energies();
504            let mo_energy_length = mo_energies_opt
505                .map(|mo_energies| {
506                    mo_energies
507                        .iter()
508                        .flat_map(|spin_mo_energies| {
509                            spin_mo_energies.map(|v| format!("{v:+.7}").chars().count())
510                        })
511                        .max()
512                        .unwrap_or(0)
513                })
514                .unwrap_or(0)
515                .max(6);
516
517            let mo_eig_above_length: usize = self
518                .mo_symmetries_thresholds
519                .as_ref()
520                .map(|mo_symmetries_thresholds| {
521                    mo_symmetries_thresholds
522                        .iter()
523                        .flat_map(|spin_mo_symmetries_thresholds| {
524                            spin_mo_symmetries_thresholds.iter().map(|(above, _)| {
525                                above
526                                    .as_ref()
527                                    .map(|eig| format!("{eig:+.3e}"))
528                                    .unwrap_or("--".to_string())
529                                    .chars()
530                                    .count()
531                            })
532                        })
533                        .max()
534                        .unwrap_or(10)
535                        .max(10)
536                })
537                .unwrap_or(10);
538            let mo_eig_below_length: usize = self
539                .mo_symmetries_thresholds
540                .as_ref()
541                .map(|mo_symmetries_thresholds| {
542                    mo_symmetries_thresholds
543                        .iter()
544                        .flat_map(|spin_mo_symmetries_thresholds| {
545                            spin_mo_symmetries_thresholds.iter().map(|(_, below)| {
546                                below
547                                    .as_ref()
548                                    .map(|eig| format!("{eig:+.3e}"))
549                                    .unwrap_or("--".to_string())
550                                    .chars()
551                                    .count()
552                            })
553                        })
554                        .max()
555                        .unwrap_or(10)
556                        .max(10)
557                })
558                .unwrap_or(10);
559
560            let precision = Real::ceil(ComplexFloat::abs(ComplexFloat::log10(
561                self.parameters.linear_independence_threshold,
562            )))
563            .to_usize()
564            .expect("Unable to convert the linear independence threshold exponent to `usize`.");
565            let mo_sym_projections_str_opt =
566                self.mo_symmetry_projections
567                    .as_ref()
568                    .map(|mo_sym_projectionss| {
569                        mo_sym_projectionss
570                            .iter()
571                            .enumerate()
572                            .map(|(ispin, mo_sym_projections)| {
573                                mo_sym_projections
574                                    .iter()
575                                    .enumerate()
576                                    .map(|(imo, mo_sym_projection)| {
577                                        mo_sym_projection
578                                            .as_ref()
579                                            .map(|sym_proj| {
580                                                if let Some(mo_symmetriess) = self.mo_symmetries() {
581                                                    if let Some(sym) = &mo_symmetriess[ispin][imo] {
582                                                        let sym_proj_hashmap = sym_proj
583                                                            .iter()
584                                                            .cloned()
585                                                            .collect::<HashMap<_, _>>();
586                                                        sym.subspaces()
587                                                            .iter()
588                                                            .map(|(subspace, _)| {
589                                                                format!(
590                                                                    "{subspace}: {}",
591                                                                    sym_proj_hashmap
592                                                                        .get(subspace)
593                                                                        .map(|composition| {
594                                                                            if abs_diff_eq!(
595                                                                                composition.im,
596                                                                                0.0,
597                                                                                epsilon = self.parameters.linear_independence_threshold.to_f64().expect("Unable to convert the linear independence threshold to `f64`.")
598                                                                            ) {
599                                                                                format!(
600                                                                                    "{:+.precision$}",
601                                                                                    composition.re
602                                                                                )
603                                                                            } else {
604                                                                                format!(
605                                                                            "{composition:+.precision$}")
606                                                                            }
607                                                                        })
608                                                                        .unwrap_or(
609                                                                            "--".to_string()
610                                                                        )
611                                                                )
612                                                            })
613                                                            .join(", ")
614                                                    } else {
615                                                        "--".to_string()
616                                                    }
617                                                } else {
618                                                    "--".to_string()
619                                                }
620                                            })
621                                            .unwrap_or("--".to_string())
622                                    })
623                                    .collect::<Vec<_>>()
624                            })
625                            .collect::<Vec<_>>()
626                    });
627            let mo_sym_projections_length_opt =
628                mo_sym_projections_str_opt
629                    .as_ref()
630                    .map(|mo_sym_projectionss| {
631                        mo_sym_projectionss
632                            .iter()
633                            .flat_map(|mo_sym_projections| {
634                                mo_sym_projections
635                                    .iter()
636                                    .map(|mo_sym_projection| mo_sym_projection.chars().count())
637                            })
638                            .max()
639                            .unwrap_or(18)
640                            .max(18)
641                    });
642            let mo_sym_projections_length = mo_sym_projections_length_opt.unwrap_or(0);
643            let mo_sym_projections_gap = mo_sym_projections_length_opt.map(|_| 2).unwrap_or(0);
644            let mo_sym_projections_heading = mo_sym_projections_length_opt
645                .map(|_| "MO sym. projection")
646                .unwrap_or("");
647
648            let mirrors = self
649                .group
650                .filter_cc_symbols(|cc| cc.is_spatial_reflection());
651            let mo_mirror_parities_length_opt = self.mo_mirror_parities.as_ref().map(|_| {
652                let mirror_heading = mirrors.iter().map(|sigma| format!("p[{sigma}]")).join("  ");
653                let length = mirror_heading.chars().count();
654                (mirror_heading, length)
655            });
656            let mo_mirror_parities_gap = mo_mirror_parities_length_opt
657                .as_ref()
658                .map(|_| 2)
659                .unwrap_or(0);
660            let (mo_mirror_parities_heading, mo_mirror_parities_length) =
661                mo_mirror_parities_length_opt.unwrap_or((String::new(), 0));
662
663            let mo_den_symss_str_opt = self.mo_density_symmetries.as_ref().map(|mo_den_symss| {
664                mo_den_symss
665                    .iter()
666                    .map(|mo_den_syms| {
667                        mo_den_syms
668                            .iter()
669                            .map(|mo_den_sym| {
670                                mo_den_sym
671                                    .as_ref()
672                                    .map(|sym| sym.to_string())
673                                    .unwrap_or("--".to_string())
674                            })
675                            .collect::<Vec<_>>()
676                    })
677                    .collect::<Vec<_>>()
678            });
679            let mo_density_length_opt = mo_den_symss_str_opt.as_ref().map(|mo_den_symss| {
680                mo_den_symss
681                    .iter()
682                    .flat_map(|mo_den_syms| {
683                        mo_den_syms
684                            .iter()
685                            .map(|mo_den_sym| mo_den_sym.chars().count())
686                    })
687                    .max()
688                    .unwrap_or(12)
689                    .max(12)
690            });
691            let mo_density_length = mo_density_length_opt.unwrap_or(0);
692            let mo_density_gap = mo_density_length_opt.map(|_| 2).unwrap_or(0);
693            let mo_density_heading = mo_density_length_opt.map(|_| "Density sym.").unwrap_or("");
694
695            let table_width = 14
696                + mo_spin_index_length
697                + mo_index_length
698                + mo_occ_length
699                + mo_energy_length
700                + mo_symmetry_length
701                + mo_mirror_parities_gap
702                + mo_mirror_parities_length
703                + mo_eig_above_length
704                + mo_eig_below_length
705                + mo_sym_projections_gap
706                + mo_sym_projections_length
707                + mo_density_gap
708                + mo_density_length;
709
710            writeln!(f, "> Molecular orbital results")?;
711            writeln!(
712                f,
713                "  Structure constraint: {}",
714                self.determinant
715                    .structure_constraint()
716                    .to_string()
717                    .to_lowercase()
718            )?;
719            if self.mo_mirror_parities.as_ref().is_some() {
720                writeln!(f, "")?;
721                writeln!(
722                    f,
723                    "Column p[σ] gives the parity under the reflection class σ: {} => even, {} => odd, {} => neither.",
724                    MirrorParity::Even,
725                    MirrorParity::Odd,
726                    MirrorParity::Neither
727                )?;
728            }
729            writeln!(f, "{}", "┈".repeat(table_width))?;
730            if mo_density_length > 0 {
731                writeln!(
732                    f,
733                    " {:>mo_spin_index_length$}  {:>mo_index_length$}  {:<mo_occ_length$}  {:<mo_energy_length$}  {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$}  {:<mo_eig_above_length$}  {:<mo_eig_below_length$}{}{:mo_sym_projections_length$}{}{}",
734                    "Spin",
735                    "MO",
736                    "Occ.",
737                    "Energy",
738                    "Symmetry",
739                    " ".repeat(mo_mirror_parities_gap),
740                    mo_mirror_parities_heading,
741                    "Eig. above",
742                    "Eig. below",
743                    " ".repeat(mo_sym_projections_gap),
744                    mo_sym_projections_heading,
745                    " ".repeat(mo_density_gap),
746                    mo_density_heading
747                )?;
748            } else if mo_sym_projections_length > 0 {
749                writeln!(
750                    f,
751                    " {:>mo_spin_index_length$}  {:>mo_index_length$}  {:<mo_occ_length$}  {:<mo_energy_length$}  {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$}  {:<mo_eig_above_length$}  {:<mo_eig_below_length$}{}{}",
752                    "Spin",
753                    "MO",
754                    "Occ.",
755                    "Energy",
756                    "Symmetry",
757                    " ".repeat(mo_mirror_parities_gap),
758                    mo_mirror_parities_heading,
759                    "Eig. above",
760                    "Eig. below",
761                    " ".repeat(mo_sym_projections_gap),
762                    mo_sym_projections_heading,
763                )?;
764            } else {
765                writeln!(
766                    f,
767                    " {:>mo_spin_index_length$}  {:>mo_index_length$}  {:<mo_occ_length$}  {:<mo_energy_length$}  {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$}  {:<mo_eig_above_length$}  {}",
768                    "Spin",
769                    "MO",
770                    "Occ.",
771                    "Energy",
772                    "Symmetry",
773                    " ".repeat(mo_mirror_parities_gap),
774                    mo_mirror_parities_heading,
775                    "Eig. above",
776                    "Eig. below",
777                )?;
778            };
779            writeln!(f, "{}", "┈".repeat(table_width))?;
780
781            let empty_string = String::new();
782            for (spini, spin_mo_symmetries) in mo_symmetries.iter().enumerate() {
783                writeln!(f, " Spin {spini}")?;
784                for (moi, mo_sym) in spin_mo_symmetries.iter().enumerate() {
785                    let occ_str = self
786                        .determinant
787                        .occupations()
788                        .get(spini)
789                        .and_then(|spin_occs| spin_occs.get(moi))
790                        .map(|occ| format!("{occ:>.3}"))
791                        .unwrap_or("--".to_string());
792                    let mo_energy_str = mo_energies_opt
793                        .and_then(|mo_energies| mo_energies.get(spini))
794                        .and_then(|spin_mo_energies| spin_mo_energies.get(moi))
795                        .map(|mo_energy| format!("{mo_energy:>+mo_energy_length$.7}"))
796                        .unwrap_or("--".to_string());
797                    let mo_sym_str = mo_sym
798                        .as_ref()
799                        .map(|sym| sym.to_string())
800                        .unwrap_or("--".to_string());
801
802                    let mo_mirror_parities_str = self
803                        .mo_mirror_parities
804                        .as_ref()
805                        .and_then(|mo_mirror_paritiess| {
806                            mo_mirror_paritiess
807                                .get(spini)
808                                .and_then(|spin_mo_mirror_parities| {
809                                    spin_mo_mirror_parities
810                                        .get(moi)
811                                        .map(|mo_mirror_parities_opt| {
812                                            mo_mirror_parities_opt
813                                                .as_ref()
814                                                .map(|mo_mirror_parities| {
815                                                    mirrors
816                                                        .iter()
817                                                        .map(|sigma| {
818                                                            let sigma_length =
819                                                                sigma.to_string().chars().count()
820                                                                    + 3;
821                                                            mo_mirror_parities
822                                                                .get(sigma)
823                                                                .map(|parity| {
824                                                                    format!(
825                                                                        "{:^sigma_length$}",
826                                                                        parity.to_string()
827                                                                    )
828                                                                })
829                                                                .unwrap_or_else(|| {
830                                                                    format!(
831                                                                        "{:^sigma_length$}",
832                                                                        "--"
833                                                                    )
834                                                                })
835                                                        })
836                                                        .join("  ")
837                                                })
838                                                .unwrap_or(String::new())
839                                        })
840                                })
841                        })
842                        .unwrap_or(String::new());
843
844                    let (eig_above_str, eig_below_str) = self
845                        .mo_symmetries_thresholds
846                        .as_ref()
847                        .map(|mo_symmetries_thresholds| {
848                            mo_symmetries_thresholds
849                                .get(spini)
850                                .and_then(|spin_mo_symmetries_thresholds| {
851                                    spin_mo_symmetries_thresholds.get(moi)
852                                })
853                                .map(|(eig_above_opt, eig_below_opt)| {
854                                    (
855                                        eig_above_opt
856                                            .map(|eig_above| format!("{eig_above:>+.3e}"))
857                                            .unwrap_or("--".to_string()),
858                                        eig_below_opt
859                                            .map(|eig_below| format!("{eig_below:>+.3e}"))
860                                            .unwrap_or("--".to_string()),
861                                    )
862                                })
863                                .unwrap_or(("--".to_string(), "--".to_string()))
864                        })
865                        .unwrap_or(("--".to_string(), "--".to_string()));
866
867                    let mo_symmetry_projections_str = mo_sym_projections_str_opt
868                        .as_ref()
869                        .and_then(|mo_symmetry_projectionss| {
870                            mo_symmetry_projectionss.get(spini).and_then(
871                                |spin_mo_symmetry_projections| {
872                                    spin_mo_symmetry_projections.get(moi)
873                                },
874                            )
875                        })
876                        .unwrap_or(&empty_string);
877
878                    let mo_density_symmetries_str = mo_den_symss_str_opt
879                        .as_ref()
880                        .and_then(|mo_density_symmetriess| {
881                            mo_density_symmetriess.get(spini).and_then(
882                                |spin_mo_density_symmetries| spin_mo_density_symmetries.get(moi),
883                            )
884                        })
885                        .unwrap_or(&empty_string);
886
887                    match (mo_density_length, mo_sym_projections_length) {
888                        (0, 0) => {
889                            writeln!(
890                                f,
891                                " {spini:>mo_spin_index_length$}  \
892                                {moi:>mo_index_length$}  \
893                                {occ_str:<mo_occ_length$}  \
894                                {mo_energy_str:<mo_energy_length$}  \
895                                {mo_sym_str:<mo_symmetry_length$}\
896                                {}{:mo_mirror_parities_length$}  \
897                                {eig_above_str:<mo_eig_above_length$}  \
898                                {eig_below_str}",
899                                " ".repeat(mo_mirror_parities_gap),
900                                mo_mirror_parities_str,
901                            )?;
902                        }
903                        (_, 0) => {
904                            writeln!(
905                                f,
906                                " {spini:>mo_spin_index_length$}  \
907                                {moi:>mo_index_length$}  \
908                                {occ_str:<mo_occ_length$}  \
909                                {mo_energy_str:<mo_energy_length$}  \
910                                {mo_sym_str:<mo_symmetry_length$}\
911                                {}{:mo_mirror_parities_length$}  \
912                                {eig_above_str:<mo_eig_above_length$}  \
913                                {eig_below_str:<mo_eig_below_length$}  \
914                                {mo_density_symmetries_str}",
915                                " ".repeat(mo_mirror_parities_gap),
916                                mo_mirror_parities_str,
917                            )?;
918                        }
919                        (0, _) => {
920                            writeln!(
921                                f,
922                                " {spini:>mo_spin_index_length$}  \
923                                {moi:>mo_index_length$}  \
924                                {occ_str:<mo_occ_length$}  \
925                                {mo_energy_str:<mo_energy_length$}  \
926                                {mo_sym_str:<mo_symmetry_length$}\
927                                {}{:mo_mirror_parities_length$}  \
928                                {eig_above_str:<mo_eig_above_length$}  \
929                                {eig_below_str:<mo_eig_below_length$}  \
930                                {mo_symmetry_projections_str}",
931                                " ".repeat(mo_mirror_parities_gap),
932                                mo_mirror_parities_str,
933                            )?;
934                        }
935                        (_, _) => {
936                            writeln!(
937                                f,
938                                " {spini:>mo_spin_index_length$}  \
939                                {moi:>mo_index_length$}  \
940                                {occ_str:<mo_occ_length$}  \
941                                {mo_energy_str:<mo_energy_length$}  \
942                                {mo_sym_str:<mo_symmetry_length$}\
943                                {}{:mo_mirror_parities_length$}  \
944                                {eig_above_str:<mo_eig_above_length$}  \
945                                {eig_below_str:<mo_eig_below_length$}  \
946                                {mo_symmetry_projections_str:<mo_sym_projections_length$}  \
947                                {mo_density_symmetries_str}",
948                                " ".repeat(mo_mirror_parities_gap),
949                                mo_mirror_parities_str,
950                            )?;
951                        }
952                    }
953                }
954            }
955
956            writeln!(f, "{}", "┈".repeat(table_width))?;
957            writeln!(f)?;
958        }
959
960        Ok(())
961    }
962}
963
964impl<'a, G, T, SC> fmt::Debug for SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
965where
966    G: SymmetryGroupProperties + Clone,
967    G::CharTab: SubspaceDecomposable<T>,
968    T: ComplexFloat + Lapack,
969    SC: StructureConstraint + fmt::Display,
970    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
971{
972    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
973        writeln!(f, "{self}")
974    }
975}
976
977// ------
978// Driver
979// ------
980
981// ~~~~~~~~~~~~~~~~~
982// Struct definition
983// ~~~~~~~~~~~~~~~~~
984
985/// Driver structure for performing representation analysis on Slater determinants.
986#[derive(Clone, Builder)]
987#[builder(build_fn(validate = "Self::validate"))]
988pub struct SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
989where
990    G: SymmetryGroupProperties + Clone,
991    G::CharTab: SubspaceDecomposable<T>,
992    T: ComplexFloat + Lapack,
993    SC: StructureConstraint + fmt::Display,
994    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
995{
996    /// The control parameters for Slater determinant representation analysis.
997    parameters: &'a SlaterDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
998
999    /// The Slater determinant to be analysed.
1000    determinant: &'a SlaterDeterminant<'a, T, SC>,
1001
1002    /// The result from symmetry-group detection on the underlying molecular structure of the
1003    /// Slater determinant.
1004    symmetry_group: &'a SymmetryGroupDetectionResult,
1005
1006    /// The atomic-orbital overlap matrix of the underlying basis set used to describe the
1007    /// determinant. This is either for a single component corresponding to the basis functions
1008    /// specified by the basis angular order structure in [`Self::determinant`], or for *all*
1009    /// explicit components specified by the coefficients in [`Self::determinant`].
1010    sao: &'a Array2<T>,
1011
1012    /// The complex-symmetric atomic-orbital overlap matrix of the underlying basis set used to
1013    /// describe the determinant. This is either for a single component corresponding to the basis
1014    /// functions specified by the basis angular order structure in [`Self::determinant`], or for
1015    /// *all* explicit components specified by the coefficients in [`Self::determinant`]. This is
1016    /// required if antiunitary symmetry operations are involved. If none is provided, this will be
1017    /// assumed to be the same as [`Self::sao`].
1018    #[builder(default = "None")]
1019    sao_h: Option<&'a Array2<T>>,
1020
1021    /// The atomic-orbital four-centre spatial overlap matrix of the underlying basis set used to
1022    /// describe the determinant. This is only required for density symmetry analysis.
1023    #[builder(default = "None")]
1024    sao_spatial_4c: Option<&'a Array4<T>>,
1025
1026    /// The complex-symmetric atomic-orbital four-centre spatial overlap matrix of the underlying
1027    /// basis set used to describe the determinant. This is only required for density symmetry
1028    /// analysis. This is required if antiunitary symmetry operations are involved. If none is
1029    /// provided, this will be assumed to be the same as [`Self::sao_spatial_4c`], if any.
1030    #[builder(default = "None")]
1031    sao_spatial_4c_h: Option<&'a Array4<T>>,
1032
1033    /// The control parameters for symmetry analysis of angular functions.
1034    angular_function_parameters: &'a AngularFunctionRepAnalysisParams,
1035
1036    /// The result of the Slater determinant representation analysis.
1037    #[builder(setter(skip), default = "None")]
1038    result: Option<SlaterDeterminantRepAnalysisResult<'a, G, T, SC>>,
1039}
1040
1041impl<'a, G, T, SC> SlaterDeterminantRepAnalysisDriverBuilder<'a, G, T, SC>
1042where
1043    G: SymmetryGroupProperties + Clone,
1044    G::CharTab: SubspaceDecomposable<T>,
1045    T: ComplexFloat + Lapack,
1046    SC: StructureConstraint + fmt::Display,
1047    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1048{
1049    fn validate(&self) -> Result<(), String> {
1050        let params = self
1051            .parameters
1052            .ok_or("No Slater determinant representation analysis parameters found.".to_string())?;
1053
1054        let sym_res = self
1055            .symmetry_group
1056            .ok_or("No symmetry group information found.".to_string())?;
1057
1058        let sao = self.sao.ok_or("No SAO matrix found.".to_string())?;
1059
1060        if let Some(sao_h) = self.sao_h.flatten() {
1061            if sao_h.shape() != sao.shape() {
1062                return Err("Mismatched shapes between `sao` and `sao_h`.".to_string());
1063            }
1064        }
1065
1066        match (
1067            self.sao_spatial_4c.flatten(),
1068            self.sao_spatial_4c_h.flatten(),
1069        ) {
1070            (Some(sao_spatial_4c), Some(sao_spatial_4c_h)) => {
1071                if sao_spatial_4c_h.shape() != sao_spatial_4c.shape() {
1072                    return Err(
1073                        "Mismatched shapes between `sao_spatial_4c` and `sao_spatial_4c_h`."
1074                            .to_string(),
1075                    );
1076                }
1077            }
1078            (None, Some(_)) => {
1079                return Err("`sao_spatial_4c_h` is provided without `sao_spatial_4c`.".to_string());
1080            }
1081            _ => {}
1082        }
1083
1084        let det = self
1085            .determinant
1086            .ok_or("No Slater determinant found.".to_string())?;
1087
1088        let sym = if params.use_magnetic_group.is_some() {
1089            sym_res
1090                .magnetic_symmetry
1091                .as_ref()
1092                .ok_or("Magnetic symmetry requested for representation analysis, but no magnetic symmetry found.")?
1093        } else {
1094            &sym_res.unitary_symmetry
1095        };
1096
1097        if sym.is_infinite() && params.infinite_order_to_finite.is_none() {
1098            Err(
1099                format!(
1100                    "Representation analysis cannot be performed using the entirety of the infinite group `{}`. \
1101                    Consider setting the parameter `infinite_order_to_finite` to restrict to a finite subgroup instead.",
1102                    sym.group_name.as_ref().expect("No symmetry group name found.")
1103                )
1104            )
1105        } else if (det.bao().n_funcs() != sao.nrows() || det.bao().n_funcs() != sao.ncols())
1106            && (det.bao().n_funcs()
1107                * det
1108                    .structure_constraint()
1109                    .n_explicit_comps_per_coefficient_matrix()
1110                != sao.nrows()
1111                || det.bao().n_funcs()
1112                    * det
1113                        .structure_constraint()
1114                        .n_explicit_comps_per_coefficient_matrix()
1115                    != sao.ncols())
1116        {
1117            Err("The dimensions of the SAO matrix do not match either the number of spatial AO basis functions or the number of spatial AO basis functions multiplied by the number of explicit components per coefficient matrix.".to_string())
1118        } else {
1119            Ok(())
1120        }
1121    }
1122}
1123
1124// ~~~~~~~~~~~~~~~~~~~~~~
1125// Struct implementations
1126// ~~~~~~~~~~~~~~~~~~~~~~
1127
1128// Generic for all symmetry groups G, determinant numeric type T, and structure constraint SC
1129// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1130
1131impl<'a, G, T, SC> SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1132where
1133    G: SymmetryGroupProperties + Clone,
1134    G::CharTab: SubspaceDecomposable<T>,
1135    T: ComplexFloat + Lapack,
1136    SC: StructureConstraint + Clone + fmt::Display,
1137    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1138{
1139    /// Returns a builder to construct a [`SlaterDeterminantRepAnalysisDriver`] structure.
1140    pub fn builder() -> SlaterDeterminantRepAnalysisDriverBuilder<'a, G, T, SC> {
1141        SlaterDeterminantRepAnalysisDriverBuilder::default()
1142    }
1143
1144    /// Constructs the appropriate atomic-orbital overlap matrix based on the structure constraint
1145    /// of the determinant and the specified overlap matrix.
1146    fn construct_sao(&self) -> Result<(Array2<T>, Option<Array2<T>>), anyhow::Error> {
1147        let nbas = self.determinant.bao().n_funcs();
1148        let ncomps = self
1149            .determinant
1150            .structure_constraint()
1151            .n_explicit_comps_per_coefficient_matrix();
1152        let provided_dim = self.sao.nrows();
1153
1154        if provided_dim == nbas {
1155            let sao = {
1156                let mut sao_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
1157                (0..ncomps).for_each(|icomp| {
1158                    let start = icomp * nbas;
1159                    let end = (icomp + 1) * nbas;
1160                    sao_mut
1161                        .slice_mut(s![start..end, start..end])
1162                        .assign(self.sao);
1163                });
1164                sao_mut
1165            };
1166
1167            let sao_h = self.sao_h.map(|sao_h| {
1168                let mut sao_h_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
1169                (0..ncomps).for_each(|icomp| {
1170                    let start = icomp * nbas;
1171                    let end = (icomp + 1) * nbas;
1172                    sao_h_mut
1173                        .slice_mut(s![start..end, start..end])
1174                        .assign(sao_h);
1175                });
1176                sao_h_mut
1177            });
1178
1179            Ok((sao, sao_h))
1180        } else {
1181            assert_eq!(provided_dim, nbas * ncomps);
1182            Ok((self.sao.clone(), self.sao_h.cloned()))
1183        }
1184    }
1185}
1186
1187// Specific for unitary-represented symmetry groups, but generic for determinant numeric type T and structure constraint SC
1188// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1189
1190impl<'a, T, SC> SlaterDeterminantRepAnalysisDriver<'a, UnitaryRepresentedSymmetryGroup, T, SC>
1191where
1192    T: ComplexFloat + Lapack + Sync + Send,
1193    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + Sync + Send,
1194    SC: StructureConstraint + fmt::Display,
1195    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
1196{
1197    fn_construct_unitary_group!(
1198        /// Constructs the unitary-represented group (which itself can be unitary or magnetic) ready
1199        /// for Slater determinant representation analysis.
1200        construct_unitary_group
1201    );
1202}
1203
1204// Specific for magnetic-represented symmetry groups, but generic for determinant numeric type T and structure constraint SC
1205// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1206
1207impl<'a, T, SC> SlaterDeterminantRepAnalysisDriver<'a, MagneticRepresentedSymmetryGroup, T, SC>
1208where
1209    T: ComplexFloat + Lapack + Sync + Send,
1210    <T as ComplexFloat>::Real: From<f64> + Sync + Send + fmt::LowerExp + fmt::Debug,
1211    SC: StructureConstraint + fmt::Display,
1212    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
1213{
1214    fn_construct_magnetic_group!(
1215        /// Constructs the magnetic-represented group (which itself can only be magnetic) ready for
1216        /// Slater determinant corepresentation analysis.
1217        construct_magnetic_group
1218    );
1219}
1220
1221// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric types f64 and C128
1222// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1223
1224#[duplicate_item(
1225    duplicate!{
1226        [
1227            dtype_nested sctype_nested;
1228            [ f64 ] [ SpinConstraint ];
1229            [ Complex<f64> ] [ SpinConstraint ];
1230            [ Complex<f64> ] [ SpinOrbitCoupled ];
1231        ]
1232        [
1233            gtype_ [ UnitaryRepresentedSymmetryGroup ]
1234            dtype_ [ dtype_nested ]
1235            sctype_ [ sctype_nested ]
1236            doc_sub_ [ "Performs representation analysis using a unitary-represented group and stores the result." ]
1237            analyse_fn_ [ analyse_representation ]
1238            construct_group_ [ self.construct_unitary_group()? ]
1239            calc_projections_ [
1240                log_subtitle("Slater determinant projection decompositions");
1241                qsym2_output!("");
1242                qsym2_output!("  Projections are defined w.r.t. the following inner product:");
1243                qsym2_output!("    {}", det_orbit.origin().overlap_definition());
1244                qsym2_output!("");
1245                det_orbit
1246                    .projections_to_string(
1247                        &det_orbit.calc_projection_compositions()?,
1248                        params.integrality_threshold,
1249                    )
1250                    .log_output_display();
1251                qsym2_output!("");
1252            ]
1253            calc_mo_projections_ [
1254                let mo_projectionss = mo_orbitss.iter().map(|mo_orbits| {
1255                    mo_orbits
1256                        .iter()
1257                        .map(|mo_orbit| mo_orbit.calc_projection_compositions().ok())
1258                        .collect::<Vec<_>>()
1259                }).collect::<Vec<_>>();
1260                Some(mo_projectionss)
1261            ]
1262        ]
1263    }
1264    duplicate!{
1265        [
1266            dtype_nested sctype_nested;
1267            [ f64 ] [ SpinConstraint ];
1268            [ Complex<f64> ] [ SpinConstraint ];
1269            [ Complex<f64> ] [ SpinOrbitCoupled ];
1270        ]
1271        [
1272            gtype_ [ MagneticRepresentedSymmetryGroup ]
1273            dtype_ [ dtype_nested ]
1274            sctype_ [ sctype_nested ]
1275            doc_sub_ [ "Performs corepresentation analysis using a magnetic-represented group and stores the result." ]
1276            analyse_fn_ [ analyse_corepresentation ]
1277            construct_group_ [ self.construct_magnetic_group()? ]
1278            calc_projections_ [ ]
1279            calc_mo_projections_ [
1280                None::<Vec<Vec<Option<Vec<(MullikenIrcorepSymbol, Complex<f64>)>>>>>
1281            ]
1282        ]
1283    }
1284)]
1285impl<'a> SlaterDeterminantRepAnalysisDriver<'a, gtype_, dtype_, sctype_> {
1286    #[doc = doc_sub_]
1287    fn analyse_fn_(&mut self) -> Result<(), anyhow::Error> {
1288        let params = self.parameters;
1289        let (sao, sao_h) = self.construct_sao()?;
1290        let group = construct_group_;
1291        log_cc_transversal(&group);
1292        let _ = find_angular_function_representation(&group, self.angular_function_parameters);
1293        if group.is_double_group() {
1294            let _ = find_spinor_function_representation(&group, self.angular_function_parameters);
1295        }
1296        log_bao(self.determinant.bao());
1297
1298        // Determinant and orbital symmetries
1299        let (
1300            det_symmetry,
1301            mo_symmetries,
1302            mo_symmetry_projections,
1303            mo_mirror_parities,
1304            mo_symmetries_thresholds,
1305        ) = if params.analyse_mo_symmetries {
1306            let mos = self.determinant.to_orbitals();
1307            let (mut det_orbit, mut mo_orbitss) = generate_det_mo_orbits(
1308                self.determinant,
1309                &mos,
1310                &group,
1311                &sao,
1312                sao_h.as_ref(),
1313                params.integrality_threshold,
1314                params.linear_independence_threshold,
1315                params.symmetry_transformation_kind.clone(),
1316                params.eigenvalue_comparison_mode.clone(),
1317                params.use_cayley_table,
1318            )?;
1319            det_orbit.calc_xmat(false)?;
1320            if params.write_overlap_eigenvalues {
1321                if let Some(smat_eigvals) = det_orbit.smat_eigvals.as_ref() {
1322                    log_overlap_eigenvalues(
1323                        "Determinant orbit overlap eigenvalues",
1324                        smat_eigvals,
1325                        params.linear_independence_threshold,
1326                        &params.eigenvalue_comparison_mode,
1327                    );
1328                    qsym2_output!("");
1329                }
1330            }
1331
1332            let det_symmetry = det_orbit.analyse_rep().map_err(|err| err.to_string());
1333
1334            {
1335                calc_projections_
1336            }
1337
1338            let mo_symmetries = mo_orbitss
1339                .iter_mut()
1340                .map(|mo_orbits| {
1341                    mo_orbits
1342                        .par_iter_mut()
1343                        .map(|mo_orbit| {
1344                            mo_orbit.calc_xmat(false).ok()?;
1345                            mo_orbit.analyse_rep().ok()
1346                        })
1347                        .collect::<Vec<_>>()
1348                })
1349                .collect::<Vec<_>>();
1350
1351            let mo_symmetry_projections_opt = if params.analyse_mo_symmetry_projections {
1352                calc_mo_projections_
1353            } else {
1354                None
1355            };
1356
1357            let mo_mirror_parities_opt = if params.analyse_mo_mirror_parities {
1358                Some(
1359                    mo_symmetries
1360                        .iter()
1361                        .map(|spin_mo_symmetries| {
1362                            spin_mo_symmetries
1363                                .iter()
1364                                .map(|mo_sym_opt| {
1365                                    mo_sym_opt.as_ref().map(|mo_sym| {
1366                                        deduce_mirror_parities(det_orbit.group(), mo_sym)
1367                                    })
1368                                })
1369                                .collect::<Vec<_>>()
1370                        })
1371                        .collect::<Vec<_>>(),
1372                )
1373            } else {
1374                None
1375            };
1376
1377            let mo_symmetries_thresholds = mo_orbitss
1378                .iter_mut()
1379                .map(|mo_orbits| {
1380                    mo_orbits
1381                        .par_iter_mut()
1382                        .map(|mo_orbit| {
1383                            mo_orbit
1384                                .smat_eigvals
1385                                .as_ref()
1386                                .map(|eigvals| {
1387                                    let mut eigvals_vec = eigvals.iter().collect::<Vec<_>>();
1388                                    match mo_orbit.eigenvalue_comparison_mode {
1389                                        EigenvalueComparisonMode::Modulus => {
1390                                            eigvals_vec.sort_by(|a, b| {
1391                                                a.abs().partial_cmp(&b.abs()).expect("Unable to compare two eigenvalues based on their moduli.")
1392                                            });
1393                                        }
1394                                        EigenvalueComparisonMode::Real => {
1395                                            eigvals_vec.sort_by(|a, b| {
1396                                                a.re().partial_cmp(&b.re()).expect("Unable to compare two eigenvalues based on their real parts.")
1397                                            });
1398                                        }
1399                                    }
1400                                    let eigval_above = match mo_orbit.eigenvalue_comparison_mode {
1401                                        EigenvalueComparisonMode::Modulus => eigvals_vec
1402                                            .iter()
1403                                            .find(|val| {
1404                                                val.abs() >= mo_orbit.linear_independence_threshold
1405                                            })
1406                                            .copied()
1407                                            .copied(),
1408                                        EigenvalueComparisonMode::Real => eigvals_vec
1409                                            .iter()
1410                                            .find(|val| {
1411                                                val.re() >= mo_orbit.linear_independence_threshold
1412                                            })
1413                                            .copied()
1414                                            .copied(),
1415                                    };
1416                                    eigvals_vec.reverse();
1417                                    let eigval_below = match mo_orbit.eigenvalue_comparison_mode {
1418                                        EigenvalueComparisonMode::Modulus => eigvals_vec
1419                                            .iter()
1420                                            .find(|val| {
1421                                                val.abs() < mo_orbit.linear_independence_threshold
1422                                            })
1423                                            .copied()
1424                                            .copied(),
1425                                        EigenvalueComparisonMode::Real => eigvals_vec
1426                                            .iter()
1427                                            .find(|val| {
1428                                                val.re() < mo_orbit.linear_independence_threshold
1429                                            })
1430                                            .copied()
1431                                            .copied(),
1432                                    };
1433                                    (eigval_above, eigval_below)
1434                                })
1435                                .unwrap_or((None, None))
1436                        })
1437                        .collect::<Vec<_>>()
1438                })
1439                .collect::<Vec<_>>();
1440            (
1441                det_symmetry,
1442                Some(mo_symmetries),
1443                mo_symmetry_projections_opt,
1444                mo_mirror_parities_opt,
1445                Some(mo_symmetries_thresholds),
1446            ) // det_symmetry, mo_symmetries, mo_symmetry_projections, mo_mirror_parities, mo_symmetries_thresholds
1447        } else {
1448            let mut det_orbit = SlaterDeterminantSymmetryOrbit::builder()
1449                .group(&group)
1450                .origin(self.determinant)
1451                .integrality_threshold(params.integrality_threshold)
1452                .linear_independence_threshold(params.linear_independence_threshold)
1453                .symmetry_transformation_kind(params.symmetry_transformation_kind.clone())
1454                .eigenvalue_comparison_mode(params.eigenvalue_comparison_mode.clone())
1455                .build()?;
1456            let det_symmetry = det_orbit
1457                .calc_smat(Some(&sao), sao_h.as_ref(), params.use_cayley_table)
1458                .and_then(|det_orb| det_orb.normalise_smat())
1459                .map_err(|err| err.to_string())
1460                .and_then(|det_orb| {
1461                    det_orb.calc_xmat(false).map_err(|err| err.to_string())?;
1462                    if params.write_overlap_eigenvalues {
1463                        if let Some(smat_eigvals) = det_orb.smat_eigvals.as_ref() {
1464                            log_overlap_eigenvalues(
1465                                "Determinant orbit overlap eigenvalues",
1466                                smat_eigvals,
1467                                params.linear_independence_threshold,
1468                                &params.eigenvalue_comparison_mode,
1469                            );
1470                            qsym2_output!("");
1471                        }
1472                    }
1473                    det_orb.analyse_rep().map_err(|err| err.to_string())
1474                });
1475
1476            {
1477                calc_projections_
1478            }
1479
1480            (det_symmetry, None, None, None, None) // det_symmetry, mo_symmetries, mo_symmetry_projections, mo_mirror_parities, mo_symmetries_thresholds
1481        };
1482
1483        // Density and orbital density symmetries
1484        let (den_symmetries, mo_den_symmetries) = if params.analyse_density_symmetries {
1485            let den_syms = self.determinant.to_densities().map(|densities| {
1486                let mut spin_den_syms = densities
1487                    .iter()
1488                    .enumerate()
1489                    .map(|(ispin, den)| {
1490                        let den_sym_res = || {
1491                            let mut den_orbit = DensitySymmetryOrbit::builder()
1492                                .group(&group)
1493                                .origin(den)
1494                                .integrality_threshold(params.integrality_threshold)
1495                                .linear_independence_threshold(params.linear_independence_threshold)
1496                                .symmetry_transformation_kind(
1497                                    params.symmetry_transformation_kind.clone(),
1498                                )
1499                                .eigenvalue_comparison_mode(
1500                                    params.eigenvalue_comparison_mode.clone(),
1501                                )
1502                                .build()?;
1503                            den_orbit
1504                                .calc_smat(
1505                                    self.sao_spatial_4c,
1506                                    self.sao_spatial_4c_h,
1507                                    params.use_cayley_table,
1508                                )?
1509                                .normalise_smat()?
1510                                .calc_xmat(false)?;
1511                            den_orbit.analyse_rep().map_err(|err| format_err!(err))
1512                        };
1513                        (
1514                            format!("Spin-{ispin} density"),
1515                            den_sym_res().map_err(|err| err.to_string()),
1516                        )
1517                    })
1518                    .collect::<Vec<_>>();
1519                let mut extra_den_syms = densities
1520                    .calc_extra_densities()
1521                    .expect("Unable to calculate extra densities.")
1522                    .iter()
1523                    .map(|(label, den)| {
1524                        let den_sym_res = || {
1525                            let mut den_orbit = DensitySymmetryOrbit::builder()
1526                                .group(&group)
1527                                .origin(den)
1528                                .integrality_threshold(params.integrality_threshold)
1529                                .linear_independence_threshold(params.linear_independence_threshold)
1530                                .symmetry_transformation_kind(
1531                                    params.symmetry_transformation_kind.clone(),
1532                                )
1533                                .eigenvalue_comparison_mode(
1534                                    params.eigenvalue_comparison_mode.clone(),
1535                                )
1536                                .build()?;
1537                            den_orbit
1538                                .calc_smat(
1539                                    self.sao_spatial_4c,
1540                                    self.sao_spatial_4c_h,
1541                                    params.use_cayley_table,
1542                                )?
1543                                .calc_xmat(false)?;
1544                            den_orbit.analyse_rep().map_err(|err| format_err!(err))
1545                        };
1546                        (label.clone(), den_sym_res().map_err(|err| err.to_string()))
1547                    })
1548                    .collect_vec();
1549                // let mut extra_den_syms = match self.determinant.structure_constraint() {
1550                //     SpinConstraint::Restricted(_) => {
1551                //         vec![("Total density".to_string(), spin_den_syms[0].1.clone())]
1552                //     }
1553                //     SpinConstraint::Unrestricted(nspins, _)
1554                //     | SpinConstraint::Generalised(nspins, _) => {
1555                //         let total_den_sym_res = || {
1556                //             let nspatial = self.determinant.bao().n_funcs();
1557                //             let zero_den = Density::<dtype_>::builder()
1558                //                 .density_matrix(Array2::<dtype_>::zeros((nspatial, nspatial)))
1559                //                 .bao(self.determinant.bao())
1560                //                 .mol(self.determinant.mol())
1561                //                 .complex_symmetric(self.determinant.complex_symmetric())
1562                //                 .threshold(self.determinant.threshold())
1563                //                 .build()?;
1564                //             let total_den =
1565                //                 densities.iter().fold(zero_den, |acc, denmat| acc + denmat);
1566                //             let mut total_den_orbit = DensitySymmetryOrbit::builder()
1567                //                 .group(&group)
1568                //                 .origin(&total_den)
1569                //                 .integrality_threshold(params.integrality_threshold)
1570                //                 .linear_independence_threshold(params.linear_independence_threshold)
1571                //                 .symmetry_transformation_kind(
1572                //                     params.symmetry_transformation_kind.clone(),
1573                //                 )
1574                //                 .eigenvalue_comparison_mode(
1575                //                     params.eigenvalue_comparison_mode.clone(),
1576                //                 )
1577                //                 .build()?;
1578                //             total_den_orbit
1579                //                 .calc_smat(
1580                //                     self.sao_spatial_4c,
1581                //                     self.sao_spatial_4c_h,
1582                //                     params.use_cayley_table,
1583                //                 )?
1584                //                 .calc_xmat(false)?;
1585                //             total_den_orbit
1586                //                 .analyse_rep()
1587                //                 .map_err(|err| format_err!(err))
1588                //         };
1589                //         let mut extra_syms = vec![(
1590                //             "Total density".to_string(),
1591                //             total_den_sym_res().map_err(|err| err.to_string()),
1592                //         )];
1593                //         extra_syms.extend((0..usize::from(*nspins)).combinations(2).map(
1594                //             |indices| {
1595                //                 let i = indices[0];
1596                //                 let j = indices[1];
1597                //                 let den_ij = &densities[i] - &densities[j];
1598                //                 let den_ij_sym_res = || {
1599                //                     let mut den_ij_orbit = DensitySymmetryOrbit::builder()
1600                //                         .group(&group)
1601                //                         .origin(&den_ij)
1602                //                         .integrality_threshold(params.integrality_threshold)
1603                //                         .linear_independence_threshold(
1604                //                             params.linear_independence_threshold,
1605                //                         )
1606                //                         .symmetry_transformation_kind(
1607                //                             params.symmetry_transformation_kind.clone(),
1608                //                         )
1609                //                         .eigenvalue_comparison_mode(
1610                //                             params.eigenvalue_comparison_mode.clone(),
1611                //                         )
1612                //                         .build()?;
1613                //                     den_ij_orbit
1614                //                         .calc_smat(
1615                //                             self.sao_spatial_4c,
1616                //                             self.sao_spatial_4c_h,
1617                //                             params.use_cayley_table,
1618                //                         )?
1619                //                         .calc_xmat(false)?;
1620                //                     den_ij_orbit.analyse_rep().map_err(|err| format_err!(err))
1621                //                 };
1622                //                 (
1623                //                     format!("Spin-polarised density {i} - {j}"),
1624                //                     den_ij_sym_res().map_err(|err| err.to_string()),
1625                //                 )
1626                //             },
1627                //         ));
1628                //         extra_syms
1629                //     }
1630                // };
1631                spin_den_syms.append(&mut extra_den_syms);
1632                spin_den_syms
1633            });
1634
1635            let mo_den_syms = if params.analyse_mo_symmetries {
1636                let mo_den_symmetries = self
1637                    .determinant
1638                    .to_orbitals()
1639                    .iter()
1640                    .map(|mos| {
1641                        mos.par_iter()
1642                            .map(|mo| {
1643                                let mo_den = mo.to_total_density().ok()?;
1644                                let mut mo_den_orbit = DensitySymmetryOrbit::builder()
1645                                    .group(&group)
1646                                    .origin(&mo_den)
1647                                    .integrality_threshold(params.integrality_threshold)
1648                                    .linear_independence_threshold(
1649                                        params.linear_independence_threshold,
1650                                    )
1651                                    .symmetry_transformation_kind(
1652                                        params.symmetry_transformation_kind.clone(),
1653                                    )
1654                                    .eigenvalue_comparison_mode(
1655                                        params.eigenvalue_comparison_mode.clone(),
1656                                    )
1657                                    .build()
1658                                    .ok()?;
1659                                log::debug!("Computing overlap matrix for an MO density orbit...");
1660                                mo_den_orbit
1661                                    .calc_smat(
1662                                        self.sao_spatial_4c,
1663                                        self.sao_spatial_4c_h,
1664                                        params.use_cayley_table,
1665                                    )
1666                                    .ok()?
1667                                    .normalise_smat()
1668                                    .ok()?
1669                                    .calc_xmat(false)
1670                                    .ok()?;
1671                                log::debug!(
1672                                    "Computing overlap matrix for an MO density orbit... Done."
1673                                );
1674                                mo_den_orbit.analyse_rep().ok()
1675                            })
1676                            .collect::<Vec<_>>()
1677                    })
1678                    .collect::<Vec<_>>();
1679                Some(mo_den_symmetries)
1680            } else {
1681                None
1682            };
1683
1684            (den_syms.ok(), mo_den_syms)
1685        } else {
1686            (None, None)
1687        };
1688
1689        let result = SlaterDeterminantRepAnalysisResult::builder()
1690            .parameters(params)
1691            .determinant(self.determinant)
1692            .group(group)
1693            .determinant_symmetry(det_symmetry)
1694            .determinant_density_symmetries(den_symmetries)
1695            .mo_symmetries(mo_symmetries)
1696            .mo_symmetry_projections(mo_symmetry_projections)
1697            .mo_mirror_parities(mo_mirror_parities)
1698            .mo_symmetries_thresholds(mo_symmetries_thresholds)
1699            .mo_density_symmetries(mo_den_symmetries)
1700            .build()?;
1701        self.result = Some(result);
1702
1703        Ok(())
1704    }
1705}
1706
1707// ~~~~~~~~~~~~~~~~~~~~~
1708// Trait implementations
1709// ~~~~~~~~~~~~~~~~~~~~~
1710
1711// Generic for all symmetry groups G and determinant numeric type T
1712// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1713
1714impl<'a, G, T, SC> fmt::Display for SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1715where
1716    G: SymmetryGroupProperties + Clone,
1717    G::CharTab: SubspaceDecomposable<T>,
1718    T: ComplexFloat + Lapack,
1719    SC: StructureConstraint + fmt::Display,
1720    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1721{
1722    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1723        write_title(f, "Slater Determinant Symmetry Analysis")?;
1724        writeln!(f)?;
1725        writeln!(f, "{}", self.parameters)?;
1726        Ok(())
1727    }
1728}
1729
1730impl<'a, G, T, SC> fmt::Debug for SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1731where
1732    G: SymmetryGroupProperties + Clone,
1733    G::CharTab: SubspaceDecomposable<T>,
1734    T: ComplexFloat + Lapack,
1735    SC: StructureConstraint + fmt::Display,
1736    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1737{
1738    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1739        writeln!(f, "{self}")
1740    }
1741}
1742
1743// Specific for unitary/magnetic-represented groups and determinant numeric type f64/Complex<f64>
1744// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1745
1746#[duplicate_item(
1747    duplicate!{
1748        [
1749            dtype_nested sctype_nested;
1750            [ f64 ] [ SpinConstraint ];
1751            [ Complex<f64> ] [ SpinConstraint ];
1752            [ Complex<f64> ] [ SpinOrbitCoupled ];
1753        ]
1754        [
1755            gtype_ [ UnitaryRepresentedSymmetryGroup ]
1756            dtype_ [ dtype_nested ]
1757            sctype_ [ sctype_nested ]
1758            analyse_fn_ [ analyse_representation ]
1759        ]
1760    }
1761    duplicate!{
1762        [
1763            dtype_nested sctype_nested;
1764            [ f64 ] [ SpinConstraint ];
1765            [ Complex<f64> ] [ SpinConstraint ];
1766            [ Complex<f64> ] [ SpinOrbitCoupled ];
1767        ]
1768        [
1769            gtype_ [ MagneticRepresentedSymmetryGroup ]
1770            dtype_ [ dtype_nested ]
1771            sctype_ [ sctype_nested ]
1772            analyse_fn_ [ analyse_corepresentation ]
1773        ]
1774    }
1775)]
1776impl<'a> QSym2Driver for SlaterDeterminantRepAnalysisDriver<'a, gtype_, dtype_, sctype_> {
1777    type Params = SlaterDeterminantRepAnalysisParams<f64>;
1778
1779    type Outcome = SlaterDeterminantRepAnalysisResult<'a, gtype_, dtype_, sctype_>;
1780
1781    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
1782        self.result.as_ref().ok_or_else(|| {
1783            format_err!("No Slater determinant representation analysis results found.")
1784        })
1785    }
1786
1787    fn run(&mut self) -> Result<(), anyhow::Error> {
1788        self.log_output_display();
1789        self.analyse_fn_()?;
1790        self.result()?.log_output_display();
1791        Ok(())
1792    }
1793}