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

2015-10-20 Thread Jose Manuel Gally 
Dear Paolo,

thank you for your explanation, very useful as always!

Another question if I'm allowed, I would like to generate multiple rings 
conformations per molecule.
However, from what I understand from the protocol implemented in RDKit, 
I must first generate "n" starting conformers.
Those are created using a distance geometry matrix and not in a force field.

Would it be possible to add conformations to the molecule while keeping 
the constraints on linkers and side chains?

Then I would be clustering the conformers minimized with constraints to 
keep only conformations that are actually different

Best,
Jose Manuel

On 10/20/2015 04:01 PM, Paolo Tosco wrote:
> Dear Jose Manuel,
>
> actually each constraint you set results in an additional term in the
> force-field equation; furthermore, the more constraints, the more
> trouble for the minimizer to converge. In general, constraints should be
> kept to a minimum; you won't need constraints on aromatics, since they
> are already kept rigid and planar by torsion and oop terms.
>
> Cheers,
> p.
>
> On 20/10/2015 14:38, Jose Manuel Gally wrote:
>> 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
> --
> ___
> Rdkit-discuss mailing list
> Rdkit-disc

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

2015-10-20 Thread Paolo Tosco 
Dear Jose Manuel,

actually each constraint you set results in an additional term in the 
force-field equation; furthermore, the more constraints, the more 
trouble for the minimizer to converge. In general, constraints should be 
kept to a minimum; you won't need constraints on aromatics, since they 
are already kept rigid and planar by torsion and oop terms.

Cheers,
p.

On 20/10/2015 14:38, Jose Manuel Gally wrote:
> 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


--
___
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

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("[!$(*#*)&X3&R]-!@[!$(*#*)&X3&R]")
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 Paolo Tosco 
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('CC')
   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=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

[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