qsym2/target/density/
density_analysis.rs

1//! Implementation of symmetry analysis for electron densities.
2
3use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, ensure, format_err, Context};
7use approx;
8use derive_builder::Builder;
9use itertools::Itertools;
10use log;
11use ndarray::{Array1, Array2, Array4, Axis, Ix4};
12use ndarray_einsum_beta::*;
13use ndarray_linalg::{
14    eig::Eig,
15    eigh::Eigh,
16    norm::Norm,
17    types::{Lapack, Scalar},
18    UPLO,
19};
20use num_complex::{Complex, ComplexFloat};
21use num_traits::{Float, ToPrimitive, Zero};
22
23use crate::analysis::{
24    fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
25    Overlap, RepAnalysis,
26};
27use crate::auxiliary::misc::complex_modified_gram_schmidt;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
31use crate::symmetry::symmetry_group::SymmetryGroupProperties;
32use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
33use crate::target::density::Density;
34
35// =======
36// Overlap
37// =======
38
39impl<'a, T> Overlap<T, Ix4> for Density<'a, T>
40where
41    T: Lapack
42        + ComplexFloat<Real = <T as Scalar>::Real>
43        + fmt::Debug
44        + Mul<<T as ComplexFloat>::Real, Output = T>,
45    <T as ComplexFloat>::Real: fmt::Debug
46        + approx::RelativeEq<<T as ComplexFloat>::Real>
47        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
48{
49    fn complex_symmetric(&self) -> bool {
50        self.complex_symmetric
51    }
52
53    /// Computes the overlap between two densities.
54    ///
55    /// When one or both of the densities have been acted on by an antiunitary operation, the
56    /// correct Hermitian or complex-symmetric metric will be chosen in the evalulation of the
57    /// overlap.
58    ///
59    /// # Panics
60    ///
61    /// Panics if `self` and `other` have mismatched spin constraints.
62    fn overlap(
63        &self,
64        other: &Self,
65        metric: Option<&Array4<T>>,
66        metric_h: Option<&Array4<T>>,
67    ) -> Result<T, anyhow::Error> {
68        ensure!(
69            self.density_matrix.shape() == other.density_matrix.shape(),
70            "Inconsistent shapes of density matrices between `self` and `other`."
71        );
72        ensure!(
73            self.bao == other.bao,
74            "Inconsistent basis angular order between `self` and `other`."
75        );
76        let dennorm = self.density_matrix().norm_l2();
77        ensure!(
78            approx::abs_diff_ne!(
79                dennorm,
80                <T as ndarray_linalg::Scalar>::Real::zero(),
81                epsilon = self.threshold()
82            ),
83            "Zero density (density matrix has Frobenius norm {dennorm:.3e})"
84        );
85
86        let sao_4c = metric
87            .ok_or_else(|| format_err!("No atomic-orbital four-centre overlap tensor found for density overlap calculation."))?;
88        let sao_4c_h = metric_h.unwrap_or(sao_4c);
89
90        if self.complex_symmetric() {
91            match (self.complex_conjugated, other.complex_conjugated) {
92                (false, false) => einsum(
93                    "ijkl,ji,lk->",
94                    &[
95                        &sao_4c_h.view(),
96                        &self.density_matrix().view(),
97                        &other.density_matrix().view(),
98                    ],
99                )
100                .map_err(|err| format_err!(err))?
101                .into_iter()
102                .next()
103                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
104                (true, false) => einsum(
105                    "ijkl,ji,lk->",
106                    &[
107                        &sao_4c.view(),
108                        &self.density_matrix().view(),
109                        &other.density_matrix().view(),
110                    ],
111                )
112                .map_err(|err| format_err!(err))?
113                .into_iter()
114                .next()
115                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
116                (false, true) => einsum(
117                    "ijkl,ji,lk->",
118                    &[
119                        &sao_4c.view(),
120                        &other.density_matrix().view(),
121                        &self.density_matrix().view(),
122                    ],
123                )
124                .map_err(|err| format_err!(err))?
125                .into_iter()
126                .next()
127                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
128                (true, true) => einsum(
129                    "klij,ji,lk->",
130                    &[
131                        &sao_4c_h.view(),
132                        &self.density_matrix().view(),
133                        &other.density_matrix().view(),
134                    ],
135                )
136                .map_err(|err| format_err!(err))?
137                .into_iter()
138                .next()
139                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
140            }
141        } else {
142            match (self.complex_conjugated, other.complex_conjugated) {
143                (false, false) => einsum(
144                    "ijkl,ji,lk->",
145                    &[
146                        &sao_4c.view(),
147                        &self.density_matrix().mapv(|x| x.conj()).view(),
148                        &other.density_matrix().view(),
149                    ],
150                )
151                .map_err(|err| format_err!(err))?
152                .into_iter()
153                .next()
154                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
155                (true, false) => einsum(
156                    "ijkl,ji,lk->",
157                    &[
158                        &sao_4c_h.view(),
159                        &self.density_matrix().mapv(|x| x.conj()).view(),
160                        &other.density_matrix().view(),
161                    ],
162                )
163                .map_err(|err| format_err!(err))?
164                .into_iter()
165                .next()
166                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
167                (false, true) => einsum(
168                    "ijkl,ji,lk->",
169                    &[
170                        &sao_4c_h.view(),
171                        &other.density_matrix().mapv(|x| x.conj()).view(),
172                        &self.density_matrix().view(),
173                    ],
174                )
175                .map_err(|err| format_err!(err))?
176                .into_iter()
177                .next()
178                .ok_or(format_err!("Unable to extract the density overlap scalar."))
179                .map(|x| x.conj()),
180                (true, true) => einsum(
181                    "klij,ji,lk->",
182                    &[
183                        &sao_4c.view(),
184                        &self.density_matrix().mapv(|x| x.conj()).view(),
185                        &other.density_matrix().view(),
186                    ],
187                )
188                .map_err(|err| format_err!(err))?
189                .into_iter()
190                .next()
191                .ok_or(format_err!("Unable to extract the density overlap scalar.")),
192            }
193        }
194    }
195
196    /// Returns the mathematical definition of the overlap between two densities.
197    fn overlap_definition(&self) -> String {
198        let k = if self.complex_symmetric() { "κ " } else { "" };
199        format!("⟨{k}ρ_1|ρ_2⟩ = ∫ [{k}ρ_1(r)]* ρ_2(r) dr")
200    }
201}
202
203// ====================
204// DensitySymmetryOrbit
205// ====================
206
207// -----------------
208// Struct definition
209// -----------------
210
211/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of
212/// densities.
213#[derive(Builder, Clone)]
214pub struct DensitySymmetryOrbit<'a, G, T>
215where
216    G: SymmetryGroupProperties,
217    T: ComplexFloat + fmt::Debug + Lapack,
218    Density<'a, T>: SymmetryTransformable,
219{
220    /// The generating symmetry group.
221    group: &'a G,
222
223    /// The origin density of the orbit.
224    origin: &'a Density<'a, T>,
225
226    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
227    pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
228
229    /// The threshold for determining if calculated multiplicities in representation analysis are
230    /// integral.
231    integrality_threshold: <T as ComplexFloat>::Real,
232
233    /// The kind of transformation determining the way the symmetry operations in `group` act on
234    /// [`Self::origin`].
235    symmetry_transformation_kind: SymmetryTransformationKind,
236
237    /// The overlap matrix between the symmetry-equivalent densities in the orbit.
238    #[builder(setter(skip), default = "None")]
239    smat: Option<Array2<T>>,
240
241    /// The eigenvalues of the overlap matrix between the symmetry-equivalent densities in
242    /// the orbit.
243    #[builder(setter(skip), default = "None")]
244    pub(crate) smat_eigvals: Option<Array1<T>>,
245
246    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent densities
247    /// in the orbit.
248    ///
249    /// See [`RepAnalysis::xmat`] for further information.
250    #[builder(setter(skip), default = "None")]
251    xmat: Option<Array2<T>>,
252
253    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
254    /// eigenvalues.
255    eigenvalue_comparison_mode: EigenvalueComparisonMode,
256}
257
258// ----------------------------
259// Struct method implementation
260// ----------------------------
261
262impl<'a, G, T> DensitySymmetryOrbit<'a, G, T>
263where
264    G: SymmetryGroupProperties + Clone,
265    T: ComplexFloat + fmt::Debug + Lapack,
266    Density<'a, T>: SymmetryTransformable,
267{
268    /// Returns a builder for constructing a new density symmetry orbit.
269    pub fn builder() -> DensitySymmetryOrbitBuilder<'a, G, T> {
270        DensitySymmetryOrbitBuilder::default()
271    }
272}
273
274impl<'a, G> DensitySymmetryOrbit<'a, G, f64>
275where
276    G: SymmetryGroupProperties,
277{
278    fn_calc_xmat_real!(
279        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
280        /// $`\mathbf{S}`$ between the symmetry-equivalent densities in the orbit.
281        ///
282        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
283        ///
284        /// # Arguments
285        ///
286        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
287        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit densities. If
288        /// `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
289        /// full rank.
290        pub calc_xmat
291    );
292}
293
294impl<'a, G, T> DensitySymmetryOrbit<'a, G, Complex<T>>
295where
296    G: SymmetryGroupProperties,
297    T: Float + Scalar<Complex = Complex<T>>,
298    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
299    Density<'a, Complex<T>>: SymmetryTransformable + Overlap<Complex<T>, Ix4>,
300{
301    fn_calc_xmat_complex!(
302        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
303        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent densities in the orbit.
304        ///
305        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
306        ///
307        /// # Arguments
308        ///
309        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
310        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit densities. If
311        /// `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
312        /// full rank.
313        pub calc_xmat
314    );
315}
316
317// ---------------------
318// Trait implementations
319// ---------------------
320
321// ~~~~~
322// Orbit
323// ~~~~~
324
325impl<'a, G, T> Orbit<G, Density<'a, T>> for DensitySymmetryOrbit<'a, G, T>
326where
327    G: SymmetryGroupProperties,
328    T: ComplexFloat + fmt::Debug + Lapack,
329    Density<'a, T>: SymmetryTransformable,
330{
331    type OrbitIter = OrbitIterator<'a, G, Density<'a, T>>;
332
333    fn group(&self) -> &G {
334        self.group
335    }
336
337    fn origin(&self) -> &Density<'a, T> {
338        self.origin
339    }
340
341    fn iter(&self) -> Self::OrbitIter {
342        OrbitIterator::new(
343            self.group,
344            self.origin,
345            match self.symmetry_transformation_kind {
346                SymmetryTransformationKind::Spatial => |op, det| {
347                    det.sym_transform_spatial(op).with_context(|| {
348                        format!("Unable to apply `{op}` spatially on the origin density")
349                    })
350                },
351                SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, det| {
352                    det.sym_transform_spatial_with_spintimerev(op).with_context(|| {
353                        format!("Unable to apply `{op}` spatially (with spin-including time-reversal) on the origin density")
354                    })
355                },
356                SymmetryTransformationKind::Spin => |op, det| {
357                    det.sym_transform_spin(op).with_context(|| {
358                        format!("Unable to apply `{op}` spin-wise on the origin density")
359                    })
360                },
361                SymmetryTransformationKind::SpinSpatial => |op, det| {
362                    det.sym_transform_spin_spatial(op).with_context(|| {
363                        format!("Unable to apply `{op}` spin-spatially on the origin density")
364                    })
365                },
366            },
367        )
368    }
369}
370
371// ~~~~~~~~~~~
372// RepAnalysis
373// ~~~~~~~~~~~
374
375impl<'a, G, T> RepAnalysis<G, Density<'a, T>, T, Ix4> for DensitySymmetryOrbit<'a, G, T>
376where
377    G: SymmetryGroupProperties,
378    G::CharTab: SubspaceDecomposable<T>,
379    T: Lapack
380        + ComplexFloat<Real = <T as Scalar>::Real>
381        + fmt::Debug
382        + Mul<<T as ComplexFloat>::Real, Output = T>,
383    <T as ComplexFloat>::Real: fmt::Debug
384        + Zero
385        + From<u16>
386        + ToPrimitive
387        + approx::RelativeEq<<T as ComplexFloat>::Real>
388        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
389    Density<'a, T>: SymmetryTransformable,
390{
391    fn set_smat(&mut self, smat: Array2<T>) {
392        self.smat = Some(smat)
393    }
394
395    fn smat(&self) -> Option<&Array2<T>> {
396        self.smat.as_ref()
397    }
398
399    fn xmat(&self) -> &Array2<T> {
400        self.xmat
401            .as_ref()
402            .expect("Orbit overlap orthogonalisation matrix not found.")
403    }
404
405    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
406        if self.origin.complex_symmetric {
407            Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
408        } else {
409            if self
410                .group
411                .get_index(i)
412                .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
413                .contains_time_reversal()
414            {
415                Ok(ComplexFloat::conj)
416            } else {
417                Ok(|x| x)
418            }
419        }
420    }
421
422    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
423        self.integrality_threshold
424    }
425
426    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
427        &self.eigenvalue_comparison_mode
428    }
429
430    /// Reduces the representation or corepresentation spanned by the densities in the orbit to
431    /// a direct sum of the irreducible representations or corepresentations of the generating
432    /// symmetry group.
433    ///
434    /// # Returns
435    ///
436    /// The decomposed result.
437    ///
438    /// # Errors
439    ///
440    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
441    /// are non-integral.
442    fn analyse_rep(
443        &self,
444    ) -> Result<
445        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
446        DecompositionError,
447    > {
448        log::debug!("Analysing representation symmetry for a density...");
449        let chis = self
450            .calc_characters()
451            .map_err(|err| DecompositionError(err.to_string()))?;
452        log::debug!("Characters calculated.");
453        let res = self.group().character_table().reduce_characters(
454            &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
455            self.integrality_threshold(),
456        );
457        log::debug!("Characters reduced.");
458        log::debug!("Analysing representation symmetry for a density... Done.");
459        res
460    }
461}