RMSD (with sub-atom selection) (resolves #28)

atomium/structures.py

@@ -164,6 +164,27 @@ def pairing_with(self, structure, *args, **kwargs):
                 a.element, a.name, a.id, a.distance_to(center), id(a)
             ))
         return {atom1: atom2 for atom1, atom2 in zip(self_atoms, other_atoms)}
+
+
+    def rmsd_with(self, structure, *args, **kwargs):
+        pairing = self.pairing_with(structure, *args, **kwargs)
+        coords1, coords2 = [[
+            a.location for a in atoms
+        ] for atoms in zip(*pairing.items())]
+        coords1 = coords1 - np.array(self.center_of_mass)
+        coords2 = coords2 - np.array(structure.center_of_mass)
+        cov_matrix = np.dot(np.transpose(coords1), coords2)
+        V, S, W = np.linalg.svd(cov_matrix)
+        d = (np.linalg.det(V) * np.linalg.det(W)) < 0.0
+        if d: S[-1] = -S[-1]
+        if d: V[:, -1] = -V[:, -1]
+        rot_matrix = np.dot(V, W)
+        coords1_rotated = np.dot(coords1, rot_matrix)
+        diff = coords1_rotated - coords2
+        N = len(coords2)
+        return np.sqrt((diff * diff).sum() / N)
+
+
 
 
 class Model(AtomStructure, metaclass=StructureClass):

requirements.txt

@@ -1,6 +1,5 @@
 numpy
 requests
-rmsd
 msgpack==0.6.1
 valerius==0.2
 python-coveralls

tests/integration/test_interactions.py

@@ -159,16 +159,17 @@ def test_can_get_xmp_pairing(self):
         for a1, a2 in pairing.items():
             self.assertEqual(a1.element, a2.element)
             self.assertEqual(a1.name, a2.name)
+        self.assertEqual(round(mol1.rmsd_with(mol2), 3), 0.48)
 
 
-
     def test_can_get_bu2_pairing(self):
         mol1 = self.pdb.model.non_polymer("D")
         mol2 = self.pdb.model.non_polymer("F")
         pairing = mol1.pairing_with(mol2)
         for a1, a2 in pairing.items():
             self.assertEqual(a1.element, a2.element)
             self.assertEqual(a1.name, a2.name)
+        self.assertEqual(round(mol1.rmsd_with(mol2), 3), 0.133)
 
 
     def test_can_get_polymer_pairing(self):
@@ -178,9 +179,7 @@ def test_can_get_polymer_pairing(self):
         sorted_by_first_id = sorted(pairing.items(), key=lambda p: p[0].id)
         second_ids = [a[1].id for a in sorted_by_first_id]
         self.assertEqual(second_ids, sorted(second_ids))
-
-
-
+        self.assertEqual(round(mol1.rmsd_with(mol2, name="CA", residue__number__lt=187), 3), 1.133)
 
 
     def test_can_handle_mispairing(self):