Memory overflow with GetBestRMS when removeHs=False

[jandom]

By memory overflow I don't mean an actual exception: my laptop just stops responding when its 20GB of ram become exhausted ;]

The problem is super-simple: i load a molecule from mol file,

mol = Chem.MolFromMolFile(f, removeHs=False)
ref = Chem.MolFromMolFile(f, removeHs=False)

no hydrogens removed, and try to rms align the two mol objects.

rms = AllChem.GetBestRMS(ref, mol) 

This explodes - memory is just filled up until my machine crashes. I looked at the molecules visually and they seem totally fine. The problem disappears when i set removeHs=True.

Example script and inputs to reproduce the problem are accessible via link below. This test case was created with rdkit from 94aef1f

https://drive.google.com/file/d/0BzI3NK6qw0lJV0RGUU1EM3ZORm8/edit?usp=sharing

[greglandrum]

GetBestRMS() has a "small" problem with combinatorial explosions when there are Hs in the molecule. There are some notes about this, and suggestions about what to do, here:
https://www.mail-archive.com/rdkit-discuss@lists.sourceforge.net/msg03417.html

The short description of the problem is that every CH3 increases the number of possible alignments by a factor of 3, and every CH2 group doubles the work.

Something as simple as butane has 288 possible mappings that need to be explored when Hs are taken into account:

In [2]: m= Chem.AddHs(Chem.MolFromSmiles('CCCC'))

In [3]: len(m.GetSubstructMatches(m,uniquify=False))
Out[3]: 288

[jandom]

Oh, ok. So the RMS in rdkit attempts all the atom-order permutations? My expectation (from protein alignment) was that even structures with multiple thousands of atoms can be aligned in a flash but that's provided atom order is preserved.

[greglandrum]

Yeah, AllChem.GetBestRMS() tries all possible atom-order permutations that are compatible with the topological symmetry. If you know the atom order you can just use AllChem.AlignMol(); this only does one alignment and is much faster.