qsym2/target/vibration/
mod.rs

1//! Vibrational coordinates for normal modes.
2
3use std::fmt;
4
5use approx;
6use derive_builder::Builder;
7use ndarray::{Array1, Array2};
8use ndarray_linalg::types::Lapack;
9use num_complex::{Complex, ComplexFloat};
10use num_traits::float::{Float, FloatConst};
11
12use crate::auxiliary::molecule::Molecule;
13
14#[cfg(test)]
15mod vibration_tests;
16
17pub mod vibration_analysis;
18mod vibration_transformation;
19
20// ==================
21// Struct definitions
22// ==================
23
24// ---------------------
25// VibrationalCoordinate
26// ---------------------
27
28/// Structure to manage vibrational coordinates.
29#[derive(Builder, Clone)]
30#[builder(build_fn(validate = "Self::validate"))]
31pub struct VibrationalCoordinate<'a, T>
32where
33    T: ComplexFloat + Lapack,
34{
35    /// The associated molecule.
36    mol: &'a Molecule,
37
38    /// The frequency of this vibration in $`\mathrm{cm}^{-1}`$.
39    frequency: T,
40
41    /// The coefficients describing this vibrational coordinate in terms of the local Cartesian
42    /// coordinates.
43    coefficients: Array1<T>,
44
45    /// The threshold for comparing vibrational coordinates.
46    threshold: <T as ComplexFloat>::Real,
47}
48
49impl<'a, T> VibrationalCoordinateBuilder<'a, T>
50where
51    T: ComplexFloat + Lapack,
52{
53    fn validate(&self) -> Result<(), String> {
54        let coefficients = self
55            .coefficients
56            .as_ref()
57            .ok_or("No coefficients found.".to_string())?;
58        let mol = self.mol.ok_or("No molecule found.".to_string())?;
59        let natoms = 3 * mol.atoms.len() == coefficients.len();
60        if !natoms {
61            log::error!("The number of coefficients for this vibrational coordinate does not match the number of atoms in the molecule.");
62        }
63
64        if natoms {
65            Ok(())
66        } else {
67            Err("Vibrational coordinate validation failed.".to_string())
68        }
69    }
70}
71
72impl<'a, T> VibrationalCoordinate<'a, T>
73where
74    T: ComplexFloat + Clone + Lapack,
75{
76    /// Returns a builder to construct a new [`VibrationalCoordinate`].
77    pub fn builder() -> VibrationalCoordinateBuilder<'a, T> {
78        VibrationalCoordinateBuilder::default()
79    }
80
81    /// Returns the frequency of the vibration.
82    pub fn frequency(&self) -> &T {
83        &self.frequency
84    }
85
86    /// Returns a shared reference to the coefficient array.
87    pub fn coefficients(&self) -> &Array1<T> {
88        &self.coefficients
89    }
90}
91
92// =====================
93// Trait implementations
94// =====================
95
96// ----
97// From
98// ----
99impl<'a, T> From<VibrationalCoordinate<'a, T>> for VibrationalCoordinate<'a, Complex<T>>
100where
101    T: Float + FloatConst + Lapack,
102    Complex<T>: Lapack,
103{
104    fn from(value: VibrationalCoordinate<'a, T>) -> Self {
105        VibrationalCoordinate::<'a, Complex<T>>::builder()
106            .coefficients(value.coefficients.map(Complex::from))
107            .mol(value.mol)
108            .frequency(Complex::from(value.frequency))
109            .threshold(value.threshold)
110            .build()
111            .expect("Unable to construct a complex vibrational coordinate.")
112    }
113}
114
115// ---------
116// PartialEq
117// ---------
118impl<'a, T> PartialEq for VibrationalCoordinate<'a, T>
119where
120    T: ComplexFloat<Real = f64> + Lapack,
121{
122    fn eq(&self, other: &Self) -> bool {
123        let thresh = (self.threshold * other.threshold).sqrt();
124        let coefficients_eq = approx::relative_eq!(
125            (&self.coefficients - &other.coefficients)
126                .map(|x| ComplexFloat::abs(*x).powi(2))
127                .sum()
128                .sqrt(),
129            0.0,
130            epsilon = thresh,
131            max_relative = thresh,
132        );
133        let frequency_eq = approx::relative_eq!(
134            ComplexFloat::abs(self.frequency - other.frequency),
135            0.0,
136            epsilon = thresh,
137            max_relative = thresh,
138        );
139        self.mol == other.mol && coefficients_eq && frequency_eq
140    }
141}
142
143// --
144// Eq
145// --
146impl<'a, T> Eq for VibrationalCoordinate<'a, T> where T: ComplexFloat<Real = f64> + Lapack {}
147
148// -----
149// Debug
150// -----
151impl<'a, T> fmt::Debug for VibrationalCoordinate<'a, T>
152where
153    T: fmt::Debug + ComplexFloat + Lapack,
154{
155    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
156        write!(
157            f,
158            "VibrationalCoordinate[frequency = {:?}, coefficient array of length {}]",
159            self.frequency,
160            self.coefficients.len()
161        )?;
162        Ok(())
163    }
164}
165
166// -------
167// Display
168// -------
169impl<'a, T> fmt::Display for VibrationalCoordinate<'a, T>
170where
171    T: fmt::Display + ComplexFloat + Lapack,
172{
173    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
174        write!(
175            f,
176            "VibrationalCoordinate[frequency = {}, coefficient array of length {}]",
177            self.frequency,
178            self.coefficients.len()
179        )?;
180        Ok(())
181    }
182}
183
184// -------------------------------
185// VibrationalCoordinateCollection
186// -------------------------------
187
188/// Structure to manage multiple vibrational coordinates of a single molecule.
189#[derive(Builder, Clone)]
190#[builder(build_fn(validate = "Self::validate"))]
191pub struct VibrationalCoordinateCollection<'a, T>
192where
193    T: ComplexFloat + Lapack,
194{
195    /// The associated molecule.
196    mol: &'a Molecule,
197
198    /// The frequencies of the vibrations in $`\mathrm{cm}^{-1}`$.
199    frequencies: Array1<T>,
200
201    /// The coefficients describing this vibrational coordinates in terms of the local Cartesian
202    /// coordinates. Each column corresponds to one vibrational coordinate.
203    coefficients: Array2<T>,
204
205    /// The threshold for comparing vibrational coordinates.
206    threshold: <T as ComplexFloat>::Real,
207}
208
209impl<'a, T> VibrationalCoordinateCollectionBuilder<'a, T>
210where
211    T: ComplexFloat + Lapack,
212{
213    fn validate(&self) -> Result<(), String> {
214        let coefficients = self
215            .coefficients
216            .as_ref()
217            .ok_or("No coefficients found.".to_string())?;
218        let mol = self.mol.ok_or("No molecule found.".to_string())?;
219        let natoms = 3 * mol.atoms.len() == coefficients.nrows();
220        if !natoms {
221            log::error!("The number of coefficient components does not match the number of atoms in the molecule.");
222        }
223
224        let frequencies = self
225            .frequencies
226            .as_ref()
227            .ok_or("No frequencies found.".to_string())?;
228        let nfreqs = frequencies.len() == coefficients.ncols();
229
230        if natoms && nfreqs {
231            Ok(())
232        } else {
233            Err("Vibrational coordinate collection validation failed.".to_string())
234        }
235    }
236}
237
238impl<'a, T> VibrationalCoordinateCollection<'a, T>
239where
240    T: ComplexFloat + Clone + Lapack,
241{
242    /// Returns a builder to construct a new [`VibrationalCoordinateCollection`].
243    pub fn builder() -> VibrationalCoordinateCollectionBuilder<'a, T> {
244        VibrationalCoordinateCollectionBuilder::default()
245    }
246
247    /// Returns the number of vibrational modes in this collection.
248    pub fn n_modes(&self) -> usize {
249        self.frequencies.len()
250    }
251
252    /// Returns the frequencies of the vibrations.
253    pub fn frequencies(&self) -> &Array1<T> {
254        &self.frequencies
255    }
256
257    /// Returns a shared reference to the coefficient array.
258    pub fn coefficients(&self) -> &Array2<T> {
259        &self.coefficients
260    }
261
262    /// Returns a vector of separate vibrational coordinates from this collection.
263    pub fn to_vibrational_coordinates(&self) -> Vec<VibrationalCoordinate<'a, T>> {
264        self.frequencies
265            .iter()
266            .zip(self.coefficients.columns())
267            .map(|(freq, vib_coord)| {
268                VibrationalCoordinate::builder()
269                    .coefficients(vib_coord.to_owned())
270                    .frequency(*freq)
271                    .mol(self.mol)
272                    .threshold(self.threshold)
273                    .build()
274                    .expect("Unable to construct a vibrational coordinate.")
275            })
276            .collect::<Vec<_>>()
277    }
278}
279
280// =====================
281// Trait implementations
282// =====================
283
284// ----
285// From
286// ----
287impl<'a, T> From<VibrationalCoordinateCollection<'a, T>>
288    for VibrationalCoordinateCollection<'a, Complex<T>>
289where
290    T: Float + FloatConst + Lapack,
291    Complex<T>: Lapack,
292{
293    fn from(value: VibrationalCoordinateCollection<'a, T>) -> Self {
294        VibrationalCoordinateCollection::<'a, Complex<T>>::builder()
295            .coefficients(value.coefficients.map(Complex::from))
296            .mol(value.mol)
297            .frequencies(value.frequencies.map(Complex::from))
298            .threshold(value.threshold)
299            .build()
300            .expect("Unable to construct a complex vibrational coordinate collection.")
301    }
302}
303
304// -----
305// Debug
306// -----
307impl<'a, T> fmt::Debug for VibrationalCoordinateCollection<'a, T>
308where
309    T: fmt::Debug + ComplexFloat + Lapack,
310{
311    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
312        let n = self.frequencies.len();
313        write!(
314            f,
315            "VibrationalCoordinateCollection[{} {}, coefficient array of shape {} × {}]",
316            n,
317            if n == 1 { "mode" } else { "modes" },
318            self.coefficients.nrows(),
319            self.coefficients.ncols(),
320        )?;
321        Ok(())
322    }
323}
324
325// -------
326// Display
327// -------
328impl<'a, T> fmt::Display for VibrationalCoordinateCollection<'a, T>
329where
330    T: fmt::Display + ComplexFloat + Lapack,
331{
332    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
333        let n = self.frequencies.len();
334        write!(
335            f,
336            "VibrationalCoordinateCollection[{} {}, coefficient array of shape {} × {}]",
337            n,
338            if n == 1 { "mode" } else { "modes" },
339            self.coefficients.nrows(),
340            self.coefficients.ncols(),
341        )?;
342        Ok(())
343    }
344}