qsym2/target/determinant/
determinant_analysis.rs

1//! Implementation of symmetry analysis for Slater determinants.
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::{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;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::group::GroupType;
31use crate::io::format::{log_subtitle, qsym2_output, QSym2Output};
32use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
33use crate::symmetry::symmetry_group::SymmetryGroupProperties;
34use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
35use crate::target::determinant::SlaterDeterminant;
36
37// =======
38// Overlap
39// =======
40
41impl<'a, T, SC> Overlap<T, Ix2> for SlaterDeterminant<'a, T, SC>
42where
43    T: Lapack
44        + ComplexFloat<Real = <T as Scalar>::Real>
45        + fmt::Debug
46        + Mul<<T as ComplexFloat>::Real, Output = T>,
47    <T as ComplexFloat>::Real: fmt::Debug
48        + approx::RelativeEq<<T as ComplexFloat>::Real>
49        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
50    SC: StructureConstraint + Eq + fmt::Display,
51{
52    fn complex_symmetric(&self) -> bool {
53        self.complex_symmetric
54    }
55
56    /// Computes the overlap between two Slater determinants.
57    ///
58    /// Determinants with fractional electrons are currently not supported.
59    ///
60    /// When one or both of the Slater determinants have been acted on by an antiunitary operation,
61    /// the correct Hermitian or complex-symmetric metric will be chosen in the evalulation of the
62    /// overlap.
63    ///
64    /// # Arguments
65    ///
66    /// * `metric` - The atomic-orbital overlap matrix with respect to the conventional sesquilinear
67    /// inner product.
68    /// * `metric_h` - The atomic-orbital overlap matrix with respect to the bilinear inner product.
69    ///
70    /// # Panics
71    ///
72    /// Panics if `self` and `other` have mismatched spin constraints or numbers of coefficient
73    /// matrices, or if fractional occupation numbers are detected.
74    fn overlap(
75        &self,
76        other: &Self,
77        metric: Option<&Array2<T>>,
78        metric_h: Option<&Array2<T>>,
79    ) -> Result<T, anyhow::Error> {
80        ensure!(
81            self.structure_constraint == other.structure_constraint,
82            "Inconsistent structure constraints between `self` and `other`."
83        );
84        ensure!(
85            self.coefficients.len() == other.coefficients.len(),
86            "Inconsistent numbers of coefficient matrices between `self` and `other`."
87        );
88        ensure!(
89            self.bao == other.bao,
90            "Inconsistent basis angular order between `self` and `other`."
91        );
92
93        let thresh = Float::sqrt(self.threshold * other.threshold);
94        ensure!(self
95            .occupations
96            .iter()
97            .chain(other.occupations.iter())
98            .all(|occs| occs.iter().all(|&occ| approx::relative_eq!(
99                occ,
100                occ.round(),
101                epsilon = thresh,
102                max_relative = thresh
103            ))),
104            "Overlaps between determinants with fractional occupation numbers are currently not supported."
105        );
106
107        let sao = metric.ok_or_else(|| format_err!("No atomic-orbital metric found."))?;
108        let sao_h = metric_h.unwrap_or(sao);
109
110        let ov = izip!(
111            &self.coefficients,
112            &self.occupations,
113            &other.coefficients,
114            &other.occupations
115        )
116        .map(|(cw, occw, cx, occx)| {
117            let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
118            let cw_o = cw.select(Axis(1), &nonzero_occ_w);
119            let nonzero_occ_x = occx.iter().positions(|&occ| occ > thresh).collect_vec();
120            let cx_o = cx.select(Axis(1), &nonzero_occ_x);
121
122            let mo_ov_mat = if self.complex_symmetric() {
123                match (self.complex_conjugated, other.complex_conjugated) {
124                    (false, false) => cw_o.t().dot(sao_h).dot(&cx_o),
125                    (true, false) => cw_o.t().dot(sao).dot(&cx_o),
126                    (false, true) => cx_o.t().dot(sao).dot(&cw_o),
127                    (true, true) => cw_o.t().dot(&sao_h.t()).dot(&cx_o),
128                }
129            } else {
130                match (self.complex_conjugated, other.complex_conjugated) {
131                    (false, false) => cw_o.t().mapv(|x| x.conj()).dot(sao).dot(&cx_o),
132                    (true, false) => cw_o.t().mapv(|x| x.conj()).dot(sao_h).dot(&cx_o),
133                    (false, true) => cx_o
134                        .t()
135                        .mapv(|x| x.conj())
136                        .dot(sao_h)
137                        .dot(&cw_o)
138                        .mapv(|x| x.conj()),
139                    (true, true) => cw_o.t().mapv(|x| x.conj()).dot(&sao.t()).dot(&cx_o),
140                }
141            };
142            mo_ov_mat
143                .det()
144                .expect("The determinant of the MO overlap matrix could not be found.")
145        })
146        .fold(T::one(), |acc, x| acc * x);
147
148        let implicit_factor = self.structure_constraint.implicit_factor()?;
149        if implicit_factor > 1 {
150            let p_i32 = i32::try_from(implicit_factor)?;
151            Ok(ComplexFloat::powi(ov, p_i32))
152        } else {
153            Ok(ov)
154        }
155    }
156
157    /// Returns the mathematical definition of the overlap between two Slater determinants.
158    fn overlap_definition(&self) -> String {
159        let k = if self.complex_symmetric() { "κ " } else { "" };
160        format!("⟨{k}Ψ_1|Ψ_2⟩ = ∫ [{k}Ψ_1(x^Ne)]* Ψ_2(x^Ne) dx^Ne")
161    }
162}
163
164// ==============================
165// SlaterDeterminantSymmetryOrbit
166// ==============================
167
168// -----------------
169// Struct definition
170// -----------------
171
172/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of Slater
173/// determinants.
174#[derive(Builder, Clone)]
175pub struct SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
176where
177    G: SymmetryGroupProperties,
178    T: ComplexFloat + fmt::Debug + Lapack,
179    SC: StructureConstraint + fmt::Display,
180    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
181{
182    /// The generating symmetry group.
183    group: &'a G,
184
185    /// The origin Slater determinant of the orbit.
186    origin: &'a SlaterDeterminant<'a, T, SC>,
187
188    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
189    linear_independence_threshold: <T as ComplexFloat>::Real,
190
191    /// The threshold for determining if calculated multiplicities in representation analysis are
192    /// integral.
193    integrality_threshold: <T as ComplexFloat>::Real,
194
195    /// The kind of transformation determining the way the symmetry operations in `group` act on
196    /// [`Self::origin`].
197    symmetry_transformation_kind: SymmetryTransformationKind,
198
199    /// The overlap matrix between the symmetry-equivalent Slater determinants in the orbit.
200    #[builder(setter(skip), default = "None")]
201    smat: Option<Array2<T>>,
202
203    /// The eigenvalues of the overlap matrix between the symmetry-equivalent Slater determinants in
204    /// the orbit.
205    #[builder(setter(skip), default = "None")]
206    pub(crate) smat_eigvals: Option<Array1<T>>,
207
208    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent Slater
209    /// determinants in the orbit.
210    ///
211    /// See [`RepAnalysis::xmat`] for further information.
212    #[builder(setter(skip), default = "None")]
213    xmat: Option<Array2<T>>,
214
215    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
216    /// eigenvalues.
217    eigenvalue_comparison_mode: EigenvalueComparisonMode,
218}
219
220// ----------------------------
221// Struct method implementation
222// ----------------------------
223
224impl<'a, G, T, SC> SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
225where
226    G: SymmetryGroupProperties + Clone,
227    T: ComplexFloat + fmt::Debug + Lapack,
228    SC: StructureConstraint + Clone + fmt::Display,
229    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
230{
231    /// Returns a builder for constructing a new Slater determinant symmetry orbit.
232    pub fn builder() -> SlaterDeterminantSymmetryOrbitBuilder<'a, G, T, SC> {
233        SlaterDeterminantSymmetryOrbitBuilder::default()
234    }
235}
236
237impl<'a, G, SC> SlaterDeterminantSymmetryOrbit<'a, G, f64, SC>
238where
239    G: SymmetryGroupProperties,
240    SC: StructureConstraint + fmt::Display,
241    SlaterDeterminant<'a, f64, SC>: SymmetryTransformable,
242{
243    fn_calc_xmat_real!(
244        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
245        /// $`\mathbf{S}`$ between the symmetry-equivalent Slater determinants in the orbit.
246        ///
247        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
248        ///
249        /// # Arguments
250        ///
251        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
252        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit determinants.
253        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
254        /// full rank.
255        pub calc_xmat
256    );
257}
258
259impl<'a, G, T, SC> SlaterDeterminantSymmetryOrbit<'a, G, Complex<T>, SC>
260where
261    G: SymmetryGroupProperties,
262    T: Float + Scalar<Complex = Complex<T>>,
263    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
264    SC: StructureConstraint + fmt::Display,
265    SlaterDeterminant<'a, Complex<T>, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
266{
267    fn_calc_xmat_complex!(
268        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
269        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent Slater determinants in the orbit.
270        ///
271        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
272        ///
273        /// # Arguments
274        ///
275        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
276        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit determinants.
277        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
278        /// full rank.
279        pub calc_xmat
280    );
281}
282
283// ---------------------
284// Trait implementations
285// ---------------------
286
287// ~~~~~
288// Orbit
289// ~~~~~
290
291impl<'a, G, T, SC> Orbit<G, SlaterDeterminant<'a, T, SC>>
292    for SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
293where
294    G: SymmetryGroupProperties,
295    T: ComplexFloat + fmt::Debug + Lapack,
296    SC: StructureConstraint + fmt::Display,
297    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
298{
299    type OrbitIter = OrbitIterator<'a, G, SlaterDeterminant<'a, T, SC>>;
300
301    fn group(&self) -> &G {
302        self.group
303    }
304
305    fn origin(&self) -> &SlaterDeterminant<'a, T, SC> {
306        self.origin
307    }
308
309    fn iter(&self) -> Self::OrbitIter {
310        OrbitIterator::new(
311            self.group,
312            self.origin,
313            match self.symmetry_transformation_kind {
314                SymmetryTransformationKind::Spatial => |op, det| {
315                    det.sym_transform_spatial(op).with_context(|| {
316                        format!("Unable to apply `{op}` spatially on the origin determinant")
317                    })
318                },
319                SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, det| {
320                    det.sym_transform_spatial_with_spintimerev(op).with_context(|| {
321                        format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin determinant")
322                    })
323                },
324                SymmetryTransformationKind::Spin => |op, det| {
325                    det.sym_transform_spin(op).with_context(|| {
326                        format!("Unable to apply `{op}` spin-wise on the origin determinant")
327                    })
328                },
329                SymmetryTransformationKind::SpinSpatial => |op, det| {
330                    det.sym_transform_spin_spatial(op).with_context(|| {
331                        format!("Unable to apply `{op}` spin-spatially on the origin determinant")
332                    })
333                },
334            },
335        )
336    }
337}
338
339// ~~~~~~~~~~~
340// RepAnalysis
341// ~~~~~~~~~~~
342
343impl<'a, G, T, SC> RepAnalysis<G, SlaterDeterminant<'a, T, SC>, T, Ix2>
344    for SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
345where
346    G: SymmetryGroupProperties,
347    G::CharTab: SubspaceDecomposable<T>,
348    T: Lapack
349        + ComplexFloat<Real = <T as Scalar>::Real>
350        + fmt::Debug
351        + Mul<<T as ComplexFloat>::Real, Output = T>,
352    <T as ComplexFloat>::Real: fmt::Debug
353        + Zero
354        + From<u16>
355        + ToPrimitive
356        + approx::RelativeEq<<T as ComplexFloat>::Real>
357        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
358    SC: StructureConstraint + Eq + fmt::Display,
359    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
360{
361    fn set_smat(&mut self, smat: Array2<T>) {
362        self.smat = Some(smat)
363    }
364
365    fn smat(&self) -> Option<&Array2<T>> {
366        self.smat.as_ref()
367    }
368
369    fn xmat(&self) -> &Array2<T> {
370        self.xmat
371            .as_ref()
372            .expect("Orbit overlap orthogonalisation matrix not found.")
373    }
374
375    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
376        if self.origin.complex_symmetric {
377            Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
378        } else {
379            if self
380                .group
381                .get_index(i)
382                .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
383                .contains_time_reversal()
384            {
385                Ok(ComplexFloat::conj)
386            } else {
387                Ok(|x| x)
388            }
389        }
390    }
391
392    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
393        self.integrality_threshold
394    }
395
396    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
397        &self.eigenvalue_comparison_mode
398    }
399
400    /// Reduces the representation or corepresentation spanned by the determinants in the orbit to
401    /// a direct sum of the irreducible representations or corepresentations of the generating
402    /// symmetry group.
403    ///
404    /// # Returns
405    ///
406    /// The decomposed result.
407    ///
408    /// # Errors
409    ///
410    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
411    /// are non-integral, or also because the combination of group type, transformation type, and
412    /// oddity of the number of electrons would not give sensible symmetry results. In particular,
413    /// spin or spin-spatial symmetry analysis of odd-electron systems in unitary-represented
414    /// magnetic groups is not valid.
415    fn analyse_rep(
416        &self,
417    ) -> Result<
418        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
419        DecompositionError,
420    > {
421        log::debug!("Analysing representation symmetry for a Slater determinant...");
422        let nelectrons_float = self.origin().nelectrons();
423        if approx::relative_eq!(
424            nelectrons_float.round(),
425            nelectrons_float,
426            epsilon = self.integrality_threshold,
427            max_relative = self.integrality_threshold
428        ) {
429            let nelectrons_usize = nelectrons_float.round().to_usize().unwrap_or_else(|| {
430                panic!(
431                    "Unable to convert the number of electrons `{nelectrons_float:.7}` to `usize`."
432                );
433            });
434            let (valid_symmetry, err_str) = if nelectrons_usize.rem_euclid(2) == 0 {
435                // Even number of electrons; always valid
436                (true, String::new())
437            } else {
438                // Odd number of electrons; validity depends on group and orbit type
439                match self.symmetry_transformation_kind {
440                    SymmetryTransformationKind::Spatial => (true, String::new()),
441                    SymmetryTransformationKind::SpatialWithSpinTimeReversal
442                        | SymmetryTransformationKind::Spin
443                        | SymmetryTransformationKind::SpinSpatial => {
444                        match self.group().group_type() {
445                            GroupType::Ordinary(_) => (true, String::new()),
446                            GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
447                                (!self.group().unitary_represented(),
448                                "Unitary-represented magnetic groups cannot be used for symmetry analysis of odd-electron systems where spin is treated explicitly.".to_string())
449                            }
450                        }
451                    }
452                }
453            };
454
455            if valid_symmetry {
456                let chis = self
457                    .calc_characters()
458                    .map_err(|err| DecompositionError(err.to_string()))?;
459                log::debug!("Characters calculated.");
460
461                log_subtitle("Determinant orbit characters");
462                qsym2_output!("");
463                self.characters_to_string(&chis, self.integrality_threshold)
464                    .log_output_display();
465                qsym2_output!("");
466
467                let res = self.group().character_table().reduce_characters(
468                    &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
469                    self.integrality_threshold(),
470                );
471                log::debug!("Characters reduced.");
472                log::debug!("Analysing representation symmetry for a Slater determinant... Done.");
473                res
474            } else {
475                Err(DecompositionError(err_str))
476            }
477        } else {
478            Err(DecompositionError(format!(
479                "Symmetry analysis for determinant with non-integer number of electrons `{nelectrons_float:.7}` (threshold = {:.3e}) not supported.",
480                self.integrality_threshold
481            )))
482        }
483    }
484}