Dear James,
>>> for atom in mh.GetAtoms():
>>> if not atom.GetAtomicNum() == 1:
>>> idx = atom.GetIdx()
>>> ff.MMFFAddPositionConstraint(idx, maxDispl=0.5,
forceConstant=100)
>>> ff.Minimize(maxIts=1)
>>> ff.CalcEnergy()
This will give you the energy *including* the constraint term, which is
probably not what you want.
If you want the clean MMFF energy term with no contraints, I'd rebuild
the force field after the constrained minimization without the
constraint and then invoke CalcEnergy(). This will give you the energy
of the partially relaxed system without any constraint terms.
If you wish to break down the energy contributions to the total energy,
currently the procedure is really inefficient from Python, as it
involves a separate force-field rebuild and energy calculation for each
MMFF energy term. I hope this is not a problem given the small number of
molecules involved; I'll look into making this more efficient.
Please find below a Python snippet that will do what you need on molecule m:
from rdkit import Chem
from rdkit.Chem import rdForceFieldHelpers
from rdkit.Chem import ChemicalForceFields
def energyTermBreakdown(m):
mp = ChemicalForceFields.MMFFGetMoleculeProperties(m)
eTotal = 0.0
ffTerms = ('Bond', 'Angle', 'StretchBend', 'Torsion', 'Oop', 'VdW',
'Ele')
for iTerm in ffTerms:
for jTerm in ffTerms:
state = (iTerm == jTerm)
setMethod = getattr(mp, 'SetMMFF' + jTerm + 'Term')
setMethod(state)
ff = rdForceFieldHelpers.MMFFGetMoleculeForceField(m, mp)
e = ff.CalcEnergy()
print ('{0:s} energy: {1:.4f}'.format(iTerm, e))
eTotal += e
print ('Total energy: {0:.4f}'.format(eTotal))
Cheers,
p.
On 03/22/18 08:27, James Davidson wrote:
Dear All,
Recently I have been assessing some ligand conformations from crystal
structures to identify any non-ideal bond lengths, angles, torsions,
or non-bonded contacts.
What I am doing at the moment is adding some positional constraints to
the crystallographic heavy atom positions, and calculating the energy
before and after minimisation:
>>> m = Chem.MolFromMolFile(‘input.mol’)
>>> mh = AllChem.AddHs(m, addCoords=True)
>>> mp = AllChem.MMFFGetMoleculeProperties(mh, mmffVariant='MMFF94s')
>>> ff = AllChem.MMFFGetMoleculeForceField(mh, mp)
>>> ff.CalcEnergy()
This gives the ‘raw’ energy.
>>> for atom in mh.GetAtoms():
>>> if not atom.GetAtomicNum() == 1:
>>> idx = atom.GetIdx()
>>> ff.MMFFAddPositionConstraint(idx, maxDispl=0.5,
forceConstant=100)
>>> ff.Minimize(maxIts=1)
>>> ff.CalcEnergy()
And this gives the energy after applying a moderate restraint (100
kcal/mol, with a maximum displacement of 0.5 A).
So I think this is ok(?), and I can compare the two energies and
inspect the conformations visually.
What I was wondering was whether there is a way of obtaining the
individual energy terms (ie each bonded and non-bonded term, angle,
and torsion)?
Because what I’d really like to do is identify the areas of the
molecule that contribute the most to the pre- and post- minimisation
energy difference.
Any suggestions would be greatly appreciated!
Kind regards
James
__
PLEASE READ: This email is confidential and may be privileged. It is
intended for the named addressee(s) only and access to it by anyone
else is unauthorised. If you are not an addressee, any disclosure or
copying of the contents of this email or any action taken (or not
taken) in reliance on it is unauthorised and may be unlawful. If you
have received this email in error, please notify the sender or
postmas...@vernalis.com. Email is not a secure method of communication
and the Company cannot accept responsibility for the accuracy or
completeness of this message or any attachment(s). Please check this
email for virus infection for which the Company accepts no
responsibility. If verification of this email is sought then please
request a hard copy. Unless otherwise stated, any views or opinions
presented are solely those of the author and do not represent those of
the Company.
The Vernalis Group of Companies
100 Berkshire Place
Wharfedale Road
Winnersh, Berkshire
RG41 5RD, England
Tel: +44 (0)118 938
To access trading company registration and address details, please go
to the Vernalis website at www.vernalis.com and click on the "Company
address and registration details" link at the bottom of the page..
__
--
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
