qsym2/target/noci/
basis_metric.rs

1//! Traits for defining metrics in bases.
2
3use std::fmt;
4
5use anyhow::{self, format_err};
6use itertools::{Itertools, iproduct, izip};
7use ndarray::{Array2, Array3, Axis, Ix2, s};
8use ndarray_linalg::solve::Determinant;
9use ndarray_linalg::types::Lapack;
10use num_complex::ComplexFloat;
11use rayon::prelude::*;
12
13use crate::analysis::{Orbit, Overlap};
14use crate::angmom::spinor_rotation_3d::StructureConstraint;
15use crate::symmetry::symmetry_element::SpecialSymmetryTransformation;
16use crate::symmetry::symmetry_group::SymmetryGroupProperties;
17use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
18use crate::target::noci::basis::basis_orbit::BasisSymmetryOrbit;
19use crate::target::noci::basis::{Basis, EagerBasis, FCIBasis, OrbitBasis};
20
21// ========================
22// BasisSymmetryOrbitMetric
23// ========================
24
25// -----------------
26// Trait definitions
27// -----------------
28
29/// Trait specifying the symmetry-orbit metric for bases of Slater determinants.
30///
31/// For a group $`\mathcal{G}`$, this metric is defined as
32///
33/// ```math
34/// \braket{\hat{g} \mathbf{w} | \mathbf{x}}
35/// ```
36///
37/// where $`g \in \mathcal{G}`$, and $`\mathbf{w}`$ and $`\mathbf{x}`$ are members of the basis.
38pub trait BasisSymmetryOrbitMetric<G, T>: Basis
39where
40    G: SymmetryGroupProperties,
41    Self: SymmetryTransformable,
42{
43    /// Computes the basis symmetry-orbit metric for this basis.
44    ///
45    /// For a group $`\mathcal{G}`$, this metric is defined as
46    ///
47    /// ```math
48    /// \braket{\hat{g} \mathbf{w} | \mathbf{x}}
49    /// ```
50    ///
51    /// where $`g \in \mathcal{G}`$, and $`\mathbf{w}`$ and $`\mathbf{x}`$ are members of the basis.
52    ///
53    /// # Arguments
54    ///
55    /// * `group` - The group $`\mathcal{G}`$.
56    /// * `symmetry_transformation_kind` - Enum specifying how the elements of $`\mathcal{G}`$ act
57    ///   on the elements of the basis.
58    /// * `metric` - The atomic-orbital overlap matrix with respect to the conventional sesquilinear
59    ///   inner product.
60    /// * `metric_h` - The atomic-orbital overlap matrix with respect to the bilinear inner product.
61    ///
62    /// # Returns
63    ///
64    /// The basis symmetry-orbit metric as a rank-3 array indexed by $`g`$, $`\mathbf{w}`$, and
65    /// $`\mathbf{x}`$.
66    fn basis_symmetry_orbit_metric(
67        &self,
68        group: &G,
69        symmetry_transformation_kind: &SymmetryTransformationKind,
70        metric: Option<&Array2<T>>,
71        metric_h: Option<&Array2<T>>,
72    ) -> Result<Array3<T>, anyhow::Error>;
73}
74
75// ---------------------
76// Trait implementations
77// ---------------------
78
79// ~~~~~~~~~~
80// OrbitBasis
81// ~~~~~~~~~~
82
83impl<'g, G, I, T> BasisSymmetryOrbitMetric<G, T> for OrbitBasis<'g, G, I>
84where
85    G: SymmetryGroupProperties + Clone,
86    I: Overlap<T, Ix2> + Clone,
87    T: ComplexFloat + fmt::Debug + Lapack,
88    OrbitBasis<'g, G, I>: SymmetryTransformable,
89{
90    /// Computes the basis symmetry-orbit metric for this orbit basis.
91    ///
92    /// Since this orbit basis already contains information about its group and group action,
93    /// `_group` and `_symmetry_transformation_kind` will be ignored.
94    ///
95    /// Given that the associated group of the orbit basis is $`\mathcal{G}`$ and that its origins
96    /// are $`\{ \mathbf{w}_{I} \}`$, the metric is defined as
97    ///
98    /// ```math
99    /// \braket{\hat{g}_k \hat{g}_i \mathbf{w}_I | \hat{g}_j \mathbf{w}_J},
100    /// ```
101    ///
102    /// arranged in a rank-3 array indexed by $`k`$, $`iI`$, and $`jJ`$.
103    ///
104    /// # Arguments
105    ///
106    /// * `metric` - The atomic-orbital overlap matrix with respect to the conventional sesquilinear
107    ///   inner product.
108    /// * `metric_h` - The atomic-orbital overlap matrix with respect to the bilinear inner product.
109    ///
110    /// # Returns
111    ///
112    /// The basis symmetry-orbit metric as a rank-3 array indexed as described above.
113    fn basis_symmetry_orbit_metric(
114        &self,
115        _group: &G,
116        _symmetry_transformation_kind: &SymmetryTransformationKind,
117        metric: Option<&Array2<T>>,
118        metric_h: Option<&Array2<T>>,
119    ) -> Result<Array3<T>, anyhow::Error> {
120        log::debug!("Computing basis symmetry-orbit metric for `OrbitBasis`...");
121        let order = self.group.order();
122        let origins = self.origins();
123        let n_origins = origins.len();
124
125        if let Some(ctb) = self.group.cayley_table() {
126            log::debug!(
127                "Cayley table available. Group closure will be used to speed up `OrbitBasis` metric computation."
128            );
129            let detov_kpwx_vec = self
130                .iter()
131                .collect::<Result<Vec<_>, _>>()?
132                .into_iter()
133                .cartesian_product(origins.iter())
134                .map(|(w, x)| w.overlap(x, metric, metric_h))
135                .collect::<Result<Vec<_>, _>>()?;
136            let detov_kpwx = Array3::from_shape_vec((order, n_origins, n_origins), detov_kpwx_vec)?;
137
138            let norm_preserving_scalar_map = |i, v: T| {
139                if origins.iter().any(|origin| origin.complex_symmetric()) {
140                    Err(format_err!(
141                        "`norm_preserving_scalar_map` is currently not implemented for complex-symmetric overlaps. This thus precludes the use of the Cayley table to speed up the computation of the orbit overlap matrix."
142                    ))
143                } else if self
144                    .group
145                    .get_index(i)
146                    .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
147                    .contains_time_reversal()
148                {
149                    Ok(v.conj())
150                } else {
151                    Ok(v)
152                }
153            };
154            let ukipwjpx_vec = [
155                (0..order),
156                (0..order),
157                (0..n_origins),
158                (0..order),
159                (0..n_origins),
160            ]
161            .into_iter()
162            .multi_cartesian_product()
163            .map(|v| {
164                let k = v[0];
165                let ip = v[1];
166                let w = v[2];
167                let jp = v[3];
168                let x = v[4];
169
170                let jpinv =
171                    ctb.slice(s![.., jp])
172                        .iter()
173                        .position(|&x| x == 0)
174                        .ok_or(format_err!(
175                            "Unable to find the inverse of group element `{jp}`."
176                        ))?;
177                let kp = ctb[(jpinv, ctb[(k, ip)])];
178                norm_preserving_scalar_map(jp, detov_kpwx[(kp, w, x)])
179            })
180            .collect::<Result<Vec<_>, _>>()?;
181            Ok(Array3::from_shape_vec(
182                (order, order * n_origins, order * n_origins),
183                ukipwjpx_vec,
184            )?)
185        } else {
186            let ukipwjpx_vec = iproduct!(
187                self.group.elements().clone().into_iter(),
188                self.iter().collect::<Result<Vec<_>, _>>()?.into_iter(),
189                self.iter().collect::<Result<Vec<_>, _>>()?.into_iter()
190            )
191            .map(|(gk, gipww, gjpwx)| {
192                let gkgipww = (self.action)(&gk, &gipww)?;
193                gkgipww.overlap(&gjpwx, metric, metric_h)
194            })
195            .collect::<Result<Vec<_>, _>>()?;
196            Ok(Array3::from_shape_vec(
197                (order, order * n_origins, order * n_origins),
198                ukipwjpx_vec,
199            )?)
200        }
201    }
202}
203
204// ~~~~~~~~~~
205// EagerBasis
206// ~~~~~~~~~~
207
208impl<G, I, T> BasisSymmetryOrbitMetric<G, T> for EagerBasis<I>
209where
210    G: SymmetryGroupProperties + Clone,
211    I: Clone + Overlap<T, Ix2> + Sync + Send,
212    T: ComplexFloat + fmt::Debug + Lapack + Sync + Send,
213    <T as ComplexFloat>::Real: Sync + Send,
214    EagerBasis<I>: SymmetryTransformable,
215{
216    fn basis_symmetry_orbit_metric(
217        &self,
218        group: &G,
219        symmetry_tranformation_kind: &SymmetryTransformationKind,
220        metric: Option<&Array2<T>>,
221        metric_h: Option<&Array2<T>>,
222    ) -> Result<Array3<T>, anyhow::Error> {
223        log::debug!("Computing basis symmetry-orbit metric for `EagerBasis`...");
224        let g_self_orbit = BasisSymmetryOrbit::builder()
225            .group(group)
226            .origin(self)
227            .symmetry_transformation_kind(symmetry_tranformation_kind.clone())
228            .build()?;
229        let ovss_vec = g_self_orbit
230            .iter()
231            .map(|g_basis_res| {
232                g_basis_res.and_then(|g_basis| {
233                    let mut ovs_vec = g_basis
234                        .elements
235                        .iter()
236                        .cartesian_product(self.elements.iter())
237                        .enumerate()
238                        .par_bridge()
239                        .map(|(i, (g_item, item))| {
240                            let ov = g_item.overlap(item, metric, metric_h)?;
241                            Ok((i, ov))
242                        })
243                        .collect::<Result<Vec<_>, anyhow::Error>>()?;
244                    ovs_vec.sort_by_key(|v| v.0);
245                    Ok(ovs_vec.into_iter().map(|v| v.1).collect_vec())
246                })
247            })
248            .collect::<Result<Vec<_>, _>>()?
249            .into_iter()
250            .flatten()
251            .collect_vec();
252        Ok(Array3::from_shape_vec(
253            (group.order(), self.n_items(), self.n_items()),
254            ovss_vec,
255        )?)
256    }
257}
258
259// ~~~~~~~~
260// FCIBasis
261// ~~~~~~~~
262
263impl<'a, G, T, SC> BasisSymmetryOrbitMetric<G, T> for FCIBasis<'a, T, SC>
264where
265    G: SymmetryGroupProperties + Clone,
266    T: ComplexFloat + Lapack + Sync + Send,
267    <T as ComplexFloat>::Real: Sync + Send,
268    SC: 'a + StructureConstraint + fmt::Display + Clone + PartialEq + Sync + Send,
269    FCIBasis<'a, T, SC>: SymmetryTransformable,
270{
271    /// Computes the basis symmetry-orbit metric for this configuration-interaction basis.
272    ///
273    /// For a group $`\mathcal{G}`$ and a Slater-determinant reference $`\Psi_0`$, this metric is
274    /// defined as
275    ///
276    /// ```math
277    /// \braket{\hat{g} \Psi_I | \Psi_J}
278    /// ```
279    ///
280    /// where $`g \in \mathcal{G}`$, and $`\Psi_I`$ and $`\Psi_J`$ are replacement Slater
281    /// determinants constructed from $`\Psi_0`$.
282    /// These determinants are specified by the occupation patterns with respect to the set of
283    /// molecular orbitals defined by $`\Psi_0`$.
284    ///
285    /// # Arguments
286    ///
287    /// * `group` - The group $`\mathcal{G}`$.
288    /// * `symmetry_transformation_kind` - Enum specifying how the elements of $`\mathcal{G}`$ act
289    ///   on the elements of the basis.
290    /// * `metric` - The atomic-orbital overlap matrix with respect to the conventional sesquilinear
291    ///   inner product.
292    /// * `metric_h` - The atomic-orbital overlap matrix with respect to the bilinear inner product.
293    ///
294    /// # Returns
295    ///
296    /// The basis symmetry-orbit metric as a rank-3 array indexed by $`g`$, $`I`$, and $`J`$.
297    fn basis_symmetry_orbit_metric(
298        &self,
299        group: &G,
300        symmetry_tranformation_kind: &SymmetryTransformationKind,
301        metric: Option<&Array2<T>>,
302        metric_h: Option<&Array2<T>>,
303    ) -> Result<Array3<T>, anyhow::Error> {
304        log::debug!("Computing basis symmetry-orbit metric for `FCIBasis`...");
305        let g_self_orbit = BasisSymmetryOrbit::builder()
306            .group(group)
307            .origin(self)
308            .symmetry_transformation_kind(symmetry_tranformation_kind.clone())
309            .build()?;
310        let sao = metric.ok_or_else(|| format_err!("No atomic-orbital metric found."))?;
311        let sao_h = metric_h.unwrap_or(sao);
312
313        let s_ref = self.reference();
314        let thresh = s_ref.threshold();
315
316        let ovss_vec = g_self_orbit
317            .iter()
318            .map(|g_self_res| {
319                g_self_res.and_then(|g_self| {
320                    let g_s_ref = g_self.reference();
321                    let mo_ov_mats = g_s_ref
322                        .coefficients()
323                        .iter()
324                        .zip(s_ref.coefficients().iter())
325                        .map(|(cw, cx)| {
326                            if s_ref.complex_symmetric() {
327                                match (g_s_ref.complex_conjugated(), s_ref.complex_conjugated()) {
328                                    (false, false) => cw.t().dot(sao_h).dot(cx),
329                                    (true, false) => cw.t().dot(sao).dot(cx),
330                                    (false, true) => cx.t().dot(sao).dot(cw),
331                                    (true, true) => cw.t().dot(&sao_h.t()).dot(cx),
332                                }
333                            } else {
334                                match (g_s_ref.complex_conjugated(), s_ref.complex_conjugated()) {
335                                    (false, false) => cw.t().mapv(|x| x.conj()).dot(sao).dot(cx),
336                                    (true, false) => cw.t().mapv(|x| x.conj()).dot(sao_h).dot(cx),
337                                    (false, true) => cx
338                                        .t()
339                                        .mapv(|x| x.conj())
340                                        .dot(sao_h)
341                                        .dot(cw)
342                                        .mapv(|x| x.conj()),
343                                    (true, true) => cw.t().mapv(|x| x.conj()).dot(&sao.t()).dot(cx),
344                                }
345                            }
346                        })
347                        .collect_vec();
348
349                    let mut ovs_vec = g_self
350                        .occupation_patterns
351                        .iter()
352                        .cartesian_product(self.occupation_patterns.iter())
353                        .enumerate()
354                        .par_bridge()
355                        .map(|(i, (g_s_occs, s_occs))| {
356                            let ov = izip!(g_s_occs, s_occs, &mo_ov_mats)
357                                .map(|(g_s_occ, s_occ, mo_ov_mat)| {
358                                    let nonzero_g_s_occ =
359                                        g_s_occ.iter().positions(|&occ| occ > thresh).collect_vec();
360                                    let nonzero_s_occ =
361                                        s_occ.iter().positions(|&occ| occ > thresh).collect_vec();
362                                    let mo_ov_mat_occ = mo_ov_mat
363                                        .select(Axis(0), &nonzero_g_s_occ)
364                                        .select(Axis(1), &nonzero_s_occ);
365                                    mo_ov_mat_occ.det().expect(
366                                    "The determinant of the MO overlap matrix could not be found.",
367                                )
368                                })
369                                .fold(T::one(), |acc, x| acc * x);
370
371                            let implicit_factor = s_ref.structure_constraint().implicit_factor()?;
372                            if implicit_factor > 1 {
373                                let p_i32 = i32::try_from(implicit_factor)?;
374                                Ok((i, ComplexFloat::powi(ov, p_i32)))
375                            } else {
376                                Ok((i, ov))
377                            }
378                        })
379                        .collect::<Result<Vec<_>, anyhow::Error>>()?;
380                    ovs_vec.sort_by_key(|v| v.0);
381                    let ovs_vec = ovs_vec.into_iter().map(|v| v.1).collect_vec();
382                    Ok(ovs_vec)
383                })
384            })
385            .collect::<Result<Vec<_>, _>>()?
386            .into_iter()
387            .flatten()
388            .collect_vec();
389
390        Ok(Array3::from_shape_vec(
391            (group.order(), self.n_items(), self.n_items()),
392            ovss_vec,
393        )?)
394    }
395}
qsym2/target/noci/basis_metric.rs

qsym2/target/noci/
basis_metric.rs
