1 //
2 // Copyright (C) 2016 Greg Landrum
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 #include <GraphMol/RDKitBase.h>
11 #include <GraphMol/MolTransforms/MolTransforms.h>
12 #include <Geometry/point.h>
13
14 #include "PMI.h"
15
16 #include <Eigen/Dense>
17
18 namespace RDKit {
19 namespace Descriptors {
20 namespace {
21
getMoments(const ROMol & mol,int confId,bool useAtomicMasses,double & pm1,double & pm2,double & pm3,bool force)22 bool getMoments(const ROMol& mol, int confId, bool useAtomicMasses, double& pm1,
23 double& pm2, double& pm3, bool force) {
24 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
25 const char* pn1 = useAtomicMasses ? "_PMI1_mass" : "_PMI1";
26 const char* pn2 = useAtomicMasses ? "_PMI2_mass" : "_PMI2";
27 const char* pn3 = useAtomicMasses ? "_PMI3_mass" : "_PMI3";
28
29 if (!force && mol.hasProp(pn1) && mol.hasProp(pn2) && mol.hasProp(pn3)) {
30 mol.getProp(pn1, pm1);
31 mol.getProp(pn2, pm2);
32 mol.getProp(pn3, pm3);
33 return true;
34 }
35
36 const Conformer& conf = mol.getConformer(confId);
37
38 Eigen::Matrix3d axes;
39 Eigen::Vector3d moments;
40 bool res;
41 bool ignoreHs = false;
42 if (useAtomicMasses) {
43 std::vector<double> weights;
44 weights.resize(mol.getNumAtoms());
45 for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
46 ++cai) {
47 weights[(*cai)->getIdx()] = (*cai)->getMass();
48 }
49 res = MolTransforms::computePrincipalAxesAndMoments(
50 conf, axes, moments, ignoreHs, force, &weights);
51 } else {
52 res = MolTransforms::computePrincipalAxesAndMoments(conf, axes, moments,
53 ignoreHs, force);
54 }
55 if (res) {
56 pm1 = moments(0);
57 pm2 = moments(1);
58 pm3 = moments(2);
59 mol.setProp(pn1, pm1, true);
60 mol.setProp(pn2, pm2, true);
61 mol.setProp(pn3, pm3, true);
62 }
63 return res;
64 }
getMomentsFromGyration(const ROMol & mol,int confId,bool useAtomicMasses,double & pm1,double & pm2,double & pm3,bool force)65 bool getMomentsFromGyration(const ROMol& mol, int confId, bool useAtomicMasses,
66 double& pm1, double& pm2, double& pm3, bool force) {
67 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
68 const char* pn1 = useAtomicMasses ? "_PMI1_mass_cov" : "_PMI1_cov";
69 const char* pn2 = useAtomicMasses ? "_PMI2_mass_cov" : "_PMI2_cov";
70 const char* pn3 = useAtomicMasses ? "_PMI3_mass_cov" : "_PMI3_cov";
71
72 if (!force && mol.hasProp(pn1) && mol.hasProp(pn2) && mol.hasProp(pn3)) {
73 mol.getProp(pn1, pm1);
74 mol.getProp(pn2, pm2);
75 mol.getProp(pn3, pm3);
76 return true;
77 }
78
79 const Conformer& conf = mol.getConformer(confId);
80
81 Eigen::Matrix3d axes;
82 Eigen::Vector3d moments;
83 bool res;
84 bool ignoreHs = false;
85 if (useAtomicMasses) {
86 std::vector<double> weights;
87 weights.resize(mol.getNumAtoms());
88 for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
89 ++cai) {
90 weights[(*cai)->getIdx()] = (*cai)->getMass();
91 }
92 res = MolTransforms::computePrincipalAxesAndMomentsFromGyrationMatrix(
93 conf, axes, moments, ignoreHs, force, &weights);
94 } else {
95 res = MolTransforms::computePrincipalAxesAndMomentsFromGyrationMatrix(
96 conf, axes, moments, ignoreHs, force);
97 }
98 if (res) {
99 pm1 = moments(0);
100 pm2 = moments(1);
101 pm3 = moments(2);
102 mol.setProp(pn1, pm1, true);
103 mol.setProp(pn2, pm2, true);
104 mol.setProp(pn3, pm3, true);
105 }
106 return res;
107 }
108
109 } // end of anonymous namespace
110
NPR1(const ROMol & mol,int confId,bool useAtomicMasses,bool force)111 double NPR1(const ROMol& mol, int confId, bool useAtomicMasses, bool force) {
112 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
113 double pm1, pm2, pm3;
114 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
115 // the eigenvector calculation failed
116 return 0.0; // FIX: throw an exception here?
117 }
118 if (pm3 < 1e-8) {
119 return 0.0;
120 }
121 return pm1 / pm3;
122 }
NPR2(const ROMol & mol,int confId,bool useAtomicMasses,bool force)123 double NPR2(const ROMol& mol, int confId, bool useAtomicMasses, bool force) {
124 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
125 double pm1, pm2, pm3;
126 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
127 // the eigenvector calculation failed
128 return 0.0; // FIX: throw an exception here?
129 }
130 if (pm3 < 1e-8) {
131 return 0.0;
132 }
133 return pm2 / pm3;
134 }
PMI1(const ROMol & mol,int confId,bool useAtomicMasses,bool force)135 double PMI1(const ROMol& mol, int confId, bool useAtomicMasses, bool force) {
136 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
137 double pm1, pm2, pm3;
138 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
139 // the eigenvector calculation failed
140 return 0.0; // FIX: throw an exception here?
141 }
142 return pm1;
143 }
PMI2(const ROMol & mol,int confId,bool useAtomicMasses,bool force)144 double PMI2(const ROMol& mol, int confId, bool useAtomicMasses, bool force) {
145 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
146 double pm1, pm2, pm3;
147 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
148 // the eigenvector calculation failed
149 return 0.0; // FIX: throw an exception here?
150 }
151 return pm2;
152 }
PMI3(const ROMol & mol,int confId,bool useAtomicMasses,bool force)153 double PMI3(const ROMol& mol, int confId, bool useAtomicMasses, bool force) {
154 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
155 double pm1, pm2, pm3;
156 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
157 // the eigenvector calculation failed
158 return 0.0; // FIX: throw an exception here?
159 }
160 return pm3;
161 }
162
radiusOfGyration(const ROMol & mol,int confId,bool useAtomicMasses,bool force)163 double radiusOfGyration(const ROMol& mol, int confId, bool useAtomicMasses,
164 bool force) {
165 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
166 double pm1, pm2, pm3;
167 if (!getMomentsFromGyration(mol, confId, useAtomicMasses, pm1, pm2, pm3,
168 force)) {
169 // the eigenvector calculation failed
170 return 0.0; // FIX: throw an exception here?
171 }
172 return sqrt(pm1 + pm2 + pm3);
173 }
174
inertialShapeFactor(const ROMol & mol,int confId,bool useAtomicMasses,bool force)175 double inertialShapeFactor(const ROMol& mol, int confId, bool useAtomicMasses,
176 bool force) {
177 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
178 double pm1, pm2, pm3;
179 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
180 // the eigenvector calculation failed
181 return 0.0; // FIX: throw an exception here?
182 }
183 if (pm1 < 1e-4 || pm3 < 1e-4) {
184 // planar or no coordinates
185 return 0.0;
186 } else {
187 return pm2 / (pm1 * pm3);
188 }
189 }
eccentricity(const ROMol & mol,int confId,bool useAtomicMasses,bool force)190 double eccentricity(const ROMol& mol, int confId, bool useAtomicMasses,
191 bool force) {
192 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
193 double pm1, pm2, pm3;
194 if (!getMoments(mol, confId, useAtomicMasses, pm1, pm2, pm3, force)) {
195 // the eigenvector calculation failed
196 return 0.0; // FIX: throw an exception here?
197 }
198 if (pm3 < 1e-4 || (pm3 * pm3 - pm1 * pm1) < 1e-4) {
199 // no coordinates or very close to degeneracy
200 return 0.0;
201 } else {
202 return sqrt(pm3 * pm3 - pm1 * pm1) / pm3;
203 }
204 }
205
asphericity(const ROMol & mol,int confId,bool useAtomicMasses,bool force)206 double asphericity(const ROMol& mol, int confId, bool useAtomicMasses,
207 bool force) {
208 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
209 double pm1, pm2, pm3;
210 if (!getMomentsFromGyration(mol, confId, useAtomicMasses, pm1, pm2, pm3,
211 force)) {
212 // the eigenvector calculation failed
213 return 0.0; // FIX: throw an exception here?
214 }
215 if (pm3 < 1e-4) {
216 // no coordinates
217 return 0.0;
218 } else {
219 double denom = pm1 + pm2 + pm3;
220
221 return 0.5 * (pow(pm1 - pm2, 2) + pow(pm1 - pm3, 2) + pow(pm2 - pm3, 2)) /
222 (denom * denom);
223 }
224 }
spherocityIndex(const ROMol & mol,int confId,bool force)225 double spherocityIndex(const ROMol& mol, int confId, bool force) {
226 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers");
227 bool useAtomicMasses = false;
228 double pm1, pm2, pm3;
229 if (!getMomentsFromGyration(mol, confId, useAtomicMasses, pm1, pm2, pm3,
230 force)) {
231 // the eigenvector calculation failed
232 return 0.0; // FIX: throw an exception here?
233 }
234 if (pm3 < 1e-4) {
235 // no coordinates
236 return 0.0;
237 } else {
238 return 3. * pm1 / (pm1 + pm2 + pm3);
239 }
240 }
241
242 } // namespace Descriptors
243 } // namespace RDKit
244