qsym2/bindings/python/representation_analysis/
density.rs

1//! Python bindings for QSym² symmetry analysis of electron densities.
2
3use std::path::PathBuf;
4
5use anyhow::format_err;
6use ndarray::{Array2, Array4};
7use num_complex::Complex;
8use numpy::{PyArray2, PyArrayMethods};
9use pyo3::exceptions::{PyIOError, PyRuntimeError};
10use pyo3::prelude::*;
11
12use crate::analysis::EigenvalueComparisonMode;
13use crate::auxiliary::molecule::Molecule;
14use crate::basis::ao::BasisAngularOrder;
15use crate::bindings::python::integrals::PyBasisAngularOrder;
16use crate::bindings::python::representation_analysis::PyArray4RC;
17use crate::drivers::representation_analysis::angular_function::AngularFunctionRepAnalysisParams;
18use crate::drivers::representation_analysis::density::{
19    DensityRepAnalysisDriver, DensityRepAnalysisParams,
20};
21use crate::drivers::representation_analysis::{
22    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
23};
24use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
25use crate::drivers::QSym2Driver;
26use crate::io::format::qsym2_output;
27use crate::io::{read_qsym2_binary, QSym2FileType};
28use crate::symmetry::symmetry_group::{
29    MagneticRepresentedSymmetryGroup, UnitaryRepresentedSymmetryGroup,
30};
31use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
32use crate::target::density::Density;
33
34type C128 = Complex<f64>;
35
36// ==================
37// Struct definitions
38// ==================
39
40/// Python-exposed structure to marshall real electron density information between Rust and
41/// Python.
42///
43/// # Constructor arguments
44///
45/// * `complex_symmetric` - A boolean indicating if inner products involving this density
46/// are complex-symmetric. Python type: `bool`.
47/// * `density_matrix` - The real density matrix describing this density.
48/// Python type: `numpy.2darray[float]`.
49/// * `threshold` - The threshold for comparisons. Python type: `float`.
50#[pyclass]
51#[derive(Clone)]
52pub struct PyDensityReal {
53    /// A boolean indicating if inner products involving this density should be the
54    /// complex-symmetric bilinear form, rather than the conventional Hermitian sesquilinear form.
55    ///
56    /// Python type: `bool`.
57    complex_symmetric: bool,
58
59    /// The real density matrix describing this density.
60    ///
61    /// Python type: `numpy.2darray[float]`.
62    density_matrix: Array2<f64>,
63
64    /// The threshold for comparing densities.
65    ///
66    /// Python type: `float`.
67    threshold: f64,
68}
69
70#[pymethods]
71impl PyDensityReal {
72    /// Constructs a real Python-exposed electron density.
73    ///
74    /// # Arguments
75    ///
76    /// * `complex_symmetric` - A boolean indicating if inner products involving this density
77    /// are complex-symmetric. Python type: `bool`.
78    /// * `density_matrix` - The real density matrix describing this density.
79    /// Python type: `numpy.2darray[float]`.
80    /// * `threshold` - The threshold for comparisons. Python type: `float`.
81    #[new]
82    fn new(
83        complex_symmetric: bool,
84        density_matrix: Bound<'_, PyArray2<f64>>,
85        threshold: f64,
86    ) -> Self {
87        let det = Self {
88            complex_symmetric,
89            density_matrix: density_matrix.to_owned_array(),
90            threshold,
91        };
92        det
93    }
94}
95
96impl PyDensityReal {
97    /// Extracts the information in the [`PyDensityReal`] structure into `QSym2`'s native
98    /// [`Density`] structure.
99    ///
100    /// # Arguments
101    ///
102    /// * `bao` - The [`BasisAngularOrder`] for the basis set in which the density is given.
103    /// * `mol` - The molecule with which the density is associated.
104    ///
105    /// # Returns
106    ///
107    /// The [`Density`] structure with the same information.
108    ///
109    /// # Errors
110    ///
111    /// Errors if the [`Density`] fails to build.
112    fn to_qsym2<'b, 'a: 'b>(
113        &'b self,
114        bao: &'a BasisAngularOrder,
115        mol: &'a Molecule,
116    ) -> Result<Density<'b, f64>, anyhow::Error> {
117        let den = Density::<f64>::builder()
118            .bao(bao)
119            .complex_symmetric(self.complex_symmetric)
120            .mol(mol)
121            .density_matrix(self.density_matrix.clone())
122            .threshold(self.threshold)
123            .build()
124            .map_err(|err| format_err!(err));
125        den
126    }
127}
128
129/// Python-exposed structure to marshall complex electron density information between Rust and
130/// Python.
131///
132/// # Constructor arguments
133///
134/// * `complex_symmetric` - A boolean indicating if inner products involving this density
135/// are complex-symmetric. Python type: `bool`.
136/// * `density_matrix` - The complex density matrix describing this density.
137/// Python type: `numpy.2darray[complex]`.
138/// * `threshold` - The threshold for comparisons. Python type: `float`.
139#[pyclass]
140#[derive(Clone)]
141pub struct PyDensityComplex {
142    /// A boolean indicating if inner products involving this density should be the
143    /// complex-symmetric bilinear form, rather than the conventional Hermitian sesquilinear form.
144    ///
145    /// Python type: `bool`.
146    complex_symmetric: bool,
147
148    /// The complex density matrix describing this density.
149    ///
150    /// Python type: `numpy.2darray[complex]`.
151    density_matrix: Array2<C128>,
152
153    /// The threshold for comparing densities.
154    ///
155    /// Python type: `float`.
156    threshold: f64,
157}
158
159#[pymethods]
160impl PyDensityComplex {
161    /// Constructs a complex Python-exposed electron density.
162    ///
163    /// # Arguments
164    ///
165    /// * `complex_symmetric` - A boolean indicating if inner products involving this density
166    /// are complex-symmetric. Python type: `bool`.
167    /// * `density_matrix` - The complex density matrix describing this density.
168    /// Python type: `numpy.2darray[complex]`.
169    /// * `threshold` - The threshold for comparisons. Python type: `float`.
170    #[new]
171    fn new(
172        complex_symmetric: bool,
173        density_matrix: Bound<'_, PyArray2<C128>>,
174        threshold: f64,
175    ) -> Self {
176        let det = Self {
177            complex_symmetric,
178            density_matrix: density_matrix.to_owned_array(),
179            threshold,
180        };
181        det
182    }
183}
184
185impl PyDensityComplex {
186    /// Extracts the information in the [`PyDensityComplex`] structure into `QSym2`'s native
187    /// [`Density`] structure.
188    ///
189    /// # Arguments
190    ///
191    /// * `bao` - The [`BasisAngularOrder`] for the basis set in which the density is given.
192    /// * `mol` - The molecule with which the density is associated.
193    ///
194    /// # Returns
195    ///
196    /// The [`Density`] structure with the same information.
197    ///
198    /// # Errors
199    ///
200    /// Errors if the [`Density`] fails to build.
201    fn to_qsym2<'b, 'a: 'b>(
202        &'b self,
203        bao: &'a BasisAngularOrder,
204        mol: &'a Molecule,
205    ) -> Result<Density<'b, C128>, anyhow::Error> {
206        let den = Density::<C128>::builder()
207            .bao(bao)
208            .complex_symmetric(self.complex_symmetric)
209            .mol(mol)
210            .density_matrix(self.density_matrix.clone())
211            .threshold(self.threshold)
212            .build()
213            .map_err(|err| format_err!(err));
214        den
215    }
216}
217
218// ================
219// Enum definitions
220// ================
221
222/// Python-exposed enumerated type to handle the union type `PyDensityReal | PyDensityComplex` in
223/// Python.
224#[derive(FromPyObject)]
225pub enum PyDensity {
226    /// Variant for real Python-exposed electron density.
227    Real(PyDensityReal),
228
229    /// Variant for complex Python-exposed electron density.
230    Complex(PyDensityComplex),
231}
232
233// =====================
234// Functions definitions
235// =====================
236
237/// Python-exposed function to perform representation symmetry analysis for real and complex
238/// electron densities and log the result via the `qsym2-output` logger at the `INFO` level.
239///
240/// # Arguments
241///
242/// * `inp_sym` - A path to the [`QSym2FileType::Sym`] file containing the symmetry-group detection
243/// result for the system. This will be used to construct abstract groups and character tables for
244/// representation analysis. Python type: `str`.
245/// * `pydens` - A sequence of Python-exposed electron densities whose density matrices are of type
246/// `float64` or `complex128`. Each density is accompanied by a description string.
247/// Python type: `list[tuple[str, PyDensityReal | PyDensityComplex]]`.
248/// * `sao_spatial_4c` - The atomic-orbital four-centre overlap matrix whose elements are of type
249/// `float64` or `complex128`. Python type: `numpy.4darray[float] | numpy.4darray[complex]`.
250/// * `sao_spatial_4c_h` - The optional complex-symmetric atomic-orbital four-centre overlap matrix
251/// whose elements are of type `float64` or `complex128`. This is required if antiunitary symmetry
252/// operations are involved. Python type: `numpy.2darray[float] | numpy.2darray[complex] | None`.
253/// * `integrality_threshold` - The threshold for verifying if subspace multiplicities are
254/// integral. Python type: `float`.
255/// * `linear_independence_threshold` - The threshold for determining the linear independence
256/// subspace via the non-zero eigenvalues of the orbit overlap matrix. Python type: `float`.
257/// * `use_magnetic_group` - An option indicating if the magnetic group is to be used for symmetry
258/// analysis, and if so, whether unitary representations or unitary-antiunitary corepresentations
259/// should be used. Python type: `None | MagneticSymmetryAnalysisKind`.
260/// * `use_double_group` - A boolean indicating if the double group of the prevailing symmetry
261/// group is to be used for representation analysis instead. Python type: `bool`.
262/// * `use_cayley_table` - A boolean indicating if the Cayley table for the group, if available,
263/// should be used to speed up the calculation of orbit overlap matrices. Python type: `bool`.
264/// * `symmetry_transformation_kind` - An enumerated type indicating the type of symmetry
265/// transformations to be performed on the origin electron density to generate the orbit. Python
266/// type: `SymmetryTransformationKind`.
267/// * `eigenvalue_comparison_mode` - An enumerated type indicating the mode of comparison of orbit
268/// overlap eigenvalues with the specified `linear_independence_threshold`.
269/// Python type: `EigenvalueComparisonMode`.
270/// * `write_character_table` - A boolean indicating if the character table of the prevailing
271/// symmetry group is to be printed out. Python type: `bool`.
272/// * `infinite_order_to_finite` - The finite order with which infinite-order generators are to be
273/// interpreted to form a finite subgroup of the prevailing infinite group. This finite subgroup
274/// will be used for symmetry analysis. Python type: `Optional[int]`.
275/// * `angular_function_integrality_threshold` - The threshold for verifying if subspace
276/// multiplicities are integral for the symmetry analysis of angular functions. Python type:
277/// `float`.
278/// * `angular_function_linear_independence_threshold` - The threshold for determining the linear
279/// independence subspace via the non-zero eigenvalues of the orbit overlap matrix for the symmetry
280/// analysis of angular functions. Python type: `float`.
281/// * `angular_function_max_angular_momentum` - The maximum angular momentum order to be used in
282/// angular function symmetry analysis. Python type: `int`.
283#[pyfunction]
284#[pyo3(signature = (
285    inp_sym,
286    pydens,
287    pybao,
288    integrality_threshold,
289    linear_independence_threshold,
290    use_magnetic_group,
291    use_double_group,
292    use_cayley_table,
293    symmetry_transformation_kind,
294    eigenvalue_comparison_mode,
295    sao_spatial_4c,
296    sao_spatial_4c_h=None,
297    write_character_table=true,
298    infinite_order_to_finite=None,
299    angular_function_integrality_threshold=1e-7,
300    angular_function_linear_independence_threshold=1e-7,
301    angular_function_max_angular_momentum=2
302))]
303pub fn rep_analyse_densities(
304    py: Python<'_>,
305    inp_sym: PathBuf,
306    pydens: Vec<(String, PyDensity)>,
307    pybao: &PyBasisAngularOrder,
308    integrality_threshold: f64,
309    linear_independence_threshold: f64,
310    use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
311    use_double_group: bool,
312    use_cayley_table: bool,
313    symmetry_transformation_kind: SymmetryTransformationKind,
314    eigenvalue_comparison_mode: EigenvalueComparisonMode,
315    sao_spatial_4c: PyArray4RC,
316    sao_spatial_4c_h: Option<PyArray4RC>,
317    write_character_table: bool,
318    infinite_order_to_finite: Option<u32>,
319    angular_function_integrality_threshold: f64,
320    angular_function_linear_independence_threshold: f64,
321    angular_function_max_angular_momentum: u32,
322) -> PyResult<()> {
323    let pd_res: SymmetryGroupDetectionResult =
324        read_qsym2_binary(inp_sym.clone(), QSym2FileType::Sym)
325            .map_err(|err| PyIOError::new_err(err.to_string()))?;
326
327    let mut file_name = inp_sym.to_path_buf();
328    file_name.set_extension(QSym2FileType::Sym.ext());
329    qsym2_output!(
330        "Symmetry-group detection results read in from {}.",
331        file_name.display(),
332    );
333    qsym2_output!("");
334
335    let mol = &pd_res.pre_symmetry.recentred_molecule;
336    let bao = pybao
337        .to_qsym2(mol)
338        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
339    let afa_params = AngularFunctionRepAnalysisParams::builder()
340        .integrality_threshold(angular_function_integrality_threshold)
341        .linear_independence_threshold(angular_function_linear_independence_threshold)
342        .max_angular_momentum(angular_function_max_angular_momentum)
343        .build()
344        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
345    let sda_params = DensityRepAnalysisParams::<f64>::builder()
346        .integrality_threshold(integrality_threshold)
347        .linear_independence_threshold(linear_independence_threshold)
348        .use_magnetic_group(use_magnetic_group.clone())
349        .use_double_group(use_double_group)
350        .use_cayley_table(use_cayley_table)
351        .symmetry_transformation_kind(symmetry_transformation_kind)
352        .eigenvalue_comparison_mode(eigenvalue_comparison_mode)
353        .write_character_table(if write_character_table {
354            Some(CharacterTableDisplay::Symbolic)
355        } else {
356            None
357        })
358        .infinite_order_to_finite(infinite_order_to_finite)
359        .build()
360        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
361
362    let any_complex = pydens
363        .iter()
364        .any(|(_, pyden)| matches!(pyden, PyDensity::Complex(_)));
365
366    match (any_complex, &sao_spatial_4c) {
367        (false, PyArray4RC::Real(pysao4c_r)) => {
368            // Both coefficients and sao_4c are real.
369            let dens = pydens
370                .iter()
371                .map(|(_, pyden)| match pyden {
372                    PyDensity::Real(pyden_r) => pyden_r.to_qsym2(&bao, mol),
373                    PyDensity::Complex(_) => panic!("Unexpected complex density."),
374                })
375                .collect::<Result<Vec<_>, _>>()
376                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
377            let dens_ref = dens
378                .iter()
379                .zip(pydens.iter())
380                .map(|(den, (desc, _))| (desc.clone(), den))
381                .collect::<Vec<_>>();
382
383            let sao_spatial_4c = pysao4c_r.to_owned_array();
384
385            match &use_magnetic_group {
386                Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
387                    let mut da_driver =
388                        DensityRepAnalysisDriver::<MagneticRepresentedSymmetryGroup, f64>::builder(
389                        )
390                        .parameters(&sda_params)
391                        .angular_function_parameters(&afa_params)
392                        .densities(dens_ref)
393                        .sao_spatial_4c(&sao_spatial_4c)
394                        .sao_spatial_4c_h(None)
395                        .symmetry_group(&pd_res)
396                        .build()
397                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
398                    py.allow_threads(|| {
399                        da_driver
400                            .run()
401                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
402                    })?
403                }
404                Some(MagneticSymmetryAnalysisKind::Representation) | None => {
405                    let mut da_driver =
406                        DensityRepAnalysisDriver::<UnitaryRepresentedSymmetryGroup, f64>::builder()
407                            .parameters(&sda_params)
408                            .angular_function_parameters(&afa_params)
409                            .densities(dens_ref)
410                            .sao_spatial_4c(&sao_spatial_4c)
411                            .sao_spatial_4c_h(None)
412                            .symmetry_group(&pd_res)
413                            .build()
414                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
415                    py.allow_threads(|| {
416                        da_driver
417                            .run()
418                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
419                    })?
420                }
421            };
422        }
423        (_, _) => {
424            // At least one of coefficients or sao_4c are not real.
425            let dens: Vec<Density<C128>> = pydens
426                .iter()
427                .map(|(_, pyden)| match pyden {
428                    PyDensity::Real(pyden_r) => {
429                        pyden_r.to_qsym2(&bao, mol).map(|den_r| den_r.into())
430                    }
431                    PyDensity::Complex(pyden_c) => pyden_c.to_qsym2(&bao, mol),
432                })
433                .collect::<Result<Vec<_>, _>>()
434                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
435            let dens_ref = dens
436                .iter()
437                .zip(pydens.iter())
438                .map(|(den, (desc, _))| (desc.clone(), den))
439                .collect::<Vec<_>>();
440
441            let (sao_spatial_4c_c, sao_spatial_4c_h_c): (Array4<C128>, Option<Array4<C128>>) =
442                match sao_spatial_4c {
443                    PyArray4RC::Real(pysao4c_r) => {
444                        (pysao4c_r.to_owned_array().mapv(Complex::from), None)
445                    }
446                    PyArray4RC::Complex(pysao4c_c) => (
447                        pysao4c_c.to_owned_array(),
448                        sao_spatial_4c_h.map(|v| match v {
449                            PyArray4RC::Real(pysao4c_h_r) => {
450                                pysao4c_h_r.to_owned_array().mapv(Complex::from)
451                            }
452                            PyArray4RC::Complex(pysao4c_h_c) => pysao4c_h_c.to_owned_array(),
453                        }),
454                    ),
455                };
456
457            match &use_magnetic_group {
458                Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
459                    let mut da_driver = DensityRepAnalysisDriver::<
460                        MagneticRepresentedSymmetryGroup,
461                        C128,
462                    >::builder()
463                    .parameters(&sda_params)
464                    .angular_function_parameters(&afa_params)
465                    .densities(dens_ref)
466                    .sao_spatial_4c(&sao_spatial_4c_c)
467                    .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
468                    .symmetry_group(&pd_res)
469                    .build()
470                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
471                    py.allow_threads(|| {
472                        da_driver
473                            .run()
474                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
475                    })?
476                }
477                Some(MagneticSymmetryAnalysisKind::Representation) | None => {
478                    let mut da_driver =
479                        DensityRepAnalysisDriver::<UnitaryRepresentedSymmetryGroup, C128>::builder(
480                        )
481                        .parameters(&sda_params)
482                        .angular_function_parameters(&afa_params)
483                        .densities(dens_ref)
484                        .sao_spatial_4c(&sao_spatial_4c_c)
485                        .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
486                        .symmetry_group(&pd_res)
487                        .build()
488                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
489                    py.allow_threads(|| {
490                        da_driver
491                            .run()
492                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
493                    })?
494                }
495            };
496        }
497    }
498    Ok(())
499}