qsym2/target/density/
density_transformation.rs

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

qsym2/target/density/
density_transformation.rs
