qsym2/target/density/
density_projection.rs

1//! Implementation of symmetry projection for electron densities.
2
3use anyhow::{self, format_err};
4use ndarray::Array2;
5use num::{Complex, ToPrimitive};
6
7use super::density_analysis::DensitySymmetryOrbit;
8use crate::analysis::Orbit;
9use crate::chartab::chartab_symbols::LinearSpaceSymbol;
10use crate::projection::Projectable;
11use crate::symmetry::symmetry_group::SymmetryGroupProperties;
12use crate::target::density::Density;
13
14impl<'a, G> Projectable<G, Density<'a, f64>> for DensitySymmetryOrbit<'a, G, f64>
15where
16    G: SymmetryGroupProperties,
17{
18    type Projected<'p>
19        = Result<Density<'p, f64>, anyhow::Error>
20    where
21        Self: 'p;
22
23    fn project_onto(&self, row: &G::RowSymbol) -> Self::Projected<'_> {
24        let group_order = self
25            .group()
26            .order()
27            .to_f64()
28            .ok_or_else(|| format_err!("Unable to convert the group order to `f64`."))?;
29        let dim_f64 = row
30            .dimensionality()
31            .to_f64()
32            .ok_or_else(|| format_err!("Unable to convert the degeneracy to `f64`."))?;
33        let denmat_res = self.generate_orbit_algebra_terms(row).fold(
34            Ok(Array2::<f64>::zeros(self.origin().density_matrix().raw_dim())),
35            |acc_res, item_res| acc_res.and_then(|acc| item_res.and_then(|(chr, den)| {
36                let chr_complex = chr.complex_conjugate().complex_value();
37                if chr_complex.im > self.origin().threshold() {
38                    Err(format_err!("Complex characters encountered. Density projection fails over the reals."))
39                } else {
40                    Ok(acc + dim_f64 * chr_complex.re / group_order * den.density_matrix())
41                }
42            })),
43        );
44        denmat_res.and_then(|denmat| {
45            Density::builder()
46                .bao(self.origin().bao)
47                .complex_symmetric(self.origin().complex_symmetric)
48                .complex_conjugated(self.origin().complex_conjugated)
49                .mol(self.origin().mol)
50                .density_matrix(denmat)
51                .threshold(self.origin().threshold)
52                .build()
53                .map_err(|err| format_err!(err))
54        })
55    }
56}
57
58impl<'a, G> Projectable<G, Density<'a, Complex<f64>>> for DensitySymmetryOrbit<'a, G, Complex<f64>>
59where
60    G: SymmetryGroupProperties,
61{
62    type Projected<'p>
63        = Result<Density<'p, Complex<f64>>, anyhow::Error>
64    where
65        Self: 'p;
66
67    fn project_onto(&self, row: &G::RowSymbol) -> Self::Projected<'_> {
68        let group_order = Complex::from(
69            self.group()
70                .order()
71                .to_f64()
72                .ok_or_else(|| format_err!("Unable to convert the group order to `f64`."))?,
73        );
74        let dim_f64 = row
75            .dimensionality()
76            .to_f64()
77            .ok_or_else(|| format_err!("Unable to convert the degeneracy to `f64`."))?;
78        let denmat_res = self.generate_orbit_algebra_terms(row).fold(
79            Ok(Array2::<Complex<f64>>::zeros(
80                self.origin().density_matrix().raw_dim(),
81            )),
82            |acc_res, item_res| {
83                acc_res.and_then(|acc| {
84                    item_res.and_then(|(chr, den)| {
85                        let chr_complex = chr.complex_conjugate().complex_value();
86                        Ok(acc + dim_f64 * chr_complex / group_order * den.density_matrix())
87                    })
88                })
89            },
90        );
91        denmat_res.and_then(|denmat| {
92            Density::builder()
93                .bao(self.origin().bao)
94                .complex_symmetric(self.origin().complex_symmetric)
95                .complex_conjugated(self.origin().complex_conjugated)
96                .mol(self.origin().mol)
97                .density_matrix(denmat)
98                .threshold(self.origin().threshold)
99                .build()
100                .map_err(|err| format_err!(err))
101        })
102    }
103}
qsym2/target/density/density_projection.rs

qsym2/target/density/
density_projection.rs
