Re: [Rdkit-discuss] detect dihedral angles in a conformation

[Paolo T. <pao...@un...>]

Dear Jose Manuel,

the enumerateTorsions() function in this Python script should allow to 
enumerate all torsions you might want to constrain in a molecule; it'll 
be up to you to carry out further filtering, though:

#!/usr/bin/env python


import sys
import rdkit
from rdkit import Chem


def enumerateTorsions(mol):
   torsionSmarts = '[!$(*#*)&!D1]~[!$(*#*)&!D1]'
   torsionQuery = Chem.MolFromSmarts(torsionSmarts)
   matches = mol.GetSubstructMatches(torsionQuery)
   torsionList = []
   for match in matches:
     idx2 = match[0]
     idx3 = match[1]
     bond = mol.GetBondBetweenAtoms(idx2, idx3)
     jAtom = mol.GetAtomWithIdx(idx2)
     kAtom = mol.GetAtomWithIdx(idx3)
     if (((jAtom.GetHybridization() != Chem.HybridizationType.SP2)
       and (jAtom.GetHybridization() != Chem.HybridizationType.SP3))
       or ((kAtom.GetHybridization() != Chem.HybridizationType.SP2)
       and (kAtom.GetHybridization() != Chem.HybridizationType.SP3))):
       continue
     for b1 in jAtom.GetBonds():
       if (b1.GetIdx() == bond.GetIdx()):
         continue
       idx1 = b1.GetOtherAtomIdx(idx2)
       for b2 in kAtom.GetBonds():
         if ((b2.GetIdx() == bond.GetIdx())
           or (b2.GetIdx() == b1.GetIdx())):
           continue
         idx4 = b2.GetOtherAtomIdx(idx3)
         # skip 3-membered rings
         if (idx4 == idx1):
           continue
         torsionList.append((idx1, idx2, idx3, idx4))
   return torsionList

if (__name__ == "__main__"):
   mol = Chem.MolFromSmiles('CCCCCC')
   mol = Chem.AddHs(mol)
   torsionList = enumerateTorsions(mol)
   for torsion in torsionList:
     sys.stdout.write(str(torsion) + '\n')

Cheers,
p.

On 20/10/2015 09:41, Jose Manuel Gally wrote:
> Hi RDKitters,
>
> I would like to consider parts of a conformation rigid (fixed dihedral
> angles) during minimization
> My end goal would be to generate only ring conformations starting with
> valid 3D molecules.
>
> I can already consider a specific dihedral angle as rigid:
>
> from rdkit import Chem
> from rdkit.Chem import AllChem, rdMolTransforms
>
> # create test mol
> s = 'COCCN1CCOCC1'
> m = Chem.MolFromSmiles(s)
> m = Chem.AddHs(m)
>
> # add 3D coordinates
> AllChem.EmbedMolecule(m)
>
> # freeze one dihedral angle (composed of atoms 0-3)
> MMFFs_MP = AllChem.MMFFGetMoleculeProperties(m, mmffVariant='MMFF94s')
> MMFFs_FF = AllChem.MMFFGetMoleculeForceField(m, MMFFs_MP)
> MMFFs_FF.MMFFAddTorsionConstraint(0, 1, 2, 3, relative=True,
> minDihedralDeg=0.0, maxDihedralDeg=0.0,forceConstant=99999999999999.0)
> c = m.GetConformer()
> print "before min", rdMolTransforms.GetDihedralDeg(c, 0,1,2,3) #
> -53.0873064656
>
> # minimize molecule with constrained dihedral angle
> MMFFs_FF.Minimize(maxIts=10)
> print "after first min", rdMolTransforms.GetDihedralDeg(c,0,1,2,3) #
> -53.0873064656
> MMFFs_FF.Minimize(maxIts=10)
> print "after second min", rdMolTransforms.GetDihedralDeg(c,0,1,2,3) #
> -53.0873064656
>
> However, I have difficulties to find all dihedral angles to consider
> rigid...
> I would like to detect dihedral angles with 4 atoms where:
>       - none is hydrogen
>       - no more than 2 atoms are in different rings
>
> First I looked for a function to return me the list of dihedral angles
> and iterate over it, but could not find any.
> My other alternative would be to iterate over atoms to get their
> neighbors, and then get their neighbor' neighbors, but that looks very
> very slow.
> Any other way to do this?
>
> Thank you!
>
> Jose Manuel
>
> ------------------------------------------------------------------------------
> _______________________________________________
> Rdkit-discuss mailing list
> Rdk...@li...
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss