qsym2/target/noci/
basis.rs

1//! Basis for configuration interaction of Slater determinants.
2
3use std::collections::VecDeque;
4use std::fmt;
5
6use anyhow::{self, format_err};
7use derive_builder::Builder;
8use itertools::Itertools;
9use itertools::structs::Product;
10use ndarray::Array1;
11use ndarray_linalg::types::Lapack;
12use num_complex::ComplexFloat;
13
14use crate::angmom::spinor_rotation_3d::StructureConstraint;
15use crate::group::GroupProperties;
16use crate::target::determinant::SlaterDeterminant;
17
18#[path = "basis_transformation.rs"]
19mod basis_transformation;
20
21#[path = "basis_metric.rs"]
22pub mod basis_metric;
23
24#[path = "basis_orbit.rs"]
25pub mod basis_orbit;
26
27#[cfg(test)]
28#[path = "basis_tests.rs"]
29mod basis_tests;
30
31// =====
32// Basis
33// =====
34
35// -----------------
36// Trait definitions
37// -----------------
38
39/// Trait defining behaviours of a basis consisting of linear-space items.
40pub trait Basis {
41    /// Type of the item in the basis.
42    type Item;
43
44    /// Type of the iterator over items in the basis.
45    type BasisIter<'b>: Iterator<Item = Result<Self::Item, anyhow::Error>>
46    where
47        Self: 'b;
48
49    /// Returns the number of items in the basis.
50    fn n_items(&self) -> usize;
51
52    /// An iterator over items in the basis.
53    fn iter(&self) -> Self::BasisIter<'_>;
54
55    /// Shared reference to the first item in the basis.
56    fn first(&self) -> Option<Self::Item>;
57}
58
59// --------------------------------------
60// Struct definitions and implementations
61// --------------------------------------
62
63// ~~~~~~~~~~~~~~~~~~~~~~
64// Lazy basis from orbits
65// ~~~~~~~~~~~~~~~~~~~~~~
66
67/// Basis defined by a set of origins $`\{\mathbf{w}_I\}`$ and a group $`\mathcal{G} = \{g_i\}`$.
68///
69/// Each element in the basis is given by the action of the group on the origins, $`\mathcal{G}
70/// \cdot \{ \mathbf{w}_I \} = \{ \hat{g}_i \mathbf{w}_I \}`$.
71#[derive(Builder, Clone)]
72pub struct OrbitBasis<'g, G, I>
73where
74    G: GroupProperties,
75{
76    /// The origins from which orbits are generated.
77    origins: Vec<I>,
78
79    /// The group acting on the origins to generate orbits, the concatenation of which forms the
80    /// basis.
81    group: &'g G,
82
83    /// Additional operators acting on the entire orbit basis (right-most operator acts first). Each
84    /// operator has an associated action that defines how it operates on the elements in the
85    /// orbit basis.
86    #[builder(default = "None")]
87    #[allow(clippy::type_complexity)]
88    prefactors: Option<
89        VecDeque<(
90            G::GroupElement,
91            fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
92        )>,
93    >,
94
95    /// A function defining the action of each group element on the origin.
96    action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
97}
98
99impl<'g, G, I> OrbitBasis<'g, G, I>
100where
101    G: GroupProperties + Clone,
102    I: Clone,
103{
104    pub fn builder() -> OrbitBasisBuilder<'g, G, I> {
105        OrbitBasisBuilder::<'g, G, I>::default()
106    }
107
108    /// The origins from which orbits are generated.
109    pub fn origins(&self) -> &Vec<I> {
110        &self.origins
111    }
112
113    /// The group acting on the origins to generate orbits, the concatenation of which forms the
114    /// basis.
115    pub fn group(&self) -> &G {
116        self.group
117    }
118
119    /// Additional operators acting on the entire orbit basis (right-most operator acts first). Each
120    /// operator has an associated action that defines how it operatres on the elements in the
121    /// orbit basis.
122    #[allow(clippy::type_complexity)]
123    pub fn prefactors(
124        &self,
125    ) -> Option<
126        &VecDeque<(
127            G::GroupElement,
128            fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
129        )>,
130    > {
131        self.prefactors.as_ref()
132    }
133}
134
135impl<'g, G, I> OrbitBasis<'g, G, I>
136where
137    G: GroupProperties,
138    I: Clone,
139{
140    /// Converts this orbit basis into the equivalent eager basis.
141    pub fn to_eager(&self) -> Result<EagerBasis<I>, anyhow::Error> {
142        let elements = self.iter().collect::<Result<Vec<_>, _>>()?;
143        EagerBasis::builder()
144            .elements(elements)
145            .build()
146            .map_err(|err| format_err!(err))
147    }
148}
149
150impl<'g, G, I> Basis for OrbitBasis<'g, G, I>
151where
152    G: GroupProperties,
153    I: Clone,
154{
155    type Item = I;
156
157    type BasisIter<'b>
158        = OrbitBasisIterator<G, I>
159    where
160        Self: 'b;
161
162    fn n_items(&self) -> usize {
163        self.origins.len() * self.group.order()
164    }
165
166    /// Iterates over the elements $`\hat{g}_i \mathbf{w}_I`$ of the [`OrbitBasis`]. Each element is
167    /// indexed by `iI` where `i` enumerates the group elements and `I` enumerates the origins. `I`
168    /// is the fast index and `i` the slow one.
169    fn iter(&self) -> Self::BasisIter<'_> {
170        OrbitBasisIterator::new(
171            self.prefactors.clone(),
172            self.group,
173            self.origins.clone(),
174            self.action,
175        )
176    }
177
178    fn first(&self) -> Option<Self::Item> {
179        if let Some(prefactors) = self.prefactors.as_ref() {
180            prefactors
181                .iter()
182                .rev()
183                .try_fold(self.origins.first()?.clone(), |acc, (symop, action)| {
184                    (action)(symop, &acc).ok()
185                })
186        } else {
187            self.origins.first().cloned()
188        }
189    }
190}
191
192// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
193// Iterator for lazy basis from orbits
194// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
195
196/// Lazy iterator for basis constructed from the concatenation of orbits generated from multiple
197/// origins.
198pub struct OrbitBasisIterator<G, I>
199where
200    G: GroupProperties,
201{
202    #[allow(clippy::type_complexity)]
203    prefactors: Option<
204        VecDeque<(
205            G::GroupElement,
206            fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
207        )>,
208    >,
209
210    /// A mutable iterator over the Cartesian product between the group elements and the origins.
211    group_origin_iter: Product<
212        <<G as GroupProperties>::ElementCollection as IntoIterator>::IntoIter,
213        std::vec::IntoIter<I>,
214    >,
215
216    /// A function defining the action of each group element on the origin.
217    action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
218}
219
220impl<G, I> OrbitBasisIterator<G, I>
221where
222    G: GroupProperties,
223    I: Clone,
224{
225    /// Creates and returns a new orbit basis iterator.
226    ///
227    /// Each element returned by the iterator is indexed by `iI` where `i` enumerates the group
228    /// elements and `I` enumerates the origins. `I` is the fast index and `i` the slow one.
229    ///
230    /// # Arguments
231    ///
232    /// * `group` - A group.
233    /// * `origins` - A slice of origins.
234    /// * `action` - A function or closure defining the action of each group element on the origins.
235    ///
236    /// # Returns
237    ///
238    /// An orbit basis iterator.
239    #[allow(clippy::type_complexity)]
240    fn new(
241        prefactors: Option<
242            VecDeque<(
243                G::GroupElement,
244                fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
245            )>,
246        >,
247        group: &G,
248        origins: Vec<I>,
249        action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
250    ) -> Self {
251        Self {
252            prefactors,
253            group_origin_iter: group
254                .elements()
255                .clone()
256                .into_iter()
257                .cartesian_product(origins),
258            action,
259        }
260    }
261}
262
263impl<G, I> Iterator for OrbitBasisIterator<G, I>
264where
265    G: GroupProperties,
266    I: Clone,
267{
268    type Item = Result<I, anyhow::Error>;
269
270    fn next(&mut self) -> Option<Self::Item> {
271        if let Some(prefactors) = self.prefactors.as_ref() {
272            self.group_origin_iter.next().and_then(|(op, origin)| {
273                prefactors
274                    .iter()
275                    .rev()
276                    .try_fold((self.action)(&op, &origin), |acc_res, (symop, action)| {
277                        acc_res.map(|acc| (action)(symop, &acc))
278                    })
279                    .ok()
280            })
281        } else {
282            self.group_origin_iter
283                .next()
284                .map(|(op, origin)| (self.action)(&op, &origin))
285        }
286    }
287}
288
289// ~~~~~~~~~~~
290// Eager basis
291// ~~~~~~~~~~~
292
293/// Basis defined by specifying all of its elements explicitly.
294#[derive(Builder, Clone)]
295pub struct EagerBasis<I: Clone> {
296    /// The elements in the basis.
297    elements: Vec<I>,
298}
299
300impl<I: Clone> EagerBasis<I> {
301    pub fn builder() -> EagerBasisBuilder<I> {
302        EagerBasisBuilder::default()
303    }
304}
305
306impl<I: Clone> Basis for EagerBasis<I> {
307    type Item = I;
308
309    type BasisIter<'b>
310        = std::vec::IntoIter<Result<I, anyhow::Error>>
311    where
312        Self: 'b;
313
314    fn n_items(&self) -> usize {
315        self.elements.len()
316    }
317
318    fn iter(&self) -> Self::BasisIter<'_> {
319        self.elements
320            .iter()
321            .cloned()
322            .map(Ok)
323            .collect_vec()
324            .into_iter()
325    }
326
327    fn first(&self) -> Option<Self::Item> {
328        self.elements.first().cloned()
329    }
330}
331
332// ~~~~~~~~~
333// FCI basis
334// ~~~~~~~~~
335
336/// Basis defined by a reference Slater determinant $`\Psi_0`$ and its replacement determinants
337/// specified by occupation patterns.
338#[derive(Builder, Clone)]
339#[builder(build_fn(validate = "Self::validate"))]
340pub struct FCIBasis<'a, T, SC>
341where
342    T: ComplexFloat + Lapack,
343    SC: StructureConstraint + fmt::Display,
344{
345    /// The reference determinant of the full-configuration-interaction basis. This determinant
346    /// carries both occupied and virtual molecular orbitals.
347    reference: SlaterDeterminant<'a, T, SC>,
348
349    /// Vector of full-configuration-interaction occupation patterns.
350    occupation_patterns: Vec<Vec<Array1<<T as ComplexFloat>::Real>>>,
351}
352
353impl<'a, T, SC> FCIBasisBuilder<'a, T, SC>
354where
355    T: ComplexFloat + Lapack,
356    SC: StructureConstraint + fmt::Display,
357{
358    fn validate(&self) -> Result<(), String> {
359        if let Some(ref reference) = self.reference
360            && let Some(ref occupation_patterns) = self.occupation_patterns
361        {
362            let reference_occs = reference.occupations();
363            if occupation_patterns.iter().all(|occs| {
364                occs.len() == reference_occs.len()
365                    && occs
366                        .iter()
367                        .zip(reference_occs)
368                        .all(|(occ, reference_occ)| occ.shape() == reference_occ.shape())
369            }) {
370                Ok(())
371            } else {
372                Err("Inconsistent occupation patterns.".to_string())
373            }
374        } else {
375            Err("Missing `reference` or `occupation_patterns`.".to_string())
376        }
377    }
378}
379
380impl<'a, T, SC> FCIBasis<'a, T, SC>
381where
382    T: ComplexFloat + Lapack,
383    SC: StructureConstraint + fmt::Display + Clone,
384{
385    pub fn builder() -> FCIBasisBuilder<'a, T, SC> {
386        FCIBasisBuilder::<T, SC>::default()
387    }
388
389    /// The reference determinant of the full-configuration-interaction basis.
390    pub fn reference(&self) -> &SlaterDeterminant<'a, T, SC> {
391        &self.reference
392    }
393
394    /// The full-configuration-interaction occupation patterns.
395    pub fn occupation_patterns(&self) -> &Vec<Vec<Array1<<T as ComplexFloat>::Real>>> {
396        &self.occupation_patterns
397    }
398}
399
400impl<'a, T, SC> Basis for FCIBasis<'a, T, SC>
401where
402    T: ComplexFloat + Lapack,
403    SC: 'a + StructureConstraint + fmt::Display + Clone,
404{
405    type Item = SlaterDeterminant<'a, T, SC>;
406
407    type BasisIter<'b>
408        = FCIBasisIterator<'b, 'a, T, SC>
409    where
410        Self: 'b;
411
412    fn n_items(&self) -> usize {
413        self.occupation_patterns.len()
414    }
415
416    /// Iterates over the elements of the [`FCIBasis`]. Each element is a Slater determinant
417    /// obtained by replacing the occupation pattern in the reference with the corresponding one
418    /// from the [`Self::occupation_patterns`] list.
419    fn iter(&self) -> Self::BasisIter<'_> {
420        FCIBasisIterator::new(self)
421    }
422
423    fn first(&self) -> Option<SlaterDeterminant<'a, T, SC>> {
424        Some(self.reference.clone())
425    }
426}
427
428// ~~~~~~~~~~~~~~~~~~
429// FCIBasis interator
430// ~~~~~~~~~~~~~~~~~~
431
432/// Lazy iterator for full-configuration-interaction basis.
433#[derive(Builder, Clone)]
434pub struct FCIBasisIterator<'b, 'a, T, SC>
435where
436    T: ComplexFloat + Lapack,
437    SC: StructureConstraint + fmt::Display,
438{
439    /// The associated full-configuration-interaction basis.
440    fci_basis: &'b FCIBasis<'a, T, SC>,
441
442    /// The current index of the occupation patterns.
443    current_occupation_index: usize,
444}
445
446impl<'b, 'a, T, SC> FCIBasisIterator<'b, 'a, T, SC>
447where
448    T: ComplexFloat + Lapack,
449    SC: StructureConstraint + fmt::Display,
450{
451    /// Creates and returns a new full-configuration-interaction basis iterator.
452    ///
453    /// # Arguments
454    ///
455    /// * `fci_basis` - The associated full-configuration-interaction basis.
456    ///
457    /// # Returns
458    ///
459    /// An full-configuration-interaction basis iterator.
460    fn new(fci_basis: &'b FCIBasis<'a, T, SC>) -> Self {
461        Self {
462            fci_basis,
463            current_occupation_index: 0,
464        }
465    }
466}
467
468impl<'b, 'a, T, SC> Iterator for FCIBasisIterator<'b, 'a, T, SC>
469where
470    T: ComplexFloat + Lapack,
471    SC: StructureConstraint + fmt::Display + Clone,
472{
473    type Item = Result<SlaterDeterminant<'a, T, SC>, anyhow::Error>;
474
475    fn next(&mut self) -> Option<Self::Item> {
476        let reference = &self.fci_basis.reference;
477        let index = self.current_occupation_index;
478        self.current_occupation_index += 1;
479        if index >= self.fci_basis.occupation_patterns.len() {
480            None
481        } else {
482            let det = SlaterDeterminant::builder()
483                .structure_constraint(reference.structure_constraint().clone())
484                .baos(reference.baos().clone())
485                .complex_symmetric(reference.complex_symmetric())
486                .complex_conjugated(reference.complex_conjugated())
487                .mol(reference.mol())
488                .coefficients(reference.coefficients())
489                .occupations(self.fci_basis.occupation_patterns.get(index)?)
490                .mo_energies(reference.mo_energies().cloned())
491                .threshold(reference.threshold())
492                .build()
493                .map_err(|err| format_err!(err));
494            Some(det)
495        }
496    }
497}
qsym2/target/noci/basis.rs

qsym2/target/noci/
basis.rs
