Dear Joseph,
please find attached a commented Python script which illustrates how to
achieve what you describe. You can do the same from C++. In the
C++/Python API documentation you may also find that accessor functions
are available to retrieve force constants/equilibrium values for the
various energy terms.
Kind regards,
Paolo
On 05/07/15 21:22, Joseph Katigbak wrote:
Hi Guys,
I'm fairly new at using rdkit so everything is still new to me. I
would like to get the individual energy terms from the MMFF force
field, ie the energy for the stretch, bending, torsion, electrostatic,
etc.... Any help would be greatly appreciated and thank you so much
for this amazing tool
-Joseph
------------------------------------------------------------------------------
One dashboard for servers and applications across Physical-Virtual-Cloud
Widest out-of-the-box monitoring support with 50+ applications
Performance metrics, stats and reports that give you Actionable Insights
Deep dive visibility with transaction tracing using APM Insight.
http://ad.doubleclick.net/ddm/clk/290420510;117567292;y
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
#!/usr/bin/env python
import rdkit
from rdkit import Chem
from rdkit.Chem import rdmolops
from rdkit.Chem import ChemicalForceFields
molBlock = """\
bilastine
RDKit 3D
34 37 0 0 0 0 0 0 0 0999 V2000
9.6369 -1.5816 -0.4538 O 0 0 0 0 0 0 0 0 0 0 0 0
9.2303 -0.3807 -0.2633 C 0 0 0 0 0 0 0 0 0 0 0 0
9.5332 0.6344 -0.9854 O 0 0 0 0 0 0 0 0 0 0 0 0
8.2967 -0.1411 0.9097 C 0 0 0 0 0 0 0 0 0 0 0 0
6.8663 -0.2160 0.3704 C 0 0 0 0 0 0 0 0 0 0 0 0
5.9989 0.8500 0.4803 C 0 0 0 0 0 0 0 0 0 0 0 0
4.7201 0.7725 -0.0198 C 0 0 0 0 0 0 0 0 0 0 0 0
4.2834 -0.3754 -0.6406 C 0 0 0 0 0 0 0 0 0 0 0 0
5.1405 -1.4474 -0.7519 C 0 0 0 0 0 0 0 0 0 0 0 0
6.4200 -1.3648 -0.2510 C 0 0 0 0 0 0 0 0 0 0 0 0
2.8696 -0.4259 -1.1981 C 0 0 0 0 0 0 0 0 0 0 0 0
1.8257 -1.1216 -0.2911 C 0 0 0 0 0 0 0 0 0 0 0 0
0.4579 -0.4915 -0.4046 N 0 0 0 0 0 0 0 0 0 0 0 0
0.3051 0.7102 0.4947 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.0550 1.3967 0.2648 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.2149 0.4047 0.4938 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.5516 1.1154 0.3398 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.6597 2.3384 -0.1560 N 0 0 0 0 0 0 0 0 0 0 0 0
-4.9822 2.7066 -0.1222 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.6589 3.8678 -0.4994 C 0 0 0 0 0 0 0 0 0 0 0 0
-7.0104 3.8882 -0.3066 C 0 0 0 0 0 0 0 0 0 0 0 0
-7.6867 2.7902 0.2442 C 0 0 0 0 0 0 0 0 0 0 0 0
-7.0204 1.6533 0.6115 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.6449 1.6396 0.4122 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.7772 0.6249 0.7039 N 0 0 0 0 0 0 0 0 0 0 0 0
-5.0969 -0.6936 1.2805 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.7544 -1.6506 0.2687 C 0 0 0 0 0 0 0 0 0 0 0 0
-7.1570 -1.3508 0.1409 O 0 0 0 0 0 0 0 0 0 0 0 0
-7.8157 -2.1358 -0.8719 C 0 0 0 0 0 0 0 0 0 0 0 0
-8.0664 -3.5888 -0.4305 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.0347 -0.8020 -0.4527 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.6674 -1.4778 -0.2106 C 0 0 0 0 0 0 0 0 0 0 0 0
8.6674 1.2130 1.5547 C 0 0 0 0 0 0 0 0 0 0 0 0
8.5347 -1.1968 2.0090 C 0 0 0 0 0 0 0 0 0 0 0 0
2 1 1 0
3 2 2 0
4 2 1 0
5 4 1 0
6 5 2 0
7 6 1 0
8 7 2 0
9 8 1 0
10 9 2 0
5 10 1 0
11 8 1 0
12 11 1 0
13 12 1 0
14 13 1 0
15 14 1 0
16 15 1 0
17 16 1 0
18 17 2 0
19 18 1 0
20 19 2 0
21 20 1 0
22 21 2 0
23 22 1 0
24 23 2 0
19 24 1 0
25 24 1 0
17 25 1 0
26 25 1 0
27 26 1 0
28 27 1 0
29 28 1 0
30 29 1 0
31 16 1 0
32 31 1 0
13 32 1 0
33 4 1 0
34 4 1 0
M CHG 2 1 -1 13 1
M END
$$$$
"""
mol = Chem.MolFromMolBlock(molBlock)
mol = rdmolops.AddHs(mol, addCoords = True)
mp = ChemicalForceFields.MMFFGetMoleculeProperties(mol)
ff = ChemicalForceFields.MMFFGetMoleculeForceField(mol, mp)
# first we minimize the molecule with all terms included
ff.Minimize()
# setting verbosity to 2 (high) allows printing out all individual
# energy contributions, split by energy term
mp = ChemicalForceFields.MMFFGetMoleculeProperties(mol, mmffVerbosity = 2)
ff = ChemicalForceFields.MMFFGetMoleculeForceField(mol, mp)
# using the various setter methods on MMFFGetMoleculeProperties
# you can include/exclude each individual energy term; the force-field
# equation will be computed using only the energy terms which
# have been included BEFORE creating the force-field
# the following example includes one energy term at a time
# and computes the energy according that term only
mp = ChemicalForceFields.MMFFGetMoleculeProperties(mol)
print '\n{0:>16s} energy: {1:12.4f} kcal/mol'.format('Total',
ff.CalcEnergy())
for i in range(7):
termList = [['BondStretch', False], ['AngleBend', False],
['StretchBend', False], ['OopBend', False], ['Torsion', False],
['VdW', False], ['Electrostatic', False]]
termList[i][1] = True
mp.SetMMFFBondTerm(termList[0][1])
mp.SetMMFFAngleTerm(termList[1][1])
mp.SetMMFFStretchBendTerm(termList[2][1])
mp.SetMMFFOopTerm(termList[3][1])
mp.SetMMFFTorsionTerm(termList[4][1])
mp.SetMMFFVdWTerm(termList[5][1])
mp.SetMMFFEleTerm(termList[6][1])
ff = ChemicalForceFields.MMFFGetMoleculeForceField(mol, mp)
print '{0:>16s} energy: {1:12.4f} kcal/mol'.format(termList[i][0],
ff.CalcEnergy())
------------------------------------------------------------------------------
One dashboard for servers and applications across Physical-Virtual-Cloud
Widest out-of-the-box monitoring support with 50+ applications
Performance metrics, stats and reports that give you Actionable Insights
Deep dive visibility with transaction tracing using APM Insight.
http://ad.doubleclick.net/ddm/clk/290420510;117567292;y
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
