qsym2/bindings/python/projection/
density.rs

1//! Python bindings for QSym² symmetry projection of electron densities.
2
3use std::path::PathBuf;
4
5use indexmap::IndexMap;
6use itertools::Itertools;
7use num_complex::Complex;
8use pyo3::exceptions::{PyIOError, PyRuntimeError};
9use pyo3::{IntoPyObjectExt, prelude::*};
10
11use crate::bindings::python::integrals::PyBasisAngularOrder;
12use crate::bindings::python::projection::PyProjectionTarget;
13use crate::bindings::python::representation_analysis::density::{
14    PyDensity, PyDensityComplex, PyDensityReal,
15};
16use crate::drivers::QSym2Driver;
17use crate::drivers::projection::density::{DensityProjectionDriver, DensityProjectionParams};
18use crate::drivers::representation_analysis::{
19    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
20};
21use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
22use crate::io::format::qsym2_output;
23use crate::io::{QSym2FileType, read_qsym2_binary};
24use crate::symmetry::symmetry_group::UnitaryRepresentedSymmetryGroup;
25use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
26use crate::target::density::Density;
27
28type C128 = Complex<f64>;
29
30// =====================
31// Functions definitions
32// =====================
33
34/// Python-exposed function to perform symmetry projection for real and complex electron densities.
35///
36/// # Arguments
37///
38/// * `inp_sym` - A path to the [`QSym2FileType::Sym`] file containing the symmetry-group detection
39/// result for the system. This will be used to construct abstract groups and character tables for
40/// symmetry projection.
41/// * `pydens` - A sequence of Python-exposed electron densities whose density matrices are of type
42/// `float64` or `complex128`. Each density is accompanied by a description string.
43/// * `projection_targets` - A sequence of subspace labels for projection. Each label is either a
44/// symbolic string or a numerical index for the subspace in the character table of the prevailing
45/// group.
46/// * `pybao` - Python-exposed structure containing basis angular order information for the density
47/// matrices.
48/// * `use_magnetic_group` - An option indicating if the magnetic group is to be used for symmetry
49/// analysis, and if so, whether unitary representations or unitary-antiunitary corepresentations
50/// should be used.
51/// * `use_double_group` - A boolean indicating if the double group of the prevailing symmetry
52/// group is to be used for representation analysis instead.
53/// * `symmetry_transformation_kind` - An enumerated type indicating the type of symmetry
54/// transformations to be performed on the origin electron density to generate the orbit.
55/// * `write_character_table` - A boolean indicating if the character table of the prevailing
56/// symmetry group is to be printed out.
57/// * `infinite_order_to_finite` - The finite order with which infinite-order generators are to be
58/// interpreted to form a finite subgroup of the prevailing infinite group. This finite subgroup
59/// will be used for symmetry analysis.
60///
61/// # Returns
62///
63/// The result will be returned as a list of tuples, each of which contains the name/description of
64/// an original density and a dictionary in which the keys are the subspace labels and the values
65/// are the corresponding projected density.
66#[allow(clippy::type_complexity, clippy::too_many_arguments)]
67#[pyfunction]
68#[pyo3(signature = (
69    inp_sym,
70    pydens,
71    projection_targets,
72    pybao,
73    use_magnetic_group,
74    use_double_group,
75    symmetry_transformation_kind,
76    write_character_table=true,
77    infinite_order_to_finite=None,
78))]
79pub fn project_densities(
80    py: Python<'_>,
81    inp_sym: PathBuf,
82    pydens: Vec<(String, PyDensity)>,
83    projection_targets: Vec<PyProjectionTarget>,
84    pybao: &PyBasisAngularOrder,
85    use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
86    use_double_group: bool,
87    symmetry_transformation_kind: SymmetryTransformationKind,
88    write_character_table: bool,
89    infinite_order_to_finite: Option<u32>,
90) -> PyResult<Vec<(String, IndexMap<String, Py<PyAny>>)>> {
91    let pd_res: SymmetryGroupDetectionResult =
92        read_qsym2_binary(inp_sym.clone(), QSym2FileType::Sym)
93            .map_err(|err| PyIOError::new_err(err.to_string()))?;
94
95    let mut file_name = inp_sym.to_path_buf();
96    file_name.set_extension(QSym2FileType::Sym.ext());
97    qsym2_output!(
98        "Symmetry-group detection results read in from {}.",
99        file_name.display(),
100    );
101    qsym2_output!("");
102
103    let mol = &pd_res.pre_symmetry.recentred_molecule;
104    let bao = pybao
105        .to_qsym2(mol)
106        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
107
108    let symbolic_projection_targets = projection_targets
109        .iter()
110        .filter_map(|pt| match pt {
111            PyProjectionTarget::Symbolic(symbolic_pt) => Some(symbolic_pt.clone()),
112            PyProjectionTarget::Numeric(_) => None,
113        })
114        .collect_vec();
115    let numeric_projection_targets = projection_targets
116        .iter()
117        .filter_map(|pt| match pt {
118            PyProjectionTarget::Symbolic(_) => None,
119            PyProjectionTarget::Numeric(numeric_pt) => Some(*numeric_pt),
120        })
121        .collect_vec();
122
123    let dp_params = DensityProjectionParams::builder()
124        .symbolic_projection_targets(Some(symbolic_projection_targets))
125        .numeric_projection_targets(Some(numeric_projection_targets))
126        .use_magnetic_group(use_magnetic_group.clone())
127        .use_double_group(use_double_group)
128        .symmetry_transformation_kind(symmetry_transformation_kind)
129        .write_character_table(if write_character_table {
130            Some(CharacterTableDisplay::Symbolic)
131        } else {
132            None
133        })
134        .infinite_order_to_finite(infinite_order_to_finite)
135        .build()
136        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
137
138    let any_complex = pydens
139        .iter()
140        .any(|(_, pyden)| matches!(pyden, PyDensity::Complex(_)));
141
142    let projected_densities = if !any_complex {
143        // All density matrices are real.
144        let dens = pydens
145            .iter()
146            .map(|(_, pyden)| match pyden {
147                PyDensity::Real(pyden_r) => pyden_r.to_qsym2(&bao, mol),
148                PyDensity::Complex(_) => panic!("Unexpected complex density."),
149            })
150            .collect::<Result<Vec<_>, _>>()
151            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
152        let dens_ref = dens
153            .iter()
154            .zip(pydens.iter())
155            .map(|(den, (desc, _))| (desc.clone(), den))
156            .collect::<Vec<_>>();
157
158        let projected_densities = match &use_magnetic_group {
159            Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
160                return Err(PyRuntimeError::new_err(
161                    "Projection using corepresentations is not yet supported.",
162                ));
163            }
164            Some(MagneticSymmetryAnalysisKind::Representation) | None => {
165                let mut dp_driver =
166                    DensityProjectionDriver::<UnitaryRepresentedSymmetryGroup, f64>::builder()
167                        .parameters(&dp_params)
168                        .densities(dens_ref)
169                        .symmetry_group(&pd_res)
170                        .build()
171                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
172                py.detach(|| {
173                    dp_driver
174                        .run()
175                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))
176                })?;
177                dp_driver
178                    .result()
179                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
180                    .projected_densities()
181                    .clone()
182            }
183        };
184        projected_densities
185            .into_iter()
186            .map(|(name, projected_dens)| {
187                (
188                    name,
189                    projected_dens
190                        .into_iter()
191                        .map(|(row, den_res)| {
192                            (
193                                row.to_string(),
194                                den_res
195                                    .and_then(|den| {
196                                        let pyden = PyDensityReal {
197                                            complex_symmetric: den.complex_symmetric(),
198                                            density_matrix: den.density_matrix().clone(),
199                                            threshold: den.threshold(),
200                                        };
201                                        pyden.into_py_any(py).map_err(|err| err.to_string())
202                                    })
203                                    .expect("Unable to convert a projected density into a Python object."),
204                            )
205                        })
206                        .collect::<IndexMap<_, _>>(),
207                )
208            })
209            .collect::<Vec<_>>()
210    } else {
211        // At least one of coefficients or sao_4c are not real.
212        let dens: Vec<Density<C128>> = pydens
213            .iter()
214            .map(|(_, pyden)| match pyden {
215                PyDensity::Real(pyden_r) => pyden_r.to_qsym2(&bao, mol).map(|den_r| den_r.into()),
216                PyDensity::Complex(pyden_c) => pyden_c.to_qsym2(&bao, mol),
217            })
218            .collect::<Result<Vec<_>, _>>()
219            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
220        let dens_ref = dens
221            .iter()
222            .zip(pydens.iter())
223            .map(|(den, (desc, _))| (desc.clone(), den))
224            .collect::<Vec<_>>();
225
226        let projected_densities = match &use_magnetic_group {
227            Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
228                return Err(PyRuntimeError::new_err(
229                    "Projection using corepresentations is not yet supported.",
230                ));
231            }
232            Some(MagneticSymmetryAnalysisKind::Representation) | None => {
233                let mut dp_driver =
234                    DensityProjectionDriver::<UnitaryRepresentedSymmetryGroup, C128>::builder()
235                        .parameters(&dp_params)
236                        .densities(dens_ref)
237                        .symmetry_group(&pd_res)
238                        .build()
239                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
240                py.detach(|| {
241                    dp_driver
242                        .run()
243                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))
244                })?;
245                dp_driver
246                    .result()
247                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
248                    .projected_densities()
249                    .clone()
250            }
251        };
252        projected_densities
253            .into_iter()
254            .map(|(name, projected_dens)| {
255                (
256                    name,
257                    projected_dens
258                        .into_iter()
259                        .map(|(row, den_res)| {
260                            (
261                                row.to_string(),
262                                den_res
263                                    .and_then(|den| {
264                                        let pyden = PyDensityComplex {
265                                            complex_symmetric: den.complex_symmetric(),
266                                            density_matrix: den.density_matrix().clone(),
267                                            threshold: den.threshold(),
268                                        };
269                                        pyden.into_py_any(py).map_err(|err| err.to_string())
270                                    })
271                                    .expect("Unable to convert a projected density into a Python object."),
272                            )
273                        })
274                        .collect::<IndexMap<_, _>>(),
275                )
276            })
277            .collect::<Vec<_>>()
278    };
279    Ok(projected_densities)
280}
qsym2/bindings/python/projection/density.rs

qsym2/bindings/python/projection/
density.rs
