qsym2/target/noci/
multideterminant_analysis.rs

1//! Implementation of symmetry analysis for multi-determinantal wavefunctions.
2
3use std::collections::HashSet;
4use std::fmt;
5use std::hash::Hash;
6use std::ops::Mul;
7
8use anyhow::{self, Context, ensure, format_err};
9use approx;
10use derive_builder::Builder;
11use itertools::Itertools;
12use log;
13use ndarray::{Array1, Array2, Array3, Axis, Ix2, s};
14use ndarray_linalg::{
15    UPLO,
16    eig::Eig,
17    eigh::Eigh,
18    types::{Lapack, Scalar},
19};
20use num_complex::{Complex, ComplexFloat};
21use num_traits::{Float, ToPrimitive, Zero};
22
23use crate::analysis::{
24    EigenvalueComparisonMode, Orbit, OrbitIterator, Overlap, RepAnalysis, fn_calc_xmat_complex,
25    fn_calc_xmat_real,
26};
27use crate::angmom::spinor_rotation_3d::StructureConstraint;
28use crate::auxiliary::misc::complex_modified_gram_schmidt;
29use crate::chartab::chartab_group::CharacterProperties;
30use crate::chartab::{DecompositionError, SubspaceDecomposable};
31use crate::group::GroupType;
32use crate::io::format::{QSym2Output, log_subtitle, qsym2_output};
33use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
34use crate::symmetry::symmetry_group::SymmetryGroupProperties;
35use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
36use crate::target::determinant::SlaterDeterminant;
37use crate::target::noci::backend::solver::check_complex_matrix_symmetry;
38use crate::target::noci::basis::{Basis, OrbitBasis};
39use crate::target::noci::multideterminant::MultiDeterminant;
40
41// =======
42// Overlap
43// =======
44
45impl<'a, T, B, SC> Overlap<T, Ix2> for MultiDeterminant<'a, T, B, SC>
46where
47    T: Lapack
48        + ComplexFloat<Real = <T as Scalar>::Real>
49        + fmt::Debug
50        + Mul<<T as ComplexFloat>::Real, Output = T>,
51    <T as ComplexFloat>::Real: fmt::Debug
52        + approx::RelativeEq<<T as ComplexFloat>::Real>
53        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
54    SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
55    B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
56{
57    fn complex_symmetric(&self) -> bool {
58        self.complex_symmetric
59    }
60
61    /// Computes the overlap between two multi-determinantal wavefunctions.
62    ///
63    /// Determinants with fractional electrons are currently not supported.
64    ///
65    /// When one or both of the multi-determinantal wavefunctions have been acted on by an
66    /// antiunitary operation, the correct Hermitian or complex-symmetric metric will be chosen in
67    /// the evalulation of the overlap.
68    ///
69    /// # Arguments
70    ///
71    /// * `metric` - The atomic-orbital overlap matrix with respect to the conventional sesquilinear
72    ///   inner product.
73    /// * `metric_h` - The atomic-orbital overlap matrix with respect to the bilinear inner product.
74    ///
75    /// # Panics
76    ///
77    /// Panics if `self` and `other` have mismatched spin constraints or numbers of coefficient
78    /// matrices, or if fractional occupation numbers are detected.
79    fn overlap(
80        &self,
81        other: &Self,
82        metric: Option<&Array2<T>>,
83        metric_h: Option<&Array2<T>>,
84    ) -> Result<T, anyhow::Error> {
85        ensure!(
86            self.structure_constraint() == other.structure_constraint(),
87            "Inconsistent structure constraints between `self` and `other`."
88        );
89        ensure!(
90            self.coefficients.len() == other.coefficients.len(),
91            "Inconsistent numbers of coefficients between `self` and `other`."
92        );
93
94        let s_dets = self.basis.iter().collect::<Result<Vec<_>, _>>()?;
95        let o_dets = other.basis.iter().collect::<Result<Vec<_>, _>>()?;
96        let swx_vec = s_dets
97            .iter()
98            .cartesian_product(o_dets.iter())
99            .map(|(w, x)| w.overlap(x, metric, metric_h))
100            .collect::<Result<Vec<_>, _>>()?;
101        let d = self.coefficients.len();
102        let swx = Array2::from_shape_vec((d, d), swx_vec)?;
103
104        if self.complex_symmetric {
105            Ok(self.coefficients.t().dot(&swx).dot(&other.coefficients))
106        } else {
107            Ok(self
108                .coefficients
109                .t()
110                .mapv(|x| x.conj())
111                .dot(&swx)
112                .dot(&other.coefficients))
113        }
114    }
115
116    /// Returns the mathematical definition of the overlap between two multi-determinantal
117    /// wavefunctions.
118    fn overlap_definition(&self) -> String {
119        let k = if self.complex_symmetric() { "κ " } else { "" };
120        format!("⟨{k}Ψ_1|Ψ_2⟩ = ∫ [{k}Ψ_1(x^Ne)]* Ψ_2(x^Ne) dx^Ne")
121    }
122}
123
124// =============================
125// MultiDeterminantSymmetryOrbit
126// =============================
127
128// -----------------
129// Struct definition
130// -----------------
131
132/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of
133/// multi-determinantal wavefunctions.
134#[derive(Builder, Clone)]
135pub struct MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
136where
137    G: SymmetryGroupProperties,
138    T: ComplexFloat + fmt::Debug + Lapack,
139    SC: StructureConstraint + Hash + Eq + fmt::Display,
140    B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
141    MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
142{
143    /// The generating symmetry group.
144    group: &'g G,
145
146    /// The origin multi-determinantal wavefunction of the orbit.
147    origin: &'a MultiDeterminant<'a, T, B, SC>,
148
149    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
150    pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
151
152    /// The threshold for determining if calculated multiplicities in representation analysis are
153    /// integral.
154    integrality_threshold: <T as ComplexFloat>::Real,
155
156    /// The kind of transformation determining the way the symmetry operations in `group` act on
157    /// [`Self::origin`].
158    symmetry_transformation_kind: SymmetryTransformationKind,
159
160    /// The overlap matrix between the symmetry-equivalent multi-eterminantal wavefunctions in the
161    /// orbit.
162    #[builder(setter(skip), default = "None")]
163    smat: Option<Array2<T>>,
164
165    /// The eigenvalues of the overlap matrix between the symmetry-equivalent multi-determinantal
166    /// wavefunctions in the orbit.
167    #[builder(setter(skip), default = "None")]
168    pub(crate) smat_eigvals: Option<Array1<T>>,
169
170    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent
171    /// multi-determinantal wavefunctions in the orbit.
172    ///
173    /// See [`RepAnalysis::xmat`] for further information.
174    #[builder(setter(skip), default = "None")]
175    xmat: Option<Array2<T>>,
176
177    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
178    /// eigenvalues.
179    eigenvalue_comparison_mode: EigenvalueComparisonMode,
180}
181
182// ----------------------------
183// Struct method implementation
184// ----------------------------
185
186impl<'a, 'g, G, T, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
187where
188    G: SymmetryGroupProperties + Clone,
189    T: ComplexFloat + fmt::Debug + Lapack,
190    SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
191    B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
192    MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
193{
194    /// Returns a builder for constructing a new multi-determinantal wavefunction symmetry orbit.
195    pub fn builder() -> MultiDeterminantSymmetryOrbitBuilder<'a, 'g, G, T, B, SC> {
196        MultiDeterminantSymmetryOrbitBuilder::default()
197    }
198
199    /// Returns the origin of the multi-determinantal wavefunction symmetry orbit.
200    pub fn origin(&self) -> &MultiDeterminant<'a, T, B, SC> {
201        self.origin
202    }
203}
204
205impl<'a, 'g, G, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, f64, B, SC>
206where
207    G: SymmetryGroupProperties,
208    SC: StructureConstraint + Hash + Eq + fmt::Display,
209    B: 'a + Basis<SlaterDeterminant<'a, f64, SC>> + Clone,
210    MultiDeterminant<'a, f64, B, SC>: SymmetryTransformable,
211{
212    fn_calc_xmat_real!(
213        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
214        /// $`\mathbf{S}`$ between the symmetry-equivalent Slater determinants in the orbit.
215        ///
216        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
217        ///
218        /// # Arguments
219        ///
220        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
221        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit determinants.
222        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
223        /// full rank.
224        pub calc_xmat
225    );
226}
227
228impl<'a, 'g, G, T, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, Complex<T>, B, SC>
229where
230    G: SymmetryGroupProperties,
231    T: Float + Scalar<Complex = Complex<T>>,
232    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
233    SC: StructureConstraint + Hash + Eq + fmt::Display,
234    B: 'a + Basis<SlaterDeterminant<'a, Complex<T>, SC>> + Clone,
235    MultiDeterminant<'a, Complex<T>, B, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
236{
237    fn_calc_xmat_complex!(
238        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
239        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent Slater determinants in the orbit.
240        ///
241        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
242        ///
243        /// # Arguments
244        ///
245        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
246        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit determinants.
247        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
248        /// full rank.
249        pub calc_xmat
250    );
251}
252
253// ---------------------
254// Trait implementations
255// ---------------------
256
257// ~~~~~
258// Orbit
259// ~~~~~
260
261impl<'a, 'g, G, T, B, SC> Orbit<G, MultiDeterminant<'a, T, B, SC>>
262    for MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
263where
264    G: SymmetryGroupProperties,
265    T: ComplexFloat + fmt::Debug + Lapack,
266    SC: StructureConstraint + Hash + Eq + fmt::Display,
267    B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
268    MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
269{
270    type OrbitIter = OrbitIterator<'a, G, MultiDeterminant<'a, T, B, SC>>;
271
272    fn group(&self) -> &G {
273        self.group
274    }
275
276    fn origin(&self) -> &MultiDeterminant<'a, T, B, SC> {
277        self.origin
278    }
279
280    fn iter(&self) -> Self::OrbitIter {
281        OrbitIterator::new(
282            self.group,
283            self.origin,
284            match self.symmetry_transformation_kind {
285                SymmetryTransformationKind::Spatial => |op, multidet| {
286                    multidet.sym_transform_spatial(op).with_context(|| {
287                        format!("Unable to apply `{op}` spatially on the origin multi-determinantal wavefunction")
288                    })
289                },
290                SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, multidet| {
291                    multidet.sym_transform_spatial_with_spintimerev(op).with_context(|| {
292                        format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin multi-determinantal wavefunction")
293                    })
294                },
295                SymmetryTransformationKind::Spin => |op, multidet| {
296                    multidet.sym_transform_spin(op).with_context(|| {
297                        format!("Unable to apply `{op}` spin-wise on the origin multi-determinantal wavefunction")
298                    })
299                },
300                SymmetryTransformationKind::SpinSpatial => |op, multidet| {
301                    multidet.sym_transform_spin_spatial(op).with_context(|| {
302                        format!("Unable to apply `{op}` spin-spatially on the origin multi-determinantal wavefunction")
303                    })
304                },
305            },
306        )
307    }
308}
309
310// ~~~~~~~~~~~
311// RepAnalysis
312// ~~~~~~~~~~~
313
314impl<'a, 'g, G, T, B, SC> RepAnalysis<G, MultiDeterminant<'a, T, B, SC>, T, Ix2>
315    for MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
316where
317    G: SymmetryGroupProperties,
318    G::CharTab: SubspaceDecomposable<T>,
319    T: Lapack
320        + ComplexFloat<Real = <T as Scalar>::Real>
321        + fmt::Debug
322        + Mul<<T as ComplexFloat>::Real, Output = T>,
323    <T as ComplexFloat>::Real: fmt::Debug
324        + Zero
325        + From<u16>
326        + ToPrimitive
327        + approx::RelativeEq<<T as ComplexFloat>::Real>
328        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
329    SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
330    B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
331    MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
332{
333    fn set_smat(&mut self, smat: Array2<T>) {
334        self.smat = Some(smat)
335    }
336
337    fn smat(&self) -> Option<&Array2<T>> {
338        self.smat.as_ref()
339    }
340
341    fn xmat(&self) -> &Array2<T> {
342        self.xmat
343            .as_ref()
344            .expect("Orbit overlap orthogonalisation matrix not found.")
345    }
346
347    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
348        if self.origin.complex_symmetric {
349            Err(format_err!(
350                "`norm_preserving_scalar_map` is currently not implemented for complex-symmetric overlaps. This thus precludes the use of the Cayley table to speed up the computation of the orbit overlap matrix."
351            ))
352        } else if self
353            .group
354            .get_index(i)
355            .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
356            .contains_time_reversal()
357        {
358            Ok(ComplexFloat::conj)
359        } else {
360            Ok(|x| x)
361        }
362    }
363
364    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
365        self.integrality_threshold
366    }
367
368    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
369        &self.eigenvalue_comparison_mode
370    }
371
372    /// Reduces the representation or corepresentation spanned by the multi-determinantal
373    /// wavefunctions in the orbit to a direct sum of the irreducible representations or
374    /// corepresentations of the generating symmetry group.
375    ///
376    /// # Returns
377    ///
378    /// The decomposed result.
379    ///
380    /// # Errors
381    ///
382    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
383    /// are non-integral, or also because the combination of group type, transformation type, and
384    /// oddity of the number of electrons would not give sensible symmetry results. In particular,
385    /// spin or spin-spatial symmetry analysis of odd-electron systems in unitary-represented
386    /// magnetic groups is not valid.
387    fn analyse_rep(
388        &self,
389    ) -> Result<
390        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
391        DecompositionError,
392    > {
393        log::debug!("Analysing representation symmetry for a multi-determinantal wavefunction...");
394        let mut nelectrons_set = self
395            .origin()
396            .basis
397            .iter()
398            .map(|det_res| det_res.and_then(|det| {
399                let det_nelectrons_float = det.nelectrons();
400                if approx::relative_eq!(
401                    det_nelectrons_float.round(),
402                    det_nelectrons_float,
403                    epsilon = self.integrality_threshold,
404                    max_relative = self.integrality_threshold
405                ) {
406                    det_nelectrons_float.round().to_usize().ok_or(format_err!("Unable to convert the number of electrons `{det_nelectrons_float:.7}` to `usize`."))
407                } else {
408                    Err(format_err!("Fractional number of electrons encountered: `{det_nelectrons_float:.7}`"))
409                }
410            }))
411            .collect::<Result<HashSet<_>, _>>()
412            .map_err(|err| DecompositionError(err.to_string()))?;
413        if nelectrons_set.len() != 1 {
414            Err(DecompositionError("Symmetry analysis for multi-determinantal wavefunctions with multiple numbers of electrons is not yet supported.".to_string()))
415        } else {
416            let nelectrons = nelectrons_set.drain().next().ok_or(DecompositionError(
417                "Unable to obtain the number of electrons.".to_string(),
418            ))?;
419            let (valid_symmetry, err_str) = if nelectrons.rem_euclid(2) == 0 {
420                // Even number of electrons; always valid
421                (true, String::new())
422            } else {
423                // Odd number of electrons; validity depends on group and orbit type
424                match self.symmetry_transformation_kind {
425                    SymmetryTransformationKind::Spatial => (true, String::new()),
426                    SymmetryTransformationKind::SpatialWithSpinTimeReversal
427                        | SymmetryTransformationKind::Spin
428                        | SymmetryTransformationKind::SpinSpatial => {
429                        match self.group().group_type() {
430                            GroupType::Ordinary(_) => (true, String::new()),
431                            GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
432                                (!self.group().unitary_represented(),
433                                "Unitary-represented magnetic groups cannot be used for symmetry analysis of odd-electron systems where spin is treated explicitly.".to_string())
434                            }
435                        }
436                    }
437                }
438            };
439
440            if valid_symmetry {
441                let chis = self
442                    .calc_characters()
443                    .map_err(|err| DecompositionError(err.to_string()))?;
444                log::debug!("Characters calculated.");
445
446                log_subtitle("Multi-determinantal wavefunction orbit characters");
447                qsym2_output!("");
448                self.characters_to_string(&chis, self.integrality_threshold)
449                    .log_output_display();
450                qsym2_output!("");
451
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!("Characters reduced.");
457                log::debug!(
458                    "Analysing representation symmetry for a multi-determinantal wavefunction... Done."
459                );
460                res
461            } else {
462                Err(DecompositionError(err_str))
463            }
464        }
465    }
466}
467
468// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
469// Optimised implementation for multi-determinantal wavefunctions constructed from orbits
470// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471impl<'a, 'go, 'g, G, T, SC>
472    MultiDeterminantSymmetryOrbit<
473        'a,
474        'go,
475        G,
476        T,
477        OrbitBasis<'g, G, SlaterDeterminant<'a, T, SC>>,
478        SC,
479    >
480where
481    G: SymmetryGroupProperties + Clone,
482    G::CharTab: SubspaceDecomposable<T>,
483    T: Lapack
484        + ComplexFloat<Real = <T as Scalar>::Real>
485        + fmt::Debug
486        + Mul<<T as ComplexFloat>::Real, Output = T>,
487    <T as ComplexFloat>::Real: fmt::Debug
488        + Zero
489        + From<u16>
490        + ToPrimitive
491        + approx::RelativeEq<<T as ComplexFloat>::Real>
492        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
493    SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
494    MultiDeterminant<'a, T, OrbitBasis<'g, G, SlaterDeterminant<'a, T, SC>>, SC>:
495        SymmetryTransformable,
496{
497    /// Calculates and stores the overlap matrix between multi-determinantal wavefunctions in the
498    /// orbit, with respect to a metric of the basis in which the constituting Slater determinants
499    /// are expressed.
500    ///
501    /// This function is particularly optimised for multi-determinantal wavefunctions constructed
502    /// from orbits of origin Slater determinants such that the multi-determinantal wavefunctions
503    /// are never explicitly transformed by group operations.
504    ///
505    /// # Arguments
506    ///
507    /// * `metric` - The metric of the basis in which the orbit items are expressed.
508    /// * `metric_h` - The complex-symmetric metric of the basis in which the orbit items are
509    ///   expressed. This is required if antiunitary operations are involved.
510    /// * `use_cayley_table` - A boolean indicating if the Cayley table of the group should be used
511    ///   to speed up the computation of the overlap matrix. If `false`, this will revert back to the
512    ///   non-optimised overlap matrix calculation.
513    pub(crate) fn calc_smat_optimised(
514        &mut self,
515        metric: Option<&Array2<T>>,
516        metric_h: Option<&Array2<T>>,
517        use_cayley_table: bool,
518    ) -> Result<&mut Self, anyhow::Error> {
519        ensure!(
520            self.group.name() == self.origin.basis.group().name(),
521            "Multi-determinantal wavefunction orbit-generating group does not match symmetry analysis group."
522        );
523
524        if let (Some(ctb), true) = (self.group().cayley_table(), use_cayley_table) {
525            log::debug!(
526                "Cayley table available. Group closure will be used to speed up overlap matrix computation."
527            );
528
529            let order = self.group.order();
530            let mut smat = Array2::<T>::zeros((order, order));
531            let multidet_0 = self.origin();
532            let det_origins = multidet_0.basis.origins();
533            let n_det_origins = det_origins.len();
534
535            let detov_kpwx_vec = multidet_0
536                .basis
537                .iter()
538                .collect::<Result<Vec<_>, _>>()?
539                .into_iter()
540                .cartesian_product(det_origins.iter())
541                .map(|(w, x)| w.overlap(x, metric, metric_h))
542                .collect::<Result<Vec<_>, _>>()?;
543            let detov_kpwx =
544                Array3::from_shape_vec((order, n_det_origins, n_det_origins), detov_kpwx_vec)?;
545            let ovs = (0..order)
546                .map(|k| {
547                    [
548                        (0..order),
549                        (0..order),
550                        (0..n_det_origins),
551                        (0..n_det_origins),
552                    ]
553                    .into_iter()
554                    .multi_cartesian_product()
555                    .try_fold(T::zero(), |acc, v| {
556                        let ip = v[0];
557                        let jp = v[1];
558                        let w = v[2];
559                        let x = v[3];
560                        let aipw = multidet_0.coefficients[ip * n_det_origins + w];
561                        let ajpx = multidet_0.coefficients[jp * n_det_origins + x];
562
563                        let jpinv = ctb.slice(s![.., jp]).iter().position(|&x| x == 0).ok_or(
564                            format_err!("Unable to find the inverse of group element `{jp}`."),
565                        )?;
566                        let kp = ctb[(jpinv, ctb[(k, ip)])];
567                        let ukipwjpx = self.norm_preserving_scalar_map(jp)?(detov_kpwx[(kp, w, x)]);
568
569                        Ok::<_, anyhow::Error>(
570                            acc + self.norm_preserving_scalar_map(k)?(aipw.conj())
571                                * ukipwjpx
572                                * ajpx,
573                        )
574                    })
575                })
576                .collect::<Result<Vec<_>, _>>()?;
577            for (i, j) in (0..order).cartesian_product(0..order) {
578                let jinv = ctb
579                    .slice(s![.., j])
580                    .iter()
581                    .position(|&x| x == 0)
582                    .ok_or(format_err!(
583                        "Unable to find the inverse of group element `{j}`."
584                    ))?;
585                let jinv_i = ctb[(jinv, i)];
586                smat[(i, j)] = self.norm_preserving_scalar_map(jinv)?(ovs[jinv_i]);
587            }
588            if self.origin().complex_symmetric() {
589                let _ = check_complex_matrix_symmetry(
590                    &smat.view(),
591                    true,
592                    self.linear_independence_threshold,
593                    "Orbit overlap",
594                    "S_orbit",
595                );
596                self.set_smat(
597                    (smat.clone() + smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())),
598                )
599            } else {
600                let _ = check_complex_matrix_symmetry(
601                    &smat.view(),
602                    false,
603                    self.linear_independence_threshold,
604                    "Orbit overlap",
605                    "S_orbit",
606                );
607                self.set_smat(
608                    (smat.clone() + smat.t().to_owned().mapv(|x| x.conj()))
609                        .mapv(|x| x / (T::one() + T::one())),
610                )
611            }
612            Ok(self)
613        } else {
614            self.calc_smat(metric, metric_h, use_cayley_table)
615        }
616    }
617}
qsym2/target/noci/multideterminant_analysis.rs

qsym2/target/noci/
multideterminant_analysis.rs
