qsym2/drivers/representation_analysis/
angular_function.rs

1//! Driver for symmetry analysis of angular functions.
2
3use 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// ==================
30// Struct definitions
31// ==================
32
33// ----------
34// Parameters
35// ----------
36
37/// Structure containing control parameters for angular function representation analysis.
38#[derive(Clone, Builder, Debug)]
39pub struct AngularFunctionRepAnalysisParams {
40    /// Threshold for checking if subspace multiplicities are integral.
41    #[builder(default = "1e-7")]
42    pub integrality_threshold: f64,
43
44    /// Threshold for determining zero eigenvalues in the orbit overlap matrix.
45    #[builder(default = "1e-7")]
46    pub linear_independence_threshold: f64,
47
48    /// The maximum angular momentum degree to be analysed.
49    #[builder(default = "2")]
50    pub max_angular_momentum: u32,
51}
52
53impl AngularFunctionRepAnalysisParams {
54    /// Returns a builder to construct a [`AngularFunctionRepAnalysisParams`] structure.
55    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
68// =========
69// Functions
70// =========
71
72/// Determines the (co)representations of a group spanned by the angular functions (spherical/solid
73/// harmonics or Cartesian functions).
74///
75/// # Arguments
76///
77/// * `group` - A symmetry group.
78/// * `params` - A parameter structure controlling the determination of angular function
79/// symmetries.
80pub(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, // Is this right for complex spherical harmonics?
152                            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                    // Rot sym to follow.
370                    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                    // No rot sym to follow.
379                    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
431/// Determines the (co)representations of a group spanned by the spinor functions.
432///
433/// # Arguments
434///
435/// * `group` - A symmetry group.
436/// * `params` - A parameter structure controlling the determination of spinor function
437/// symmetries.
438pub(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, // Is this right for complex spherical harmonics?
499                        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}