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 = 'OCCCC1=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=9999999999999.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
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
