1 //
2 // Copyright (c) 2012, Institue of Cancer Research.
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are
7 // met:
8 //
9 // * Redistributions of source code must retain the above copyright
10 // notice, this list of conditions and the following disclaimer.
11 // * Redistributions in binary form must reproduce the above
12 // copyright notice, this list of conditions and the following
13 // disclaimer in the documentation and/or other materials provided
14 // with the distribution.
15 // * Neither the name of Institue of Cancer Research.
16 // nor the names of its contributors may be used to endorse or promote
17 // products derived from this software without specific prior written
18 // permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 //
32 // For more information on the Plane of Best Fit please see
33 // http://pubs.acs.org/doi/abs/10.1021/ci300293f
34 //
35 // If this code has been useful to you, please include the reference
36 // in any work which has made use of it:
37
38 // Plane of Best Fit: A Novel Method to Characterize the Three-Dimensionality
39 // of Molecules, Nicholas C. Firth, Nathan Brown, and Julian Blagg, Journal of
40 // Chemical Information and Modeling 2012 52 (10), 2516-2525
41
42 //
43 //
44 // Created by Nicholas Firth, November 2011
45 // Modified by Greg Landrum for inclusion in the RDKit distribution November
46 // 2012
47 //
48
49 #include "PBFRDKit.h"
50 #include <Numerics/Matrix.h>
51 #include <Numerics/SquareMatrix.h>
52 #include <Numerics/SymmMatrix.h>
53
54 #include <Eigen/Dense>
55 using namespace RDKit;
56
57 void getSmallestEigenVector(double fSumXX, double fSumXY, double fSumXZ,
58 double fSumYY, double fSumYZ, double fSumZZ,
59 double &x, double &y, double &z);
60
distanceFromAPlane(const RDGeom::Point3D & pt,const std::vector<double> & plane,double denom)61 double distanceFromAPlane(const RDGeom::Point3D &pt,
62 const std::vector<double> &plane, double denom) {
63 double numer = 0.0;
64 numer =
65 std::fabs(pt.x * plane[0] + pt.y * plane[1] + pt.z * plane[2] + plane[3]);
66
67 return numer / denom;
68 }
69
getBestFitPlane(const std::vector<RDGeom::Point3D> & points,std::vector<double> & plane,const std::vector<double> * weights)70 bool getBestFitPlane(const std::vector<RDGeom::Point3D> &points,
71 std::vector<double> &plane,
72 const std::vector<double> *weights) {
73 PRECONDITION((!weights || weights->size() >= points.size()),
74 "bad weights vector");
75 RDGeom::Point3D origin(0, 0, 0);
76 double wSum = 0.0;
77
78 for (unsigned int i = 0; i < points.size(); ++i) {
79 if (weights) {
80 double w = (*weights)[i];
81 wSum += w;
82 origin += points[i] * w;
83 } else {
84 wSum += 1;
85 origin += points[i];
86 }
87 }
88 origin /= wSum;
89
90 double sumXX = 0, sumXY = 0, sumXZ = 0, sumYY = 0, sumYZ = 0, sumZZ = 0;
91 for (unsigned int i = 0; i < points.size(); ++i) {
92 RDGeom::Point3D delta = points[i] - origin;
93 if (weights) {
94 double w = (*weights)[i];
95 delta *= w;
96 }
97 sumXX += delta.x * delta.x;
98 sumXY += delta.x * delta.y;
99 sumXZ += delta.x * delta.z;
100 sumYY += delta.y * delta.y;
101 sumYZ += delta.y * delta.z;
102 sumZZ += delta.z * delta.z;
103 }
104 sumXX /= wSum;
105 sumXY /= wSum;
106 sumXZ /= wSum;
107 sumYY /= wSum;
108 sumYZ /= wSum;
109 sumZZ /= wSum;
110
111 Eigen::Matrix3d mat;
112 mat << sumXX, sumXY, sumXZ, sumXY, sumYY, sumYZ, sumXZ, sumYZ, sumZZ;
113 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(mat);
114 if (eigensolver.info() != Eigen::Success) {
115 BOOST_LOG(rdErrorLog) << "eigenvalue calculation did not converge"
116 << std::endl;
117 return 0.0;
118 }
119 RDGeom::Point3D normal;
120 normal.x = eigensolver.eigenvectors()(0, 0);
121 normal.y = eigensolver.eigenvectors()(1, 0);
122 normal.z = eigensolver.eigenvectors()(2, 0);
123
124 plane[0] = normal.x;
125 plane[1] = normal.y;
126 plane[2] = normal.z;
127 plane[3] = -1 * normal.dotProduct(origin);
128 }
129
PBFRD(ROMol & mol,int confId)130 double PBFRD(ROMol &mol, int confId) {
131 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
132 int numAtoms = mol.getNumAtoms();
133 if (numAtoms < 4) return 0;
134
135 const Conformer &conf = mol.getConformer(confId);
136 if (!conf.is3D()) return 0;
137
138 std::vector<RDGeom::Point3D> points;
139 points.reserve(numAtoms);
140 for (unsigned int i = 0; i < numAtoms; ++i) {
141 points.push_back(conf.getAtomPos(i));
142 }
143
144 std::vector<double> plane(4);
145 getBestFitPlane(points, plane, 0);
146
147 double denom = 0.0;
148 for (unsigned int i = 0; i < 3; ++i) {
149 denom += plane[i] * plane[i];
150 }
151 denom = pow(denom, 0.5);
152
153 double res = 0.0;
154 for (unsigned int i = 0; i < numAtoms; ++i) {
155 res += distanceFromAPlane(points[i], plane, denom);
156 }
157 res /= numAtoms;
158
159 return res;
160 }
161