[Rdkit-discuss] Calculating strain energy of conformers

[Lewis M. <lew...@gm...>]

Hi RDKit,
I'm exploring strain energies in the context of virtual screening,
something that has been considered for a while[1] and is still being
explored today[2].

There may not be a canonical way, but is this a valid/good way to calculate
strain energy? I'm just not sure if I'm using the MMFF correctly.

from rdkit import Chem
from rdkit.Chem import AllChem
testmol =
Chem.MolFromSmiles('Cc1c(c2cc(ccc2n1C(=O)c3ccc(cc3)Cl)OC)CC(=O)O')
#indomethacin
testmol = Chem.AddHs(testmol)
AllChem.EmbedMultipleConfs(testmol, 50)

mp = AllChem.MMFFGetMoleculeProperties(testmol, mmffVariant='MMFF94s')

# keep angle-related terms
mp.SetMMFFOopTerm(True)
mp.SetMMFFAngleTerm(True)
mp.SetMMFFTorsionTerm(True)
# turn off anything unrelated to angles
mp.SetMMFFStretchBendTerm(False)
mp.SetMMFFBondTerm(False)
mp.SetMMFFVdWTerm(False)
mp.SetMMFFEleTerm(False)

for c in range(50):
    print(AllChem.MMFFGetMoleculeForceField(testmol, mp,
confId=c).CalcEnergy())
    MMFFOptimizeMolecule(testmol, confId=c, mmffVariant='MMFF94s')
    print('\t'+str(AllChem.MMFFGetMoleculeForceField(testmol, mp,
confId=c).CalcEnergy()))

Thanks for any comments or advice!
Lewis


[1]Conformational Analysis of Drug-Like Molecules Bound to Proteins: An
Extensive Study of Ligand Reorganization upon Binding, 2003
[2]Ligand Strain Energy in Large Library Docking, 2021