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