A while ago I had a project which needed PMI descriptors (specifically NPR1
and NPR2) which were not available in the main branch of RDKit at the time.
At the time I used the fork by 'hahnda6' which provided the
calcPMIDescriptors() function, and this worked well. Now that PMI
descriptors are available in the main RDKit distrubution I thought I'd
rewrite my code to use the official version.
Building the new RDKit was no problem, but things went downhill shortly
after that. There's every chance that I've missed the relevant
documentation (I hope someone can point me in the right direction if so)
and done something stupid!
The issues are -
1) I can't find any documentation of the C++ API - the only reference to
PMI in the online RDKit documentation appears to be to the PMI.h file
2) Having written a program using the PMI and/or NPR functions, I
couldn't get it to compile until I added the -DRDK_BUILD_DESCRIPTORS3D
g++ -o sdf_pmi_blob sdf_pmi.cpp -I/packages/rdkit/include/rdkit
-L/packages/rdkit/lib -lDescriptors -lGraphMol -lFileParsers
-Wno-deprecated -O2 -DRDK_BUILD_DESCRIPTORS3D
This seems a bit odd...
3) Is it necessary to make separate calls to the individual PMI() and/or
NPR() functions? Surely this results in duplication of some of the heavier
calculations? I can't find any equivalent of calcPMIDescriptors() which
returned a 'Moments' struct containing all the PMI and NPR values in one go.
4) The big one! The returned results look very odd. They appear to relate
more to the dimensions of the molecule than the moments of inertia. For a
rod-like molecule (dimethylacetylene) I'd expect two large and one small
PMI (e.g. PMI1: 6.61651 PMI2: 150.434 PMI3: 150.434 NPR1: 0.0439828
NPR2: 0.999998) but actually get PMI1: 0.061647 PMI2: 0.061652 PMI3:
25.3699 NPR1: 0.002430 NPR2: 0.002430.
For disk-like (benzene) the result should be one large and two medium (e.g.
PMI1: 89.1448 PMI2: 89.1495 PMI3: 178.294 NPR1: 0.499987 NPR2:
0.500013) but get PMI1: 2.37457e-10 PMI2: 11.0844 PMI3: 11.0851 NPR1:
2.14213e-11 NPR2: 0.999933.
Finally for a roughly spherical molecule (neopentane) the NPR values look
reasonable (no great surprise) but the absolute PMI values may be too
small: old program - PMI1: 114.795 PMI2: 114.797 PMI3: 114.799
NPR1: 0.999966 NPR2: 0.999988, new program - PMI1: 6.59466 PMI2: 6.59488
PMI3: 6.59531 NPR1: 0.999902 NPR2: 0.999935
As I say, it's entirely likely that I'm doing something stupid here so any
pointers will be gratefully received. FWIW, the core of my program is -
mol = MolBlockToMol(ctab, true, false);
double pmi1 = RDKit::Descriptors::PMI1(*mol);
double pmi2 = RDKit::Descriptors::PMI2(*mol);
double pmi3 = RDKit::Descriptors::PMI3(*mol);
double npr1 = RDKit::Descriptors::NPR1(*mol);
double npr2 = RDKit::Descriptors::NPR2(*mol);
Thanks for any help!
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