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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
|
// $Id$
//
// Copyright (C) 2005-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include "ShapeUtils.h"
#include "ShapeEncoder.h"
#include <Geometry/UniformGrid3D.h>
#include <GraphMol/RDKitBase.h>
#include <Geometry/Transform3D.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include <Geometry/GridUtils.h>
namespace RDKit {
namespace MolShapes {
void computeConfBox(const Conformer &conf, RDGeom::Point3D &leftBottom,
RDGeom::Point3D &rightTop, const RDGeom::Transform3D *trans,
double padding) {
double xmin, xmax, ymin, ymax, zmin, zmax;
xmin = ymin = zmin = 1.e8;
xmax = ymax = zmax = -1.e8;
unsigned int i, nAtms = conf.getNumAtoms();
for (i = 0; i < nAtms; ++i) {
RDGeom::Point3D loc = conf.getAtomPos(i);
if (trans) {
trans->TransformPoint(loc);
}
xmax = std::max(xmax, loc.x);
xmin = std::min(xmin, loc.x);
ymax = std::max(ymax, loc.y);
ymin = std::min(ymin, loc.y);
zmax = std::max(zmax, loc.z);
zmin = std::min(zmin, loc.z);
}
RDGeom::Point3D padPt(padding, padding, padding);
leftBottom.x = xmin;
leftBottom.y = ymin;
leftBottom.z = zmin;
rightTop.x = xmax;
rightTop.y = ymax;
rightTop.z = zmax;
leftBottom -= padPt;
rightTop += padPt;
}
void computeConfDimsAndOffset(const Conformer &conf, RDGeom::Point3D &dims,
RDGeom::Point3D &offSet,
const RDGeom::Transform3D *trans,
double padding) {
// RDGeom::Point3D lb, rb;
computeConfBox(conf, offSet, dims, trans, padding);
dims -= offSet;
}
std::vector<double> getConfDimensions(const Conformer &conf, double padding,
const RDGeom::Point3D *center,
bool ignoreHs) {
RDGeom::Point3D lb, rb;
computeConfBox(conf, lb, rb, nullptr, padding);
if (!center) {
RDGeom::Point3D cpt = MolTransforms::computeCentroid(conf, ignoreHs);
rb -= cpt;
lb -= cpt;
} else {
rb -= (*center);
lb -= (*center);
}
lb *= -1.0;
double dimX = 2.0 * std::max(rb.x, lb.x);
double dimY = 2.0 * std::max(rb.y, lb.y);
double dimZ = 2.0 * std::max(rb.z, lb.z);
std::vector<double> res;
res.reserve(3);
res.push_back(dimX);
res.push_back(dimY);
res.push_back(dimZ);
return res;
}
void computeUnionBox(const RDGeom::Point3D &leftBottom1,
const RDGeom::Point3D &rightTop1,
const RDGeom::Point3D &leftBottom2,
const RDGeom::Point3D &rightTop2,
RDGeom::Point3D &uLeftBottom, RDGeom::Point3D &uRightTop) {
uLeftBottom.x = std::min(leftBottom1.x, leftBottom2.x);
uLeftBottom.y = std::min(leftBottom1.y, leftBottom2.y);
uLeftBottom.z = std::min(leftBottom1.z, leftBottom2.z);
uRightTop.x = std::max(rightTop1.x, rightTop2.x);
uRightTop.y = std::max(rightTop1.y, rightTop2.y);
uRightTop.z = std::max(rightTop1.z, rightTop2.z);
}
double tverskyIndex(const ROMol &mol1, const ROMol &mol2, double alpha, double beta, int confId1,
int confId2, double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs) {
const Conformer &conf1 = mol1.getConformer(confId1);
const Conformer &conf2 = mol2.getConformer(confId2);
return tverskyIndex(conf1, conf2, alpha, beta, gridSpacing, bitsPerPoint,
vdwScale, stepSize, maxLayers, ignoreHs);
}
double tverskyIndex(const Conformer &conf1, const Conformer &conf2, double alpha, double beta,
double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs) {
RDGeom::Transform3D *trans = MolTransforms::computeCanonicalTransform(conf1);
// now use this transform and figure out what size grid we will need
// find the lower-left and upper-right corners for each of the conformers
// and take a union of these boxes - we will use this fo grid dimensions
RDGeom::Point3D leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop;
computeConfBox(conf1, leftBottom1, rightTop1, trans);
computeConfBox(conf2, leftBottom2, rightTop2, trans);
computeUnionBox(leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop);
// make the grid object to store the encoding
uRightTop -= uLeftBottom; // uRightTop now has grid dimensions
RDGeom::UniformGrid3D grd1(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
RDGeom::UniformGrid3D grd2(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
EncodeShape(conf1, grd1, trans, vdwScale, stepSize, maxLayers, ignoreHs);
EncodeShape(conf2, grd2, trans, vdwScale, stepSize, maxLayers, ignoreHs);
return RDGeom::tverskyIndex(grd1, grd2, alpha, beta);
}
double tanimotoDistance(const ROMol &mol1, const ROMol &mol2, int confId1,
int confId2, double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs) {
const Conformer &conf1 = mol1.getConformer(confId1);
const Conformer &conf2 = mol2.getConformer(confId2);
return tanimotoDistance(conf1, conf2, gridSpacing = 0.5, bitsPerPoint,
vdwScale, stepSize, maxLayers, ignoreHs);
}
double tanimotoDistance(const Conformer &conf1, const Conformer &conf2,
double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs) {
RDGeom::Transform3D *trans = MolTransforms::computeCanonicalTransform(conf1);
// now use this transform and figure out what size grid we will need
// find the lower-left and upper-right corners for each of the conformers
// and take a union of these boxes - we will use this fo grid dimensions
RDGeom::Point3D leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop;
computeConfBox(conf1, leftBottom1, rightTop1, trans);
computeConfBox(conf2, leftBottom2, rightTop2, trans);
computeUnionBox(leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop);
// make the grid object to store the encoding
uRightTop -= uLeftBottom; // uRightTop now has grid dimensions
RDGeom::UniformGrid3D grd1(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
RDGeom::UniformGrid3D grd2(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
EncodeShape(conf1, grd1, trans, vdwScale, stepSize, maxLayers, ignoreHs);
EncodeShape(conf2, grd2, trans, vdwScale, stepSize, maxLayers, ignoreHs);
return RDGeom::tanimotoDistance(grd1, grd2);
}
double protrudeDistance(const ROMol &mol1, const ROMol &mol2, int confId1,
int confId2, double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs, bool allowReordering) {
const Conformer &conf1 = mol1.getConformer(confId1);
const Conformer &conf2 = mol2.getConformer(confId2);
return protrudeDistance(conf1, conf2, gridSpacing = 0.5, bitsPerPoint,
vdwScale, stepSize, maxLayers, ignoreHs,
allowReordering);
}
double protrudeDistance(const Conformer &conf1, const Conformer &conf2,
double gridSpacing,
DiscreteValueVect::DiscreteValueType bitsPerPoint,
double vdwScale, double stepSize, int maxLayers,
bool ignoreHs, bool allowReordering) {
//
// FIX: all this duplicated code needs to be refactored out.
//
RDGeom::Transform3D *trans = MolTransforms::computeCanonicalTransform(conf1);
// now use this transform and figure out what size grid we will need
// find the lower-left and upper-right corners for each of the conformers
// and take a union of these boxes - we will use this fo grid dimensions
RDGeom::Point3D leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop;
computeConfBox(conf1, leftBottom1, rightTop1, trans);
computeConfBox(conf2, leftBottom2, rightTop2, trans);
computeUnionBox(leftBottom1, rightTop1, leftBottom2, rightTop2, uLeftBottom,
uRightTop);
// make the grid object to store the encoding
uRightTop -= uLeftBottom; // uRightTop now has grid dimensions
RDGeom::UniformGrid3D grd1(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
RDGeom::UniformGrid3D grd2(uRightTop.x, uRightTop.y, uRightTop.z, gridSpacing,
bitsPerPoint, &uLeftBottom);
EncodeShape(conf1, grd1, trans, vdwScale, stepSize, maxLayers, ignoreHs);
EncodeShape(conf2, grd2, trans, vdwScale, stepSize, maxLayers, ignoreHs);
double res;
if (allowReordering && (grd2.getOccupancyVect()->getTotalVal() <
grd1.getOccupancyVect()->getTotalVal())) {
res = RDGeom::protrudeDistance(grd2, grd1);
} else {
res = RDGeom::protrudeDistance(grd1, grd2);
}
return res;
}
}
}
|
