Re: [Rdkit-discuss] out-of-plane bends

2015-04-16 Thread Michal Krompiec
Hi Paolo,
Thanks, this is an interesting feature - may come useful one day!
I found a different solution - via torsion constraints:

mp = ChemicalForceFields.MMFFGetMoleculeProperties(mol)
ff = ChemicalForceFields.MMFFGetMoleculeForceField(mol, mp)
sp=Chem.MolFromSmarts(c1ccsc1!@c) #thiophene connected to another
aromatic ring

atoms=(1,0,4,5) #out-of-plane bend, atoms 1, 0 and 4 are in the thiophene
ring, atom 5 is the first atom of the next ring

maplist = mol.GetSubstructMatches(sp)

if (len(maplist)0):
for match in maplist :

angle=180.0 #desired dihedral angle
a=[]
for i in range (0,4) :
a.append(match[atoms[i]])
cur_angle=rdMolTransforms.GetDihedralDeg(mol.GetConformer(),
a[0],a[1],a[2], a[3])

if abs(abs(cur_angle)-abs(angle))45.0: # closer to 0 than to
180?
angle=0.0 # freeze at 0
   ff.MMFFAddTorsionConstraint(a[0],a[1],a[2],a[3], False, angle, angle,
1e4)

ff.Minimize()


Best wishes,

Michal

On 15 April 2015 at 23:06, Paolo Tosco paolo.to...@unito.it wrote:

  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 
 exerciseshttp://www.bonitasoft.com/be-part-of-it/events/bpm-camp-virtual- 
 event?utm_
 source=Sourceforge_BPM_Camp_5_6_15utm_medium=emailutm_campaign=VA_SF



 ___
 Rdkit-discuss mailing 
 listRdkit-discuss@lists.sourceforge.nethttps://lists.sourceforge.net/lists/listinfo/rdkit-discuss




 --
 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_15utm_medium=emailutm_campaign=VA_SF
 ___
 Rdkit-discuss mailing list
 Rdkit-discuss@lists.sourceforge.net
 https://lists.sourceforge.net/lists/listinfo/rdkit-discuss


--
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_15utm_medium=emailutm_campaign=VA_SF___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss


Re: [Rdkit-discuss] out-of-plane bends

2015-04-15 Thread Paolo Tosco

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_15utm_medium=emailutm_campaign=VA_SF


___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
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.73161.23670.7532 S   0  0  0  0  0  0  0  0  0  0  0  0
   -2.67580.5095   -0.4904 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.41760.8536   -0.3142 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.20021.72110.8189 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.29992.01451.4780 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.37542.07483.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.41160.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.61960.5172   -0.9595 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.22622.09411.0990 H   0  0  0  0  0  0  0  0  0  0  0  0
   -3.40912.08633.4416 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.88341.20763.5436 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.88332.97013.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.31651.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.55920.0279   -2.2113 H   0  0  0  0  0  0  0  0  0  0  0  0
   -5.73260.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 

[Rdkit-discuss] out-of-plane bends

2015-04-15 Thread Michal Krompiec
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_15utm_medium=emailutm_campaign=VA_SF___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss