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(&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 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 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 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 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 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 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 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 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
794impl<'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 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#[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#[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}