1 #include <cstdlib>
2 #include <iostream>
3 #include "MolData3Ddescriptors.h"
4 #include <GraphMol/RDKitBase.h>
5
6 #include "GraphMol/PartialCharges/GasteigerCharges.h"
7 #include "GraphMol/PartialCharges/GasteigerParams.h"
8
9 using namespace std;
10
MolData3Ddescriptors()11 MolData3Ddescriptors::MolData3Ddescriptors() {}
12
GetUn(int numAtoms)13 std::vector<double> MolData3Ddescriptors::GetUn(int numAtoms) {
14 std::vector<double> u(numAtoms, 1.0);
15
16 return u;
17 }
18
GetRelativeMW(const RDKit::ROMol & mol)19 std::vector<double> MolData3Ddescriptors::GetRelativeMW(
20 const RDKit::ROMol& mol) {
21 double* relativeMw = data3D.getMW();
22 int numAtoms = mol.getNumAtoms();
23
24 std::vector<double> pol(numAtoms, 0.0);
25 for (int i = 0; i < numAtoms; ++i) {
26 pol[i] = relativeMw[mol.getAtomWithIdx(i)->getAtomicNum() - 1];
27 }
28 return pol;
29 }
30
GetRelativePol(const RDKit::ROMol & mol)31 std::vector<double> MolData3Ddescriptors::GetRelativePol(
32 const RDKit::ROMol& mol) {
33 int numAtoms = mol.getNumAtoms();
34 double* relativePol = data3D.getPOL();
35
36 std::vector<double> pol(numAtoms, 0.0);
37 for (int i = 0; i < numAtoms; ++i) {
38 pol[i] = relativePol[mol.getAtomWithIdx(i)->getAtomicNum() - 1];
39 }
40 return pol;
41 }
42
GetRelativeVdW(const RDKit::ROMol & mol)43 std::vector<double> MolData3Ddescriptors::GetRelativeVdW(
44 const RDKit::ROMol& mol) {
45 int numAtoms = mol.getNumAtoms();
46 double* relativeVdW = data3D.getVDW();
47
48 std::vector<double> vdw(numAtoms, 0.0);
49 for (int i = 0; i < numAtoms; ++i) {
50 vdw[i] = relativeVdW[mol.getAtomWithIdx(i)->getAtomicNum() - 1];
51 }
52 return vdw;
53 }
54
GetRelativeRcov(const RDKit::ROMol & mol)55 std::vector<double> MolData3Ddescriptors::GetRelativeRcov(
56 const RDKit::ROMol& mol) {
57 int numAtoms = mol.getNumAtoms();
58 double* rcov = data3D.getRCOV();
59
60 std::vector<double> wroc(numAtoms, 0.0);
61 for (int i = 0; i < numAtoms; ++i) {
62 wroc[i] = rcov[mol.getAtomWithIdx(i)->getAtomicNum() - 1] / rcov[5];
63 }
64 return wroc;
65 }
66
GetRelativeENeg(const RDKit::ROMol & mol)67 std::vector<double> MolData3Ddescriptors::GetRelativeENeg(
68 const RDKit::ROMol& mol) {
69 int numAtoms = mol.getNumAtoms();
70 double* relativeNeg = data3D.getNEG();
71
72 std::vector<double> neg(numAtoms, 0.0);
73 for (int i = 0; i < numAtoms; ++i) {
74 neg[i] = relativeNeg[mol.getAtomWithIdx(i)->getAtomicNum() - 1];
75 }
76 return neg;
77 }
78
GetRelativeIonPol(const RDKit::ROMol & mol)79 std::vector<double> MolData3Ddescriptors::GetRelativeIonPol(
80 const RDKit::ROMol& mol) {
81 int numAtoms = mol.getNumAtoms();
82 double* absionpol = data3D.getIonPOL();
83
84 std::vector<double> ionpols(numAtoms, 0.0);
85 for (int i = 0; i < numAtoms; ++i) {
86 ionpols[i] = absionpol[mol.getAtomWithIdx(i)->getAtomicNum() - 1];
87 }
88 return ionpols;
89 }
90
GetCustomAtomProp(const RDKit::ROMol & mol,const std::string & customAtomPropName)91 std::vector<double> MolData3Ddescriptors::GetCustomAtomProp(
92 const RDKit::ROMol& mol, const std::string& customAtomPropName) {
93 int numAtoms = mol.getNumAtoms();
94
95 std::vector<double> customAtomArray(numAtoms, 1.0);
96 for (auto& atom : mol.atoms()) {
97 atom->getPropIfPresent(customAtomPropName, customAtomArray[atom->getIdx()]);
98 }
99 return customAtomArray;
100 }
101
GetCharges(const RDKit::ROMol & mol)102 std::vector<double> MolData3Ddescriptors::GetCharges(const RDKit::ROMol& mol) {
103 std::vector<double> charges(mol.getNumAtoms(), 0);
104 // use 12 iterations... can be more
105 RDKit::computeGasteigerCharges(mol, charges, 12, true);
106 return charges;
107 }
108
GetPrincipalQuantumNumber(int AtomicNum)109 int MolData3Ddescriptors::GetPrincipalQuantumNumber(int AtomicNum) {
110 if (AtomicNum <= 2) {
111 return 1;
112 } else if (AtomicNum <= 10) {
113 return 2;
114 } else if (AtomicNum <= 18) {
115 return 3;
116 } else if (AtomicNum <= 36) {
117 return 4;
118 } else if (AtomicNum <= 54) {
119 return 5;
120 } else if (AtomicNum <= 86) {
121 return 6;
122 } else {
123 return 7;
124 }
125 }
126
GetIState(const RDKit::ROMol & mol)127 std::vector<double> MolData3Ddescriptors::GetIState(const RDKit::ROMol& mol) {
128 int numAtoms = mol.getNumAtoms();
129 std::vector<double> Is(numAtoms, 1.0); // values set to 1 for Hs
130
131 for (int i = 0; i < numAtoms; ++i) {
132 const RDKit::Atom* atom = mol.getAtomWithIdx(i);
133 int atNum = atom->getAtomicNum();
134 int degree = atom->getDegree(); // number of substituants (heavy of not?)
135 if (degree > 0 && atNum > 1) {
136 int h = atom->getTotalNumHs(
137 true); // caution getTotalNumHs(true) to count h !!!!
138 int dv = RDKit::PeriodicTable::getTable()->getNouterElecs(atNum) -
139 h; // number of valence (explicit with Hs)
140 int N = GetPrincipalQuantumNumber(atNum); // principal quantum number
141 double d = (double)degree - h; // degree-h
142 if (d > 0) {
143 Is[i] = std::round(1000 * (4.0 / (N * N) * dv + 1.0) / d) / 1000;
144 }
145 }
146 }
147
148 return Is;
149 }
150
GetIStateDrag(const RDKit::ROMol & mol)151 std::vector<double> MolData3Ddescriptors::GetIStateDrag(
152 const RDKit::ROMol& mol) {
153 int numAtoms = mol.getNumAtoms();
154 std::vector<double> Is(numAtoms, 1.0);
155
156 for (int i = 0; i < numAtoms; ++i) {
157 const RDKit::Atom* atom = mol.getAtomWithIdx(i);
158 int atNum = atom->getAtomicNum();
159 int degree = atom->getDegree(); // number of substituants
160 if (degree > 0 && atNum > 1) {
161 int h = atom->getTotalNumHs(true);
162 int Zv = RDKit::PeriodicTable::getTable()->getNouterElecs(
163 atNum); // number of valence (explicit with Hs)
164 double dv = (double)Zv - h; // number of valence electron without Hs
165 int N = GetPrincipalQuantumNumber(atNum); // principal quantum number
166 double d = (double)degree - h; // degree-h
167 if (d > 0) {
168 Is[i] = std::round(1000 * (4.0 / (N * N) * dv + 1.0) / d) / 1000;
169 }
170 }
171 }
172
173 return Is;
174 }
175
176 // adaptation from EState.py
177 // we need the Is value only there
GetEState(const RDKit::ROMol & mol)178 std::vector<double> MolData3Ddescriptors::GetEState(const RDKit::ROMol& mol) {
179 int numAtoms = mol.getNumAtoms();
180
181 std::vector<double> Is = GetIState(mol);
182
183 double tmp, p;
184 double* dist = RDKit::MolOps::getDistanceMat(mol, false, false);
185 std::vector<double> accum(numAtoms, 0.0);
186
187 for (int i = 0; i < numAtoms; i++) {
188 for (int j = i + 1; j < numAtoms; j++) {
189 p = dist[i * numAtoms + j] + 1;
190 if (p < 1e6) {
191 tmp = (Is[i] - Is[j]) / (p * p);
192 accum[i] += tmp;
193 accum[j] -= tmp;
194 }
195 }
196 }
197
198 for (int i = 0; i < numAtoms; i++) {
199 Is[i] += accum[i];
200 }
201
202 return Is;
203 }
204
205 // modification of previous code to follow documentation from Padel code
GetEState2(const RDKit::ROMol & mol)206 std::vector<double> MolData3Ddescriptors::GetEState2(const RDKit::ROMol& mol) {
207 int numAtoms = mol.getNumAtoms();
208
209 std::vector<double> Si = GetIState(mol);
210
211 // in WHIM definition it's write:
212 double tmp, p, d;
213 double* dist = RDKit::MolOps::getDistanceMat(mol, false, false);
214 std::vector<double> accum(numAtoms, 0.0);
215
216 for (int i = 0; i < numAtoms; i++) {
217 for (int j = i + 1; j < numAtoms; j++) {
218 d = dist[i * numAtoms + j];
219 p = dist[i * numAtoms + j] + 1;
220 if (d == 1) {
221 tmp = (Si[i] - Si[j]) / (p * p);
222 accum[i] += tmp;
223 accum[j] -= tmp;
224 }
225 }
226 }
227
228 // add the Accum to the Si
229 // WHIM Si values
230 // electrotopological indices are scaled thus: Si'=Si + 7 => Si' > 0
231 // In this case, only the nonhydrogen atoms are considered,
232 // and the atomic electrotopological charge of each atom depends on its atom
233 // neighbor.
234 // So we should not use all the terms in the sum but only Adj matrix cases!
235
236 // Correct the Si adding the rescaling parameter for WHIM only
237 for (int i = 0; i < numAtoms; i++) {
238 Si[i] += accum[i] + 7.0;
239 }
240
241 return Si;
242 }
243