1use 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#[derive(Builder, Clone)]
30#[builder(build_fn(validate = "Self::validate"))]
31pub struct VibrationalCoordinate<'a, T>
32where
33 T: ComplexFloat + Lapack,
34{
35 mol: &'a Molecule,
37
38 frequency: T,
40
41 coefficients: Array1<T>,
44
45 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 pub fn builder() -> VibrationalCoordinateBuilder<'a, T> {
78 VibrationalCoordinateBuilder::default()
79 }
80
81 pub fn frequency(&self) -> &T {
83 &self.frequency
84 }
85
86 pub fn coefficients(&self) -> &Array1<T> {
88 &self.coefficients
89 }
90}
91
92impl<'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
115impl<'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
143impl<'a, T> Eq for VibrationalCoordinate<'a, T> where T: ComplexFloat<Real = f64> + Lapack {}
147
148impl<'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
166impl<'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#[derive(Builder, Clone)]
190#[builder(build_fn(validate = "Self::validate"))]
191pub struct VibrationalCoordinateCollection<'a, T>
192where
193 T: ComplexFloat + Lapack,
194{
195 mol: &'a Molecule,
197
198 frequencies: Array1<T>,
200
201 coefficients: Array2<T>,
204
205 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 pub fn builder() -> VibrationalCoordinateCollectionBuilder<'a, T> {
244 VibrationalCoordinateCollectionBuilder::default()
245 }
246
247 pub fn n_modes(&self) -> usize {
249 self.frequencies.len()
250 }
251
252 pub fn frequencies(&self) -> &Array1<T> {
254 &self.frequencies
255 }
256
257 pub fn coefficients(&self) -> &Array2<T> {
259 &self.coefficients
260 }
261
262 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
280impl<'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
304impl<'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
325impl<'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}