From e110a5c47939762c9106c3c93d01acd98e48081c Mon Sep 17 00:00:00 2001 From: Greg Landrum Date: Mon, 3 Oct 2022 10:27:22 +0200 Subject: [PATCH] Add GetBestTFDBetweenMolecules() (#5577) * a bit of cleanup * add GetBestTFDBetweenMolecules * a bit more cleanup * changes in response to review --- rdkit/Chem/TorsionFingerprints.py | 70 ++++++++++++++++++++++++------- rdkit/Chem/UnitTestMol3D.py | 9 ++++ 2 files changed, 64 insertions(+), 15 deletions(-) diff --git a/rdkit/Chem/TorsionFingerprints.py b/rdkit/Chem/TorsionFingerprints.py index 3f535e2e4..bd19ac0f0 100644 --- a/rdkit/Chem/TorsionFingerprints.py +++ b/rdkit/Chem/TorsionFingerprints.py @@ -566,20 +566,18 @@ def GetTFDBetweenConformers(mol, confIds1, confIds2, useWeights=True, maxDev='eq tfd = [] if useWeights: weights = CalculateTorsionWeights(mol, ignoreColinearBonds=ignoreColinearBonds) - for t1 in torsions1: - for t2 in torsions2: - tfd.append(CalculateTFD(t1, t2, weights=weights)) else: - for t1 in torsions1: - for t2 in torsions2: - tfd.append(CalculateTFD(t1, t2)) + weights = None + for t1 in torsions1: + for t2 in torsions2: + tfd.append(CalculateTFD(t1, t2, weights=weights)) return tfd def GetTFDBetweenMolecules(mol1, mol2, confId1=-1, confId2=-1, useWeights=True, maxDev='equal', symmRadius=2, ignoreColinearBonds=True): """ Wrapper to calculate the TFD between two molecules. - Important: The two molecules must be instances of the same molecule + Important: The two molecules must be isomorphic Arguments: - mol1: first instance of the molecule of interest @@ -611,11 +609,54 @@ def GetTFDBetweenMolecules(mol1, mol2, confId1=-1, confId2=-1, useWeights=True, torsion2 = CalculateTorsionAngles(mol2, tl, tlr, confId=confId2) if useWeights: weights = CalculateTorsionWeights(mol1, ignoreColinearBonds=ignoreColinearBonds) - tfd = CalculateTFD(torsion1, torsion2, weights=weights) else: - tfd = CalculateTFD(torsion1, torsion2) + weights = None + tfd = CalculateTFD(torsion1, torsion2, weights=weights) return tfd +def GetBestTFDBetweenMolecules(mol1, mol2, confId1=-1, useWeights=True, maxDev='equal', + symmRadius=2, ignoreColinearBonds=True): + """ Wrapper to calculate the best TFD between a single conformer of mol1 and all the conformers of mol2 + Important: The two molecules must be isomorphic + + Arguments: + - mol1: first instance of the molecule of interest + - mol2: second instance the molecule of interest + - confId1: conformer index for mol1 (default: first conformer) + - useWeights: flag for using torsion weights in the TFD calculation + - maxDev: maximal deviation used for normalization + 'equal': all torsions are normalized using 180.0 (default) + 'spec': each torsion is normalized using its specific + maximal deviation as given in the paper + - symmRadius: radius used for calculating the atom invariants + (default: 2) + - ignoreColinearBonds: if True (default), single bonds adjacent to + triple bonds are ignored + if False, alternative not-covalently bound + atoms are used to define the torsion + + Return: TFD value + """ + if (Chem.MolToSmiles(mol1) != Chem.MolToSmiles(mol2)): + raise ValueError("The two molecules must be instances of the same molecule!") + mol2 = _getSameAtomOrder(mol1, mol2) + tl, tlr = CalculateTorsionLists(mol1, maxDev=maxDev, symmRadius=symmRadius, + ignoreColinearBonds=ignoreColinearBonds) + # first molecule + torsion1 = CalculateTorsionAngles(mol1, tl, tlr, confId=confId1) + if useWeights: + weights = CalculateTorsionWeights(mol1, ignoreColinearBonds=ignoreColinearBonds) + else: + weights = None + best = 1e8 + for conf in mol2.GetConformers(): + # second molecule + torsion2 = CalculateTorsionAngles(mol2, tl, tlr, confId=conf.GetId()) + tfd = CalculateTFD(torsion1, torsion2, weights=weights) + best = min(best,tfd) + return best + + def GetTFDMatrix(mol, useWeights=True, maxDev='equal', symmRadius=2, ignoreColinearBonds=True): """ Wrapper to calculate the matrix of TFD values for the @@ -652,11 +693,10 @@ def GetTFDMatrix(mol, useWeights=True, maxDev='equal', symmRadius=2, ignoreColin tfdmat = [] if useWeights: weights = CalculateTorsionWeights(mol, ignoreColinearBonds=ignoreColinearBonds) - for i in range(0, numconf): - for j in range(0, i): - tfdmat.append(CalculateTFD(torsions[i], torsions[j], weights=weights)) else: - for i in range(0, numconf): - for j in range(0, i): - tfdmat.append(CalculateTFD(torsions[i], torsions[j])) + weights = None + for i in range(0, numconf): + for j in range(0, i): + tfdmat.append(CalculateTFD(torsions[i], torsions[j], weights=weights)) + return tfdmat diff --git a/rdkit/Chem/UnitTestMol3D.py b/rdkit/Chem/UnitTestMol3D.py index 3a8b90495..e191f737e 100755 --- a/rdkit/Chem/UnitTestMol3D.py +++ b/rdkit/Chem/UnitTestMol3D.py @@ -100,6 +100,11 @@ class TestCase(unittest.TestCase): tfd = TorsionFingerprints.GetTFDBetweenMolecules(mol, mol2) self.assertAlmostEqual(tfd, 0.0691, 4) + # exactly equivalent to the above since the mols each only have one conformer + tfd = TorsionFingerprints.GetBestTFDBetweenMolecules(mol2,mol) + self.assertAlmostEqual(tfd,0.0691,4) + + mol.AddConformer(mol2.GetConformer(), assignId=True) mol.AddConformer(mol2.GetConformer(), assignId=True) tfd = TorsionFingerprints.GetTFDBetweenConformers(mol, confIds1=[0], confIds2=[1, 2]) @@ -109,6 +114,10 @@ class TestCase(unittest.TestCase): tfdmat = TorsionFingerprints.GetTFDMatrix(mol) self.assertEqual(len(tfdmat), 3) + tfd = TorsionFingerprints.GetBestTFDBetweenMolecules(mol2,mol) + self.assertAlmostEqual(tfd,0,4) + + def testTorsionFingerprintsAtomReordering(self): # we use the xray structure from the paper (JCIM, 52, 1499, 2012): 1DWD refFile = os.path.join(RDConfig.RDCodeDir, 'Chem', 'test_data', '1DWD_ligand.pdb')