1use 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
56lazy_static! {
61 static ref SP_PATH_RE: Regex =
62 Regex::new(r"(.*sp)\\energy_function$").expect("Regex pattern invalid.");
63}
64
65#[derive(Clone, Builder)]
73pub struct QChemSlaterDeterminantH5Driver<'a, T>
74where
75 T: Clone,
76{
77 filename: PathBuf,
79
80 symmetry_group_detection_input: &'a SymmetryGroupDetectionInputKind,
82
83 angular_function_analysis_parameters: &'a AngularFunctionRepAnalysisParams,
85
86 rep_analysis_parameters: &'a SlaterDeterminantRepAnalysisParams<f64>,
88
89 #[builder(default = "None")]
93 result: Option<Vec<(String, String)>>,
94
95 #[builder(setter(skip), default = "PhantomData")]
97 numerical_type: PhantomData<T>,
98}
99
100impl<'a, T> QChemSlaterDeterminantH5Driver<'a, T>
105where
106 T: Clone,
107{
108 pub fn builder() -> QChemSlaterDeterminantH5DriverBuilder<'a, T> {
110 QChemSlaterDeterminantH5DriverBuilder::default()
111 }
112}
113
114impl<'a> QChemSlaterDeterminantH5Driver<'a, f64> {
121 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
316pub enum OrbitalType {
326 Canonical,
328
329 Localised,
331}
332
333#[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 sp_group: &'a hdf5::Group,
349
350 #[builder(setter(custom))]
353 energy_function_index: String,
354
355 symmetry_group_detection_input: &'a SymmetryGroupDetectionInputKind,
357
358 angular_function_analysis_parameters: &'a AngularFunctionRepAnalysisParams,
360
361 rep_analysis_parameters: &'a SlaterDeterminantRepAnalysisParams<f64>,
363
364 #[builder(default = "None")]
366 result: Option<(
367 Symmetry,
368 Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>,
369 )>,
370}
371
372impl<'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 pub fn builder() -> QChemSlaterDeterminantH5SinglePointDriverBuilder<'a, G, T> {
398 QChemSlaterDeterminantH5SinglePointDriverBuilder::default()
399 }
400
401 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 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 vec![BasisShell::new(0, ShellOrder::Cart(CartOrder::qchem(0)))]
461 } else if *shell_type == 1 {
462 vec![BasisShell::new(1, ShellOrder::Cart(CartOrder::qchem(1)))]
464 } else if *shell_type == -1 {
465 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 let l = shell_type.unsigned_abs();
473 vec![BasisShell::new(l, ShellOrder::Cart(CartOrder::qchem(l)))]
474 } else {
475 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 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 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 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 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 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 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 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 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 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 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
797impl<'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 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#[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#[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}