qsym2/auxiliary/
atom.rs

1//! Atoms and chemical elements.
2
3use std::collections::HashMap;
4use std::fmt;
5use std::hash::{Hash, Hasher};
6
7use approx;
8use nalgebra::{Matrix3, Point3, Rotation3, Translation3, UnitVector3, Vector3};
9use num_traits::ToPrimitive;
10use periodic_table;
11use serde::{Deserialize, Serialize};
12
13use crate::auxiliary::geometry::{self, ImproperRotationKind, Transform};
14use crate::auxiliary::misc::{self, HashableFloat};
15
16// https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
17/// Constant for converting $`a_0`$ to $`Å`$.
18pub(crate) const BOHR_TO_ANGSTROM: f64 = 0.529177210903;
19
20/// Constant for converting $`Å`$ to $`a_0`$.
21pub(crate) const ANGSTROM_TO_BOHR: f64 = 1.889726124626;
22
23/// Structure storing a look-up of element symbols to give atomic numbers
24/// and atomic masses.
25pub struct ElementMap<'a> {
26    /// A [`HashMap`] from a symbol string to a tuple of atomic number and atomic mass.
27    map: HashMap<&'a str, (u32, f64)>,
28}
29
30impl Default for ElementMap<'static> {
31    fn default() -> Self {
32        Self::new()
33    }
34}
35
36impl ElementMap<'static> {
37    /// Creates a new [`ElementMap`] for all elements in the periodic table.
38    #[must_use]
39    pub fn new() -> ElementMap<'static> {
40        let elements = periodic_table::periodic_table();
41        let map = elements
42            .iter()
43            .map(|element| {
44                let mass = parse_atomic_mass(element.atomic_mass);
45                (element.symbol, (element.atomic_number, mass))
46            })
47            .collect::<HashMap<_, _>>();
48        ElementMap { map }
49    }
50
51    /// Returns a tuple of atomic number and atomic mass for the requested element.
52    pub fn get(&self, element: &str) -> Option<&(u32, f64)> {
53        self.map.get(element)
54    }
55}
56
57/// Parses the atomic mass string in the format of [`periodic_table`] to a single float value.
58///
59/// # Arguments
60///
61/// * `mass_str` - A string of mass value that is either `x.y(z)` where the
62///     uncertain digit `z` is enclosed in parentheses, or `[x]` where `x`
63///     is the mass number in place of precise experimental values.
64///
65/// # Returns
66///
67/// The numeric mass value.
68fn parse_atomic_mass(mass_str: &str) -> f64 {
69    let mass = mass_str.replace(&['(', ')', '[', ']'][..], "");
70    mass.parse::<f64>()
71        .unwrap_or_else(|_| panic!("Unable to parse atomic mass string {mass}."))
72}
73
74/// Structure representing an atom.
75#[derive(Clone, Serialize, Deserialize)]
76pub struct Atom {
77    /// The atom kind.
78    pub kind: AtomKind,
79
80    /// The atomic number of the atom.
81    pub atomic_number: u32,
82
83    /// The atomic symbol of the atom.
84    pub atomic_symbol: String,
85
86    /// The weighted-average atomic mass for all naturally occuring isotopes.
87    pub atomic_mass: f64,
88
89    /// The position of the atom.
90    pub coordinates: Point3<f64>,
91
92    /// A threshold for approximate equality comparisons.
93    pub threshold: f64,
94}
95
96impl Atom {
97    /// Parses an atom line in an `xyz` file to construct an [`Atom`].
98    ///
99    /// # Arguments
100    ///
101    /// * `line` - A line in an `xyz` file containing an atomic symbol and
102    ///     three Cartesian coordinates.
103    /// * `emap` - A hash map between atomic symbols and atomic numbers and
104    ///     masses.
105    /// * `thresh` - A threshold for approximate equality comparisons.
106    ///
107    /// # Returns
108    ///
109    /// The parsed [`Atom`] struct if the line has the correct format,
110    /// otherwise [`None`].
111    #[must_use]
112    pub(crate) fn from_xyz(line: &str, emap: &ElementMap, thresh: f64) -> Option<Atom> {
113        let split: Vec<&str> = line.split_whitespace().collect();
114        if split.len() != 4 {
115            return None;
116        };
117        let atomic_symbol = split.first().expect("Unable to get the element symbol.");
118        let (atomic_number, atomic_mass) = emap
119            .map
120            .get(atomic_symbol)
121            .expect("Invalid atomic symbol encountered.");
122        let coordinates = Point3::new(
123            split
124                .get(1)
125                .expect("Unable to get the x coordinate.")
126                .parse::<f64>()
127                .expect("Unable to parse the x coordinate."),
128            split
129                .get(2)
130                .expect("Unable to get the y coordinate.")
131                .parse::<f64>()
132                .expect("Unable to parse the y coordinate."),
133            split
134                .get(3)
135                .expect("Unable to get the z coordinate.")
136                .parse::<f64>()
137                .expect("Unable to parse the z coordinate."),
138        );
139        let atom = Atom {
140            kind: AtomKind::Ordinary,
141            atomic_number: *atomic_number,
142            atomic_symbol: (*atomic_symbol).to_string(),
143            atomic_mass: *atomic_mass,
144            coordinates,
145            threshold: thresh,
146        };
147        Some(atom)
148    }
149
150    /// Writes the atom to an XYZ line.
151    #[must_use]
152    pub(crate) fn to_xyz(&self) -> String {
153        let precision = self.precision();
154        let length = (precision + precision.div_euclid(2)).max(6);
155        if let AtomKind::Ordinary = self.kind {
156            format!(
157                "{:<3} {:+length$.precision$} {:+length$.precision$} {:+length$.precision$}",
158                self.atomic_symbol, self.coordinates[0], self.coordinates[1], self.coordinates[2],
159            )
160        } else {
161            String::new()
162        }
163    }
164
165    /// Creates an ordinary atom.
166    ///
167    /// # Arguments
168    ///
169    /// * `atomic_symbol` - The element symbol of the atom.
170    /// * `coordinates` - The coordinates of the atom.
171    /// * `emap` - A hash map between atomic symbols and atomic numbers and
172    ///     masses.
173    /// * `thresh` - The threshold for comparing atoms.
174    ///
175    /// # Returns
176    ///
177    /// The required ordinary atom.
178    #[must_use]
179    pub fn new_ordinary(
180        atomic_symbol: &str,
181        coordinates: Point3<f64>,
182        emap: &ElementMap,
183        thresh: f64,
184    ) -> Atom {
185        let (atomic_number, atomic_mass) = emap
186            .map
187            .get(atomic_symbol)
188            .expect("Invalid atomic symbol encountered.");
189        Atom {
190            kind: AtomKind::Ordinary,
191            atomic_number: *atomic_number,
192            atomic_symbol: atomic_symbol.to_string(),
193            atomic_mass: *atomic_mass,
194            coordinates,
195            threshold: thresh,
196        }
197    }
198
199    /// Creates a special atom.
200    ///
201    /// # Arguments
202    ///
203    /// * `kind` - The required special kind.
204    /// * `coordinates` - The coordinates of the special atom.
205    /// * `thresh` - The threshold for comparing atoms.
206    ///
207    /// # Returns
208    ///
209    /// `None` if `kind` is not one of the special atom kinds, `Some(Atom)`
210    /// otherwise.
211    #[must_use]
212    pub fn new_special(kind: AtomKind, coordinates: Point3<f64>, thresh: f64) -> Option<Atom> {
213        match kind {
214            AtomKind::Magnetic(_) | AtomKind::Electric(_) => Some(Atom {
215                kind,
216                atomic_number: 0,
217                atomic_symbol: String::new(),
218                atomic_mass: 100.0,
219                coordinates,
220                threshold: thresh,
221            }),
222            AtomKind::Ordinary => None,
223        }
224    }
225
226    /// Computes the precision for printing out atom coordinates in a manner compatible with the
227    /// comparison threshold of the atom.
228    ///
229    /// # Returns
230    ///
231    /// The print-out precision.
232    fn precision(&self) -> usize {
233        self.threshold.log10().abs().round().to_usize().unwrap_or(7) + 1
234    }
235}
236
237impl fmt::Display for Atom {
238    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
239        let precision = self.precision();
240        let length = (precision + precision.div_euclid(2)).max(6);
241        if let AtomKind::Ordinary = self.kind {
242            write!(
243                f,
244                "{:>9} {:>3} {:+length$.precision$} {:+length$.precision$} {:+length$.precision$}",
245                "Atom",
246                self.atomic_symbol,
247                self.coordinates[0],
248                self.coordinates[1],
249                self.coordinates[2],
250            )
251        } else {
252            write!(
253                f,
254                "{:>13} {:+length$.precision$} {:+length$.precision$} {:+length$.precision$}",
255                self.kind, self.coordinates[0], self.coordinates[1], self.coordinates[2],
256            )
257        }
258    }
259}
260
261impl fmt::Debug for Atom {
262    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
263        write!(f, "{self}")
264    }
265}
266
267/// Enumerated type describing the atom kind.
268#[derive(Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
269pub enum AtomKind {
270    /// An ordinary atom.
271    Ordinary,
272
273    /// A fictitious atom representing a magnetic field. This variant contains
274    /// a flag to indicate if the fictitious atom is of positive type or not.
275    Magnetic(bool),
276
277    /// A fictitious atom representing an electric field. This variant contains
278    /// a flag to indicate if the fictitious atom is of positive type or not.
279    Electric(bool),
280}
281
282impl fmt::Display for AtomKind {
283    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
284        match self {
285            Self::Ordinary => write!(f, "{}", "Atom".to_owned()),
286            Self::Magnetic(pos) => {
287                if *pos {
288                    write!(f, "{}", "MagneticAtom+".to_owned())
289                } else {
290                    write!(f, "{}", "MagneticAtom-".to_owned())
291                }
292            }
293            Self::Electric(pos) => {
294                if *pos {
295                    write!(f, "{}", "ElectricAtom+".to_owned())
296                } else {
297                    write!(f, "{}", "ElectricAtom-".to_owned())
298                }
299            }
300        }
301    }
302}
303
304impl Transform for Atom {
305    fn transform_mut(&mut self, mat: &Matrix3<f64>) {
306        let det = mat.determinant();
307        assert!(
308            approx::relative_eq!(
309                det,
310                1.0,
311                epsilon = self.threshold,
312                max_relative = self.threshold
313            ) || approx::relative_eq!(
314                det,
315                -1.0,
316                epsilon = self.threshold,
317                max_relative = self.threshold
318            )
319        );
320        self.coordinates = mat * self.coordinates;
321        if approx::relative_eq!(
322            det,
323            -1.0,
324            epsilon = self.threshold,
325            max_relative = self.threshold
326        ) {
327            if let AtomKind::Magnetic(pos) = self.kind {
328                self.kind = AtomKind::Magnetic(!pos);
329            }
330        };
331    }
332
333    fn rotate_mut(&mut self, angle: f64, axis: &Vector3<f64>) {
334        let normalised_axis = UnitVector3::new_normalize(*axis);
335        let rotation = Rotation3::from_axis_angle(&normalised_axis, angle);
336        self.coordinates = rotation.transform_point(&self.coordinates);
337    }
338
339    fn improper_rotate_mut(
340        &mut self,
341        angle: f64,
342        axis: &Vector3<f64>,
343        kind: &ImproperRotationKind,
344    ) {
345        let mat = geometry::improper_rotation_matrix(angle, axis, 1, kind);
346        self.transform_mut(&mat);
347    }
348
349    fn translate_mut(&mut self, tvec: &Vector3<f64>) {
350        let translation = Translation3::from(*tvec);
351        self.coordinates = translation.transform_point(&self.coordinates);
352    }
353
354    fn recentre_mut(&mut self) {
355        self.coordinates = Point3::origin();
356    }
357
358    fn reverse_time_mut(&mut self) {
359        if let AtomKind::Magnetic(polarity) = self.kind {
360            self.kind = AtomKind::Magnetic(!polarity);
361        }
362    }
363
364    fn transform(&self, mat: &Matrix3<f64>) -> Self {
365        let mut transformed_atom = self.clone();
366        transformed_atom.transform_mut(mat);
367        transformed_atom
368    }
369
370    fn rotate(&self, angle: f64, axis: &Vector3<f64>) -> Self {
371        let mut rotated_atom = self.clone();
372        rotated_atom.rotate_mut(angle, axis);
373        rotated_atom
374    }
375
376    fn improper_rotate(
377        &self,
378        angle: f64,
379        axis: &Vector3<f64>,
380        kind: &ImproperRotationKind,
381    ) -> Self {
382        let mut improper_rotated_atom = self.clone();
383        improper_rotated_atom.improper_rotate_mut(angle, axis, kind);
384        improper_rotated_atom
385    }
386
387    fn translate(&self, tvec: &Vector3<f64>) -> Self {
388        let mut translated_atom = self.clone();
389        translated_atom.translate_mut(tvec);
390        translated_atom
391    }
392
393    fn recentre(&self) -> Self {
394        let mut recentred_atom = self.clone();
395        recentred_atom.recentre_mut();
396        recentred_atom
397    }
398
399    fn reverse_time(&self) -> Self {
400        let mut time_reversed_atom = self.clone();
401        time_reversed_atom.reverse_time_mut();
402        time_reversed_atom
403    }
404}
405
406impl PartialEq for Atom {
407    /// The `[Self::threshold]` value for each atom defines a discrete grid over
408    /// the real number field. All real numbers (*e.g.* atomic mass, coordinates)
409    /// are rounded to take on the values on this discrete grid which are then
410    /// used for hashing and comparisons.
411    fn eq(&self, other: &Self) -> bool {
412        let result = self.atomic_number == other.atomic_number
413            && self.kind == other.kind
414            && approx::relative_eq!(
415                self.atomic_mass.round_factor(self.threshold),
416                other.atomic_mass.round_factor(other.threshold),
417            )
418            && approx::relative_eq!(
419                self.coordinates[0].round_factor(self.threshold),
420                other.coordinates[0].round_factor(other.threshold),
421            )
422            && approx::relative_eq!(
423                self.coordinates[1].round_factor(self.threshold),
424                other.coordinates[1].round_factor(other.threshold),
425            )
426            && approx::relative_eq!(
427                self.coordinates[2].round_factor(self.threshold),
428                other.coordinates[2].round_factor(other.threshold),
429            );
430        // if result {
431        //     assert_eq!(misc::calculate_hash(self), misc::calculate_hash(other));
432        // }
433        result && (misc::calculate_hash(self) == misc::calculate_hash(other))
434    }
435}
436
437impl Eq for Atom {}
438
439impl Hash for Atom {
440    fn hash<H: Hasher>(&self, state: &mut H) {
441        self.atomic_number.hash(state);
442        self.kind.hash(state);
443        self.atomic_mass
444            .round_factor(self.threshold)
445            .integer_decode()
446            .hash(state);
447        self.coordinates[0]
448            .round_factor(self.threshold)
449            .integer_decode()
450            .hash(state);
451        self.coordinates[1]
452            .round_factor(self.threshold)
453            .integer_decode()
454            .hash(state);
455        self.coordinates[2]
456            .round_factor(self.threshold)
457            .integer_decode()
458            .hash(state);
459    }
460}