On Aug 31, 2018, at 07:41, Paolo Tosco <[email protected]> wrote:
> this gist should do what you need:
Unless I misinterpreted what Jim is looking for, I don't think that returns the
contiguous rotatable bonds in a small molecule.
In the following there are only two rotatable bonds:
>>> mol = Chem.MolFromSmiles("NCC1CCC1C(=O)O")
>>> mol.GetSubstructMatches(RotatableBondSmarts)
((1, 2), (5, 6))
The code from the gist identifies three bonds:
>>> find_bond_groups(mol)
((<rdkit.Chem.rdchem.Bond object at 0x108604a80>, <rdkit.Chem.rdchem.Bond
object at 0x108604bc0>, <rdkit.Chem.rdchem.Bond object at 0x1086049e0>),)
>>> [(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in _[0]]
[(5, 2), (5, 6), (1, 2)]
In this case, (2,5) is part of a ring, and not rotatable.
I think the problem is that:
nbrs = [nbr.GetIdx() for nbr in a.GetNeighbors() if (
(nbr.GetIdx() in rot_atom_set_tmp) and (not (nbr.GetIdx() in
connected_atom_set)))]
finds that both atoms of the bond are in connected bond groups, while the bond
itself is not part of the match.
I have put an alternative implementation of this function at the bottom of this
email. For my test case it returns:
>>> find_bond_groups2(mol)
((<rdkit.Chem.rdchem.Bond object at 0x1086048a0>,), (<rdkit.Chem.rdchem.Bond
object at 0x108604800>,))
>>> [[(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in x] for x in _]
[[(5, 6)], [(1, 2)]]
Cheers,
Andrew
[email protected]
from collections import defaultdict
def find_bond_groups_dalke(mol):
rot_atom_pairs = mol.GetSubstructMatches(RotatableBondSmarts)
bond_groups = defaultdict(set)
for (left, right) in rot_atom_pairs:
# Ensure they are in order (I'm not sure if this is required)
if right < left:
left, right = right, left
pair = (left, right)
# Add the pair to the group associated with the left atom
bond_groups[left].add(pair)
# Merge the left atom group with the right atom group
bond_groups[left].update(bond_groups[right])
bond_groups[right] = bond_groups[left]
# Get just the unique bond groups, in order from largest to smallest
unique_bond_groups = set(frozenset(v) for v in bond_groups.values())
sorted_bond_groups = sorted(unique_bond_groups, reverse=True, key=len)
# Return as bond objects, to match Paolo's code
return tuple(
tuple(mol.GetBondBetweenAtoms(left, right) for (left, right) in
bond_group)
for bond_group in sorted_bond_groups)
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
