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