qsym2/drivers/representation_analysis/
angular_function.rs1use anyhow::{self, ensure, format_err};
4use derive_builder::Builder;
5use nalgebra::{Point3, Vector3};
6use ndarray::{Array1, Array2};
7use num_complex::{Complex, ComplexFloat};
8use rayon::prelude::*;
9
10use crate::analysis::{EigenvalueComparisonMode, RepAnalysis};
11use crate::angmom::sh_conversion::sh_cart2rl_mat;
12use crate::angmom::spinor_rotation_3d::{SpinConstraint, SpinOrbitCoupled};
13use crate::auxiliary::atom::{Atom, ElementMap};
14use crate::auxiliary::molecule::Molecule;
15use crate::basis::ao::{
16 cart_tuple_to_str, BasisAngularOrder, BasisAtom, BasisShell, CartOrder, PureOrder, ShellOrder,
17 SpinorOrder,
18};
19use crate::chartab::chartab_group::CharacterProperties;
20use crate::chartab::SubspaceDecomposable;
21use crate::io::format::{log_subtitle, qsym2_output};
22use crate::symmetry::symmetry_group::SymmetryGroupProperties;
23use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
24use crate::target::determinant::SlaterDeterminant;
25use crate::target::orbital::orbital_analysis::generate_det_mo_orbits;
26use crate::target::tensor::axialvector::axialvector_analysis::AxialVector3SymmetryOrbit;
27use crate::target::tensor::axialvector::{AxialVector3, TimeParity};
28
29#[derive(Clone, Builder, Debug)]
39pub struct AngularFunctionRepAnalysisParams {
40 #[builder(default = "1e-7")]
42 pub integrality_threshold: f64,
43
44 #[builder(default = "1e-7")]
46 pub linear_independence_threshold: f64,
47
48 #[builder(default = "2")]
50 pub max_angular_momentum: u32,
51}
52
53impl AngularFunctionRepAnalysisParams {
54 pub fn builder() -> AngularFunctionRepAnalysisParamsBuilder {
56 AngularFunctionRepAnalysisParamsBuilder::default()
57 }
58}
59
60impl Default for AngularFunctionRepAnalysisParams {
61 fn default() -> Self {
62 Self::builder()
63 .build()
64 .expect("Unable to build as default `AngularFunctionRepAnalysisParams`.")
65 }
66}
67
68pub(crate) fn find_angular_function_representation<G>(
81 group: &G,
82 params: &AngularFunctionRepAnalysisParams,
83) -> Result<(), anyhow::Error>
84where
85 G: SymmetryGroupProperties + Clone + Send + Sync,
86 G::CharTab: SubspaceDecomposable<f64>,
87 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<f64>>::Decomposition: Send + Sync,
88{
89 let emap = ElementMap::new();
90 let thresh = group
91 .elements()
92 .clone()
93 .into_iter()
94 .next()
95 .ok_or_else(|| format_err!("No symmetry operation found."))?
96 .generating_element
97 .threshold();
98 let mol = Molecule::from_atoms(
99 &[Atom::new_ordinary("H", Point3::origin(), &emap, thresh)],
100 thresh,
101 );
102 let lmax = params.max_angular_momentum;
103
104 let (pure_symss, cart_symss) = (0..=lmax).fold(
105 (
106 Vec::with_capacity(usize::try_from(lmax)?),
107 Vec::with_capacity(usize::try_from(lmax)?),
108 ),
109 |mut acc, l| {
110 [
111 ShellOrder::Pure(PureOrder::increasingm(l)),
112 ShellOrder::Cart(CartOrder::lex(l)),
113 ]
114 .iter()
115 .for_each(|shell_order| {
116 let bao = BasisAngularOrder::new(&[BasisAtom::new(
117 &mol.atoms[0],
118 &[BasisShell::new(l, shell_order.clone())],
119 )]);
120 let nbas = bao.n_funcs();
121 let cs = vec![Array2::<f64>::eye(nbas)];
122 let occs = vec![Array1::<f64>::ones(nbas)];
123 let sao = match shell_order {
124 ShellOrder::Pure(_) | ShellOrder::Spinor(_) => {
125 Array2::<f64>::eye(bao.n_funcs())
126 }
127 ShellOrder::Cart(cartorder) => {
128 let cart2rl =
129 sh_cart2rl_mat(l, l, cartorder, true, &PureOrder::increasingm(l));
130 cart2rl.mapv(ComplexFloat::conj).t().dot(&cart2rl)
131 }
132 };
133
134 let mo_symmetries = SlaterDeterminant::<f64, SpinConstraint>::builder()
135 .structure_constraint(SpinConstraint::Restricted(1))
136 .bao(&bao)
137 .complex_symmetric(false)
138 .mol(&mol)
139 .coefficients(&cs)
140 .occupations(&occs)
141 .threshold(params.linear_independence_threshold)
142 .build()
143 .map_err(|err| format_err!(err))
144 .and_then(|det| {
145 let mos = det.to_orbitals();
146 generate_det_mo_orbits(
147 &det,
148 &mos,
149 group,
150 &sao,
151 None, params.integrality_threshold,
153 params.linear_independence_threshold,
154 SymmetryTransformationKind::Spatial,
155 EigenvalueComparisonMode::Real,
156 true,
157 )
158 .map(|(_, mut mo_orbitss)| {
159 mo_orbitss[0]
160 .par_iter_mut()
161 .map(|mo_orbit| {
162 mo_orbit.calc_xmat(false)?;
163 mo_orbit.analyse_rep().map_err(|err| format_err!(err))
164 })
165 .collect::<Vec<_>>()
166 })
167 });
168 match shell_order {
169 ShellOrder::Pure(_) => acc.0.push(mo_symmetries),
170 ShellOrder::Cart(_) => acc.1.push(mo_symmetries),
171 ShellOrder::Spinor(_) => todo!(),
172 }
173 });
174 acc
175 },
176 );
177
178 let pure_sym_strss = pure_symss
179 .into_iter()
180 .map(|l_pure_syms_opt| {
181 l_pure_syms_opt
182 .map(|l_pure_syms| {
183 l_pure_syms
184 .into_iter()
185 .map(|sym_opt| {
186 sym_opt
187 .map(|sym| sym.to_string())
188 .unwrap_or_else(|err| err.to_string())
189 })
190 .collect::<Vec<_>>()
191 })
192 .unwrap_or_else(|err| vec![err.to_string()])
193 })
194 .collect::<Vec<_>>();
195 let cart_sym_strss = cart_symss
196 .into_iter()
197 .map(|l_cart_syms_opt| {
198 l_cart_syms_opt
199 .map(|l_cart_syms| {
200 l_cart_syms
201 .into_iter()
202 .map(|sym_opt| {
203 sym_opt
204 .map(|sym| sym.to_string())
205 .unwrap_or_else(|err| err.to_string())
206 })
207 .collect::<Vec<_>>()
208 })
209 .unwrap_or_else(|err| vec![err.to_string()])
210 })
211 .collect::<Vec<_>>();
212 let rot_sym_strs = [Vector3::x(), Vector3::y(), Vector3::z()]
213 .into_iter()
214 .map(|v| {
215 let rv = AxialVector3::<f64>::builder()
216 .components(v)
217 .time_parity(TimeParity::Even)
218 .threshold(params.linear_independence_threshold)
219 .build()
220 .map_err(|err| format_err!(err))?;
221 let mut orbit_rv = AxialVector3SymmetryOrbit::builder()
222 .group(group)
223 .origin(&rv)
224 .integrality_threshold(params.integrality_threshold)
225 .linear_independence_threshold(params.linear_independence_threshold)
226 .symmetry_transformation_kind(SymmetryTransformationKind::Spatial)
227 .eigenvalue_comparison_mode(EigenvalueComparisonMode::Real)
228 .build()
229 .map_err(|err| format_err!(err))?;
230 let _ = orbit_rv
231 .calc_smat(None, None, true)
232 .unwrap()
233 .calc_xmat(false);
234 orbit_rv
235 .analyse_rep()
236 .map(|sym| sym.to_string())
237 .map_err(|err| format_err!(err))
238 })
239 .map(|sym| sym.unwrap_or_else(|err| err.to_string()))
240 .collect::<Vec<_>>();
241
242 let pure_sym_width = pure_sym_strss
243 .iter()
244 .flat_map(|syms| syms.iter().map(|sym| sym.chars().count()))
245 .max()
246 .unwrap_or(9)
247 .max(9);
248 let cart_sym_width = cart_sym_strss
249 .iter()
250 .flat_map(|syms| syms.iter().map(|sym| sym.chars().count()))
251 .max()
252 .unwrap_or(9)
253 .max(9);
254 let rot_sym_width = rot_sym_strs
255 .iter()
256 .map(|sym| sym.chars().count())
257 .max()
258 .unwrap_or(8)
259 .max(8);
260
261 let l_width = lmax.to_string().chars().count();
262 let pure_width = (l_width + 1).max(4);
263 let cart_width = usize::try_from(lmax)?.max(4);
264 let rot_width = 3;
265
266 log_subtitle(&format!(
267 "Space-fixed spatial angular function symmetries in {}",
268 group.finite_subgroup_name().unwrap_or(&group.name())
269 ));
270 qsym2_output!("");
271 if lmax < 1 {
272 qsym2_output!(
273 "{}",
274 "┈".repeat(l_width + pure_width + pure_sym_width + cart_width + cart_sym_width + 11)
275 );
276 qsym2_output!(
277 " {:>l_width$} {:>pure_width$} {:<pure_sym_width$} {:>cart_width$} {:<}",
278 "l",
279 "Pure",
280 "Pure sym.",
281 "Cart",
282 "Cart sym."
283 );
284 qsym2_output!(
285 "{}",
286 "┈".repeat(l_width + pure_width + pure_sym_width + cart_width + cart_sym_width + 11)
287 );
288 } else {
289 qsym2_output!(
290 "{}",
291 "┈".repeat(
292 l_width
293 + pure_width
294 + pure_sym_width
295 + cart_width
296 + cart_sym_width
297 + rot_width
298 + rot_sym_width
299 + 15
300 )
301 );
302 qsym2_output!(
303 " {:>l_width$} {:>pure_width$} {:<pure_sym_width$} {:>cart_width$} {:<cart_sym_width$} {:>rot_width$} {:<}",
304 "l",
305 "Pure",
306 "Pure sym.",
307 "Cart",
308 "Cart sym.",
309 "Rot",
310 "Rot sym."
311 );
312 qsym2_output!(
313 "{}",
314 "┈".repeat(
315 l_width
316 + pure_width
317 + pure_sym_width
318 + cart_width
319 + cart_sym_width
320 + rot_width
321 + rot_sym_width
322 + 15
323 )
324 );
325 }
326
327 let empty_str = String::new();
328 (0..=usize::try_from(lmax)?).for_each(|l| {
329 if l > 0 {
330 qsym2_output!("");
331 }
332 let n_pure = 2 * l + 1;
333 let mut i_pure = 0;
334 let mut i_rot = 0;
335
336 let l_u32 = u32::try_from(l).unwrap_or_else(|err| panic!("{err}"));
337 let cartorder = CartOrder::lex(l_u32);
338 cartorder
339 .iter()
340 .enumerate()
341 .for_each(|(i_cart, cart_tuple)| {
342 let l_str = if i_cart == 0 {
343 format!("{l:>l_width$}")
344 } else {
345 " ".repeat(l_width)
346 };
347
348 let pure_str = if i_pure < n_pure {
349 let pure_str = format!(
350 "{:>pure_width$} {:<pure_sym_width$}",
351 if i_pure < l {
352 format!("-{}", i_pure.abs_diff(l))
353 } else {
354 format!("+{}", i_pure.abs_diff(l))
355 },
356 pure_sym_strss
357 .get(l)
358 .and_then(|l_pure_sym_strs| l_pure_sym_strs.get(i_pure))
359 .unwrap_or(&empty_str)
360 );
361 i_pure += 1;
362 pure_str
363 } else {
364 " ".repeat(pure_width + pure_sym_width + 2)
365 };
366
367 let cart_symbol = cart_tuple_to_str(cart_tuple, true);
368 let cart_str = if l == 1 {
369 format!(
371 "{cart_symbol:>cart_width$} {:<cart_sym_width$}",
372 cart_sym_strss
373 .get(l)
374 .and_then(|l_cart_sym_strs| l_cart_sym_strs.get(i_cart))
375 .unwrap_or(&empty_str)
376 )
377 } else {
378 format!(
380 "{cart_symbol:>cart_width$} {:<}",
381 cart_sym_strss
382 .get(l)
383 .and_then(|l_cart_sym_strs| l_cart_sym_strs.get(i_cart))
384 .unwrap_or(&empty_str)
385 )
386 };
387
388 if l == 1 && i_rot < 3 {
389 let rot_str = format!(
390 "{:>rot_width$} {:<}",
391 match i_rot {
392 0 => "Rx",
393 1 => "Ry",
394 2 => "Rz",
395 _ => "--",
396 },
397 rot_sym_strs.get(i_rot).unwrap_or(&"--".to_string())
398 );
399 i_rot += 1;
400 qsym2_output!(" {l_str} {pure_str} {cart_str} {rot_str}");
401 } else {
402 qsym2_output!(" {l_str} {pure_str} {cart_str}");
403 }
404 });
405 });
406 if lmax < 1 {
407 qsym2_output!(
408 "{}",
409 "┈".repeat(l_width + pure_width + pure_sym_width + cart_width + cart_sym_width + 11)
410 );
411 } else {
412 qsym2_output!(
413 "{}",
414 "┈".repeat(
415 l_width
416 + pure_width
417 + pure_sym_width
418 + cart_width
419 + cart_sym_width
420 + rot_width
421 + rot_sym_width
422 + 15
423 )
424 );
425 }
426 qsym2_output!("");
427
428 Ok(())
429}
430
431pub(crate) fn find_spinor_function_representation<G>(
439 group: &G,
440 params: &AngularFunctionRepAnalysisParams,
441) -> Result<(), anyhow::Error>
442where
443 G: SymmetryGroupProperties + Clone + Send + Sync,
444 G::CharTab: SubspaceDecomposable<Complex<f64>>,
445 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<Complex<f64>>>::Decomposition:
446 Send + Sync,
447{
448 ensure!(
449 group.is_double_group(),
450 "The specified group is not a double group."
451 );
452
453 let emap = ElementMap::new();
454 let thresh = group
455 .elements()
456 .clone()
457 .into_iter()
458 .next()
459 .ok_or_else(|| format_err!("No symmetry operation found."))?
460 .generating_element
461 .threshold();
462 let mol = Molecule::from_atoms(
463 &[Atom::new_ordinary("H", Point3::origin(), &emap, thresh)],
464 thresh,
465 );
466 let lmax = params.max_angular_momentum;
467
468 let spinor_symss = (1..2 * lmax).step_by(2).fold(
469 Vec::with_capacity(usize::try_from(lmax)?),
470 |mut acc, two_j| {
471 let shell_order = ShellOrder::Spinor(SpinorOrder::increasingm(two_j, true));
472 let bao = BasisAngularOrder::new(&[BasisAtom::new(
473 &mol.atoms[0],
474 &[BasisShell::new(two_j, shell_order.clone())],
475 )]);
476 let nbas = bao.n_funcs();
477 let cs = vec![Array2::<Complex<f64>>::eye(nbas)];
478 let occs = vec![Array1::<f64>::ones(nbas)];
479 let sao = Array2::<Complex<f64>>::eye(bao.n_funcs());
480
481 let mo_symmetries = SlaterDeterminant::<Complex<f64>, SpinOrbitCoupled>::builder()
482 .structure_constraint(SpinOrbitCoupled::JAdapted(1))
483 .bao(&bao)
484 .complex_symmetric(false)
485 .mol(&mol)
486 .coefficients(&cs)
487 .occupations(&occs)
488 .threshold(params.linear_independence_threshold)
489 .build()
490 .map_err(|err| format_err!(err))
491 .and_then(|det| {
492 let mos = det.to_orbitals();
493 generate_det_mo_orbits(
494 &det,
495 &mos,
496 group,
497 &sao,
498 None, params.integrality_threshold,
500 params.linear_independence_threshold,
501 SymmetryTransformationKind::SpinSpatial,
502 EigenvalueComparisonMode::Modulus,
503 true,
504 )
505 .map(|(_, mut mo_orbitss)| {
506 mo_orbitss[0]
507 .par_iter_mut()
508 .map(|mo_orbit| {
509 mo_orbit.calc_xmat(false)?;
510 mo_orbit.analyse_rep().map_err(|err| format_err!(err))
511 })
512 .collect::<Vec<_>>()
513 })
514 });
515 acc.push(mo_symmetries);
516 acc
517 },
518 );
519
520 let spinor_sym_strss = spinor_symss
521 .into_iter()
522 .map(|l_spinor_syms_opt| {
523 l_spinor_syms_opt
524 .map(|l_spinor_syms| {
525 l_spinor_syms
526 .into_iter()
527 .map(|sym_opt| {
528 sym_opt
529 .map(|sym| sym.to_string())
530 .unwrap_or_else(|err| err.to_string())
531 })
532 .collect::<Vec<_>>()
533 })
534 .unwrap_or_else(|err| vec![err.to_string()])
535 })
536 .collect::<Vec<_>>();
537
538 let spinor_sym_width = spinor_sym_strss
539 .iter()
540 .flat_map(|syms| syms.iter().map(|sym| sym.chars().count()))
541 .max()
542 .unwrap_or(11)
543 .max(11);
544
545 let j_width = (2 * lmax).to_string().chars().count() + 2;
546 let mj_width = (j_width + 1).max(6);
547
548 log_subtitle(&format!(
549 "Space-fixed spinor function symmetries in {}",
550 group.finite_subgroup_name().unwrap_or(&group.name())
551 ));
552 qsym2_output!("");
553 qsym2_output!("{}", "┈".repeat(j_width + mj_width + spinor_sym_width + 6));
554 qsym2_output!(
555 " {:>j_width$} {:>mj_width$} {:<}",
556 "j",
557 "Spinor",
558 "Spinor sym.",
559 );
560 qsym2_output!("{}", "┈".repeat(j_width + mj_width + spinor_sym_width + 6));
561
562 let empty_str = String::new();
563 (1..usize::try_from(2 * lmax)?)
564 .step_by(2)
565 .enumerate()
566 .for_each(|(i_two_j, two_j)| {
567 if two_j > 1 {
568 qsym2_output!("");
569 }
570 let n_spinor = two_j + 1;
571
572 let two_j_u32 = u32::try_from(two_j).unwrap_or_else(|err| panic!("{err}"));
573 let spinororder = SpinorOrder::increasingm(two_j_u32, true);
574 spinororder
575 .iter()
576 .enumerate()
577 .for_each(|(i_spinor, two_mj)| {
578 let j_str = if i_spinor == 0 {
579 let j_str_temp = format!("{two_j}/2");
580 format!("{j_str_temp:>j_width$}")
581 } else {
582 " ".repeat(j_width)
583 };
584
585 let spinor_str = if i_spinor < n_spinor {
586 let spinor_str_temp = format!("{two_mj:+}/2");
587 let spinor_str = format!(
588 "{spinor_str_temp:>mj_width$} {:<}",
589 spinor_sym_strss
590 .get(i_two_j)
591 .and_then(|l_spinor_sym_strs| l_spinor_sym_strs.get(i_spinor))
592 .unwrap_or(&empty_str)
593 );
594 spinor_str
595 } else {
596 " ".repeat(mj_width + spinor_sym_width + 2)
597 };
598
599 qsym2_output!(" {j_str} {spinor_str}");
600 });
601 });
602 qsym2_output!("{}", "┈".repeat(j_width + mj_width + spinor_sym_width + 6));
603 qsym2_output!("");
604
605 Ok(())
606}