Dear Rafal,
On 22/04/2016 17:15, Rafal Roszak wrote:
> Hello,
>
> I want to find:
> A) "global minimum" [*] for given compound(s) and
> B) (for the same compound) minimum with constraint(s) (e.g. frozen
> angle)
>
> For problem A) I tried folowing code:
>
> self.MOLEC=Chem.MolFromSmiles(SMILES)
> self.MOLEC=Chem.AddHs(self.MOLEC)
> allconf=Chem.rdDistGeom.EmbedMultipleConfs(self.MOLEC,
> numConfs=HowMany, enforceChirality=True, pruneRmsThresh=limitRMS) for
> confId in range(len(allconf)): AllChem.UFFOptimizeMolecule(self.MOLEC,
> maxIters=12540, confId=confId) self.ff =
> AllChem.UFFGetMoleculeForceField(self.MOLEC) energy_value = str
> (self.ff.CalcEnergy()) print (energy_value)
>
Here you do generate and optimize many conformations, but then you
retrieve the force-field with UFFGetMoleculeForceField() always for
conformation -1 (the default), hence you get the same energy.
I slightly modified your script (see below) to actually generate a few
dozens of diverse conformations:
#!/usr/bin/env python
import sys
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
SMILES = 'CCOCCn1c(C2CC[NH+](CCc3ccc(C(C)(C)C(=O)[O-])cc3)CC2)nc2ccccc21'
HowMany = 100
limitRMS = 2.0
mol = Chem.MolFromSmiles(SMILES)
mol = AllChem.AddHs(mol)
allconf = AllChem.EmbedMultipleConfs(mol,
numConfs = HowMany, enforceChirality=True, pruneRmsThresh = limitRMS)
w = Chem.SDWriter(sys.stdout)
for confId in allconf:
ff = AllChem.UFFGetMoleculeForceField(mol, confId = confId)
ff.Minimize()
energy_value = ff.CalcEnergy()
mol.SetProp('ENERGY', '{0:.2f}'.format(energy_value))
w.write(mol, confId = confId)
w.close()
> No matter what values are used for the numConfs or pruneRmsThresh after
> optimization I end up in conformer(s) with the same energy. The energy
> value vary from one run to another. In other words EmbedMultipleConfs()
> generate random structure (that make sense to me) but all conformation
> are very similar - after optimization all fall into one structure. Is
> it bug or feature? :) How effectively generate set of _different_
> conformation? I found two solution: a) execute presented above code
> multiple times b) set dihedral angle to unnatural value (using
> Chem.rdMolTransforms.SetDihedralDeg() function) for selected dihedral
> angle and then minimize structure Both works but seems to me that both
> are not optimal.
>
>
> For problem B:
>
> I want to minimize structure with frozen dihedral angle. I set the
> angle value using Chem.rdMolTransforms.SetDihedralDeg I could not find
> any possibility to freeze angle so I used
> rdForceField.ForceField.MMFFAddTorsionConstraint function. with very
> high forceConstant (relative=True minDihedralDeg=-5, maxDihedralDeg=5,
> forceConstant=9999999.9). The code is as follow:
>
>
> conformers=Chem.rdDistGeom.EmbedMultipleConfs(self.MOLEC,numConfs=confNum,
> enforceChirality=enforceChirality,numThreads=numThreads,
> pruneRmsThresh=RMSThreshold)
> for confId in range(len(conformers)):
> ff = AllChem.UFFGetMoleculeForceField(self.MOLEC, confId=confId)
> rdForceField.ForceField.MMFFAddTorsionConstraint( ff, idx0, idx1,
> idx2, idx3, True, minDihedralDeg=-5, maxDihedralDeg=5,
> forceConstant=9999999.9)
> Chem.rdMolTransforms.SetDihedralDeg(self.MOLEC.GetConformer(confId),
> idx0, idx1, idx2, idx3, 0.0 )
> AllChem.UFFOptimizeMolecule(self.MOLEC, maxIters=12540, confId=confId)
>
> but it does not work - i end up in structure with dihedral value far from
> this which was set. How can I force optimization with defined value (or some
> range of values)? In other words what I am doing wrong?
The problem here is that you are setting the constraint in MMFF and then
use UFF for minimization; furthermore MMFFAddTorsionConstraint() and
UFFAddTorsionConstraint() are member functions of the respective
force-fields. Please find below a sample script which should do what you
describe; here I read the value of the dihedral before and after
minimization, to verify that it is actually frozen.
#!/usr/bin/env python
import sys
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
SMILES = 'CCOCCn1c(C2CC[NH+](CCc3ccc(C(C)(C)C(=O)[O-])cc3)CC2)nc2ccccc21'
HowMany = 100
limitRMS = 2.0
mol = Chem.MolFromSmiles(SMILES)
mol = AllChem.AddHs(mol)
tors = mol.GetSubstructMatch(Chem.MolFromSmarts('cCC[NH+]'))
allconf = AllChem.EmbedMultipleConfs(mol,
numConfs = HowMany, enforceChirality=True, pruneRmsThresh = limitRMS)
w = Chem.SDWriter(sys.stdout)
for confId in allconf:
ff = AllChem.UFFGetMoleculeForceField(mol, confId = confId)
before = AllChem.GetDihedralDeg(mol.GetConformer(confId),
tors[0], tors[1], tors[2], tors[3])
ff.UFFAddTorsionConstraint(tors[0], tors[1], tors[2], tors[3],
True, -5, 5, 9999)
ff.Minimize()
energy_value = ff.CalcEnergy()
after = AllChem.GetDihedralDeg(mol.GetConformer(confId),
tors[0], tors[1], tors[2], tors[3])
mol.SetProp('DIHEDRAL', '{0:.2f},{1:.2f}'.format \
(before, after))
mol.SetProp('ENERGY', '{0:.2f}'.format(energy_value))
w.write(mol, confId = confId)
w.close()
Kind regards,
Paolo
>
> Regards,
>
> Rafal
>
>
>
>
> [*] Yes I aware that it is very complex task to find a true global
> minima but I only want to find reasonably good/stable conformer, that
> why I wrote "global minimum"
>
> ------------------------------------------------------------------------------
> Find and fix application performance issues faster with Applications Manager
> Applications Manager provides deep performance insights into multiple tiers of
> your business applications. It resolves application problems quickly and
> reduces your MTTR. Get your free trial!
> https://ad.doubleclick.net/ddm/clk/302982198;130105516;z
> _______________________________________________
> Rdkit-discuss mailing list
> Rdkit-discuss@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
------------------------------------------------------------------------------
Find and fix application performance issues faster with Applications Manager
Applications Manager provides deep performance insights into multiple tiers of
your business applications. It resolves application problems quickly and
reduces your MTTR. Get your free trial!
https://ad.doubleclick.net/ddm/clk/302982198;130105516;z
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
