qsym2/bindings/python/
molecule_symmetrisation.rs

1//! Python bindings for QSym² molecule symmetrisation by bootstrapping.
2//!
3//! See [`crate::drivers::molecule_symmetrisation_bootstrap`] for more information.
4
5use std::path::PathBuf;
6
7use anyhow::format_err;
8use pyo3::exceptions::PyRuntimeError;
9use pyo3::prelude::*;
10
11use crate::auxiliary::atom::AtomKind;
12use crate::auxiliary::molecule::Molecule;
13use crate::bindings::python::symmetry_group_detection::PyMolecule;
14use crate::drivers::molecule_symmetrisation_bootstrap::{
15    MoleculeSymmetrisationBootstrapDriver, MoleculeSymmetrisationBootstrapParams,
16};
17use crate::drivers::QSym2Driver;
18use crate::io::QSym2FileType;
19
20/// Python-exposed function to perform molecule symmetrisation by bootstrapping and log the result
21/// via the `qsym2-output` logger at the `INFO` level.
22///
23/// See [`crate::drivers::molecule_symmetrisation_bootstrap`] for more information.
24///
25/// # Arguments
26///
27/// * `inp_xyz` - An optional string providing the path to an XYZ file containing the molecule to
28/// be symmetrised. Only one of `inp_xyz` or `inp_mol` can be specified. Python type: `Optional[str]`.
29/// * `inp_mol` - An optional `PyMolecule` structure containing the molecule to be symmetrised. Only
30/// one of `inp_xyz` or `inp_mol` can be specified. Python type: `PyMolecule`.
31/// * `out_target_sym` - An optional path for a [`QSym2FileType::Sym`] file to be saved that
32/// contains the symmetry-group detection results of the symmetrised molecule at the target
33/// thresholds. Python type: `Optional[str]`.
34/// * `loose_moi_threshold` - The loose MoI threshold. Python type: `float`.
35/// * `loose_distance_threshold` - The loose distance threshold. Python type: `float`.
36/// * `target_moi_threshold` - The target (tight) MoI threshold. Python type: `float`.
37/// * `target_distance_threshold` - The target (tight) distance threshold. Python type: `float`.
38/// * `use_magnetic_group` - A boolean indicating if the magnetic group (*i.e.* the group including
39/// time-reversed operations) is to be used for the symmetrisation. Python type: `bool`.
40/// * `reorientate_molecule` - A boolean indicating if the molecule is also reoriented to align its
41/// principal axes with the Cartesian axes. Python type: `bool`.
42/// * `max_iterations` - The maximum number of iterations for the symmetrisation process. Python
43/// type: `int`.
44/// * `consistent_target_symmetry_iterations` - The number of consecutive iterations during which
45/// the symmetry group at the target level of threshold must be consistently found for convergence
46/// to be reached, if this group cannot become identical to the symmetry group at the loose level
47/// of threshold. Python type: `int`.
48/// * `verbose` - The print-out level. Python type: `int`.
49/// * `infinite_order_to_finite` - The finite order with which infinite-order generators are to be
50/// interpreted to form a finite subgroup of the prevailing infinite group. This finite subgroup
51/// will be used for the symmetrisation. Python type: `Optional[int]`.
52///
53/// # Returns
54///
55/// The symmetrised molecule.
56///
57/// Python type: `PyMolecule`
58///
59/// # Errors
60///
61/// Errors if any intermediate step in the symmetrisation procedure fails.
62#[pyfunction]
63#[pyo3(signature = (
64    inp_xyz,
65    inp_mol,
66    out_target_sym=None,
67    loose_moi_threshold=1e-2,
68    loose_distance_threshold=1e-2,
69    target_moi_threshold=1e-7,
70    target_distance_threshold=1e-7,
71    use_magnetic_group=false,
72    reorientate_molecule=true,
73    max_iterations=50,
74    consistent_target_symmetry_iterations=10,
75    verbose=0,
76    infinite_order_to_finite=None
77))]
78pub fn symmetrise_molecule(
79    py: Python<'_>,
80    inp_xyz: Option<PathBuf>,
81    inp_mol: Option<PyMolecule>,
82    out_target_sym: Option<PathBuf>,
83    loose_moi_threshold: f64,
84    loose_distance_threshold: f64,
85    target_moi_threshold: f64,
86    target_distance_threshold: f64,
87    use_magnetic_group: bool,
88    reorientate_molecule: bool,
89    max_iterations: usize,
90    consistent_target_symmetry_iterations: usize,
91    verbose: u8,
92    infinite_order_to_finite: Option<u32>,
93) -> PyResult<PyMolecule> {
94    py.allow_threads(|| {
95        let mol = match (inp_xyz, inp_mol) {
96            (Some(xyz_path), None) => Molecule::from_xyz(xyz_path, 1e-7),
97            (None, Some(pymol)) => Molecule::from(pymol),
98            _ => {
99                return Err(PyRuntimeError::new_err(
100                    "One and only one of `inp_xyz` or `inp_mol` must be specified.",
101                ))
102            }
103        };
104
105        let msb_params = MoleculeSymmetrisationBootstrapParams::builder()
106            .reorientate_molecule(reorientate_molecule)
107            .use_magnetic_group(use_magnetic_group)
108            .loose_moi_threshold(loose_moi_threshold)
109            .loose_distance_threshold(loose_distance_threshold)
110            .target_moi_threshold(target_moi_threshold)
111            .target_distance_threshold(target_distance_threshold)
112            .infinite_order_to_finite(infinite_order_to_finite)
113            .max_iterations(max_iterations)
114            .consistent_target_symmetry_iterations(consistent_target_symmetry_iterations)
115            .verbose(verbose)
116            .symmetrised_result_save_name(out_target_sym)
117            .build()
118            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
119
120        let mut msb_driver = MoleculeSymmetrisationBootstrapDriver::builder()
121            .parameters(&msb_params)
122            .molecule(&mol)
123            .build()
124            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
125        msb_driver
126            .run()
127            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
128
129        let symmol = &msb_driver
130            .result()
131            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?
132            .symmetrised_molecule;
133
134        // Note that the magnetic field here will not have the original magnitude as it is back-deduded
135        // from the magnetic atoms which have been added using the normalised version of the original
136        // magnetic field.
137        let magnetic_field = symmol
138            .magnetic_atoms
139            .as_ref()
140            .map(|mag_atoms| {
141                if mag_atoms.len() != 2 {
142                    Err(format_err!("Only a uniform magnetic field is supported."))
143                } else {
144                    match (&mag_atoms[0].kind, &mag_atoms[1].kind) {
145                        (AtomKind::Magnetic(true), AtomKind::Magnetic(false)) => {
146                            let bvec = mag_atoms[0].coordinates - mag_atoms[1].coordinates;
147                            Ok([bvec[0], bvec[1], bvec[2]])
148                        }
149                        (AtomKind::Magnetic(false), AtomKind::Magnetic(true)) => {
150                            let bvec = mag_atoms[1].coordinates - mag_atoms[0].coordinates;
151                            Ok([bvec[0], bvec[1], bvec[2]])
152                        }
153                        _ => Err(format_err!("Invalid fictitious magnetic atoms detected.")),
154                    }
155                }
156            })
157            .transpose()
158            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
159
160        // Note that the electric field here will not have the original magnitude as it is back-deduded
161        // from the electric atoms which have been added using the normalised version of the original
162        // electric field.
163        let electric_field = symmol
164            .electric_atoms
165            .as_ref()
166            .map(|elec_atoms| {
167                if elec_atoms.len() != 1 {
168                    Err(format_err!("Only a uniform electric field is supported."))
169                } else {
170                    match &elec_atoms[0].kind {
171                        AtomKind::Electric(pos) => {
172                            let evec = if *pos {
173                                elec_atoms[0].coordinates
174                            } else {
175                                -elec_atoms[0].coordinates
176                            };
177                            Ok([evec[0], evec[1], evec[2]])
178                        }
179                        _ => Err(format_err!("Invalid fictitious electric atoms detected.")),
180                    }
181                }
182            })
183            .transpose()
184            .map_err(|err| PyRuntimeError::new_err(err.to_string()))?;
185
186        let pymol = PyMolecule::new(
187            symmol
188                .atoms
189                .iter()
190                .map(|atom| {
191                    (
192                        atom.atomic_symbol.clone(),
193                        [
194                            atom.coordinates[0],
195                            atom.coordinates[1],
196                            atom.coordinates[2],
197                        ],
198                    )
199                })
200                .collect::<Vec<_>>(),
201            symmol.threshold,
202            magnetic_field,
203            electric_field,
204        );
205        Ok(pymol)
206    })
207}