qsym2/interfaces/qchem/hdf5/
slater_determinant.rs

1//! Slater determinants from Q-Chem HDF5 archives.
2
3use std::fmt;
4use std::marker::PhantomData;
5use std::path::PathBuf;
6
7use anyhow::{self, bail, format_err, Context};
8use derive_builder::Builder;
9use duplicate::duplicate_item;
10use factorial::DoubleFactorial;
11use hdf5::{self, H5Type};
12use lazy_static::lazy_static;
13use log;
14use nalgebra::Point3;
15use ndarray::{s, Array1, Array2, Array4, Axis, Ix3};
16use ndarray_linalg::types::Lapack;
17use num_complex::ComplexFloat;
18use num_traits::{One, ToPrimitive, Zero};
19use numeric_sort;
20use periodic_table::periodic_table;
21use regex::Regex;
22
23use crate::angmom::spinor_rotation_3d::SpinConstraint;
24use crate::auxiliary::atom::{Atom, ElementMap};
25use crate::auxiliary::molecule::Molecule;
26use crate::basis::ao::*;
27use crate::basis::ao_integrals::*;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::SubspaceDecomposable;
30use crate::drivers::representation_analysis::angular_function::AngularFunctionRepAnalysisParams;
31use crate::drivers::representation_analysis::slater_determinant::{
32    SlaterDeterminantRepAnalysisDriver, SlaterDeterminantRepAnalysisParams,
33};
34use crate::drivers::representation_analysis::MagneticSymmetryAnalysisKind;
35use crate::drivers::symmetry_group_detection::{
36    SymmetryGroupDetectionDriver, SymmetryGroupDetectionResult,
37};
38use crate::drivers::QSym2Driver;
39#[cfg(feature = "integrals")]
40use crate::integrals::shell_tuple::build_shell_tuple_collection;
41use crate::interfaces::input::SymmetryGroupDetectionInputKind;
42use crate::io::format::{
43    log_macsec_begin, log_macsec_end, log_micsec_begin, log_micsec_end, qsym2_error, qsym2_output,
44};
45use crate::io::{read_qsym2_binary, QSym2FileType};
46use crate::symmetry::symmetry_core::Symmetry;
47use crate::symmetry::symmetry_group::{
48    MagneticRepresentedSymmetryGroup, SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup,
49};
50use crate::target::determinant::SlaterDeterminant;
51
52#[cfg(test)]
53#[path = "slater_determinant_tests.rs"]
54mod slater_determinant_tests;
55
56// =====================
57// Full Q-Chem H5 Driver
58// =====================
59
60lazy_static! {
61    static ref SP_PATH_RE: Regex =
62        Regex::new(r"(.*sp)\\energy_function$").expect("Regex pattern invalid.");
63}
64
65// -----------------
66// Struct definition
67// -----------------
68
69/// Driver to perform symmetry-group detection and Slater determinant representation symmetry
70/// analysis for all discoverable single-point calculation data stored in a Q-Chem's `qarchive.h5`
71/// file.
72#[derive(Clone, Builder)]
73pub struct QChemSlaterDeterminantH5Driver<'a, T>
74where
75    T: Clone,
76{
77    /// The `qarchive.h5` file name.
78    filename: PathBuf,
79
80    /// The input specification controlling symmetry-group detection.
81    symmetry_group_detection_input: &'a SymmetryGroupDetectionInputKind,
82
83    /// The parameters controlling representation analysis of standard angular functions.
84    angular_function_analysis_parameters: &'a AngularFunctionRepAnalysisParams,
85
86    /// The parameters controlling representation analysis.
87    rep_analysis_parameters: &'a SlaterDeterminantRepAnalysisParams<f64>,
88
89    /// The simplified result of the analysis. Each element in the vector is a tuple containing the
90    /// group name and the representation symmetry of the Slater determinant for one single-point
91    /// calculation.
92    #[builder(default = "None")]
93    result: Option<Vec<(String, String)>>,
94
95    /// The numerical type of the Slater determinant.
96    #[builder(setter(skip), default = "PhantomData")]
97    numerical_type: PhantomData<T>,
98}
99
100// ----------------------
101// Struct implementations
102// ----------------------
103
104impl<'a, T> QChemSlaterDeterminantH5Driver<'a, T>
105where
106    T: Clone,
107{
108    /// Returns a builder to construct a [`QChemSlaterDeterminantH5Driver`].
109    pub fn builder() -> QChemSlaterDeterminantH5DriverBuilder<'a, T> {
110        QChemSlaterDeterminantH5DriverBuilder::default()
111    }
112}
113
114// ~~~~~~~~~~~~~~~~~~~
115// Slater determinants
116// ~~~~~~~~~~~~~~~~~~~
117
118// Specific for Slater determinant numeric type f64
119// ''''''''''''''''''''''''''''''''''''''''''''''''
120impl<'a> QChemSlaterDeterminantH5Driver<'a, f64> {
121    /// Performs analysis for all real-valued single-point determinants.
122    fn analyse(&mut self) -> Result<(), anyhow::Error> {
123        let f = hdf5::File::open(&self.filename)?;
124        let mut sp_paths = f
125            .group(".counters")?
126            .member_names()?
127            .iter()
128            .filter_map(|path| {
129                if SP_PATH_RE.is_match(path) {
130                    let path = path.replace("\\", "/");
131                    let mut energy_function_indices = f
132                        .group(&path)
133                        .and_then(|sp_energy_function_group| {
134                            sp_energy_function_group.member_names()
135                        })
136                        .ok()?;
137                    numeric_sort::sort(&mut energy_function_indices);
138                    Some((
139                        path.replace("/energy_function", ""),
140                        energy_function_indices,
141                    ))
142                } else {
143                    None
144                }
145            })
146            .collect::<Vec<_>>();
147        sp_paths.sort_by(|(path_a, _), (path_b, _)| numeric_sort::cmp(path_a, path_b));
148
149        let pd_input = self.symmetry_group_detection_input;
150        let afa_params = self.angular_function_analysis_parameters;
151        let sda_params = self.rep_analysis_parameters;
152        let result = sp_paths
153            .iter()
154            .flat_map(|(sp_path, energy_function_indices)| {
155                energy_function_indices.iter().map(|energy_function_index| {
156                    log_macsec_begin(&format!(
157                        "Analysis for {} (energy function {energy_function_index})",
158                        sp_path.clone()
159                    ));
160                    qsym2_output!("");
161                    let sp = f.group(sp_path)?;
162                    let sp_driver_result = match sda_params.use_magnetic_group {
163                        Some(MagneticSymmetryAnalysisKind::Corepresentation) => {
164                            let mut sp_driver = QChemSlaterDeterminantH5SinglePointDriver::<
165                                MagneticRepresentedSymmetryGroup,
166                                f64,
167                            >::builder()
168                            .sp_group(&sp)
169                            .energy_function_index(energy_function_index)
170                            .symmetry_group_detection_input(pd_input)
171                            .angular_function_analysis_parameters(afa_params)
172                            .rep_analysis_parameters(sda_params)
173                            .build()?;
174                            let _ = sp_driver.run();
175                            sp_driver.result().map(|(sym, rep)| {
176                                (
177                                    sym.group_name
178                                        .as_ref()
179                                        .unwrap_or(&String::new())
180                                        .to_string(),
181                                    rep.as_ref()
182                                        .map(|rep| rep.to_string())
183                                        .unwrap_or_else(|err| err.to_string()),
184                                )
185                            })
186                        }
187                        Some(MagneticSymmetryAnalysisKind::Representation) | None => {
188                            let mut sp_driver = QChemSlaterDeterminantH5SinglePointDriver::<
189                                UnitaryRepresentedSymmetryGroup,
190                                f64,
191                            >::builder()
192                            .sp_group(&sp)
193                            .energy_function_index(energy_function_index)
194                            .symmetry_group_detection_input(pd_input)
195                            .angular_function_analysis_parameters(afa_params)
196                            .rep_analysis_parameters(sda_params)
197                            .build()?;
198                            let _ = sp_driver.run();
199                            sp_driver.result().map(|(sym, rep)| {
200                                (
201                                    sym.group_name
202                                        .as_ref()
203                                        .unwrap_or(&String::new())
204                                        .to_string(),
205                                    rep.as_ref()
206                                        .map(|rep| rep.to_string())
207                                        .unwrap_or_else(|err| err.to_string()),
208                                )
209                            })
210                        }
211                    };
212                    qsym2_output!("");
213                    log_macsec_end(&format!(
214                        "Analysis for {} (energy function {energy_function_index})",
215                        sp_path.clone()
216                    ));
217                    qsym2_output!("");
218                    sp_driver_result
219                })
220            })
221            .map(|res| {
222                res.unwrap_or_else(|err| {
223                    (
224                        "Unidentified symmetry group".to_string(),
225                        format!("Unidentified (co)representation: {err}"),
226                    )
227                })
228            })
229            .collect::<Vec<_>>();
230
231        log_macsec_begin("Q-Chem HDF5 Archive Summary");
232        qsym2_output!("");
233        let path_length = sp_paths
234            .iter()
235            .map(|(path, _)| path.chars().count())
236            .max()
237            .unwrap_or(18)
238            .max(18);
239        let energy_function_length = sp_paths
240            .iter()
241            .map(|(_, energy_function_indices)| {
242                energy_function_indices
243                    .iter()
244                    .map(|index| index.chars().count())
245                    .max()
246                    .unwrap_or(1)
247                    .max(1)
248            })
249            .max()
250            .unwrap_or(7)
251            .max(7);
252        let group_length = result
253            .iter()
254            .map(|(group, _)| group.chars().count())
255            .max()
256            .unwrap_or(5)
257            .max(5);
258        let sym_length = result
259            .iter()
260            .map(|(_, sym)| sym.chars().count())
261            .max()
262            .unwrap_or(13)
263            .max(13);
264        let table_width = path_length + energy_function_length + group_length + sym_length + 8;
265        qsym2_output!("{}", "┈".repeat(table_width));
266        qsym2_output!(
267            " {:<path_length$}  {:<energy_function_length$}  {:<group_length$}  {:<}",
268            "Single-point calc.",
269            "E func.",
270            "Group",
271            "Det. symmetry"
272        );
273        qsym2_output!("{}", "┈".repeat(table_width));
274        sp_paths
275            .iter()
276            .flat_map(|(sp_path, energy_function_indices)| {
277                energy_function_indices
278                    .iter()
279                    .map(|index| (sp_path.clone(), index))
280            })
281            .zip(result.iter())
282            .for_each(|((path, index), (group, sym))| {
283                qsym2_output!(
284                    " {:<path_length$}  {:<energy_function_length$}  {:<group_length$}  {:<#}",
285                    path,
286                    index,
287                    group,
288                    sym
289                );
290            });
291        qsym2_output!("{}", "┈".repeat(table_width));
292        qsym2_output!("");
293        log_macsec_end("Q-Chem HDF5 Archive Summary");
294
295        self.result = Some(result);
296        Ok(())
297    }
298}
299
300impl<'a> QSym2Driver for QChemSlaterDeterminantH5Driver<'a, f64> {
301    type Params = SlaterDeterminantRepAnalysisParams<f64>;
302
303    type Outcome = Vec<(String, String)>;
304
305    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
306        self.result.as_ref().ok_or(format_err!(
307            "No Q-Chem HDF5 analysis results for a real Slater determinant found."
308        ))
309    }
310
311    fn run(&mut self) -> Result<(), anyhow::Error> {
312        self.analyse()
313    }
314}
315
316// ==================
317// SinglePoint Driver
318// ==================
319
320// ---------------
321// Enum definition
322// ---------------
323
324/// Enumerated type to distinguish different kinds of molecular orbitals.
325pub enum OrbitalType {
326    /// Canonical molecular orbitals as obtained by diagonalising Fock matrices.
327    Canonical,
328
329    /// Localised molecular orbitals as obtained by a localisation method.
330    Localised,
331}
332
333// -----------------
334// Struct definition
335// -----------------
336
337/// Driver to perform symmetry-group detection and representation analysis for a single-point
338/// calculation result in a Q-Chem's `qarchive.h5` file.
339#[derive(Clone, Builder)]
340pub struct QChemSlaterDeterminantH5SinglePointDriver<'a, G, T>
341where
342    G: SymmetryGroupProperties + Clone,
343    G::CharTab: SubspaceDecomposable<T>,
344    T: ComplexFloat + Lapack,
345    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
346{
347    /// A H5 group containing data from a single-point calculation.
348    sp_group: &'a hdf5::Group,
349
350    /// The index of the energy function whose results are to be considered for this single-point
351    /// symmetry analysis.
352    #[builder(setter(custom))]
353    energy_function_index: String,
354
355    /// The parameters controlling symmetry-group detection.
356    symmetry_group_detection_input: &'a SymmetryGroupDetectionInputKind,
357
358    /// The parameters controlling representation analysis of standard angular functions.
359    angular_function_analysis_parameters: &'a AngularFunctionRepAnalysisParams,
360
361    /// The parameters controlling representation analysis of Slater determinants.
362    rep_analysis_parameters: &'a SlaterDeterminantRepAnalysisParams<f64>,
363
364    /// The symmetry of the system and the representation of the Slater determinant.
365    #[builder(default = "None")]
366    result: Option<(
367        Symmetry,
368        Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
369    )>,
370}
371
372// ----------------------
373// Struct implementations
374// ----------------------
375
376impl<'a, G, T> QChemSlaterDeterminantH5SinglePointDriverBuilder<'a, G, T>
377where
378    G: SymmetryGroupProperties + Clone,
379    G::CharTab: SubspaceDecomposable<T>,
380    T: ComplexFloat + Lapack,
381    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
382{
383    pub fn energy_function_index(&mut self, idx: &str) -> &mut Self {
384        self.energy_function_index = Some(idx.to_string());
385        self
386    }
387}
388
389impl<'a, G, T> QChemSlaterDeterminantH5SinglePointDriver<'a, G, T>
390where
391    G: SymmetryGroupProperties + Clone,
392    G::CharTab: SubspaceDecomposable<T>,
393    T: ComplexFloat + Lapack + H5Type,
394    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
395{
396    /// Returns a builder to construct a [`QChemSlaterDeterminantH5SinglePointDriver`].
397    pub fn builder() -> QChemSlaterDeterminantH5SinglePointDriverBuilder<'a, G, T> {
398        QChemSlaterDeterminantH5SinglePointDriverBuilder::default()
399    }
400
401    /// Extracts the molecular structure from the single-point H5 group.
402    pub fn extract_molecule(&self) -> Result<Molecule, anyhow::Error> {
403        let emap = ElementMap::new();
404        let coordss = self
405            .sp_group
406            .dataset("structure/coordinates")?
407            .read_2d::<f64>()?;
408        let atomic_numbers = self
409            .sp_group
410            .dataset("structure/nuclei")?
411            .read_1d::<usize>()?;
412        let atoms = coordss
413            .rows()
414            .into_iter()
415            .zip(atomic_numbers.iter())
416            .map(|(coords, atomic_number)| {
417                let element = periodic_table()
418                    .get(*atomic_number - 1)
419                    .ok_or(hdf5::Error::from(
420                        format!(
421                            "Element with atomic number {atomic_number} could not be identified."
422                        )
423                        .as_str(),
424                    ))?
425                    .symbol;
426                let coordinates = Point3::new(coords[0], coords[1], coords[2]);
427                Ok::<_, hdf5::Error>(Atom::new_ordinary(element, coordinates, &emap, 1e-8))
428            })
429            .collect::<Result<Vec<Atom>, _>>()?;
430        let mol = Molecule::from_atoms(&atoms, 1e-14);
431        Ok(mol)
432    }
433
434    /// Extracts the basis angular order information from the single-point H5 group.
435    pub fn extract_bao(&self, mol: &'a Molecule) -> Result<BasisAngularOrder<'a>, anyhow::Error> {
436        let shell_types = self
437            .sp_group
438            .dataset("aobasis/shell_types")?
439            .read_1d::<i32>()?;
440        let shell_to_atom_map = self
441            .sp_group
442            .dataset("aobasis/shell_to_atom_map")?
443            .read_1d::<usize>()?
444            .iter()
445            .zip(shell_types.iter())
446            .flat_map(|(&idx, shell_type)| {
447                if *shell_type == -1 {
448                    vec![idx, idx]
449                } else {
450                    vec![idx]
451                }
452            })
453            .collect::<Vec<_>>();
454
455        let bss: Vec<BasisShell> = shell_types
456            .iter()
457            .flat_map(|shell_type| {
458                if *shell_type == 0 {
459                    // S shell
460                    vec![BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)))]
461                } else if *shell_type == 1 {
462                    // P shell
463                    vec![BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)))]
464                } else if *shell_type == -1 {
465                    // SP shell
466                    vec![
467                        BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0))),
468                        BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1))),
469                    ]
470                } else if *shell_type < 0 {
471                    // Cartesian D shell or higher
472                    let l = shell_type.unsigned_abs();
473                    vec![BasisShell::new(l, ShellOrder::Cart(CartOrder::qchem(l)))]
474                } else {
475                    // Pure D shell or higher
476                    let l = shell_type.unsigned_abs();
477                    vec![BasisShell::new(
478                        l,
479                        ShellOrder::Pure(PureOrder::increasingm(l)),
480                    )]
481                }
482            })
483            .collect::<Vec<BasisShell>>();
484
485        let batms = mol
486            .atoms
487            .iter()
488            .enumerate()
489            .map(|(atom_i, atom)| {
490                let shells = bss
491                    .iter()
492                    .zip(shell_to_atom_map.iter())
493                    .filter_map(|(bs, atom_index)| {
494                        if *atom_index == atom_i {
495                            Some(bs.clone())
496                        } else {
497                            None
498                        }
499                    })
500                    .collect::<Vec<_>>();
501                BasisAtom::new(atom, &shells)
502            })
503            .collect::<Vec<BasisAtom>>();
504        Ok(BasisAngularOrder::new(&batms))
505    }
506
507    /// Extracts the spatial atomic-orbital overlap matrix from the single-point H5 group.
508    ///
509    /// Note that the overlap matrix in the HDF5 file uses lexicographic order for Cartesian
510    /// functions. This is inconsistent with the conventional Q-Chem ordering used for molecular
511    /// orbital coefficients. See [`Self::recompute_sao`] for a way to get the atomic-orbital
512    /// overlap matrix with the consistent Cartesian ordering.
513    pub fn extract_sao(&self) -> Result<Array2<T>, anyhow::Error> {
514        self.sp_group
515            .dataset("aobasis/overlap_matrix")?
516            .read_2d::<T>()
517            .map_err(|err| err.into())
518    }
519
520    /// Extracts the full basis set information from the single-point H5 group.
521    pub fn extract_basis_set(
522        &self,
523        mol: &'a Molecule,
524    ) -> Result<BasisSet<f64, f64>, anyhow::Error> {
525        let shell_types = self
526            .sp_group
527            .dataset("aobasis/shell_types")?
528            .read_1d::<i32>()?;
529        let shell_to_atom_map = self
530            .sp_group
531            .dataset("aobasis/shell_to_atom_map")?
532            .read_1d::<usize>()?
533            .iter()
534            .zip(shell_types.iter())
535            .flat_map(|(&idx, shell_type)| {
536                if *shell_type == -1 {
537                    vec![idx, idx]
538                } else {
539                    vec![idx]
540                }
541            })
542            .collect::<Vec<_>>();
543
544        let primitives_per_shell = self
545            .sp_group
546            .dataset("aobasis/primitives_per_shell")?
547            .read_1d::<usize>()?;
548        let contraction_coefficients = self
549            .sp_group
550            .dataset("aobasis/contraction_coefficients")?
551            .read_1d::<f64>()?;
552        let sp_contraction_coefficients = self
553            .sp_group
554            .dataset("aobasis/sp_contraction_coefficients")?
555            .read_1d::<f64>()?;
556        let primitive_exponents = self
557            .sp_group
558            .dataset("aobasis/primitive_exponents")?
559            .read_1d::<f64>()?;
560        // `shell_coordinates` in Bohr radius.
561        let shell_coordinates = self
562            .sp_group
563            .dataset("aobasis/shell_coordinates")?
564            .read_2d::<f64>()?;
565
566        let bscs: Vec<BasisShellContraction<f64, f64>> = primitives_per_shell
567            .iter()
568            .scan(0, |end, n_prims| {
569                let start = *end;
570                *end += n_prims;
571                Some((start, *end))
572            })
573            .zip(shell_types.iter())
574            .zip(shell_coordinates.rows())
575            .flat_map(|(((start, end), shell_type), centre)| {
576                if *shell_type == 0 {
577                    // S shell
578                    let basis_shell = BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)));
579                    let primitives = primitive_exponents
580                        .slice(s![start..end])
581                        .iter()
582                        .cloned()
583                        .zip(
584                            contraction_coefficients
585                                .slice(s![start..end])
586                                .iter()
587                                .cloned(),
588                        )
589                        .collect::<Vec<_>>();
590                    let contraction = GaussianContraction { primitives };
591                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
592                    vec![BasisShellContraction {
593                        basis_shell,
594                        contraction,
595                        cart_origin,
596                        k: None,
597                    }]
598                } else if *shell_type == 1 {
599                    // P shell
600                    let basis_shell = BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)));
601                    let primitives = primitive_exponents
602                        .slice(s![start..end])
603                        .iter()
604                        .cloned()
605                        .zip(
606                            contraction_coefficients
607                                .slice(s![start..end])
608                                .iter()
609                                .cloned(),
610                        )
611                        .collect::<Vec<_>>();
612                    let contraction = GaussianContraction { primitives };
613                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
614                    vec![BasisShellContraction {
615                        basis_shell,
616                        contraction,
617                        cart_origin,
618                        k: None,
619                    }]
620                } else if *shell_type == -1 {
621                    // SP shell
622                    let basis_shell_s = BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)));
623                    let primitives_s = primitive_exponents
624                        .slice(s![start..end])
625                        .iter()
626                        .cloned()
627                        .zip(
628                            contraction_coefficients
629                                .slice(s![start..end])
630                                .iter()
631                                .cloned(),
632                        )
633                        .collect::<Vec<_>>();
634                    let contraction_s = GaussianContraction {
635                        primitives: primitives_s,
636                    };
637
638                    let basis_shell_p = BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)));
639                    let primitives_p = primitive_exponents
640                        .slice(s![start..end])
641                        .iter()
642                        .cloned()
643                        .zip(
644                            sp_contraction_coefficients
645                                .slice(s![start..end])
646                                .iter()
647                                .cloned(),
648                        )
649                        .collect::<Vec<_>>();
650                    let contraction_p = GaussianContraction {
651                        primitives: primitives_p,
652                    };
653
654                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
655                    vec![
656                        BasisShellContraction {
657                            basis_shell: basis_shell_s,
658                            contraction: contraction_s,
659                            cart_origin: cart_origin.clone(),
660                            k: None,
661                        },
662                        BasisShellContraction {
663                            basis_shell: basis_shell_p,
664                            contraction: contraction_p,
665                            cart_origin,
666                            k: None,
667                        },
668                    ]
669                } else if *shell_type < 0 {
670                    // Cartesian D shell or higher
671                    let l = shell_type.unsigned_abs();
672                    let basis_shell = BasisShell::new(l, ShellOrder::Cart(CartOrder::qchem(l)));
673                    let primitives = primitive_exponents
674                        .slice(s![start..end])
675                        .iter()
676                        .cloned()
677                        .zip(
678                            contraction_coefficients
679                                .slice(s![start..end])
680                                .iter()
681                                .cloned(),
682                        )
683                        .collect::<Vec<_>>();
684                    let contraction = GaussianContraction { primitives };
685                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
686                    vec![BasisShellContraction {
687                        basis_shell,
688                        contraction,
689                        cart_origin,
690                        k: None,
691                    }]
692                } else {
693                    // Pure D shell or higher
694                    let l = shell_type.unsigned_abs();
695                    let basis_shell =
696                        BasisShell::new(l, ShellOrder::Pure(PureOrder::increasingm(l)));
697                    let primitives = primitive_exponents
698                        .slice(s![start..end])
699                        .iter()
700                        .cloned()
701                        .zip(
702                            contraction_coefficients
703                                .slice(s![start..end])
704                                .iter()
705                                .cloned(),
706                        )
707                        .collect::<Vec<_>>();
708                    let contraction = GaussianContraction { primitives };
709                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
710                    vec![BasisShellContraction {
711                        basis_shell,
712                        contraction,
713                        cart_origin,
714                        k: None,
715                    }]
716                }
717            })
718            .collect::<Vec<BasisShellContraction<f64, f64>>>();
719
720        let basis_atoms = mol
721            .atoms
722            .iter()
723            .enumerate()
724            .map(|(atom_i, _)| {
725                let shells = bscs
726                    .iter()
727                    .zip(shell_to_atom_map.iter())
728                    .filter_map(|(bs, atom_index)| {
729                        if *atom_index == atom_i {
730                            Some(bs.clone())
731                        } else {
732                            None
733                        }
734                    })
735                    .collect::<Vec<_>>();
736                shells
737            })
738            .collect::<Vec<Vec<BasisShellContraction<f64, f64>>>>();
739
740        let mut basis_set = BasisSet::<f64, f64>::new(basis_atoms);
741
742        // Q-Chem renormalises each Gaussian primitive, but this is not the convention used in
743        // QSym². We therefore un-normalise the Gaussian primitives to restore the original
744        // contraction coefficients.
745        let prefactor = (2.0 / std::f64::consts::PI).powf(0.75);
746        basis_set.all_shells_mut().for_each(|bsc| {
747            let l = bsc.basis_shell().l;
748            let l_i32 = l
749                .to_i32()
750                .unwrap_or_else(|| panic!("Unable to convert `{l}` to `i32`."));
751            let l_f64 = l
752                .to_f64()
753                .unwrap_or_else(|| panic!("Unable to convert `{l}` to `f64`."));
754            let doufac_sqrt = if l == 0 {
755                1
756            } else {
757                ((2 * l) - 1)
758                    .checked_double_factorial()
759                    .unwrap_or_else(|| panic!("Unable to obtain `{}!!`.", 2 * l - 1))
760            }
761            .to_f64()
762            .unwrap_or_else(|| panic!("Unable to convert `{}!!` to `f64`.", 2 * l - 1))
763            .sqrt();
764            bsc.contraction.primitives.iter_mut().for_each(|(a, c)| {
765                let n = prefactor * 2.0.powi(l_i32) * a.powf(l_f64 / 2.0 + 0.75) / doufac_sqrt;
766                *c /= n;
767            });
768        });
769        Ok(basis_set)
770    }
771
772    /// Recomputes the spatial atomic-orbital overlap matrix.
773    ///
774    /// The overlap matrix stored in the H5 group unfortunately uses lexicographic order for
775    /// Cartesian functions, which is inconsistent with that used in the coefficients. We thus
776    /// recompute the overlap matrix from the basis set information using the conventional Q-Chem
777    /// order for Cartesian functions.
778    pub fn recompute_sao(&self) -> Result<Array2<f64>, anyhow::Error> {
779        log::debug!("Recomputing atomic-orbital overlap matrix...");
780        let mol = self.extract_molecule()?;
781        let basis_set = self.extract_basis_set(&mol)?;
782        let stc = build_shell_tuple_collection![
783            <s1, s2>;
784            false, false;
785            &basis_set, &basis_set;
786            f64
787        ];
788        let sao_res = stc
789            .overlap([0, 0])
790            .pop()
791            .ok_or(format_err!("Unable to compute the AO overlap matrix."));
792        log::debug!("Recomputing atomic-orbital overlap matrix... Done.");
793        sao_res
794    }
795}
796
797// ~~~~~~~~~~~~~~~~~~~
798// Slater determinants
799// ~~~~~~~~~~~~~~~~~~~
800
801// Generic for all symmetry groups G and determinant numeric type T
802// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
803impl<'a, G, T> QChemSlaterDeterminantH5SinglePointDriver<'a, G, T>
804where
805    G: SymmetryGroupProperties + Clone,
806    G::CharTab: SubspaceDecomposable<T>,
807    T: ComplexFloat + Lapack + H5Type,
808    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
809{
810    /// Extracts the Slater determinant from the single-point H5 group.
811    ///
812    /// # Arguments
813    ///
814    /// * `mol` - The molecule to be associated with the extracted determinant.
815    /// * `bao` - The basis angular order information to be associated with the extracted determinant.
816    /// * `threshold` - The comparison threshold to be associated with the extracted determinant.
817    ///
818    /// # Returns
819    ///
820    /// The extracted Slater determinant.
821    pub fn extract_determinant(
822        &self,
823        mol: &'a Molecule,
824        bao: &'a BasisAngularOrder,
825        threshold: <T as ComplexFloat>::Real,
826        orbital_type: OrbitalType,
827    ) -> Result<SlaterDeterminant<'a, T, SpinConstraint>, anyhow::Error> {
828        let energy = self
829            .sp_group
830            .dataset(&format!(
831                "energy_function/{}/energy",
832                self.energy_function_index
833            ))?
834            .read_scalar::<T>()
835            .map_err(|err| err.to_string());
836        let orbital_path = match orbital_type {
837            OrbitalType::Canonical => format!(
838                "energy_function/{}/method/scf/molecular_orbitals",
839                self.energy_function_index
840            ),
841            OrbitalType::Localised => format!(
842                "energy_function/{}/analysis/localized_orbitals/{}/molecular_orbitals",
843                self.energy_function_index, self.energy_function_index
844            ),
845        };
846        let nspins = self
847            .sp_group
848            .dataset(&format!("{orbital_path}/nsets"))?
849            .read_scalar::<usize>()?;
850        let nmo = self
851            .sp_group
852            .dataset(&format!("{orbital_path}/norb",))?
853            .read_scalar::<usize>()?;
854        let (spincons, occs) = match nspins {
855            1 => {
856                log::warn!(
857                    "The number of spin spaces detected is 1. \
858                    It will be assumed that this implies an RHF calculation. \
859                    However, it must be noted that, if the calculation is GHF instead, then the \
860                    following symmetry analysis will be wrong, because Q-Chem does not archive \
861                    GHF MO coefficients correctly."
862                );
863                let nalpha = self
864                    .sp_group
865                    .dataset("structure/nalpha")?
866                    .read_scalar::<usize>()?;
867                let occ_a = Array1::from_vec(
868                    (0..nmo)
869                        .map(|i| {
870                            if i < nalpha {
871                                <T as ComplexFloat>::Real::one()
872                            } else {
873                                <T as ComplexFloat>::Real::zero()
874                            }
875                        })
876                        .collect::<Vec<_>>(),
877                );
878                (SpinConstraint::Restricted(2), vec![occ_a])
879            }
880            2 => {
881                let nalpha = self
882                    .sp_group
883                    .dataset("structure/nalpha")?
884                    .read_scalar::<usize>()?;
885                let nbeta = self
886                    .sp_group
887                    .dataset("structure/nbeta")?
888                    .read_scalar::<usize>()?;
889                let occ_a = Array1::from_vec(
890                    (0..nmo)
891                        .map(|i| {
892                            if i < nalpha {
893                                <T as ComplexFloat>::Real::one()
894                            } else {
895                                <T as ComplexFloat>::Real::zero()
896                            }
897                        })
898                        .collect::<Vec<_>>(),
899                );
900                let occ_b = Array1::from_vec(
901                    (0..nmo)
902                        .map(|i| {
903                            if i < nbeta {
904                                <T as ComplexFloat>::Real::one()
905                            } else {
906                                <T as ComplexFloat>::Real::zero()
907                            }
908                        })
909                        .collect::<Vec<_>>(),
910                );
911                (SpinConstraint::Unrestricted(2, true), vec![occ_a, occ_b])
912            }
913            _ => {
914                bail!("Unexpected number of spin spaces from Q-Chem.")
915            }
916        };
917        let cs = self
918            .sp_group
919            .dataset(&format!("{orbital_path}/mo_coefficients"))?
920            .read::<T, Ix3>()?
921            .axis_iter(Axis(0))
922            .map(|c| c.to_owned())
923            .collect::<Vec<_>>();
924        let mo_energies = self
925            .sp_group
926            .dataset(&format!("{orbital_path}/mo_energies"))
927            .and_then(|mo_energies_dataset| {
928                mo_energies_dataset.read_2d::<T>().map(|mo_energies_arr| {
929                    mo_energies_arr
930                        .columns()
931                        .into_iter()
932                        .map(|c| c.to_owned())
933                        .collect::<Vec<_>>()
934                })
935            })
936            .ok();
937
938        SlaterDeterminant::builder()
939            .structure_constraint(spincons)
940            .bao(bao)
941            .complex_symmetric(false)
942            .mol(mol)
943            .coefficients(&cs)
944            .occupations(&occs)
945            .mo_energies(mo_energies)
946            .energy(energy)
947            .threshold(threshold)
948            .build()
949            .map_err(|err| err.into())
950    }
951}
952
953// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric type f64
954// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
955#[duplicate_item(
956    [
957        gtype_ [ UnitaryRepresentedSymmetryGroup ]
958        doc_ [ "Performs symmetry-group detection and unitary-represented representation analysis." ]
959    ]
960    [
961        gtype_ [ MagneticRepresentedSymmetryGroup ]
962        doc_ [ "Performs symmetry-group detection and magnetic-represented corepresentation analysis." ]
963    ]
964)]
965impl<'a> QChemSlaterDeterminantH5SinglePointDriver<'a, gtype_, f64> {
966    #[doc = doc_]
967    pub fn analyse(&mut self) -> Result<(), anyhow::Error> {
968        let mol = self.extract_molecule()
969            .with_context(|| "Unable to extract the molecule from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")?;
970        log::debug!("Performing symmetry-group detection...");
971        let pd_res = match self.symmetry_group_detection_input {
972            SymmetryGroupDetectionInputKind::Parameters(pd_params) => {
973                let mut pd_driver = SymmetryGroupDetectionDriver::builder()
974                    .parameters(pd_params)
975                    .molecule(Some(&mol))
976                    .build()
977                    .unwrap();
978                let pd_run = pd_driver.run();
979                if let Err(err) = pd_run {
980                    qsym2_error!("Symmetry-group detection has failed with error:");
981                    qsym2_error!("  {err:#}");
982                }
983                let pd_res = pd_driver.result()?;
984                pd_res.clone()
985            }
986            SymmetryGroupDetectionInputKind::FromFile(path) => {
987                read_qsym2_binary::<SymmetryGroupDetectionResult, _>(path, QSym2FileType::Sym)
988                    .with_context(|| "Unable to read the specified .qsym2.sym file while performing symmetry analysis for a single-point Q-Chem calculation")?
989            }
990        };
991        let recentred_mol = &pd_res.pre_symmetry.recentred_molecule;
992        let sym = if self.rep_analysis_parameters.use_magnetic_group.is_some() {
993            pd_res.magnetic_symmetry.clone()
994        } else {
995            Some(pd_res.unitary_symmetry.clone())
996        }
997        .ok_or(format_err!("Symmetry not found."))?;
998        log::debug!("Performing symmetry-group detection... Done.");
999
1000        let rep = || {
1001            log::debug!("Extracting AO basis information for representation analysis...");
1002            let sao = self.recompute_sao()
1003                .with_context(|| "Unable to extract the SAO matrix from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")
1004                .map_err(|err| err.to_string())?;
1005            let bao = self.extract_bao(recentred_mol)
1006                .with_context(|| "Unable to extract the basis angular order information from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")
1007                .map_err(|err| err.to_string())?;
1008            let basis_set_opt = if self.rep_analysis_parameters.analyse_density_symmetries {
1009                self.extract_basis_set(recentred_mol).ok()
1010            } else {
1011                None
1012            };
1013            log::debug!("Extracting AO basis information for representation analysis... Done.");
1014
1015            #[cfg(feature = "integrals")]
1016            let sao_4c: Option<Array4<f64>> = basis_set_opt.map(|basis_set| {
1017                log::debug!("Computing four-centre overlap integrals for density symmetry analysis...");
1018                let stc = build_shell_tuple_collection![
1019                    <s1, s2, s3, s4>;
1020                    false, false, false, false;
1021                    &basis_set, &basis_set, &basis_set, &basis_set;
1022                    f64
1023                ];
1024                let sao_4c = stc.overlap([0, 0, 0, 0])
1025                    .pop()
1026                    .expect("Unable to retrieve the four-centre overlap tensor.");
1027                log::debug!("Computing four-centre overlap integrals for density symmetry analysis... Done.");
1028                sao_4c
1029            });
1030
1031            #[cfg(not(feature = "integrals"))]
1032            let sao_4c: Option<Array4<f64>> = None;
1033
1034            log::debug!(
1035                "Extracting canonical determinant information for representation analysis..."
1036            );
1037            let det = self.extract_determinant(
1038                recentred_mol,
1039                &bao,
1040                self.rep_analysis_parameters
1041                    .linear_independence_threshold,
1042                OrbitalType::Canonical,
1043            )
1044            .with_context(|| "Unable to extract the determinant from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")
1045            .map_err(|err| err.to_string())?;
1046            log::debug!(
1047                "Extracting canonical determinant information for representation analysis... Done."
1048            );
1049
1050            log::debug!("Running representation analysis on canonical determinant...");
1051            let mut sda_driver =
1052                SlaterDeterminantRepAnalysisDriver::<gtype_, f64, SpinConstraint>::builder()
1053                    .parameters(self.rep_analysis_parameters)
1054                    .angular_function_parameters(self.angular_function_analysis_parameters)
1055                    .determinant(&det)
1056                    .sao(&sao)
1057                    .sao_spatial_4c(sao_4c.as_ref())
1058                    .symmetry_group(&pd_res)
1059                    .build()
1060                    .with_context(|| "Unable to construct a Slater determinant representation analysis driver while performing symmetry analysis for a single-point Q-Chem calculation")
1061                    .map_err(|err| err.to_string())?;
1062            log_micsec_begin("Canonical orbital representation analysis");
1063            let sda_run = sda_driver.run();
1064            log_micsec_end("Canonical orbital representation analysis");
1065            qsym2_output!("");
1066            log::debug!("Running representation analysis on canonical determinant... Done.");
1067            if let Err(err) = sda_run {
1068                qsym2_error!("Representation analysis has failed with error:");
1069                qsym2_error!("  {err:#}");
1070            }
1071
1072            let _ = self
1073                .extract_determinant(
1074                    recentred_mol,
1075                    &bao,
1076                    self.rep_analysis_parameters.linear_independence_threshold,
1077                    OrbitalType::Localised,
1078                )
1079                .and_then(|loc_det| {
1080                    log::debug!("Running representation analysis on localised determinant...");
1081                    let mut loc_sda_driver = SlaterDeterminantRepAnalysisDriver::<
1082                        UnitaryRepresentedSymmetryGroup,
1083                        f64,
1084                        SpinConstraint,
1085                    >::builder()
1086                    .parameters(self.rep_analysis_parameters)
1087                    .angular_function_parameters(self.angular_function_analysis_parameters)
1088                    .determinant(&loc_det)
1089                    .sao(&sao)
1090                    .sao_spatial_4c(sao_4c.as_ref())
1091                    .symmetry_group(&pd_res)
1092                    .build()?;
1093                    log_micsec_begin("Localised orbital representation analysis");
1094                    let res = loc_sda_driver.run();
1095                    log_micsec_end("Localised orbital representation analysis");
1096                    qsym2_output!("");
1097                    log::debug!(
1098                        "Running representation analysis on localised determinant... Done."
1099                    );
1100                    res
1101                });
1102
1103            sda_driver
1104                .result()
1105                .map_err(|err| err.to_string())
1106                .and_then(|sda_res| sda_res.determinant_symmetry().clone())
1107        };
1108        self.result = Some((sym, rep()));
1109        Ok(())
1110    }
1111}
1112
1113// ---------------------
1114// Trait implementations
1115// ---------------------
1116
1117// ~~~~~~~~~~~~~~~~~~~
1118// Slater determinants
1119// ~~~~~~~~~~~~~~~~~~~
1120
1121// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric type f64
1122// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1123#[duplicate_item(
1124    [
1125        gtype_ [ UnitaryRepresentedSymmetryGroup ]
1126        err_ [ "No Q-Chem single-point analysis results (unitary-represented group, real determinant) found." ]
1127    ]
1128    [
1129        gtype_ [ MagneticRepresentedSymmetryGroup ]
1130        err_ [ "No Q-Chem single-point analysis results (magnetic-represented group, real determinant) found." ]
1131    ]
1132)]
1133impl<'a> QSym2Driver for QChemSlaterDeterminantH5SinglePointDriver<'a, gtype_, f64> {
1134    type Params = SlaterDeterminantRepAnalysisParams<f64>;
1135
1136    type Outcome = (
1137        Symmetry,
1138        Result<
1139            <<gtype_ as CharacterProperties>::CharTab as SubspaceDecomposable<f64>>::Decomposition,
1140            String,
1141        >,
1142    );
1143
1144    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
1145        self.result.as_ref().ok_or(format_err!(err_))
1146    }
1147
1148    fn run(&mut self) -> Result<(), anyhow::Error> {
1149        self.analyse()
1150    }
1151}