1 //
2 // Copyright (C) 2019 Greg Landrum
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 #define PY_ARRAY_UNIQUE_SYMBOL rdeht_array_API
11 #include <RDBoost/python.h>
12 #include <boost/python/list.hpp>
13 #include <RDGeneral/Exceptions.h>
14 #include <GraphMol/RDKitBase.h>
15 #include "EHTTools.h"
16
17 #include <RDBoost/boost_numpy.h>
18 #include <RDBoost/PySequenceHolder.h>
19 #include <RDBoost/Wrap.h>
20 #include <RDBoost/import_array.h>
21
22 namespace python = boost::python;
23
24 namespace RDKit {
25
26 namespace {
27
28 // from: https://stackoverflow.com/a/35960614
29 /// @brief Transfer ownership to a Python object. If the transfer fails,
30 /// then object will be destroyed and an exception is thrown.
31 template <typename T>
transfer_to_python(T * t)32 boost::python::object transfer_to_python(T *t) {
33 // Transfer ownership to a smart pointer, allowing for proper cleanup
34 // incase Boost.Python throws.
35 std::unique_ptr<T> ptr(t);
36
37 // Use the manage_new_object generator to transfer ownership to Python.
38 namespace python = boost::python;
39 typename python::manage_new_object::apply<T *>::type converter;
40
41 // Transfer ownership to the Python handler and release ownership
42 // from C++.
43 python::handle<> handle(converter(*ptr));
44 ptr.release();
45
46 return python::object(handle);
47 }
48
getMatrixProp(const double * mat,unsigned int dim1,unsigned int dim2)49 PyObject *getMatrixProp(const double *mat, unsigned int dim1,
50 unsigned int dim2) {
51 if (!mat) {
52 throw_value_error("matrix has not been initialized");
53 }
54 npy_intp dims[2];
55 dims[0] = dim1;
56 dims[1] = dim2;
57
58 auto *res = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
59
60 memcpy(PyArray_DATA(res), static_cast<const void *>(mat),
61 dim1 * dim2 * sizeof(double));
62
63 return PyArray_Return(res);
64 }
getSymmMatrixProp(const double * mat,unsigned int sz)65 PyObject *getSymmMatrixProp(const double *mat, unsigned int sz) {
66 if (!mat) {
67 throw_value_error("matrix has not been initialized");
68 }
69 npy_intp dims[1];
70 dims[0] = sz * (sz + 1) / 2;
71
72 auto *res = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
73
74 memcpy(PyArray_DATA(res), static_cast<const void *>(mat),
75 dims[0] * sizeof(double));
76
77 return PyArray_Return(res);
78 }
getVectorProp(const double * mat,unsigned int sz)79 PyObject *getVectorProp(const double *mat, unsigned int sz) {
80 if (!mat) {
81 throw_value_error("vector has not been initialized");
82 }
83 npy_intp dims[1];
84 dims[0] = sz;
85
86 auto *res = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
87
88 memcpy(PyArray_DATA(res), static_cast<const void *>(mat),
89 sz * sizeof(double));
90
91 return PyArray_Return(res);
92 }
getChargeMatrix(EHTTools::EHTResults & self)93 PyObject *getChargeMatrix(EHTTools::EHTResults &self) {
94 return getMatrixProp(self.reducedChargeMatrix.get(), self.numAtoms,
95 self.numOrbitals);
96 }
97
getOPMatrix(EHTTools::EHTResults & self)98 PyObject *getOPMatrix(EHTTools::EHTResults &self) {
99 return getSymmMatrixProp(self.reducedOverlapPopulationMatrix.get(),
100 self.numAtoms);
101 }
102
getCharges(EHTTools::EHTResults & self)103 PyObject *getCharges(EHTTools::EHTResults &self) {
104 return getVectorProp(self.atomicCharges.get(), self.numAtoms);
105 }
106
getOrbitalEnergies(EHTTools::EHTResults & self)107 PyObject *getOrbitalEnergies(EHTTools::EHTResults &self) {
108 return getVectorProp(self.orbitalEnergies.get(), self.numOrbitals);
109 }
110
runCalc(const RDKit::ROMol & mol,int confId,bool keepMatrices)111 python::tuple runCalc(const RDKit::ROMol &mol, int confId, bool keepMatrices) {
112 auto eRes = new RDKit::EHTTools::EHTResults();
113 bool ok = RDKit::EHTTools::runMol(mol, *eRes, confId, keepMatrices);
114 return python::make_tuple(ok, transfer_to_python(eRes));
115 }
getHamiltonian(EHTTools::EHTResults & self)116 PyObject *getHamiltonian(EHTTools::EHTResults &self) {
117 if (!self.hamiltonianMatrix) {
118 throw_value_error(
119 "Hamiltonian not available, set keepOverlapAndHamiltonianMatrices=True "
120 "to preserve it.");
121 }
122 return getMatrixProp(self.hamiltonianMatrix.get(), self.numOrbitals,
123 self.numOrbitals);
124 }
getOverlapMatrix(EHTTools::EHTResults & self)125 PyObject *getOverlapMatrix(EHTTools::EHTResults &self) {
126 if (!self.overlapMatrix) {
127 throw_value_error(
128 "Overlap matrix not available, set "
129 "keepOverlapAndHamiltonianMatrices=True "
130 "to preserve it.");
131 }
132 return getMatrixProp(self.overlapMatrix.get(), self.numOrbitals,
133 self.numOrbitals);
134 }
135
136 } // end of anonymous namespace
137
138 struct EHT_wrapper {
wrapRDKit::EHT_wrapper139 static void wrap() {
140 std::string docString = "";
141
142 python::class_<RDKit::EHTTools::EHTResults, boost::noncopyable>(
143 "EHTResults", docString.c_str(), python::no_init)
144 .def_readonly("numOrbitals", &RDKit::EHTTools::EHTResults::numOrbitals)
145 .def_readonly("numElectrons",
146 &RDKit::EHTTools::EHTResults::numElectrons)
147 .def_readonly("fermiEnergy", &RDKit::EHTTools::EHTResults::fermiEnergy)
148 .def_readonly("totalEnergy", &RDKit::EHTTools::EHTResults::totalEnergy)
149 .def("GetReducedChargeMatrix", getChargeMatrix,
150 "returns the reduced charge matrix")
151 .def("GetReducedOverlapPopulationMatrix", getOPMatrix,
152 "returns the reduced overlap population matrix")
153 .def("GetAtomicCharges", getCharges,
154 "returns the calculated atomic charges")
155 .def("GetHamiltonian", getHamiltonian,
156 "returns the symmetric Hamiltonian matrix")
157 .def("GetOverlapMatrix", getOverlapMatrix,
158 "returns the symmetric overlap matrix")
159 .def("GetOrbitalEnergies", getOrbitalEnergies,
160 "returns the energies of the molecular orbitals as a vector");
161
162 docString =
163 R"DOC(Runs an extended Hueckel calculation for a molecule.
164 The molecule should have at least one conformation
165
166 ARGUMENTS:
167 - mol: molecule to use
168 - confId: (optional) conformation to use
169 - keepOverlapAndHamiltonianMatrices: (optional) triggers storing the overlap
170 and hamiltonian matrices in the EHTResults object
171
172 RETURNS: a 2-tuple:
173 - a boolean indicating whether or not the calculation succeeded
174 - an EHTResults object with the results
175 )DOC";
176 python::def("RunMol", runCalc,
177 (python::arg("mol"), python::arg("confId") = -1,
178 python::arg("keepOverlapAndHamiltonianMatrices") = false),
179 docString.c_str());
180 }
181 };
182
183 } // end of namespace RDKit
BOOST_PYTHON_MODULE(rdEHTTools)184 BOOST_PYTHON_MODULE(rdEHTTools) {
185 python::scope().attr("__doc__") =
186 R"DOC(Module containing interface to the YAeHMOP extended Hueckel library.
187 Please note that this interface should still be considered experimental and may
188 change from one release to the next.)DOC";
189 rdkit_import_array();
190
191 RDKit::EHT_wrapper::wrap();
192 }