Added bonds infer function for #137

src/structs/pdb.rs

@@ -1073,6 +1073,84 @@ impl<'a> PDB {
         }
         names
     }
+
+    /// This function infers bonds between protein components. Most macromolecule PDB/CIF files don't include
+    /// explicit bond information. This functions infers bond lengths by comparing each interactomic
+    /// bond distance, and matching against known amino acid bond lengths.
+    ///
+    /// Some info here: https://www.ruppweb.org/Xray/tutorial/protein_structure.htm
+    /// https://itp.uni-frankfurt.de/~engel/amino.html
+    pub fn connect_atoms(&mut self) {
+        // Written by David-OConnor in issue #137
+
+        // Peptide
+        // Double bond len of C' to N.
+        const LEN_CP_N: f64 = 1.33;
+        const LEN_N_CALPHA: f64 = 1.46;
+        const LEN_CALPHA_CP: f64 = 1.53;
+
+        // Single bond
+        const LEN_C_C: f64 = 1.54;
+        const LEN_C_N: f64 = 1.48;
+        const LEN_C_O: f64 = 1.43;
+
+        // Hydrogen
+        const LEN_OH_OH: f64 = 2.8;
+        const LEN_NH_OC: f64 = 2.9;
+        const LEN_OH_OC: f64 = 2.8;
+
+        // Bonds to H. Mostly ~1
+        const LEN_N_H: f64 = 1.00;
+        const LEN_C_H: f64 = 1.10;
+        const LEN_O_H: f64 = 1.0; // In water molecules. What is it in proteins?
+
+        // If interatomic distance is within this distance of one of our known bond lengths, consider it to be a bond.
+        const BOND_LEN_THRESH: f64 = 0.04; // todo: Adjust A/R based on performance.
+        const GRID_SIZE: f64 = 3.0; // Slightly larger than the largest bond threshold
+
+        // Infer bonds from atom distances. Uses spacial partitioning for efficiency.
+        let lens_covalent = vec![
+            LEN_CP_N,
+            LEN_N_CALPHA,
+            LEN_CALPHA_CP,
+            LEN_C_C,
+            LEN_C_N,
+            LEN_C_O,
+            LEN_N_H,
+            LEN_C_H,
+            LEN_O_H,
+        ];
+
+        let lens_hydrogen = vec![LEN_OH_OH, LEN_NH_OC, LEN_OH_OC];
+        let mut bonds = Vec::new();
+
+        // We use spacial partitioning, so as not to compare every pair of atoms.
+        let grid = self.create_atom_rtree();
+
+        for atom in &grid {
+            for other_atom in
+                grid.locate_within_distance((atom.x(), atom.y(), atom.z()), 2.9_f64.sqrt())
+            {
+                if atom.counter() >= other_atom.counter() {
+                    continue;
+                }
+                let dist = atom.distance(other_atom);
+
+                for (lens, bond_type) in [
+                    (&lens_covalent, Bond::Covalent),
+                    (&lens_hydrogen, Bond::Hydrogen),
+                ] {
+                    for &bond_len in lens {
+                        if (dist - bond_len).abs() < BOND_LEN_THRESH {
+                            bonds.push((atom.counter(), other_atom.counter(), bond_type));
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+        self.bonds.extend(bonds);
+    }
 }
 
 impl fmt::Display for PDB {