Hi Thomas,
I think there are a couple of problems with the code here:
1) you aren't storing the confirmation of the conformer of mol that
produces the best alignment (cid in the above code)
2) you aren't keeping the best alignment since you repeatedly run the same
molecule instances through the alignment code. You could fix this by, after
the double conformer loop, adding something like:
AllChem.AlignMol(refMol, mol, prbCid=refConformerId,
refCid=molConformerId, atomMap=zip(refMatch,match))
where I have assumed that molConformerId is the variable you use to solve
problem 1.
The loop would become something like this:
minRmsd = 1000;
for refCid in refCids:
for cid in cids:
rmsd = AllChem.AlignMol(refMol, mol, prbCid=refCid,
refCid=cid, atomMap=zip(refMatch,match))
if rmsd minRmsd:
minRmsd = rmsd
refConformerId = refCid
molConformerId = cid
rmsd = AllChem.AlignMol(refMol, mol, prbCid=refConformerId,
refCid=molConformerId,
atomMap=zip(refMatch,match))
It's also not really right to compare the results of this, which uses the
MCS code to find a common set of atoms and then does an RMSD alignment of
those, and the Open3DAlign results, which use a fuzzier scheme to identify
the atom-atom mapping between the molecules.
-greg
On Thu, Jun 26, 2014 at 2:32 PM, Thomas Strunz beginn...@hotmail.de wrote:
I'm trying to align all conformers of 2 molecules (and keep the best ones)
using the python api by following some of the tutorials:
http://nbviewer.ipython.org/gist/greglandrum/4316435/Working%20in%203D.ipynb
http://nbviewer.ipython.org/github/greglandrum/rdkit_blog/blob/master/notebooks/Using%20ConstrainedEmbed.ipynb
However whatever I try the alignment is considerably different than the
one created using the Open 3D alignment node in Knime.
Currently the issue is that the molecules do not seem to be aligned
properly. There is some sort of small shift. See code below.
Best Regards,
Thomas
refCids = generateConformers(refMol, numConformers)
mcs = MCS.FindMCS([refMol,mol],
ringMatchesRingOnly=matchesRingOnly)
if mcs.completed == 1 and mcs.numAtoms 0:
core = Chem.MolFromSmarts(mcs.smarts)
logger.info('MCS: %s', Chem.MolToSmiles(core))
refMatch = refMol.GetSubstructMatch(core)
match = mol.GetSubstructMatch(core)
# conformers for current target
cids = generateConformers(mol, numConformers,
coordMap=coordMap)
minRmsd = 1000;
for refCid in refCids:
for cid in cids:
rmsd = AllChem.AlignMol(refMol, mol, prbCid=refCid,
refCid=cid, atomMap=zip(refMatch,match))
logger.debug('RMSD: %.2f', rmsd)
if rmsd minRmsd:
logger.debug('New min RMSD: %.2f', rmsd)
minRmsd = rmsd
refConformerId = refCid
def generateConformers(mol, numConformers):
AllChem.EmbedMolecule(mol)
AllChem.MMFFOptimizeMolecule(mol)
cids=AllChem.EmbedMultipleConfs(mol, numConfs=numConformers,
maxAttempts=50, pruneRmsThresh=0.5, coordMap=coordMap)
for cid in cids: AllChem.MMFFOptimizeMolecule(mol,confId=cid)
return cids
--
Open source business process management suite built on Java and Eclipse
Turn processes into business applications with Bonita BPM Community Edition
Quickly connect people, data, and systems into organized workflows
Winner of BOSSIE, CODIE, OW2 and Gartner awards
http://p.sf.net/sfu/Bonitasoft
___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
--
Open source business process management suite built on Java and Eclipse
Turn processes into business applications with Bonita BPM Community Edition
Quickly connect people, data, and systems into organized workflows
Winner of BOSSIE, CODIE, OW2 and Gartner awards
http://p.sf.net/sfu/Bonitasoft___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
