Dear Méliné,
if you wish to keep input hydrogens you need to add a "removeHs = False"
parameter to Chem.SDMolSupplier(); hydrogens are removed by default.
As Peter pointed out, the original coordinates may give large energy
values just because of small deviations from ideal values of bond
length, angles and dihedrals according to the specific force-field you
are using, which could be fixed by a geometry optimizazion.
Additionally, the absolute force-field energy value is not particularly
meaningful per se: you will get high values from molecules which contain
small strained rings, for instance, which does not indicate that their
geometry is not reasonable.
Finally, you should consider that the RDKit implementation of UFF does
not include the electrostatic term, which may give rise to significantly
different geometries/energies, particularly with molecules characterized
by opposite formal charges or intramolecular hydrogen bonds.
Cheers,
p.
On 02/09/17 16:48, Méliné Simsir wrote:
Dear all,
I'm still a beginner with Rdkit, and I’m having a hard time
calculating the energy of conformers.
I think I have read all the topics on the subject, but I might have
missed something.
Here is a shorten version of a code I wrote thanks to all the
information I could find here.
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdmolops
mols = Chem.SDMolSupplier('test_3c.sdf')
for mol in mols:
# UFF
ffu = AllChem.UFFGetMoleculeForceField(mol)
energy_value_U = ffu.CalcEnergy()
print(energy_value_U)
# MMFF
mol = rdmolops.AddHs(mol, addCoords = True)
mp = AllChem.MMFFGetMoleculeProperties(mol)
ffm = AllChem.MMFFGetMoleculeForceField(mol, mp)
energy_value_M = ffm.CalcEnergy()
print(energy_value_M)
My problem is that I get weird energy values (not all the time I
think, but too often). Either too high or too low. And i don't get why.
I would like to get the energy of my ligand as he is. Without changing
it's conformation at all. That's why I don't use the minimization
function, or optimization. First I also didn't add the hydrogen, but
thanks to one of Paolo Tosco's code (mmffEnergyTerms.py), I saw that H
are not included even if they are in the sdf file. But, I still have
some weird values i think.
Here three examples:
GNP_1CIP : -305.2896 kcal/mol (MMFF, with H added)
-46.4812 kcal/mol (MMFF, without H add)
278.6606 kcal/mol (UFF)
SAM_2IGT: 1675.9951 (MMFF, with H added)
151.6481 (MMFF, without H add)
217.1217 (UFF)
ACP_3A1C: 222.9396 (MMFF, with H added)
269.151 (MMFF, without H add)
217.1217 (UFF)
Regards,
Méliné
------------------------------------------------------------------------------
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
------------------------------------------------------------------------------
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
