[Rdkit-discuss] detect dihedral angles in a conformation

2015-10-20 Thread Jose Manuel Gally 
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=99.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
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss

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

2015-10-20 Thread Michal Krompiec 
Hi Jose
I have a similar problem, but I look for dihedrals between
aromatic/unsaturated rings:

rotatable_ring_bonds=Chem.MolFromSmarts("[!$(*#*)]-!@[!$(*#*)]")
rotatable_matches=mol.GetSubstructMatches(rotatable_ring_bonds)
rotatable=set() #set of sorted indices of the rotatable bonds
for bond in rotatable_matches:
if bond[0] 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=99.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
> Rdkit-discuss@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
>
--
___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss

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

2015-10-20 Thread Jose Manuel Gally 
Hi Paolo, Michal,

thanks for your replies!

I don't want to freeze only linkers between saturated rings, but also 
'side chains' so I went with Paolo's snippet to get the maximum amount 
of torsions.
Then I filtered the torsions having either:
-  at least 2 atoms in the same ring
-  at least 1 hydrogen

This works for linkers and "side chains"!
Would it make sense to also freeze torsions in aromatic rings to speed 
up minimization, given that they should not move that much anyway?

Here is how I (sub-optimally) did, using Paolo's function:

import sys
from rdkit import Chem
from rdkit.Chem import AllChem, rdMolTransforms

# create test mol
s = 'O1=CC(NCC2=CC=CC=C2)=CC2=C1C=C(OCC1=CC=CC=C1)C=C2'
m = Chem.MolFromSmiles(s)
m = Chem.AddHs(m)
# add 3D coordinates
AllChem.EmbedMolecule(m)

# get full torsion list
torsionList = enumerateTorsions(m)
# detect Hs
p = '[#1&!#7]'# match all hs
query = Chem.MolFromSmarts(p)
temp_Hs = m.GetSubstructMatches(query)
Hs=[]
for H in temp_Hs:
 Hs.append(H[0])

# get rings info
RI=m.GetRingInfo()
rings = RI.AtomRings()

# Filter full torsion list
torsionList_noH=[]
for torsion in torsionList:
 ok=True
 print "\ncurrent torsion:", torsion
 print "read as : ", [x +1 for x in torsion]
 # check if 2 atoms in same ring
 for ring in rings:
 print "ring:", ring
 notInRing = [ x for x in torsion if x not in ring ]
 if len(notInRing) < 3:
 ok=False
 print "ring check : rejected (2 atoms in same ring)"
 else:
 print "ring check : passed"
 if ok:
 noHs = [ x for x in torsion if x not in Hs ]
 print "noHs", noHs
 if len(noHs)==4:
 print "hydrogen : passed"
 torsionList_noH.append(noHs)
 else:
 print "hydrogen : rejected (hydrogen found)"
 for t in torsion:
 print m.GetAtomWithIdx(t).GetSymbol()
 else:
 print "rejected"

# remove duplicates
torsionList_noH  = [list(i) for i in set(tuple(i) for i in 
torsionList_noH)]# (2)

### RESULTS
print "\nFINAL: ANGLES TO CONSTRAINT:"
for t in torsionList_noH:
 print [x +1 for x in t]

# freeze one dihedral angle
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=9.0)
c = m.GetConformer()
print "before min"
# check values
for t in torsionList_noH:
 print t, rdMolTransforms.GetDihedralDeg(c, t[0], t[1], t[2], t[3])
# minimize molecule with constrained dihedral angle?
MMFFs_FF.Minimize(maxIts=10)
print "after first min"
# check values
for t in torsionList_noH:
 print t, rdMolTransforms.GetDihedralDeg(c, t[0], t[1], t[2], t[3])
MMFFs_FF.Minimize(maxIts=10)
print "after second min"
# check values
for t in torsionList_noH:
 print t, rdMolTransforms.GetDihedralDeg(c, t[0], t[1], t[2], t[3])

Jose Manuel

Hs=[]
for H in temp_Hs:
 Hs.append(H[0])

--
___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss