Dear Greg,
I found my mistake I need the bond angle not the torsion.
I add this function for that...
from numpy import (array, dot, arccos, clip,zeros,degrees)
from numpy.linalg import norm
def cangle(v1,v2):
c=dot(v1,v2)/norm(v1)/norm(v2) # -> cosine of the angle
return arccos(clip(c, -1, 1)) # if you really want the angle
# THIS CODE IS WORKING BUT IT'S SLOW!
def AnglesBond(mol):
angles = mol.GetSubstructMatches(Chem.MolFromSmarts('*~*~*'))
conf = mol.GetConformer()
A60=0
A90=0
A102=0
for ang in angles:
a1 = conf.GetAtomPosition(ang[0])
c = conf.GetAtomPosition(ang[1])
a2 = conf.GetAtomPosition(ang[2])
v1=array([a1.x-c.x,a1.y-c.y,a1.z-c.z])
v2=array([a2.x-c.x,a2.y-c.y,a2.z-c.z])
Angle = abs(degrees(cangle(v1,v2)))
if Angle<=60:
A60+=1
elif Angle>60 and Angle<=90:
A90+=1
elif Angle>90 and Angle<=102:
A102+=1
return (A60,A90,A102)
the stats are not perfect but in progress.
[cid:046142a7-c63a-41fe-9e15-0ddd4ea26d62]
BR,
Guillaume
________________________________
De : Greg Landrum <greg.land...@gmail.com>
Envoyé : mercredi 14 septembre 2016 14:14
À : Guillaume GODIN
Cc : RDKit Discuss
Objet : Re: [Rdkit-discuss] Angstroms Hydrogen bonding
On Wed, Sep 14, 2016 at 4:16 AM, Guillaume GODIN
<guillaume.go...@firmenich.com<mailto:guillaume.go...@firmenich.com>> wrote:
Your solution is perfect!
glad it worked
I am currently implementing this article
"http://www.mdpi.com/1420-3049/20/10/18279"; using RDKit.
Interesting, I hadn't seen that one.
It's now almost done, I need to check the results on Heat of Formation now and
understand why my code to get Angle in torsions is not working.
def Angles(mol):
tors =
mol.GetSubstructMatches(Chem.MolFromSmarts('[C]-[C;O;S]-[C;O;S]-[C]'))
conf = mol.GetConformer()
A60=0
A90=0
A102=0
for tor in tors:
Angle = abs(AllChem.GetDihedralDeg(conf,tor[0], tor[1], tor[2], tor[3]))
if Angle<=60:
A60+=1
elif Angle>60 and Angle<=90:
A90+=1
elif Angle>90 and Angle<=102:
A102+=1
return (A60,A90,A102)
Comparing to the publication on Heat of Formation page 18294, I have to much
more occurences (ie. a factor of 50-100 more then expected).
I suspect that the problem is that the SMARTS you specify:
'[C]-[C;O;S]-[C;O;S]-[C]' won't actually match anything. [C;O;S] means
'aliphatic C AND aliphatic O AND aliphatic S'. It think you like mean
'[C,O,S]', which is 'aliphatic C OR aliphatic O OR aliphatic S'.
You might also want to think about whether or not you should replace C, O, S
(all of which only match aliphatic elements) with #6, #8, #16 (which can match
aliphatic or aromatic) throughout the SMARTS
When it will be done I will share the code with you guys.
Cool, I look forward to trying that out.
Best,
-greg
BR,
Guillaume
________________________________
De : Greg Landrum <greg.land...@gmail.com<mailto:greg.land...@gmail.com>>
Envoyé : mercredi 14 septembre 2016 12:16
À : Guillaume GODIN
Cc : RDKit Discuss
Objet : Re: [Rdkit-discuss] Angstroms Hydrogen bonding
Hi Guillaume,
On Tue, Sep 13, 2016 at 10:12 PM, Guillaume GODIN
<guillaume.go...@firmenich.com<mailto:guillaume.go...@firmenich.com>> wrote:
1 Does 3D coordinates of a conformer is in Angstroms ?
If you read the conformer from a file, for example a mol file, then the 3D
coordinates are in whatever units they were in in that file. This is usually
Angstrsom.
If you generate the conformation using one of the RDKit embedding functions,
then they are certainly in Angstroms.
2 How to enumerate all HBonding to determine the bond length ?
Interesting question. Here's a python function that might be a starting point:
def
findHBonds(m,confId=-1,possiblePartners='[#8,#7]',possibleHs='[#1][#8,#7]',distThresh=2.5):
conf = m.GetConformer(confId)
partners =[x[0] for x in
m.GetSubstructMatches(Chem.MolFromSmarts(possiblePartners))]
hs= [x[0] for x in m.GetSubstructMatches(Chem.MolFromSmarts(possibleHs))]
res = []
for h in hs:
ph = conf.GetAtomPosition(h)
for partner in partners:
if m.GetBondBetweenAtoms(h,partner) is not None:
continue
d = conf.GetAtomPosition(partner).Distance(ph)
if d<=distThresh:
res.append((h,partner,d))
return res
In order to allow flexibility about what an H bond is, I left the definitions
of the acceptors (partners in the above code) and polar Hs (just Hs in the
above code) as SMARTS definitions so that they can be customized.
**********************************************************************
DISCLAIMER
This email and any files transmitted with it, including replies and forwarded
copies (which may contain alterations) subsequently transmitted from Firmenich,
are confidential and solely for the use of the intended recipient. The contents
do not represent the opinion of Firmenich except to the extent that it relates
to their official business.
**********************************************************************
------------------------------------------------------------------------------
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
