Hello all,
I am trying to calculate 3D Descriptors following this publication:
"*Beyond Rotatable Bond Counts: Capturing 3D Conformational Flexibility in
a Single Descriptor"*, Jerome G. P. Wicker and Richard I. Cooper. J. Chem.
Inf. Model. 2016, 56, 2347−2352
I am essentially using the same script as they have in the supporting
information and i have attached it here as well. In Table 2 from the above
calculation, the value of the descriptor (nConf20) for ZINC000290539224
molecule is listed as 10. However, when I run the exact code as the one
they used, I get different value at each run.
I have already contacted the authors but got no response. I am wondering
if the code they have in the supporting information is not right or the
value they listed in the table is wrong?
The SMILES string for this particular molecule is:
'CC(C)N2CC(NCc1cnc(C(C)O)s1)CC2=O'
Thanks in advance for your help!
Ali
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from collections import OrderedDict
def GenerateConformers(smil, numConfs):
mol=Chem.MolFromSmiles(smil)
molecule = Chem.AddHs(mol)
conformerIntegers = []
conformers = AllChem.EmbedMultipleConfs (molecule,numConfs,pruneRmsThresh=0.5, numThreads =3)
optimised_and_energies = AllChem.MMFFOptimizeMoleculeConfs(molecule , maxIters=600, numThreads =3, nonBondedThresh =100.0)
EnergyDictionaryWithIDAsKey = {}
FinalConformersToUse = {}
for conformer in conformers :
optimised , energy = optimised_and_energies[conformer]
if optimised == 0:
EnergyDictionaryWithIDAsKey [ conformer ] = energy
conformerIntegers.append(conformer)
lowestenergy = min( EnergyDictionaryWithIDAsKey.values () )
for k, v in (EnergyDictionaryWithIDAsKey.items()):
if v == lowestenergy:
lowestEnergyConformerID = k
FinalConformersToUse [lowestEnergyConformerID] = lowestenergy
molecule = AllChem.RemoveHs(molecule)
matches = molecule.GetSubstructMatches ( molecule , uniquify=False )
maps = [ list (enumerate(match)) for match in matches]
for conformerID in EnergyDictionaryWithIDAsKey.keys ():
okayToAdd = True
for finalconformerID in FinalConformersToUse.keys ():
RMS = AllChem.GetBestRMS ( molecule , molecule , finalconformerID , conformerID , maps)
if RMS< 1.0:
okayToAdd = False
break
if okayToAdd:
FinalConformersToUse [conformerID] = EnergyDictionaryWithIDAsKey [conformerID]
sortedDictionary =OrderedDict (sorted ((FinalConformersToUse.items()),key=lambda t: t[1] ))
energies = [val for val in (sortedDictionary.values())]
energy_descriptor = 0
relative_energies=np.array(energies)-energies[0]
for energy in relative_energies [1:]:
if 0 <= energy < 20:
energy_descriptor += 1
return energy_descriptor
sm='CC(C)N2CC(NCc1cnc(C(C)O)s1)CC2=O'
non_SMILES=['C(F)C(F)(F)', 'C(OC1=C2)=NC1=CC=C2C(OC3=C4)=NC3=CC=C4','C(C=C1)=CC=C1OC(C=C2)=CC=C2O',
'C(OC1=C2)=NC1=CC=C2C(OC3=C4)=NC3=CC=C4','CC(C(#N))CC(C(#N))','C(C=C1)=CC=C1OC(C=C2)=CC=C2O',
'CC(C=C1)=CC=C1C(C=C2)=CC=C2O','NC(=O)NC(=O)']
print (GenerateConformers('CC(C)N2CC(NCc1cnc(C(C)O)s1)CC2=O',50))------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
