1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
//! # Dihedral //! //! [![crates.io](https://img.shields.io/crates/d/dihedral.svg)](https://crates.io/crates/dihedral) //! [![crates.io](https://img.shields.io/crates/v/dihedral.svg)](https://crates.io/crates/dihedral) //! [![crates.io](https://img.shields.io/crates/l/dihedral.svg)](https://crates.io/crates/dihedral) //! //! This crate provides functions for working with dihedral angles. Currently, there are two functions: //! //! - `dihedral` calculates the dihedral angle in the range -π to π in accordance with biochemistry textbooks (see also: https://en.wikipedia.org/wiki/Dihedral_angle#In_stereochemistry) //! - `dihedral_unsigned` ignores the direction of rotation and outputs the angle within the range 0 to π. This function is faster than the above signed version. //! //! If you want to use `f32` instead of `f64` for calculation, you can add `dihedral = {version = "*", features = ["f32"]}` to your `Cargo.toml`. //! //! # References //! //! - https://math.stackexchange.com/a/47084 //! - https://en.wikipedia.org/wiki/Dihedral_angle#In_stereochemistry #[cfg(not(feature = "f32"))] type Float = f64; #[cfg(feature = "f32")] type Float = f32; type Vector3D = [Float; 3]; type Point3D = [Float; 3]; /// Calculates the dihedral angle, in the range -π to π, of the four ordered coordinates. /// /// This is in accordance with the definition of dihedral angle in biochemistry textbooks. /// /// # Examples /// /// ``` /// use dihedral::dihedral; /// /// let P0 = [24.969, 13.428, 30.692]; // N /// let P1 = [24.044, 12.661, 29.808]; // CA /// let P2 = [22.785, 13.482, 29.543]; // C /// let P3 = [21.951, 13.670, 30.431]; // O /// let P4 = [23.672, 11.328, 30.466]; // CB /// let P5 = [22.881, 10.326, 29.620]; // CG /// let P6 = [23.691, 9.935, 28.389]; // CD1 /// let P7 = [22.557, 9.096, 30.459]; // CD2 /// /// assert!((dihedral([P0, P1, P2, P3]).to_degrees() - (-71.21515)).abs() < 1E-2); /// assert!((dihedral([P0, P1, P4, P5]).to_degrees() - (-171.94319)).abs() < 1E-2); /// assert!((dihedral([P1, P4, P5, P6]).to_degrees() - (60.82226)).abs() < 1E-2); /// assert!((dihedral([P1, P4, P5, P7]).to_degrees() - (-177.63641)).abs() < 1E-2); /// ``` #[inline] pub fn dihedral([a, b, c, d]: [Point3D; 4]) -> Float { let (a, b, c) = (v(a, b), v(b, c), v(c, d)); let (r, s) = (u(cross(a, b)), u(cross(b, c))); let t = cross(r, u(b)); let (x, y) = (dot(r, s), dot(s, t)); -y.atan2(x) } /// Calculates the unsigned dihedral angle, in the range 0 to π, of the four ordered coordinates /// /// # Examples /// /// ``` /// use dihedral::dihedral_unsigned; /// /// let P0 = [24.969, 13.428, 30.692]; // N /// let P1 = [24.044, 12.661, 29.808]; // CA /// let P2 = [22.785, 13.482, 29.543]; // C /// let P3 = [21.951, 13.670, 30.431]; // O /// let P4 = [23.672, 11.328, 30.466]; // CB /// let P5 = [22.881, 10.326, 29.620]; // CG /// let P6 = [23.691, 9.935, 28.389]; // CD1 /// let P7 = [22.557, 9.096, 30.459]; // CD2 /// /// assert!((dihedral_unsigned([P0, P1, P2, P3]).to_degrees() - (71.21515)).abs() < 1E-2); /// assert!((dihedral_unsigned([P0, P1, P4, P5]).to_degrees() - (171.94319)).abs() < 1E-2); /// assert!((dihedral_unsigned([P1, P4, P5, P6]).to_degrees() - (60.82226)).abs() < 1E-2); /// assert!((dihedral_unsigned([P1, P4, P5, P7]).to_degrees() - (177.63641)).abs() < 1E-2); /// ``` #[inline] pub fn dihedral_unsigned([a, b, c, d]: [Point3D; 4]) -> Float { let (a, b, c) = (v(a, b), v(b, c), v(c, d)); let (r, s) = (cross(a, b), cross(b, c)); (dot(r, s) / (norm(r) * norm(s))).acos() } #[inline(always)] fn norm(a: Vector3D) -> Float { dot(a, a).sqrt() } #[inline(always)] fn u(v: Vector3D) -> Vector3D { let norm = norm(v); [v[0] / norm, v[1] / norm, v[2] / norm] } #[inline(always)] fn v(p1: Point3D, p2: Point3D) -> Vector3D { [p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]] } #[inline(always)] fn dot(m: Vector3D, n: Vector3D) -> Float { m[0] * n[0] + m[1] * n[1] + m[2] * n[2] } #[inline(always)] fn cross(m: Vector3D, n: Vector3D) -> Vector3D { [ m[1] * n[2] - m[2] * n[1], m[2] * n[0] - m[0] * n[2], m[0] * n[1] - m[1] * n[0], ] } #[cfg(test)] mod tests { use super::*; const P0: Point3D = [24.969, 13.428, 30.692]; // N const P1: Point3D = [24.044, 12.661, 29.808]; // CA const P2: Point3D = [22.785, 13.482, 29.543]; // C const P3: Point3D = [21.951, 13.670, 30.431]; // O const P4: Point3D = [23.672, 11.328, 30.466]; // CB const P5: Point3D = [22.881, 10.326, 29.620]; // CG const P6: Point3D = [23.691, 9.935, 28.389]; // CD1 const P7: Point3D = [22.557, 9.096, 30.459]; // CD2 #[test] fn test_dihedral() { assert!((dihedral([P0, P1, P2, P3]).to_degrees() - (-71.21515)).abs() < 1E-2); assert!((dihedral([P0, P1, P4, P5]).to_degrees() - (-171.94319)).abs() < 1E-2); assert!((dihedral([P1, P4, P5, P6]).to_degrees() - (60.82226)).abs() < 1E-2); assert!((dihedral([P1, P4, P5, P7]).to_degrees() - (-177.63641)).abs() < 1E-2); } #[test] fn test_dihedral_unsigned() { println!("{}", dihedral_unsigned([P0, P1, P4, P5]).to_degrees()); assert!((dihedral_unsigned([P0, P1, P2, P3]).to_degrees() - (71.21515)).abs() < 1E-2); assert!((dihedral_unsigned([P0, P1, P4, P5]).to_degrees() - (171.94319)).abs() < 1E-2); assert!((dihedral_unsigned([P1, P4, P5, P6]).to_degrees() - (60.82226)).abs() < 1E-2); assert!((dihedral_unsigned([P1, P4, P5, P7]).to_degrees() - (177.63641)).abs() < 1E-2); } }