Conformers generation - RMSD cutoff limit not applied?

[Fransu86]

I noticed that when running conformer generation in Datamol, the rms_cutoff option doesn't seem to be applied correctly.
Specifically:

import datamol as dm
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

smiles = "N(c1n[nH]c(C2CC2)c1)c1nc(Nc2ccc(CC#N)cc2)ncc1"
mol = dm.to_mol(smiles)

mol = dm.conformers.generate(
    mol,
    n_confs=200,
    rms_cutoff=0.25,
    minimize_energy=True,
    align_conformers=True,
    ewindow=5.0
)

n_confs = mol.GetNumConformers()
print(f"Number of Conformations: {n_confs}")

Should generate 145 aligned conformations and make sure that no two conformations have a pairwise RMSD equal or below 0.25.
Looking at them in 3D, I noticed that several seem to be very very similar and indeed with the following debug code:

rmsd_matrix = dm.conformers.rmsd(mol)
threshold = 0.25
mask_below_threshold = rmsd_matrix < threshold
total_count = np.sum(mask_below_threshold)
count_excluding_diagonal = total_count - n_confs
unique_pairs_count = count_excluding_diagonal // 2

It tells me that out of the 145 conformations, 50 of them have pairwise RMSD <0.25.
Just one example from the compounds above; the four conformations generated via Datamol have pairwise RMSD<<0.25.

Am I missing something on the role or the rms_cutoff, particularly together with the align_conformers option?
I have seen the same behavior with several other molecules.

[zhu0619]

@Fransu86
The clustering approach might not be reliable.
A quick work-around would be iterating the rsm matrix .

def remove_duplicate_conformers(mol, rmsd_threshold=0.25):
    """
    Removes duplicate conformers from a molecule based on an RMSD threshold.

    :param mol: The RDKit molecule object with multiple conformers.
    :param rmsd_threshold: The RMSD threshold for considering two conformers as duplicates (in Angstroms).
    :return: A new molecule object with only unique conformers.
    """
    if mol.GetNumConformers() <= 1:
        return mol
    
    # Calculate pairwise RMSD values
    num_confs = mol.GetNumConformers()
    # Use a symmetric matrix for storing RMSD to avoid redundant calculations
    rmsd_matrix = np.ones((num_confs, num_confs)) * -1.0
    for i in range(num_confs):
        for j in range(i + 1, num_confs):
            # The GetBestRMS function aligns the conformers optimally before calculating RMSD
            rmsd = AllChem.GetBestRMS(mol, mol, i, j)
            rmsd_matrix[i, j] = rmsd
            rmsd_matrix[j, i] = rmsd

    # Identify and keep unique conformers
    kept_conf_ids = []
    removed_conf_ids = set()

    for i in range(num_confs):
        if i in removed_conf_ids:
            continue
        
        kept_conf_ids.append(i)
        
        for j in range(i + 1, num_confs):
            if j in removed_conf_ids:
                continue
            
            if rmsd_matrix[i, j] < rmsd_threshold:
                removed_conf_ids.add(j)

    # Create a new molecule with only the kept conformers
    unique_mol = Chem.Mol(mol)
    unique_mol.RemoveAllConformers()
    for conf_id in kept_conf_ids:
        unique_mol.AddConformer(mol.GetConformer(conf_id), assignId=True)
        
    return unique_mol

[Fransu86]

The code does the trick! It removes the RMSD duplicates as I would have expected from the rms_cutoff flag.

[j-wags]

A complete drive-by comment here, but one frequent source of confusion I've observed with RMSDs are whether they're calculated for all atoms or just heavy atoms, and whether the molecule in question has implicit or explicit hydrogens when the various RMSDs are calculated.