qsym2/target/orbital/
orbital_analysis.rs

1//! Implementation of symmetry analysis for orbitals.
2
3use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, ensure, format_err, Context};
7use approx;
8use derive_builder::Builder;
9use itertools::{izip, Itertools};
10use log;
11use ndarray::{s, Array1, Array2, Axis, Ix2};
12use ndarray_linalg::{
13    eig::Eig,
14    eigh::Eigh,
15    solve::Determinant,
16    types::{Lapack, Scalar},
17    UPLO,
18};
19use num_complex::{Complex, ComplexFloat};
20use num_traits::{Float, ToPrimitive, Zero};
21
22use crate::analysis::{
23    fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
24    Overlap, RepAnalysis,
25};
26use crate::angmom::spinor_rotation_3d::StructureConstraint;
27use crate::auxiliary::misc::{complex_modified_gram_schmidt, ProductRepeat};
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::group::GroupType;
31use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
32use crate::symmetry::symmetry_group::SymmetryGroupProperties;
33use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
34use crate::target::determinant::determinant_analysis::SlaterDeterminantSymmetryOrbit;
35use crate::target::determinant::SlaterDeterminant;
36use crate::target::orbital::MolecularOrbital;
37
38// -------
39// Overlap
40// -------
41
42impl<'a, T, SC> Overlap<T, Ix2> for MolecularOrbital<'a, T, SC>
43where
44    T: Lapack
45        + ComplexFloat<Real = <T as Scalar>::Real>
46        + fmt::Debug
47        + Mul<<T as ComplexFloat>::Real, Output = T>,
48    <T as ComplexFloat>::Real: fmt::Debug
49        + approx::RelativeEq<<T as ComplexFloat>::Real>
50        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
51    SC: StructureConstraint + Eq + fmt::Display,
52{
53    fn complex_symmetric(&self) -> bool {
54        self.complex_symmetric
55    }
56
57    /// Computes the overlap between two molecular orbitals.
58    ///
59    /// When one or both of the orbitals have been acted on by an antiunitary operation, the correct
60    /// Hermitian or complex-symmetric metric will be chosen in the evalulation of the overlap.
61    ///
62    /// # Panics
63    ///
64    /// Panics if `self` and `other` have mismatched spin constraints or coefficient array lengths.
65    fn overlap(
66        &self,
67        other: &Self,
68        metric: Option<&Array2<T>>,
69        metric_h: Option<&Array2<T>>,
70    ) -> Result<T, anyhow::Error> {
71        ensure!(
72            self.structure_constraint == other.structure_constraint,
73            "Inconsistent spin constraints between `self` and `other`."
74        );
75        ensure!(
76            self.coefficients.len() == other.coefficients.len(),
77            "Inconsistent numbers of coefficient matrices between `self` and `other`."
78        );
79
80        let sao = metric.ok_or_else(|| format_err!("No atomic-orbital metric found."))?;
81        let sao_h = metric_h.unwrap_or(sao);
82        let ov = if self.component_index != other.component_index {
83            T::zero()
84        } else if self.complex_symmetric() {
85            match (self.complex_conjugated, other.complex_conjugated) {
86                (false, false) => self.coefficients.t().dot(sao_h).dot(&other.coefficients),
87                (true, false) => self.coefficients.t().dot(sao).dot(&other.coefficients),
88                (false, true) => other.coefficients.t().dot(sao).dot(&self.coefficients),
89                (true, true) => self
90                    .coefficients
91                    .t()
92                    .dot(&sao_h.t())
93                    .dot(&other.coefficients),
94            }
95        } else {
96            match (self.complex_conjugated, other.complex_conjugated) {
97                (false, false) => self
98                    .coefficients
99                    .t()
100                    .mapv(|x| x.conj()) // This conjugation is still needed as it comes from the sesquilinear form.
101                    .dot(sao)
102                    .dot(&other.coefficients),
103                (true, false) => self
104                    .coefficients
105                    .t()
106                    .mapv(|x| x.conj()) // This conjugation is still needed as it comes from the sesquilinear form.
107                    .dot(sao_h)
108                    .dot(&other.coefficients),
109                (false, true) => other
110                    .coefficients
111                    .t()
112                    .mapv(|x| x.conj()) // This conjugation is still needed as it comes from the sesquilinear form.
113                    .dot(sao_h)
114                    .dot(&self.coefficients)
115                    .conj(),
116                (true, true) => self
117                    .coefficients
118                    .t()
119                    .mapv(|x| x.conj()) // This conjugation is still needed as it comes from the sesquilinear form.
120                    .dot(&sao.t())
121                    .dot(&other.coefficients),
122            }
123        };
124        Ok(ov)
125    }
126
127    /// Returns the mathematical definition of the overlap between two orbitals.
128    fn overlap_definition(&self) -> String {
129        let k = if self.complex_symmetric() { "κ " } else { "" };
130        format!("⟨{k}ψ_1|ψ_2⟩ = ∫ [{k}ψ_1(x)]* ψ_2(x) dx")
131    }
132}
133
134// =============================
135// MolecularOrbitalSymmetryOrbit
136// =============================
137
138// -----------------
139// Struct definition
140// -----------------
141
142/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of molecular
143/// orbitals.
144#[derive(Builder, Clone)]
145pub struct MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
146where
147    G: SymmetryGroupProperties,
148    T: ComplexFloat + fmt::Debug + Lapack,
149    SC: StructureConstraint + fmt::Display,
150    MolecularOrbital<'a, T, SC>: SymmetryTransformable,
151{
152    /// The generating symmetry group.
153    group: &'a G,
154
155    /// The origin molecular orbital of the orbit.
156    origin: &'a MolecularOrbital<'a, T, SC>,
157
158    /// The threshold for determining if calculated multiplicities in representation analysis are
159    /// integral.
160    integrality_threshold: <T as ComplexFloat>::Real,
161
162    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
163    pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
164
165    /// The kind of transformation determining the way the symmetry operations in `group` act on
166    /// [`Self::origin`].
167    symmetry_transformation_kind: SymmetryTransformationKind,
168
169    /// The overlap matrix between the symmetry-equivalent molecular orbitals in the orbit.
170    #[builder(setter(skip), default = "None")]
171    smat: Option<Array2<T>>,
172
173    /// The eigenvalues of the overlap matrix between the symmetry-equivalent molecular orbitals in
174    /// the orbit.
175    #[builder(setter(skip), default = "None")]
176    pub(crate) smat_eigvals: Option<Array1<T>>,
177
178    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent molecular
179    /// orbitals in the orbit.
180    ///
181    /// See [`RepAnalysis::xmat`] for further information.
182    #[builder(setter(skip), default = "None")]
183    xmat: Option<Array2<T>>,
184
185    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
186    /// eigenvalues.
187    pub(crate) eigenvalue_comparison_mode: EigenvalueComparisonMode,
188}
189
190// ----------------------------
191// Struct method implementation
192// ----------------------------
193
194impl<'a, G, T, SC> MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
195where
196    G: SymmetryGroupProperties + Clone,
197    T: ComplexFloat + fmt::Debug + Lapack,
198    SC: StructureConstraint + Clone + fmt::Display,
199    MolecularOrbital<'a, T, SC>: SymmetryTransformable,
200{
201    /// Returns a builder to construct a new [`MolecularOrbitalSymmetryOrbit`] structure.
202    pub fn builder() -> MolecularOrbitalSymmetryOrbitBuilder<'a, G, T, SC> {
203        MolecularOrbitalSymmetryOrbitBuilder::default()
204    }
205
206    /// Constructs multiple molecular orbital orbits, each from one of the supplied orbitals.
207    ///
208    /// # Arguments
209    ///
210    /// * `group` - The orbit-generating group.
211    /// * `orbitals` - The origin orbitals, each of which generates its own orbit.
212    /// * `sym_kind` - The symmetry transformation kind.
213    /// * `integrality_thresh` - The threshold of integrality check of multiplicity coefficients in
214    /// each orbit.
215    /// * `linear_independence_thresh` - The threshold of linear independence for each orbit.
216    ///
217    /// # Returns
218    ///
219    /// A vector of molecular orbital orbits.
220    pub fn from_orbitals(
221        group: &'a G,
222        orbitals: &'a [Vec<MolecularOrbital<'a, T, SC>>],
223        sym_kind: SymmetryTransformationKind,
224        eig_comp_mode: EigenvalueComparisonMode,
225        integrality_thresh: <T as ComplexFloat>::Real,
226        linear_independence_thresh: <T as ComplexFloat>::Real,
227    ) -> Vec<Vec<Self>> {
228        orbitals
229            .iter()
230            .map(|orbs_spin| {
231                orbs_spin
232                    .iter()
233                    .map(|orb| {
234                        MolecularOrbitalSymmetryOrbit::builder()
235                            .group(group)
236                            .origin(orb)
237                            .integrality_threshold(integrality_thresh)
238                            .linear_independence_threshold(linear_independence_thresh)
239                            .symmetry_transformation_kind(sym_kind.clone())
240                            .eigenvalue_comparison_mode(eig_comp_mode.clone())
241                            .build()
242                            .expect("Unable to construct a molecular orbital symmetry orbit.")
243                    })
244                    .collect_vec()
245            })
246            .collect_vec()
247    }
248}
249
250impl<'a, G, SC> MolecularOrbitalSymmetryOrbit<'a, G, f64, SC>
251where
252    G: SymmetryGroupProperties,
253    SC: StructureConstraint + fmt::Display,
254    MolecularOrbital<'a, f64, SC>: SymmetryTransformable,
255{
256    fn_calc_xmat_real!(
257        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
258        /// $`\mathbf{S}`$ between the symmetry-equivalent molecular orbitals in the orbit.
259        ///
260        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
261        ///
262        /// # Arguments
263        ///
264        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
265        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit molecular
266        /// orbitals. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
267        /// already of full rank.
268        pub calc_xmat
269    );
270}
271
272impl<'a, G, T, SC> MolecularOrbitalSymmetryOrbit<'a, G, Complex<T>, SC>
273where
274    G: SymmetryGroupProperties,
275    T: Float + Scalar<Complex = Complex<T>>,
276    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
277    SC: StructureConstraint + fmt::Display,
278    MolecularOrbital<'a, Complex<T>, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
279{
280    fn_calc_xmat_complex!(
281        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
282        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent molecular orbitals in the orbit.
283        ///
284        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
285        ///
286        /// # Arguments
287        ///
288        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
289        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit molecular
290        /// orbitals. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
291        /// already of full rank.
292        pub calc_xmat
293    );
294}
295
296// ---------------------
297// Trait implementations
298// ---------------------
299
300// ~~~~~
301// Orbit
302// ~~~~~
303
304impl<'a, G, T, SC> Orbit<G, MolecularOrbital<'a, T, SC>>
305    for MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
306where
307    G: SymmetryGroupProperties,
308    T: ComplexFloat + fmt::Debug + Lapack,
309    SC: StructureConstraint + fmt::Display,
310    MolecularOrbital<'a, T, SC>: SymmetryTransformable,
311{
312    type OrbitIter = OrbitIterator<'a, G, MolecularOrbital<'a, T, SC>>;
313
314    fn group(&self) -> &G {
315        self.group
316    }
317
318    fn origin(&self) -> &MolecularOrbital<'a, T, SC> {
319        self.origin
320    }
321
322    fn iter(&self) -> Self::OrbitIter {
323        OrbitIterator::new(
324            self.group,
325            self.origin,
326            match self.symmetry_transformation_kind {
327                SymmetryTransformationKind::Spatial => |op, orb| {
328                    orb.sym_transform_spatial(op).with_context(|| {
329                        format!("Unable to apply `{op}` spatially on the origin orbital")
330                    })
331                },
332                SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, orb| {
333                    orb.sym_transform_spatial_with_spintimerev(op).with_context(|| {
334                        format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin orbital")
335                    })
336                },
337                SymmetryTransformationKind::Spin => |op, orb| {
338                    orb.sym_transform_spin(op).with_context(|| {
339                        format!("Unable to apply `{op}` spin-wise on the origin orbital")
340                    })
341                },
342                SymmetryTransformationKind::SpinSpatial => |op, orb| {
343                    orb.sym_transform_spin_spatial(op).with_context(|| {
344                        format!("Unable to apply `{op}` spin-spatially on the origin orbital",)
345                    })
346                },
347            },
348        )
349    }
350}
351
352// ~~~~~~~~~~~
353// RepAnalysis
354// ~~~~~~~~~~~
355
356impl<'a, G, T, SC> RepAnalysis<G, MolecularOrbital<'a, T, SC>, T, Ix2>
357    for MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
358where
359    G: SymmetryGroupProperties,
360    G::CharTab: SubspaceDecomposable<T>,
361    T: Lapack
362        + ComplexFloat<Real = <T as Scalar>::Real>
363        + fmt::Debug
364        + Mul<<T as ComplexFloat>::Real, Output = T>,
365    <T as ComplexFloat>::Real: fmt::Debug
366        + Zero
367        + approx::RelativeEq<<T as ComplexFloat>::Real>
368        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
369    SC: StructureConstraint + Eq + fmt::Display,
370    MolecularOrbital<'a, T, SC>: SymmetryTransformable,
371{
372    fn set_smat(&mut self, smat: Array2<T>) {
373        self.smat = Some(smat)
374    }
375
376    fn smat(&self) -> Option<&Array2<T>> {
377        self.smat.as_ref()
378    }
379
380    fn xmat(&self) -> &Array2<T> {
381        self.xmat
382            .as_ref()
383            .expect("Orbit overlap orthogonalisation matrix not found.")
384    }
385
386    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
387        if self.origin.complex_symmetric {
388            Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
389        } else {
390            if self
391                .group
392                .get_index(i)
393                .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
394                .contains_time_reversal()
395            {
396                Ok(ComplexFloat::conj)
397            } else {
398                Ok(|x| x)
399            }
400        }
401    }
402
403    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
404        self.integrality_threshold
405    }
406
407    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
408        &self.eigenvalue_comparison_mode
409    }
410
411    /// Reduces the representation or corepresentation spanned by the molecular orbitals in the
412    /// orbit to a direct sum of the irreducible representations or corepresentations of the
413    /// generating symmetry group.
414    ///
415    /// # Returns
416    ///
417    /// The decomposed result.
418    ///
419    /// # Errors
420    ///
421    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
422    /// are non-integral, or also because the combination of group type and transformation type
423    /// would not give sensible symmetry results for a single-electron orbital. In particular, spin
424    /// or spin-spatial symmetry analysis in unitary-represented magnetic groups is not valid for
425    /// one-electron orbitals.
426    fn analyse_rep(
427        &self,
428    ) -> Result<
429        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
430        DecompositionError,
431    > {
432        // A single electron; validity depends on group and orbit type
433        let (valid_symmetry, err_str) = match self.symmetry_transformation_kind {
434                SymmetryTransformationKind::Spatial => (true, String::new()),
435                SymmetryTransformationKind::SpatialWithSpinTimeReversal
436                    | SymmetryTransformationKind::Spin
437                    | SymmetryTransformationKind::SpinSpatial => {
438                    match self.group().group_type() {
439                        GroupType::Ordinary(_) => (true, String::new()),
440                        GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
441                            (!self.group().unitary_represented(),
442                            "Unitary-represented magnetic groups cannot be used for symmetry analysis of a one-electron molecular orbital where spin is treated explicitly.".to_string())
443                        }
444                    }
445                }
446        };
447        if valid_symmetry {
448            log::debug!("Analysing representation symmetry for an MO...");
449            let chis = self
450                .calc_characters()
451                .map_err(|err| DecompositionError(err.to_string()))?;
452            let res = self.group().character_table().reduce_characters(
453                &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
454                self.integrality_threshold(),
455            );
456            log::debug!("Analysing representation symmetry for an MO... Done.");
457            res
458        } else {
459            Err(DecompositionError(err_str))
460        }
461    }
462}
463
464// ---------
465// Functions
466// ---------
467
468/// Given an origin determinant, generates the determinant orbit and all molecular-orbital orbits
469/// in tandem while populating their $`\mathbf{S}`$ matrices at the same time.
470///
471/// The evaluation of the $`\mathbf{S}`$ matrix for the determinant orbit passes through
472/// intermediate values that can be used to populate the $`\mathbf{S}`$ matrices of the
473/// molecular-orbital orbits, so it saves time significantly to construct all orbits together.
474///
475/// # Arguments
476///
477/// * `det` - An origin determinant.
478/// * `mos` - The molecular orbitals of `det`.
479/// * `group` - An orbit-generating symmetry group.
480/// * `metric` - The metric of the basis in which the coefficients of `det` and `mos` are written.
481/// * `integrality_threshold` - The threshold of integrality check of multiplicity coefficients in
482/// each orbit.
483/// * `linear_independence_threshold` - The threshold of linear independence for each orbit.
484/// * `symmetry_transformation_kind` - The kind of symmetry transformation to be applied to the
485/// * `eigenvalue_comparison_mode` - The mode of comparing the overlap eigenvalues to the specified
486/// `linear_independence_threshold`.
487/// * `use_cayley_table` - A boolean indicating if the Cayley table of the group, if available,
488/// should be used to speed up the computation of the orbit overlap matrix.
489///
490/// # Returns
491///
492/// A tuple consisting of:
493/// - the determinant orbit, and
494/// - a vector of vectors of molecular-orbital orbits, where each element of the outer vector is
495/// for one spin space, and each element of an inner vector is for one molecular orbital.
496pub fn generate_det_mo_orbits<'a, G, T, SC>(
497    det: &'a SlaterDeterminant<'a, T, SC>,
498    mos: &'a [Vec<MolecularOrbital<'a, T, SC>>],
499    group: &'a G,
500    metric: &Array2<T>,
501    metric_h: Option<&Array2<T>>,
502    integrality_threshold: <T as ComplexFloat>::Real,
503    linear_independence_threshold: <T as ComplexFloat>::Real,
504    symmetry_transformation_kind: SymmetryTransformationKind,
505    eigenvalue_comparison_mode: EigenvalueComparisonMode,
506    use_cayley_table: bool,
507) -> Result<
508    (
509        SlaterDeterminantSymmetryOrbit<'a, G, T, SC>,
510        Vec<Vec<MolecularOrbitalSymmetryOrbit<'a, G, T, SC>>>,
511    ),
512    anyhow::Error,
513>
514where
515    G: SymmetryGroupProperties + Clone,
516    G::CharTab: SubspaceDecomposable<T>,
517    T: Lapack
518        + ComplexFloat<Real = <T as Scalar>::Real>
519        + fmt::Debug
520        + Mul<<T as ComplexFloat>::Real, Output = T>,
521    <T as ComplexFloat>::Real: fmt::Debug
522        + Zero
523        + From<u16>
524        + ToPrimitive
525        + approx::RelativeEq<<T as ComplexFloat>::Real>
526        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
527    SC: StructureConstraint + Clone + Eq + fmt::Display,
528    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
529    MolecularOrbital<'a, T, SC>: SymmetryTransformable,
530{
531    log::debug!("Constructing determinant and MO orbits in tandem...");
532    let order = group.order();
533    let mut det_orbit = SlaterDeterminantSymmetryOrbit::<G, T, SC>::builder()
534        .group(group)
535        .origin(det)
536        .integrality_threshold(integrality_threshold)
537        .linear_independence_threshold(linear_independence_threshold)
538        .symmetry_transformation_kind(symmetry_transformation_kind.clone())
539        .eigenvalue_comparison_mode(eigenvalue_comparison_mode.clone())
540        .build()
541        .map_err(|err| format_err!(err))?;
542    let mut mo_orbitss = MolecularOrbitalSymmetryOrbit::from_orbitals(
543        group,
544        mos,
545        symmetry_transformation_kind,
546        eigenvalue_comparison_mode,
547        integrality_threshold,
548        linear_independence_threshold,
549    );
550
551    let mut det_smat = Array2::<T>::zeros((order, order));
552    let mut mo_smatss = mo_orbitss
553        .iter()
554        .map(|mo_orbits| {
555            mo_orbits
556                .iter()
557                .map(|_| Array2::<T>::zeros((order, order)))
558                .collect::<Vec<_>>()
559        })
560        .collect::<Vec<_>>();
561
562    let thresh = det.threshold();
563
564    let sao = metric;
565    let sao_h = metric_h.unwrap_or(sao);
566
567    if let (Some(ctb), true) = (group.cayley_table(), use_cayley_table) {
568        log::debug!("Cayley table available. Group closure will be used to speed up overlap matrix computation.");
569        let mut det_smatw0 = Array1::<T>::zeros(order);
570        let mut mo_smatw0ss = mo_orbitss
571            .iter()
572            .map(|mo_orbits| {
573                mo_orbits
574                    .iter()
575                    .map(|_| Array1::<T>::zeros(order))
576                    .collect::<Vec<_>>()
577            })
578            .collect::<Vec<_>>();
579        let det_0 = det_orbit.origin();
580        for (w, det_w_res) in det_orbit.iter().enumerate() {
581            let det_w = det_w_res?;
582            let w0_ov = izip!(
583                det_w.coefficients(),
584                det_w.occupations(),
585                det_0.coefficients(),
586                det_0.occupations(),
587            )
588            .enumerate()
589            .map(|(ispin, (cw, occw, c0, occ0))| {
590                let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
591                let nonzero_occ_0 = occ0.iter().positions(|&occ| occ > thresh).collect_vec();
592
593                let all_mo_w0_ov_mat = if det.complex_symmetric() {
594                    match (det_w.complex_conjugated(), det_0.complex_conjugated()) {
595                        (false, false) => cw.t().dot(sao_h).dot(c0),
596                        (true, false) => cw.t().dot(sao).dot(c0),
597                        (false, true) => c0.t().dot(sao).dot(cw),
598                        (true, true) => cw.t().dot(&sao_h.t()).dot(c0),
599                    }
600                } else {
601                    match (det_w.complex_conjugated(), det_0.complex_conjugated()) {
602                        (false, false) => cw.t().mapv(|x| x.conj()).dot(sao).dot(c0),
603                        (true, false) => cw.t().mapv(|x| x.conj()).dot(sao_h).dot(c0),
604                        (false, true) => c0
605                            .t()
606                            .mapv(|x| x.conj())
607                            .dot(sao_h)
608                            .dot(cw)
609                            .mapv(|x| x.conj()),
610                        (true, true) => cw.t().mapv(|x| x.conj()).dot(&sao.t()).dot(c0),
611                    }
612                };
613
614                all_mo_w0_ov_mat
615                    .diag()
616                    .iter()
617                    .enumerate()
618                    .for_each(|(imo, mo_w0_ov)| {
619                        mo_smatw0ss[ispin][imo][w] = *mo_w0_ov;
620                    });
621
622                let occ_mo_w0_ov_mat = all_mo_w0_ov_mat
623                    .select(Axis(0), &nonzero_occ_w)
624                    .select(Axis(1), &nonzero_occ_0);
625
626                occ_mo_w0_ov_mat
627                    .det()
628                    .expect("The determinant of the MO overlap matrix `w0` could not be found.")
629            })
630            .fold(T::one(), |acc, x| acc * x);
631
632            let implicit_factor = det.structure_constraint().implicit_factor()?;
633            let w0_ov = if implicit_factor > 1 {
634                let p_i32 = i32::try_from(implicit_factor)?;
635                ComplexFloat::powi(w0_ov, p_i32)
636            } else {
637                w0_ov
638            };
639            det_smatw0[w] = w0_ov;
640        }
641
642        for (i, j) in (0..order).cartesian_product(0..order) {
643            let jinv = ctb
644                .slice(s![.., j])
645                .iter()
646                .position(|&x| x == 0)
647                .ok_or(format_err!(
648                    "Unable to find the inverse of group element `{j}`."
649                ))?;
650            let jinv_i = ctb[(jinv, i)];
651
652            for (ispin, mo_smatw0s) in mo_smatw0ss.iter().enumerate() {
653                for (imo, mo_smat_w0) in mo_smatw0s.iter().enumerate() {
654                    mo_smatss[ispin][imo][(i, j)] =
655                        mo_orbitss[ispin][imo].norm_preserving_scalar_map(jinv)?(mo_smat_w0[jinv_i])
656                }
657            }
658
659            det_smat[(i, j)] = det_orbit.norm_preserving_scalar_map(jinv)?(det_smatw0[jinv_i]);
660        }
661    } else {
662        log::debug!("Cayley table not available or the use of Cayley table not requested. Overlap matrix will be constructed without group-closure speed-up.");
663        let indexed_dets = det_orbit
664            .iter()
665            .map(|det_res| det_res.map_err(|err| err.to_string()))
666            .enumerate()
667            .collect::<Vec<_>>();
668        for det_pair in indexed_dets.iter().product_repeat(2) {
669            let (w, det_w_res) = &det_pair[0];
670            let (x, det_x_res) = &det_pair[1];
671            let det_w = det_w_res
672                .as_ref()
673                .map_err(|err| format_err!(err.to_owned()))
674                .with_context(|| "One of the determinants in the orbit is not available")?;
675            let det_x = det_x_res
676                .as_ref()
677                .map_err(|err| format_err!(err.to_owned()))
678                .with_context(|| "One of the determinants in the orbit is not available")?;
679
680            let wx_ov = izip!(
681                det_w.coefficients(),
682                det_w.occupations(),
683                det_x.coefficients(),
684                det_x.occupations(),
685            )
686            .enumerate()
687            .map(|(ispin, (cw, occw, cx, occx))| {
688                let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
689                let nonzero_occ_x = occx.iter().positions(|&occ| occ > thresh).collect_vec();
690
691                // let all_mo_wx_ov_mat = if det.complex_symmetric() {
692                //     cw.t().dot(metric).dot(cx)
693                // } else {
694                //     cw.t().mapv(|x| x.conj()).dot(metric).dot(cx)
695                // };
696                let all_mo_wx_ov_mat = if det.complex_symmetric() {
697                    match (det_w.complex_conjugated(), det_x.complex_conjugated()) {
698                        (false, false) => cw.t().dot(sao_h).dot(cx),
699                        (true, false) => cw.t().dot(sao).dot(cx),
700                        (false, true) => cx.t().dot(sao).dot(cw),
701                        (true, true) => cw.t().dot(&sao_h.t()).dot(cx),
702                    }
703                } else {
704                    match (det_w.complex_conjugated(), det_x.complex_conjugated()) {
705                        (false, false) => cw.t().mapv(|x| x.conj()).dot(sao).dot(cx),
706                        (true, false) => cw.t().mapv(|x| x.conj()).dot(sao_h).dot(cx),
707                        (false, true) => cx
708                            .t()
709                            .mapv(|x| x.conj())
710                            .dot(sao_h)
711                            .dot(cw)
712                            .mapv(|x| x.conj()),
713                        (true, true) => cw.t().mapv(|x| x.conj()).dot(&sao.t()).dot(cx),
714                    }
715                };
716
717                all_mo_wx_ov_mat
718                    .diag()
719                    .iter()
720                    .enumerate()
721                    .for_each(|(imo, mo_wx_ov)| {
722                        mo_smatss[ispin][imo][(*w, *x)] = *mo_wx_ov;
723                    });
724
725                let occ_mo_wx_ov_mat = all_mo_wx_ov_mat
726                    .select(Axis(0), &nonzero_occ_w)
727                    .select(Axis(1), &nonzero_occ_x);
728
729                occ_mo_wx_ov_mat
730                    .det()
731                    .expect("The determinant of the MO overlap matrix could not be found.")
732            })
733            .fold(T::one(), |acc, x| acc * x);
734
735            let implicit_factor = det.structure_constraint().implicit_factor()?;
736            let wx_ov = if implicit_factor > 1 {
737                let p_i32 = i32::try_from(implicit_factor)?;
738                ComplexFloat::powi(wx_ov, p_i32)
739            } else {
740                wx_ov
741            };
742            det_smat[(*w, *x)] = wx_ov;
743        }
744    }
745
746    if det_orbit.origin().complex_symmetric() {
747        det_orbit.set_smat(
748            (det_smat.clone() + det_smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())),
749        );
750    } else {
751        det_orbit.set_smat(
752            (det_smat.clone() + det_smat.t().to_owned().mapv(|x| x.conj()))
753                .mapv(|x| x / (T::one() + T::one())),
754        );
755    };
756
757    mo_orbitss
758        .iter_mut()
759        .enumerate()
760        .for_each(|(ispin, mo_orbits)| {
761            mo_orbits
762                .iter_mut()
763                .enumerate()
764                .for_each(|(imo, mo_orbit)| {
765                    if mo_orbit.origin().complex_symmetric() {
766                        mo_orbit.set_smat(
767                            (mo_smatss[ispin][imo].clone() + mo_smatss[ispin][imo].t().to_owned())
768                                .mapv(|x| x / (T::one() + T::one())),
769                        )
770                    } else {
771                        mo_orbit.set_smat(
772                            (mo_smatss[ispin][imo].clone()
773                                + mo_smatss[ispin][imo].t().to_owned().mapv(|x| x.conj()))
774                            .mapv(|x| x / (T::one() + T::one())),
775                        )
776                    }
777                })
778        });
779
780    log::debug!("Constructing determinant and MO orbits in tandem... Done.");
781    Ok((det_orbit, mo_orbitss))
782}