Hi Nitzan,
you should try to modify the molecule conformation as little as possible
between the two calculations; see the example below.
rom rdkit import Chem
from rdkit.Chem import AllChem
test_ion1 = Chem.MolFromSmiles('COC(=O)C(c1ccccc1)C1CCCC[NH2+]1')
test_ion1 = Chem.AddHs(test_ion1)
AllChem.EmbedMolecule(test_ion1, randomSeed=2)
prop1 = AllChem.MMFFGetMoleculeProperties(test_ion1,
mmffVariant="MMFF94")
ff1 = AllChem.MMFFGetMoleculeForceField(test_ion1, prop1)
ff1.Minimize(maxIts=1000)
print(ff1.CalcEnergy())
50.90705167464979
mn = test_ion1.GetSubstructMatch(Chem.MolFromSmarts("[NH2+]([H])[H]"))
n = mn[0]
na = test_ion1.GetAtomWithIdx(n)
mo = test_ion1.GetSubstructMatch(Chem.MolFromSmarts("O=CO"))
o = mo[0]
for h in mn[1:]:
test_ion2 = Chem.RWMol(test_ion1)
na.SetFormalCharge(0)
test_ion2.RemoveBond(n, h)
test_ion2.RemoveAtom(h)
oa = test_ion2.GetAtomWithIdx(o)
oa.SetFormalCharge(1)
oa.SetNoImplicit(True)
oa.SetNumExplicitHs(1)
Chem.SanitizeMol(test_ion2)
test_ion2 = Chem.AddHs(test_ion2, addCoords=True)
prop2 = AllChem.MMFFGetMoleculeProperties(test_ion2,
mmffVariant="MMFF94")
ff2 = AllChem.MMFFGetMoleculeForceField(test_ion2, prop2)
print(ff2.CalcEnergy())
53.20633094793681
58.41577310354003
In this example I built two conformers with axial and equatorial N-H,
and both their MMFF94 energies are larger than the NH2+ tautomer, as
expected.
To make the calculation more accurate one should actually take into
consideration multiple conformers and orientations of the mobile
protons. Also, molecular mechanics are most likely not accurate enough
for this kind of exercise - you should look into using a QM method to
get more accurate figures.
Cheers,
p.
On 08/05/19 19:26, Nitzan Tzanani wrote:
Dear rdkitters,
I am trying to use rdkit descriptors to investigate protonation of
molecules in the gas phase
while testing my procedure on Ritalin, I try to calculate the
protonation energy on different functional groups in the molecule
although my understanding tells me protonation should occur on the
secondary amine, energy calculations point at the keto group as the
favorable site.
Bellow is my code
Is it wrong to use this approach? Is it a problem of using the wrong
parameters?
Any help would be greatly appreciated
Yours
Nitzan
test_ion1 = Chem.MolFromSmiles('COC(=O)C(c1ccccc1)C1CCCC[NH2+]1')
test_ion1 = Chem.Add Hs(test_ion1)
AllChem.Embed Molecule( test_ion1,randomSeed=2)
AllChem.MMFFOptimizeMolecule(test_ion1,maxlters-1000,mmffVariant="MMFF94")
Prop1 = AllChem.MMFFGetMoleculeProperties(test_ion1,mmffVariant="MMFF94")
ff1 = AllChem.MMFFGetMoleculeForceField (test_ion1,Prop1)
print ff1.CalcEnergy()
test_ion2 = Chem.MolFromSmiles('COC(=[OH+])C(c1ccccc1)C1CCCCN1')
test_ion2 = Chem.Add Hs(test_ion2)
AllChem.Embed Molecule(test_ion2,randomSeed=2)
AllChem.MMFFOptimizeMolecule( test_ion2,maxlters-1000,mmffVariant="MMFF94")
Prop2 = AllChem.MMFFGetMoleculeProperties(test_ion2,mmffVariant="MMFF94")
ff2 = AllChem.MMFFGetMoleculeForceField (test_ion2,Prop2)
print ff2.CalcEnergy()
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail&utm_term=icon>
Virus-free. www.avast.com
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail&utm_term=link>
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
