1 // $Id$
2 //
3 // Copyright (C) 2005-2008 Greg Landrum and Rational Discovery LLC
4 //
5 // @@ All Rights Reserved @@
6 // This file is part of the RDKit.
7 // The contents are covered by the terms of the BSD license
8 // which is included in the file license.txt, found at the root
9 // of the RDKit source tree.
10 //
11 #define PY_ARRAY_UNIQUE_SYMBOL rdmoltransforms_array_API
12 #include <RDBoost/python.h>
13 #include <RDBoost/import_array.h>
14 #include "numpy/arrayobject.h"
15 #include <GraphMol/ROMol.h>
16 #include <RDBoost/Wrap.h>
17 #include <GraphMol/Conformer.h>
18 #include <GraphMol/MolTransforms/MolTransforms.h>
19 #include <Geometry/Transform3D.h>
20 #include <Geometry/point.h>
21
22 namespace python = boost::python;
23
24 namespace RDKit {
computeCanonTrans(const Conformer & conf,const RDGeom::Point3D * center=nullptr,bool normalizeCovar=false,bool ignoreHs=true)25 PyObject *computeCanonTrans(const Conformer &conf,
26 const RDGeom::Point3D *center = nullptr,
27 bool normalizeCovar = false, bool ignoreHs = true) {
28 RDGeom::Transform3D *trans;
29 trans = MolTransforms::computeCanonicalTransform(conf, center, normalizeCovar,
30 ignoreHs);
31 npy_intp dims[2];
32 dims[0] = 4;
33 dims[1] = 4;
34 auto *res = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
35 auto *resData = reinterpret_cast<double *>(PyArray_DATA(res));
36 const double *tdata = trans->getData();
37 memcpy(static_cast<void *>(resData), static_cast<const void *>(tdata),
38 4 * 4 * sizeof(double));
39 delete trans;
40 return PyArray_Return(res);
41 }
42
43 #ifdef RDK_HAS_EIGEN3
computePrincAxesMomentsHelper(bool func (const Conformer &,Eigen::Matrix3d &,Eigen::Vector3d &,bool,bool,const std::vector<double> *),const Conformer & conf,bool ignoreHs,const python::object & weights)44 PyObject *computePrincAxesMomentsHelper(bool func(const Conformer &,
45 Eigen::Matrix3d &, Eigen::Vector3d &,
46 bool, bool, const std::vector<double> *),
47 const Conformer &conf,
48 bool ignoreHs, const python::object &weights) {
49 Eigen::Matrix3d axes;
50 Eigen::Vector3d moments;
51 std::vector<double> *weightsVecPtr = nullptr;
52 std::vector<double> weightsVec;
53 size_t i;
54 if (weights != python::object()) {
55 size_t numElements = python::extract<int>(weights.attr("__len__")());
56 if (numElements != conf.getNumAtoms()) {
57 throw ValueErrorException("The Python container must have length equal to conf.GetNumAtoms()");
58 }
59 weightsVec.resize(numElements);
60 for (i = 0; i < numElements; ++i) {
61 weightsVec[i] = python::extract<double>(weights[i]);
62 }
63 weightsVecPtr = &weightsVec;
64 }
65 PyObject *res = PyTuple_New(2);
66 bool success = func(conf, axes, moments, ignoreHs, true, weightsVecPtr);
67 if (success) {
68 npy_intp dims[2];
69 dims[0] = 3;
70 dims[1] = 3;
71 auto *axesNpy = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
72 auto *axesNpyData = reinterpret_cast<double *>(PyArray_DATA(axesNpy));
73 i = 0;
74 for (size_t y = 0; y < 3; ++y) {
75 for (size_t x = 0; x < 3; ++x) {
76 axesNpyData[i++] = axes(y, x);
77 }
78 }
79 auto *momentsNpy = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
80 auto *momentsNpyData = reinterpret_cast<double *>(PyArray_DATA(momentsNpy));
81 for (size_t y = 0; y < 3; ++y) {
82 momentsNpyData[y] = moments(y);
83 }
84 res = PyTuple_New(2);
85 PyTuple_SetItem(res, 0, (PyObject *)axesNpy);
86 PyTuple_SetItem(res, 1, (PyObject *)momentsNpy);
87 }
88 else {
89 PyTuple_SetItem(res, 0, Py_None);
90 PyTuple_SetItem(res, 1, Py_None);
91 }
92 return res;
93 }
94
computePrincAxesMoments(const Conformer & conf,bool ignoreHs,const python::object & weights)95 PyObject *computePrincAxesMoments(const Conformer &conf,
96 bool ignoreHs, const python::object &weights) {
97 return computePrincAxesMomentsHelper(
98 MolTransforms::computePrincipalAxesAndMoments, conf, ignoreHs, weights);
99 }
100
computePrincAxesMomentsFromGyrationMatrix(const Conformer & conf,bool ignoreHs,const python::object & weights)101 PyObject *computePrincAxesMomentsFromGyrationMatrix(const Conformer &conf,
102 bool ignoreHs, const python::object &weights) {
103 return computePrincAxesMomentsHelper(
104 MolTransforms::computePrincipalAxesAndMomentsFromGyrationMatrix, conf, ignoreHs, weights);
105 }
106 #endif
107
transConformer(Conformer & conf,python::object trans)108 void transConformer(Conformer &conf, python::object trans) {
109 PyObject *transObj = trans.ptr();
110 if (!PyArray_Check(transObj)) {
111 throw_value_error("Expecting a numeric array for transformation");
112 }
113 auto *transMat = reinterpret_cast<PyArrayObject *>(transObj);
114 unsigned int nrows = PyArray_DIM(transMat, 0);
115 unsigned int dSize = nrows * nrows;
116 auto *inData = reinterpret_cast<double *>(PyArray_DATA(transMat));
117 RDGeom::Transform3D transform;
118 double *tData = transform.getData();
119 memcpy(static_cast<void *>(tData), static_cast<void *>(inData),
120 dSize * sizeof(double));
121 MolTransforms::transformConformer(conf, transform);
122 }
123 }
124
BOOST_PYTHON_MODULE(rdMolTransforms)125 BOOST_PYTHON_MODULE(rdMolTransforms) {
126 python::scope().attr("__doc__") =
127 "Module containing functions to perform 3D operations like rotate and "
128 "translate conformations";
129
130 rdkit_import_array();
131
132 std::string docString =
133 "Compute the centroid of the conformation - hydrogens are ignored and no attention\n\
134 if paid to the difference in sizes of the heavy atoms\n";
135 python::def("ComputeCentroid", MolTransforms::computeCentroid,
136 (python::arg("conf"), python::arg("ignoreHs") = true),
137 docString.c_str());
138
139 docString =
140 "Compute the transformation required aligna conformer so that\n\
141 the principal axes align up with the x,y, z axes\n\
142 The conformer itself is left unchanged\n\
143 ARGUMENTS:\n\
144 - conf : the conformer of interest\n\
145 - center : optional center point to compute the principal axes around (defaults to the centroid)\n\
146 - normalizeCovar : optionally normalize the covariance matrix by the number of atoms\n";
147 python::def(
148 "ComputeCanonicalTransform", RDKit::computeCanonTrans,
149 (python::arg("conf"), python::arg("center") = (RDGeom::Point3D *)nullptr,
150 python::arg("normalizeCovar") = false, python::arg("ignoreHs") = true),
151 docString.c_str());
152
153 #ifdef RDK_HAS_EIGEN3
154 docString =
155 "Compute principal axes and moments of inertia for a conformer\n\
156 These values are calculated from the inertia tensor:\n\
157 Iij = - sum_{s=1..N}(w_s * r_{si} * r_{sj}) i != j\n\
158 Iii = sum_{s=1..N} sum_{j!=i} (w_s * r_{sj} * r_{sj})\n\
159 where the coordinates are relative to the center of mass.\n\
160 \n\
161 ARGUMENTS:\n\
162 - conf : the conformer of interest\n\
163 - ignoreHs : if True, ignore hydrogen atoms\n\
164 - weights : if present, used to weight the atomic coordinates\n\n\
165 Returns a (principal axes, principal moments) tuple\n";
166 python::def(
167 "ComputePrincipalAxesAndMoments", RDKit::computePrincAxesMoments,
168 (python::arg("conf"), python::arg("ignoreHs") = true,
169 python::arg("weights") = python::object()),
170 docString.c_str());
171
172 docString =
173 "Compute principal axes and moments from the gyration matrix of a conformer\n\
174 These values are calculated from the gyration matrix/tensor:\n\
175 Iij = sum_{s=1..N}(w_s * r_{si} * r_{sj}) i != j\n\
176 Iii = sum_{s=1..N} sum_{t!=s}(w_s * r_{si} * r_{ti})\n\
177 where the coordinates are relative to the center of mass.\n\
178 \n\
179 ARGUMENTS:\n\
180 - conf : the conformer of interest\n\
181 - ignoreHs : if True, ignore hydrogen atoms\n\
182 - weights : if present, used to weight the atomic coordinates\n\n\
183 Returns a (principal axes, principal moments) tuple\n";
184 python::def(
185 "ComputePrincipalAxesAndMomentsFromGyrationMatrix", RDKit::computePrincAxesMomentsFromGyrationMatrix,
186 (python::arg("conf"), python::arg("ignoreHs") = true,
187 python::arg("weights") = python::object()),
188 docString.c_str());
189 #endif
190
191 python::def("TransformConformer", RDKit::transConformer,
192 "Transform the coordinates of a conformer");
193
194 docString =
195 "Canonicalize the orientation of a conformer so that its principal axes\n\
196 around the specified center point coincide with the x, y, z axes\n\
197 \n\
198 ARGUMENTS:\n\
199 - conf : conformer of interest \n\
200 - center : optionally center point about which the principal axes are computed \n\
201 if not specified the centroid of the conformer will be used\n\
202 - normalizeCovar : Optionally normalize the covariance matrix by the number of atoms\n";
203 python::def(
204 "CanonicalizeConformer", MolTransforms::canonicalizeConformer,
205 (python::arg("conf"), python::arg("center") = (RDGeom::Point3D *)nullptr,
206 python::arg("normalizeCovar") = false, python::arg("ignoreHs") = true),
207 docString.c_str());
208
209 python::def("CanonicalizeMol", MolTransforms::canonicalizeMol,
210 (python::arg("mol"), python::arg("normalizeCovar") = false,
211 python::arg("ignoreHs") = true),
212 "Loop over the conformers in a molecule and canonicalize their "
213 "orientation");
214 python::def(
215 "GetBondLength", &MolTransforms::getBondLength,
216 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId")),
217 "Returns the bond length in angstrom between atoms i, j\n");
218 python::def("SetBondLength", &MolTransforms::setBondLength,
219 (python::arg("conf"), python::arg("iAtomId"),
220 python::arg("jAtomId"), python::arg("value")),
221 "Sets the bond length in angstrom between atoms i, j; "
222 "all atoms bonded to atom j are moved\n");
223 python::def("GetAngleRad", &MolTransforms::getAngleRad,
224 (python::arg("conf"), python::arg("iAtomId"),
225 python::arg("jAtomId"), python::arg("kAtomId")),
226 "Returns the angle in radians between atoms i, j, k\n");
227 python::def("GetAngleDeg", &MolTransforms::getAngleDeg,
228 (python::arg("conf"), python::arg("iAtomId"),
229 python::arg("jAtomId"), python::arg("kAtomId")),
230 "Returns the angle in degrees between atoms i, j, k\n");
231 python::def(
232 "SetAngleRad", &MolTransforms::setAngleRad,
233 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId"),
234 python::arg("kAtomId"), python::arg("value")),
235 "Sets the angle in radians between atoms i, j, k; "
236 "all atoms bonded to atom k are moved\n");
237 python::def(
238 "SetAngleDeg", &MolTransforms::setAngleDeg,
239 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId"),
240 python::arg("kAtomId"), python::arg("value")),
241 "Sets the angle in degrees between atoms i, j, k; "
242 "all atoms bonded to atom k are moved\n");
243 python::def(
244 "GetDihedralRad", &MolTransforms::getDihedralRad,
245 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId"),
246 python::arg("kAtomId"), python::arg("lAtomId")),
247 "Returns the dihedral angle in radians between atoms i, j, k, l\n");
248 python::def(
249 "GetDihedralDeg", &MolTransforms::getDihedralDeg,
250 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId"),
251 python::arg("kAtomId"), python::arg("lAtomId")),
252 "Returns the dihedral angle in degrees between atoms i, j, k, l\n");
253 python::def(
254 "SetDihedralRad", &MolTransforms::setDihedralRad,
255 (python::arg("conf"), python::arg("iAtomId"), python::arg("jAtomId"),
256 python::arg("kAtomId"), python::arg("lAtomId"), python::arg("value")),
257 "Sets the dihedral angle in radians between atoms i, j, k, l; "
258 "all atoms bonded to atom l are moved\n");
259 python::def("SetDihedralDeg", &MolTransforms::setDihedralDeg,
260 "Sets the dihedral angle in degrees between atoms i, j, k, l; "
261 "all atoms bonded to atom l are moved\n");
262 }
263