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}