1use std::collections::HashSet;
4use std::fmt;
5
6use approx;
7use derive_builder::Builder;
8use fraction;
9use itertools::{self, Itertools};
10use nalgebra::{ClosedMulAssign, Matrix3, Point3, Rotation3, Scalar, UnitVector3, Vector3};
11use num_traits::{One, ToPrimitive};
12use serde::{Deserialize, Serialize};
13
14use crate::auxiliary::atom::Atom;
15use crate::auxiliary::misc::HashableFloat;
16
17type F32 = fraction::GenericFraction<u32>;
18
19#[cfg(test)]
20#[path = "geometry_tests.rs"]
21mod geometry_tests;
22
23pub enum ImproperRotationKind {
29 MirrorPlane,
32
33 InversionCentre,
36}
37
38pub const IMSIG: ImproperRotationKind = ImproperRotationKind::MirrorPlane;
40
41pub const IMINV: ImproperRotationKind = ImproperRotationKind::InversionCentre;
43
44#[must_use]
60pub fn normalise_rotation_angle(rot_ang: f64, thresh: f64) -> (f64, u32) {
61 let frac_1_2 = 1.0 / 2.0;
62 let fraction = rot_ang / (2.0 * std::f64::consts::PI);
63 if fraction > frac_1_2 + thresh {
64 let integer_part = fraction.trunc().to_u32().unwrap_or_else(|| {
65 panic!("Unable to convert the integer part of `{fraction}` to `u32`.")
66 });
67 let x = if fraction.fract() <= frac_1_2 + thresh {
68 integer_part
69 } else {
70 integer_part + 1
71 };
72 (rot_ang - 2.0 * std::f64::consts::PI * f64::from(x), x)
73 } else if fraction <= -frac_1_2 + thresh {
74 let integer_part = (-fraction).trunc().to_u32().unwrap_or_else(|| {
75 panic!("Unable to convert the integer part of `{fraction}` to `u32`.")
76 });
77 let x = if (-fraction).fract() < frac_1_2 - thresh {
78 integer_part
79 } else {
80 integer_part + 1
81 };
82 (rot_ang + 2.0 * std::f64::consts::PI * f64::from(x), x)
83 } else {
84 (rot_ang, 0)
85 }
86}
87
88#[must_use]
99pub fn normalise_rotation_fraction(fraction: F32) -> (F32, u32) {
100 let frac_1_2 = F32::new(1u32, 2u32);
118 if fraction > frac_1_2 {
119 let integer_part = fraction.trunc();
120 let x = if fraction.fract() <= frac_1_2 {
121 integer_part
122 } else {
123 integer_part + F32::one()
124 };
125 (
126 fraction - x,
127 x.to_u32()
128 .expect("Unable to convert the 2π-turn number to `u32`."),
129 )
130 } else if fraction <= -frac_1_2 {
131 let integer_part = (-fraction).trunc();
132 let x = if (-fraction).fract() < frac_1_2 {
133 integer_part
134 } else {
135 integer_part + F32::one()
136 };
137 (
138 fraction + x,
139 x.to_u32()
140 .expect("Unable to convert the 2π-turn number to `u32`."),
141 )
142 } else {
143 (fraction, 0)
144 }
145}
146
147#[must_use]
164pub fn get_proper_fraction(angle: f64, thresh: f64, max_trial_power: u32) -> Option<F32> {
165 let (normalised_angle, _) = normalise_rotation_angle(angle, thresh);
166 let rational_order = (2.0 * std::f64::consts::PI) / normalised_angle.abs();
167 let mut power: u32 = 1;
168 while approx::relative_ne!(
169 rational_order * (f64::from(power)),
170 (rational_order * (f64::from(power))).round(),
171 max_relative = thresh,
172 epsilon = thresh
173 ) && power < max_trial_power
174 {
175 power += 1;
176 }
177 if approx::relative_eq!(
178 rational_order * (f64::from(power)),
179 (rational_order * (f64::from(power))).round(),
180 max_relative = thresh,
181 epsilon = thresh
182 ) {
183 let orderf64 = (rational_order * (f64::from(power))).round();
184 assert!(orderf64.is_sign_positive());
185 assert!(orderf64 <= f64::from(u32::MAX));
186 #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
187 let order = orderf64 as u32;
188 if normalised_angle > 0.0 {
189 Some(F32::new(power, order))
190 } else {
191 Some(F32::new_neg(power, order))
192 }
193 } else {
194 None
195 }
196}
197
198fn outer<T: Scalar + ClosedMulAssign + Copy>(vec1: &Vector3<T>, vec2: &Vector3<T>) -> Matrix3<T> {
209 let outer_product_iter: Vec<T> = vec2
210 .iter()
211 .flat_map(|&item_x| vec1.iter().map(move |&item_y| item_x * item_y))
212 .collect();
213 Matrix3::from_iterator(outer_product_iter)
214}
215
216#[must_use]
229pub fn proper_rotation_matrix(angle: f64, axis: &Vector3<f64>, power: i8) -> Matrix3<f64> {
230 let normalised_axis = UnitVector3::new_normalize(*axis);
231 Rotation3::from_axis_angle(&normalised_axis, (f64::from(power)) * angle).into_inner()
232}
233
234#[must_use]
248pub fn improper_rotation_matrix(
249 angle: f64,
250 axis: &Vector3<f64>,
251 power: i8,
252 kind: &ImproperRotationKind,
253) -> Matrix3<f64> {
254 let rotmat = proper_rotation_matrix(angle, axis, power);
255 let normalised_axis = UnitVector3::new_normalize(*axis);
256 match kind {
257 ImproperRotationKind::MirrorPlane => {
258 let refmat = Matrix3::identity()
259 - 2.0 * (f64::from(power % 2)) * outer(&normalised_axis, &normalised_axis);
260 refmat * rotmat
261 }
262 ImproperRotationKind::InversionCentre => {
263 if power % 2 == 1 {
264 -rotmat
265 } else {
266 rotmat
267 }
268 }
269 }
270}
271
272#[must_use]
286pub fn check_regular_polygon(atoms: &[&Atom]) -> bool {
287 assert!(
288 atoms.len() >= 3,
289 "Polygons can only be formed by three atoms or more."
290 );
291
292 let tot_m: f64 = atoms.iter().fold(0.0, |acc, atom| acc + atom.atomic_mass);
293 let com: Point3<f64> = atoms.iter().fold(Point3::origin(), |acc, atom| {
294 acc + (atom.coordinates * atom.atomic_mass - Point3::origin())
295 }) / tot_m;
296
297 let radial_dists: HashSet<(u64, i16, i8)> = atoms
298 .iter()
299 .map(|atom| {
300 (atom.coordinates - com)
301 .norm()
302 .round_factor(atom.threshold)
303 .integer_decode()
304 })
305 .collect();
306
307 if radial_dists.len() == 1 {
309 let regular_angle = 2.0 * std::f64::consts::PI
310 / atoms
311 .len()
312 .to_f64()
313 .unwrap_or_else(|| panic!("Unable to convert `{}` to `f64`.", atoms.len()));
314 let thresh = atoms
315 .iter()
316 .fold(0.0_f64, |acc, atom| acc.max(atom.threshold));
317 let mut rad_vectors: Vec<Vector3<f64>> =
318 atoms.iter().map(|atom| atom.coordinates - com).collect();
319 let (vec_i, vec_j) = itertools::iproduct!(rad_vectors.iter(), rad_vectors.iter())
320 .max_by(|&(v_i1, v_j1), &(v_i2, v_j2)| {
321 v_i1.cross(v_j1)
322 .norm()
323 .partial_cmp(&v_i2.cross(v_j2).norm())
324 .expect("Unable to compare the cross products of two vector pairs.")
325 })
326 .expect("Unable to find the vector pair with the largest norm cross product.");
327 let normal = UnitVector3::new_normalize(vec_i.cross(vec_j));
328 if normal.norm() < thresh {
329 return false;
330 }
331
332 let vec0 = atoms[0].coordinates - com;
333 rad_vectors.sort_by(|a, b| {
334 get_anticlockwise_angle(&vec0, a, &normal, thresh)
335 .partial_cmp(&get_anticlockwise_angle(&vec0, b, &normal, thresh))
336 .unwrap_or_else(|| {
337 panic!(
338 "Unable to compare anticlockwise angles of {a} and {b} relative to {vec0}."
339 )
340 })
341 });
342 let vector_pairs: Vec<(&Vector3<f64>, &Vector3<f64>)> =
343 rad_vectors.iter().circular_tuple_windows().collect();
344 let mut angles: HashSet<(u64, i16, i8)> = vector_pairs
345 .iter()
346 .map(|(v1, v2)| {
347 get_anticlockwise_angle(v1, v2, &normal, thresh)
348 .round_factor(thresh)
349 .integer_decode()
350 })
351 .collect();
352 angles.insert(regular_angle.round_factor(thresh).integer_decode());
353
354 angles.len() == 1
355 } else {
356 false
357 }
358}
359
360fn get_anticlockwise_angle(
377 vec1: &Vector3<f64>,
378 vec2: &Vector3<f64>,
379 normal: &UnitVector3<f64>,
380 thresh: f64,
381) -> f64 {
382 assert!(thresh >= std::f64::EPSILON);
383 assert!(vec1.norm() >= thresh);
384 assert!(vec2.norm() >= thresh);
385 let dot = vec1.dot(vec2);
386 let det = normal.into_inner().dot(&vec1.cross(vec2));
387 let mut angle = det.atan2(dot);
388 while angle < -thresh {
389 angle += 2.0 * std::f64::consts::PI;
390 }
391 angle
392}
393
394pub trait Transform {
396 fn transform_mut(&mut self, mat: &Matrix3<f64>);
403
404 fn rotate_mut(&mut self, angle: f64, axis: &Vector3<f64>);
411
412 fn improper_rotate_mut(&mut self, angle: f64, axis: &Vector3<f64>, kind: &ImproperRotationKind);
420
421 fn translate_mut(&mut self, tvec: &Vector3<f64>);
428
429 fn recentre_mut(&mut self);
431
432 fn reverse_time_mut(&mut self);
434
435 #[must_use]
446 fn transform(&self, mat: &Matrix3<f64>) -> Self;
447
448 #[must_use]
459 fn rotate(&self, angle: f64, axis: &Vector3<f64>) -> Self;
460
461 #[must_use]
473 fn improper_rotate(&self, angle: f64, axis: &Vector3<f64>, kind: &ImproperRotationKind)
474 -> Self;
475
476 #[must_use]
487 fn translate(&self, tvec: &Vector3<f64>) -> Self;
488
489 #[must_use]
495 fn recentre(&self) -> Self;
496
497 #[must_use]
504 fn reverse_time(&self) -> Self;
505}
506
507#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
517enum ImproperOrdering {
518 Greater,
519 GreaterEqual,
520 Less,
521 LessEqual,
522 Equal,
523}
524
525impl fmt::Display for ImproperOrdering {
526 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
527 match self {
528 ImproperOrdering::Greater => write!(f, ">"),
529 ImproperOrdering::GreaterEqual => write!(f, "≥"),
530 ImproperOrdering::Less => write!(f, "<"),
531 ImproperOrdering::LessEqual => write!(f, "≤"),
532 ImproperOrdering::Equal => write!(f, "="),
533 }
534 }
535}
536
537#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
547enum CartesianCoordinate {
548 X,
549 Y,
550 Z,
551}
552
553impl CartesianCoordinate {
554 fn to_index(&self) -> usize {
556 match self {
557 CartesianCoordinate::X => 0,
558 CartesianCoordinate::Y => 1,
559 CartesianCoordinate::Z => 2,
560 }
561 }
562}
563
564impl fmt::Display for CartesianCoordinate {
565 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
566 match self {
567 CartesianCoordinate::X => write!(f, "x"),
568 CartesianCoordinate::Y => write!(f, "y"),
569 CartesianCoordinate::Z => write!(f, "z"),
570 }
571 }
572}
573
574#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
580pub struct CartesianConditions {
581 conditions: Vec<Vec<(CartesianCoordinate, ImproperOrdering, f64)>>,
584}
585
586impl CartesianConditions {
587 fn check(&self, vec: &Vector3<f64>, thresh: f64) -> bool {
599 self.conditions.iter().any(|condition_set| {
600 condition_set.iter().all(|(i, order, target)| match order {
601 ImproperOrdering::Greater => vec[i.to_index()] > target + thresh,
602 ImproperOrdering::GreaterEqual => vec[i.to_index()] > target - thresh,
603 ImproperOrdering::Less => vec[i.to_index()] < target - thresh,
604 ImproperOrdering::LessEqual => vec[i.to_index()] < target + thresh,
605 ImproperOrdering::Equal => approx::relative_eq!(
606 vec[i.to_index()],
607 target,
608 max_relative = thresh,
609 epsilon = thresh
610 ),
611 })
612 })
613 }
614}
615
616impl Default for CartesianConditions {
617 fn default() -> Self {
618 Self {
619 conditions: vec![
620 vec![(CartesianCoordinate::Z, ImproperOrdering::Greater, 0.0)],
621 vec![
622 (CartesianCoordinate::Z, ImproperOrdering::Equal, 0.0),
623 (CartesianCoordinate::X, ImproperOrdering::Greater, 0.0),
624 ],
625 vec![
626 (CartesianCoordinate::Z, ImproperOrdering::Equal, 0.0),
627 (CartesianCoordinate::X, ImproperOrdering::Equal, 0.0),
628 (CartesianCoordinate::Y, ImproperOrdering::Greater, 0.0),
629 ],
630 ],
631 }
632 }
633}
634
635impl fmt::Display for CartesianConditions {
636 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
637 writeln!(f, "Cartesian conditions:")?;
638 let conditions = self
639 .conditions
640 .iter()
641 .map(|condition_set| {
642 condition_set
643 .iter()
644 .map(|(i, order, target)| format!("{i} {order} {target}"))
645 .join(", ")
646 })
647 .join("\n or\n");
648 writeln!(f, "{conditions}")?;
649 Ok(())
650 }
651}
652
653#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
663pub enum SphericalCoordinate {
664 Theta,
665 Phi,
666}
667
668impl fmt::Display for SphericalCoordinate {
669 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
670 match self {
671 SphericalCoordinate::Theta => write!(f, "θ"),
672 SphericalCoordinate::Phi => write!(f, "φ"),
673 }
674 }
675}
676
677#[derive(Debug, Clone, Builder, PartialEq, Serialize, Deserialize)]
683pub struct SphericalConditions {
684 #[builder(setter(custom))]
686 z_basis: Vector3<f64>,
687
688 #[builder(setter(custom))]
690 x_basis: Vector3<f64>,
691
692 #[builder(setter(custom))]
695 conditions: Vec<Vec<(SphericalCoordinate, ImproperOrdering, f64)>>,
696}
697
698impl SphericalConditionsBuilder {
699 fn z_basis(&mut self, z_bas: Vector3<f64>) -> &mut Self {
700 self.z_basis = Some(z_bas.normalize());
701 self
702 }
703
704 fn x_basis(&mut self, x_bas: Vector3<f64>) -> &mut Self {
705 self.x_basis = Some(x_bas.normalize());
706 self
707 }
708
709 fn conditions(
710 &mut self,
711 conds: &[Vec<(SphericalCoordinate, ImproperOrdering, f64)>],
712 ) -> &mut Self {
713 self.conditions = Some(conds.to_vec());
714 self
715 }
716}
717
718impl SphericalConditions {
719 fn builder() -> SphericalConditionsBuilder {
721 SphericalConditionsBuilder::default()
722 }
723
724 fn get_component(&self, vec: &Vector3<f64>, coord: &SphericalCoordinate, thresh: f64) -> f64 {
737 match coord {
738 SphericalCoordinate::Theta => {
739 let cos_theta = vec.dot(&self.z_basis) / (vec.norm() * self.z_basis.norm());
740 if approx::relative_eq!(cos_theta, 1.0, epsilon = thresh, max_relative = thresh) {
741 1.0f64.acos()
742 } else if approx::relative_eq!(
743 cos_theta,
744 -1.0,
745 epsilon = thresh,
746 max_relative = thresh
747 ) {
748 (-1.0f64).acos()
749 } else {
750 cos_theta.acos()
751 }
752 }
753 SphericalCoordinate::Phi => {
754 let y_vector = self.z_basis.cross(&self.x_basis).normalize();
755 let xy_vec = vec - vec.dot(&self.z_basis) / self.z_basis.norm() * self.z_basis;
756 let sgn_y = xy_vec.dot(&y_vector).signum();
757 let cos_phi = xy_vec.dot(&self.x_basis) / (xy_vec.norm() * self.x_basis.norm());
758 if approx::relative_eq!(cos_phi, 1.0, epsilon = thresh, max_relative = thresh) {
759 sgn_y * 1.0f64.acos()
760 } else if approx::relative_eq!(
761 cos_phi,
762 -1.0,
763 epsilon = thresh,
764 max_relative = thresh
765 ) {
766 sgn_y * (-1.0f64).acos()
767 } else {
768 sgn_y * cos_phi.acos()
769 }
770 }
771 }
772 }
773
774 fn check(&self, vec: &Vector3<f64>, thresh: f64) -> bool {
786 self.conditions.iter().any(|condition_set| {
787 condition_set.iter().all(|(i, order, target)| {
788 let component = self.get_component(vec, i, thresh);
789 match order {
790 ImproperOrdering::Greater => component > target + thresh,
791 ImproperOrdering::GreaterEqual => component > target - thresh,
792 ImproperOrdering::Less => component < target - thresh,
793 ImproperOrdering::LessEqual => component < target + thresh,
794 ImproperOrdering::Equal => approx::relative_eq!(
795 component,
796 target,
797 max_relative = thresh,
798 epsilon = thresh
799 ),
800 }
801 })
802 })
803 }
804
805 fn new_disjoint_equatorial_arcs(
840 z_basis: Vector3<f64>,
841 x_basis: Vector3<f64>,
842 n: usize,
843 ) -> Self {
844 assert!(n > 0 && n.rem_euclid(2) == 1);
845 let n_f64 = n
846 .to_f64()
847 .expect("Unable to convert the number of arcs to `f64`.");
848 let half_arc = std::f64::consts::PI / (2.0 * n_f64);
849 let sep = 2.0 * std::f64::consts::PI / n_f64;
850 let half_pi = 0.5 * std::f64::consts::PI;
851
852 let mut conditions = vec![vec![
853 (
854 SphericalCoordinate::Theta,
855 ImproperOrdering::GreaterEqual,
856 0.0,
857 ),
858 (SphericalCoordinate::Theta, ImproperOrdering::Less, half_pi),
859 ]];
860
861 let phi_conditions = (0..n)
862 .map(|i| {
863 let (centre, _) = normalise_rotation_angle(i.to_f64().unwrap() * sep, f64::EPSILON);
864 let min_exc = centre - half_arc;
865 let max_inc = centre + half_arc;
866 vec![
867 (SphericalCoordinate::Theta, ImproperOrdering::Equal, half_pi),
868 (SphericalCoordinate::Phi, ImproperOrdering::Greater, min_exc),
869 (
870 SphericalCoordinate::Phi,
871 ImproperOrdering::LessEqual,
872 max_inc,
873 ),
874 ]
875 })
876 .collect_vec();
877
878 conditions.extend(phi_conditions.into_iter());
879 Self::builder()
880 .z_basis(z_basis)
881 .x_basis(x_basis)
882 .conditions(&conditions)
883 .build()
884 .expect("Unable to construct a set of spherical-coordinate conditions.")
885 }
886}
887
888impl Default for SphericalConditions {
889 fn default() -> Self {
890 let half_pi = 0.5 * std::f64::consts::PI;
891 let conditions = vec![
892 vec![
893 (
894 SphericalCoordinate::Theta,
895 ImproperOrdering::GreaterEqual,
896 0.0,
897 ),
898 (SphericalCoordinate::Theta, ImproperOrdering::Less, half_pi),
899 ],
900 vec![
901 (SphericalCoordinate::Theta, ImproperOrdering::Equal, half_pi),
902 (
903 SphericalCoordinate::Phi,
904 ImproperOrdering::Greater,
905 -half_pi,
906 ),
907 (
908 SphericalCoordinate::Phi,
909 ImproperOrdering::LessEqual,
910 half_pi,
911 ),
912 ],
913 ];
914 Self::builder()
915 .z_basis(Vector3::z())
916 .x_basis(Vector3::x())
917 .conditions(&conditions)
918 .build()
919 .expect("Unable to construct a set of spherical-coordinate conditions.")
920 }
921}
922
923impl fmt::Display for SphericalConditions {
924 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
925 writeln!(f, "Spherical conditions:")?;
926 writeln!(f, "Polar axis (z-like): {:?}", self.z_basis)?;
927 writeln!(f, "Azimuthal axis (x-like): {:?}", self.x_basis)?;
928 let conditions = self
929 .conditions
930 .iter()
931 .map(|condition_set| {
932 condition_set
933 .iter()
934 .map(|(i, order, target)| format!("{i} {order} {target}"))
935 .join(", ")
936 })
937 .join("\n or\n");
938 writeln!(f, "{conditions}")?;
939 Ok(())
940 }
941}
942
943#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
949pub enum PositiveHemisphere {
950 Cartesian(CartesianConditions),
951 Spherical(SphericalConditions),
952}
953
954impl PositiveHemisphere {
955 pub fn new_standard_cartesian() -> Self {
957 Self::Cartesian(CartesianConditions::default())
958 }
959
960 pub fn new_standard_spherical() -> Self {
962 Self::Spherical(SphericalConditions::default())
963 }
964
965 pub fn new_spherical_disjoint_equatorial_arcs(
979 z_basis: Vector3<f64>,
980 x_basis: Vector3<f64>,
981 n: usize,
982 ) -> Self {
983 Self::Spherical(SphericalConditions::new_disjoint_equatorial_arcs(
984 z_basis, x_basis, n,
985 ))
986 }
987
988 pub fn check_positive_pole(&self, axis: &Vector3<f64>, thresh: f64) -> bool {
1003 let normalised_axis = axis.normalize();
1004 match self {
1005 PositiveHemisphere::Cartesian(cart_conditions) => {
1006 cart_conditions.check(&normalised_axis, thresh)
1007 }
1008 PositiveHemisphere::Spherical(sph_conditions) => {
1009 sph_conditions.check(&normalised_axis, thresh)
1010 }
1011 }
1012 }
1013
1014 pub fn get_positive_pole(&self, axis: &Vector3<f64>, thresh: f64) -> Vector3<f64> {
1030 let normalised_axis = axis.normalize();
1031 if self.check_positive_pole(&normalised_axis, thresh) {
1032 normalised_axis
1033 } else {
1034 -normalised_axis
1035 }
1036 }
1037}
1038
1039impl Default for PositiveHemisphere {
1040 fn default() -> Self {
1041 Self::Cartesian(CartesianConditions::default())
1042 }
1043}
1044
1045impl fmt::Display for PositiveHemisphere {
1046 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
1047 match self {
1048 PositiveHemisphere::Cartesian(cart_conds) => write!(f, "{cart_conds}"),
1049 PositiveHemisphere::Spherical(sph_conds) => write!(f, "{sph_conds}"),
1050 }
1051 }
1052}
1053
1054#[must_use]
1072pub fn get_standard_positive_pole(axis: &Vector3<f64>, thresh: f64) -> Vector3<f64> {
1073 let poshem = PositiveHemisphere::new_standard_cartesian();
1074 poshem.get_positive_pole(axis, thresh)
1075}
1076
1077#[must_use]
1095pub fn check_standard_positive_pole(axis: &Vector3<f64>, thresh: f64) -> bool {
1096 let poshem = PositiveHemisphere::new_standard_cartesian();
1097 poshem.check_positive_pole(axis, thresh)
1098}