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(&self, mol: &'a Molecule) -> Result<BasisSet<f64, f64>, anyhow::Error> {
522        let shell_types = self
523            .sp_group
524            .dataset("aobasis/shell_types")?
525            .read_1d::<i32>()?;
526        let shell_to_atom_map = self
527            .sp_group
528            .dataset("aobasis/shell_to_atom_map")?
529            .read_1d::<usize>()?
530            .iter()
531            .zip(shell_types.iter())
532            .flat_map(|(&idx, shell_type)| {
533                if *shell_type == -1 {
534                    vec![idx, idx]
535                } else {
536                    vec![idx]
537                }
538            })
539            .collect::<Vec<_>>();
540
541        let primitives_per_shell = self
542            .sp_group
543            .dataset("aobasis/primitives_per_shell")?
544            .read_1d::<usize>()?;
545        let contraction_coefficients = self
546            .sp_group
547            .dataset("aobasis/contraction_coefficients")?
548            .read_1d::<f64>()?;
549        let sp_contraction_coefficients = self
550            .sp_group
551            .dataset("aobasis/sp_contraction_coefficients")?
552            .read_1d::<f64>()?;
553        let primitive_exponents = self
554            .sp_group
555            .dataset("aobasis/primitive_exponents")?
556            .read_1d::<f64>()?;
557        // `shell_coordinates` in Bohr radius.
558        let shell_coordinates = self
559            .sp_group
560            .dataset("aobasis/shell_coordinates")?
561            .read_2d::<f64>()?;
562
563        let bscs: Vec<BasisShellContraction<f64, f64>> = primitives_per_shell
564            .iter()
565            .scan(0, |end, n_prims| {
566                let start = *end;
567                *end += n_prims;
568                Some((start, *end))
569            })
570            .zip(shell_types.iter())
571            .zip(shell_coordinates.rows())
572            .flat_map(|(((start, end), shell_type), centre)| {
573                if *shell_type == 0 {
574                    // S shell
575                    let basis_shell = BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)));
576                    let primitives = primitive_exponents
577                        .slice(s![start..end])
578                        .iter()
579                        .cloned()
580                        .zip(
581                            contraction_coefficients
582                                .slice(s![start..end])
583                                .iter()
584                                .cloned(),
585                        )
586                        .collect::<Vec<_>>();
587                    let contraction = GaussianContraction { primitives };
588                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
589                    vec![BasisShellContraction {
590                        basis_shell,
591                        contraction,
592                        cart_origin,
593                        k: None,
594                    }]
595                } else if *shell_type == 1 {
596                    // P shell
597                    let basis_shell = BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)));
598                    let primitives = primitive_exponents
599                        .slice(s![start..end])
600                        .iter()
601                        .cloned()
602                        .zip(
603                            contraction_coefficients
604                                .slice(s![start..end])
605                                .iter()
606                                .cloned(),
607                        )
608                        .collect::<Vec<_>>();
609                    let contraction = GaussianContraction { primitives };
610                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
611                    vec![BasisShellContraction {
612                        basis_shell,
613                        contraction,
614                        cart_origin,
615                        k: None,
616                    }]
617                } else if *shell_type == -1 {
618                    // SP shell
619                    let basis_shell_s = BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)));
620                    let primitives_s = primitive_exponents
621                        .slice(s![start..end])
622                        .iter()
623                        .cloned()
624                        .zip(
625                            contraction_coefficients
626                                .slice(s![start..end])
627                                .iter()
628                                .cloned(),
629                        )
630                        .collect::<Vec<_>>();
631                    let contraction_s = GaussianContraction {
632                        primitives: primitives_s,
633                    };
634
635                    let basis_shell_p = BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)));
636                    let primitives_p = primitive_exponents
637                        .slice(s![start..end])
638                        .iter()
639                        .cloned()
640                        .zip(
641                            sp_contraction_coefficients
642                                .slice(s![start..end])
643                                .iter()
644                                .cloned(),
645                        )
646                        .collect::<Vec<_>>();
647                    let contraction_p = GaussianContraction {
648                        primitives: primitives_p,
649                    };
650
651                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
652                    vec![
653                        BasisShellContraction {
654                            basis_shell: basis_shell_s,
655                            contraction: contraction_s,
656                            cart_origin: cart_origin.clone(),
657                            k: None,
658                        },
659                        BasisShellContraction {
660                            basis_shell: basis_shell_p,
661                            contraction: contraction_p,
662                            cart_origin,
663                            k: None,
664                        },
665                    ]
666                } else if *shell_type < 0 {
667                    // Cartesian D shell or higher
668                    let l = shell_type.unsigned_abs();
669                    let basis_shell = BasisShell::new(l, ShellOrder::Cart(CartOrder::qchem(l)));
670                    let primitives = primitive_exponents
671                        .slice(s![start..end])
672                        .iter()
673                        .cloned()
674                        .zip(
675                            contraction_coefficients
676                                .slice(s![start..end])
677                                .iter()
678                                .cloned(),
679                        )
680                        .collect::<Vec<_>>();
681                    let contraction = GaussianContraction { primitives };
682                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
683                    vec![BasisShellContraction {
684                        basis_shell,
685                        contraction,
686                        cart_origin,
687                        k: None,
688                    }]
689                } else {
690                    // Pure D shell or higher
691                    let l = shell_type.unsigned_abs();
692                    let basis_shell =
693                        BasisShell::new(l, ShellOrder::Pure(PureOrder::increasingm(l)));
694                    let primitives = primitive_exponents
695                        .slice(s![start..end])
696                        .iter()
697                        .cloned()
698                        .zip(
699                            contraction_coefficients
700                                .slice(s![start..end])
701                                .iter()
702                                .cloned(),
703                        )
704                        .collect::<Vec<_>>();
705                    let contraction = GaussianContraction { primitives };
706                    let cart_origin = Point3::new(centre[0], centre[1], centre[2]);
707                    vec![BasisShellContraction {
708                        basis_shell,
709                        contraction,
710                        cart_origin,
711                        k: None,
712                    }]
713                }
714            })
715            .collect::<Vec<BasisShellContraction<f64, f64>>>();
716
717        let basis_atoms = mol
718            .atoms
719            .iter()
720            .enumerate()
721            .map(|(atom_i, _)| {
722                let shells = bscs
723                    .iter()
724                    .zip(shell_to_atom_map.iter())
725                    .filter_map(|(bs, atom_index)| {
726                        if *atom_index == atom_i {
727                            Some(bs.clone())
728                        } else {
729                            None
730                        }
731                    })
732                    .collect::<Vec<_>>();
733                shells
734            })
735            .collect::<Vec<Vec<BasisShellContraction<f64, f64>>>>();
736
737        let mut basis_set = BasisSet::<f64, f64>::new(basis_atoms);
738
739        // Q-Chem renormalises each Gaussian primitive, but this is not the convention used in
740        // QSym². We therefore un-normalise the Gaussian primitives to restore the original
741        // contraction coefficients.
742        let prefactor = (2.0 / std::f64::consts::PI).powf(0.75);
743        basis_set.all_shells_mut().for_each(|bsc| {
744            let l = bsc.basis_shell().l;
745            let l_i32 = l
746                .to_i32()
747                .unwrap_or_else(|| panic!("Unable to convert `{l}` to `i32`."));
748            let l_f64 = l
749                .to_f64()
750                .unwrap_or_else(|| panic!("Unable to convert `{l}` to `f64`."));
751            let doufac_sqrt = if l == 0 {
752                1
753            } else {
754                ((2 * l) - 1)
755                    .checked_double_factorial()
756                    .unwrap_or_else(|| panic!("Unable to obtain `{}!!`.", 2 * l - 1))
757            }
758            .to_f64()
759            .unwrap_or_else(|| panic!("Unable to convert `{}!!` to `f64`.", 2 * l - 1))
760            .sqrt();
761            bsc.contraction.primitives.iter_mut().for_each(|(a, c)| {
762                let n = prefactor * 2.0.powi(l_i32) * a.powf(l_f64 / 2.0 + 0.75) / doufac_sqrt;
763                *c /= n;
764            });
765        });
766        Ok(basis_set)
767    }
768
769    /// Recomputes the spatial atomic-orbital overlap matrix.
770    ///
771    /// The overlap matrix stored in the H5 group unfortunately uses lexicographic order for
772    /// Cartesian functions, which is inconsistent with that used in the coefficients. We thus
773    /// recompute the overlap matrix from the basis set information using the conventional Q-Chem
774    /// order for Cartesian functions.
775    pub fn recompute_sao(&self) -> Result<Array2<f64>, anyhow::Error> {
776        log::debug!("Recomputing atomic-orbital overlap matrix...");
777        let mol = self.extract_molecule()?;
778        let basis_set = self.extract_basis_set(&mol)?;
779        let stc = build_shell_tuple_collection![
780            <s1, s2>;
781            false, false;
782            &basis_set, &basis_set;
783            f64
784        ];
785        let sao_res = stc
786            .overlap([0, 0])
787            .pop()
788            .ok_or(format_err!("Unable to compute the AO overlap matrix."));
789        log::debug!("Recomputing atomic-orbital overlap matrix... Done.");
790        sao_res
791    }
792}
793
794// ~~~~~~~~~~~~~~~~~~~
795// Slater determinants
796// ~~~~~~~~~~~~~~~~~~~
797
798// Generic for all symmetry groups G and determinant numeric type T
799// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
800impl<'a, G, T> QChemSlaterDeterminantH5SinglePointDriver<'a, G, T>
801where
802    G: SymmetryGroupProperties + Clone,
803    G::CharTab: SubspaceDecomposable<T>,
804    T: ComplexFloat + Lapack + H5Type,
805    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
806{
807    /// Extracts the Slater determinant from the single-point H5 group.
808    ///
809    /// # Arguments
810    ///
811    /// * `mol` - The molecule to be associated with the extracted determinant.
812    /// * `bao` - The basis angular order information to be associated with the extracted determinant.
813    /// * `threshold` - The comparison threshold to be associated with the extracted determinant.
814    ///
815    /// # Returns
816    ///
817    /// The extracted Slater determinant.
818    pub fn extract_determinant(
819        &self,
820        mol: &'a Molecule,
821        bao: &'a BasisAngularOrder,
822        threshold: <T as ComplexFloat>::Real,
823        orbital_type: OrbitalType,
824    ) -> Result<SlaterDeterminant<'a, T, SpinConstraint>, anyhow::Error> {
825        let energy = self
826            .sp_group
827            .dataset(&format!(
828                "energy_function/{}/energy",
829                self.energy_function_index
830            ))?
831            .read_scalar::<T>()
832            .map_err(|err| err.to_string());
833        let orbital_path = match orbital_type {
834            OrbitalType::Canonical => format!(
835                "energy_function/{}/method/scf/molecular_orbitals",
836                self.energy_function_index
837            ),
838            OrbitalType::Localised => format!(
839                "energy_function/{}/analysis/localized_orbitals/{}/molecular_orbitals",
840                self.energy_function_index, self.energy_function_index
841            ),
842        };
843        let nspins = self
844            .sp_group
845            .dataset(&format!("{orbital_path}/nsets"))?
846            .read_scalar::<usize>()?;
847        let nmo = self
848            .sp_group
849            .dataset(&format!("{orbital_path}/norb",))?
850            .read_scalar::<usize>()?;
851        let (spincons, occs) = match nspins {
852            1 => {
853                log::warn!(
854                    "The number of spin spaces detected is 1. \
855                    It will be assumed that this implies an RHF calculation. \
856                    However, it must be noted that, if the calculation is GHF instead, then the \
857                    following symmetry analysis will be wrong, because Q-Chem does not archive \
858                    GHF MO coefficients correctly."
859                );
860                let nalpha = self
861                    .sp_group
862                    .dataset("structure/nalpha")?
863                    .read_scalar::<usize>()?;
864                let occ_a = Array1::from_vec(
865                    (0..nmo)
866                        .map(|i| {
867                            if i < nalpha {
868                                <T as ComplexFloat>::Real::one()
869                            } else {
870                                <T as ComplexFloat>::Real::zero()
871                            }
872                        })
873                        .collect::<Vec<_>>(),
874                );
875                (SpinConstraint::Restricted(2), vec![occ_a])
876            }
877            2 => {
878                let nalpha = self
879                    .sp_group
880                    .dataset("structure/nalpha")?
881                    .read_scalar::<usize>()?;
882                let nbeta = self
883                    .sp_group
884                    .dataset("structure/nbeta")?
885                    .read_scalar::<usize>()?;
886                let occ_a = Array1::from_vec(
887                    (0..nmo)
888                        .map(|i| {
889                            if i < nalpha {
890                                <T as ComplexFloat>::Real::one()
891                            } else {
892                                <T as ComplexFloat>::Real::zero()
893                            }
894                        })
895                        .collect::<Vec<_>>(),
896                );
897                let occ_b = Array1::from_vec(
898                    (0..nmo)
899                        .map(|i| {
900                            if i < nbeta {
901                                <T as ComplexFloat>::Real::one()
902                            } else {
903                                <T as ComplexFloat>::Real::zero()
904                            }
905                        })
906                        .collect::<Vec<_>>(),
907                );
908                (SpinConstraint::Unrestricted(2, true), vec![occ_a, occ_b])
909            }
910            _ => {
911                bail!("Unexpected number of spin spaces from Q-Chem.")
912            }
913        };
914        let cs = self
915            .sp_group
916            .dataset(&format!("{orbital_path}/mo_coefficients"))?
917            .read::<T, Ix3>()?
918            .axis_iter(Axis(0))
919            .map(|c| c.to_owned())
920            .collect::<Vec<_>>();
921        let mo_energies = self
922            .sp_group
923            .dataset(&format!("{orbital_path}/mo_energies"))
924            .and_then(|mo_energies_dataset| {
925                mo_energies_dataset.read_2d::<T>().map(|mo_energies_arr| {
926                    mo_energies_arr
927                        .columns()
928                        .into_iter()
929                        .map(|c| c.to_owned())
930                        .collect::<Vec<_>>()
931                })
932            })
933            .ok();
934
935        SlaterDeterminant::builder()
936            .structure_constraint(spincons)
937            .bao(bao)
938            .complex_symmetric(false)
939            .mol(mol)
940            .coefficients(&cs)
941            .occupations(&occs)
942            .mo_energies(mo_energies)
943            .energy(energy)
944            .threshold(threshold)
945            .build()
946            .map_err(|err| err.into())
947    }
948}
949
950// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric type f64
951// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
952#[duplicate_item(
953    [
954        gtype_ [ UnitaryRepresentedSymmetryGroup ]
955        doc_ [ "Performs symmetry-group detection and unitary-represented representation analysis." ]
956    ]
957    [
958        gtype_ [ MagneticRepresentedSymmetryGroup ]
959        doc_ [ "Performs symmetry-group detection and magnetic-represented corepresentation analysis." ]
960    ]
961)]
962impl<'a> QChemSlaterDeterminantH5SinglePointDriver<'a, gtype_, f64> {
963    #[doc = doc_]
964    pub fn analyse(&mut self) -> Result<(), anyhow::Error> {
965        let mol = self.extract_molecule()
966            .with_context(|| "Unable to extract the molecule from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")?;
967        log::debug!("Performing symmetry-group detection...");
968        let pd_res = match self.symmetry_group_detection_input {
969            SymmetryGroupDetectionInputKind::Parameters(pd_params) => {
970                let mut pd_driver = SymmetryGroupDetectionDriver::builder()
971                    .parameters(pd_params)
972                    .molecule(Some(&mol))
973                    .build()
974                    .unwrap();
975                let pd_run = pd_driver.run();
976                if let Err(err) = pd_run {
977                    qsym2_error!("Symmetry-group detection has failed with error:");
978                    qsym2_error!("  {err:#}");
979                }
980                let pd_res = pd_driver.result()?;
981                pd_res.clone()
982            }
983            SymmetryGroupDetectionInputKind::FromFile(path) => {
984                read_qsym2_binary::<SymmetryGroupDetectionResult, _>(path, QSym2FileType::Sym)
985                    .with_context(|| "Unable to read the specified .qsym2.sym file while performing symmetry analysis for a single-point Q-Chem calculation")?
986            }
987        };
988        let recentred_mol = &pd_res.pre_symmetry.recentred_molecule;
989        let sym = if self.rep_analysis_parameters.use_magnetic_group.is_some() {
990            pd_res.magnetic_symmetry.clone()
991        } else {
992            Some(pd_res.unitary_symmetry.clone())
993        }
994        .ok_or(format_err!("Symmetry not found."))?;
995        log::debug!("Performing symmetry-group detection... Done.");
996
997        let rep = || {
998            log::debug!("Extracting AO basis information for representation analysis...");
999            let sao = self.recompute_sao()
1000                .with_context(|| "Unable to extract the SAO matrix from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")
1001                .map_err(|err| err.to_string())?;
1002            let bao = self.extract_bao(recentred_mol)
1003                .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")
1004                .map_err(|err| err.to_string())?;
1005            let basis_set_opt = if self.rep_analysis_parameters.analyse_density_symmetries {
1006                self.extract_basis_set(recentred_mol).ok()
1007            } else {
1008                None
1009            };
1010            log::debug!("Extracting AO basis information for representation analysis... Done.");
1011
1012            #[cfg(feature = "integrals")]
1013            let sao_4c: Option<Array4<f64>> = basis_set_opt.map(|basis_set| {
1014                log::debug!("Computing four-centre overlap integrals for density symmetry analysis...");
1015                let stc = build_shell_tuple_collection![
1016                    <s1, s2, s3, s4>;
1017                    false, false, false, false;
1018                    &basis_set, &basis_set, &basis_set, &basis_set;
1019                    f64
1020                ];
1021                let sao_4c = stc.overlap([0, 0, 0, 0])
1022                    .pop()
1023                    .expect("Unable to retrieve the four-centre overlap tensor.");
1024                log::debug!("Computing four-centre overlap integrals for density symmetry analysis... Done.");
1025                sao_4c
1026            });
1027
1028            #[cfg(not(feature = "integrals"))]
1029            let sao_4c: Option<Array4<f64>> = None;
1030
1031            log::debug!(
1032                "Extracting canonical determinant information for representation analysis..."
1033            );
1034            let det = self.extract_determinant(
1035                recentred_mol,
1036                &bao,
1037                self.rep_analysis_parameters
1038                    .linear_independence_threshold,
1039                OrbitalType::Canonical,
1040            )
1041            .with_context(|| "Unable to extract the determinant from the HDF5 file while performing symmetry analysis for a single-point Q-Chem calculation")
1042            .map_err(|err| err.to_string())?;
1043            log::debug!(
1044                "Extracting canonical determinant information for representation analysis... Done."
1045            );
1046
1047            log::debug!("Running representation analysis on canonical determinant...");
1048            let mut sda_driver =
1049                SlaterDeterminantRepAnalysisDriver::<gtype_, f64, SpinConstraint>::builder()
1050                    .parameters(self.rep_analysis_parameters)
1051                    .angular_function_parameters(self.angular_function_analysis_parameters)
1052                    .determinant(&det)
1053                    .sao(&sao)
1054                    .sao_spatial_4c(sao_4c.as_ref())
1055                    .symmetry_group(&pd_res)
1056                    .build()
1057                    .with_context(|| "Unable to construct a Slater determinant representation analysis driver while performing symmetry analysis for a single-point Q-Chem calculation")
1058                    .map_err(|err| err.to_string())?;
1059            log_micsec_begin("Canonical orbital representation analysis");
1060            let sda_run = sda_driver.run();
1061            log_micsec_end("Canonical orbital representation analysis");
1062            qsym2_output!("");
1063            log::debug!("Running representation analysis on canonical determinant... Done.");
1064            if let Err(err) = sda_run {
1065                qsym2_error!("Representation analysis has failed with error:");
1066                qsym2_error!("  {err:#}");
1067            }
1068
1069            let _ = self
1070                .extract_determinant(
1071                    recentred_mol,
1072                    &bao,
1073                    self.rep_analysis_parameters.linear_independence_threshold,
1074                    OrbitalType::Localised,
1075                )
1076                .and_then(|loc_det| {
1077                    log::debug!("Running representation analysis on localised determinant...");
1078                    let mut loc_sda_driver = SlaterDeterminantRepAnalysisDriver::<
1079                        UnitaryRepresentedSymmetryGroup,
1080                        f64,
1081                        SpinConstraint,
1082                    >::builder()
1083                    .parameters(self.rep_analysis_parameters)
1084                    .angular_function_parameters(self.angular_function_analysis_parameters)
1085                    .determinant(&loc_det)
1086                    .sao(&sao)
1087                    .sao_spatial_4c(sao_4c.as_ref())
1088                    .symmetry_group(&pd_res)
1089                    .build()?;
1090                    log_micsec_begin("Localised orbital representation analysis");
1091                    let res = loc_sda_driver.run();
1092                    log_micsec_end("Localised orbital representation analysis");
1093                    qsym2_output!("");
1094                    log::debug!(
1095                        "Running representation analysis on localised determinant... Done."
1096                    );
1097                    res
1098                });
1099
1100            sda_driver
1101                .result()
1102                .map_err(|err| err.to_string())
1103                .and_then(|sda_res| sda_res.determinant_symmetry().clone())
1104        };
1105        self.result = Some((sym, rep()));
1106        Ok(())
1107    }
1108}
1109
1110// ---------------------
1111// Trait implementations
1112// ---------------------
1113
1114// ~~~~~~~~~~~~~~~~~~~
1115// Slater determinants
1116// ~~~~~~~~~~~~~~~~~~~
1117
1118// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric type f64
1119// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
1120#[duplicate_item(
1121    [
1122        gtype_ [ UnitaryRepresentedSymmetryGroup ]
1123        err_ [ "No Q-Chem single-point analysis results (unitary-represented group, real determinant) found." ]
1124    ]
1125    [
1126        gtype_ [ MagneticRepresentedSymmetryGroup ]
1127        err_ [ "No Q-Chem single-point analysis results (magnetic-represented group, real determinant) found." ]
1128    ]
1129)]
1130impl<'a> QSym2Driver for QChemSlaterDeterminantH5SinglePointDriver<'a, gtype_, f64> {
1131    type Params = SlaterDeterminantRepAnalysisParams<f64>;
1132
1133    type Outcome = (
1134        Symmetry,
1135        Result<
1136            <<gtype_ as CharacterProperties>::CharTab as SubspaceDecomposable<f64>>::Decomposition,
1137            String,
1138        >,
1139    );
1140
1141    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
1142        self.result.as_ref().ok_or(format_err!(err_))
1143    }
1144
1145    fn run(&mut self) -> Result<(), anyhow::Error> {
1146        self.analyse()
1147    }
1148}