1use 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
16pub(crate) const BOHR_TO_ANGSTROM: f64 = 0.529177210903;
19
20pub(crate) const ANGSTROM_TO_BOHR: f64 = 1.889726124626;
22
23pub struct ElementMap<'a> {
26 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 #[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 pub fn get(&self, element: &str) -> Option<&(u32, f64)> {
53 self.map.get(element)
54 }
55}
56
57fn 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#[derive(Clone, Serialize, Deserialize)]
76pub struct Atom {
77 pub kind: AtomKind,
79
80 pub atomic_number: u32,
82
83 pub atomic_symbol: String,
85
86 pub atomic_mass: f64,
88
89 pub coordinates: Point3<f64>,
91
92 pub threshold: f64,
94}
95
96impl Atom {
97 #[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 #[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 #[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 #[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 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#[derive(Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
269pub enum AtomKind {
270 Ordinary,
272
273 Magnetic(bool),
276
277 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 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 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}