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
|
// $Id: ShapeUtils.cpp 1528 2010-09-26 17:04:37Z glandrum $
//
// 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, 0, 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 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;
}
}
}
|
