1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
//
// Copyright (c) 2016, Guillaume GODIN
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above
// copyright notice, this list of conditions and the following
// disclaimer in the documentation and/or other materials provided
// with the distribution.
// * Neither the name of Institue of Cancer Research.
// nor the names of its contributors may be used to endorse or promote
// products derived from this software without specific prior written
// permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Adding RBF descriptors to 3D descriptors by Guillaume Godin
#include <GraphMol/RDKitBase.h>
#include "RDF.h"
#include "MolData3Ddescriptors.h"
#include <math.h>
namespace RDKit {
namespace Descriptors {
namespace {
MolData3Ddescriptors moldata3D;
std::vector<double> getG(int n) {
std::vector<double> res(n);
for (int i = 0; i < n; i++) {
res[i] = 1 + i * 0.5;
}
return res;
}
std::vector<double> prepareIState(const ROMol& mol) {
std::vector<double> IState = moldata3D.GetIStateDrag(
mol); // get the real IState value not the EState!
return IState;
}
void getRDFDesc(double* DM, const ROMol& mol, const Conformer& conf,
std::vector<double>& res) {
// std::vector<double> reserror(210, 0); need to avoid the res NaNs
// if (numAtoms < 4) return reserror;
// if (!conf.is3D()) return reserror;
int numAtoms = conf.getNumAtoms();
std::vector<double> R = getG(30);
std::vector<double> R1(30);
std::vector<double> R2(30);
std::vector<double> R3(30);
std::vector<double> R4(30);
std::vector<double> R5(30);
std::vector<double> R6(30);
std::vector<double> R7(30);
std::vector<double> Mass = moldata3D.GetRelativeMW(mol);
std::vector<double> RelativePol = moldata3D.GetRelativePol(mol);
std::vector<double> IonPol = moldata3D.GetRelativeIonPol(mol);
std::vector<double> RelativeElectroNeg = moldata3D.GetRelativeENeg(mol);
std::vector<double> RelativeVdW = moldata3D.GetRelativeVdW(mol);
std::vector<double> IState = prepareIState(mol);
double p;
for (size_t i = 0; i < R.size(); i++) {
double res1 = 0.0;
double res2 = 0.0;
double res3 = 0.0;
double res4 = 0.0;
double res5 = 0.0;
double res6 = 0.0;
double res7 = 0.0;
for (int j = 0; j < numAtoms - 1; j++) {
for (int k = j + 1; k < numAtoms; k++) {
p = exp(-100 * pow(R[i] - DM[j * numAtoms + k], 2));
res1 += p; // "u"
res2 += Mass[j] * Mass[k] * p; // "m"
res3 += RelativeVdW[j] * RelativeVdW[k] * p; // "v"
res4 += RelativeElectroNeg[j] * RelativeElectroNeg[k] * p; //"e"
res5 += RelativePol[j] * RelativePol[k] * p; // "p"
res6 += IonPol[j] * IonPol[k] * p; // "i"
res7 += IState[j] * IState[k] * p; // "s"
}
}
R1[i] = round(1000 * res1) / 1000;
R2[i] = round(1000 * res2) / 1000;
R3[i] = round(1000 * res3) / 1000;
R4[i] = round(1000 * res4) / 1000;
R5[i] = round(1000 * res5) / 1000;
R6[i] = round(1000 * res6) / 1000;
R7[i] = round(1000 * res7) / 1000;
}
R1.insert(R1.end(), R2.begin(), R2.end());
R1.insert(R1.end(), R3.begin(), R3.end());
R1.insert(R1.end(), R4.begin(), R4.end());
R1.insert(R1.end(), R5.begin(), R5.end());
R1.insert(R1.end(), R6.begin(), R6.end());
R1.insert(R1.end(), R7.begin(), R7.end());
res = R1;
}
void getRDFDescCustom(double* DM, const ROMol& mol, const Conformer& conf,
std::vector<double>& res,
const std::string& customAtomPropName) {
int numAtoms = conf.getNumAtoms();
std::vector<double> R = getG(30);
std::vector<double> R1(30);
std::vector<double> customAtomArray =
moldata3D.GetCustomAtomProp(mol, customAtomPropName);
double p;
for (size_t i = 0; i < R.size(); i++) {
double res = 0.0;
for (int j = 0; j < numAtoms - 1; j++) {
for (int k = j + 1; k < numAtoms; k++) {
p = exp(-100 * pow(R[i] - DM[j * numAtoms + k], 2));
res += customAtomArray[j] * customAtomArray[k] * p; // "custom"
}
}
R1[i] = round(1000 * res) / 1000;
}
res = R1;
}
void GetRDF(double* dist3D, const ROMol& mol, const Conformer& conf,
std::vector<double>& res) {
getRDFDesc(dist3D, mol, conf, res);
}
void GetRDFone(double* dist3D, const ROMol& mol, const Conformer& conf,
std::vector<double>& res, const std::string customAtomPropName) {
getRDFDescCustom(dist3D, mol, conf, res, customAtomPropName);
}
} // end of anonymous namespace
void RDF(const ROMol& mol, std::vector<double>& res, int confId,
const std::string& customAtomPropName) {
// RDF010u RDF015u RDF020u RDF025u RDF030u RDF035u RDF040u RDF045u RDF050u
// RDF055u RDF060u RDF065u RDF070u RDF075u RDF080u RDF085u RDF090u RDF095u
// RDF100u RDF105u RDF110u RDF115u RDF120u RDF125u RDF130u RDF135u RDF140u
// RDF145u RDF150u RDF155u
// RDF010m RDF015m RDF020m RDF025m RDF030m RDF035m RDF040m RDF045m RDF050m
// RDF055m RDF060m RDF065m RDF070m RDF075m RDF080m RDF085m RDF090m RDF095m
// RDF100m RDF105m RDF110m RDF115m RDF120m RDF125m RDF130m RDF135m RDF140m
// RDF145m RDF150m RDF155m
// RDF010v RDF015v RDF020v RDF025v RDF030v RDF035v RDF040v RDF045v RDF050v
// RDF055v RDF060v RDF065v RDF070v RDF075v RDF080v RDF085v RDF090v RDF095v
// RDF100v RDF105v RDF110v RDF115v RDF120v RDF125v RDF130v RDF135v RDF140v
// RDF145v RDF150v RDF155v
// RDF010e RDF015e RDF020e RDF025e RDF030e RDF035e RDF040e RDF045e RDF050e
// RDF055e RDF060e RDF065e RDF070e RDF075e RDF080e RDF085e RDF090e RDF095e
// RDF100e RDF105e RDF110e RDF115e RDF120e RDF125e RDF130e RDF135e RDF140e
// RDF145e RDF150e RDF155e
// RDF010p RDF015p RDF020p RDF025p RDF030p RDF035p RDF040p RDF045p RDF050p
// RDF055p RDF060p RDF065p RDF070p RDF075p RDF080p RDF085p RDF090p RDF095p
// RDF100p RDF105p RDF110p RDF115p RDF120p RDF125p RDF130p RDF135p RDF140p
// RDF145p RDF150p RDF155p
// RDF010i RDF015i RDF020i RDF025i RDF030i RDF035i RDF040i RDF045i RDF050i
// RDF055i RDF060i RDF065i RDF070i RDF075i RDF080i RDF085i RDF090i RDF095i
// RDF100i RDF105i RDF110i RDF115i RDF120i RDF125i RDF130i RDF135i RDF140i
// RDF145i RDF150i RDF155i
// RDF010s RDF015s RDF020s RDF025s RDF030s RDF035s RDF040s RDF045s RDF050s
// RDF055s RDF060s RDF065s RDF070s RDF075s RDF080s RDF085s RDF090s RDF095s
// RDF100s RDF105s RDF110s RDF115s RDF120s RDF125s RDF130s RDF135s RDF140s
// RDF145s RDF150s RDF155s
PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
// int numAtoms = mol.getNumAtoms();
// if (numAtoms < 4) return reserror;
const Conformer& conf = mol.getConformer(confId);
// if (!conf.is3D()) return reserror;
double* dist3D =
MolOps::get3DDistanceMat(mol, confId, false, true); // 3D distance matrix
if (customAtomPropName != "") {
// std::cout << " Using CustomAtomPropertie\n";
// do something
res.clear();
res.resize(30); // 7 * 30
GetRDFone(dist3D, mol, conf, res, customAtomPropName);
} else {
res.clear();
res.resize(210); // 7 * 30
GetRDF(dist3D, mol, conf, res);
}
}
} // namespace Descriptors
} // namespace RDKit
|
