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

Reply via email to