1use std::collections::HashMap;
4use std::fmt;
5use std::ops::Mul;
6
7use anyhow::{self, bail, format_err};
8use approx::abs_diff_eq;
9use derive_builder::Builder;
10use duplicate::duplicate_item;
11use indexmap::IndexMap;
12use itertools::Itertools;
13use ndarray::{s, Array2, Array4};
14use ndarray_linalg::types::Lapack;
15use num::ToPrimitive;
16use num_complex::{Complex, ComplexFloat};
17use num_traits::real::Real;
18use num_traits::Float;
19use rayon::prelude::*;
20use serde::{Deserialize, Serialize};
21
22use crate::analysis::{
23 log_overlap_eigenvalues, EigenvalueComparisonMode, Orbit, Overlap, ProjectionDecomposition,
24 RepAnalysis,
25};
26use crate::angmom::spinor_rotation_3d::{SpinConstraint, SpinOrbitCoupled, StructureConstraint};
27use crate::chartab::chartab_group::CharacterProperties;
28use crate::chartab::chartab_symbols::ReducibleLinearSpaceSymbol;
29use crate::chartab::{CharacterTable, SubspaceDecomposable};
30use crate::drivers::representation_analysis::angular_function::{
31 find_angular_function_representation, find_spinor_function_representation,
32 AngularFunctionRepAnalysisParams,
33};
34use crate::drivers::representation_analysis::{
35 fn_construct_magnetic_group, fn_construct_unitary_group, log_bao, log_cc_transversal,
36 CharacterTableDisplay, MagneticSymmetryAnalysisKind,
37};
38use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
39use crate::drivers::QSym2Driver;
40use crate::group::{GroupProperties, MagneticRepresentedGroup, UnitaryRepresentedGroup};
41use crate::io::format::{
42 log_subtitle, nice_bool, qsym2_output, write_subtitle, write_title, QSym2Output,
43};
44use crate::symmetry::symmetry_element::symmetry_operation::{
45 SpecialSymmetryTransformation, SymmetryOperation,
46};
47use crate::symmetry::symmetry_group::{
48 MagneticRepresentedSymmetryGroup, SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup,
49};
50use crate::symmetry::symmetry_symbols::{
51 deduce_mirror_parities, MirrorParity, MullikenIrcorepSymbol, SymmetryClassSymbol,
52};
53use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
54use crate::target::density::density_analysis::DensitySymmetryOrbit;
55use crate::target::determinant::determinant_analysis::SlaterDeterminantSymmetryOrbit;
56use crate::target::determinant::SlaterDeterminant;
57use crate::target::orbital::orbital_analysis::generate_det_mo_orbits;
58
59#[cfg(test)]
60#[path = "slater_determinant_tests.rs"]
61mod slater_determinant_tests;
62
63const fn default_true() -> bool {
72 true
73}
74const fn default_symbolic() -> Option<CharacterTableDisplay> {
75 Some(CharacterTableDisplay::Symbolic)
76}
77
78#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
80pub struct SlaterDeterminantRepAnalysisParams<T: From<f64>> {
81 pub integrality_threshold: T,
83
84 pub linear_independence_threshold: T,
86
87 #[builder(default = "true")]
90 #[serde(default = "default_true")]
91 pub analyse_mo_symmetries: bool,
92
93 #[builder(default = "false")]
96 #[serde(default)]
97 pub analyse_mo_mirror_parities: bool,
98
99 #[builder(default = "false")]
101 #[serde(default)]
102 pub analyse_mo_symmetry_projections: bool,
103
104 #[builder(default = "false")]
107 #[serde(default)]
108 pub analyse_density_symmetries: bool,
109
110 #[builder(default = "None")]
113 #[serde(default)]
114 pub use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
115
116 #[builder(default = "false")]
118 #[serde(default)]
119 pub use_double_group: bool,
120
121 #[builder(default = "true")]
124 #[serde(default = "default_true")]
125 pub use_cayley_table: bool,
126
127 #[builder(default = "SymmetryTransformationKind::Spatial")]
130 #[serde(default)]
131 pub symmetry_transformation_kind: SymmetryTransformationKind,
132
133 #[builder(default = "Some(CharacterTableDisplay::Symbolic)")]
136 #[serde(default = "default_symbolic")]
137 pub write_character_table: Option<CharacterTableDisplay>,
138
139 #[builder(default = "true")]
141 #[serde(default = "default_true")]
142 pub write_overlap_eigenvalues: bool,
143
144 #[builder(default = "EigenvalueComparisonMode::Modulus")]
146 #[serde(default)]
147 pub eigenvalue_comparison_mode: EigenvalueComparisonMode,
148
149 #[builder(default = "None")]
152 #[serde(default)]
153 pub infinite_order_to_finite: Option<u32>,
154}
155
156impl<T> SlaterDeterminantRepAnalysisParams<T>
157where
158 T: Float + From<f64>,
159{
160 pub fn builder() -> SlaterDeterminantRepAnalysisParamsBuilder<T> {
162 SlaterDeterminantRepAnalysisParamsBuilder::default()
163 }
164}
165
166impl Default for SlaterDeterminantRepAnalysisParams<f64> {
167 fn default() -> Self {
168 Self::builder()
169 .integrality_threshold(1e-7)
170 .linear_independence_threshold(1e-7)
171 .build()
172 .expect("Unable to construct a default `SlaterDeterminantRepAnalysisParams<f64>`.")
173 }
174}
175
176impl<T> fmt::Display for SlaterDeterminantRepAnalysisParams<T>
177where
178 T: From<f64> + fmt::LowerExp + fmt::Debug,
179{
180 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
181 writeln!(
182 f,
183 "Integrality threshold: {:.3e}",
184 self.integrality_threshold
185 )?;
186 writeln!(
187 f,
188 "Linear independence threshold: {:.3e}",
189 self.linear_independence_threshold
190 )?;
191 writeln!(
192 f,
193 "Orbit eigenvalue comparison mode: {}",
194 self.eigenvalue_comparison_mode
195 )?;
196 writeln!(
197 f,
198 "Write overlap eigenvalues: {}",
199 nice_bool(self.write_overlap_eigenvalues)
200 )?;
201 writeln!(f)?;
202 writeln!(
203 f,
204 "Analyse molecular orbital symmetry: {}",
205 nice_bool(self.analyse_mo_symmetries)
206 )?;
207 writeln!(
208 f,
209 "Analyse molecular orbital symmetry projection: {}",
210 nice_bool(self.analyse_mo_symmetry_projections)
211 )?;
212 writeln!(
213 f,
214 "Analyse molecular orbital mirror parity: {}",
215 nice_bool(self.analyse_mo_mirror_parities)
216 )?;
217 writeln!(
218 f,
219 "Analyse density symmetry: {}",
220 nice_bool(self.analyse_density_symmetries)
221 )?;
222 writeln!(f)?;
223 writeln!(
224 f,
225 "Use magnetic group for analysis: {}",
226 match self.use_magnetic_group {
227 None => "no",
228 Some(MagneticSymmetryAnalysisKind::Representation) =>
229 "yes, using unitary representations",
230 Some(MagneticSymmetryAnalysisKind::Corepresentation) =>
231 "yes, using magnetic corepresentations",
232 }
233 )?;
234 writeln!(
235 f,
236 "Use double group for analysis: {}",
237 nice_bool(self.use_double_group)
238 )?;
239 writeln!(
240 f,
241 "Use Cayley table for orbit overlap matrices: {}",
242 nice_bool(self.use_cayley_table)
243 )?;
244 if let Some(finite_order) = self.infinite_order_to_finite {
245 writeln!(f, "Infinite order to finite: {finite_order}")?;
246 }
247 writeln!(
248 f,
249 "Symmetry transformation kind: {}",
250 self.symmetry_transformation_kind
251 )?;
252 writeln!(f)?;
253 writeln!(
254 f,
255 "Write character table: {}",
256 if let Some(chartab_display) = self.write_character_table.as_ref() {
257 format!("yes, {}", chartab_display.to_string().to_lowercase())
258 } else {
259 "no".to_string()
260 }
261 )?;
262
263 Ok(())
264 }
265}
266
267#[derive(Clone, Builder)]
273pub struct SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
274where
275 G: SymmetryGroupProperties + Clone,
276 G::CharTab: SubspaceDecomposable<T>,
277 T: ComplexFloat + Lapack,
278 SC: StructureConstraint + fmt::Display,
279 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
280{
281 parameters: &'a SlaterDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
284
285 determinant: &'a SlaterDeterminant<'a, T, SC>,
287
288 group: G,
290
291 determinant_symmetry: Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
293
294 mo_symmetries: Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>>,
296
297 mo_symmetry_projections: Option<
299 Vec<
300 Vec<
301 Option<
302 Vec<(
303 <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
304 Complex<f64>,
305 )>,
306 >,
307 >,
308 >,
309 >,
310
311 mo_mirror_parities:
313 Option<Vec<Vec<Option<IndexMap<SymmetryClassSymbol<SymmetryOperation>, MirrorParity>>>>>,
314
315 mo_symmetries_thresholds: Option<Vec<Vec<(Option<T>, Option<T>)>>>,
318
319 determinant_density_symmetries: Option<
323 Vec<(
324 String,
325 Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
326 )>,
327 >,
328
329 mo_density_symmetries:
332 Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>>,
333}
334
335impl<'a, G, T, SC> SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
336where
337 G: SymmetryGroupProperties + Clone,
338 G::CharTab: SubspaceDecomposable<T>,
339 T: ComplexFloat + Lapack,
340 SC: StructureConstraint + Clone + fmt::Display,
341 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
342{
343 fn builder() -> SlaterDeterminantRepAnalysisResultBuilder<'a, G, T, SC> {
346 SlaterDeterminantRepAnalysisResultBuilder::default()
347 }
348}
349
350impl<'a, G, T, SC> SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
351where
352 G: SymmetryGroupProperties + Clone,
353 G::CharTab: SubspaceDecomposable<T>,
354 T: ComplexFloat + Lapack,
355 SC: StructureConstraint + fmt::Display,
356 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
357{
358 pub fn group(&self) -> &G {
360 &self.group
361 }
362
363 pub fn determinant_symmetry(
365 &self,
366 ) -> &Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String> {
367 &self.determinant_symmetry
368 }
369
370 pub fn mo_symmetries(
372 &self,
373 ) -> &Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>> {
374 &self.mo_symmetries
375 }
376
377 pub fn mo_symmetry_projections(
379 &self,
380 ) -> &Option<
381 Vec<
382 Vec<
383 Option<
384 Vec<(
385 <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
386 Complex<f64>,
387 )>,
388 >,
389 >,
390 >,
391 > {
392 &self.mo_symmetry_projections
393 }
394
395 pub fn determinant_density_symmetries(
399 &self,
400 ) -> &Option<
401 Vec<(
402 String,
403 Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
404 )>,
405 > {
406 &self.determinant_density_symmetries
407 }
408
409 pub fn mo_density_symmetries(
412 &self,
413 ) -> &Option<Vec<Vec<Option<<G::CharTab as SubspaceDecomposable<T>>::Decomposition>>>> {
414 &self.mo_density_symmetries
415 }
416}
417
418impl<'a, G, T, SC> fmt::Display for SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
419where
420 G: SymmetryGroupProperties + Clone,
421 G::CharTab: SubspaceDecomposable<T>,
422 T: ComplexFloat + Lapack,
423 SC: StructureConstraint + fmt::Display,
424 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
425{
426 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
427 write_subtitle(f, "Orbit-based symmetry analysis results")?;
428 writeln!(f)?;
429 writeln!(
430 f,
431 "> Group: {} ({})",
432 self.group
433 .finite_subgroup_name()
434 .map(|subgroup_name| format!("{} > {}", self.group.name(), subgroup_name))
435 .unwrap_or(self.group.name()),
436 self.group.group_type().to_string().to_lowercase()
437 )?;
438 writeln!(f)?;
439 writeln!(f, "> Overall determinantal result")?;
440 writeln!(
441 f,
442 " Energy : {}",
443 self.determinant
444 .energy()
445 .map(|e| e.to_string())
446 .unwrap_or_else(|err| format!("-- ({err})"))
447 )?;
448 writeln!(
449 f,
450 " Symmetry: {}",
451 self.determinant_symmetry
452 .as_ref()
453 .map(|s| s.to_string())
454 .unwrap_or_else(|err| format!("-- ({err})"))
455 )?;
456 writeln!(f)?;
457
458 if let Some(den_syms) = self.determinant_density_symmetries.as_ref() {
459 writeln!(f, "> Overall determinantal density result")?;
460 let den_type_width = den_syms
461 .iter()
462 .map(|(den_type, _)| den_type.chars().count())
463 .max()
464 .unwrap_or(7)
465 .max(7);
466 for (den_type, den_sym_res) in den_syms.iter() {
467 writeln!(
468 f,
469 " {den_type:<den_type_width$}: {}",
470 den_sym_res
471 .as_ref()
472 .map(|e| e.to_string())
473 .unwrap_or_else(|err| format!("-- ({err})"))
474 )?;
475 }
476 writeln!(f)?;
477 }
478
479 if let Some(mo_symmetries) = self.mo_symmetries.as_ref() {
480 let mo_spin_index_length = 4;
481 let mo_index_length = mo_symmetries
482 .iter()
483 .map(|spin_mo_symmetries| spin_mo_symmetries.len())
484 .max()
485 .and_then(|max_mo_length| usize::try_from(max_mo_length.ilog10() + 2).ok())
486 .unwrap_or(4);
487 let mo_occ_length = 5;
488 let mo_symmetry_length = mo_symmetries
489 .iter()
490 .flat_map(|spin_mo_symmetries| {
491 spin_mo_symmetries.iter().map(|mo_sym| {
492 mo_sym
493 .as_ref()
494 .map(|sym| sym.to_string())
495 .unwrap_or("--".to_string())
496 .chars()
497 .count()
498 })
499 })
500 .max()
501 .unwrap_or(0)
502 .max(8);
503 let mo_energies_opt = self.determinant.mo_energies();
504 let mo_energy_length = mo_energies_opt
505 .map(|mo_energies| {
506 mo_energies
507 .iter()
508 .flat_map(|spin_mo_energies| {
509 spin_mo_energies.map(|v| format!("{v:+.7}").chars().count())
510 })
511 .max()
512 .unwrap_or(0)
513 })
514 .unwrap_or(0)
515 .max(6);
516
517 let mo_eig_above_length: usize = self
518 .mo_symmetries_thresholds
519 .as_ref()
520 .map(|mo_symmetries_thresholds| {
521 mo_symmetries_thresholds
522 .iter()
523 .flat_map(|spin_mo_symmetries_thresholds| {
524 spin_mo_symmetries_thresholds.iter().map(|(above, _)| {
525 above
526 .as_ref()
527 .map(|eig| format!("{eig:+.3e}"))
528 .unwrap_or("--".to_string())
529 .chars()
530 .count()
531 })
532 })
533 .max()
534 .unwrap_or(10)
535 .max(10)
536 })
537 .unwrap_or(10);
538 let mo_eig_below_length: usize = self
539 .mo_symmetries_thresholds
540 .as_ref()
541 .map(|mo_symmetries_thresholds| {
542 mo_symmetries_thresholds
543 .iter()
544 .flat_map(|spin_mo_symmetries_thresholds| {
545 spin_mo_symmetries_thresholds.iter().map(|(_, below)| {
546 below
547 .as_ref()
548 .map(|eig| format!("{eig:+.3e}"))
549 .unwrap_or("--".to_string())
550 .chars()
551 .count()
552 })
553 })
554 .max()
555 .unwrap_or(10)
556 .max(10)
557 })
558 .unwrap_or(10);
559
560 let precision = Real::ceil(ComplexFloat::abs(ComplexFloat::log10(
561 self.parameters.linear_independence_threshold,
562 )))
563 .to_usize()
564 .expect("Unable to convert the linear independence threshold exponent to `usize`.");
565 let mo_sym_projections_str_opt =
566 self.mo_symmetry_projections
567 .as_ref()
568 .map(|mo_sym_projectionss| {
569 mo_sym_projectionss
570 .iter()
571 .enumerate()
572 .map(|(ispin, mo_sym_projections)| {
573 mo_sym_projections
574 .iter()
575 .enumerate()
576 .map(|(imo, mo_sym_projection)| {
577 mo_sym_projection
578 .as_ref()
579 .map(|sym_proj| {
580 if let Some(mo_symmetriess) = self.mo_symmetries() {
581 if let Some(sym) = &mo_symmetriess[ispin][imo] {
582 let sym_proj_hashmap = sym_proj
583 .iter()
584 .cloned()
585 .collect::<HashMap<_, _>>();
586 sym.subspaces()
587 .iter()
588 .map(|(subspace, _)| {
589 format!(
590 "{subspace}: {}",
591 sym_proj_hashmap
592 .get(subspace)
593 .map(|composition| {
594 if abs_diff_eq!(
595 composition.im,
596 0.0,
597 epsilon = self.parameters.linear_independence_threshold.to_f64().expect("Unable to convert the linear independence threshold to `f64`.")
598 ) {
599 format!(
600 "{:+.precision$}",
601 composition.re
602 )
603 } else {
604 format!(
605 "{composition:+.precision$}")
606 }
607 })
608 .unwrap_or(
609 "--".to_string()
610 )
611 )
612 })
613 .join(", ")
614 } else {
615 "--".to_string()
616 }
617 } else {
618 "--".to_string()
619 }
620 })
621 .unwrap_or("--".to_string())
622 })
623 .collect::<Vec<_>>()
624 })
625 .collect::<Vec<_>>()
626 });
627 let mo_sym_projections_length_opt =
628 mo_sym_projections_str_opt
629 .as_ref()
630 .map(|mo_sym_projectionss| {
631 mo_sym_projectionss
632 .iter()
633 .flat_map(|mo_sym_projections| {
634 mo_sym_projections
635 .iter()
636 .map(|mo_sym_projection| mo_sym_projection.chars().count())
637 })
638 .max()
639 .unwrap_or(18)
640 .max(18)
641 });
642 let mo_sym_projections_length = mo_sym_projections_length_opt.unwrap_or(0);
643 let mo_sym_projections_gap = mo_sym_projections_length_opt.map(|_| 2).unwrap_or(0);
644 let mo_sym_projections_heading = mo_sym_projections_length_opt
645 .map(|_| "MO sym. projection")
646 .unwrap_or("");
647
648 let mirrors = self
649 .group
650 .filter_cc_symbols(|cc| cc.is_spatial_reflection());
651 let mo_mirror_parities_length_opt = self.mo_mirror_parities.as_ref().map(|_| {
652 let mirror_heading = mirrors.iter().map(|sigma| format!("p[{sigma}]")).join(" ");
653 let length = mirror_heading.chars().count();
654 (mirror_heading, length)
655 });
656 let mo_mirror_parities_gap = mo_mirror_parities_length_opt
657 .as_ref()
658 .map(|_| 2)
659 .unwrap_or(0);
660 let (mo_mirror_parities_heading, mo_mirror_parities_length) =
661 mo_mirror_parities_length_opt.unwrap_or((String::new(), 0));
662
663 let mo_den_symss_str_opt = self.mo_density_symmetries.as_ref().map(|mo_den_symss| {
664 mo_den_symss
665 .iter()
666 .map(|mo_den_syms| {
667 mo_den_syms
668 .iter()
669 .map(|mo_den_sym| {
670 mo_den_sym
671 .as_ref()
672 .map(|sym| sym.to_string())
673 .unwrap_or("--".to_string())
674 })
675 .collect::<Vec<_>>()
676 })
677 .collect::<Vec<_>>()
678 });
679 let mo_density_length_opt = mo_den_symss_str_opt.as_ref().map(|mo_den_symss| {
680 mo_den_symss
681 .iter()
682 .flat_map(|mo_den_syms| {
683 mo_den_syms
684 .iter()
685 .map(|mo_den_sym| mo_den_sym.chars().count())
686 })
687 .max()
688 .unwrap_or(12)
689 .max(12)
690 });
691 let mo_density_length = mo_density_length_opt.unwrap_or(0);
692 let mo_density_gap = mo_density_length_opt.map(|_| 2).unwrap_or(0);
693 let mo_density_heading = mo_density_length_opt.map(|_| "Density sym.").unwrap_or("");
694
695 let table_width = 14
696 + mo_spin_index_length
697 + mo_index_length
698 + mo_occ_length
699 + mo_energy_length
700 + mo_symmetry_length
701 + mo_mirror_parities_gap
702 + mo_mirror_parities_length
703 + mo_eig_above_length
704 + mo_eig_below_length
705 + mo_sym_projections_gap
706 + mo_sym_projections_length
707 + mo_density_gap
708 + mo_density_length;
709
710 writeln!(f, "> Molecular orbital results")?;
711 writeln!(
712 f,
713 " Structure constraint: {}",
714 self.determinant
715 .structure_constraint()
716 .to_string()
717 .to_lowercase()
718 )?;
719 if self.mo_mirror_parities.as_ref().is_some() {
720 writeln!(f, "")?;
721 writeln!(
722 f,
723 "Column p[σ] gives the parity under the reflection class σ: {} => even, {} => odd, {} => neither.",
724 MirrorParity::Even,
725 MirrorParity::Odd,
726 MirrorParity::Neither
727 )?;
728 }
729 writeln!(f, "{}", "┈".repeat(table_width))?;
730 if mo_density_length > 0 {
731 writeln!(
732 f,
733 " {:>mo_spin_index_length$} {:>mo_index_length$} {:<mo_occ_length$} {:<mo_energy_length$} {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$} {:<mo_eig_above_length$} {:<mo_eig_below_length$}{}{:mo_sym_projections_length$}{}{}",
734 "Spin",
735 "MO",
736 "Occ.",
737 "Energy",
738 "Symmetry",
739 " ".repeat(mo_mirror_parities_gap),
740 mo_mirror_parities_heading,
741 "Eig. above",
742 "Eig. below",
743 " ".repeat(mo_sym_projections_gap),
744 mo_sym_projections_heading,
745 " ".repeat(mo_density_gap),
746 mo_density_heading
747 )?;
748 } else if mo_sym_projections_length > 0 {
749 writeln!(
750 f,
751 " {:>mo_spin_index_length$} {:>mo_index_length$} {:<mo_occ_length$} {:<mo_energy_length$} {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$} {:<mo_eig_above_length$} {:<mo_eig_below_length$}{}{}",
752 "Spin",
753 "MO",
754 "Occ.",
755 "Energy",
756 "Symmetry",
757 " ".repeat(mo_mirror_parities_gap),
758 mo_mirror_parities_heading,
759 "Eig. above",
760 "Eig. below",
761 " ".repeat(mo_sym_projections_gap),
762 mo_sym_projections_heading,
763 )?;
764 } else {
765 writeln!(
766 f,
767 " {:>mo_spin_index_length$} {:>mo_index_length$} {:<mo_occ_length$} {:<mo_energy_length$} {:<mo_symmetry_length$}{}{:mo_mirror_parities_length$} {:<mo_eig_above_length$} {}",
768 "Spin",
769 "MO",
770 "Occ.",
771 "Energy",
772 "Symmetry",
773 " ".repeat(mo_mirror_parities_gap),
774 mo_mirror_parities_heading,
775 "Eig. above",
776 "Eig. below",
777 )?;
778 };
779 writeln!(f, "{}", "┈".repeat(table_width))?;
780
781 let empty_string = String::new();
782 for (spini, spin_mo_symmetries) in mo_symmetries.iter().enumerate() {
783 writeln!(f, " Spin {spini}")?;
784 for (moi, mo_sym) in spin_mo_symmetries.iter().enumerate() {
785 let occ_str = self
786 .determinant
787 .occupations()
788 .get(spini)
789 .and_then(|spin_occs| spin_occs.get(moi))
790 .map(|occ| format!("{occ:>.3}"))
791 .unwrap_or("--".to_string());
792 let mo_energy_str = mo_energies_opt
793 .and_then(|mo_energies| mo_energies.get(spini))
794 .and_then(|spin_mo_energies| spin_mo_energies.get(moi))
795 .map(|mo_energy| format!("{mo_energy:>+mo_energy_length$.7}"))
796 .unwrap_or("--".to_string());
797 let mo_sym_str = mo_sym
798 .as_ref()
799 .map(|sym| sym.to_string())
800 .unwrap_or("--".to_string());
801
802 let mo_mirror_parities_str = self
803 .mo_mirror_parities
804 .as_ref()
805 .and_then(|mo_mirror_paritiess| {
806 mo_mirror_paritiess
807 .get(spini)
808 .and_then(|spin_mo_mirror_parities| {
809 spin_mo_mirror_parities
810 .get(moi)
811 .map(|mo_mirror_parities_opt| {
812 mo_mirror_parities_opt
813 .as_ref()
814 .map(|mo_mirror_parities| {
815 mirrors
816 .iter()
817 .map(|sigma| {
818 let sigma_length =
819 sigma.to_string().chars().count()
820 + 3;
821 mo_mirror_parities
822 .get(sigma)
823 .map(|parity| {
824 format!(
825 "{:^sigma_length$}",
826 parity.to_string()
827 )
828 })
829 .unwrap_or_else(|| {
830 format!(
831 "{:^sigma_length$}",
832 "--"
833 )
834 })
835 })
836 .join(" ")
837 })
838 .unwrap_or(String::new())
839 })
840 })
841 })
842 .unwrap_or(String::new());
843
844 let (eig_above_str, eig_below_str) = self
845 .mo_symmetries_thresholds
846 .as_ref()
847 .map(|mo_symmetries_thresholds| {
848 mo_symmetries_thresholds
849 .get(spini)
850 .and_then(|spin_mo_symmetries_thresholds| {
851 spin_mo_symmetries_thresholds.get(moi)
852 })
853 .map(|(eig_above_opt, eig_below_opt)| {
854 (
855 eig_above_opt
856 .map(|eig_above| format!("{eig_above:>+.3e}"))
857 .unwrap_or("--".to_string()),
858 eig_below_opt
859 .map(|eig_below| format!("{eig_below:>+.3e}"))
860 .unwrap_or("--".to_string()),
861 )
862 })
863 .unwrap_or(("--".to_string(), "--".to_string()))
864 })
865 .unwrap_or(("--".to_string(), "--".to_string()));
866
867 let mo_symmetry_projections_str = mo_sym_projections_str_opt
868 .as_ref()
869 .and_then(|mo_symmetry_projectionss| {
870 mo_symmetry_projectionss.get(spini).and_then(
871 |spin_mo_symmetry_projections| {
872 spin_mo_symmetry_projections.get(moi)
873 },
874 )
875 })
876 .unwrap_or(&empty_string);
877
878 let mo_density_symmetries_str = mo_den_symss_str_opt
879 .as_ref()
880 .and_then(|mo_density_symmetriess| {
881 mo_density_symmetriess.get(spini).and_then(
882 |spin_mo_density_symmetries| spin_mo_density_symmetries.get(moi),
883 )
884 })
885 .unwrap_or(&empty_string);
886
887 match (mo_density_length, mo_sym_projections_length) {
888 (0, 0) => {
889 writeln!(
890 f,
891 " {spini:>mo_spin_index_length$} \
892 {moi:>mo_index_length$} \
893 {occ_str:<mo_occ_length$} \
894 {mo_energy_str:<mo_energy_length$} \
895 {mo_sym_str:<mo_symmetry_length$}\
896 {}{:mo_mirror_parities_length$} \
897 {eig_above_str:<mo_eig_above_length$} \
898 {eig_below_str}",
899 " ".repeat(mo_mirror_parities_gap),
900 mo_mirror_parities_str,
901 )?;
902 }
903 (_, 0) => {
904 writeln!(
905 f,
906 " {spini:>mo_spin_index_length$} \
907 {moi:>mo_index_length$} \
908 {occ_str:<mo_occ_length$} \
909 {mo_energy_str:<mo_energy_length$} \
910 {mo_sym_str:<mo_symmetry_length$}\
911 {}{:mo_mirror_parities_length$} \
912 {eig_above_str:<mo_eig_above_length$} \
913 {eig_below_str:<mo_eig_below_length$} \
914 {mo_density_symmetries_str}",
915 " ".repeat(mo_mirror_parities_gap),
916 mo_mirror_parities_str,
917 )?;
918 }
919 (0, _) => {
920 writeln!(
921 f,
922 " {spini:>mo_spin_index_length$} \
923 {moi:>mo_index_length$} \
924 {occ_str:<mo_occ_length$} \
925 {mo_energy_str:<mo_energy_length$} \
926 {mo_sym_str:<mo_symmetry_length$}\
927 {}{:mo_mirror_parities_length$} \
928 {eig_above_str:<mo_eig_above_length$} \
929 {eig_below_str:<mo_eig_below_length$} \
930 {mo_symmetry_projections_str}",
931 " ".repeat(mo_mirror_parities_gap),
932 mo_mirror_parities_str,
933 )?;
934 }
935 (_, _) => {
936 writeln!(
937 f,
938 " {spini:>mo_spin_index_length$} \
939 {moi:>mo_index_length$} \
940 {occ_str:<mo_occ_length$} \
941 {mo_energy_str:<mo_energy_length$} \
942 {mo_sym_str:<mo_symmetry_length$}\
943 {}{:mo_mirror_parities_length$} \
944 {eig_above_str:<mo_eig_above_length$} \
945 {eig_below_str:<mo_eig_below_length$} \
946 {mo_symmetry_projections_str:<mo_sym_projections_length$} \
947 {mo_density_symmetries_str}",
948 " ".repeat(mo_mirror_parities_gap),
949 mo_mirror_parities_str,
950 )?;
951 }
952 }
953 }
954 }
955
956 writeln!(f, "{}", "┈".repeat(table_width))?;
957 writeln!(f)?;
958 }
959
960 Ok(())
961 }
962}
963
964impl<'a, G, T, SC> fmt::Debug for SlaterDeterminantRepAnalysisResult<'a, G, T, SC>
965where
966 G: SymmetryGroupProperties + Clone,
967 G::CharTab: SubspaceDecomposable<T>,
968 T: ComplexFloat + Lapack,
969 SC: StructureConstraint + fmt::Display,
970 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
971{
972 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
973 writeln!(f, "{self}")
974 }
975}
976
977#[derive(Clone, Builder)]
987#[builder(build_fn(validate = "Self::validate"))]
988pub struct SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
989where
990 G: SymmetryGroupProperties + Clone,
991 G::CharTab: SubspaceDecomposable<T>,
992 T: ComplexFloat + Lapack,
993 SC: StructureConstraint + fmt::Display,
994 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
995{
996 parameters: &'a SlaterDeterminantRepAnalysisParams<<T as ComplexFloat>::Real>,
998
999 determinant: &'a SlaterDeterminant<'a, T, SC>,
1001
1002 symmetry_group: &'a SymmetryGroupDetectionResult,
1005
1006 sao: &'a Array2<T>,
1011
1012 #[builder(default = "None")]
1019 sao_h: Option<&'a Array2<T>>,
1020
1021 #[builder(default = "None")]
1024 sao_spatial_4c: Option<&'a Array4<T>>,
1025
1026 #[builder(default = "None")]
1031 sao_spatial_4c_h: Option<&'a Array4<T>>,
1032
1033 angular_function_parameters: &'a AngularFunctionRepAnalysisParams,
1035
1036 #[builder(setter(skip), default = "None")]
1038 result: Option<SlaterDeterminantRepAnalysisResult<'a, G, T, SC>>,
1039}
1040
1041impl<'a, G, T, SC> SlaterDeterminantRepAnalysisDriverBuilder<'a, G, T, SC>
1042where
1043 G: SymmetryGroupProperties + Clone,
1044 G::CharTab: SubspaceDecomposable<T>,
1045 T: ComplexFloat + Lapack,
1046 SC: StructureConstraint + fmt::Display,
1047 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1048{
1049 fn validate(&self) -> Result<(), String> {
1050 let params = self
1051 .parameters
1052 .ok_or("No Slater determinant representation analysis parameters found.".to_string())?;
1053
1054 let sym_res = self
1055 .symmetry_group
1056 .ok_or("No symmetry group information found.".to_string())?;
1057
1058 let sao = self.sao.ok_or("No SAO matrix found.".to_string())?;
1059
1060 if let Some(sao_h) = self.sao_h.flatten() {
1061 if sao_h.shape() != sao.shape() {
1062 return Err("Mismatched shapes between `sao` and `sao_h`.".to_string());
1063 }
1064 }
1065
1066 match (
1067 self.sao_spatial_4c.flatten(),
1068 self.sao_spatial_4c_h.flatten(),
1069 ) {
1070 (Some(sao_spatial_4c), Some(sao_spatial_4c_h)) => {
1071 if sao_spatial_4c_h.shape() != sao_spatial_4c.shape() {
1072 return Err(
1073 "Mismatched shapes between `sao_spatial_4c` and `sao_spatial_4c_h`."
1074 .to_string(),
1075 );
1076 }
1077 }
1078 (None, Some(_)) => {
1079 return Err("`sao_spatial_4c_h` is provided without `sao_spatial_4c`.".to_string());
1080 }
1081 _ => {}
1082 }
1083
1084 let det = self
1085 .determinant
1086 .ok_or("No Slater determinant found.".to_string())?;
1087
1088 let sym = if params.use_magnetic_group.is_some() {
1089 sym_res
1090 .magnetic_symmetry
1091 .as_ref()
1092 .ok_or("Magnetic symmetry requested for representation analysis, but no magnetic symmetry found.")?
1093 } else {
1094 &sym_res.unitary_symmetry
1095 };
1096
1097 if sym.is_infinite() && params.infinite_order_to_finite.is_none() {
1098 Err(
1099 format!(
1100 "Representation analysis cannot be performed using the entirety of the infinite group `{}`. \
1101 Consider setting the parameter `infinite_order_to_finite` to restrict to a finite subgroup instead.",
1102 sym.group_name.as_ref().expect("No symmetry group name found.")
1103 )
1104 )
1105 } else if (det.bao().n_funcs() != sao.nrows() || det.bao().n_funcs() != sao.ncols())
1106 && (det.bao().n_funcs()
1107 * det
1108 .structure_constraint()
1109 .n_explicit_comps_per_coefficient_matrix()
1110 != sao.nrows()
1111 || det.bao().n_funcs()
1112 * det
1113 .structure_constraint()
1114 .n_explicit_comps_per_coefficient_matrix()
1115 != sao.ncols())
1116 {
1117 Err("The dimensions of the SAO matrix do not match either the number of spatial AO basis functions or the number of spatial AO basis functions multiplied by the number of explicit components per coefficient matrix.".to_string())
1118 } else {
1119 Ok(())
1120 }
1121 }
1122}
1123
1124impl<'a, G, T, SC> SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1132where
1133 G: SymmetryGroupProperties + Clone,
1134 G::CharTab: SubspaceDecomposable<T>,
1135 T: ComplexFloat + Lapack,
1136 SC: StructureConstraint + Clone + fmt::Display,
1137 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1138{
1139 pub fn builder() -> SlaterDeterminantRepAnalysisDriverBuilder<'a, G, T, SC> {
1141 SlaterDeterminantRepAnalysisDriverBuilder::default()
1142 }
1143
1144 fn construct_sao(&self) -> Result<(Array2<T>, Option<Array2<T>>), anyhow::Error> {
1147 let nbas = self.determinant.bao().n_funcs();
1148 let ncomps = self
1149 .determinant
1150 .structure_constraint()
1151 .n_explicit_comps_per_coefficient_matrix();
1152 let provided_dim = self.sao.nrows();
1153
1154 if provided_dim == nbas {
1155 let sao = {
1156 let mut sao_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
1157 (0..ncomps).for_each(|icomp| {
1158 let start = icomp * nbas;
1159 let end = (icomp + 1) * nbas;
1160 sao_mut
1161 .slice_mut(s![start..end, start..end])
1162 .assign(self.sao);
1163 });
1164 sao_mut
1165 };
1166
1167 let sao_h = self.sao_h.map(|sao_h| {
1168 let mut sao_h_mut = Array2::zeros((ncomps * nbas, ncomps * nbas));
1169 (0..ncomps).for_each(|icomp| {
1170 let start = icomp * nbas;
1171 let end = (icomp + 1) * nbas;
1172 sao_h_mut
1173 .slice_mut(s![start..end, start..end])
1174 .assign(sao_h);
1175 });
1176 sao_h_mut
1177 });
1178
1179 Ok((sao, sao_h))
1180 } else {
1181 assert_eq!(provided_dim, nbas * ncomps);
1182 Ok((self.sao.clone(), self.sao_h.cloned()))
1183 }
1184 }
1185}
1186
1187impl<'a, T, SC> SlaterDeterminantRepAnalysisDriver<'a, UnitaryRepresentedSymmetryGroup, T, SC>
1191where
1192 T: ComplexFloat + Lapack + Sync + Send,
1193 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + Sync + Send,
1194 SC: StructureConstraint + fmt::Display,
1195 for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
1196{
1197 fn_construct_unitary_group!(
1198 construct_unitary_group
1201 );
1202}
1203
1204impl<'a, T, SC> SlaterDeterminantRepAnalysisDriver<'a, MagneticRepresentedSymmetryGroup, T, SC>
1208where
1209 T: ComplexFloat + Lapack + Sync + Send,
1210 <T as ComplexFloat>::Real: From<f64> + Sync + Send + fmt::LowerExp + fmt::Debug,
1211 SC: StructureConstraint + fmt::Display,
1212 for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
1213{
1214 fn_construct_magnetic_group!(
1215 construct_magnetic_group
1218 );
1219}
1220
1221#[duplicate_item(
1225 duplicate!{
1226 [
1227 dtype_nested sctype_nested;
1228 [ f64 ] [ SpinConstraint ];
1229 [ Complex<f64> ] [ SpinConstraint ];
1230 [ Complex<f64> ] [ SpinOrbitCoupled ];
1231 ]
1232 [
1233 gtype_ [ UnitaryRepresentedSymmetryGroup ]
1234 dtype_ [ dtype_nested ]
1235 sctype_ [ sctype_nested ]
1236 doc_sub_ [ "Performs representation analysis using a unitary-represented group and stores the result." ]
1237 analyse_fn_ [ analyse_representation ]
1238 construct_group_ [ self.construct_unitary_group()? ]
1239 calc_projections_ [
1240 log_subtitle("Slater determinant projection decompositions");
1241 qsym2_output!("");
1242 qsym2_output!(" Projections are defined w.r.t. the following inner product:");
1243 qsym2_output!(" {}", det_orbit.origin().overlap_definition());
1244 qsym2_output!("");
1245 det_orbit
1246 .projections_to_string(
1247 &det_orbit.calc_projection_compositions()?,
1248 params.integrality_threshold,
1249 )
1250 .log_output_display();
1251 qsym2_output!("");
1252 ]
1253 calc_mo_projections_ [
1254 let mo_projectionss = mo_orbitss.iter().map(|mo_orbits| {
1255 mo_orbits
1256 .iter()
1257 .map(|mo_orbit| mo_orbit.calc_projection_compositions().ok())
1258 .collect::<Vec<_>>()
1259 }).collect::<Vec<_>>();
1260 Some(mo_projectionss)
1261 ]
1262 ]
1263 }
1264 duplicate!{
1265 [
1266 dtype_nested sctype_nested;
1267 [ f64 ] [ SpinConstraint ];
1268 [ Complex<f64> ] [ SpinConstraint ];
1269 [ Complex<f64> ] [ SpinOrbitCoupled ];
1270 ]
1271 [
1272 gtype_ [ MagneticRepresentedSymmetryGroup ]
1273 dtype_ [ dtype_nested ]
1274 sctype_ [ sctype_nested ]
1275 doc_sub_ [ "Performs corepresentation analysis using a magnetic-represented group and stores the result." ]
1276 analyse_fn_ [ analyse_corepresentation ]
1277 construct_group_ [ self.construct_magnetic_group()? ]
1278 calc_projections_ [ ]
1279 calc_mo_projections_ [
1280 None::<Vec<Vec<Option<Vec<(MullikenIrcorepSymbol, Complex<f64>)>>>>>
1281 ]
1282 ]
1283 }
1284)]
1285impl<'a> SlaterDeterminantRepAnalysisDriver<'a, gtype_, dtype_, sctype_> {
1286 #[doc = doc_sub_]
1287 fn analyse_fn_(&mut self) -> Result<(), anyhow::Error> {
1288 let params = self.parameters;
1289 let (sao, sao_h) = self.construct_sao()?;
1290 let group = construct_group_;
1291 log_cc_transversal(&group);
1292 let _ = find_angular_function_representation(&group, self.angular_function_parameters);
1293 if group.is_double_group() {
1294 let _ = find_spinor_function_representation(&group, self.angular_function_parameters);
1295 }
1296 log_bao(self.determinant.bao());
1297
1298 let (
1300 det_symmetry,
1301 mo_symmetries,
1302 mo_symmetry_projections,
1303 mo_mirror_parities,
1304 mo_symmetries_thresholds,
1305 ) = if params.analyse_mo_symmetries {
1306 let mos = self.determinant.to_orbitals();
1307 let (mut det_orbit, mut mo_orbitss) = generate_det_mo_orbits(
1308 self.determinant,
1309 &mos,
1310 &group,
1311 &sao,
1312 sao_h.as_ref(),
1313 params.integrality_threshold,
1314 params.linear_independence_threshold,
1315 params.symmetry_transformation_kind.clone(),
1316 params.eigenvalue_comparison_mode.clone(),
1317 params.use_cayley_table,
1318 )?;
1319 det_orbit.calc_xmat(false)?;
1320 if params.write_overlap_eigenvalues {
1321 if let Some(smat_eigvals) = det_orbit.smat_eigvals.as_ref() {
1322 log_overlap_eigenvalues(
1323 "Determinant orbit overlap eigenvalues",
1324 smat_eigvals,
1325 params.linear_independence_threshold,
1326 ¶ms.eigenvalue_comparison_mode,
1327 );
1328 qsym2_output!("");
1329 }
1330 }
1331
1332 let det_symmetry = det_orbit.analyse_rep().map_err(|err| err.to_string());
1333
1334 {
1335 calc_projections_
1336 }
1337
1338 let mo_symmetries = mo_orbitss
1339 .iter_mut()
1340 .map(|mo_orbits| {
1341 mo_orbits
1342 .par_iter_mut()
1343 .map(|mo_orbit| {
1344 mo_orbit.calc_xmat(false).ok()?;
1345 mo_orbit.analyse_rep().ok()
1346 })
1347 .collect::<Vec<_>>()
1348 })
1349 .collect::<Vec<_>>();
1350
1351 let mo_symmetry_projections_opt = if params.analyse_mo_symmetry_projections {
1352 calc_mo_projections_
1353 } else {
1354 None
1355 };
1356
1357 let mo_mirror_parities_opt = if params.analyse_mo_mirror_parities {
1358 Some(
1359 mo_symmetries
1360 .iter()
1361 .map(|spin_mo_symmetries| {
1362 spin_mo_symmetries
1363 .iter()
1364 .map(|mo_sym_opt| {
1365 mo_sym_opt.as_ref().map(|mo_sym| {
1366 deduce_mirror_parities(det_orbit.group(), mo_sym)
1367 })
1368 })
1369 .collect::<Vec<_>>()
1370 })
1371 .collect::<Vec<_>>(),
1372 )
1373 } else {
1374 None
1375 };
1376
1377 let mo_symmetries_thresholds = mo_orbitss
1378 .iter_mut()
1379 .map(|mo_orbits| {
1380 mo_orbits
1381 .par_iter_mut()
1382 .map(|mo_orbit| {
1383 mo_orbit
1384 .smat_eigvals
1385 .as_ref()
1386 .map(|eigvals| {
1387 let mut eigvals_vec = eigvals.iter().collect::<Vec<_>>();
1388 match mo_orbit.eigenvalue_comparison_mode {
1389 EigenvalueComparisonMode::Modulus => {
1390 eigvals_vec.sort_by(|a, b| {
1391 a.abs().partial_cmp(&b.abs()).expect("Unable to compare two eigenvalues based on their moduli.")
1392 });
1393 }
1394 EigenvalueComparisonMode::Real => {
1395 eigvals_vec.sort_by(|a, b| {
1396 a.re().partial_cmp(&b.re()).expect("Unable to compare two eigenvalues based on their real parts.")
1397 });
1398 }
1399 }
1400 let eigval_above = match mo_orbit.eigenvalue_comparison_mode {
1401 EigenvalueComparisonMode::Modulus => eigvals_vec
1402 .iter()
1403 .find(|val| {
1404 val.abs() >= mo_orbit.linear_independence_threshold
1405 })
1406 .copied()
1407 .copied(),
1408 EigenvalueComparisonMode::Real => eigvals_vec
1409 .iter()
1410 .find(|val| {
1411 val.re() >= mo_orbit.linear_independence_threshold
1412 })
1413 .copied()
1414 .copied(),
1415 };
1416 eigvals_vec.reverse();
1417 let eigval_below = match mo_orbit.eigenvalue_comparison_mode {
1418 EigenvalueComparisonMode::Modulus => eigvals_vec
1419 .iter()
1420 .find(|val| {
1421 val.abs() < mo_orbit.linear_independence_threshold
1422 })
1423 .copied()
1424 .copied(),
1425 EigenvalueComparisonMode::Real => eigvals_vec
1426 .iter()
1427 .find(|val| {
1428 val.re() < mo_orbit.linear_independence_threshold
1429 })
1430 .copied()
1431 .copied(),
1432 };
1433 (eigval_above, eigval_below)
1434 })
1435 .unwrap_or((None, None))
1436 })
1437 .collect::<Vec<_>>()
1438 })
1439 .collect::<Vec<_>>();
1440 (
1441 det_symmetry,
1442 Some(mo_symmetries),
1443 mo_symmetry_projections_opt,
1444 mo_mirror_parities_opt,
1445 Some(mo_symmetries_thresholds),
1446 ) } else {
1448 let mut det_orbit = SlaterDeterminantSymmetryOrbit::builder()
1449 .group(&group)
1450 .origin(self.determinant)
1451 .integrality_threshold(params.integrality_threshold)
1452 .linear_independence_threshold(params.linear_independence_threshold)
1453 .symmetry_transformation_kind(params.symmetry_transformation_kind.clone())
1454 .eigenvalue_comparison_mode(params.eigenvalue_comparison_mode.clone())
1455 .build()?;
1456 let det_symmetry = det_orbit
1457 .calc_smat(Some(&sao), sao_h.as_ref(), params.use_cayley_table)
1458 .and_then(|det_orb| det_orb.normalise_smat())
1459 .map_err(|err| err.to_string())
1460 .and_then(|det_orb| {
1461 det_orb.calc_xmat(false).map_err(|err| err.to_string())?;
1462 if params.write_overlap_eigenvalues {
1463 if let Some(smat_eigvals) = det_orb.smat_eigvals.as_ref() {
1464 log_overlap_eigenvalues(
1465 "Determinant orbit overlap eigenvalues",
1466 smat_eigvals,
1467 params.linear_independence_threshold,
1468 ¶ms.eigenvalue_comparison_mode,
1469 );
1470 qsym2_output!("");
1471 }
1472 }
1473 det_orb.analyse_rep().map_err(|err| err.to_string())
1474 });
1475
1476 {
1477 calc_projections_
1478 }
1479
1480 (det_symmetry, None, None, None, None) };
1482
1483 let (den_symmetries, mo_den_symmetries) = if params.analyse_density_symmetries {
1485 let den_syms = self.determinant.to_densities().map(|densities| {
1486 let mut spin_den_syms = densities
1487 .iter()
1488 .enumerate()
1489 .map(|(ispin, den)| {
1490 let den_sym_res = || {
1491 let mut den_orbit = DensitySymmetryOrbit::builder()
1492 .group(&group)
1493 .origin(den)
1494 .integrality_threshold(params.integrality_threshold)
1495 .linear_independence_threshold(params.linear_independence_threshold)
1496 .symmetry_transformation_kind(
1497 params.symmetry_transformation_kind.clone(),
1498 )
1499 .eigenvalue_comparison_mode(
1500 params.eigenvalue_comparison_mode.clone(),
1501 )
1502 .build()?;
1503 den_orbit
1504 .calc_smat(
1505 self.sao_spatial_4c,
1506 self.sao_spatial_4c_h,
1507 params.use_cayley_table,
1508 )?
1509 .normalise_smat()?
1510 .calc_xmat(false)?;
1511 den_orbit.analyse_rep().map_err(|err| format_err!(err))
1512 };
1513 (
1514 format!("Spin-{ispin} density"),
1515 den_sym_res().map_err(|err| err.to_string()),
1516 )
1517 })
1518 .collect::<Vec<_>>();
1519 let mut extra_den_syms = densities
1520 .calc_extra_densities()
1521 .expect("Unable to calculate extra densities.")
1522 .iter()
1523 .map(|(label, den)| {
1524 let den_sym_res = || {
1525 let mut den_orbit = DensitySymmetryOrbit::builder()
1526 .group(&group)
1527 .origin(den)
1528 .integrality_threshold(params.integrality_threshold)
1529 .linear_independence_threshold(params.linear_independence_threshold)
1530 .symmetry_transformation_kind(
1531 params.symmetry_transformation_kind.clone(),
1532 )
1533 .eigenvalue_comparison_mode(
1534 params.eigenvalue_comparison_mode.clone(),
1535 )
1536 .build()?;
1537 den_orbit
1538 .calc_smat(
1539 self.sao_spatial_4c,
1540 self.sao_spatial_4c_h,
1541 params.use_cayley_table,
1542 )?
1543 .calc_xmat(false)?;
1544 den_orbit.analyse_rep().map_err(|err| format_err!(err))
1545 };
1546 (label.clone(), den_sym_res().map_err(|err| err.to_string()))
1547 })
1548 .collect_vec();
1549 spin_den_syms.append(&mut extra_den_syms);
1632 spin_den_syms
1633 });
1634
1635 let mo_den_syms = if params.analyse_mo_symmetries {
1636 let mo_den_symmetries = self
1637 .determinant
1638 .to_orbitals()
1639 .iter()
1640 .map(|mos| {
1641 mos.par_iter()
1642 .map(|mo| {
1643 let mo_den = mo.to_total_density().ok()?;
1644 let mut mo_den_orbit = DensitySymmetryOrbit::builder()
1645 .group(&group)
1646 .origin(&mo_den)
1647 .integrality_threshold(params.integrality_threshold)
1648 .linear_independence_threshold(
1649 params.linear_independence_threshold,
1650 )
1651 .symmetry_transformation_kind(
1652 params.symmetry_transformation_kind.clone(),
1653 )
1654 .eigenvalue_comparison_mode(
1655 params.eigenvalue_comparison_mode.clone(),
1656 )
1657 .build()
1658 .ok()?;
1659 log::debug!("Computing overlap matrix for an MO density orbit...");
1660 mo_den_orbit
1661 .calc_smat(
1662 self.sao_spatial_4c,
1663 self.sao_spatial_4c_h,
1664 params.use_cayley_table,
1665 )
1666 .ok()?
1667 .normalise_smat()
1668 .ok()?
1669 .calc_xmat(false)
1670 .ok()?;
1671 log::debug!(
1672 "Computing overlap matrix for an MO density orbit... Done."
1673 );
1674 mo_den_orbit.analyse_rep().ok()
1675 })
1676 .collect::<Vec<_>>()
1677 })
1678 .collect::<Vec<_>>();
1679 Some(mo_den_symmetries)
1680 } else {
1681 None
1682 };
1683
1684 (den_syms.ok(), mo_den_syms)
1685 } else {
1686 (None, None)
1687 };
1688
1689 let result = SlaterDeterminantRepAnalysisResult::builder()
1690 .parameters(params)
1691 .determinant(self.determinant)
1692 .group(group)
1693 .determinant_symmetry(det_symmetry)
1694 .determinant_density_symmetries(den_symmetries)
1695 .mo_symmetries(mo_symmetries)
1696 .mo_symmetry_projections(mo_symmetry_projections)
1697 .mo_mirror_parities(mo_mirror_parities)
1698 .mo_symmetries_thresholds(mo_symmetries_thresholds)
1699 .mo_density_symmetries(mo_den_symmetries)
1700 .build()?;
1701 self.result = Some(result);
1702
1703 Ok(())
1704 }
1705}
1706
1707impl<'a, G, T, SC> fmt::Display for SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1715where
1716 G: SymmetryGroupProperties + Clone,
1717 G::CharTab: SubspaceDecomposable<T>,
1718 T: ComplexFloat + Lapack,
1719 SC: StructureConstraint + fmt::Display,
1720 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1721{
1722 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1723 write_title(f, "Slater Determinant Symmetry Analysis")?;
1724 writeln!(f)?;
1725 writeln!(f, "{}", self.parameters)?;
1726 Ok(())
1727 }
1728}
1729
1730impl<'a, G, T, SC> fmt::Debug for SlaterDeterminantRepAnalysisDriver<'a, G, T, SC>
1731where
1732 G: SymmetryGroupProperties + Clone,
1733 G::CharTab: SubspaceDecomposable<T>,
1734 T: ComplexFloat + Lapack,
1735 SC: StructureConstraint + fmt::Display,
1736 <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
1737{
1738 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1739 writeln!(f, "{self}")
1740 }
1741}
1742
1743#[duplicate_item(
1747 duplicate!{
1748 [
1749 dtype_nested sctype_nested;
1750 [ f64 ] [ SpinConstraint ];
1751 [ Complex<f64> ] [ SpinConstraint ];
1752 [ Complex<f64> ] [ SpinOrbitCoupled ];
1753 ]
1754 [
1755 gtype_ [ UnitaryRepresentedSymmetryGroup ]
1756 dtype_ [ dtype_nested ]
1757 sctype_ [ sctype_nested ]
1758 analyse_fn_ [ analyse_representation ]
1759 ]
1760 }
1761 duplicate!{
1762 [
1763 dtype_nested sctype_nested;
1764 [ f64 ] [ SpinConstraint ];
1765 [ Complex<f64> ] [ SpinConstraint ];
1766 [ Complex<f64> ] [ SpinOrbitCoupled ];
1767 ]
1768 [
1769 gtype_ [ MagneticRepresentedSymmetryGroup ]
1770 dtype_ [ dtype_nested ]
1771 sctype_ [ sctype_nested ]
1772 analyse_fn_ [ analyse_corepresentation ]
1773 ]
1774 }
1775)]
1776impl<'a> QSym2Driver for SlaterDeterminantRepAnalysisDriver<'a, gtype_, dtype_, sctype_> {
1777 type Params = SlaterDeterminantRepAnalysisParams<f64>;
1778
1779 type Outcome = SlaterDeterminantRepAnalysisResult<'a, gtype_, dtype_, sctype_>;
1780
1781 fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
1782 self.result.as_ref().ok_or_else(|| {
1783 format_err!("No Slater determinant representation analysis results found.")
1784 })
1785 }
1786
1787 fn run(&mut self) -> Result<(), anyhow::Error> {
1788 self.log_output_display();
1789 self.analyse_fn_()?;
1790 self.result()?.log_output_display();
1791 Ok(())
1792 }
1793}