Re: [Rdkit-discuss] UFF and MMFF conformers energy

[Paolo Tosco]

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

Re: [Rdkit-discuss] UFF and MMFF conformers energy

[Peter S. Shenkin]

Small atomic displacements can cause large forcefield energy differences.
Computing molecular-mechanics energies from exactly the same coordinates
using two different force-fields is probably not a reasonable procedure.

It would be better to do an energy minimization with the two force fields
separately and live with the differing coordinates. Often the coordinate
differences will be small, but there are situations in which forcefield X
gives no conformation that is close to a minimized conformation obtained
from forcefield Y.

-P.

On Thu, Feb 9, 2017 at 11:48 AM, 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