qsym2/drivers/representation_analysis/
multideterminant.rs

1//! Driver for symmetry analysis of multi-determinantal wavefunctions.
2
3use std::collections::HashSet;
4use std::fmt;
5use std::hash::Hash;
6use std::ops::Mul;
7
8use anyhow::{self, bail, format_err};
9use derive_builder::Builder;
10use duplicate::duplicate_item;
11use ndarray::{s, Array2};
12use ndarray_linalg::types::Lapack;
13use num_complex::{Complex, ComplexFloat};
14use num_traits::Float;
15use serde::{Deserialize, Serialize};
16
17use crate::analysis::{
18    log_overlap_eigenvalues, EigenvalueComparisonMode, Overlap, ProjectionDecomposition,
19    RepAnalysis,
20};
21use crate::angmom::spinor_rotation_3d::{SpinConstraint, StructureConstraint};
22use crate::chartab::chartab_group::CharacterProperties;
23use crate::chartab::SubspaceDecomposable;
24use crate::drivers::representation_analysis::angular_function::{
25    find_angular_function_representation, find_spinor_function_representation,
26    AngularFunctionRepAnalysisParams,
27};
28use crate::drivers::representation_analysis::{
29    fn_construct_magnetic_group, fn_construct_unitary_group, log_bao, log_cc_transversal,
30    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
31};
32use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
33use crate::drivers::QSym2Driver;
34use crate::group::{GroupProperties, MagneticRepresentedGroup, UnitaryRepresentedGroup};
35use crate::io::format::{
36    log_micsec_begin, log_micsec_end, log_subtitle, nice_bool, qsym2_output, write_subtitle,
37    write_title, QSym2Output,
38};
39use crate::symmetry::symmetry_group::{
40    MagneticRepresentedSymmetryGroup, SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup,
41};
42use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
43use crate::target::determinant::SlaterDeterminant;
44use crate::target::noci::basis::{Basis, EagerBasis, OrbitBasis};
45use crate::target::noci::multideterminant::multideterminant_analysis::MultiDeterminantSymmetryOrbit;
46use crate::target::noci::multideterminant::MultiDeterminant;
47
48#[cfg(test)]
49#[path = "multideterminant_tests.rs"]
50mod multideterminant_tests;
51
52// ==================
53// Struct definitions
54// ==================
55
56// ----------
57// Parameters
58// ----------
59
60const fn default_true() -> bool {
61    true
62}
63const fn default_symbolic() -> Option<CharacterTableDisplay> {
64    Some(CharacterTableDisplay::Symbolic)
65}
66
67/// Structure containing control parameters for multi-determinantal wavefunction representation
68/// analysis.
69#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
70pub struct MultiDeterminantRepAnalysisParams<T: From<f64>> {
71    /// Threshold for checking if subspace multiplicities are integral.
72    pub integrality_threshold: T,
73
74    /// Threshold for determining zero eigenvalues in the orbit overlap matrix.
75    pub linear_independence_threshold: T,
76
77    /// Option indicating if the magnetic group is to be used for symmetry analysis, and if so,
78    /// whether unitary representations or unitary-antiunitary corepresentations should be used.
79    #[builder(default = "None")]
80    #[serde(default)]
81    pub use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
82
83    /// Boolean indicating if the double group is to be used for symmetry analysis.
84    #[builder(default = "false")]
85    #[serde(default)]
86    pub use_double_group: bool,
87
88    /// Boolean indicating if the Cayley table of the group, if available, should be used to speed
89    /// up the computation of orbit overlap matrices.
90    #[builder(default = "true")]
91    #[serde(default = "default_true")]
92    pub use_cayley_table: bool,
93
94    /// The kind of symmetry transformation to be applied on the reference multi-determinantal
95    /// wavefunction to generate the orbit for symmetry analysis.
96    #[builder(default = "SymmetryTransformationKind::Spatial")]
97    #[serde(default)]
98    pub symmetry_transformation_kind: SymmetryTransformationKind,
99
100    /// Option indicating if the character table of the group used for symmetry analysis is to be
101    /// printed out.
102    #[builder(default = "Some(CharacterTableDisplay::Symbolic)")]
103    #[serde(default = "default_symbolic")]
104    pub write_character_table: Option<CharacterTableDisplay>,
105
106    /// Boolean indicating if the eigenvalues of the orbit overlap matrix are to be printed out.
107    #[builder(default = "true")]
108    #[serde(default = "default_true")]
109    pub write_overlap_eigenvalues: bool,
110
111    /// The comparison mode for filtering out orbit overlap eigenvalues.
112    #[builder(default = "EigenvalueComparisonMode::Modulus")]
113    #[serde(default)]
114    pub eigenvalue_comparison_mode: EigenvalueComparisonMode,
115
116    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
117    /// subgroup of an infinite group can be used for the symmetry analysis.
118    #[builder(default = "None")]
119    #[serde(default)]
120    pub infinite_order_to_finite: Option<u32>,
121}
122
123impl<T> MultiDeterminantRepAnalysisParams<T>
124where
125    T: Float + From<f64>,
126{
127    /// Returns a builder to construct a [`MultiDeterminantRepAnalysisParams`] structure.
128    pub fn builder() -> MultiDeterminantRepAnalysisParamsBuilder<T> {
129        MultiDeterminantRepAnalysisParamsBuilder::default()
130    }
131}
132
133impl Default for MultiDeterminantRepAnalysisParams<f64> {
134    fn default() -> Self {
135        Self::builder()
136            .integrality_threshold(1e-7)
137            .linear_independence_threshold(1e-7)
138            .build()
139            .expect("Unable to construct a default `MultiDeterminantRepAnalysisParams<f64>`.")
140    }
141}
142
143impl<T> fmt::Display for MultiDeterminantRepAnalysisParams<T>
144where
145    T: From<f64> + fmt::LowerExp + fmt::Debug,
146{
147    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
148        writeln!(
149            f,
150            "Integrality threshold: {:.3e}",
151            self.integrality_threshold
152        )?;
153        writeln!(
154            f,
155            "Linear independence threshold: {:.3e}",
156            self.linear_independence_threshold
157        )?;
158        writeln!(
159            f,
160            "Orbit eigenvalue comparison mode: {}",
161            self.eigenvalue_comparison_mode
162        )?;
163        writeln!(
164            f,
165            "Write overlap eigenvalues: {}",
166            nice_bool(self.write_overlap_eigenvalues)
167        )?;
168        writeln!(f)?;
169        writeln!(
170            f,
171            "Use magnetic group for analysis: {}",
172            match self.use_magnetic_group {
173                None => "no",
174                Some(MagneticSymmetryAnalysisKind::Representation) =>
175                    "yes, using unitary representations",
176                Some(MagneticSymmetryAnalysisKind::Corepresentation) =>
177                    "yes, using magnetic corepresentations",
178            }
179        )?;
180        writeln!(
181            f,
182            "Use double group for analysis: {}",
183            nice_bool(self.use_double_group)
184        )?;
185        writeln!(
186            f,
187            "Use Cayley table for orbit overlap matrices: {}",
188            nice_bool(self.use_cayley_table)
189        )?;
190        if let Some(finite_order) = self.infinite_order_to_finite {
191            writeln!(f, "Infinite order to finite: {finite_order}")?;
192        }
193        writeln!(
194            f,
195            "Symmetry transformation kind: {}",
196            self.symmetry_transformation_kind
197        )?;
198        writeln!(f)?;
199        writeln!(
200            f,
201            "Write character table: {}",
202            if let Some(chartab_display) = self.write_character_table.as_ref() {
203                format!("yes, {}", chartab_display.to_string().to_lowercase())
204            } else {
205                "no".to_string()
206            }
207        )?;
208
209        Ok(())
210    }
211}
212
213// ------
214// Result
215// ------
216
217/// Structure to contain multi-determinantal representation analysis results.
218#[derive(Clone, Builder)]
219pub struct MultiDeterminantRepAnalysisResult<'a, G, T, B, SC>
220where
221    G: SymmetryGroupProperties + Clone,
222    G::CharTab: SubspaceDecomposable<T>,
223    T: ComplexFloat + Lapack,
224    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
225    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
226    SC: StructureConstraint + Hash + Eq + fmt::Display,
227{
228    /// The control parameters used to obtain this set of multi-determinantal wavefunction
229    /// representation analysis results.
230    parameters: &'a MultiDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
231
232    /// The multi-determinantal wavefunctions being analysed.
233    multidets: Vec<&'a MultiDeterminant<'a, T, B, SC>>,
234
235    /// The group used for the representation analysis.
236    group: G,
237
238    /// The deduced symmetries of the multi-determinantal wavefunctions.
239    multidet_symmetries:
240        Vec<Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>>,
241
242    /// The overlap eigenvalues above and below the linear independence threshold for each
243    /// multi-determinantal wavefunction symmetry deduction.
244    multidet_symmetries_thresholds: Vec<(Option<T>, Option<T>)>,
245}
246
247impl<'a, G, T, B, SC> MultiDeterminantRepAnalysisResult<'a, G, T, B, SC>
248where
249    G: SymmetryGroupProperties + Clone,
250    G::CharTab: SubspaceDecomposable<T>,
251    T: ComplexFloat + Lapack,
252    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
253    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
254    SC: StructureConstraint + Clone + Hash + Eq + fmt::Display,
255{
256    /// Returns a builder to construct a new [`MultiDeterminantRepAnalysisResultBuilder`]
257    /// structure.
258    fn builder() -> MultiDeterminantRepAnalysisResultBuilder<'a, G, T, B, SC> {
259        MultiDeterminantRepAnalysisResultBuilder::default()
260    }
261
262    /// Returns the multi-determinantal wavefunction symmetries obtained from the analysis result.
263    pub fn multidet_symmetries(
264        &self,
265    ) -> &Vec<Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>> {
266        &self.multidet_symmetries
267    }
268}
269
270impl<'a, G, T, B, SC> fmt::Display for MultiDeterminantRepAnalysisResult<'a, G, T, B, SC>
271where
272    G: SymmetryGroupProperties + Clone,
273    G::CharTab: SubspaceDecomposable<T>,
274    T: ComplexFloat + Lapack,
275    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
276    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
277    SC: StructureConstraint + Clone + Hash + Eq + fmt::Display,
278{
279    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
280        write_subtitle(f, "Orbit-based symmetry analysis results")?;
281        writeln!(f)?;
282        writeln!(
283            f,
284            "> Group: {} ({})",
285            self.group
286                .finite_subgroup_name()
287                .map(|subgroup_name| format!("{} > {}", self.group.name(), subgroup_name))
288                .unwrap_or(self.group.name()),
289            self.group.group_type().to_string().to_lowercase()
290        )?;
291        writeln!(f)?;
292
293        let multidet_index_length = usize::try_from(self.multidets.len().ilog10() + 1).unwrap_or(4);
294        let multidet_symmetry_length = self
295            .multidet_symmetries
296            .iter()
297            .map(|multidet_sym| {
298                multidet_sym
299                    .as_ref()
300                    .map(|sym| sym.to_string())
301                    .unwrap_or_else(|err| err.clone())
302                    .chars()
303                    .count()
304            })
305            .max()
306            .unwrap_or(0)
307            .max(8);
308        let multidet_energy_length = self
309            .multidets
310            .iter()
311            .map(|multidet| {
312                multidet
313                    .energy()
314                    .map(|v| format!("{v:+.7}").chars().count())
315                    .unwrap_or(2)
316            })
317            .max()
318            .unwrap_or(6)
319            .max(6);
320
321        let multidet_eig_above_length: usize = self
322            .multidet_symmetries_thresholds
323            .iter()
324            .map(|(above, _)| {
325                above
326                    .as_ref()
327                    .map(|eig| format!("{eig:+.3e}"))
328                    .unwrap_or("--".to_string())
329                    .chars()
330                    .count()
331            })
332            .max()
333            .unwrap_or(10)
334            .max(10);
335
336        let multidet_eig_below_length: usize = self
337            .multidet_symmetries_thresholds
338            .iter()
339            .map(|(_, below)| {
340                below
341                    .as_ref()
342                    .map(|eig| format!("{eig:+.3e}"))
343                    .unwrap_or("--".to_string())
344                    .chars()
345                    .count()
346            })
347            .max()
348            .unwrap_or(10)
349            .max(10);
350
351        let table_width = 10
352            + multidet_index_length
353            + multidet_energy_length
354            + multidet_symmetry_length
355            + multidet_eig_above_length
356            + multidet_eig_below_length;
357
358        writeln!(f, "> Multi-determinantal results")?;
359        writeln!(
360            f,
361            "  Structure constraint: {}",
362            self.multidets
363                .get(0)
364                .map(|multidet_0| multidet_0.structure_constraint().to_string().to_lowercase())
365                .unwrap_or("--".to_string())
366        )?;
367        writeln!(f, "{}", "┈".repeat(table_width))?;
368        writeln!(
369            f,
370            " {:>multidet_index_length$}  {:<multidet_energy_length$}  {:<multidet_symmetry_length$}  {:<multidet_eig_above_length$}  Eig. below",
371            "#",
372            "Energy",
373            "Symmetry",
374            "Eig. above",
375        )?;
376        writeln!(f, "{}", "┈".repeat(table_width))?;
377
378        for (multidet_i, multidet) in self.multidets.iter().enumerate() {
379            let multidet_energy_str = multidet
380                .energy()
381                .map(|multidet_energy| format!("{multidet_energy:>+multidet_energy_length$.7}"))
382                .unwrap_or("--".to_string());
383            let multidet_sym_str = self
384                .multidet_symmetries
385                .get(multidet_i)
386                .ok_or_else(|| format!("Unable to retrieve the symmetry of multideterminantal wavefunction index `{multidet_i}`."))
387                .and_then(|sym_res| sym_res.as_ref().map(|sym| sym.to_string()).map_err(|err| err.to_string()))
388                .unwrap_or_else(|err| err);
389
390            let (eig_above_str, eig_below_str) = self
391                .multidet_symmetries_thresholds
392                .get(multidet_i)
393                .map(|(eig_above_opt, eig_below_opt)| {
394                    (
395                        eig_above_opt
396                            .map(|eig_above| format!("{eig_above:>+.3e}"))
397                            .unwrap_or("--".to_string()),
398                        eig_below_opt
399                            .map(|eig_below| format!("{eig_below:>+.3e}"))
400                            .unwrap_or("--".to_string()),
401                    )
402                })
403                .unwrap_or(("--".to_string(), "--".to_string()));
404            writeln!(
405                f,
406                " {multidet_i:>multidet_index_length$}  \
407                {multidet_energy_str:<multidet_energy_length$}  \
408                {multidet_sym_str:<multidet_symmetry_length$}  \
409                {eig_above_str:<multidet_eig_above_length$}  \
410                {eig_below_str}"
411            )?;
412        }
413
414        writeln!(f, "{}", "┈".repeat(table_width))?;
415        writeln!(f)?;
416
417        Ok(())
418    }
419}
420
421impl<'a, G, T, B, SC> fmt::Debug for MultiDeterminantRepAnalysisResult<'a, G, T, B, SC>
422where
423    G: SymmetryGroupProperties + Clone,
424    G::CharTab: SubspaceDecomposable<T>,
425    T: ComplexFloat + Lapack,
426    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
427    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
428    SC: StructureConstraint + Clone + Hash + Eq + fmt::Display,
429{
430    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
431        writeln!(f, "{self}")
432    }
433}
434
435// ------
436// Driver
437// ------
438
439// ~~~~~~~~~~~~~~~~~
440// Struct definition
441// ~~~~~~~~~~~~~~~~~
442
443/// Driver structure for performing representation analysis on multi-determinantal wavefunctions.
444#[derive(Clone, Builder)]
445#[builder(build_fn(validate = "Self::validate"))]
446pub struct MultiDeterminantRepAnalysisDriver<'a, G, T, B, SC>
447where
448    G: SymmetryGroupProperties + Clone,
449    G::CharTab: SubspaceDecomposable<T>,
450    T: ComplexFloat + Lapack,
451    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
452    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
453    SC: StructureConstraint + Hash + Eq + fmt::Display,
454{
455    /// The control parameters for multi-determinantal wavefunction representation analysis.
456    parameters: &'a MultiDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
457
458    /// The multi-determinantal wavefunctions to be analysed.
459    multidets: Vec<&'a MultiDeterminant<'a, T, B, SC>>,
460
461    /// The result from symmetry-group detection on the underlying molecular structure of the
462    /// multi-determinantal wavefunctions.
463    symmetry_group: &'a SymmetryGroupDetectionResult,
464
465    /// The atomic-orbital overlap matrix of the underlying basis set used to describe the
466    /// wavefunctions. This is either for a single component corresponding to the basis functions
467    /// specified by the basis angular order structure in the determinants in [`Self::multidets`],
468    /// or for *all* explicit components specified by the coefficients in the determinants.
469    sao: &'a Array2<T>,
470
471    /// The complex-symmetric atomic-orbital overlap matrix of the underlying basis set used to
472    /// describe the wavefunctions. This is either for a single component corresponding to the basis
473    /// functions specified by the basis angular order structure in the determinants, or for *all*
474    /// explicit components specified by the coefficients in the determinants.This is required if
475    /// antiunitary symmetry operations are involved. If none is provided, this will be assumed to
476    /// be the same as [`Self::sao`].
477    #[builder(default = "None")]
478    sao_h: Option<&'a Array2<T>>,
479
480    /// The control parameters for symmetry analysis of angular functions.
481    angular_function_parameters: &'a AngularFunctionRepAnalysisParams,
482
483    /// The result of the multi-determinantal wavefunction representation analysis.
484    #[builder(setter(skip), default = "None")]
485    result: Option<MultiDeterminantRepAnalysisResult<'a, G, T, B, SC>>,
486}
487
488impl<'a, G, T, B, SC> MultiDeterminantRepAnalysisDriverBuilder<'a, G, T, B, SC>
489where
490    G: SymmetryGroupProperties + Clone,
491    G::CharTab: SubspaceDecomposable<T>,
492    T: ComplexFloat + Lapack,
493    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
494    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
495    SC: StructureConstraint + Hash + Eq + fmt::Display,
496{
497    fn validate(&self) -> Result<(), String> {
498        let params = self.parameters.ok_or(
499            "No multi-determinantal wavefunction representation analysis parameters found."
500                .to_string(),
501        )?;
502
503        let sym_res = self
504            .symmetry_group
505            .ok_or("No symmetry group information found.".to_string())?;
506
507        let sao = self.sao.ok_or("No SAO matrix found.".to_string())?;
508
509        if let Some(sao_h) = self.sao_h.flatten() {
510            if sao_h.shape() != sao.shape() {
511                return Err("Mismatched shapes between `sao` and `sao_h`.".to_string());
512            }
513        }
514
515        let multidets = self
516            .multidets
517            .as_ref()
518            .ok_or("No multi-determinantal wavefunctions found.".to_string())?;
519        let mut n_spatial_set = multidets
520            .iter()
521            .flat_map(|multidet| {
522                multidet.basis().iter().map(|det_res| {
523                    det_res.map(|det| {
524                        (
525                            det.bao().n_funcs(),
526                            det.structure_constraint()
527                                .n_explicit_comps_per_coefficient_matrix(),
528                        )
529                    })
530                })
531            })
532            .collect::<Result<HashSet<(usize, usize)>, _>>()
533            .map_err(|err| err.to_string())?;
534        let (n_spatial, n_comps) = if n_spatial_set.len() == 1 {
535            n_spatial_set
536                .drain()
537                .next()
538                .ok_or("Unable to retrieve the number of spatial AO basis functions and the number of explicit components.".to_string())
539        } else {
540            Err("Inconsistent numbers of spatial AO basis functions and/or explicit components across multi-determinantal wavefunctions.".to_string())
541        }?;
542
543        let sym = if params.use_magnetic_group.is_some() {
544            sym_res
545                .magnetic_symmetry
546                .as_ref()
547                .ok_or("Magnetic symmetry requested for representation analysis, but no magnetic symmetry found.")?
548        } else {
549            &sym_res.unitary_symmetry
550        };
551
552        if sym.is_infinite() && params.infinite_order_to_finite.is_none() {
553            Err(
554                format!(
555                    "Representation analysis cannot be performed using the entirety of the infinite group `{}`. \
556                    Consider setting the parameter `infinite_order_to_finite` to restrict to a finite subgroup instead.",
557                    sym.group_name.as_ref().expect("No symmetry group name found.")
558                )
559            )
560        } else if (n_spatial != sao.nrows() || n_spatial != sao.ncols())
561            || (n_spatial * n_comps != sao.nrows() || n_spatial * n_comps != sao.ncols())
562        {
563            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())
564        } else {
565            Ok(())
566        }
567    }
568}
569
570// ~~~~~~~~~~~~~~~~~~~~~~
571// Struct implementations
572// ~~~~~~~~~~~~~~~~~~~~~~
573
574// Generic for all symmetry groups G and wavefunction numeric type T
575// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
576
577impl<'a, G, T, B, SC> MultiDeterminantRepAnalysisDriver<'a, G, T, B, SC>
578where
579    G: SymmetryGroupProperties + Clone,
580    G::CharTab: SubspaceDecomposable<T>,
581    T: ComplexFloat + Lapack,
582    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
583    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
584    SC: StructureConstraint + Clone + Hash + Eq + fmt::Display,
585{
586    /// Returns a builder to construct a [`MultiDeterminantRepAnalysisDriver`] structure.
587    pub fn builder() -> MultiDeterminantRepAnalysisDriverBuilder<'a, G, T, B, SC> {
588        MultiDeterminantRepAnalysisDriverBuilder::default()
589    }
590
591    /// Constructs the appropriate atomic-orbital overlap matrix based on the structure constraint of
592    /// the multi-determinantal wavefunctions and the provided overlap matrix.
593    fn construct_sao(&self) -> Result<(Array2<T>, Option<Array2<T>>), anyhow::Error> {
594        let mut structure_constraint_set = self
595            .multidets
596            .iter()
597            .map(|multidet| multidet.structure_constraint())
598            .collect::<HashSet<_>>();
599        let structure_constraint = if structure_constraint_set.len() == 1 {
600            structure_constraint_set.drain().next().ok_or(format_err!(
601                "Unable to retrieve the structure constraint of the multi-determinantal wavefunctions."
602            ))
603        } else {
604            Err(format_err!(
605                "Inconsistent structure constraints across multi-determinantal wavefunctions."
606            ))
607        }?;
608
609        let mut nbas_set = self
610            .multidets
611            .iter()
612            .map(|multidet| {
613                multidet
614                    .basis()
615                    .first()
616                    .expect("Unable to obtain the first determinant in the basis.")
617                    .bao()
618                    .n_funcs()
619            })
620            .collect::<HashSet<_>>();
621        let nbas = if nbas_set.len() == 1 {
622            nbas_set.drain().next().ok_or(format_err!(
623                "Unable to retrieve the number of basis functions describing the multi-determinantal wavefunctions."
624            ))
625        } else {
626            Err(format_err!(
627                "Inconsistent numbers of basis functions across multi-determinantal wavefunctions."
628            ))
629        }?;
630        let ncomps = structure_constraint.n_explicit_comps_per_coefficient_matrix();
631        let provided_dim = self.sao.nrows();
632
633        if provided_dim == nbas {
634            let sao = {
635                let mut sao_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
636                (0..ncomps).for_each(|icomp| {
637                    let start = icomp * nbas;
638                    let end = (icomp + 1) * nbas;
639                    sao_mut
640                        .slice_mut(s![start..end, start..end])
641                        .assign(self.sao);
642                });
643                sao_mut
644            };
645
646            let sao_h = self.sao_h.map(|sao_h| {
647                let mut sao_h_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
648                (0..ncomps).for_each(|icomp| {
649                    let start = icomp * nbas;
650                    let end = (icomp + 1) * nbas;
651                    sao_h_mut
652                        .slice_mut(s![start..end, start..end])
653                        .assign(sao_h);
654                });
655                sao_h_mut
656            });
657
658            Ok((sao, sao_h))
659        } else {
660            assert_eq!(provided_dim, nbas * ncomps);
661            Ok((self.sao.clone(), self.sao_h.cloned()))
662        }
663    }
664}
665
666// Specific for unitary-represented symmetry groups, but generic for wavefunction numeric type T
667// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
668
669impl<'a, T, B, SC> MultiDeterminantRepAnalysisDriver<'a, UnitaryRepresentedSymmetryGroup, T, B, SC>
670where
671    T: ComplexFloat + Lapack + Sync + Send,
672    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + Sync + Send,
673    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
674    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
675    SC: StructureConstraint + Hash + Eq + fmt::Display,
676{
677    fn_construct_unitary_group!(
678        /// Constructs the unitary-represented group (which itself can be unitary or magnetic) ready
679        /// for Slater determinant representation analysis.
680        construct_unitary_group
681    );
682}
683
684// Specific for magnetic-represented symmetry groups, but generic for wavefunction numeric type T
685// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
686
687impl<'a, T, B, SC> MultiDeterminantRepAnalysisDriver<'a, MagneticRepresentedSymmetryGroup, T, B, SC>
688where
689    T: ComplexFloat + Lapack + Sync + Send,
690    <T as ComplexFloat>::Real: From<f64> + Sync + Send + fmt::LowerExp + fmt::Debug,
691    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
692    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
693    SC: StructureConstraint + Hash + Eq + fmt::Display,
694{
695    fn_construct_magnetic_group!(
696        /// Constructs the magnetic-represented group (which itself can only be magnetic) ready for
697        /// Slater determinant corepresentation analysis.
698        construct_magnetic_group
699    );
700}
701
702// Specific for unitary-represented and magnetic-represented symmetry groups and wavefunction numeric types f64 and C128
703// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
704
705#[duplicate_item(
706    duplicate!{
707        [ dtype_nested; [f64]; [Complex<f64>] ]
708        duplicate!{
709            [ sctype_nested; [SpinConstraint] ]
710            duplicate!{
711                [
712                    [
713                        btype_nested [OrbitBasis<'a, UnitaryRepresentedSymmetryGroup, SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
714                        calc_smat_nested [calc_smat_optimised]
715                    ]
716                    [
717                        btype_nested [EagerBasis<SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
718                        calc_smat_nested [calc_smat]
719                    ]
720                ]
721                [
722                    gtype_ [ UnitaryRepresentedSymmetryGroup ]
723                    dtype_ [ dtype_nested ]
724                    btype_ [ btype_nested ]
725                    sctype_ [ sctype_nested ]
726                    doc_sub_ [ "Performs representation analysis using a unitary-represented group and stores the result." ]
727                    analyse_fn_ [ analyse_representation ]
728                    construct_group_ [ self.construct_unitary_group()? ]
729                    calc_smat_ [ calc_smat_nested ]
730                    calc_projections_ [
731                        log_subtitle("Multi-determinantal wavefunction projection decompositions");
732                        qsym2_output!("");
733                        qsym2_output!("  Projections are defined w.r.t. the following inner product:");
734                        qsym2_output!("    {}", multidet_orbit.origin().overlap_definition());
735                        qsym2_output!("");
736                        multidet_orbit
737                            .projections_to_string(
738                                &multidet_orbit.calc_projection_compositions()?,
739                                params.integrality_threshold,
740                            )
741                            .log_output_display();
742                        qsym2_output!("");
743                    ]
744                ]
745            }
746        }
747    }
748    duplicate!{
749        [ dtype_nested; [f64]; [Complex<f64>] ]
750        duplicate!{
751            [ sctype_nested; [SpinConstraint] ]
752            duplicate!{
753                [
754                    [
755                        btype_nested [OrbitBasis<'a, MagneticRepresentedSymmetryGroup, SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
756                        calc_smat_nested [calc_smat_optimised]
757                    ]
758                    [
759                        btype_nested [EagerBasis<SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
760                        calc_smat_nested [calc_smat]
761                    ]
762                ]
763                [
764                    gtype_ [ MagneticRepresentedSymmetryGroup ]
765                    dtype_ [ dtype_nested ]
766                    btype_ [ btype_nested ]
767                    sctype_ [ sctype_nested ]
768                    doc_sub_ [ "Performs corepresentation analysis using a magnetic-represented group and stores the result." ]
769                    analyse_fn_ [ analyse_corepresentation ]
770                    construct_group_ [ self.construct_magnetic_group()? ]
771                    calc_smat_ [ calc_smat_nested ]
772                    calc_projections_ [ ]
773                ]
774            }
775    }
776    }
777)]
778impl<'a> MultiDeterminantRepAnalysisDriver<'a, gtype_, dtype_, btype_, sctype_> {
779    #[doc = doc_sub_]
780    fn analyse_fn_(&mut self) -> Result<(), anyhow::Error> {
781        let params = self.parameters;
782        let (sao, sao_h) = self.construct_sao()?;
783        let group = construct_group_;
784        log_cc_transversal(&group);
785        let _ = find_angular_function_representation(&group, self.angular_function_parameters);
786        if group.is_double_group() {
787            let _ = find_spinor_function_representation(&group, self.angular_function_parameters);
788        }
789        if let Some(det) = self
790            .multidets
791            .get(0)
792            .and_then(|multidet| multidet.basis().first())
793        {
794            log_bao(det.bao());
795        }
796
797        let (multidet_symmetries, multidet_symmetries_thresholds): (Vec<_>, Vec<_>) = self.multidets
798            .iter()
799            .enumerate()
800            .map(|(i, multidet)| {
801                log_micsec_begin(&format!("Multi-determinantal wavefunction {i}"));
802                qsym2_output!("");
803                let res = MultiDeterminantSymmetryOrbit::builder()
804                    .group(&group)
805                    .origin(multidet)
806                    .integrality_threshold(params.integrality_threshold)
807                    .linear_independence_threshold(params.linear_independence_threshold)
808                    .symmetry_transformation_kind(params.symmetry_transformation_kind.clone())
809                    .eigenvalue_comparison_mode(params.eigenvalue_comparison_mode.clone())
810                    .build()
811                    .map_err(|err| format_err!(err))
812                    .and_then(|mut multidet_orbit| {
813                        multidet_orbit
814                            .calc_smat_(Some(&sao), sao_h.as_ref(), params.use_cayley_table)?
815                            .normalise_smat()?
816                            .calc_xmat(false)?;
817                        log_overlap_eigenvalues(
818                            "Overlap eigenvalues",
819                            multidet_orbit.smat_eigvals.as_ref().ok_or(format_err!("Orbit overlap eigenvalues not found."))?,
820                            params.linear_independence_threshold,
821                            &params.eigenvalue_comparison_mode
822                        );
823                        qsym2_output!("");
824                        let multidet_symmetry_thresholds = multidet_orbit
825                            .smat_eigvals
826                            .as_ref()
827                            .map(|eigvals| {
828                                let mut eigvals_vec = eigvals.iter().collect::<Vec<_>>();
829                                match multidet_orbit.eigenvalue_comparison_mode() {
830                                    EigenvalueComparisonMode::Modulus => {
831                                        eigvals_vec.sort_by(|a, b| {
832                                            a.abs().partial_cmp(&b.abs()).expect("Unable to compare two eigenvalues based on their moduli.")
833                                        });
834                                    }
835                                    EigenvalueComparisonMode::Real => {
836                                        eigvals_vec.sort_by(|a, b| {
837                                            a.re().partial_cmp(&b.re()).expect("Unable to compare two eigenvalues based on their real parts.")
838                                        });
839                                    }
840                                }
841                                let eigval_above = match multidet_orbit.eigenvalue_comparison_mode() {
842                                    EigenvalueComparisonMode::Modulus => eigvals_vec
843                                        .iter()
844                                        .find(|val| {
845                                            val.abs() >= multidet_orbit.linear_independence_threshold
846                                        })
847                                        .copied()
848                                        .copied(),
849                                    EigenvalueComparisonMode::Real => eigvals_vec
850                                        .iter()
851                                        .find(|val| {
852                                            val.re() >= multidet_orbit.linear_independence_threshold
853                                        })
854                                        .copied()
855                                        .copied(),
856                                };
857                                eigvals_vec.reverse();
858                                let eigval_below = match multidet_orbit.eigenvalue_comparison_mode() {
859                                    EigenvalueComparisonMode::Modulus => eigvals_vec
860                                        .iter()
861                                        .find(|val| {
862                                            val.abs() < multidet_orbit.linear_independence_threshold
863                                        })
864                                        .copied()
865                                        .copied(),
866                                    EigenvalueComparisonMode::Real => eigvals_vec
867                                        .iter()
868                                        .find(|val| {
869                                            val.re() < multidet_orbit.linear_independence_threshold
870                                        })
871                                        .copied()
872                                        .copied(),
873                                };
874                                (eigval_above, eigval_below)
875                            })
876                            .unwrap_or((None, None));
877                        let multidet_sym = multidet_orbit.analyse_rep().map_err(|err| err.to_string());
878                        { calc_projections_ }
879                        Ok((multidet_sym, multidet_symmetry_thresholds))
880                    })
881                    .unwrap_or_else(|err| (Err(err.to_string()), (None, None)));
882                log_micsec_end(&format!("Multi-determinantal wavefunction {i}"));
883                qsym2_output!("");
884                res
885            }).unzip();
886
887        let result = MultiDeterminantRepAnalysisResult::builder()
888            .parameters(params)
889            .multidets(self.multidets.clone())
890            .group(group)
891            .multidet_symmetries(multidet_symmetries)
892            .multidet_symmetries_thresholds(multidet_symmetries_thresholds)
893            .build()?;
894        self.result = Some(result);
895
896        Ok(())
897    }
898}
899
900// ~~~~~~~~~~~~~~~~~~~~~
901// Trait implementations
902// ~~~~~~~~~~~~~~~~~~~~~
903
904// Generic for all symmetry groups G, basis B, and wavefunction numeric type T
905// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
906
907impl<'a, G, T, B, SC> fmt::Display for MultiDeterminantRepAnalysisDriver<'a, G, T, B, SC>
908where
909    G: SymmetryGroupProperties + Clone,
910    G::CharTab: SubspaceDecomposable<T>,
911    T: ComplexFloat + Lapack,
912    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
913    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
914    SC: StructureConstraint + Hash + Eq + fmt::Display,
915{
916    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
917        write_title(f, "Multi-determinantal Wavefunction Symmetry Analysis")?;
918        writeln!(f)?;
919        writeln!(f, "{}", self.parameters)?;
920        Ok(())
921    }
922}
923
924impl<'a, G, T, B, SC> fmt::Debug for MultiDeterminantRepAnalysisDriver<'a, G, T, B, SC>
925where
926    G: SymmetryGroupProperties + Clone,
927    G::CharTab: SubspaceDecomposable<T>,
928    T: ComplexFloat + Lapack,
929    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
930    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
931    SC: StructureConstraint + Hash + Eq + fmt::Display,
932{
933    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
934        writeln!(f, "{self}")
935    }
936}
937
938// Specific for unitary/magnetic-represented groups and determinant numeric type f64/Complex<f64>
939// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
940
941#[duplicate_item(
942    duplicate!{
943        [ dtype_nested; [f64]; [Complex<f64>] ]
944        duplicate!{
945            [ sctype_nested; [SpinConstraint] ]
946            duplicate!{
947                [
948                    btype_nested;
949                    [OrbitBasis<'a, UnitaryRepresentedSymmetryGroup, SlaterDeterminant<'a, dtype_nested, sctype_nested>>];
950                    [EagerBasis<SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
951                ]
952                [
953                    gtype_ [ UnitaryRepresentedSymmetryGroup ]
954                    dtype_ [ dtype_nested ]
955                    btype_ [ btype_nested ]
956                    sctype_ [ sctype_nested ]
957                    analyse_fn_ [ analyse_representation ]
958                ]
959            }
960        }
961    }
962    duplicate!{
963        [ dtype_nested; [f64]; [Complex<f64>] ]
964        duplicate!{
965            [ sctype_nested; [SpinConstraint] ]
966            duplicate!{
967                [
968                    btype_nested;
969                    [OrbitBasis<'a, MagneticRepresentedSymmetryGroup, SlaterDeterminant<'a, dtype_nested, sctype_nested>>];
970                    [EagerBasis<SlaterDeterminant<'a, dtype_nested, sctype_nested>>]
971                ]
972                [
973                    gtype_ [ MagneticRepresentedSymmetryGroup ]
974                    dtype_ [ dtype_nested ]
975                    btype_ [ btype_nested ]
976                    sctype_ [ sctype_nested ]
977                    analyse_fn_ [ analyse_corepresentation ]
978                ]
979            }
980        }
981    }
982)]
983impl<'a> QSym2Driver for MultiDeterminantRepAnalysisDriver<'a, gtype_, dtype_, btype_, sctype_> {
984    type Params = MultiDeterminantRepAnalysisParams<f64>;
985
986    type Outcome = MultiDeterminantRepAnalysisResult<'a, gtype_, dtype_, btype_, sctype_>;
987
988    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
989        self.result.as_ref().ok_or_else(|| {
990            format_err!(
991                "No multi-determinantal wavefunction representation analysis results found."
992            )
993        })
994    }
995
996    fn run(&mut self) -> Result<(), anyhow::Error> {
997        self.log_output_display();
998        self.analyse_fn_()?;
999        self.result()?.log_output_display();
1000        Ok(())
1001    }
1002}