qsym2/target/density/
density_transformation.rs

1//! Implementation of symmetry transformations for electron densities.
2
3use ndarray::{concatenate, s, Array2, Axis, LinalgScalar, ScalarOperand};
4use ndarray_linalg::types::Lapack;
5use num_complex::{Complex, ComplexFloat};
6
7use crate::permutation::{IntoPermutation, PermutableCollection, Permutation};
8use crate::symmetry::symmetry_element::SymmetryOperation;
9use crate::symmetry::symmetry_transformation::{
10    assemble_sh_rotation_3d_matrices, permute_array_by_atoms, ComplexConjugationTransformable,
11    DefaultTimeReversalTransformable, SpatialUnitaryTransformable, SpinUnitaryTransformable,
12    SymmetryTransformable, TimeReversalTransformable, TransformationError,
13};
14use crate::target::density::Density;
15
16// ---------------------------
17// SpatialUnitaryTransformable
18// ---------------------------
19impl<'a, T> SpatialUnitaryTransformable for Density<'a, T>
20where
21    T: ComplexFloat + LinalgScalar + ScalarOperand + Copy + Lapack,
22    f64: Into<T>,
23{
24    fn transform_spatial_mut(
25        &mut self,
26        rmat: &Array2<f64>,
27        perm: Option<&Permutation<usize>>,
28    ) -> Result<&mut Self, TransformationError> {
29        let tmats: Vec<Array2<T>> = assemble_sh_rotation_3d_matrices(self.bao, rmat, perm)
30            .map_err(|err| TransformationError(err.to_string()))?
31            .iter()
32            .map(|tmat| tmat.map(|&x| x.into()))
33            .collect();
34        let pbao = if let Some(p) = perm {
35            self.bao
36                .permute(p)
37                .map_err(|err| TransformationError(err.to_string()))?
38        } else {
39            self.bao.clone()
40        };
41        let old_denmat = &self.density_matrix;
42        let p_coeff = if let Some(p) = perm {
43            permute_array_by_atoms(old_denmat, p, &[Axis(0), Axis(1)], self.bao)
44        } else {
45            old_denmat.clone()
46        };
47        let trow_p_blocks = pbao
48            .shell_boundary_indices()
49            .into_iter()
50            .zip(tmats.iter())
51            .map(|((shl_start, shl_end), tmat)| {
52                tmat.dot(&p_coeff.slice(s![shl_start..shl_end, ..]))
53            })
54            .collect::<Vec<_>>();
55        let trow_p_coeff = concatenate(
56            Axis(0),
57            &trow_p_blocks
58                .iter()
59                .map(|trow_p_block| trow_p_block.view())
60                .collect::<Vec<_>>(),
61        )
62        .expect("Unable to concatenate the transformed rows for the various shells.");
63
64        let tcol_trow_p_blocks = pbao
65            .shell_boundary_indices()
66            .into_iter()
67            .zip(tmats.iter())
68            .map(|((shl_start, shl_end), tmat)| {
69                // tmat is real-valued, so there is no need for tmat.t().map(|x| x.conj()).
70                trow_p_coeff
71                    .slice(s![.., shl_start..shl_end])
72                    .dot(&tmat.t())
73            })
74            .collect::<Vec<_>>();
75        let new_denmat = concatenate(
76            Axis(1),
77            &tcol_trow_p_blocks
78                .iter()
79                .map(|tcol_trow_p_block| tcol_trow_p_block.view())
80                .collect::<Vec<_>>(),
81        )
82        .expect("Unable to concatenate the transformed columns for the various shells.");
83        self.density_matrix = new_denmat;
84        Ok(self)
85    }
86}
87
88// ------------------------
89// SpinUnitaryTransformable
90// ------------------------
91
92impl<'a, T> SpinUnitaryTransformable for Density<'a, T>
93where
94    T: ComplexFloat + Lapack,
95{
96    /// Performs a spin transformation in-place.
97    ///
98    /// Since densities are entirely spatial, spin transformations have no effect on them. This
99    /// thus simply returns `self` without modification.
100    fn transform_spin_mut(
101        &mut self,
102        _: &Array2<Complex<f64>>,
103    ) -> Result<&mut Self, TransformationError> {
104        Ok(self)
105    }
106}
107
108// -------------------------------
109// ComplexConjugationTransformable
110// -------------------------------
111
112impl<'a, T> ComplexConjugationTransformable for Density<'a, T>
113where
114    T: ComplexFloat + Lapack,
115{
116    fn transform_cc_mut(&mut self) -> Result<&mut Self, TransformationError> {
117        self.density_matrix.mapv_inplace(|x| x.conj());
118        self.complex_conjugated = !self.complex_conjugated;
119        Ok(self)
120    }
121}
122
123// --------------------------------
124// DefaultTimeReversalTransformable
125// --------------------------------
126impl<'a, T> DefaultTimeReversalTransformable for Density<'a, T> where T: ComplexFloat + Lapack {}
127
128// ---------------------
129// SymmetryTransformable
130// ---------------------
131impl<'a, T> SymmetryTransformable for Density<'a, T>
132where
133    T: ComplexFloat + Lapack,
134    Density<'a, T>: SpatialUnitaryTransformable + TimeReversalTransformable,
135{
136    fn sym_permute_sites_spatial(
137        &self,
138        symop: &SymmetryOperation,
139    ) -> Result<Permutation<usize>, TransformationError> {
140        if (symop.generating_element.threshold().log10() - self.mol.threshold.log10()).abs() >= 3.0
141        {
142            log::warn!(
143                "Symmetry operation threshold ({:.3e}) and molecule threshold ({:.3e}) \
144                differ by more than three orders of magnitudes.",
145                symop.generating_element.threshold(),
146                self.mol.threshold
147            )
148        }
149        symop
150            .act_permute(&self.mol.molecule_ordinary_atoms())
151            .ok_or(TransformationError(format!(
152            "Unable to determine the atom permutation corresponding to the operation `{symop}`.",
153        )))
154    }
155}
qsym2/target/density/density_transformation.rs

qsym2/target/density/
density_transformation.rs
