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

qsym2/target/noci/
multideterminant_analysis.rs
