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
