bug fix in mcs_rmsd

useful_rdkit_utils/__init__.py

@@ -15,4 +15,4 @@
 from .split_utils import *
 from .scaffold_finder import *
 
-__version__ = "0.63"
+__version__ = "0.64"

useful_rdkit_utils/geometry.py

@@ -8,6 +8,7 @@
 from rdkit.Chem.rdFMCS import FindMCS
 from rdkit.Chem.rdMolTransforms import ComputeCentroid
 from rdkit.Chem.rdchem import Mol
+from typing import Tuple
 
 
 # ----------- Molecular geometry
@@ -62,27 +63,32 @@ def gen_conformers(mol, num_confs=10):
     return mol
 
 
-def mcs_rmsd(mol_1: Mol, mol_2: Mol) -> float:
+def mcs_rmsd(mol_1: Mol, mol_2: Mol, id_1: int = 0, id_2: int = 0) -> Tuple[int, float]:
     """
-    Calculate the minimum Root Mean Square Deviation (RMSD) between the Maximum Common Substructure (MCS)
-    of two RDKit molecule objects.
+    Calculate the RMSD (Root Mean Square Deviation) between the MCS (Maximum Common Substructure) of two molecules.
 
-    :param mol_1: The first RDKit molecule object.
-    :param mol_2: The second RDKit molecule object.
-    :return: The minimum RMSD as a float.
+    :param mol_1: First RDKit molecule
+    :param mol_2: Second RDKit molecule
+    :param id_1: Conformer ID for the first molecule
+    :param id_2: Conformer ID for the second molecule
+    :return: A tuple containing the number of MCS atoms and the RMSD value
     """
     mcs_res = FindMCS([mol_1, mol_2])
+    num_mcs_atoms = mcs_res.numAtoms
     pat = Chem.MolFromSmarts(mcs_res.smartsString)
     match_1 = mol_1.GetSubstructMatches(pat)
     match_2 = mol_2.GetSubstructMatches(pat)
     min_rmsd = 1e6
     for m1 in match_1:
         for m2 in match_2:
-            crd_1 = mol_1.GetConformer(0).GetPositions()[list(m1)]
-            crd_2 = mol_2.GetConformer(0).GetPositions()[list(m2)]
-            rmsd = np.sqrt((((crd_1 - crd_2) ** 2).sum()).mean())
+            crd_1 = mol_1.GetConformer(id_1).GetPositions()[list(m1)]
+            crd_2 = mol_2.GetConformer(id_2).GetPositions()[list(m2)]
+            diff = crd_1 - crd_2
+            squared_dist = np.sum(diff ** 2, axis=1)
+            msd = np.mean(squared_dist)
+            rmsd = np.sqrt(msd)
             min_rmsd = min(min_rmsd, rmsd)
-    return min_rmsd
+    return num_mcs_atoms, float(min_rmsd)
 
 
 # from https://birdlet.github.io/2019/10/02/py3dmol_example/