datamol.isomers

canonical_tautomer(mol)

Get the canonical tautomer of the current molecule.

Parameters:

Name	Type	Description	Default
mol	Mol	A molecule.	required

Source code in datamol/isomers/_enumerate.py

def canonical_tautomer(mol: dm.Mol):
    """Get the canonical tautomer of the current molecule.

    Args:
        mol: A molecule.
    """
    enumerator = rdMolStandardize.TautomerEnumerator()
    return enumerator.Canonicalize(mol)

count_stereoisomers(mol, n_variants=20, undefined_only=False, rationalise=True, timeout_seconds=None, clean_it=True, precise=False)

Get the number of possible stereoisomers for a molecule.

Warning: By default, this function compute an estimtion number based on the stereo bonds which gives an upper bound of possible stereoisomers. By setting precise=True, the number is computed by enumrerating the stereoisomers. However, it can be computationnaly intensive.

Parameters:

Name	Type	Description	Default
mol	Mol	The molecule whose state we should enumerate.	required
n_variants	int	The maximum amount of molecules that should be returned.	20
undefined_only	bool	If we should enumerate all stereocenters and bonds or only those with undefined stereochemistry.	False
rationalise	bool	If we should try to build and rationalise the molecule to ensure it can exist.	True
timeout_seconds	int	The maximum amount of time to spend on enumeration. None will disable the timeout. Note that the timeout might be inaccurate as a running single variant computation is not stopped when the duration is reached.	None
clean_it	bool	A flag for assigning stereochemistry. If True, it will remove previous stereochemistry markings on the bonds.	True
precise	bool	Whether compute counts by enumerate the stereoisomers using enumerate_stereoisomers.	False

Source code in datamol/isomers/_enumerate.py

def count_stereoisomers(
    mol: dm.Mol,
    n_variants: int = 20,
    undefined_only: bool = False,
    rationalise: bool = True,
    timeout_seconds: int = None,
    clean_it: bool = True,
    precise: bool = False,
):
    """Get the number of possible stereoisomers for a molecule.

    Warning: By default, this function compute an estimtion number based on the stereo bonds which
    gives an upper bound of possible stereoisomers. By setting `precise=True`, the number is computed
    by enumrerating the stereoisomers. However, it can be computationnaly intensive.

    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.
        timeout_seconds: The maximum amount of time to spend on enumeration. None
            will disable the timeout. Note that the timeout might be inaccurate as a running single variant
            computation is not stopped when the duration is reached.
        clean_it: A flag for assigning stereochemistry. If True, it will remove previous stereochemistry
            markings on the bonds.
        precise: Whether compute counts by enumerate the stereoisomers using `enumerate_stereoisomers`.

    """
    if precise:
        num_variants = len(
            enumerate_stereoisomers(
                mol=mol,
                n_variants=n_variants,
                undefined_only=undefined_only,
                rationalise=rationalise,
                timeout_seconds=timeout_seconds,
                clean_it=clean_it,
            )
        )
    else:
        # 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

enumerate_stereoisomers(mol, n_variants=20, undefined_only=False, rationalise=True, timeout_seconds=None, clean_it=True)

Enumerate the stereocenters and bonds of the current molecule.

Original source: the openff-toolkit lib.

Warning: this function can be computationnaly intensive.

Parameters:

Name	Type	Description	Default
mol	Mol	The molecule whose state we should enumerate.	required
n_variants	int	The maximum amount of molecules that should be returned.	20
undefined_only	bool	If we should enumerate all stereocenters and bonds or only those with undefined stereochemistry.	False
rationalise	bool	If we should try to build and rationalise the molecule to ensure it can exist.	True
timeout_seconds	int	The maximum amount of time to spend on enumeration. None will disable the timeout. Note that the timeout might be inaccurate as a running single variant computation is not stopped when the duration is reached.	None
clean_it	bool	A flag for assigning stereochemistry. If True, it will remove previous stereochemistry markings on the bonds.	True

Source code in datamol/isomers/_enumerate.py

def enumerate_stereoisomers(
    mol: dm.Mol,
    n_variants: int = 20,
    undefined_only: bool = False,
    rationalise: bool = True,
    timeout_seconds: int = None,
    clean_it: bool = True,
):
    """Enumerate the stereocenters and bonds of the current molecule.

    Original source: the `openff-toolkit` lib.

    Warning: this function can be computationnaly intensive.

    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.
        timeout_seconds: The maximum amount of time to spend on enumeration. None
            will disable the timeout. Note that the timeout might be inaccurate as a running single variant
            computation is not stopped when the duration is reached.
        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,
        maxIsomers=n_variants,
        unique=True,
    )

    isomers_iterator = EnumerateStereoisomers(mol, options=stereo_opts)

    start = time.time()
    duration = 0

    isomers = []
    while timeout_seconds is None or duration < timeout_seconds:
        try:
            isomer = next(isomers_iterator)
            isomers.append(isomer)
        except StopIteration:
            break

        duration = time.time() - start

    variants = []
    for isomer in isomers:
        # isomer has CIS/TRANS tags so convert back to E/Z
        Chem.SetDoubleBondNeighborDirections(isomer)  # type: ignore
        Chem.AssignStereochemistry(isomer, force=True, cleanIt=clean_it)  # type: ignore
        variants.append(isomer)

    return variants

enumerate_structisomers(mol, n_variants=20, allow_cycle=False, allow_double_bond=False, allow_triple_bond=False, depth=None, timeout_seconds=None)

Enumerate the structural isomers of the input molecule

Warning: this function can be computationnaly intensive.

Parameters:

Name	Type	Description	Default
mol	Mol	The molecule whose state we should enumerate.	required
n_variants	int	The maximum amount of molecules that should be returned.	20
allow_cycle	bool	whether to allow transformation involving cycle creation.	False
allow_double_bond	bool	whether to allow transformation involving at least one double bond.	False
allow_triple_bond	bool	whether to allow transformation involving at least one triple bond.	False
depth	int	depth of the search, use a sensible value to decrease search time. Will fall back to number of atoms.	None
timeout_seconds	int	The maximum amount of time to spend on enumeration. None will disable the timeout. Note that the timeout might be inaccurate as a running single variant computation is not stopped when the duration is reached.	None

Source code in datamol/isomers/_enumerate.py

def enumerate_structisomers(
    mol: dm.Mol,
    n_variants: int = 20,
    allow_cycle: bool = False,
    allow_double_bond: bool = False,
    allow_triple_bond: bool = False,
    depth: int = None,
    timeout_seconds: int = None,
):
    """Enumerate the structural isomers of the input molecule

    Warning: this function can be computationnaly intensive.

    Args:
        mol: The molecule whose state we should enumerate.
        n_variants: The maximum amount of molecules that should be returned.
        allow_cycle: whether to allow transformation involving cycle creation.
        allow_double_bond: whether to allow transformation involving at least one double bond.
        allow_triple_bond: whether to allow transformation involving at least one triple bond.
        depth: depth of the search, use a sensible value to decrease search time. Will fall back to number of atoms.
        timeout_seconds: The maximum amount of time to spend on enumeration. None
            will disable the timeout. Note that the timeout might be inaccurate as a running single variant
            computation is not stopped when the duration is reached.
    """

    mol = dm.copy_mol(mol)

    enumerator = IsomerEnumerator(
        allow_cycle=allow_cycle,
        allow_double_bond=allow_double_bond,
        allow_triple_bond=allow_triple_bond,
    )

    if depth is None:
        depth = mol.GetNumAtoms()

    isomers_iterator = enumerator(
        mol,
        max_mols=n_variants,
        include_input=False,
        depth=depth,
    )

    start = time.time()
    duration = 0

    isomers = []
    while timeout_seconds is None or duration < timeout_seconds:
        try:
            isomer = next(isomers_iterator)
            isomers.append(dm.to_mol(isomer))
        except StopIteration:
            break

        duration = time.time() - start

    return isomers

enumerate_tautomers(mol, n_variants=20, max_transforms=1000, reassign_stereo=True, remove_bond_stereo=True, remove_sp3_stereo=True)

Enumerate the possible tautomers of the current molecule.

Parameters:

Name	Type	Description	Default
mol	Mol	The molecule whose state we should enumerate.	required
n_variants	int	The maximum amount of molecules that should be returned.	20
max_transforms	int	Set the maximum number of transformations to be applied. This limit is usually hit earlier than the n_variants limit and leads to a more linear scaling of CPU time with increasing number of tautomeric centers (see Sitzmann et al.).	1000
reassign_stereo	bool	Whether to reassign stereo centers.	True
remove_bond_stereo	bool	Whether stereochemistry information is removed from double bonds involved in tautomerism. This means that enols will lose their E/Z stereochemistry after going through tautomer enumeration because of the keto-enolic tautomerism.	True
remove_sp3_stereo	bool	Whether stereochemistry information is removed from sp3 atoms involved in tautomerism. This means that S-aminoacids will lose their stereochemistry after going through tautomer enumeration because of the amido-imidol tautomerism.	True

Source code in datamol/isomers/_enumerate.py

def enumerate_tautomers(
    mol: dm.Mol,
    n_variants: int = 20,
    max_transforms: int = 1000,
    reassign_stereo: bool = True,
    remove_bond_stereo: bool = True,
    remove_sp3_stereo: bool = True,
):
    """Enumerate the possible tautomers of the current molecule.

    Args:
        mol: The molecule whose state we should enumerate.
        n_variants: The maximum amount of molecules that should be returned.
        max_transforms: Set the maximum number of transformations to be applied. This limit is usually
            hit earlier than the n_variants limit and leads to a more linear scaling
            of CPU time with increasing number of tautomeric centers (see Sitzmann et al.).
        reassign_stereo: Whether to reassign stereo centers.
        remove_bond_stereo: Whether stereochemistry information is removed from double bonds involved in tautomerism.
            This means that enols will lose their E/Z stereochemistry after going through tautomer enumeration because
            of the keto-enolic tautomerism.
        remove_sp3_stereo: Whether stereochemistry information is removed from sp3 atoms involved in tautomerism. This
            means that S-aminoacids will lose their stereochemistry after going through tautomer enumeration because
            of the amido-imidol tautomerism.
    """
    enumerator = rdMolStandardize.TautomerEnumerator()

    # Configure the enumeration
    enumerator.SetMaxTautomers(n_variants)
    enumerator.SetMaxTransforms(max_transforms)
    enumerator.SetReassignStereo(reassign_stereo)
    enumerator.SetRemoveBondStereo(remove_bond_stereo)
    enumerator.SetRemoveSp3Stereo(remove_sp3_stereo)

    tautomers = enumerator.Enumerate(mol)
    return list(tautomers)

remove_stereochemistry(mol, copy=True)

Removes all stereochemistry info from the molecule.

Parameters:

Name	Type	Description	Default
mol	Mol	a molecule	required
copy	bool	Whether to copy the molecule.	True

Source code in datamol/isomers/_enumerate.py

def remove_stereochemistry(mol: dm.Mol, copy: bool = True):
    """Removes all stereochemistry info from the molecule.

    Args:
        mol: a molecule
        copy: Whether to copy the molecule.
    """

    if copy is True:
        mol = dm.copy_mol(mol)
    rdmolops.RemoveStereochemistry(mol)
    return mol

IsomerEnumerator

Implementation of the isomer enumeration algorithm described in https://doi.org/10.1186/s13321-017-0252-9

..note:: Due to the chemical reaction used, atom mapping will not be preserved !

Source code in datamol/isomers/_structural.py

class IsomerEnumerator:
    """
    Implementation of the isomer enumeration algorithm described in  https://doi.org/10.1186/s13321-017-0252-9

    ..note::
        Due to the chemical reaction used, atom mapping will not be preserved !

    """

    def __init__(
        self,
        allow_cycle: bool = True,
        allow_double_bond: bool = False,
        allow_triple_bond: bool = False,
        rxn_list: Optional[list] = None,
    ):
        """Structural isomer enumeration

        Args:
            allow_cycle (bool, optional): whether to allow transformation involving cycle creation
            allow_double_bond (bool, optional): whether to allow transformation involving at least one double bond
            allow_triple_bond (bool, optional): whether to allow transformation involving at least one triple bond
            rxn_list (list, optional): List of str to specifiy the reactions that should be used
        """
        self.allow_cycle = allow_cycle
        self.allow_triple_bond = allow_triple_bond
        self.allow_double_bond = allow_double_bond
        self.rxn_cache = self._resolve(rxn_list)
        self._ring_smarts = dm.from_smarts("[R]")
        # there should be a lot more of this
        self._impossible_chemistry = [dm.from_smarts("[$([*;r3,r4](=*)=*)]")]
        if not self.rxn_cache:
            raise ValueError("No isomeric transformation matches your")

    def _resolve(self, rxn_list: Optional[list]):
        """Resolve map of rxn to use

        Args:
            rxn_list: list of reaction to use
        """
        if rxn_list is not None and len(rxn_list) > 0:
            rxn_list = [_rxn for _rxn in rxn_list if _rxn in ISOMERIZERS]
        else:
            rxn_list = list([x for x, val in ISOMERIZERS.items() if val.use])

        final_rxns = dict()

        for x in rxn_list:
            rxn_trans = ISOMERIZERS[x]
            # check for triple bonds
            can_add = (not rxn_trans.triplebond) or self.allow_triple_bond
            # check for double bonds
            can_add &= not (rxn_trans.doublebond) or self.allow_double_bond
            # check for cycles
            can_add &= rxn_trans.acyclic or self.allow_cycle
            if can_add:
                as_rxn = Chem.rdChemReactions.ReactionFromSmarts(rxn_trans.smarts)
                final_rxns[x] = as_rxn
        return final_rxns

    def _clean(self, mol: dm.Mol):
        """Clean and sanitize a molecule during enumeration

        Args:
            mol: molecule to clean
        """
        clean_mol = None
        clean_smiles = ""
        try:
            mol = dm.sanitize_mol(mol)
            Chem.rdmolops.AssignRadicals(mol)
            for a in mol.GetAtoms():
                a.SetNoImplicit(True)
                a.SetNumExplicitHs(a.GetTotalNumHs() + a.GetNumRadicalElectrons())
                a.SetNumRadicalElectrons(0)
            mol.UpdatePropertyCache()
            clean_smiles = dm.to_smiles(mol)
            clean_mol = dm.to_mol(clean_smiles)
        except Exception:
            pass

        return clean_mol, clean_smiles

    def _is_valid(self, new_mol, new_smiles, need_substruct: Optional[dm.Mol] = None):
        """Check whether the mol or smiles is valid"""
        # not disconnected
        is_ok = new_smiles.count(".") == 0
        if not self.allow_cycle:
            is_ok &= not new_mol.HasSubstructMatch(self._ring_smarts)
        is_ok &= not (any(new_mol.HasSubstructMatch(x) for x in self._impossible_chemistry))
        if need_substruct is not None:
            is_ok &= new_mol.HasSubstructMatch(need_substruct)
        return is_ok

    def _has_correct_size(self, new_mol, expected_size):
        """Check whether the molecule has the correct size"""
        return new_mol.GetNumHeavyAtoms() == expected_size

    def __call__(self, *args, **kwargs):
        return self.enumerate(*args, **kwargs)

    def enumerate(
        self,
        mol: dm.Mol,
        depth: Optional[int] = None,
        include_input: bool = True,
        protect_substruct: Optional[Chem.rdchem.Mol] = None,
        max_mols: Optional[int] = None,
    ):
        """Enumerate the list of isomers of the current mol

        Args:
            mol: input molecule or smiles
            depth: optional maximum depth of the breadth search (default=None)
            include_input: whether to include the input molecule (default=True)
            protect_substruct: Optional substruct to protect
            max_mols: maximum number of molecule to sample
        """

        mol = dm.to_mol(mol)
        if protect_substruct is not None:
            mol = dm.protect_atoms(mol, substruct=protect_substruct, in_place=False)
        smiles = dm.to_smiles(mol, isomeric=True, canonical=True)
        mol_size = mol.GetNumHeavyAtoms()
        source_mols = [(mol, 0)]
        seen = set([smiles])
        seen_inchi_key = set([dm.to_inchikey(smiles)])

        _depth = depth or float("inf")
        _max_mols = max_mols or float("inf")

        if include_input:
            yield smiles
        while len(source_mols) > 0 and _max_mols > 0:
            curmol, curdepth = source_mols.pop(0)
            try:
                # we first try to kekulize the molecule
                Chem.Kekulize(curmol)
                curmol = Chem.AddHs(curmol, explicitOnly=False)
            except Exception:
                pass

            if curdepth >= _depth:
                return
            with dm.without_rdkit_log():
                for k, rxn in self.rxn_cache.items():
                    if _max_mols < 1:
                        break
                    for new_mols in rxn.RunReactants((curmol,)):
                        if _max_mols < 1:
                            break
                        new_mol = new_mols[0]
                        # sanitize mol
                        new_mol, new_smiles = self._clean(new_mol)
                        # need to clean twice, not sure why
                        _, new_smiles = self._clean(new_mol)
                        new_smiles_inchikey = dm.to_inchikey(new_mol)

                        if (
                            new_mol is None
                            or new_smiles_inchikey in seen_inchi_key
                            or not self._has_correct_size(new_mol, mol_size)
                        ):
                            continue
                        if protect_substruct is not None:
                            new_mol = dm.protect_atoms(
                                new_mol,
                                substruct=protect_substruct,
                                in_place=False,
                            )
                        seen.add(new_smiles)
                        seen_inchi_key.add(new_smiles_inchikey)
                        source_mols.append((new_mol, curdepth + 1))
                        if self._is_valid(new_mol, new_smiles, protect_substruct):
                            yield new_smiles
                            _max_mols -= 1

__init__(allow_cycle=True, allow_double_bond=False, allow_triple_bond=False, rxn_list=None)

Structural isomer enumeration

Parameters:

Name	Type	Description	Default
allow_cycle	bool	whether to allow transformation involving cycle creation	True
allow_double_bond	bool	whether to allow transformation involving at least one double bond	False
allow_triple_bond	bool	whether to allow transformation involving at least one triple bond	False
rxn_list	list	List of str to specifiy the reactions that should be used	None

Source code in datamol/isomers/_structural.py

def __init__(
    self,
    allow_cycle: bool = True,
    allow_double_bond: bool = False,
    allow_triple_bond: bool = False,
    rxn_list: Optional[list] = None,
):
    """Structural isomer enumeration

    Args:
        allow_cycle (bool, optional): whether to allow transformation involving cycle creation
        allow_double_bond (bool, optional): whether to allow transformation involving at least one double bond
        allow_triple_bond (bool, optional): whether to allow transformation involving at least one triple bond
        rxn_list (list, optional): List of str to specifiy the reactions that should be used
    """
    self.allow_cycle = allow_cycle
    self.allow_triple_bond = allow_triple_bond
    self.allow_double_bond = allow_double_bond
    self.rxn_cache = self._resolve(rxn_list)
    self._ring_smarts = dm.from_smarts("[R]")
    # there should be a lot more of this
    self._impossible_chemistry = [dm.from_smarts("[$([*;r3,r4](=*)=*)]")]
    if not self.rxn_cache:
        raise ValueError("No isomeric transformation matches your")

enumerate(mol, depth=None, include_input=True, protect_substruct=None, max_mols=None)

Enumerate the list of isomers of the current mol

Parameters:

Name	Type	Description	Default
mol	Mol	input molecule or smiles	required
depth	Optional[int]	optional maximum depth of the breadth search (default=None)	None
include_input	bool	whether to include the input molecule (default=True)	True
protect_substruct	Optional[Mol]	Optional substruct to protect	None
max_mols	Optional[int]	maximum number of molecule to sample	None

Source code in datamol/isomers/_structural.py

def enumerate(
    self,
    mol: dm.Mol,
    depth: Optional[int] = None,
    include_input: bool = True,
    protect_substruct: Optional[Chem.rdchem.Mol] = None,
    max_mols: Optional[int] = None,
):
    """Enumerate the list of isomers of the current mol

    Args:
        mol: input molecule or smiles
        depth: optional maximum depth of the breadth search (default=None)
        include_input: whether to include the input molecule (default=True)
        protect_substruct: Optional substruct to protect
        max_mols: maximum number of molecule to sample
    """

    mol = dm.to_mol(mol)
    if protect_substruct is not None:
        mol = dm.protect_atoms(mol, substruct=protect_substruct, in_place=False)
    smiles = dm.to_smiles(mol, isomeric=True, canonical=True)
    mol_size = mol.GetNumHeavyAtoms()
    source_mols = [(mol, 0)]
    seen = set([smiles])
    seen_inchi_key = set([dm.to_inchikey(smiles)])

    _depth = depth or float("inf")
    _max_mols = max_mols or float("inf")

    if include_input:
        yield smiles
    while len(source_mols) > 0 and _max_mols > 0:
        curmol, curdepth = source_mols.pop(0)
        try:
            # we first try to kekulize the molecule
            Chem.Kekulize(curmol)
            curmol = Chem.AddHs(curmol, explicitOnly=False)
        except Exception:
            pass

        if curdepth >= _depth:
            return
        with dm.without_rdkit_log():
            for k, rxn in self.rxn_cache.items():
                if _max_mols < 1:
                    break
                for new_mols in rxn.RunReactants((curmol,)):
                    if _max_mols < 1:
                        break
                    new_mol = new_mols[0]
                    # sanitize mol
                    new_mol, new_smiles = self._clean(new_mol)
                    # need to clean twice, not sure why
                    _, new_smiles = self._clean(new_mol)
                    new_smiles_inchikey = dm.to_inchikey(new_mol)

                    if (
                        new_mol is None
                        or new_smiles_inchikey in seen_inchi_key
                        or not self._has_correct_size(new_mol, mol_size)
                    ):
                        continue
                    if protect_substruct is not None:
                        new_mol = dm.protect_atoms(
                            new_mol,
                            substruct=protect_substruct,
                            in_place=False,
                        )
                    seen.add(new_smiles)
                    seen_inchi_key.add(new_smiles_inchikey)
                    source_mols.append((new_mol, curdepth + 1))
                    if self._is_valid(new_mol, new_smiles, protect_substruct):
                        yield new_smiles
                        _max_mols -= 1