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.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 pub fn get(&self, element: &str) -> Option<&(u32, f64)> {
51 self.map.get(element)
52 }
53}
54
55fn 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#[derive(Clone, Serialize, Deserialize)]
74pub struct Atom {
75 pub kind: AtomKind,
77
78 pub atomic_number: u32,
80
81 pub atomic_symbol: String,
83
84 pub atomic_mass: f64,
86
87 pub coordinates: Point3<f64>,
89
90 pub threshold: f64,
92}
93
94impl Atom {
95 #[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 #[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 #[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 #[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 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#[derive(Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
267pub enum AtomKind {
268 Ordinary,
270
271 Magnetic(bool),
274
275 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 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 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}