Added function to get the number of stereoisomers

datamol/__init__.py

@@ -135,6 +135,7 @@
     "save_df": "datamol.io",
     # isomers
     "enumerate_stereoisomers": "datamol.isomers",
+    "count_stereoisomers": "datamol.isomers",
     "enumerate_tautomers": "datamol.isomers",
     "enumerate_structisomers": "datamol.isomers",
     "canonical_tautomer": "datamol.isomers",
@@ -329,6 +330,7 @@ def __dir__():
     from .io import open_df
 
     from .isomers import enumerate_stereoisomers
+    from .isomers import count_stereoisomers
     from .isomers import enumerate_tautomers
     from .isomers import enumerate_structisomers
     from .isomers import canonical_tautomer

datamol/isomers/__init__.py

@@ -3,5 +3,6 @@
 from ._enumerate import enumerate_stereoisomers
 from ._enumerate import enumerate_tautomers
 from ._enumerate import enumerate_structisomers
+from ._enumerate import count_stereoisomers
 from ._enumerate import remove_stereochemistry
 from ._enumerate import canonical_tautomer

datamol/isomers/_enumerate.py

@@ -5,6 +5,7 @@
 
 from rdkit.Chem.MolStandardize import rdMolStandardize
 from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers
+from rdkit.Chem.EnumerateStereoisomers import GetStereoisomerCount
 from rdkit.Chem.EnumerateStereoisomers import StereoEnumerationOptions
 
 import datamol as dm
@@ -82,7 +83,8 @@ def enumerate_stereoisomers(
 
     # in case any bonds/centers are missing stereo chem flag it here
     Chem.AssignStereochemistry(mol, force=False, flagPossibleStereoCenters=True, cleanIt=clean_it)  # type: ignore
-    Chem.FindPotentialStereoBonds(mol, cleanIt=clean_it)  # type: ignore
+    # lu: do not clean (overwrite bond stereo information) when set `undefined_only=Ture`
+    Chem.FindPotentialStereoBonds(mol, cleanIt=not undefined_only and clean_it)
 
     # set up the options
     stereo_opts = StereoEnumerationOptions(
@@ -117,6 +119,43 @@ def enumerate_stereoisomers(
     return variants
 
 
+def count_stereoisomers(
+    mol: dm.Mol, undefined_only: bool = False, rationalise: bool = True, clean_it: bool = True
+):
+    """Get an estimate (upper bound) of the number of possible stereoisomers for a molecule.
+
+    Warning: this function is an estimtion, therefore it might be less accurate than enumrerating the stereoisomers.
+
+    Args:
+        mol: The molecule whose state we should enumerate.
+        n_variants: The maximum amount of molecules that should be returned.
+        undefined_only: If we should enumerate all stereocenters and bonds or only those
+            with undefined stereochemistry.
+        rationalise: If we should try to build and rationalise the molecule to ensure it
+            can exist.
+        clean_it: A flag for assigning stereochemistry. If True, it will remove previous stereochemistry
+            markings on the bonds.
+    """
+    # safety first
+    mol = dm.copy_mol(mol)
+
+    # in case any bonds/centers are missing stereo chem flag it here
+    Chem.AssignStereochemistry(mol, force=False, flagPossibleStereoCenters=True, cleanIt=clean_it)  # type: ignore
+    # lu: do not clean (overwrite bond stereo information) when set `undefined_only=Ture`
+    Chem.FindPotentialStereoBonds(mol, cleanIt=not undefined_only and clean_it)
+
+    # set up the options
+    stereo_opts = StereoEnumerationOptions(
+        tryEmbedding=rationalise,
+        onlyUnassigned=undefined_only,
+        unique=True,
+    )
+
+    num_variants = GetStereoisomerCount(mol, options=stereo_opts)
+
+    return num_variants
+
+
 def enumerate_structisomers(
     mol: dm.Mol,
     n_variants: int = 20,

tests/test_isomers.py

@@ -43,6 +43,14 @@ def test_enumerate_stereo_timeout():
     dm.enumerate_stereoisomers(mol, n_variants=2, timeout_seconds=1)
 
 
+def test_count_stereoisomers():
+    num_isomers_1 = dm.count_stereoisomers(dm.to_mol("CC=CC"), undefined_only=True)
+    num_isomers_2 = dm.count_stereoisomers(dm.to_mol("CC=CC"), undefined_only=False)
+    assert num_isomers_1 == num_isomers_2
+
+    assert dm.count_stereoisomers(dm.to_mol("Br/C=C\\Br"), undefined_only=True) == 1
+
+
 def test_enumerate_structural():
     mol = dm.to_mol("CCCCC")  # pentane has only three structural isomers