qsym2/drivers/projection/
density.rs

1//! Driver for symmetry projection of electron densities.
2
3use std::collections::HashSet;
4use std::fmt;
5use std::str::FromStr;
6
7use anyhow::{bail, format_err};
8use derive_builder::Builder;
9use duplicate::duplicate_item;
10use indexmap::IndexMap;
11use ndarray_linalg::{Lapack, Norm};
12use num::Complex;
13use num_complex::ComplexFloat;
14use serde::{Deserialize, Serialize, de::DeserializeOwned};
15
16use crate::analysis::EigenvalueComparisonMode;
17use crate::chartab::CharacterTable;
18use crate::chartab::chartab_group::CharacterProperties;
19use crate::drivers::QSym2Driver;
20use crate::drivers::representation_analysis::{
21    CharacterTableDisplay, MagneticSymmetryAnalysisKind, fn_construct_unitary_group, log_bao,
22    log_cc_transversal,
23};
24use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
25use crate::group::{GroupProperties, UnitaryRepresentedGroup};
26use crate::io::format::{
27    QSym2Output, log_subtitle, nice_bool, qsym2_output, write_subtitle, write_title,
28};
29use crate::projection::Projectable;
30use crate::symmetry::symmetry_group::{SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup};
31use crate::symmetry::symmetry_symbols::MullikenIrrepSymbol;
32use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
33use crate::target::density::Density;
34use crate::target::density::density_analysis::DensitySymmetryOrbit;
35
36#[cfg(test)]
37#[path = "density_tests.rs"]
38mod density_tests;
39
40// ----------
41// Parameters
42// ----------
43
44const fn default_symbolic() -> Option<CharacterTableDisplay> {
45    Some(CharacterTableDisplay::Symbolic)
46}
47
48/// Structure containing control parameters for electron density projection.
49#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
50pub struct DensityProjectionParams {
51    /// Option indicating if the magnetic group is to be used for projection, and if so,
52    /// whether unitary representations or unitary-antiunitary corepresentations should be used.
53    #[builder(default = "None")]
54    #[serde(default)]
55    pub use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
56
57    /// Boolean indicating if the double group is to be used for symmetry projection.
58    #[builder(default = "false")]
59    #[serde(default)]
60    pub use_double_group: bool,
61
62    /// Option indicating if the character table of the group used for symmetry projection is to be
63    /// printed out.
64    #[builder(default = "Some(CharacterTableDisplay::Symbolic)")]
65    #[serde(default = "default_symbolic")]
66    pub write_character_table: Option<CharacterTableDisplay>,
67
68    /// The kind of symmetry transformation to be applied on the reference density to generate
69    /// the orbit for symmetry projection.
70    #[builder(default = "SymmetryTransformationKind::Spatial")]
71    #[serde(default)]
72    pub symmetry_transformation_kind: SymmetryTransformationKind,
73
74    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
75    /// subgroup of an infinite group can be used for the symmetry projection.
76    #[builder(default = "None")]
77    #[serde(default)]
78    pub infinite_order_to_finite: Option<u32>,
79
80    /// The projection targets supplied symbolically.
81    #[builder(default = "None")]
82    #[serde(default)]
83    pub symbolic_projection_targets: Option<Vec<String>>,
84
85    /// The projection targets supplied numerically, where each value gives the index of the
86    /// projection subspace based on the group's character table.
87    #[builder(default = "None")]
88    #[serde(default)]
89    pub numeric_projection_targets: Option<Vec<usize>>,
90}
91
92impl DensityProjectionParams {
93    /// Returns a builder to construct a [`DensityProjectionParams`] structure.
94    pub fn builder() -> DensityProjectionParamsBuilder {
95        DensityProjectionParamsBuilder::default()
96    }
97}
98
99impl fmt::Display for DensityProjectionParams {
100    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
101        if let Some(symbolic_projection_targets) = self.symbolic_projection_targets.as_ref() {
102            writeln!(f, "Projection subspaces (symbolic):")?;
103            for projection_target in symbolic_projection_targets.iter() {
104                writeln!(f, "  {projection_target}")?;
105            }
106            writeln!(f)?;
107        }
108        if let Some(numeric_projection_targets) = self.numeric_projection_targets.as_ref() {
109            writeln!(f, "Projection subspaces (numeric):")?;
110            for projection_target in numeric_projection_targets.iter() {
111                writeln!(f, "  {projection_target}")?;
112            }
113            writeln!(f)?;
114        }
115        writeln!(
116            f,
117            "Use magnetic group for projection: {}",
118            match self.use_magnetic_group {
119                None => "no",
120                Some(MagneticSymmetryAnalysisKind::Representation) =>
121                    "yes, using unitary representations",
122                Some(MagneticSymmetryAnalysisKind::Corepresentation) =>
123                    "yes, using magnetic corepresentations",
124            }
125        )?;
126        writeln!(
127            f,
128            "Use double group for projection: {}",
129            nice_bool(self.use_double_group)
130        )?;
131        if let Some(finite_order) = self.infinite_order_to_finite {
132            writeln!(f, "Infinite order to finite: {finite_order}")?;
133        }
134        writeln!(
135            f,
136            "Symmetry transformation kind: {}",
137            self.symmetry_transformation_kind
138        )?;
139        writeln!(f)?;
140        writeln!(
141            f,
142            "Write character table: {}",
143            if let Some(chartab_display) = self.write_character_table.as_ref() {
144                format!("yes, {}", chartab_display.to_string().to_lowercase())
145            } else {
146                "no".to_string()
147            }
148        )?;
149
150        Ok(())
151    }
152}
153
154// ------
155// Result
156// ------
157
158/// Structure to contain electron density projection results.
159#[derive(Clone, Builder)]
160pub struct DensityProjectionResult<'a, G, T>
161where
162    G: SymmetryGroupProperties + Clone + 'a,
163    T: ComplexFloat + Lapack,
164    Density<'a, T>: SymmetryTransformable,
165    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
166{
167    /// The control parameters used to obtain this set of electron density projection results.
168    parameters: &'a DensityProjectionParams,
169
170    /// The densities being projected and their associated names or descriptions.
171    densities: Vec<(String, &'a Density<'a, T>)>,
172
173    /// The group used for the projection.
174    group: G,
175
176    /// The projected densities. Each tuple in the vector contains the name or description of the
177    /// density being projected and an indexmap containing the projected densities indexed by the
178    /// requested subspace labels.
179    #[allow(clippy::type_complexity)]
180    projected_densities: Vec<(
181        String,
182        IndexMap<G::RowSymbol, Result<Density<'a, T>, String>>,
183    )>,
184}
185
186impl<'a, G, T> DensityProjectionResult<'a, G, T>
187where
188    G: SymmetryGroupProperties + Clone,
189    G::RowSymbol: Serialize + DeserializeOwned,
190    T: ComplexFloat + Lapack,
191    Density<'a, T>: SymmetryTransformable,
192    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
193{
194    pub fn builder() -> DensityProjectionResultBuilder<'a, G, T> {
195        DensityProjectionResultBuilder::default()
196    }
197
198    /// Returns the control parameters used to obtain this set of electron density projection results.
199    pub fn parameters(&self) -> &DensityProjectionParams {
200        self.parameters
201    }
202
203    /// Returns the densities being projected and their associated names or descriptions.
204    pub fn densities(&self) -> &Vec<(String, &'a Density<'a, T>)> {
205        &self.densities
206    }
207
208    /// Returns the projected densities.
209    #[allow(clippy::type_complexity)]
210    pub fn projected_densities(
211        &self,
212    ) -> &Vec<(
213        String,
214        IndexMap<G::RowSymbol, Result<Density<'a, T>, String>>,
215    )> {
216        &self.projected_densities
217    }
218}
219
220impl<'a, G, T> fmt::Display for DensityProjectionResult<'a, G, T>
221where
222    G: SymmetryGroupProperties + Clone,
223    G::RowSymbol: Serialize + DeserializeOwned,
224    T: ComplexFloat + Lapack,
225    Density<'a, T>: SymmetryTransformable,
226    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
227{
228    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
229        write_subtitle(f, "Orbit-based symmetry projection summary")?;
230        writeln!(f)?;
231        writeln!(
232            f,
233            "> Group: {} ({})",
234            self.group
235                .finite_subgroup_name()
236                .map(|subgroup_name| format!("{} > {}", self.group.name(), subgroup_name))
237                .unwrap_or(self.group.name()),
238            self.group.group_type().to_string().to_lowercase()
239        )?;
240        writeln!(f)?;
241
242        for (name, projected_densities) in self.projected_densities.iter() {
243            writeln!(f, ">> Density-matrix Frobenius norm for projected {name}")?;
244
245            let (rows, norms): (Vec<String>, Vec<String>) = projected_densities
246                .iter()
247                .map(|(row, den_res)| {
248                    let norm = den_res
249                        .as_ref()
250                        .map(|den| format!("{:.3e}", den.density_matrix().norm_l2()))
251                        .unwrap_or_else(|err| err.to_string());
252                    (row.to_string(), norm)
253                })
254                .unzip();
255
256            let row_length = rows
257                .iter()
258                .map(|row| row.chars().count())
259                .max()
260                .unwrap_or(8)
261                .max(8);
262            let norm_length = norms
263                .iter()
264                .map(|norm| norm.chars().count())
265                .max()
266                .unwrap_or(14)
267                .max(14);
268            let table_width = 4 + row_length + norm_length;
269            writeln!(f, "{}", "┈".repeat(table_width))?;
270            writeln!(f, " {:<row_length$}  Frobenius norm", "Subspace",)?;
271            writeln!(f, "{}", "┈".repeat(table_width))?;
272            for (row, norm) in rows.iter().zip(norms) {
273                writeln!(f, " {:<row_length$}  {:<}", row, norm)?;
274            }
275            writeln!(f, "{}", "┈".repeat(table_width))?;
276
277            writeln!(f)?;
278        }
279        Ok(())
280    }
281}
282
283impl<'a, G, T> fmt::Debug for DensityProjectionResult<'a, G, T>
284where
285    G: SymmetryGroupProperties + Clone,
286    G::RowSymbol: Serialize + DeserializeOwned,
287    T: ComplexFloat + Lapack,
288    Density<'a, T>: SymmetryTransformable,
289    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
290{
291    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
292        writeln!(f, "{self}")
293    }
294}
295
296// ------
297// Driver
298// ------
299
300// ~~~~~~~~~~~~~~~~~
301// Struct definition
302// ~~~~~~~~~~~~~~~~~
303//
304/// Driver structure for performing projection on electron densities.
305#[derive(Clone, Builder)]
306#[builder(build_fn(validate = "Self::validate"))]
307pub struct DensityProjectionDriver<'a, G, T>
308where
309    G: SymmetryGroupProperties + Clone,
310    G::RowSymbol: Serialize + DeserializeOwned,
311    T: ComplexFloat + Lapack,
312    Density<'a, T>: SymmetryTransformable,
313    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
314{
315    /// The control parameters used to obtain this set of electron density projection results.
316    parameters: &'a DensityProjectionParams,
317
318    /// The densities being projected and their associated names or descriptions.
319    densities: Vec<(String, &'a Density<'a, T>)>,
320
321    /// The result from symmetry-group detection on the underlying molecular structure of the
322    /// electron densities. Only the unitary symmetry group will be used for projection, since
323    /// magnetic-group projection is not yet formulated.
324    symmetry_group: &'a SymmetryGroupDetectionResult,
325
326    /// The result of the electron density projection.
327    #[builder(setter(skip), default = "None")]
328    result: Option<DensityProjectionResult<'a, G, T>>,
329}
330
331impl<'a, G, T> DensityProjectionDriverBuilder<'a, G, T>
332where
333    G: SymmetryGroupProperties + Clone,
334    G::RowSymbol: Serialize + DeserializeOwned,
335    T: ComplexFloat + Lapack,
336    Density<'a, T>: SymmetryTransformable,
337    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
338{
339    fn validate(&self) -> Result<(), String> {
340        let params = self
341            .parameters
342            .ok_or("No electron density projection parameters found.".to_string())?;
343
344        let sym_res = self
345            .symmetry_group
346            .ok_or("No symmetry group information found.".to_string())?;
347
348        let dens = self
349            .densities
350            .as_ref()
351            .ok_or("No electron densities found.".to_string())?;
352
353        let sym = if params.use_magnetic_group.is_some() {
354            sym_res
355                .magnetic_symmetry
356                .as_ref()
357                .ok_or("Magnetic symmetry requested for symmetry projection, but no magnetic symmetry found.")?
358        } else {
359            &sym_res.unitary_symmetry
360        };
361
362        if sym.is_infinite() && params.infinite_order_to_finite.is_none() {
363            Err(format!(
364                "Projection cannot be performed using the entirety of the infinite group `{}`. \
365                    Consider setting the parameter `infinite_order_to_finite` to restrict to a finite subgroup instead.",
366                sym.group_name
367                    .as_ref()
368                    .expect("No symmetry group name found.")
369            ))
370        } else {
371            let baos = dens
372                .iter()
373                .map(|(_, den)| den.bao())
374                .collect::<HashSet<_>>();
375            if baos.len() != 1 {
376                Err("Inconsistent basis angular order information between densities.".to_string())
377            } else {
378                Ok(())
379            }
380        }
381    }
382}
383
384// ~~~~~~~~~~~~~~~~~~~~~~
385// Struct implementations
386// ~~~~~~~~~~~~~~~~~~~~~~
387
388// Generic for all symmetry groups G and density numeric type T
389// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
390
391impl<'a, G, T> DensityProjectionDriver<'a, G, T>
392where
393    G: SymmetryGroupProperties + Clone,
394    G::RowSymbol: Serialize + DeserializeOwned,
395    T: ComplexFloat + Lapack,
396    Density<'a, T>: SymmetryTransformable,
397    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
398{
399    /// Returns a builder to construct a [`DensityProjectionDriver`] structure.
400    pub fn builder() -> DensityProjectionDriverBuilder<'a, G, T> {
401        DensityProjectionDriverBuilder::default()
402    }
403}
404
405// Specific for unitary-represented symmetry groups, but generic for density numeric type T
406// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
407
408impl<'a, T> DensityProjectionDriver<'a, UnitaryRepresentedSymmetryGroup, T>
409where
410    T: ComplexFloat + Lapack,
411    Density<'a, T>: SymmetryTransformable,
412    DensitySymmetryOrbit<'a, UnitaryRepresentedSymmetryGroup, T>:
413        Projectable<UnitaryRepresentedSymmetryGroup, Density<'a, T>>,
414{
415    fn_construct_unitary_group!(
416        /// Constructs the unitary-represented group (which itself can be unitary or magnetic) ready
417        /// for electron density projection.
418        construct_unitary_group
419    );
420}
421
422// Specific for unitary-represented symmetry groups and density numeric types f64 and C128
423// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
424
425#[duplicate_item(
426    duplicate!{
427        [ dtype_nested; [f64]; [Complex<f64>] ]
428        [
429            gtype_ [ UnitaryRepresentedSymmetryGroup ]
430            dtype_ [ dtype_nested ]
431            doc_sub_ [ "Performs projection using a unitary-represented group and stores the result." ]
432            projection_fn_ [ projection_representation ]
433            construct_group_ [ self.construct_unitary_group()? ]
434        ]
435    }
436)]
437impl<'a> DensityProjectionDriver<'a, gtype_, dtype_> {
438    #[doc = doc_sub_]
439    fn projection_fn_(&mut self) -> Result<(), anyhow::Error> {
440        let params = self.parameters;
441        let group = construct_group_;
442        log_cc_transversal(&group);
443        let bao = self
444            .densities
445            .first()
446            .map(|(_, den)| den.bao())
447            .ok_or_else(|| {
448                format_err!("Basis angular order information could not be extracted.")
449            })?;
450        log_bao(bao, None);
451
452        let all_rows = group.character_table().get_all_rows();
453        let rows = params
454            .symbolic_projection_targets
455            .as_ref()
456            .unwrap_or(&vec![])
457            .iter()
458            .map(|row_str| MullikenIrrepSymbol::from_str(row_str).map_err(|err| format_err!(err)))
459            .chain(
460                params
461                    .numeric_projection_targets
462                    .as_ref()
463                    .unwrap_or(&vec![])
464                    .iter()
465                    .map(|row_index| {
466                        all_rows.get_index(*row_index).cloned().ok_or_else(|| {
467                            format_err!(
468                                "Unable to retrieve the subspace label with index {row_index}."
469                            )
470                        })
471                    }),
472            )
473            .collect::<Result<Vec<_>, _>>()?;
474
475        let projected_densities = self
476            .densities
477            .iter()
478            .map(|(name, den)| {
479                let projections = DensitySymmetryOrbit::builder()
480                    .group(&group)
481                    .origin(den)
482                    .symmetry_transformation_kind(params.symmetry_transformation_kind.clone())
483                    .integrality_threshold(1e-14)
484                    .linear_independence_threshold(1e-14)
485                    .eigenvalue_comparison_mode(EigenvalueComparisonMode::Modulus)
486                    .build()
487                    .map_err(|err| format_err!(err))
488                    .map(|den_orbit| {
489                        rows.iter()
490                            .map(|row| {
491                                (
492                                    row.clone(),
493                                    den_orbit
494                                        .project_onto(row)
495                                        .map_err(|err| err.to_string())
496                                        .and_then(|temp_projected_den| {
497                                            Density::builder()
498                                                .bao(bao)
499                                                .complex_symmetric(den.complex_symmetric())
500                                                .complex_conjugated(den.complex_conjugated())
501                                                .mol(den.mol())
502                                                .density_matrix(
503                                                    temp_projected_den.density_matrix().clone(),
504                                                )
505                                                .threshold(den.threshold())
506                                                .build()
507                                                .map_err(|err| err.to_string())
508                                        }),
509                                )
510                            })
511                            .collect::<IndexMap<_, _>>()
512                    })?;
513                Ok((name.clone(), projections))
514            })
515            .collect::<Result<Vec<_>, anyhow::Error>>()?;
516
517        let result = DensityProjectionResult::builder()
518            .parameters(params)
519            .densities(self.densities.clone())
520            .group(group.clone())
521            .projected_densities(projected_densities)
522            .build()?;
523        self.result = Some(result);
524
525        Ok(())
526    }
527}
528
529// ~~~~~~~~~~~~~~~~~~~~~
530// Trait implementations
531// ~~~~~~~~~~~~~~~~~~~~~
532
533// Generic for all symmetry groups G and density numeric type T
534// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
535
536impl<'a, G, T> fmt::Display for DensityProjectionDriver<'a, G, T>
537where
538    G: SymmetryGroupProperties + Clone,
539    G::RowSymbol: Serialize + DeserializeOwned,
540    T: ComplexFloat + Lapack,
541    Density<'a, T>: SymmetryTransformable,
542    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
543{
544    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
545        write_title(f, "Electron Density Symmetry Projection")?;
546        writeln!(f)?;
547        writeln!(f, "{}", self.parameters)?;
548        Ok(())
549    }
550}
551
552impl<'a, G, T> fmt::Debug for DensityProjectionDriver<'a, G, T>
553where
554    G: SymmetryGroupProperties + Clone,
555    G::RowSymbol: Serialize + DeserializeOwned,
556    T: ComplexFloat + Lapack,
557    Density<'a, T>: SymmetryTransformable,
558    DensitySymmetryOrbit<'a, G, T>: Projectable<G, Density<'a, T>>,
559{
560    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
561        writeln!(f, "{self}")
562    }
563}
564
565// Specific for unitary-represented symmetry groups and density numeric types f64 and C128
566// '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
567
568#[duplicate_item(
569    duplicate!{
570        [ dtype_nested; [f64]; [Complex<f64>] ]
571        [
572            gtype_ [ UnitaryRepresentedSymmetryGroup ]
573            dtype_ [ dtype_nested ]
574            doc_sub_ [ "Performs projection using a unitary-represented group and stores the result." ]
575            projection_fn_ [ projection_representation ]
576            construct_group_ [ self.construct_unitary_group()? ]
577        ]
578    }
579)]
580impl<'a> QSym2Driver for DensityProjectionDriver<'a, gtype_, dtype_> {
581    type Params = DensityProjectionParams;
582
583    type Outcome = DensityProjectionResult<'a, gtype_, dtype_>;
584
585    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
586        self.result
587            .as_ref()
588            .ok_or_else(|| format_err!("No electron density projection results found."))
589    }
590
591    fn run(&mut self) -> Result<(), anyhow::Error> {
592        self.log_output_display();
593        self.projection_fn_()?;
594        self.result()?.log_output_display();
595        Ok(())
596    }
597}
qsym2/drivers/projection/density.rs

qsym2/drivers/projection/
density.rs
