Dear Michal,
please find attached a small script which accomplishes what you describe
by a different approach, i.e. it minimizes only the methyl group in
2-methylthiophene while keeping the rest fixed, effectively pushing it
back in plane. Would that work for you?
Best,
Paolo
On 04/15/2015 10:57 AM, Michal Krompiec wrote:
Hello,
I'm trying to manipulate out-of-plane bends of substituents attached
to aromatic rings. Obviously, they should lie in the plane of the
ring, but sometimes (after constrained optimization) they come
slightly out of plane and there doesn't seem to be a way to
push/rotate them back to the ring plane. Adding constraints to the
forcefield (which seems to work for out-of-plane bends as well) is not
a perfect solution.
For example, there is no way to manipulate the out-of-plane bend
(dihedral angle) of the methyl group in the 2-methylthiophene molecule
c1cc(C)sc1, defined by the first four atoms c1cc(C). If you try to use
rdMolTransforms.SetDihedralDeg, it will raise an exception because the
bond around it tries to rotate is part of a ring (line 375 in
MolTransforms.cpp) - and for a good reason, because what I wanted to
do is to rotate just the methyl group (atom l as defined in
setDihedralRad, and everything attached to it), not all atoms
connected to atom k being on the "right hand side" of the j-k bond.
Would it make sense to make a modified version of setDihedralRad():
setOutOfPlaneBendRad(), with the following modifications:
* omit line 375: if(queryIsBondInRing(bondJK)) throw
ValueErrorException("bond (j,k) must not belong to a ring"); instead
perhaps check if bond k,l belongs to a ring
* change line 401: _toBeMovedIdxList(mol, jAtomId, kAtomId, alist); to
_toBeMovedIdxList(mol, kAtomId, lAtomId, alist); (we want to rotate
the tree rooted in atom l, not atom k)
Best regards,
Michal
------------------------------------------------------------------------------
BPM Camp - Free Virtual Workshop May 6th at 10am PDT/1PM EDT
Develop your own process in accordance with the BPMN 2 standard
Learn Process modeling best practices with Bonita BPM through live exercises
http://www.bonitasoft.com/be-part-of-it/events/bpm-camp-virtual- event?utm_
source=Sourceforge_BPM_Camp_5_6_15&utm_medium=email&utm_campaign=VA_SF
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
#!/usr/bin/env python
import rdkit
from rdkit import Chem
from rdkit.Chem import ChemicalForceFields, rdForceFieldHelpers
oopThiophene = \
"""oopThiophene
3D
21 21 0 0 0 0 0 0 0 0999 V2000
-3.7316 1.2367 0.7532 S 0 0 0 0 0 0 0 0 0 0 0 0
-2.6758 0.5095 -0.4904 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.4176 0.8536 -0.3142 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2002 1.7211 0.8189 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.2999 2.0145 1.4780 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.3754 2.0748 3.1288 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.1713 -0.4058 -1.5990 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.4116 0.0803 -2.3745 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.7054 -0.1744 -1.5816 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.6196 0.5172 -0.9595 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.2262 2.0941 1.0990 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.4091 2.0863 3.4416 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.8834 1.2076 3.5436 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.8833 2.9701 3.4792 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.3665 -0.5309 -2.3083 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.4122 -1.3608 -1.1559 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.3165 1.1398 -2.5607 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.4667 -0.4472 -3.3153 H 0 0 0 0 0 0 0 0 0 0 0 0
-6.5592 0.0279 -2.2113 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.7326 0.4748 -0.7189 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.7325 -1.2044 -1.2581 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 2 0
3 4 1 0
4 5 2 0
5 1 1 0
5 6 1 0
2 7 1 0
7 8 1 0
8 9 1 0
3 10 1 0
4 11 1 0
6 12 1 0
6 13 1 0
6 14 1 0
7 15 1 0
7 16 1 0
8 17 1 0
8 18 1 0
9 19 1 0
9 20 1 0
9 21 1 0
M END"""
mol = Chem.MolFromMolBlock(oopThiophene, sanitize = True, removeHs = False)
methyl = Chem.MolFromSmarts('C([H])([H])([H])')
match = mol.GetSubstructMatch(methyl)
mp = ChemicalForceFields.MMFFGetMoleculeProperties(mol)
ff = ChemicalForceFields.MMFFGetMoleculeForceField(mol, mp)
for i in range(mol.GetNumAtoms()):
if (not (i in match)):
ff.AddFixedPoint(i)
ff.Minimize()
w = Chem.SDWriter('minOopThiophene.sdf')
w.write(mol)
w.close()
------------------------------------------------------------------------------
BPM Camp - Free Virtual Workshop May 6th at 10am PDT/1PM EDT
Develop your own process in accordance with the BPMN 2 standard
Learn Process modeling best practices with Bonita BPM through live exercises
http://www.bonitasoft.com/be-part-of-it/events/bpm-camp-virtual- event?utm_
source=Sourceforge_BPM_Camp_5_6_15&utm_medium=email&utm_campaign=VA_SF
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
