Re: [Rdkit-discuss] RDKit conformer handling possible bug?

[Greg L. <gre...@gm...>]

Hi Steve,

[A general request: please send RDKit questions/comments to the
rdkit-discuss mailing list so that others can see the questions and
answers. I'm CC'ing the list on my reply because I think it could be of
general interest]

For this one, what you have encountered here is a feature, but it's not
immediately obvious.

There are two things going on here:
1) an RDKit molecule can have 0-N conformers. Chem.MolToMolBlock(), by
default, picks the one with ID zero. If there are no conformers, then it
just outputs zeros for the coordinates. Mol.AddConformer() does not replace
existing conformers, but just adds to the list.
2) By default Mol.AddConformer() does not change the ID of the conformer it
adds. So if you add the default conformer (number zero) from one molecule
to another, you end up with two conformers with ID 0. You can work around
this by calling Mol.AddConformer() with assignId=True

Here's a demo:
mol_SDF_H_2=Chem.AddHs(Chem.MolFromMolBlock(molString))
print(mol_SDF_H_2.GetNumConformers())
# output is 1
mol_SDF_H_2.AddConformer(mol_SDF_H.GetConformer(0),assignId=True)
print(mol_SDF_H_2.GetNumConformers())
# output is 2
print([x.GetId() for x in mol_SDF_H_2.GetConformers()])
# output is [0,1]
print(Chem.MolToMolBlock(mol_SDF_H_2,confId=1))
# output is what is expected.

If you want to replace the existing conformer, you can do:
mol.RemoveAllConformers() before you call mol.AddConformer().

Hope this helps,
-greg




On Tue, Jan 10, 2017 at 2:58 PM, Stephen Roughley <S.R...@ve...>
wrote:

> from rdkit import Chem
> from rdkit.Chem import AllChem
>
> #From SMILES - this works!
> mol_SMI = Chem.MolFromSmiles("C")
> mol_SMI_H=Chem.AddHs(mol_SMI)
> mol_SMI_H_2 = Chem.AddHs(Chem.MolFromSmiles("C"))
>
> ids=AllChem.EmbedMultipleConfs(mol_SMI_H, numConfs=10)
> mol_SMI_H_2.AddConformer(mol_SMI_H.GetConformer(0))
> print(Chem.MolToMolBlock(mol_SMI_H_2))
>
> ## As expected, the conformer is transferred to the vanilla copy
> ##
> ##     RDKit          3D
> ##
> ##  5  4  0  0  0  0  0  0  0  0999 V2000
> ##   -0.0221    0.0032    0.0165 C   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.6665    0.8884   -0.1014 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.3778   -0.8578   -0.5883 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.0964   -0.3151    1.0638 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.9699    0.2812   -0.3906 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##  1  2  1  0
> ##  1  3  1  0
> ##  1  4  1  0
> ##  1  5  1  0
> ##M  END
>
> molString="""
>   Mrv16b2101101711412D
>
>   1  0  0  0  0  0            999 V2000
>     0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
> M  END
> $$
> """
>
> #Now try from MOL Block
> mol_SDF=Chem.MolFromMolBlock(molString)
> print(Chem.MolToMolBlock(mol_SDF))
>
> ## As expected we have a molecule with 1 atom and no conformers
> ##
> ##     RDKit          2D
> ##
> ##  1  0  0  0  0  0  0  0  0  0999 V2000
> ##    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
> ##M  END
> ##
>
> mol_SDF_H=Chem.AddHs(mol_SDF)
> Chem.MolFromMolBlock(molString)
> mol_SDF_H_2=Chem.AddHs(Chem.MolFromMolBlock(molString))
>
> ids=AllChem.EmbedMultipleConfs(mol_SDF_H, numConfs=10)
> print(Chem.MolToMolBlock(mol_SDF_H,confId=0))
>
> ## As expected, we have an molecule with Hs added, and conformers with
> coords
> ##
> ##     RDKit          3D
> ##
> ##  5  4  0  0  0  0  0  0  0  0999 V2000
> ##    0.0260   -0.0121   -0.0153 C   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.9804   -0.4982   -0.2778 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.0317    1.0327   -0.3649 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.1070    0.0168    1.0896 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.8677   -0.5393   -0.4316 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##  1  2  1  0
> ##  1  3  1  0
> ##  1  4  1  0
> ##  1  5  1  0
> ##M  END
>
> mol_SDF_H_2.AddConformer(mol_SDF_H.GetConformer(0))
> print(Chem.MolToMolBlock(mol_SDF_H_2))
> ## Whoa! This time, we have added a conformation to a vanilla copy
> ## exactly as we did for SMILES, but no conformer coords are exported:
> ##
> ##     RDKit          2D
> ##
> ##  5  4  0  0  0  0  0  0  0  0999 V2000
> ##    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##  1  2  1  0
> ##  1  3  1  0
> ##  1  4  1  0
> ##  1  5  1  0
> ##M  END
>
> mol_SMI_H_3 = Chem.AddHs(Chem.MolFromSmiles("C"))
> mol_SMI_H_3.AddConformer(mol_SDF_H.GetConformer(0))
> print(Chem.MolToMolBlock(mol_SMI_H_3))
> ## And again, this works ok...
> ##
> ##     RDKit          3D
> ##
> ##  5  4  0  0  0  0  0  0  0  0999 V2000
> ##    0.0260   -0.0121   -0.0153 C   0  0  0  0  0  0  0  0  0  0  0  0
> ##    0.9804   -0.4982   -0.2778 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.0317    1.0327   -0.3649 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.1070    0.0168    1.0896 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##   -0.8677   -0.5393   -0.4316 H   0  0  0  0  0  0  0  0  0  0  0  0
> ##  1  2  1  0
> ##  1  3  1  0
> ##  1  4  1  0
> ##  1  5  1  0
> ##M  END
>