qsym2/bindings/python/representation_analysis/

slater_determinant.rs

//! Python bindings for QSym² symmetry analysis of Slater determinants.

use std::fmt;
use std::path::PathBuf;

use anyhow::format_err;
use itertools::Itertools;
use ndarray::{Array1, Array2};
use num_complex::Complex;
use numpy::{PyArray1, PyArray2, PyArrayMethods, ToPyArray};
use pyo3::exceptions::{PyIOError, PyRuntimeError};
use pyo3::prelude::*;
use serde::{Deserialize, Serialize};

use crate::analysis::EigenvalueComparisonMode;
use crate::angmom::spinor_rotation_3d::{SpinConstraint, SpinOrbitCoupled, StructureConstraint};
use crate::auxiliary::molecule::Molecule;
use crate::basis::ao::BasisAngularOrder;
use crate::bindings::python::integrals::{PyBasisAngularOrder, PyStructureConstraint};
use crate::bindings::python::representation_analysis::{PyArray2RC, PyArray4RC};
use crate::drivers::QSym2Driver;
use crate::drivers::representation_analysis::angular_function::AngularFunctionRepAnalysisParams;
use crate::drivers::representation_analysis::slater_determinant::{
    SlaterDeterminantRepAnalysisDriver, SlaterDeterminantRepAnalysisParams,
};
use crate::drivers::representation_analysis::{
    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
};
use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
use crate::group::GroupProperties;
use crate::io::format::qsym2_output;
use crate::io::{QSym2FileType, read_qsym2_binary};
use crate::symmetry::symmetry_group::{
    MagneticRepresentedSymmetryGroup, UnitaryRepresentedSymmetryGroup,
};
use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
use crate::target::determinant::SlaterDeterminant;

type C128 = Complex<f64>;

// ==================
// Struct definitions
// ==================

// ------------------
// Slater determinant
// ------------------

// ~~~~
// Real
// ~~~~

/// Python-exposed structure to marshall real Slater determinant information between Rust and
/// Python.
#[pyclass]
#[derive(Clone, Serialize, Deserialize)]
pub struct PySlaterDeterminantReal {
    /// The structure constraint applied to the coefficients of the determinant.
    pub(crate) structure_constraint: PyStructureConstraint,

    /// A boolean indicating if inner products involving this determinant are complex-symmetric.
    #[pyo3(get)]
    pub(crate) complex_symmetric: bool,

    /// The real coefficients for the molecular orbitals of this determinant.
    pub(crate) coefficients: Vec<Array2<f64>>,

    /// The occupation patterns for the molecular orbitals.
    pub(crate) occupations: Vec<Array1<f64>>,

    /// The threshold for comparisons.
    #[pyo3(get)]
    pub(crate) threshold: f64,

    /// The optional real molecular orbital energies.
    pub(crate) mo_energies: Option<Vec<Array1<f64>>>,

    /// The optional real determinantal energy.
    #[pyo3(get)]
    pub(crate) energy: Option<f64>,
}

#[pymethods]
impl PySlaterDeterminantReal {
    /// Constructs a real Python-exposed Slater determinant.
    ///
    /// # Arguments
    ///
    /// * `structure_constraint` - The structure constraint applied to the coefficients of the determinant.
    /// * `complex_symmetric` - A boolean indicating if inner products involving this determinant are complex-symmetric.
    /// * `coefficients` - The real coefficients for the molecular orbitals of this determinant.
    /// * `occupations` - The occupation patterns for the molecular orbitals.
    /// * `threshold` - The threshold for comparisons.
    /// * `mo_energies` - The optional real molecular orbital energies.
    /// * `energy` - The optional real determinantal energy.
    #[new]
    #[pyo3(signature = (structure_constraint, complex_symmetric, coefficients, occupations, threshold, mo_energies=None, energy=None))]
    pub fn new(
        structure_constraint: PyStructureConstraint,
        complex_symmetric: bool,
        coefficients: Vec<Bound<'_, PyArray2<f64>>>,
        occupations: Vec<Bound<'_, PyArray1<f64>>>,
        threshold: f64,
        mo_energies: Option<Vec<Bound<'_, PyArray1<f64>>>>,
        energy: Option<f64>,
    ) -> Self {
        Self {
            structure_constraint,
            complex_symmetric,
            coefficients: coefficients
                .iter()
                .map(|pyarr| pyarr.to_owned_array())
                .collect::<Vec<_>>(),
            occupations: occupations
                .iter()
                .map(|pyarr| pyarr.to_owned_array())
                .collect::<Vec<_>>(),
            threshold,
            mo_energies: mo_energies.map(|energies| {
                energies
                    .iter()
                    .map(|pyarr| pyarr.to_owned_array())
                    .collect::<Vec<_>>()
            }),
            energy,
        }
    }

    /// The occupation patterns for the molecular orbitals.
    #[getter]
    pub fn occupations<'py>(&self, py: Python<'py>) -> PyResult<Vec<Bound<'py, PyArray1<f64>>>> {
        Ok(self
            .occupations
            .iter()
            .map(|occ| occ.to_pyarray(py))
            .collect::<Vec<_>>())
    }

    /// The real coefficients for the molecular orbitals of this determinant.
    #[getter]
    pub fn coefficients<'py>(&self, py: Python<'py>) -> PyResult<Vec<Bound<'py, PyArray2<f64>>>> {
        Ok(self
            .coefficients
            .iter()
            .map(|occ| occ.to_pyarray(py))
            .collect::<Vec<_>>())
    }

    /// The real molecular orbital energies, if any.
    #[getter]
    pub fn mo_energies<'py>(
        &self,
        py: Python<'py>,
    ) -> PyResult<Option<Vec<Bound<'py, PyArray1<f64>>>>> {
        Ok(self.mo_energies.as_ref().map(|mo_energies| {
            mo_energies
                .iter()
                .map(|mo_en| mo_en.to_pyarray(py))
                .collect::<Vec<_>>()
        }))
    }
}

impl PySlaterDeterminantReal {
    /// Extracts the information in the [`PySlaterDeterminantReal`] structure into `QSym2`'s native
    /// [`SlaterDeterminant`] structure.
    ///
    /// # Arguments
    ///
    /// * `baos` - The [`BasisAngularOrder`]s for the basis set in which the Slater determinant is
    ///   given, one for each explicit component per coefficient matrix.
    /// * `mol` - The molecule with which the Slater determinant is associated.
    ///
    /// # Returns
    ///
    /// The [`SlaterDeterminant`] structure with the same information.
    ///
    /// # Errors
    ///
    /// Errors if the [`SlaterDeterminant`] fails to build.
    pub fn to_qsym2<'b, 'a: 'b, SC>(
        &'b self,
        baos: &[&'a BasisAngularOrder],
        mol: &'a Molecule,
    ) -> Result<SlaterDeterminant<'b, f64, SC>, anyhow::Error>
    where
        SC: StructureConstraint
            + Clone
            + fmt::Display
            + TryFrom<PyStructureConstraint, Error = anyhow::Error>,
    {
        SlaterDeterminant::<f64, SC>::builder()
            .structure_constraint(self.structure_constraint.clone().try_into()?)
            .baos(baos.to_vec())
            .complex_symmetric(self.complex_symmetric)
            .mol(mol)
            .coefficients(&self.coefficients)
            .occupations(&self.occupations)
            .mo_energies(self.mo_energies.clone())
            .energy(
                self.energy
                    .ok_or_else(|| "No determinantal energy set.".to_string()),
            )
            .threshold(self.threshold)
            .build()
            .map_err(|err| format_err!(err))
    }

    /// Returns the Python-exposed structure constraint information.
    pub fn structure_constraint(&self) -> &PyStructureConstraint {
        &self.structure_constraint
    }
}

impl<'a, SC> SlaterDeterminant<'a, f64, SC>
where
    SC: StructureConstraint
        + Clone
        + fmt::Display
        + TryInto<PyStructureConstraint, Error = anyhow::Error>,
{
    /// Extracts the information in the real [`SlaterDeterminant`] structure into [`PySlaterDeterminantReal`].
    ///
    /// # Returns
    ///
    /// The [`PySlaterDeterminantReal`] structure with the same information.
    pub fn to_python<'py>(
        &self,
        py: Python<'py>,
    ) -> Result<PySlaterDeterminantReal, anyhow::Error> {
        Ok(PySlaterDeterminantReal::new(
            self.structure_constraint().clone().try_into()?,
            self.complex_symmetric(),
            self.coefficients()
                .iter()
                .map(|coeffs| coeffs.to_pyarray(py))
                .collect_vec(),
            self.occupations()
                .iter()
                .map(|occ| occ.to_pyarray(py))
                .collect_vec(),
            self.threshold(),
            self.mo_energies().map(|mo_energies| {
                mo_energies
                    .iter()
                    .map(|mo_e| mo_e.to_pyarray(py))
                    .collect_vec()
            }),
            self.energy().ok().copied(),
        ))
    }
}

// ~~~~~~~
// Complex
// ~~~~~~~

/// Python-exposed structure to marshall complex Slater determinant information between Rust and Python.
#[pyclass]
#[derive(Clone, Serialize, Deserialize)]
pub struct PySlaterDeterminantComplex {
    /// The structure constraint applied to the coefficients of the determinant.
    pub(crate) structure_constraint: PyStructureConstraint,

    /// A boolean indicating if inner products involving this determinant are complex-symmetric.
    #[pyo3(get)]
    pub(crate) complex_symmetric: bool,

    /// The complex coefficients for the molecular orbitals of this determinant.
    pub(crate) coefficients: Vec<Array2<C128>>,

    /// The occupation patterns for the molecular orbitals.
    pub(crate) occupations: Vec<Array1<f64>>,

    /// The threshold for comparisons.
    #[pyo3(get)]
    pub(crate) threshold: f64,

    /// The optional complex molecular orbital energies.
    pub(crate) mo_energies: Option<Vec<Array1<C128>>>,

    /// The optional complex determinantal energy.
    #[pyo3(get)]
    pub(crate) energy: Option<C128>,
}

#[pymethods]
impl PySlaterDeterminantComplex {
    /// Constructs a complex Python-exposed Slater determinant.
    ///
    /// # Arguments
    ///
    /// * `structure_constraint` - The structure constraint applied to the coefficients of the
    /// determinant.
    /// * `complex_symmetric` - A boolean indicating if inner products involving this determinant
    /// are complex-symmetric.
    /// * `coefficients` - The complex coefficients for the molecular orbitals of this determinant.
    /// * `occupations` - The occupation patterns for the molecular orbitals.
    /// * `threshold` - The threshold for comparisons.
    /// * `mo_energies` - The optional complex molecular orbital energies.
    /// * `energy` - The optional complex determinantal energy.
    #[new]
    #[pyo3(signature = (structure_constraint, complex_symmetric, coefficients, occupations, threshold, mo_energies=None, energy=None))]
    pub fn new(
        structure_constraint: PyStructureConstraint,
        complex_symmetric: bool,
        coefficients: Vec<Bound<'_, PyArray2<C128>>>,
        occupations: Vec<Bound<'_, PyArray1<f64>>>,
        threshold: f64,
        mo_energies: Option<Vec<Bound<'_, PyArray1<C128>>>>,
        energy: Option<C128>,
    ) -> Self {
        Self {
            structure_constraint,
            complex_symmetric,
            coefficients: coefficients
                .iter()
                .map(|pyarr| pyarr.to_owned_array())
                .collect::<Vec<_>>(),
            occupations: occupations
                .iter()
                .map(|pyarr| pyarr.to_owned_array())
                .collect::<Vec<_>>(),
            threshold,
            mo_energies: mo_energies.map(|energies| {
                energies
                    .iter()
                    .map(|pyarr| pyarr.to_owned_array())
                    .collect::<Vec<_>>()
            }),
            energy,
        }
    }

    /// The occupation patterns for the molecular orbitals.
    #[getter]
    pub fn occupations<'py>(&self, py: Python<'py>) -> PyResult<Vec<Bound<'py, PyArray1<f64>>>> {
        Ok(self
            .occupations
            .iter()
            .map(|occ| occ.to_pyarray(py))
            .collect::<Vec<_>>())
    }

    /// The complex coefficients for the molecular orbitals of this determinant.
    #[getter]
    pub fn coefficients<'py>(&self, py: Python<'py>) -> PyResult<Vec<Bound<'py, PyArray2<C128>>>> {
        Ok(self
            .coefficients
            .iter()
            .map(|occ| occ.to_pyarray(py))
            .collect::<Vec<_>>())
    }

    /// The complex molecular orbital energies, if any.
    #[getter]
    pub fn mo_energies<'py>(
        &self,
        py: Python<'py>,
    ) -> PyResult<Option<Vec<Bound<'py, PyArray1<C128>>>>> {
        Ok(self.mo_energies.as_ref().map(|mo_energies| {
            mo_energies
                .iter()
                .map(|mo_en| mo_en.to_pyarray(py))
                .collect::<Vec<_>>()
        }))
    }
}

impl PySlaterDeterminantComplex {
    /// Extracts the information in the [`PySlaterDeterminantComplex`] structure into `QSym2`'s native
    /// [`SlaterDeterminant`] structure.
    ///
    /// # Arguments
    ///
    /// * `baos` - The [`BasisAngularOrder`]s for the basis set in which the Slater determinant is
    ///   given, one for each explicit component per coefficient matrix.
    /// * `mol` - The molecule with which the Slater determinant is associated.
    ///
    /// # Returns
    ///
    /// The [`SlaterDeterminant`] structure with the same information.
    ///
    /// # Errors
    ///
    /// Errors if the [`SlaterDeterminant`] fails to build.
    pub fn to_qsym2<'b, 'a: 'b, SC>(
        &'b self,
        baos: &[&'a BasisAngularOrder],
        mol: &'a Molecule,
    ) -> Result<SlaterDeterminant<'b, C128, SC>, anyhow::Error>
    where
        SC: StructureConstraint
            + Clone
            + fmt::Display
            + TryFrom<PyStructureConstraint, Error = anyhow::Error>,
    {
        SlaterDeterminant::<C128, SC>::builder()
            .structure_constraint(self.structure_constraint.clone().try_into()?)
            .baos(baos.to_vec())
            .complex_symmetric(self.complex_symmetric)
            .mol(mol)
            .coefficients(&self.coefficients)
            .occupations(&self.occupations)
            .mo_energies(self.mo_energies.clone())
            .energy(
                self.energy
                    .ok_or_else(|| "No determinantal energy set.".to_string()),
            )
            .threshold(self.threshold)
            .build()
            .map_err(|err| format_err!(err))
    }

    pub fn structure_constraint(&self) -> &PyStructureConstraint {
        &self.structure_constraint
    }
}

impl<'a, SC> SlaterDeterminant<'a, C128, SC>
where
    SC: StructureConstraint
        + Clone
        + fmt::Display
        + TryInto<PyStructureConstraint, Error = anyhow::Error>,
{
    /// Extracts the information in the complex [`SlaterDeterminant`] structure into [`PySlaterDeterminantComplex`].
    ///
    /// # Returns
    ///
    /// The [`PySlaterDeterminantComplex`] structure with the same information.
    pub fn to_python<'py>(
        &self,
        py: Python<'py>,
    ) -> Result<PySlaterDeterminantComplex, anyhow::Error> {
        Ok(PySlaterDeterminantComplex::new(
            self.structure_constraint().clone().try_into()?,
            self.complex_symmetric(),
            self.coefficients()
                .iter()
                .map(|coeffs| coeffs.to_pyarray(py))
                .collect_vec(),
            self.occupations()
                .iter()
                .map(|occ| occ.to_pyarray(py))
                .collect_vec(),
            self.threshold(),
            self.mo_energies().map(|mo_energies| {
                mo_energies
                    .iter()
                    .map(|mo_e| mo_e.to_pyarray(py))
                    .collect_vec()
            }),
            self.energy().ok().copied(),
        ))
    }
}

// --------------------------------------------
// Slater determinant symmetry analysis results
// --------------------------------------------

/// Python-exposed structure storing the results of Slater determinant representation analysis.
#[pyclass]
#[derive(Clone)]
pub struct PySlaterDeterminantRepAnalysisResult {
    /// The group used for the representation analysis.
    #[pyo3(get)]
    group: String,

    /// The deduced overall symmetry of the determinant.
    #[pyo3(get)]
    determinant_symmetry: Option<String>,

    /// The deduced symmetries of the molecular orbitals constituting the determinant, if required.
    #[pyo3(get)]
    mo_symmetries: Option<Vec<Vec<Option<String>>>>,

    /// The deduced symmetries of the various densities constructible from the determinant, if
    /// required. In each tuple, the first element gives a description of the density corresponding
    /// to the symmetry result.
    #[pyo3(get)]
    determinant_density_symmetries: Option<Vec<(String, Option<String>)>>,

    /// The deduced symmetries of the total densities of the molecular orbitals constituting the
    /// determinant, if required.
    #[pyo3(get)]
    mo_density_symmetries: Option<Vec<Vec<Option<String>>>>,
}

// ================
// Enum definitions
// ================

/// Python-exposed enumerated type to handle the union type
/// `PySlaterDeterminantReal | PySlaterDeterminantComplex` in Python.
#[derive(FromPyObject)]
pub enum PySlaterDeterminant {
    /// Variant for real Python-exposed Slater determinant.
    Real(PySlaterDeterminantReal),

    /// Variant for complex Python-exposed Slater determinant.
    Complex(PySlaterDeterminantComplex),
}

// =====================
// Functions definitions
// =====================

/// Python-exposed function to perform representation symmetry analysis for real and complex
/// Slater determinants and log the result via the `qsym2-output` logger at the `INFO` level.
///
/// If `symmetry_transformation_kind` includes spin transformation, the provided determinant will
/// be augmented to generalised spin constraint automatically.
///
/// # Arguments
///
/// * `inp_sym` - A path to the [`QSym2FileType::Sym`] file containing the symmetry-group detection
/// result for the system. This will be used to construct abstract groups and character tables for
/// representation analysis.
/// * `pydet` - A Python-exposed Slater determinant whose coefficients are of type `float64` or
/// `complex128`.
/// * `pybaos` - Python-exposed structures containing basis angular order information, one for each
/// explicit component per coefficient matrix.
/// * `integrality_threshold` - The threshold for verifying if subspace multiplicities are
/// integral.
/// * `linear_independence_threshold` - The threshold for determining the linear independence
/// subspace via the non-zero eigenvalues of the orbit overlap matrix.
/// * `orbit_overlap_symmetry_threshold` - The threshold for checking the symmetry of the orbit
/// overlap matrix.
/// * `use_magnetic_group` - An option indicating if the magnetic group is to be used for symmetry
/// analysis, and if so, whether unitary representations or unitary-antiunitary corepresentations
/// should be used.
/// * `use_double_group` - A boolean indicating if the double group of the prevailing symmetry
/// group is to be used for representation analysis instead.
/// * `use_cayley_table` - A boolean indicating if the Cayley table for the group, if available,
/// should be used to speed up the calculation of orbit overlap matrices.
/// * `symmetry_transformation_kind` - An enumerated type indicating the type of symmetry
/// transformations to be performed on the origin determinant to generate the orbit. If this
/// contains spin transformation, the determinant will be augmented to generalised spin constraint
/// automatically.
/// * `eigenvalue_comparison_mode` - An enumerated type indicating the mode of comparison of orbit
/// overlap eigenvalues with the specified `linear_independence_threshold`.
/// * `sao` - The atomic-orbital overlap matrix whose elements are of type `float64` or
/// `complex128`.
/// * `sao_h` - The optional complex-symmetric atomic-orbital overlap matrix whose elements
/// are of type `float64` or `complex128`. This is required if antiunitary symmetry operations are
/// involved.
/// * `sao_spatial_4c` - The optional atomic-orbital four-centre overlap matrix whose elements are
/// of type `float64` or `complex128`.
/// * `sao_spatial_4c_h` - The optional complex-symmetric atomic-orbital four-centre overlap matrix
/// whose elements are of type `float64` or `complex128`. This is required if antiunitary symmetry
/// operations are involved.
/// * `analyse_mo_symmetries` - A boolean indicating if the symmetries of individual molecular
/// orbitals are to be analysed.
/// * `analyse_mo_symmetry_projections` - A boolean indicating if the symmetry projections of
/// individual molecular orbitals are to be analysed.
/// * `analyse_mo_mirror_parities` - A boolean indicating if the mirror parities of individual
/// molecular orbitals are to be printed.
/// * `analyse_density_symmetries` - A boolean indicating if the symmetries of densities are to be
/// analysed.
/// * `write_overlap_eigenvalues` - A boolean indicating if the eigenvalues of the determinant
/// orbit overlap matrix are to be written to the output.
/// * `write_character_table` - A boolean indicating if the character table of the prevailing
/// symmetry group is to be printed out.
/// * `infinite_order_to_finite` - The finite order with which infinite-order generators are to be
/// interpreted to form a finite subgroup of the prevailing infinite group. This finite subgroup
/// will be used for symmetry analysis.
/// * `angular_function_integrality_threshold` - The threshold for verifying if subspace
/// multiplicities are integral for the symmetry analysis of angular functions.
/// * `angular_function_linear_independence_threshold` - The threshold for determining the linear
/// independence subspace via the non-zero eigenvalues of the orbit overlap matrix for the symmetry
/// analysis of angular functions.
/// * `angular_function_orbit_overlap_symmetry_threshold` - The threshold for checking the symmetry
/// of the orbit overlap matrices for the symmetry analysis of angular functions.
/// * `angular_function_max_angular_momentum` - The maximum angular momentum order to be used in
/// angular function symmetry analysis.
///
/// # Returns
///
/// A Python-exposed [`PySlaterDeterminantRepAnalysisResult`] structure containing the results of the
/// representation analysis.
#[allow(clippy::too_many_arguments)]
#[pyfunction]
#[pyo3(signature = (
    inp_sym,
    pydet,
    pybaos,
    integrality_threshold,
    linear_independence_threshold,
    orbit_overlap_symmetry_threshold,
    use_magnetic_group,
    use_double_group,
    use_cayley_table,
    symmetry_transformation_kind,
    eigenvalue_comparison_mode,
    sao,
    sao_h=None,
    sao_spatial_4c=None,
    sao_spatial_4c_h=None,
    analyse_mo_symmetries=true,
    analyse_mo_symmetry_projections=true,
    analyse_mo_mirror_parities=false,
    analyse_density_symmetries=false,
    write_overlap_eigenvalues=true,
    write_character_table=true,
    infinite_order_to_finite=None,
    angular_function_integrality_threshold=1e-7,
    angular_function_linear_independence_threshold=1e-7,
    angular_function_orbit_overlap_symmetry_threshold=1e-7,
    angular_function_max_angular_momentum=2
))]
pub fn rep_analyse_slater_determinant(
    py: Python<'_>,
    inp_sym: PathBuf,
    pydet: PySlaterDeterminant,
    pybaos: Vec<PyBasisAngularOrder>,
    integrality_threshold: f64,
    linear_independence_threshold: f64,
    orbit_overlap_symmetry_threshold: f64,
    use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
    use_double_group: bool,
    use_cayley_table: bool,
    symmetry_transformation_kind: SymmetryTransformationKind,
    eigenvalue_comparison_mode: EigenvalueComparisonMode,
    sao: PyArray2RC,
    sao_h: Option<PyArray2RC>,
    sao_spatial_4c: Option<PyArray4RC>,
    sao_spatial_4c_h: Option<PyArray4RC>,
    analyse_mo_symmetries: bool,
    analyse_mo_symmetry_projections: bool,
    analyse_mo_mirror_parities: bool,
    analyse_density_symmetries: bool,
    write_overlap_eigenvalues: bool,
    write_character_table: bool,
    infinite_order_to_finite: Option<u32>,
    angular_function_integrality_threshold: f64,
    angular_function_linear_independence_threshold: f64,
    angular_function_orbit_overlap_symmetry_threshold: f64,
    angular_function_max_angular_momentum: u32,
) -> PyResult<PySlaterDeterminantRepAnalysisResult> {
    let pd_res: SymmetryGroupDetectionResult =
        read_qsym2_binary(inp_sym.clone(), QSym2FileType::Sym)
            .map_err(|err| PyIOError::new_err(err.to_string()))?;

    let mut file_name = inp_sym.to_path_buf();
    file_name.set_extension(QSym2FileType::Sym.ext());
    qsym2_output!(
        "Symmetry-group detection results read in from {}.",
        file_name.display(),
    );
    qsym2_output!("");

    let mol = &pd_res.pre_symmetry.recentred_molecule;

    let baos = pybaos
        .iter()
        .map(|bao| {
            bao.to_qsym2(mol)
                .map_err(|err| PyRuntimeError::new_err(err.to_string()))
        })
        .collect::<Result<Vec<_>, _>>()?;
    let baos_ref = baos.iter().collect::<Vec<_>>();
    let augment_to_generalised = match symmetry_transformation_kind {
        SymmetryTransformationKind::SpatialWithSpinTimeReversal
        | SymmetryTransformationKind::Spin
        | SymmetryTransformationKind::SpinSpatial => true,
        SymmetryTransformationKind::Spatial => false,
    };
    let afa_params = AngularFunctionRepAnalysisParams::builder()
        .integrality_threshold(angular_function_integrality_threshold)
        .linear_independence_threshold(angular_function_linear_independence_threshold)
        .orbit_overlap_symmetry_threshold(angular_function_orbit_overlap_symmetry_threshold)
        .max_angular_momentum(angular_function_max_angular_momentum)
        .build()
        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
    let sda_params = SlaterDeterminantRepAnalysisParams::<f64>::builder()
        .integrality_threshold(integrality_threshold)
        .linear_independence_threshold(linear_independence_threshold)
        .orbit_overlap_symmetry_threshold(orbit_overlap_symmetry_threshold)
        .use_magnetic_group(use_magnetic_group.clone())
        .use_double_group(use_double_group)
        .use_cayley_table(use_cayley_table)
        .symmetry_transformation_kind(symmetry_transformation_kind)
        .eigenvalue_comparison_mode(eigenvalue_comparison_mode)
        .analyse_mo_symmetries(analyse_mo_symmetries)
        .analyse_mo_symmetry_projections(analyse_mo_symmetry_projections)
        .analyse_mo_mirror_parities(analyse_mo_mirror_parities)
        .analyse_density_symmetries(analyse_density_symmetries)
        .write_overlap_eigenvalues(write_overlap_eigenvalues)
        .write_character_table(if write_character_table {
            Some(CharacterTableDisplay::Symbolic)
        } else {
            None
        })
        .infinite_order_to_finite(infinite_order_to_finite)
        .build()
        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
    let pysda_res: PySlaterDeterminantRepAnalysisResult = match (&pydet, &sao) {
        (PySlaterDeterminant::Real(pydet_r), PyArray2RC::Real(pysao_r)) => {
            if matches!(
                pydet_r.structure_constraint,
                PyStructureConstraint::SpinOrbitCoupled(_)
            ) {
                return Err(PyRuntimeError::new_err("Real determinants are not compatible with spin--orbit-coupled structure constraint.".to_string()));
            }

            let sao = pysao_r.to_owned_array();
            let sao_spatial_4c = sao_spatial_4c.and_then(|pysao4c| match pysao4c {
                // sao_spatial_4c must have the same reality as sao_spatial.
                PyArray4RC::Real(pysao4c_r) => Some(pysao4c_r.to_owned_array()),
                PyArray4RC::Complex(_) => None,
            });

            let det_r = if augment_to_generalised {
                pydet_r
                    .to_qsym2::<SpinConstraint>(&baos_ref, mol)
                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
                    .to_generalised()
            } else {
                pydet_r
                    .to_qsym2::<SpinConstraint>(&baos_ref, mol)
                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
            };
            match &use_magnetic_group {
                Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
                    let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                        MagneticRepresentedSymmetryGroup,
                        f64,
                        SpinConstraint,
                    >::builder()
                    .parameters(&sda_params)
                    .angular_function_parameters(&afa_params)
                    .determinant(&det_r)
                    .sao(&sao)
                    .sao_h(None) // Real SAO.
                    .sao_spatial_4c(sao_spatial_4c.as_ref())
                    .sao_spatial_4c_h(None) // Real SAO.
                    .symmetry_group(&pd_res)
                    .build()
                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    py.detach(|| {
                        sda_driver
                            .run()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                    })?;
                    let sda_res = sda_driver
                        .result()
                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    PySlaterDeterminantRepAnalysisResult {
                        group: sda_res.group().name().clone(),
                        determinant_symmetry: sda_res
                            .determinant_symmetry()
                            .as_ref()
                            .ok()
                            .map(|sym| sym.to_string()),
                        mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                            mo_symss
                                .iter()
                                .map(|mo_syms| {
                                    mo_syms
                                        .iter()
                                        .map(|mo_sym_opt| {
                                            mo_sym_opt.as_ref().map(|mo_sym| mo_sym.to_string())
                                        })
                                        .collect::<Vec<_>>()
                                })
                                .collect::<Vec<_>>()
                        }),
                        determinant_density_symmetries: sda_res
                            .determinant_density_symmetries()
                            .as_ref()
                            .map(|den_syms| {
                                den_syms
                                    .iter()
                                    .map(|(den_name, den_sym_res)| {
                                        (
                                            den_name.clone(),
                                            den_sym_res
                                                .as_ref()
                                                .ok()
                                                .map(|den_sym| den_sym.to_string()),
                                        )
                                    })
                                    .collect::<Vec<_>>()
                            }),
                        mo_density_symmetries: sda_res.mo_density_symmetries().as_ref().map(
                            |mo_den_symss| {
                                mo_den_symss
                                    .iter()
                                    .map(|mo_den_syms| {
                                        mo_den_syms
                                            .iter()
                                            .map(|mo_den_sym_opt| {
                                                mo_den_sym_opt
                                                    .as_ref()
                                                    .map(|mo_den_sym| mo_den_sym.to_string())
                                            })
                                            .collect::<Vec<_>>()
                                    })
                                    .collect::<Vec<_>>()
                            },
                        ),
                    }
                }
                Some(MagneticSymmetryAnalysisKind::Representation) | None => {
                    let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                        UnitaryRepresentedSymmetryGroup,
                        f64,
                        SpinConstraint,
                    >::builder()
                    .parameters(&sda_params)
                    .angular_function_parameters(&afa_params)
                    .determinant(&det_r)
                    .sao(&sao)
                    .sao_h(None) // Real SAO.
                    .sao_spatial_4c(sao_spatial_4c.as_ref())
                    .sao_spatial_4c_h(None) // Real SAO.
                    .symmetry_group(&pd_res)
                    .build()
                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    py.detach(|| {
                        sda_driver
                            .run()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                    })?;
                    let sda_res = sda_driver
                        .result()
                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    PySlaterDeterminantRepAnalysisResult {
                        group: sda_res.group().name().clone(),
                        determinant_symmetry: sda_res
                            .determinant_symmetry()
                            .as_ref()
                            .ok()
                            .map(|sym| sym.to_string()),
                        mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                            mo_symss
                                .iter()
                                .map(|mo_syms| {
                                    mo_syms
                                        .iter()
                                        .map(|mo_sym_opt| {
                                            mo_sym_opt.as_ref().map(|mo_sym| mo_sym.to_string())
                                        })
                                        .collect::<Vec<_>>()
                                })
                                .collect::<Vec<_>>()
                        }),
                        determinant_density_symmetries: sda_res
                            .determinant_density_symmetries()
                            .as_ref()
                            .map(|den_syms| {
                                den_syms
                                    .iter()
                                    .map(|(den_name, den_sym_res)| {
                                        (
                                            den_name.clone(),
                                            den_sym_res
                                                .as_ref()
                                                .ok()
                                                .map(|den_sym| den_sym.to_string()),
                                        )
                                    })
                                    .collect::<Vec<_>>()
                            }),
                        mo_density_symmetries: sda_res.mo_density_symmetries().as_ref().map(
                            |mo_den_symss| {
                                mo_den_symss
                                    .iter()
                                    .map(|mo_den_syms| {
                                        mo_den_syms
                                            .iter()
                                            .map(|mo_den_sym_opt| {
                                                mo_den_sym_opt
                                                    .as_ref()
                                                    .map(|mo_den_sym| mo_den_sym.to_string())
                                            })
                                            .collect::<Vec<_>>()
                                    })
                                    .collect::<Vec<_>>()
                            },
                        ),
                    }
                }
            }
        }
        (PySlaterDeterminant::Real(pydet_r), PyArray2RC::Complex(pysao_c)) => {
            match pydet_r.structure_constraint {
                PyStructureConstraint::SpinConstraint(_) => {
                    let det_r = if augment_to_generalised {
                        pydet_r
                            .to_qsym2(&baos_ref, mol)
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
                            .to_generalised()
                    } else {
                        pydet_r
                            .to_qsym2(&baos_ref, mol)
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
                    };
                    let det_c: SlaterDeterminant<C128, SpinConstraint> = det_r.into();
                    let sao_c = pysao_c.to_owned_array();
                    let sao_h_c = sao_h.and_then(|pysao_h| match pysao_h {
                        // sao_spatial_h must have the same reality as sao_spatial.
                        PyArray2RC::Real(_) => None,
                        PyArray2RC::Complex(pysao_h_c) => Some(pysao_h_c.to_owned_array()),
                    });
                    let sao_spatial_4c_c = sao_spatial_4c.and_then(|pysao4c| match pysao4c {
                        // sao_spatial_4c must have the same reality as sao_spatial.
                        PyArray4RC::Real(_) => None,
                        PyArray4RC::Complex(pysao4c_c) => Some(pysao4c_c.to_owned_array()),
                    });
                    let sao_spatial_4c_h_c =
                        sao_spatial_4c_h.and_then(|pysao4c_h| match pysao4c_h {
                            // sao_spatial_4c_h must have the same reality as sao_spatial.
                            PyArray4RC::Real(_) => None,
                            PyArray4RC::Complex(pysao4c_h_c) => Some(pysao4c_h_c.to_owned_array()),
                        });
                    match &use_magnetic_group {
                        Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                MagneticRepresentedSymmetryGroup,
                                C128,
                                SpinConstraint,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                        Some(MagneticSymmetryAnalysisKind::Representation) | None => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                UnitaryRepresentedSymmetryGroup,
                                C128,
                                SpinConstraint,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                    }
                }
                PyStructureConstraint::SpinOrbitCoupled(_) => {
                    let det_r = pydet_r
                        .to_qsym2::<SpinOrbitCoupled>(&baos_ref, mol)
                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    let det_c: SlaterDeterminant<C128, SpinOrbitCoupled> = det_r.into();
                    let sao_c = pysao_c.to_owned_array();
                    let sao_h_c = sao_h.and_then(|pysao_h| match pysao_h {
                        // sao_h must have the same reality as sao.
                        PyArray2RC::Real(_) => None,
                        PyArray2RC::Complex(pysao_h_c) => Some(pysao_h_c.to_owned_array()),
                    });
                    let sao_spatial_4c_c = sao_spatial_4c.and_then(|pysao4c| match pysao4c {
                        // sao_spatial_4c must have the same reality as sao.
                        PyArray4RC::Real(_) => None,
                        PyArray4RC::Complex(pysao4c_c) => Some(pysao4c_c.to_owned_array()),
                    });
                    let sao_spatial_4c_h_c =
                        sao_spatial_4c_h.and_then(|pysao4c_h| match pysao4c_h {
                            // sao_spatial_4c_h must have the same reality as sao_spatial.
                            PyArray4RC::Real(_) => None,
                            PyArray4RC::Complex(pysao4c_h_c) => Some(pysao4c_h_c.to_owned_array()),
                        });
                    match &use_magnetic_group {
                        Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                MagneticRepresentedSymmetryGroup,
                                C128,
                                SpinOrbitCoupled,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                        Some(MagneticSymmetryAnalysisKind::Representation) | None => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                UnitaryRepresentedSymmetryGroup,
                                C128,
                                SpinOrbitCoupled,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                    }
                }
            }
        }
        (PySlaterDeterminant::Complex(pydet_c), _) => {
            match pydet_c.structure_constraint {
                PyStructureConstraint::SpinConstraint(_) => {
                    let det_c = if augment_to_generalised {
                        pydet_c
                            .to_qsym2::<SpinConstraint>(&baos_ref, mol)
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
                            .to_generalised()
                    } else {
                        pydet_c
                            .to_qsym2::<SpinConstraint>(&baos_ref, mol)
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
                    };
                    let sao_c = match sao {
                        PyArray2RC::Real(pysao_r) => pysao_r.to_owned_array().mapv(Complex::from),
                        PyArray2RC::Complex(pysao_c) => pysao_c.to_owned_array(),
                    };
                    let sao_h_c = sao_h.map(|pysao_h| match pysao_h {
                        // sao_h must have the same reality as sao.
                        PyArray2RC::Real(pysao_h_r) => {
                            pysao_h_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray2RC::Complex(pysao_h_c) => pysao_h_c.to_owned_array(),
                    });
                    let sao_spatial_4c_c = sao_spatial_4c.map(|pysao4c| match pysao4c {
                        // sao_spatial_4c must have the same reality as sao.
                        PyArray4RC::Real(pysao4c_r) => {
                            pysao4c_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray4RC::Complex(pysao4c_c) => pysao4c_c.to_owned_array(),
                    });
                    let sao_spatial_4c_h_c = sao_spatial_4c_h.map(|pysao4c_h| match pysao4c_h {
                        // sao_spatial_4c_h must have the same reality as sao.
                        PyArray4RC::Real(pysao4c_h_r) => {
                            pysao4c_h_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray4RC::Complex(pysao4c_h_c) => pysao4c_h_c.to_owned_array(),
                    });
                    match &use_magnetic_group {
                        Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                MagneticRepresentedSymmetryGroup,
                                C128,
                                SpinConstraint,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                        Some(MagneticSymmetryAnalysisKind::Representation) | None => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                UnitaryRepresentedSymmetryGroup,
                                C128,
                                SpinConstraint,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                    }
                }
                PyStructureConstraint::SpinOrbitCoupled(_) => {
                    let det_c = pydet_c
                        .to_qsym2::<SpinOrbitCoupled>(&baos_ref, mol)
                        .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                    let sao_c = match sao {
                        PyArray2RC::Real(pysao_r) => pysao_r.to_owned_array().mapv(Complex::from),
                        PyArray2RC::Complex(pysao_c) => pysao_c.to_owned_array(),
                    };
                    let sao_h_c = sao_h.map(|pysao_h| match pysao_h {
                        // sao_spatial_h must have the same reality as sao.
                        PyArray2RC::Real(pysao_h_r) => {
                            pysao_h_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray2RC::Complex(pysao_h_c) => pysao_h_c.to_owned_array(),
                    });
                    let sao_spatial_4c_c = sao_spatial_4c.map(|pysao4c| match pysao4c {
                        // sao_spatial_4c must have the same reality as sao.
                        PyArray4RC::Real(pysao4c_r) => {
                            pysao4c_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray4RC::Complex(pysao4c_c) => pysao4c_c.to_owned_array(),
                    });
                    let sao_spatial_4c_h_c = sao_spatial_4c_h.map(|pysao4c_h| match pysao4c_h {
                        // sao_spatial_4c_h must have the same reality as sao.
                        PyArray4RC::Real(pysao4c_h_r) => {
                            pysao4c_h_r.to_owned_array().mapv(Complex::from)
                        }
                        PyArray4RC::Complex(pysao4c_h_c) => pysao4c_h_c.to_owned_array(),
                    });
                    match &use_magnetic_group {
                        Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                MagneticRepresentedSymmetryGroup,
                                C128,
                                SpinOrbitCoupled,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                        Some(MagneticSymmetryAnalysisKind::Representation) | None => {
                            let mut sda_driver = SlaterDeterminantRepAnalysisDriver::<
                                UnitaryRepresentedSymmetryGroup,
                                C128,
                                SpinOrbitCoupled,
                            >::builder()
                            .parameters(&sda_params)
                            .angular_function_parameters(&afa_params)
                            .determinant(&det_c)
                            .sao(&sao_c)
                            .sao_h(sao_h_c.as_ref())
                            .sao_spatial_4c(sao_spatial_4c_c.as_ref())
                            .sao_spatial_4c_h(sao_spatial_4c_h_c.as_ref())
                            .symmetry_group(&pd_res)
                            .build()
                            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            py.detach(|| {
                                sda_driver
                                    .run()
                                    .map_err(|err| PyRuntimeError::new_err(err.to_string()))
                            })?;
                            let sda_res = sda_driver
                                .result()
                                .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
                            PySlaterDeterminantRepAnalysisResult {
                                group: sda_res.group().name().clone(),
                                determinant_symmetry: sda_res
                                    .determinant_symmetry()
                                    .as_ref()
                                    .ok()
                                    .map(|sym| sym.to_string()),
                                mo_symmetries: sda_res.mo_symmetries().as_ref().map(|mo_symss| {
                                    mo_symss
                                        .iter()
                                        .map(|mo_syms| {
                                            mo_syms
                                                .iter()
                                                .map(|mo_sym_opt| {
                                                    mo_sym_opt
                                                        .as_ref()
                                                        .map(|mo_sym| mo_sym.to_string())
                                                })
                                                .collect::<Vec<_>>()
                                        })
                                        .collect::<Vec<_>>()
                                }),
                                determinant_density_symmetries: sda_res
                                    .determinant_density_symmetries()
                                    .as_ref()
                                    .map(|den_syms| {
                                        den_syms
                                            .iter()
                                            .map(|(den_name, den_sym_res)| {
                                                (
                                                    den_name.clone(),
                                                    den_sym_res
                                                        .as_ref()
                                                        .ok()
                                                        .map(|den_sym| den_sym.to_string()),
                                                )
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                                mo_density_symmetries: sda_res
                                    .mo_density_symmetries()
                                    .as_ref()
                                    .map(|mo_den_symss| {
                                        mo_den_symss
                                            .iter()
                                            .map(|mo_den_syms| {
                                                mo_den_syms
                                                    .iter()
                                                    .map(|mo_den_sym_opt| {
                                                        mo_den_sym_opt.as_ref().map(|mo_den_sym| {
                                                            mo_den_sym.to_string()
                                                        })
                                                    })
                                                    .collect::<Vec<_>>()
                                            })
                                            .collect::<Vec<_>>()
                                    }),
                            }
                        }
                    }
                }
            }
        }
    };
    Ok(pysda_res)
}