feat(rings): identify aromatic interactions (#6)

* feat(ring): add fn to detect ring center and normal with SVD

* feat(ring): detect rings during struct initialization

* fix(ring): the U matrix should be used for finding the normal vector

Our input matrix is 3xN instead of Nx3

* refactor(ring): extract center and normal into struct

* feat(ring): find cation-pi interactions

* feat(ring): let main fn find cation-pi interactions

* feat(ring): find cation-pi and pi-pi interactions

* fix(ring): wrong vectors in both angle calculations

* feat(ring): use pi-pi nomenclature by Chakrabarti and Bhattacharyya

* docs(readme): finish aromatic interactions

Cargo.toml

@@ -8,6 +8,7 @@ edition = "2021"
 
 [dependencies]
 clap = { version = "4.5.4", features = ["derive"] }
+nalgebra = "0.32.5"
 pdbtbx = { version = "0.11.0", features = ["rayon", "rstar"] }
 rayon = "1.10.0"
 rstar = "0.12.0"

README.md

@@ -12,12 +12,12 @@ This is a port of the [Arpeggio](https://github.com/PDBeurope/arpeggio/) library
   - [x] Steric clashes
   - [x] VdW interactions
   - [x] Hydrophobic interactions
-  - [ ] Aromatic interactions
-  - [ ] Cation-pi interactions
+  - [x] Aromatic interactions
+  - [x] Cation-pi interactions
   - [x] Ionic interactions
   - [x] Hydrogen bonds
   - [x] Weak hydrogen bonds
-  - [ ] Disulfide bonds
+  - [x] Disulfide bonds
   - [ ] Covalent bonds
 - [ ] Output results in various formats (e.g., JSON, CSV)
 - [ ] Bundle into a `PyO3` extension module

src/interactions/aromatic.rs

@@ -0,0 +1,87 @@
+use super::ionic::is_pos_ionizable;
+use super::structs::Interaction;
+use crate::residues::Ring;
+
+use nalgebra as na;
+use pdbtbx::*;
+
+const CATION_PI_ANGLE_THRESHOLD: f64 = 30.0;
+const CATION_PI_DIST_THRESHOLD: f64 = 4.5;
+const PI_PI_DIST_THRESHOLD: f64 = 6.0;
+
+/// Identify cation-pi interactions.
+pub fn find_cation_pi(ring: &Ring, entity: &AtomConformerResidueChainModel) -> Option<Interaction> {
+    if is_pos_ionizable(entity.residue().name().unwrap(), entity.atom().name()) {
+        let atom_coord = entity.atom().pos();
+        let atom_point = na::Vector3::new(atom_coord.0, atom_coord.1, atom_coord.2);
+        let dist = point_ring_dist(ring, &atom_coord);
+        let theta = point_ring_angle(ring, &atom_point);
+
+        if (theta <= CATION_PI_ANGLE_THRESHOLD) & (dist <= CATION_PI_DIST_THRESHOLD) {
+            Some(Interaction::CationPi)
+        } else {
+            None
+        }
+    } else {
+        None
+    }
+}
+
+/// Identify pi-pi interactions using the classification by [Chakrabarti and Bhattacharyya (2007)](https://doi.org/10.1016/j.pbiomolbio.2007.03.016).
+pub fn find_pi_pi(ring1: &Ring, ring2: &Ring) -> Option<Interaction> {
+    let angle_vec = ring1.center - ring2.center;
+    let dist = (angle_vec).norm();
+    if dist <= PI_PI_DIST_THRESHOLD {
+        let theta = point_ring_angle(ring1, &angle_vec);
+        let dihedral = ring_ring_angle(ring1, ring2);
+
+        match dihedral {
+            d if d <= 30.0 => match theta {
+                t if t <= 30.0 => Some(Interaction::PiSandwichStacking), // ff
+                t if t <= 60.0 => Some(Interaction::PiDisplacedStacking), // of
+                t if t <= 90.0 => Some(Interaction::PiParallelInPlaneStacking), // ee
+                _ => None,
+            },
+            // ft, ot, and et
+            d if d <= 60.0 => Some(Interaction::PiTiltedStacking),
+            d if d <= 90.0 => match theta {
+                t if (30.0..60.0).contains(&t) => Some(Interaction::PiLStacking), // oe
+                _ => Some(Interaction::PiTStacking),                              // fe and ef
+            },
+            _ => None,
+        }
+    } else {
+        None
+    }
+}
+
+/// Calculate the distance from the point to the ring center.
+pub fn point_ring_dist(ring: &Ring, point: &(f64, f64, f64)) -> f64 {
+    let atom_point = na::Vector3::new(point.0, point.1, point.2);
+    (atom_point - ring.center).norm()
+}
+
+/// Calculate the angle between the ring normal and the vector pointing from
+/// the ring center to the point.
+fn point_ring_angle(ring: &Ring, point: &na::Vector3<f64>) -> f64 {
+    let v = point - ring.center;
+    let mut rad = (ring.normal.dot(&v) / (ring.normal.norm() * v.norm())).acos();
+
+    // Convert to degrees
+    if rad > std::f64::consts::FRAC_PI_2 {
+        rad = std::f64::consts::PI - rad;
+    }
+    rad.to_degrees()
+}
+
+/// Calculate the angle between two rings.
+fn ring_ring_angle(ring1: &Ring, ring2: &Ring) -> f64 {
+    let mut rad =
+        (ring1.normal.dot(&ring2.normal) / (ring1.normal.norm() * ring2.normal.norm())).acos();
+
+    // Convert to degrees
+    if rad > std::f64::consts::FRAC_PI_2 {
+        rad = std::f64::consts::PI - rad;
+    }
+    rad.to_degrees()
+}