Dear Greg,
Your solution is perfect!
BTW, I was also looking for the intramolecular H-H distances, and your function
can be a modify to do this too:
def
findintramolecularHH(m,confId=-1,possiblePartners='[#1][*]',possibleHs='[#1]',distThresh=2.5):
conf = m.GetConformer(confId)
partners =[x[0] for x in
m.GetSubstructMatches(Chem.MolFromSmarts(possiblePartners))]
linkers =[x[1] 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,linker in zip(partners,linkers):
if m.GetBondBetweenAtoms(h,linker) is not None:
continue
d = conf.GetAtomPosition(partner).Distance(ph)
if d<=distThresh:
res.append((h,partner,d))
return res
Thanks for your help. two issues solved on one comment, great!
I am currently implementing this article
"http://www.mdpi.com/1420-3049/20/10/18279"; using RDKit.
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).
When it will be done I will share the code with you guys.
BR,
Guillaume
________________________________
De : Greg Landrum <[email protected]>
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
<[email protected]<mailto:[email protected]>> 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
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
