Add Molecular Interaction Fields (#7993)

* Add RealValueVect.

* Add UniformRealValueGrid3D

* Add Molecular Interaction Fields (MIFs)

* line endings

* cherry-pick f1bc94a4c8

* format

* Adapt tests for python3.

* Adapt RealValueVector pickling for python3.

* Speed-up of MIF calculations.

* Bugfix in MIFDescriptors.cpp.

* all tests pass

* clean up some memory leaks

* update copyrights

* rename

* rename the library

* complete the rename

* lost file

* another forgotten file

* cleanup

* clang-tidy

* clang-tidy

* windows DLL builds work

* python wrapper and tests cleanup

* convert to catch2 testing

* switch RealValueVect to use std::vector

* remove obsolete friend

* - Replace explicit loops with stdlib implicit equivalents
- Replace explicit types with auto where possible
- Avoid unnecessary copy operations where possible
- Replace raw pointers with exception-safe unique_ptr
- Replace C-style #define with constexpr
- Replace C-style casts with C++ casts
- Replace C-style arrays with std::vector
- Avoid code duplication with templated operators
- Replace VdWaals class taking multiple atom type definitions and force-field name as string parameter with force-field-specific classes deriving from an abstract VdWaals class
- Replace x,y,z doubles with Point3D class where possible
- Removed unused (and untested) DistanceToClosestAtom class
- Renamed some variables and functions for better clarity
- Converted tabs to spaces
- Made the mol parameter in cube read/write functions optional for convenience
- Made the Python wrappers more pythonic (e.g., avoid C++-style passing objects as parameters which are modified in place)
- Implemented alternative Python class constructors using boost::python::make_constructor rather than with external non-class functions
- The Python wrappers taking a sequence of Point3D now take a sequence of sequences, such that the output of Conformer.GetPositions() can be passed
- Made the Python wrapper sequence parsing more robust
- Removed duplicated code from Python wrappers

* - avoid an unnecessary copy

* progress

* works

* more cleanup

* all tests pass

* changes in response to review

---------

Co-authored-by: dfhahn <dfhahn@users.noreply.github.com>
Co-authored-by: ptosco <paolo.tosco@novartis.com>

.clang-tidy

@@ -1,7 +1,6 @@
 ---
 Checks:          '-*,boost-use-to-string,modernize-use-auto,modernize-use-nullptr,modernize-loop-convert,modernize-use-override,modernize-pass-by-value,modernize-shrink-to-fit,readability-braces-around-statements,readability-simplify-boolean-expr'
 HeaderFilterRegex: 'Code/'
-AnalyzeTemporaryDtors: false
 User:            glandrum
 CheckOptions:
   - key:             modernize-loop-convert.MaxCopySize

Code/DataStructs/CMakeLists.txt

@@ -7,6 +7,7 @@ rdkit_library(DataStructs
               BitVect.cpp SparseBitVect.cpp ExplicitBitVect.cpp Utils.cpp
               base64.cpp BitOps.cpp DiscreteDistMat.cpp
               DiscreteValueVect.cpp FPBReader.cpp MultiFPBReader.cpp
+              RealValueVect.cpp
               LINK_LIBRARIES RDGeneral)
 target_compile_definitions(DataStructs PRIVATE RDKIT_DATASTRUCTS_BUILD)
 
@@ -19,6 +20,7 @@ rdkit_headers(base64.h
               DatastructsStreamOps.h
               DiscreteDistMat.h
               DiscreteValueVect.h
+              RealValueVect.h
               ExplicitBitVect.h
               SparseBitVect.h
               SparseIntVect.h

Code/DataStructs/DiscreteValueVect.h

@@ -41,7 +41,7 @@ class RDKIT_DATASTRUCTS_EXPORT DiscreteValueVect {
     d_valsPerInt = BITS_PER_INT / d_bitsPerVal;
     d_numInts = (length + d_valsPerInt - 1) / d_valsPerInt;
     d_mask = ((1 << d_bitsPerVal) - 1);
-    std::uint32_t *data = new std::uint32_t[d_numInts];
+    auto *data = new std::uint32_t[d_numInts];
     memset(static_cast<void *>(data), 0, d_numInts * sizeof(std::uint32_t));
     d_data.reset(data);
   }

Code/DataStructs/RealValueVect.cpp

@@ -0,0 +1,164 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 "RealValueVect.h"
+#include <RDGeneral/Invariant.h>
+#include <RDGeneral/StreamOps.h>
+#include "DatastructsException.h"
+
+constexpr double VAL_TOL = 0.01;
+
+namespace RDKit {
+const int ci_REALVALUEVECTPICKLE_VERSION = 0x1;
+
+double RealValueVect::getVal(unsigned int i) const {
+  if (i >= d_length) {
+    throw IndexErrorException(i);
+  }
+  return d_data[i];
+}
+
+void RealValueVect::setVal(unsigned int i, double val) {
+  if (i >= d_length) {
+    throw IndexErrorException(i);
+  }
+  d_data[i] = val;
+}
+
+double RealValueVect::getTotalVal() const {
+  return std::accumulate(d_data.begin(), d_data.end(), 0.0,
+                         std::plus<double>());
+}
+
+bool RealValueVect::compareVectors(const RealValueVect &other) const {
+  if (getLength() != other.getLength()) {
+    throw ValueErrorException("Comparing vectors of different lengths");
+  }
+  return std::equal(this->d_data.begin(), this->d_data.end(),
+                    other.d_data.begin(), [](auto v1, auto v2) {
+                      return fabs((v1 - v2) / (v1 != 0 ? v1 : 1)) <= VAL_TOL;
+                    });
+}
+
+double computeL1Norm(const RealValueVect &v1, const RealValueVect &v2) {
+  if (v1.getLength() != v2.getLength()) {
+    throw ValueErrorException("Comparing vectors of different lengths");
+  }
+  double res = 0.0;
+  for (auto i = 0u; i < v1.getLength(); ++i) {
+    res += fabs(v1.getData()[i] - v2.getData()[i]);
+  }
+  return res;
+}
+
+std::string RealValueVect::toString() const {
+  std::stringstream ss(std::ios_base::binary | std::ios_base::out |
+                       std::ios_base::in);
+
+  std::int32_t tVers = ci_REALVALUEVECTPICKLE_VERSION * -1;
+  streamWrite(ss, tVers);
+  std::uint32_t tInt = d_length;
+  streamWrite(ss, tInt);
+
+#if defined(BOOST_BIG_ENDIAN)
+  std::vector<double> td(d_length);
+  std::transform(d_data.begin(), d_data.end(), td.begin(), [](const auto v) {
+    return EndianSwapBytes<HOST_ENDIAN_ORDER, LITTLE_ENDIAN_ORDER>(v);
+  });
+  ss.write(reinterpret_cast<const char *>(td.data()),
+           d_length * sizeof(double));
+  const auto *data = td.data();
+#else
+  const auto *data = d_data.data();
+#endif
+  ss.write(reinterpret_cast<const char *>(data), d_length * sizeof(double));
+  return ss.str();
+};
+
+void RealValueVect::initFromText(const char *pkl, const unsigned int len) {
+  std::stringstream ss(std::ios_base::binary | std::ios_base::in |
+                       std::ios_base::out);
+  ss.write(pkl, len);
+  std::int32_t tVers;
+  streamRead(ss, tVers);
+  tVers *= -1;
+  if (tVers != 0x1) {
+    throw ValueErrorException("bad version in RealValueVect pickle");
+  }
+  std::uint32_t tInt;
+  streamRead(ss, tInt);
+  d_length = tInt;
+  d_data.resize(d_length);
+  auto *data = d_data.data();
+  ss.read(reinterpret_cast<char *>(data), d_length * sizeof(double));
+
+#if defined(BOOST_BIG_ENDIAN)
+  std::transform(d_data.begin(), d_data.end(), d_data.begin(), [](auto &v) {
+    return EndianSwapBytes<LITTLE_ENDIAN_ORDER, HOST_ENDIAN_ORDER>(v);
+  });
+#endif
+};
+
+template <typename O>
+RealValueVect &RealValueVect::applyBinaryOp(const RealValueVect &other, O op) {
+  PRECONDITION(other.d_length == d_length, "length mismatch");
+  for (unsigned int i = 0; i < d_length; ++i) {
+    double v1 = getVal(i);
+    double v2 = other.getVal(i);
+    setVal(i, op(v1, v2));
+  }
+  return *this;
+}
+
+RealValueVect &RealValueVect::operator&=(const RealValueVect &other) {
+  static const double &(*minOp)(const double &, const double &) =
+      std::min<double>;
+  return applyBinaryOp(other, minOp);
+}
+
+RealValueVect &RealValueVect::operator|=(const RealValueVect &other) {
+  static const double &(*maxOp)(const double &, const double &) =
+      std::max<double>;
+  return applyBinaryOp(other, maxOp);
+}
+
+RealValueVect &RealValueVect::operator+=(const RealValueVect &other) {
+  return applyBinaryOp(other, std::plus<double>());
+}
+
+RealValueVect &RealValueVect::operator-=(const RealValueVect &other) {
+  return applyBinaryOp(other, std::minus<double>());
+}
+
+RealValueVect operator+(const RealValueVect &p1, const RealValueVect &p2) {
+  RealValueVect res(p1);
+  res += p2;
+  return res;
+}
+
+RealValueVect operator-(const RealValueVect &p1, const RealValueVect &p2) {
+  RealValueVect res(p1);
+  res -= p2;
+  return res;
+}
+
+RealValueVect operator|(const RealValueVect &p1, const RealValueVect &p2) {
+  RealValueVect res(p1);
+  res |= p2;
+  return res;
+};
+
+RealValueVect operator&(const RealValueVect &p1, const RealValueVect &p2) {
+  RealValueVect res(p1);
+  res &= p2;
+  return res;
+};
+
+}  // end of namespace RDKit

Code/DataStructs/RealValueVect.h

@@ -0,0 +1,137 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <RDGeneral/export.h>
+#ifndef RD_REAL_VALUE_VECT_20140407
+#define RD_REAL_VALUE_VECT_20140407
+
+#include <boost/smart_ptr.hpp>
+#include <string>
+#include <cstring>
+
+namespace RDGeom {
+class UniformRealValueGrid3D;
+}
+namespace RDKit {
+
+//! a class for efficiently storing vectors of double values
+//! Has additional features compared to std::vector<double>:
+//! construct from and write to pickle
+class RDKIT_DATASTRUCTS_EXPORT RealValueVect {
+ public:
+  RealValueVect() = default;
+  //! initialize with a particular size
+  RealValueVect(unsigned int length) : d_length(length) {
+    d_data.resize(d_length, 0.0);
+  }
+
+  //! initialize with a particular size and initial value
+  RealValueVect(unsigned int length, double val) : d_length(length) {
+    d_data.resize(d_length, val);
+  }
+
+  //! constructor from a pickle
+  RealValueVect(const std::string &pkl) {
+    initFromText(pkl.c_str(), pkl.size());
+  };
+
+  //! constructor from a pickle
+  RealValueVect(const char *pkl, unsigned int len) { initFromText(pkl, len); };
+
+  //! return the value at an index
+  double getVal(unsigned int i) const;
+
+  //! support indexing using []
+  double operator[](unsigned int idx) const { return getVal(idx); };
+
+  //! set the value at an index
+  void setVal(unsigned int i, double val);
+
+  //! returns the sum of all the elements in the vect
+  double getTotalVal() const;
+
+  //! returns the length
+  unsigned int getLength() const { return d_length; };
+
+  //! returns the length
+  unsigned int size() const { return getLength(); };
+
+  //! compares 2 vectors and returns false if different
+  bool compareVectors(const RealValueVect &other) const;
+
+  //! in-place operator&
+  RealValueVect &operator&=(const RealValueVect &other);
+
+  //! in-place operator|
+  RealValueVect &operator|=(const RealValueVect &other);
+
+  //! in-place operator+
+  RealValueVect &operator+=(const RealValueVect &other);
+
+  //! in-place operator-
+  RealValueVect &operator-=(const RealValueVect &other);
+
+  //! returns a binary string representation (pickle)
+  std::string toString() const;
+
+  void setLength(unsigned int sz) {
+    d_length = sz;
+    d_data.resize(sz);
+  }
+  void setToVal(double val) { std::fill(d_data.begin(), d_data.end(), val); }
+
+  const std::vector<double> &getData() const { return d_data; }
+  std::vector<double> &getData() { return d_data; }
+
+ private:
+  void initFromText(const char *pkl, const unsigned int len);
+  unsigned int d_length = 0;
+  std::vector<double> d_data;
+
+  template <typename O>
+  RealValueVect &applyBinaryOp(const RealValueVect &other, O op);
+};  // end of declaration of class RealValueVect
+
+//! returns L1 Norm of vectors
+RDKIT_DATASTRUCTS_EXPORT double computeL1Norm(const RealValueVect &v1,
+                                              const RealValueVect &v2);
+
+//! returns the sum of vectors
+RDKIT_DATASTRUCTS_EXPORT RealValueVect operator+(const RealValueVect &p1,
+                                                 const RealValueVect &p2);
+
+//! returns the difference of vectors
+RDKIT_DATASTRUCTS_EXPORT RealValueVect operator-(const RealValueVect &p1,
+                                                 const RealValueVect &p2);
+
+//! support rvv3 = rvv1|rvv2
+/*!
+
+   operator| returns the maximum value for each element.
+   e.g.:
+     [0,1,2,0] | [0,1,1,1] -> [0,1,2,1]
+
+ */
+RDKIT_DATASTRUCTS_EXPORT RealValueVect operator|(const RealValueVect &p1,
+                                                 const RealValueVect &p2);
+
+//! support rvv3 = rvv1&rvv2
+/*!
+   operator& returns the minimum value for each element.
+   e.g.:
+     [0,1,2,0] & [0,1,1,1] -> [0,1,1,0]
+
+ */
+RDKIT_DATASTRUCTS_EXPORT RealValueVect operator&(const RealValueVect &p1,
+                                                 const RealValueVect &p2);
+
+}  // end of namespace RDKit
+
+#endif

Code/DataStructs/Wrap/CMakeLists.txt

@@ -1,6 +1,7 @@
 remove_definitions(-DRDKIT_DATASTRUCTS_BUILD)
 rdkit_python_extension(cDataStructs
                        DataStructs.cpp DiscreteValueVect.cpp SparseIntVect.cpp
+                       RealValueVect.cpp
                        wrap_SparseBV.cpp wrap_ExplicitBV.cpp wrap_BitOps.cpp
                        wrap_FPB.cpp
                        wrap_Utils.cpp
@@ -12,6 +13,8 @@ add_pytest(pyBV
          ${CMAKE_CURRENT_SOURCE_DIR}/testBV.py)
 add_pytest(pyDiscreteValueVect
          ${CMAKE_CURRENT_SOURCE_DIR}/testDiscreteValueVect.py)
+add_pytest(pyRealValueVect
+         ${CMAKE_CURRENT_SOURCE_DIR}/testRealValueVect.py)
 add_pytest(pySparseIntVect
          ${CMAKE_CURRENT_SOURCE_DIR}/testSparseIntVect.py)
 add_pytest(pyFPB

Code/DataStructs/Wrap/DataStructs.cpp

@@ -15,7 +15,7 @@
 #include <DataStructs/BitVects.h>
 #include <DataStructs/DiscreteValueVect.h>
 #include <DataStructs/SparseIntVect.h>
-
+#include <DataStructs/RealValueVect.h>
 #include <RDBoost/boost_numpy.h>
 #include <numpy/npy_common.h>
 #include <RDBoost/import_array.h>
@@ -28,11 +28,13 @@ void wrap_EBV();
 void wrap_BitOps();
 void wrap_Utils();
 void wrap_discreteValVect();
+void wrap_realValVect();
 void wrap_sparseIntVect();
 void wrap_FPB();
 
-template <typename T>
-void convertToNumpyArray(const T &v, python::object destArray) {
+namespace {
+template <typename T, typename U>
+void converter(const T &v, python::object destArray, U func) {
   if (!PyArray_Check(destArray.ptr())) {
     throw_value_error("Expecting a Numeric array object");
   }
@@ -44,12 +46,22 @@ void convertToNumpyArray(const T &v, python::object destArray) {
   dims.len = 1;
   PyArray_Resize(destP, &dims, 0, NPY_ANYORDER);
   for (unsigned int i = 0; i < v.size(); ++i) {
-    PyObject *iItem = PyInt_FromLong(v[i]);
+    PyObject *iItem = func(v[i]);
     PyArray_SETITEM(destP, static_cast<char *>(PyArray_GETPTR1(destP, i)),
                     iItem);
     Py_DECREF(iItem);
   }
 }
+}  // namespace
+
+template <typename T>
+void convertToIntNumpyArray(const T &v, python::object destArray) {
+  converter(v, destArray, PyLong_FromLong);
+}
+template <typename T>
+void convertToDoubleNumpyArray(const T &v, python::object destArray) {
+  converter(v, destArray, PyFloat_FromDouble);
+}
 
 BOOST_PYTHON_MODULE(cDataStructs) {
   rdkit_import_array();
@@ -72,31 +84,36 @@ BOOST_PYTHON_MODULE(cDataStructs) {
   wrap_EBV();
   wrap_BitOps();
   wrap_discreteValVect();
+  wrap_realValVect();
   wrap_sparseIntVect();
   wrap_FPB();
 
   python::def(
       "ConvertToNumpyArray",
-      (void (*)(const ExplicitBitVect &, python::object))convertToNumpyArray,
+      (void (*)(const ExplicitBitVect &, python::object))convertToIntNumpyArray,
       (python::arg("bv"), python::arg("destArray")));
   python::def("ConvertToNumpyArray",
               (void (*)(const RDKit::DiscreteValueVect &,
-                        python::object))convertToNumpyArray,
+                        python::object))convertToIntNumpyArray,
               (python::arg("bv"), python::arg("destArray")));
   python::def("ConvertToNumpyArray",
               (void (*)(const RDKit::SparseIntVect<std::int32_t> &,
-                        python::object))convertToNumpyArray,
+                        python::object))convertToIntNumpyArray,
               (python::arg("bv"), python::arg("destArray")));
   python::def("ConvertToNumpyArray",
               (void (*)(const RDKit::SparseIntVect<boost::int64_t> &,
-                        python::object))convertToNumpyArray,
+                        python::object))convertToIntNumpyArray,
               (python::arg("bv"), python::arg("destArray")));
   python::def("ConvertToNumpyArray",
               (void (*)(const RDKit::SparseIntVect<std::uint32_t> &,
-                        python::object))convertToNumpyArray,
+                        python::object))convertToIntNumpyArray,
               (python::arg("bv"), python::arg("destArray")));
   python::def("ConvertToNumpyArray",
               (void (*)(const RDKit::SparseIntVect<boost::uint64_t> &,
-                        python::object))convertToNumpyArray,
+                        python::object))convertToIntNumpyArray,
               (python::arg("bv"), python::arg("destArray")));
+  python::def("ConvertToNumpyArray",
+              (void (*)(const RDKit::RealValueVect &,
+                        python::object))convertToDoubleNumpyArray,
+              (python::arg("rvv"), python::arg("destArray")));
 }

Code/DataStructs/Wrap/RealValueVect.cpp

@@ -0,0 +1,73 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <boost/python.hpp>
+
+#include <RDGeneral/types.h>
+#include <RDGeneral/Invariant.h>
+#include <RDBoost/PySequenceHolder.h>
+#include <DataStructs/RealValueVect.h>
+
+using namespace RDKit;
+
+struct rvv_pickle_suite : rdkit_pickle_suite {
+  static python::tuple getinitargs(const RealValueVect &self) {
+    std::string res = self.toString();
+    python::object retval = python::object(
+        python::handle<>(PyBytes_FromStringAndSize(res.c_str(), res.length())));
+    return python::make_tuple(retval);
+  };
+};
+
+std::string realValVectDoc =
+    "A container class for storing real\n\
+values.\n\
+\n\
+The length of the vector is set at construction time.\n\
+\n\
+As you would expect, _RealValueVects_ support a set of binary operations\n\
+so you can do things like:\n\
+  rvv3 = rvv1 & rvv2  the result contains the smallest value in each entry\n\
+  rvv3 = rvv1 | rvv2  the result contains the largest value in each entry\n\
+  rvv1 += rvv2     \n\
+  rvv3 = rvv1 + rvv2    \n\
+  rvv1 -= rvv3    \n\
+  rvv3 = rvv1 - rvv2    \n\
+\n\
+Elements can be set and read using indexing (i.e. bv[i] = 4 or val=bv[i])\n\
+\n";
+
+struct realValVec_wrapper {
+  static void wrap() {
+    python::class_<RealValueVect>("RealValueVect", realValVectDoc.c_str(),
+                                  python::init<unsigned int>("Constructor"))
+        .def(python::init<std::string>())
+        .def("__len__", &RealValueVect::getLength,
+             "Get the number of entries in the vector")
+        .def("__setitem__", &RealValueVect::setVal,
+             "Set the value at a specified location")
+        .def("__getitem__", &RealValueVect::getVal,
+             "Get the value at a specified location")
+        .def(python::self & python::self)
+        .def(python::self | python::self)
+        .def(python::self - python::self)
+        .def(python::self -= python::self)
+        .def(python::self + python::self)
+        .def(python::self += python::self)
+        .def("GetTotalVal", &RealValueVect::getTotalVal,
+             "Get the sum of the values in the vector, basically L1 norm")
+        .def_pickle(rvv_pickle_suite());
+    python::def("ComputeL1Norm", computeL1Norm,
+                "Compute the distance between two real vector values\n");
+  }
+};
+
+void wrap_realValVect() { realValVec_wrapper::wrap(); }

Code/DataStructs/Wrap/testData/rvvs.pkl

@@ -0,0 +1,6 @@
+crdkit.DataStructs.cDataStructs
+RealValueVect
+p1
+(S'\xff\xff\xff\xff\x1e\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00\xcd\xcc\xcc\xcc\xcc\xcc\xf4?\x00\x00\x00\x00\x00\x00\x00\x00'
+tRp2
+.

Code/DataStructs/Wrap/testRealValueVect.py

@@ -0,0 +1,118 @@
+# testRealValueVect.py
+#
+# 14.04.2014 by David Hahn, hahnda6
+
+from __future__ import print_function, unicode_literals
+from rdkit import RDConfig
+import os, sys, pickle
+import unittest
+from rdkit import DataStructs as ds
+
+
+class TestCase(unittest.TestCase):
+
+  def setUp(self):
+    pass
+
+  def test1Real(self):
+    v1 = ds.RealValueVect(30)
+    for i in range(30):
+      v1[i] = i / 10.0
+
+    self.assertTrue(len(v1) == 30)
+    self.assertAlmostEqual(v1.GetTotalVal(), 43.5)
+
+    for i in range(len(v1)):
+      self.assertAlmostEqual(v1[i], 0.1 * i)
+
+  def test2RealVectDistances(self):
+    v1 = ds.RealValueVect(30)
+    v2 = ds.RealValueVect(30)
+    for i in range(15):
+      v1[2 * i] = 1.3
+      v2[2 * i] = 1.3
+    self.assertAlmostEqual(ds.ComputeL1Norm(v1, v2), 0)
+    for i in range(30):
+      if (i % 3 == 0):
+        v2[i] = 1.3
+      else:
+        v2[i] = 0
+    self.assertAlmostEqual(ds.ComputeL1Norm(v1, v2), 19.5)
+
+  def test3Pickles(self):
+    #outF = file('../testData/rvvs.pkl','wb+')
+    with open(os.path.join(RDConfig.RDBaseDir, 'Code/DataStructs/Wrap/testData/rvvs.pkl'),
+              'rb') as inF:
+      v1 = ds.RealValueVect(30)
+      for i in range(15):
+        v1[2 * i] = 1.3
+      v2 = pickle.loads(pickle.dumps(v1))
+      self.assertAlmostEqual(ds.ComputeL1Norm(v1, v2), 0)
+      #pickle.dump(v1,outF)
+      v2 = pickle.load(inF, encoding='bytes')
+      self.assertAlmostEqual(ds.ComputeL1Norm(v1, v2), 0)
+      self.assertAlmostEqual(v1.GetTotalVal(), v2.GetTotalVal())
+      self.assertTrue(v2.GetTotalVal() != 0)
+    #outF.close()
+
+  def test4RealVectOps(self):
+    v1 = ds.RealValueVect(8)
+    for i in range(4):
+      v1[2 * i] = 2.1
+    self.assertAlmostEqual(v1.GetTotalVal(), 8.4)
+    v2 = ds.RealValueVect(8)
+    for i in range(4):
+      v2[2 * i + 1] = 2.1
+      v2[2 * i] = 1.1
+    self.assertAlmostEqual(v2.GetTotalVal(), 12.8)
+
+    v3 = v1 | v2
+    self.assertTrue(len(v3) == len(v2))
+    self.assertAlmostEqual(v3.GetTotalVal(), 16.8)
+
+    v3 = v1 & v2
+    self.assertTrue(len(v3) == len(v2))
+    self.assertAlmostEqual(v3.GetTotalVal(), 4.4)
+
+    v4 = v1 + v2
+    self.assertTrue(len(v4) == len(v2))
+    self.assertAlmostEqual(v4.GetTotalVal(), 21.2)
+
+    v4 = v1 - v2
+    self.assertAlmostEqual(v4.GetTotalVal(), -4.4)
+    v4 = v2 - v1
+    self.assertAlmostEqual(v4.GetTotalVal(), 4.4)
+
+    v4 = v2
+    v4 -= v1
+    self.assertAlmostEqual(v4.GetTotalVal(), 4.4)
+    v4 -= v4
+    self.assertAlmostEqual(v4.GetTotalVal(), 0)
+
+  def testIterator(self):
+    v1 = ds.RealValueVect(30)
+    for i in range(15):
+      v1[2 * i] = 1.1
+    l1 = list(v1)
+    self.assertTrue(len(l1) == len(v1))
+    for i, v in enumerate(v1):
+      self.assertAlmostEqual(l1[i], v)
+    with self.assertRaises(IndexError):
+      v1[40]
+
+  def test9ToNumpy(self):
+    import numpy
+    bv = ds.RealValueVect(32)
+    bv[0] = 1.1
+    bv[1] = 4.4
+    bv[17] = 1.2
+    bv[23] = 8.3
+    bv[31] = 12.2
+    arr = numpy.zeros((32, ), numpy.double)
+    ds.ConvertToNumpyArray(bv, arr)
+    for i in range(len(bv)):
+      self.assertAlmostEqual(bv[i], arr[i])
+
+
+if __name__ == '__main__':
+  unittest.main()

Code/DataStructs/testDatastructs.cpp

@@ -19,6 +19,7 @@
 #include "base64.h"
 #include <cmath>
 #include "DiscreteValueVect.h"
+#include "RealValueVect.h"
 #include <RDGeneral/Invariant.h>
 #include <RDGeneral/RDLog.h>
 #include <RDGeneral/Exceptions.h>
@@ -1450,6 +1451,117 @@ void test17Github3994() {
   TEST_ASSERT(siv1 == siv2);
 }
 
+void test14RealVect() {
+  RealValueVect vect1(30);
+  unsigned int i;
+  for (i = 0; i < 15; ++i) {
+    vect1.setVal(2 * i, 1.0);
+  }
+
+  TEST_ASSERT(vect1.getLength() == 30);
+  TEST_ASSERT(feq(vect1.getTotalVal(), 15.0));
+  for (i = 0; i < vect1.getLength(); ++i) {
+    TEST_ASSERT(feq(vect1.getVal(i), (i + 1) % 2));
+  }
+
+  RealValueVect vect2(30);
+  for (i = 0; i < vect2.getLength(); ++i) {
+    vect2.setVal(i, double(1.0 / (i + 1.0)));
+  }
+
+  TEST_ASSERT(vect2.getLength() == 30);
+  for (i = 0; i < vect2.getLength(); ++i) {
+    TEST_ASSERT(feq(vect2.getVal(i), double(1.0 / (i + 1.0))));
+  }
+
+  // test copy constructor and operator[]
+  RealValueVect vect3(vect2);
+  TEST_ASSERT(vect3.getLength() == 30);
+  for (i = 0; i < vect3.getLength(); ++i) {
+    TEST_ASSERT(feq(vect3[i], double(1.0 / (i + 1.0))));
+  }
+
+  double Pi = 3.141592;
+  RealValueVect vect4(60, Pi);
+  TEST_ASSERT(feq(vect4.getTotalVal(), 60 * Pi));
+}
+
+void test15RealVectDists() {
+  RealValueVect v1(30);
+  RealValueVect v2(30);
+  unsigned int i;
+  for (i = 0; i < 15; ++i) {
+    v1.setVal(2 * i, 1.0);
+    v2.setVal(2 * i, 1.0);
+  }
+  TEST_ASSERT(feq(computeL1Norm(v1, v2), 0));
+  for (i = 0; i < 30; ++i) {
+    v2.setVal(i, i % 2);
+  }
+
+  TEST_ASSERT(feq(computeL1Norm(v1, v2), 30.0));
+
+  for (i = 0; i < 30; ++i) {
+    if (i % 3 == 0) {
+      v2.setVal(i, 1.0);
+    } else {
+      v2.setVal(i, 0.0);
+    }
+  }
+
+  TEST_ASSERT(feq(computeL1Norm(v1, v2), 15.0));
+
+  for (i = 0; i < 30; ++i) {
+    v1.setVal(i, 0.0);
+    v2.setVal(i, i / 10.0);
+  }
+
+  TEST_ASSERT(feq(computeL1Norm(v1, v2), 43.5))
+}
+
+void test16RealVectPickles() {
+  RealValueVect v1(30);
+  unsigned int i;
+  for (i = 0; i < 15; ++i) {
+    v1.setVal(2 * i, 1.1);
+  }
+  RealValueVect v2(v1.toString());
+  TEST_ASSERT(feq(computeL1Norm(v1, v2), 0.0));
+}
+
+void test17RealVectOps() {
+  RealValueVect vect1(8);
+  for (unsigned int i = 0; i < 4; ++i) {
+    vect1.setVal(2 * i, 2.1);
+  }
+  TEST_ASSERT(vect1.getLength() == 8);
+  TEST_ASSERT(feq(vect1.getTotalVal(), 8.4));
+
+  RealValueVect vect2(8);
+  for (unsigned int i = 0; i < 4; ++i) {
+    vect2.setVal(2 * i + 1, 2.2);
+    vect2.setVal(2 * i, 1.1);
+  }
+  TEST_ASSERT(feq(vect2.getTotalVal(), 13.2));
+
+  RealValueVect vect3 = vect1 & vect2;
+  TEST_ASSERT(vect3.getLength() == 8);
+  TEST_ASSERT(feq(vect3.getTotalVal(), 4.4));
+
+  RealValueVect vect4 = vect1 | vect2;
+  TEST_ASSERT(vect4.getLength() == 8);
+  TEST_ASSERT(feq(vect4.getTotalVal(), 17.2));
+
+  RealValueVect vect5 = vect1 + vect2;
+  TEST_ASSERT(vect5.getLength() == 8);
+  TEST_ASSERT(feq(vect5.getTotalVal(), 21.6));
+
+  vect5 = vect1 - vect2;
+  TEST_ASSERT(feq(vect5.getTotalVal(), -4.8));
+  vect5 = vect2 - vect1;
+  TEST_ASSERT(feq(vect5.getTotalVal(), 4.8));
+}
+
 int main() {
   RDLog::InitLogs();
   boost::logging::enable_logs("rdApp.info");
@@ -1561,5 +1673,21 @@ int main() {
                        << std::endl;
   test17Github3994();
 
+  BOOST_LOG(rdInfoLog)
+      << " Test RealValue Vectors 1: Constructors and Accessing Values -------------------------------"
+      << std::endl;
+  test14RealVect();
+  BOOST_LOG(rdInfoLog)
+      << " Test RealValue Vectors 2: L1Norm -------------------------------"
+      << std::endl;
+  test15RealVectDists();
+  BOOST_LOG(rdInfoLog)
+      << " Test RealValue Vectors 3: Pickles -------------------------------"
+      << std::endl;
+  test16RealVectPickles();
+  BOOST_LOG(rdInfoLog)
+      << " Test RealValue Vectors 4: Vector Operations  -------------------------------"
+      << std::endl;
+  test17RealVectOps();
   return 0;
 }

Code/Geometry/CMakeLists.txt

@@ -1,7 +1,8 @@
 
 rdkit_library(RDGeometryLib 
               point.cpp Transform2D.cpp Transform3D.cpp 
-              UniformGrid3D.cpp GridUtils.cpp
+              UniformGrid3D.cpp UniformRealValueGrid3D.cpp
+              GridUtils.cpp
               LINK_LIBRARIES DataStructs RDGeneral)
 target_compile_definitions(RDGeometryLib PRIVATE RDKIT_RDGEOMETRYLIB_BUILD)
 
@@ -12,13 +13,16 @@ rdkit_headers(Grid3D.h
               Transform3D.h
               Transform.h
               UniformGrid3D.h
+              UniformRealValueGrid3D.h
               Utils.h DEST Geometry)
 
 rdkit_test(testTransforms testTransforms.cpp LINK_LIBRARIES RDGeometryLib )
 rdkit_test(testGrid testGrid.cpp LINK_LIBRARIES RDGeometryLib)
+rdkit_test(testRealValueGrid testRealValueGrid.cpp LINK_LIBRARIES RDGeometryLib DataStructs RDGeneral)
 rdkit_catch_test(geometryTestsCatch catch_tests.cpp
            LINK_LIBRARIES RDGeometryLib)
 
+
 if(RDK_BUILD_PYTHON_WRAPPERS)
 add_subdirectory(Wrap)
 endif()

Code/Geometry/Grid3D.h

@@ -35,20 +35,19 @@ class RDKIT_RDGEOMETRYLIB_EXPORT GridException : public std::exception {
 };
 
 //! Virtual base class for a grid object
+template <typename VectorType, typename ValueType1, typename ValueType2>
 class RDKIT_RDGEOMETRYLIB_EXPORT Grid3D {
  public:
-  virtual ~Grid3D() {}
+  virtual ~Grid3D(){};
   virtual int getGridPointIndex(const Point3D &point) const = 0;
-  virtual int getVal(const Point3D &point) const = 0;
-  virtual void setVal(const Point3D &point, unsigned int val) = 0;
+  virtual ValueType1 getVal(const Point3D &point) const = 0;
+  virtual void setVal(const Point3D &point, ValueType2 val) = 0;
 
   virtual Point3D getGridPointLoc(unsigned int pointId) const = 0;
-  virtual unsigned int getVal(unsigned int pointId) const = 0;
-  virtual void setVal(unsigned int pointId, unsigned int val) = 0;
-
+  virtual ValueType2 getVal(unsigned int pointId) const = 0;
+  virtual void setVal(unsigned int pointId, ValueType2 val) = 0;
   virtual unsigned int getSize() const = 0;
-
-  virtual const RDKit::DiscreteValueVect *getOccupancyVect() const = 0;
+  virtual const VectorType *getOccupancyVect() const = 0;
 };
 }  // namespace RDGeom
 

Code/Geometry/UniformGrid3D.cpp

@@ -73,11 +73,11 @@ void UniformGrid3D::initGrid(
 
 UniformGrid3D::UniformGrid3D(const std::string &pkl) {
   dp_storage = nullptr;
-  this->initFromText(pkl.c_str(), pkl.size());
+  initFromText(pkl.c_str(), pkl.size());
 }
 UniformGrid3D::UniformGrid3D(const char *pkl, const unsigned int len) {
   dp_storage = nullptr;
-  this->initFromText(pkl, len);
+  initFromText(pkl, len);
 }
 
 UniformGrid3D::~UniformGrid3D() {
@@ -181,7 +181,7 @@ bool UniformGrid3D::compareParams(const UniformGrid3D &other) const {
 void UniformGrid3D::setSphereOccupancy(const Point3D &center, double radius,
                                        double stepSize, int maxNumLayers,
                                        bool ignoreOutOfBound) {
-  int ptIndex = this->getGridPointIndex(center);
+  int ptIndex = getGridPointIndex(center);
   if (ptIndex == -1) {
     if (ignoreOutOfBound) {
       return;

Code/Geometry/UniformGrid3D.h

@@ -8,8 +8,8 @@
 //  of the RDKit source tree.
 //
 #include <RDGeneral/export.h>
-#ifndef _UNIFORMGRID3D_H_20050124_1703
-#define _UNIFORMGRID3D_H_20050124_1703
+#ifndef UNIFORMGRID3D_H_20050124_1703
+#define UNIFORMGRID3D_H_20050124_1703
 
 #include "point.h"
 #include <DataStructs/DiscreteValueVect.h>
@@ -17,7 +17,8 @@
 #include <iostream>
 
 namespace RDGeom {
-class RDKIT_RDGEOMETRYLIB_EXPORT UniformGrid3D : public Grid3D {
+class RDKIT_RDGEOMETRYLIB_EXPORT UniformGrid3D
+    : public Grid3D<RDKit::DiscreteValueVect, int, unsigned int> {
  public:
   //! \brief ctor
   /*
@@ -28,13 +29,13 @@ class RDKIT_RDGEOMETRYLIB_EXPORT UniformGrid3D : public Grid3D {
       \param valType: the data type of the grid (determines the number of bits
                       per point)
       \param offset:  OPTIONAL: the offset of the grid from (0,0,0), in
-     Angstroms.
+    Angstroms.
 
-      \b Note: the values of arguments such as \c dimX and \c spacing
-      don't actually need to be in Angstroms, but they should be internally
-      consistent.
+    \b Note: the values of arguments such as \c dimX and \c spacing
+    don't actually need to be in Angstroms, but they should be internally
+    consistent.
 
-  */
+*/
   UniformGrid3D(double dimX, double dimY, double dimZ, double spacing = 0.5,
                 RDKit::DiscreteValueVect::DiscreteValueType valType =
                     RDKit::DiscreteValueVect::TWOBITVALUE,

Code/Geometry/UniformRealValueGrid3D.cpp

@@ -0,0 +1,319 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 "UniformRealValueGrid3D.h"
+#include <DataStructs/RealValueVect.h>
+#include <RDGeneral/StreamOps.h>
+#include <RDGeneral/Exceptions.h>
+#include <Geometry/point.h>
+#include <fstream>
+#include <iomanip>
+
+using namespace RDKit;
+
+namespace {
+constexpr double OFFSET_TOL = 1.e-4;
+constexpr double SPACING_TOL = 1.e-4;
+constexpr const char *INCOMPATIBLE_GRIDS = "incompatible grids";
+}  // end anonymous namespace
+
+namespace RDGeom {
+std::int32_t ci_RealValueGrid3DPICKLE_VERSION = 0x1;
+UniformRealValueGrid3D::UniformRealValueGrid3D(
+    const UniformRealValueGrid3D &other) {
+  UniformRealValueGrid3D::initGrid(
+      other.d_numX * other.d_spacing, other.d_numY * other.d_spacing,
+      other.d_numZ * other.d_spacing, other.d_spacing, other.d_offSet,
+      &other.d_storage);
+}
+
+void UniformRealValueGrid3D::initGrid(double dimX, double dimY, double dimZ,
+                                      double spacing,
+                                      const RDGeom::Point3D &offSet,
+                                      const RDKit::RealValueVect *data) {
+  PRECONDITION(dimX > 0.0, "Invalid x-dimension for grid");
+  PRECONDITION(dimY > 0.0, "Invalid y-dimension for grid");
+  PRECONDITION(dimZ > 0.0, "Invalid z-dimension for grid");
+  PRECONDITION(spacing > 0.0, "Invalid spacing for grid");
+  d_numX = static_cast<unsigned int>(floor(dimX / spacing + 0.5));
+  d_numY = static_cast<unsigned int>(floor(dimY / spacing + 0.5));
+  d_numZ = static_cast<unsigned int>(floor(dimZ / spacing + 0.5));
+  // PRECONDITION((!data)||data->getValueType()==double,"grid data type
+  // mismatch");
+  PRECONDITION(!data || data->size() == d_numX * d_numY * d_numZ,
+               "grid data size mismatch");
+
+  d_spacing = spacing;
+  d_offSet = offSet;
+  if (!data) {
+    d_storage.setLength(d_numX * d_numY * d_numZ);
+    d_storage.setToVal(0.0);
+  } else {
+    d_storage = *data;
+  }
+}
+
+UniformRealValueGrid3D::UniformRealValueGrid3D(const std::string &pkl) {
+  initFromText(pkl.c_str(), pkl.size());
+}
+
+UniformRealValueGrid3D::UniformRealValueGrid3D(const char *pkl,
+                                               unsigned int len) {
+  initFromText(pkl, len);
+}
+
+int UniformRealValueGrid3D::getGridIndex(unsigned int xi, unsigned int yi,
+                                         unsigned int zi) const {
+  if (xi >= d_numX) {
+    return -1;
+  }
+  if (yi >= d_numY) {
+    return -1;
+  }
+  if (zi >= d_numZ) {
+    return -1;
+  }
+  return (zi * d_numX * d_numY + yi * d_numX + xi);
+}
+
+void UniformRealValueGrid3D::getGridIndices(unsigned int idx, unsigned int &xi,
+                                            unsigned int &yi,
+                                            unsigned int &zi) const {
+  if (idx >= d_numX * d_numY * d_numZ) {
+    throw IndexErrorException(idx);
+  }
+  xi = idx % d_numX;
+  yi = (idx % (d_numX * d_numY)) / d_numX;
+  zi = idx / (d_numX * d_numY);
+}
+
+int UniformRealValueGrid3D::getGridPointIndex(
+    const RDGeom::Point3D &point) const {
+  auto tPt = point;
+  tPt -= d_offSet;  // d_origin;
+  tPt /= d_spacing;
+  constexpr double move = 0.5;
+  auto xi = static_cast<int>(floor(tPt.x + move));
+  auto yi = static_cast<int>(floor(tPt.y + move));
+  auto zi = static_cast<int>(floor(tPt.z + move));
+
+  if ((xi < 0) || (xi >= static_cast<int>(d_numX))) {
+    return -1;
+  }
+  if ((yi < 0) || (yi >= static_cast<int>(d_numY))) {
+    return -1;
+  }
+  if ((zi < 0) || (zi >= static_cast<int>(d_numZ))) {
+    return -1;
+  }
+
+  return (zi * d_numX * d_numY + yi * d_numX + xi);
+}
+
+double UniformRealValueGrid3D::getVal(const RDGeom::Point3D &point) const {
+  auto id = getGridPointIndex(point);
+  if (id < 0) {
+    return -1;
+  }
+  return d_storage.getVal(static_cast<unsigned int>(id));
+}
+
+double UniformRealValueGrid3D::getVal(unsigned int pointId) const {
+  return d_storage.getVal(pointId);
+}
+
+void UniformRealValueGrid3D::setVal(const RDGeom::Point3D &point, double val) {
+  auto id = getGridPointIndex(point);
+  if (id < 0) {
+    return;
+  }
+  d_storage.setVal(static_cast<unsigned int>(id), val);
+}
+
+RDGeom::Point3D UniformRealValueGrid3D::getGridPointLoc(
+    unsigned int pointId) const {
+  if (pointId >= d_numX * d_numY * d_numZ) {
+    throw IndexErrorException(pointId);
+  }
+  RDGeom::Point3D res((pointId % d_numX) * d_spacing,
+                      ((pointId % (d_numX * d_numY)) / d_numX) * d_spacing,
+                      (pointId / (d_numX * d_numY)) * d_spacing);
+  res += d_offSet;  // d_origin;
+  return res;
+}
+
+void UniformRealValueGrid3D::setVal(unsigned int pointId, double val) {
+  d_storage.setVal(pointId, val);
+}
+
+bool UniformRealValueGrid3D::compareParams(
+    const UniformRealValueGrid3D &other) const {
+  if (d_numX != other.getNumX()) {
+    return false;
+  }
+  if (d_numY != other.getNumY()) {
+    return false;
+  }
+  if (d_numZ != other.getNumZ()) {
+    return false;
+  }
+  if (fabs(d_spacing - other.getSpacing()) > SPACING_TOL) {
+    return false;
+  }
+  auto dOffset = d_offSet;
+  dOffset -= other.getOffset();
+  return dOffset.lengthSq() <= OFFSET_TOL;
+}
+
+bool UniformRealValueGrid3D::compareVectors(
+    const UniformRealValueGrid3D &other) const {
+  return d_storage.compareVectors(other.d_storage);
+}
+
+bool UniformRealValueGrid3D::compareGrids(
+    const UniformRealValueGrid3D &other) const {
+  if (!compareParams(other)) {
+    return false;
+  }
+  return d_storage.compareVectors(other.d_storage);
+}
+
+UniformRealValueGrid3D &UniformRealValueGrid3D::operator=(
+    const UniformRealValueGrid3D &other) {
+  if (this == &other) {
+    return *this;
+  }
+  d_numX = other.d_numX;
+  d_numY = other.d_numY;
+  d_numZ = other.d_numZ;
+  d_spacing = other.d_spacing;
+  d_offSet = other.d_offSet;
+  d_storage = other.d_storage;
+
+  return *this;
+}
+UniformRealValueGrid3D &UniformRealValueGrid3D::operator|=(
+    const UniformRealValueGrid3D &other) {
+  PRECONDITION(compareParams(other), INCOMPATIBLE_GRIDS);
+
+  d_storage |= other.d_storage;
+  return *this;
+}
+
+UniformRealValueGrid3D &UniformRealValueGrid3D::operator&=(
+    const UniformRealValueGrid3D &other) {
+  PRECONDITION(compareParams(other), INCOMPATIBLE_GRIDS);
+
+  d_storage &= other.d_storage;
+  return *this;
+}
+
+UniformRealValueGrid3D &UniformRealValueGrid3D::operator+=(
+    const UniformRealValueGrid3D &other) {
+  PRECONDITION(compareParams(other), INCOMPATIBLE_GRIDS);
+
+  d_storage += other.d_storage;
+  return *this;
+}
+
+UniformRealValueGrid3D &UniformRealValueGrid3D::operator-=(
+    const UniformRealValueGrid3D &other) {
+  PRECONDITION(compareParams(other), INCOMPATIBLE_GRIDS);
+
+  d_storage -= other.d_storage;
+  return *this;
+}
+
+std::string UniformRealValueGrid3D::toString() const {
+  std::stringstream ss(std::ios_base::binary | std::ios_base::out |
+                       std::ios_base::in);
+  std::int32_t tVers = ci_RealValueGrid3DPICKLE_VERSION * -1;
+  streamWrite(ss, tVers);
+  std::uint32_t tInt;
+  tInt = d_numX;
+  streamWrite(ss, tInt);
+  tInt = d_numY;
+  streamWrite(ss, tInt);
+  tInt = d_numZ;
+  streamWrite(ss, tInt);
+  streamWrite(ss, d_spacing);
+  streamWrite(ss, d_offSet.x);
+  streamWrite(ss, d_offSet.y);
+  streamWrite(ss, d_offSet.z);
+  std::string storePkl = d_storage.toString();
+  std::uint32_t pklSz = storePkl.size();
+  streamWrite(ss, pklSz);
+  ss.write(storePkl.c_str(), pklSz * sizeof(char));
+
+  return ss.str();
+}
+
+void UniformRealValueGrid3D::initFromText(const char *pkl,
+                                          const unsigned int length) {
+  std::stringstream ss(std::ios_base::binary | std::ios_base::in |
+                       std::ios_base::out);
+  ss.write(pkl, length);
+  std::int32_t tVers;
+  streamRead(ss, tVers);
+  tVers = -tVers;
+  if (tVers != 0x1) {
+    throw ValueErrorException("bad version in UniformRealValueGrid3D pickle");
+  }
+  std::uint32_t tInt;
+  streamRead(ss, tInt);
+  d_numX = tInt;
+  streamRead(ss, tInt);
+  d_numY = tInt;
+  streamRead(ss, tInt);
+  d_numZ = tInt;
+  streamRead(ss, d_spacing);
+  double oX, oY, oZ;
+  streamRead(ss, oX);
+  streamRead(ss, oY);
+  streamRead(ss, oZ);
+  d_offSet = RDGeom::Point3D(oX, oY, oZ);
+
+  std::uint32_t pklSz;
+  streamRead(ss, pklSz);
+  std::vector<char> buff(pklSz);
+  ss.read(buff.data(), pklSz * sizeof(char));
+  d_storage = RDKit::RealValueVect(buff.data(), pklSz);
+}
+
+UniformRealValueGrid3D operator&(const UniformRealValueGrid3D &grd1,
+                                 const UniformRealValueGrid3D &grd2) {
+  UniformRealValueGrid3D ans(grd1);
+  ans &= grd2;
+  return ans;
+};
+
+UniformRealValueGrid3D operator|(const UniformRealValueGrid3D &grd1,
+                                 const UniformRealValueGrid3D &grd2) {
+  UniformRealValueGrid3D ans(grd1);
+  ans |= grd2;
+  return ans;
+};
+
+UniformRealValueGrid3D operator+(const UniformRealValueGrid3D &grd1,
+                                 const UniformRealValueGrid3D &grd2) {
+  UniformRealValueGrid3D ans(grd1);
+  ans += grd2;
+  return ans;
+};
+
+UniformRealValueGrid3D operator-(const UniformRealValueGrid3D &grd1,
+                                 const UniformRealValueGrid3D &grd2) {
+  UniformRealValueGrid3D ans(grd1);
+  ans -= grd2;
+  return ans;
+};
+
+}  // namespace RDGeom

Code/Geometry/UniformRealValueGrid3D.h

@@ -0,0 +1,201 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <RDGeneral/export.h>
+#ifndef UNIFORMREALVALUEGRID3D_H_20140403
+#define UNIFORMREALVALUEGRID3D_H_20140403
+
+#include <DataStructs/RealValueVect.h>
+#include <Geometry/point.h>
+#include "Grid3D.h"
+
+namespace RDGeom {
+class RDKIT_RDGEOMETRYLIB_EXPORT UniformRealValueGrid3D
+    : public Grid3D<RDKit::RealValueVect, double, double> {
+ public:
+  UniformRealValueGrid3D() { initGrid(1.0, 1.0, 1.0, 1.0, RDGeom::Point3D()); };
+  //! \brief ctor
+  /*
+      \param dimX:    the x dimension of the grid, in Angstroms
+      \param dimY:    the y dimension of the grid, in Angstroms
+      \param dimZ:    the z dimension of the grid, in Angstroms
+      \param spacing: the grid spacing, in Angstroms
+      \param offset:  OPTIONAL: the offset of the grid from (0,0,0), in
+     Angstroms.
+
+      \b Note: the values of arguments such as \c dimX and \c spacing
+      don't actually need to be in Angstroms, but they should be internally
+      consistent.
+
+   */
+  UniformRealValueGrid3D(double dimX, double dimY, double dimZ,
+                         double spacing = 0.5,
+                         const RDGeom::Point3D *offset = nullptr,
+                         const RDKit::RealValueVect *data = nullptr) {
+    if (offset == nullptr) {
+      initGrid(dimX, dimY, dimZ, spacing,
+               RDGeom::Point3D(-0.5 * dimX, -0.5 * dimY, -0.5 * dimZ), data);
+    } else {
+      initGrid(dimX, dimY, dimZ, spacing, *offset, data);
+    }
+  }
+  UniformRealValueGrid3D(const UniformRealValueGrid3D &other);
+  UniformRealValueGrid3D &operator=(const UniformRealValueGrid3D &other);
+
+  //! construct from a string pickle
+  UniformRealValueGrid3D(const std::string &pkl);
+  //! construct from a text pickle
+  UniformRealValueGrid3D(const char *pkl, unsigned int);
+
+  //! \brief Get the index of the grid point closest to point
+  //!
+  //! \return the integer index, -1 if the specified point is outside the grid
+  int getGridPointIndex(const RDGeom::Point3D &point) const override;
+
+  //! \brief Get the value at the grid point closest to the specified point
+  //!
+  //! \return the double value, -1 if the specified index is outside the grid
+  double getVal(const RDGeom::Point3D &point) const override;
+
+  //! \brief Get the value at a specified grid point
+  //!
+  //! \return the double value
+  double getVal(unsigned int pointId) const override;
+
+  //! \brief Set the value at the grid point closest to the specified point
+  //!
+  //! doesn't do anything if the point is outside the grid
+  void setVal(const RDGeom::Point3D &point, double val) override;
+
+  //! \brief get the location of the specified grid point
+  RDGeom::Point3D getGridPointLoc(unsigned int pointId) const override;
+
+  //! \brief Set the value at the specified grid point
+  void setVal(unsigned int pointId, double val) override;
+
+  //! \brief get the size of the grid (number of grid points)
+  unsigned int getSize() const override { return d_numX * d_numY * d_numZ; };
+
+  //! \brief get the index of the grid point given the x, y, z indices
+  //!
+  //! \return the integer value, -1 if the indices are outside the grid
+  int getGridIndex(unsigned int xi, unsigned int yi, unsigned int zi) const;
+
+  //! \brief get the x, y, and z indices of a grid-point index
+  //!
+  void getGridIndices(unsigned int idx, unsigned int &xi, unsigned int &yi,
+                      unsigned int &zi) const;
+
+  //! \brief get the number of grid points along x-axis
+  unsigned int getNumX() const { return d_numX; };
+
+  //! \brief get the number of grid points along y-axis
+  unsigned int getNumY() const { return d_numY; };
+
+  //! \brief get the number of grid points along z-axis
+  unsigned int getNumZ() const { return d_numZ; };
+
+  //! \brief get the grid's offset
+  const RDGeom::Point3D &getOffset() const { return d_offSet; };
+
+  //! \brief get the grid's spacing
+  double getSpacing() const { return d_spacing; };
+
+  //! \brief return a \b const pointer to our occupancy vector
+  const RDKit::RealValueVect *getOccupancyVect() const override {
+    return &d_storage;
+  };
+
+  //! brief returns raw vector
+  const std::vector<double> &getData() const { return d_storage.getData(); }
+  std::vector<double> &getData() { return d_storage.getData(); }
+
+  //! \brief returns true if the grid \c other has parameters
+  //!        compatible with ours.
+  bool compareParams(const UniformRealValueGrid3D &other) const;
+
+  //! \brief returns true if the grid \c other has the same values
+  //!        as ours.
+  bool compareVectors(const UniformRealValueGrid3D &other) const;
+
+  //! \brief returns true if the grid \c other has the same parameters and
+  //!        values as ours.
+  bool compareGrids(const UniformRealValueGrid3D &other) const;
+
+  //! \brief calculates the union between the data on this grid and
+  //!  that on \c other.
+  //!  This grid is modified.
+  //!  NOTE that the grids must have the same parameters.
+  UniformRealValueGrid3D &operator|=(const UniformRealValueGrid3D &other);
+  //! \brief calculates the intersection between the data on this grid and
+  //!  that on \c other.
+  //!  This grid is modified.
+  //!  NOTE that the grids must have the same parameters.
+  UniformRealValueGrid3D &operator&=(const UniformRealValueGrid3D &other);
+  //! \brief calculates the sum of the data on this grid and
+  //!  that on \c other.
+  //!  This grid is modified.
+  //!  NOTE that the grids must have the same parameters.
+  UniformRealValueGrid3D &operator+=(const UniformRealValueGrid3D &other);
+  //! \brief calculates the difference between the data on this grid and
+  //!  that on \c other.
+  //!  This grid is modified.
+  //!  NOTE that the grids must have the same parameters.
+  UniformRealValueGrid3D &operator-=(const UniformRealValueGrid3D &other);
+
+  //! \brief create and return a pickle
+  std::string toString() const;
+
+  /*
+  UniformRealValueGrid3D operator& (const UniformRealValueGrid3D &other) const{
+    PRECONDITION(dp_storage,"bad storage");
+    PRECONDITION(compareParams(other),"mismatched params");
+    UniformRealValueGrid3D
+  res(d_numX*d_spacing,d_numY*d_spacing,d_numZ*d_spacing, d_spacing,&d_offSet);
+    return res;
+  };*/
+
+ private:
+  //! \brief internal initialization code
+  /*
+      \param dimX:    the x dimension of the grid, in Angstroms
+      \param dimY:    the y dimension of the grid, in Angstroms
+      \param dimZ:    the z dimension of the grid, in Angstroms
+      \param spacing: the grid spacing, in Angstroms
+      \param offset:  the offset of the grid from (0,0,0), in Angstroms.
+      \param data:    (optional) a pointer to a DoubleVector to use, we take
+                      ownership of the pointer.
+   */
+  void initGrid(double dimX, double dimY, double dimZ, double spacing,
+                const RDGeom::Point3D &offSet,
+                const RDKit::RealValueVect *data = nullptr);
+
+  unsigned int d_numX, d_numY,
+      d_numZ;        //! number of grid points along x, y, z axes
+  double d_spacing;  //! grid spacing
+
+  RDGeom::Point3D d_offSet;        //! the grid offset (from the origin)
+  RDKit::RealValueVect d_storage;  //! storage for values at each grid point
+
+  //! \brief construct from a pickle
+  void initFromText(const char *pkl, const unsigned int length);
+};
+
+RDKIT_RDGEOMETRYLIB_EXPORT UniformRealValueGrid3D operator|(
+    const UniformRealValueGrid3D &grd1, const UniformRealValueGrid3D &grd2);
+RDKIT_RDGEOMETRYLIB_EXPORT UniformRealValueGrid3D operator&(
+    const UniformRealValueGrid3D &grd1, const UniformRealValueGrid3D &grd2);
+RDKIT_RDGEOMETRYLIB_EXPORT UniformRealValueGrid3D operator+(
+    const UniformRealValueGrid3D &grd1, const UniformRealValueGrid3D &grd2);
+RDKIT_RDGEOMETRYLIB_EXPORT UniformRealValueGrid3D operator-(
+    const UniformRealValueGrid3D &grd1, const UniformRealValueGrid3D &grd2);
+}  // namespace RDGeom
+
+#endif

Code/Geometry/Wrap/CMakeLists.txt

@@ -1,6 +1,6 @@
 remove_definitions(-DRDKIT_RDGEOMETRYLIB_BUILD)
 rdkit_python_extension(rdGeometry
-                       Point.cpp UniformGrid3D.cpp rdGeometry.cpp
+                       Point.cpp UniformGrid3D.cpp UniformRealValueGrid3D.cpp rdGeometry.cpp
                        DEST Geometry
                        LINK_LIBRARIES RDGeometryLib)
 

Code/Geometry/Wrap/UniformRealValueGrid3D.cpp

@@ -0,0 +1,128 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <boost/python.hpp>
+#include <RDBoost/Wrap.h>
+#include <RDGeneral/types.h>
+#include <RDGeneral/Invariant.h>
+#include <DataStructs/RealValueVect.h>
+#include <Geometry/point.h>
+#include <Geometry/UniformRealValueGrid3D.h>
+#include <Geometry/GridUtils.h>
+namespace python = boost::python;
+
+using namespace RDKit;
+
+namespace RDGeom {
+struct urvg3d_pickle_suite : rdkit_pickle_suite {
+  static python::tuple getinitargs(const UniformRealValueGrid3D &self) {
+    auto res = self.toString();
+    python::object retval(
+        python::handle<>(PyBytes_FromStringAndSize(res.c_str(), res.length())));
+    return python::make_tuple(retval);
+  }
+};
+
+UniformRealValueGrid3D *makeUniformRealValueGrid3D(
+    double dimX, double dimY, double dimZ, double spacing = 0.5,
+    const Point3D *offSet = nullptr) {
+  UniformRealValueGrid3D *grd =
+      new UniformRealValueGrid3D(dimX, dimY, dimZ, spacing, offSet);
+  return grd;
+}
+
+double getValPoint(const UniformRealValueGrid3D &grid, const Point3D &pt) {
+  return grid.getVal(pt);
+}
+
+double getValIndex(const UniformRealValueGrid3D &grid, unsigned int id) {
+  return grid.getVal(id);
+}
+
+void setValIndex(UniformRealValueGrid3D &grid, unsigned int id, double val) {
+  grid.setVal(id, val);
+}
+
+void setValPoint(UniformRealValueGrid3D &grid, const Point3D &pt, double val) {
+  grid.setVal(pt, val);
+}
+
+python::tuple getGridIndicesWrap(const UniformRealValueGrid3D &grid,
+                                 unsigned int idx) {
+  unsigned int xi, yi, zi;
+  grid.getGridIndices(idx, xi, yi, zi);
+  return python::make_tuple(xi, yi, zi);
+}
+
+std::string urvGridClassDoc =
+    "Class to represent a uniform three-dimensional\n\
+    cubic grid. Each grid point can store a floating point value. \n";
+
+struct urvGrid3D_wrapper {
+  static void wrap() {
+    python::class_<UniformRealValueGrid3D>(
+        "UniformRealValueGrid3D", urvGridClassDoc.c_str(),
+        python::init<std::string>("Pickle constructor"))
+        .def(python::init<>("Default constructor"))
+        .def(python::init<RDGeom::UniformRealValueGrid3D>("Copy constructor"))
+        .def("__init__",
+             python::make_constructor(
+                 makeUniformRealValueGrid3D, python::default_call_policies(),
+                 (python::arg("dimX"), python::arg("dimY"), python::arg("dimZ"),
+                  python::arg("spacing") = 0.5,
+                  python::arg("offSet") = python::object())),
+             "Constructor")
+        .def("GetGridPointIndex", &UniformRealValueGrid3D::getGridPointIndex,
+             "Get the index to the grid point closest to the specified point")
+        .def(
+            "GetGridIndex", &UniformRealValueGrid3D::getGridIndex,
+            "Get the index to the grid point with the three integer indices provided")
+        .def("GetGridIndices", &getGridIndicesWrap,
+             "Returns the integer indices of the grid index provided.")
+        .def("GetValPoint", getValPoint,
+             "Get the value at the closest grid point")
+        .def("GetVal", getValIndex, "Get the value at the specified grid point")
+        .def("SetVal", setValIndex, "Set the value at the specified grid point")
+        .def("SetValPoint", setValPoint,
+             "Set the value at grid point closest to the specified point")
+        .def("GetGridPointLoc", &UniformRealValueGrid3D::getGridPointLoc,
+             "Get the location of the specified grid point")
+        .def("GetSize", &UniformRealValueGrid3D::getSize,
+             "Get the size of the grid (number of grid points)")
+        .def("GetNumX", &UniformRealValueGrid3D::getNumX,
+             "Get the number of grid points along x-axis")
+        .def("GetNumY", &UniformRealValueGrid3D::getNumY,
+             "Get the number of grid points along y-axis")
+        .def("GetNumZ", &UniformRealValueGrid3D::getNumZ,
+             "Get the number of grid points along z-axis")
+        .def("GetOffset", &UniformRealValueGrid3D::getOffset,
+             python::return_value_policy<python::copy_const_reference>(),
+             "Get the location of the center of the grid")
+        .def("GetSpacing", &UniformRealValueGrid3D::getSpacing,
+             "Get the grid spacing")
+        .def("GetOccupancyVect", &UniformRealValueGrid3D::getOccupancyVect,
+             python::return_value_policy<python::reference_existing_object>(),
+             "Get the occupancy vector for the grid")
+        .def("CompareVectors", &UniformRealValueGrid3D::compareVectors,
+             "Compare the vector values between two grid objects.")
+        .def("CompareParams", &UniformRealValueGrid3D::compareParams,
+             "Compare the parameters between two grid object.")
+        .def("CompareGrids", &UniformRealValueGrid3D::compareGrids,
+             "Compare the parameters and values between two grid objects.")
+        .def(python::self &= python::self)
+        .def(python::self |= python::self)
+        .def(python::self += python::self)
+        .def(python::self -= python::self)
+        .def_pickle(RDGeom::urvg3d_pickle_suite());
+  }
+};
+}  // namespace RDGeom
+
+void wrap_uniformrealvalueGrid() { RDGeom::urvGrid3D_wrapper::wrap(); }

Code/Geometry/Wrap/rdGeometry.cpp

@@ -14,6 +14,7 @@
 namespace python = boost::python;
 void wrap_point();
 void wrap_uniformGrid();
+void wrap_uniformrealvalueGrid();
 
 BOOST_PYTHON_MODULE(rdGeometry) {
   python::scope().attr("__doc__") =
@@ -21,4 +22,5 @@ BOOST_PYTHON_MODULE(rdGeometry) {
 
   wrap_point();
   wrap_uniformGrid();
+  wrap_uniformrealvalueGrid();
 }

Code/Geometry/Wrap/testGeometry.py

@@ -406,6 +406,405 @@ class TestCase(unittest.TestCase):
     self.assertEqual(p2.x, p3.x)
     self.assertEqual(p2.y, p3.y)
 
+  def test1bPoint2D(self):
+    pt = geom.Point2D()
+    self.assertTrue(feq(pt.x, 0.0))
+    self.assertTrue(feq(pt.y, 0.0))
+
+    pt = geom.Point2D(3., 4.)
+    self.assertTrue(feq(pt.x, 3.0))
+    self.assertTrue(feq(pt.y, 4.0))
+
+    self.assertTrue(feq(pt.x, 3.0))
+    self.assertTrue(feq(pt.y, 4.0))
+    self.assertTrue(feq(pt[0], 3.0))
+    self.assertTrue(feq(pt[1], 4.0))
+    self.assertTrue(feq(pt[-2], 3.0))
+    self.assertTrue(feq(pt[-1], 4.0))
+    lst = list(pt)
+    self.assertTrue(feq(lst[0], 3.0))
+    self.assertTrue(feq(lst[1], 4.0))
+
+    pt2 = geom.Point2D(1., 1.)
+
+    pt3 = pt + pt2
+    self.assertTrue(feq(pt3.x, 4.0))
+    self.assertTrue(feq(pt3.y, 5.0))
+
+    pt += pt2
+    self.assertTrue(feq(pt.x, 4.0))
+    self.assertTrue(feq(pt.y, 5.0))
+
+    pt3 = pt - pt2
+    self.assertTrue(feq(pt3.x, 3.0))
+    self.assertTrue(feq(pt3.y, 4.0))
+
+    pt -= pt2
+    self.assertTrue(feq(pt.x, 3.0))
+    self.assertTrue(feq(pt.y, 4.0))
+
+    pt *= 2.0
+    self.assertTrue(feq(pt.x, 6.0))
+    self.assertTrue(feq(pt.y, 8.0))
+
+    pt /= 2
+    self.assertTrue(feq(pt.x, 3.0))
+    self.assertTrue(feq(pt.y, 4.0))
+    self.assertTrue(feq(pt.Length(), 5.0))
+    self.assertTrue(feq(pt.LengthSq(), 25.0))
+    pt.Normalize()
+    self.assertTrue(feq(pt.Length(), 1.0))
+
+    pt1 = geom.Point2D(1.0, 0.0)
+    pt2 = geom.Point2D(2.0 * math.cos(math.pi / 6), 2.0 * math.sin(math.pi / 6))
+    ang = pt1.AngleTo(pt2)
+    self.assertTrue(feq(ang, math.pi / 6))
+
+    prod = pt1.DotProduct(pt2)
+    self.assertTrue(feq(prod, 2.0 * math.cos(math.pi / 6)))
+
+  def test1cPointND(self):
+    dim = 4
+    pt = geom.PointND(4)
+    for i in range(dim):
+      self.assertTrue(feq(pt[i], 0.0))
+
+    pt[0] = 3
+    pt[3] = 4
+    self.assertTrue(feq(pt[0], 3.0))
+    self.assertTrue(feq(pt[3], 4.0))
+    self.assertTrue(feq(pt[-4], 3.0))
+    self.assertTrue(feq(pt[-1], 4.0))
+    lst = list(pt)
+    self.assertTrue(feq(lst[0], 3.0))
+    self.assertTrue(feq(lst[3], 4.0))
+
+    pt2 = geom.PointND(4)
+    pt2[0] = 1.
+    pt2[2] = 1.
+
+    pt3 = pt + pt2
+    self.assertTrue(feq(pt3[0], 4.0))
+    self.assertTrue(feq(pt3[2], 1.0))
+    self.assertTrue(feq(pt3[3], 4.0))
+
+    pt += pt2
+    self.assertTrue(feq(pt[0], 4.0))
+    self.assertTrue(feq(pt[2], 1.0))
+    self.assertTrue(feq(pt[3], 4.0))
+
+    pt3 = pt - pt2
+    self.assertTrue(feq(pt3[0], 3.0))
+    self.assertTrue(feq(pt3[2], 0.0))
+    self.assertTrue(feq(pt3[3], 4.0))
+
+    pt -= pt2
+    self.assertTrue(feq(pt[0], 3.0))
+    self.assertTrue(feq(pt[2], 0.0))
+    self.assertTrue(feq(pt[3], 4.0))
+
+    pt *= 2.0
+    self.assertTrue(feq(pt[0], 6.0))
+    self.assertTrue(feq(pt[1], 0.0))
+    self.assertTrue(feq(pt[2], 0.0))
+    self.assertTrue(feq(pt[3], 8.0))
+
+    pt /= 2
+    self.assertTrue(feq(pt[0], 3.0))
+    self.assertTrue(feq(pt[1], 0.0))
+    self.assertTrue(feq(pt[2], 0.0))
+    self.assertTrue(feq(pt[3], 4.0))
+
+    self.assertTrue(feq(pt.Length(), 5.0))
+    self.assertTrue(feq(pt.LengthSq(), 25.0))
+    pt.Normalize()
+    self.assertTrue(feq(pt.Length(), 1.0))
+
+    pkl = pickle.dumps(pt)
+    pt2 = pickle.loads(pkl)
+    self.assertTrue(len(pt) == len(pt2))
+    for i in range(len(pt)):
+      self.assertTrue(feq(pt2[i], pt[i]))
+
+  def test3UniformGrid(self):
+    ugrid = geom.UniformGrid3D(20, 18, 15)
+    self.assertTrue(ugrid.GetNumX() == 40)
+    self.assertTrue(ugrid.GetNumY() == 36)
+    self.assertTrue(ugrid.GetNumZ() == 30)
+    dvect = ugrid.GetOccupancyVect()
+    ugrid = geom.UniformGrid3D(20, 18, 15, 0.5, DataStructs.DiscreteValueType.TWOBITVALUE)
+    dvect = ugrid.GetOccupancyVect()
+    self.assertTrue(dvect.GetValueType() == DataStructs.DiscreteValueType.TWOBITVALUE)
+
+    grd = geom.UniformGrid3D(10.0, 10.0, 10.0, 0.5)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(2.0, 2.0, 0.0), 1.5, 0.25)
+
+    geom.WriteGridToFile(grd, "junk.grd")
+    grd2 = geom.UniformGrid3D(10.0, 10.0, 10.0, 0.5)
+    grd2.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.5, 0.25)
+    grd2.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.5, 0.25)
+    grd2.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.5, 0.25)
+
+    dist = geom.TanimotoDistance(grd, grd2)
+    self.assertTrue(dist == 0.25)
+    dist = geom.ProtrudeDistance(grd, grd2)
+    self.assertTrue(dist == 0.25)
+    dist = geom.ProtrudeDistance(grd2, grd)
+    self.assertTrue(dist == 0.0)
+
+    grd2 = geom.UniformGrid3D(10.0, 10.0, 10.0, 0.5, DataStructs.DiscreteValueType.FOURBITVALUE)
+    grd2.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.5, 0.25, 3)
+    grd2.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.5, 0.25, 3)
+    grd2.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.5, 0.25, 3)
+    self.assertRaises(ValueError, lambda: geom.TanimotoDistance(grd, grd2))
+
+    grd2 = geom.UniformGrid3D(10.0, 10.0, 10.0, 1.0)
+    self.assertRaises(ValueError, lambda: geom.TanimotoDistance(grd, grd2))
+
+    grd2 = geom.UniformGrid3D(11.0, 10.0, 10.0, 1.0)
+    self.assertRaises(ValueError, lambda: geom.TanimotoDistance(grd, grd2))
+
+  def testSymmetry(self):
+    grd = geom.UniformGrid3D(10.0, 10.0, 10.0, 0.5)
+    grd.SetSphereOccupancy(geom.Point3D(-2.2, -2.0, 0.0), 1.65, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(2.2, -2.0, 0.0), 1.65, 0.25)
+
+    bPt1 = geom.Point3D(-4.0, -2.0, -2.0)
+    bPt2 = geom.Point3D(4.0, -2.0, -2.0)
+    for k in range(8):
+      bPt1 += geom.Point3D(0.0, 0.0, 0.5)
+      bPt2 += geom.Point3D(0.0, 0.0, 0.5)
+      for j in range(8):
+        bPt1 += geom.Point3D(0.0, 0.5, 0.0)
+        bPt2 += geom.Point3D(0.0, 0.5, 0.0)
+        for i in range(8):
+          bPt1 += geom.Point3D(0.5, 0.0, 0.0)
+          bPt2 -= geom.Point3D(0.5, 0.0, 0.0)
+          self.assertTrue(grd.GetValPoint(bPt1) == grd.GetValPoint(bPt2))
+
+        bPt1.x = -4.0
+        bPt2.x = 4.0
+      bPt1.y = -2.0
+      bPt2.y = -2.0
+
+  def test4UniformRealValueGrid(self):
+    ugrid = geom.UniformRealValueGrid3D(20, 18, 15)
+    self.assertTrue(ugrid.GetNumX() == 40)
+    self.assertTrue(ugrid.GetNumY() == 36)
+    self.assertTrue(ugrid.GetNumZ() == 30)
+    dvect = ugrid.GetOccupancyVect()
+    ugrid = geom.UniformRealValueGrid3D(20, 18, 15, 0.5)
+    ugrid.SetVal(50, 2.3)
+    val = ugrid.GetVal(50)
+    self.assertTrue(feq(val, 2.3))
+
+  def test5PointPickles(self):
+    pt = geom.Point3D(2.0, -3.0, 1.0)
+    pt2 = pickle.loads(pickle.dumps(pt))
+    self.assertTrue(feq(pt.x, pt2.x, 1e-6))
+    self.assertTrue(feq(pt.y, pt2.y, 1e-6))
+    self.assertTrue(feq(pt.z, pt2.z, 1e-6))
+
+    pt = geom.Point2D(2.0, -4.0)
+    pt2 = pickle.loads(pickle.dumps(pt))
+    self.assertTrue(feq(pt.x, pt2.x, 1e-6))
+    self.assertTrue(feq(pt.y, pt2.y, 1e-6))
+
+  def test6GridPickles(self):
+    grd = geom.UniformGrid3D(10.0, 9.0, 8.0, 0.5)
+    self.assertTrue(grd.GetNumX() == 20)
+    self.assertTrue(grd.GetNumY() == 18)
+    self.assertTrue(grd.GetNumZ() == 16)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.5, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(2.0, 2.0, 0.0), 1.5, 0.25)
+
+    self.assertTrue(geom.TanimotoDistance(grd, grd) == 0.0)
+
+    grd2 = pickle.loads(pickle.dumps(grd))
+    self.assertTrue(grd2.GetNumX() == 20)
+    self.assertTrue(grd2.GetNumY() == 18)
+    self.assertTrue(grd2.GetNumZ() == 16)
+    self.assertTrue(geom.TanimotoDistance(grd, grd2) == 0.0)
+
+  def test7RealGridPickles(self):
+    grd = geom.UniformRealValueGrid3D(10.0, 9.0, 8.0, 0.5)
+    self.assertTrue(grd.GetNumX() == 20)
+    self.assertTrue(grd.GetNumY() == 18)
+    self.assertTrue(grd.GetNumZ() == 16)
+
+    grd2 = pickle.loads(pickle.dumps(grd))
+    self.assertTrue(grd2.GetNumX() == 20)
+    self.assertTrue(grd2.GetNumY() == 18)
+    self.assertTrue(grd2.GetNumZ() == 16)
+
+  def test8GridOps(self):
+    grd = geom.UniformGrid3D(10, 10, 10)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.0, 0.25)
+    grd.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.0, 0.25)
+
+    grd2 = geom.UniformGrid3D(10, 10, 10)
+    grd2.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.0, 0.25)
+    grd2.SetSphereOccupancy(geom.Point3D(2.0, 2.0, 0.0), 1.0, 0.25)
+
+    self.assertTrue(geom.TanimotoDistance(grd, grd) == 0.0)
+    self.assertTrue(geom.TanimotoDistance(grd, grd2) == 1.0)
+
+    grd3 = copy.deepcopy(grd)
+    grd3 |= grd2
+    self.assertTrue(geom.TanimotoDistance(grd3, grd) == .5)
+    self.assertTrue(geom.TanimotoDistance(grd3, grd2) == .5)
+
+    grd3 = copy.deepcopy(grd)
+    grd3 += grd2
+    self.assertTrue(geom.TanimotoDistance(grd3, grd) == .5)
+    self.assertTrue(geom.TanimotoDistance(grd3, grd2) == .5)
+
+    grd3 -= grd
+    self.assertTrue(geom.TanimotoDistance(grd3, grd) == 1.0)
+    self.assertTrue(geom.TanimotoDistance(grd3, grd2) == 0)
+
+    grd4 = geom.UniformGrid3D(10, 10, 10)
+    grd4.SetSphereOccupancy(geom.Point3D(-2.0, -2.0, 0.0), 1.0, 0.25)
+    grd4.SetSphereOccupancy(geom.Point3D(-2.0, 2.0, 0.0), 1.0, 0.25)
+    grd4.SetSphereOccupancy(geom.Point3D(2.0, -2.0, 0.0), 1.0, 0.25)
+    self.assertTrue(feq(geom.TanimotoDistance(grd4, grd), .3333))
+    self.assertTrue(feq(geom.TanimotoDistance(grd4, grd2), .75))
+
+    grd4 &= grd2
+    self.assertTrue(feq(geom.TanimotoDistance(grd4, grd), 1.0))
+    self.assertTrue(feq(geom.TanimotoDistance(grd4, grd2), .5))
+
+  def test9RealGridOps(self):
+    grd1 = geom.UniformRealValueGrid3D(5.0, 5.0, 5.0, 0.1)
+    grd2 = geom.UniformRealValueGrid3D(5.0, 5.0, 5.0, 0.1)
+
+    grd1.SetVal(50, 37.37)
+    grd2.SetVal(50, 1.03)
+    grd3 = copy.deepcopy(grd2)
+    grd4 = geom.UniformRealValueGrid3D(grd2)
+
+    grd1 |= grd2
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd2.GetVal(50), 1.03))
+
+    grd2 |= grd1
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd2.GetVal(50), 37.37))
+
+    grd2 &= grd3
+    self.assertTrue(feq(grd2.GetVal(50), 1.03))
+    self.assertTrue(feq(grd3.GetVal(50), 1.03))
+
+    grd3 &= grd1
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd3.GetVal(50), 1.03))
+
+    grd2 &= grd4
+    self.assertTrue(feq(grd2.GetVal(50), 1.03))
+    self.assertTrue(feq(grd4.GetVal(50), 1.03))
+
+    grd4 &= grd1
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd4.GetVal(50), 1.03))
+
+    grd1 += grd2
+    self.assertTrue(feq(grd1.GetVal(50), 38.40))
+    self.assertTrue(feq(grd2.GetVal(50), 1.03))
+
+    grd1 -= grd2
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd2.GetVal(50), 1.03))
+
+    grd2 -= grd1
+    self.assertTrue(feq(grd1.GetVal(50), 37.37))
+    self.assertTrue(feq(grd2.GetVal(50), -36.34))
+
+  def test10Dihedrals(self):
+    p1 = geom.Point3D(1, 0, 0)
+    p2 = geom.Point3D(0, 0, 0)
+    p3 = geom.Point3D(0, 1, 0)
+
+    p4 = geom.Point3D(.5, 1, .5)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 4, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, -math.pi / 4, 4)
+
+    p4 = geom.Point3D(-.5, 1, .5)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, 3 * math.pi / 4, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, -3 * math.pi / 4, 4)
+
+    p4 = geom.Point3D(.5, 1, -.5)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 4, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 4, 4)
+
+    p4 = geom.Point3D(-.5, 1, -.5)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, 3 * math.pi / 4, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, 3 * math.pi / 4, 4)
+
+    p4 = geom.Point3D(0, 1, 1)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 2, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, -math.pi / 2, 4)
+
+    p4 = geom.Point3D(0, 1, -1)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 2, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi / 2, 4)
+
+    p4 = geom.Point3D(1, 1, 0)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, 0, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, 0, 4)
+
+    p4 = geom.Point3D(-1, 1, 0)
+    ang = geom.ComputeDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi, 4)
+    ang = geom.ComputeSignedDihedralAngle(p1, p2, p3, p4)
+    self.assertAlmostEqual(ang, math.pi, 4)
+
+  def test11UniformGridIndices(self):
+    ugrid = geom.UniformGrid3D(20, 18, 15)
+    idx = ugrid.GetGridIndex(3, 2, 1)
+    xi, yi, zi = ugrid.GetGridIndices(idx)
+    self.assertEqual(xi, 3)
+    self.assertEqual(yi, 2)
+    self.assertEqual(zi, 1)
+
+  def test12UniformRealGridIndices(self):
+    ugrid = geom.UniformRealValueGrid3D(20, 18, 15)
+    idx = ugrid.GetGridIndex(3, 2, 1)
+    xi, yi, zi = ugrid.GetGridIndices(idx)
+    self.assertEqual(xi, 3)
+    self.assertEqual(yi, 2)
+    self.assertEqual(zi, 1)
+
+    pt = ugrid.GetOffset()
+    ugrid.SetValPoint(pt, 2.3)
+    idx = ugrid.GetGridPointIndex(pt)
+    self.assertTrue(feq(ugrid.GetValPoint(pt), 2.3))
+    self.assertTrue(idx == 0)
+    self.assertTrue(feq(ugrid.GetVal(idx), 2.3))
+
+    pt2 = ugrid.GetGridPointLoc(idx)
+    self.assertTrue(feq(pt.x, pt2.x))
+    self.assertTrue(feq(pt.y, pt2.y))
+    self.assertTrue(feq(pt.z, pt2.z))
+
 
 if __name__ == '__main__':
   print("Testing Geometry wrapper")

Code/Geometry/testRealValueGrid.cpp

@@ -0,0 +1,185 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 "UniformRealValueGrid3D.h"
+#include <GraphMol/GraphMol.h>
+#include <Geometry/point.h>
+#include <RDGeneral/Invariant.h>
+#include <RDGeneral/utils.h>
+#include <RDGeneral/RDLog.h>
+#include <RDGeneral/types.h>
+#include <algorithm>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <ios>
+#include <stdlib.h>
+
+using namespace RDGeom;
+using namespace RDKit;
+
+void test1UniformRealValueGrid3D() {
+  RDGeom::UniformRealValueGrid3D grd(6.0, 5.0, 4.0);
+  CHECK_INVARIANT(grd.getSize() == 960, "");
+  CHECK_INVARIANT(RDKit::feq(grd.getSpacing(), .5), "");
+  CHECK_INVARIANT(grd.getNumX() == 12, "");
+  CHECK_INVARIANT(grd.getNumY() == 10, "");
+  CHECK_INVARIANT(grd.getNumZ() == 8, "");
+
+  RDGeom::UniformRealValueGrid3D grd2(grd);
+  CHECK_INVARIANT(grd2.getSize() == 960, "");
+  CHECK_INVARIANT(RDKit::feq(grd2.getSpacing(), .5), "");
+  CHECK_INVARIANT(grd2.getNumX() == 12, "");
+  CHECK_INVARIANT(grd2.getNumY() == 10, "");
+  CHECK_INVARIANT(grd2.getNumZ() == 8, "");
+  CHECK_INVARIANT(grd2.getOccupancyVect()->getTotalVal() == 0, "");
+  TEST_ASSERT(grd.compareVectors(grd2));
+  TEST_ASSERT(grd.compareParams(grd2));
+
+  grd.setVal(1, 1.0);
+  TEST_ASSERT(!grd.compareVectors(grd2));
+  TEST_ASSERT(grd.compareParams(grd2))
+  // make sure the data are actually decoupled:
+  // grd.setSphereOccupancy(Point3D(1.0, 1.0, 0.0), 1.5, 0.25);
+  // CHECK_INVARIANT(grd.getOccupancyVect()->getTotalVal()>523, "" );
+  // CHECK_INVARIANT(grd2.getOccupancyVect()->getTotalVal() == 523, "" );
+}
+
+void test2UniformRealValueGrid3DPickling() {
+  RDGeom::UniformRealValueGrid3D grd(5.0, 5.0, 5.0, 0.1);
+  grd.setVal(50, 2.3);
+  const std::string pkl = grd.toString();
+
+  RDGeom::UniformRealValueGrid3D grd2(pkl);
+  TEST_ASSERT(grd.compareVectors(grd2));
+  TEST_ASSERT(grd.compareParams(grd2));
+}
+
+void test3UniformRealValueGrid3DIndexing() {
+  RDGeom::UniformRealValueGrid3D grd(5.0, 5.0, 5.0, 0.1);
+  {
+    unsigned int xi = 3, yi = 3, zi = 3;
+    unsigned int idx = grd.getGridIndex(xi, yi, zi);
+    unsigned int nxi, nyi, nzi;
+    grd.getGridIndices(idx, nxi, nyi, nzi);
+    TEST_ASSERT(nxi == xi);
+    TEST_ASSERT(nyi == yi);
+    TEST_ASSERT(nzi == zi);
+  }
+  {
+    unsigned int xi = 3, yi = 3, zi = 5;
+    unsigned int idx = grd.getGridIndex(xi, yi, zi);
+    unsigned int nxi, nyi, nzi;
+    grd.getGridIndices(idx, nxi, nyi, nzi);
+    TEST_ASSERT(nxi == xi);
+    TEST_ASSERT(nyi == yi);
+    TEST_ASSERT(nzi == zi);
+  }
+  {
+    unsigned int xi = 3, yi = 6, zi = 3;
+    unsigned int idx = grd.getGridIndex(xi, yi, zi);
+    unsigned int nxi, nyi, nzi;
+    grd.getGridIndices(idx, nxi, nyi, nzi);
+    TEST_ASSERT(nxi == xi);
+    TEST_ASSERT(nyi == yi);
+    TEST_ASSERT(nzi == zi);
+  }
+  {
+    unsigned int xi = 0, yi = 0, zi = 0;
+    unsigned int idx = grd.getGridIndex(xi, yi, zi);
+    unsigned int nxi, nyi, nzi;
+    grd.getGridIndices(idx, nxi, nyi, nzi);
+    TEST_ASSERT(nxi == xi);
+    TEST_ASSERT(nyi == yi);
+    TEST_ASSERT(nzi == zi);
+  }
+  {
+    unsigned int xi = 8, yi = 2, zi = 1;
+    unsigned int idx = grd.getGridIndex(xi, yi, zi);
+    unsigned int nxi, nyi, nzi;
+    grd.getGridIndices(idx, nxi, nyi, nzi);
+    TEST_ASSERT(nxi == xi);
+    TEST_ASSERT(nyi == yi);
+    TEST_ASSERT(nzi == zi);
+  }
+
+  RDGeom::Point3D pt = grd.getOffset();
+  grd.setVal(pt, 2.3);
+  unsigned int id = grd.getGridPointIndex(pt);
+
+  TEST_ASSERT(RDKit::feq(grd.getVal(pt), 2.3));
+  TEST_ASSERT(id == 0);
+  TEST_ASSERT(RDKit::feq(grd.getVal(id), 2.3));
+
+  RDGeom::Point3D pt2 = grd.getGridPointLoc(id);
+  TEST_ASSERT(RDKit::feq(pt.x, pt2.x));
+  TEST_ASSERT(RDKit::feq(pt.y, pt2.y));
+  TEST_ASSERT(RDKit::feq(pt.z, pt2.z));
+}
+
+void test4UniformRealValueGrid3DOps() {
+  RDGeom::UniformRealValueGrid3D grd1(5.0, 5.0, 5.0, 0.1);
+  RDGeom::UniformRealValueGrid3D grd2(5.0, 5.0, 5.0, 0.1);
+
+  grd1.setVal(50, 37.37);
+  grd2.setVal(50, 1.03);
+  RDGeom::UniformRealValueGrid3D grd3(grd2);
+
+  grd1 |= grd2;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 37.37));
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), 1.03));
+
+  grd2 = grd2 | grd1;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 37.37));
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), 37.37));
+
+  grd2 = grd2 & grd3;
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), 1.03));
+  TEST_ASSERT(RDKit::feq(grd3.getVal(50), 1.03));
+
+  grd3 &= grd1;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 37.37));
+  TEST_ASSERT(RDKit::feq(grd3.getVal(50), 1.03));
+
+  grd1 += grd2;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 38.40));
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), 1.03));
+
+  grd1 -= grd2;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 37.37));
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), 1.03));
+
+  grd2 = grd2 - grd1;
+  TEST_ASSERT(RDKit::feq(grd1.getVal(50), 37.37));
+  TEST_ASSERT(RDKit::feq(grd2.getVal(50), -36.34));
+}
+
+int main() {
+  std::cout << "***********************************************************\n";
+  std::cout << "Testing RealValueGrid3D\n";
+
+  std::cout << "\t---------------------------------\n";
+  std::cout << "\t test1UniformRealValueGrid3D \n\n";
+  test1UniformRealValueGrid3D();
+
+  std::cout << "\t---------------------------------\n";
+  std::cout << "\t test2UniformRealValueGrid3DPickling \n\n";
+  test2UniformRealValueGrid3DPickling();
+
+  std::cout << "\t---------------------------------\n";
+  std::cout << "\t test3UniformRealValueGrid3DIndexing \n\n";
+  test3UniformRealValueGrid3DIndexing();
+
+  std::cout << "\t---------------------------------\n";
+  std::cout << "\t test4UniformRealValueGrid3DOps \n\n";
+  test4UniformRealValueGrid3DOps();
+
+  return 0;
+}

Code/GraphMol/CMakeLists.txt

@@ -125,6 +125,7 @@ add_subdirectory(MolEnumerator)
 add_subdirectory(Abbreviations)
 add_subdirectory(GeneralizedSubstruct)
 add_subdirectory(MolProcessing)
+add_subdirectory(MolInteractionFields)
 
 
 rdkit_test(graphmolTest1 test1.cpp LINK_LIBRARIES FileParsers SmilesParse GraphMol

Code/GraphMol/MolInteractionFields/CMakeLists.txt

@@ -0,0 +1,16 @@
+rdkit_library(MolInteractionFields MIFDescriptors.cpp LINK_LIBRARIES ForceFieldHelpers ForceField DataStructs RDGeneral GraphMol SmilesParse)
+target_compile_definitions(MolInteractionFields PRIVATE RDKIT_MOLINTERACTIONFIELDS_BUILD)
+
+rdkit_headers(MIFDescriptors.h DEST GraphMol/MolInteractionFields)
+
+rdkit_catch_test(testMIF testMIF.cpp LINK_LIBRARIES MolInteractionFields DistGeomHelpers DistGeometry ForceFieldHelpers ForceField
+ PartialCharges SubstructMatch FileParsers SmilesParse
+MolAlign MolTransforms
+GraphMol
+Optimizer EigenSolvers
+Alignment
+DataStructs RDGeometryLib RDGeneral
+${RDKit_THREAD_LIBS} )
+
+add_subdirectory(Wrap)
+

Code/GraphMol/MolInteractionFields/MIFDescriptors.h

@@ -0,0 +1,562 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <RDGeneral/export.h>
+#ifndef RDMIF_DESCRIPTORS_H
+#define RDMIF_DESCRIPTORS_H
+/*! \file MIFDescriptors.h
+
+  \brief code for generating molecular interaction field (MIF) descriptors
+
+  \b Note that this functionality is experimental and the API and/or results
+  amay change in future releases.
+*/
+
+#include <vector>
+#include <Geometry/point.h>
+#include <Geometry/UniformRealValueGrid3D.h>
+#include <ForceField/MMFF/Nonbonded.h>
+#include <ForceField/UFF/Nonbonded.h>
+#include <ForceField/UFF/Params.h>
+#include <GraphMol/ForceFieldHelpers/UFF/AtomTyper.h>
+#include <GraphMol/RDKitBase.h>
+
+namespace RDMIF {
+//! \brief constructs a UniformRealValueGrid3D which fits to the molecule mol
+/*!
+ \returns a pointer to a UniformRealValueGrid which has a spacing of \c spacing
+ (default: 0.5 Angstrom) and a spatial extent of conformer \c confId
+ (default: -1) of molecule \c mol plus a defined margin \c margin on each side
+ (default: 5.0 Angstrom).
+
+ \param mol		Molecule for which the grid is constructed
+ \param confId		Id of Conformer which is used
+ \param margin 	distance from molecule to the grid's wall
+ \param spacing	grid's spacing
+
+ <b>Returns</b>
+ pointer to a grid
+ */
+RDKIT_MOLINTERACTIONFIELDS_EXPORT
+std::unique_ptr<RDGeom::UniformRealValueGrid3D> constructGrid(
+    const RDKit::ROMol &mol, int confId = -1, double margin = 5.0,
+    double spacing = 0.5);
+
+//! \brief calculates a descriptor at every grid point of MIF
+/*!
+ any unary function/functor can be used taking a grid point of type
+ RDKit::Point3D as parameter
+
+ \param grd			a UniformRealValueGrid3D
+ \param functor	any unary function/functor which takes four doubles as
+ parameters (3 coordinates, 1 threshold)
+ \param thres	  cuts off atom contributions if distance of interaction higher
+ than threshold, negative threshold: no threshold, defaults to -1.0
+ */
+template <typename T>
+void calculateDescriptors(RDGeom::UniformRealValueGrid3D &grd, const T &functor,
+                          double thres = -1.0) {
+  const RDGeom::Point3D &offSet = grd.getOffset();
+  auto x = offSet.x;
+  auto y = offSet.y;
+  auto z = offSet.z;
+  auto oX = x;
+  auto oY = y;
+  auto spacing = grd.getSpacing();
+  if (thres < 0) {
+    thres = spacing * grd.getSize();
+  }
+  thres *= thres;  // comparing squared distance
+
+  unsigned int id = 0;
+  auto &data = grd.getData();
+
+  for (unsigned int idZ = 0; idZ < grd.getNumZ(); ++idZ) {
+    for (unsigned int idY = 0; idY < grd.getNumY(); ++idY) {
+      for (unsigned int idX = 0; idX < grd.getNumX(); ++idX) {
+        data[id++] = functor(x, y, z, thres);
+        x += spacing;
+      }
+      y += spacing;
+      x = oX;
+    }
+    z += spacing;
+    y = oY;
+  }
+}
+
+//! \brief class for calculation of electrostatic interaction (Coulomb energy)
+//! between probe and molecule in vacuum (no dielectric)
+/*
+ d_nAtoms     No of atoms in molecule
+ d_softcore   true if softcore interaction is used, else minimum cutoff distance
+ is used
+ d_probe      charge of probe [e]
+ d_absVal     if true, absolute values of interactions are calculated
+ d_alpha      softcore interaction parameter [A^2]
+ d_cutoff     squared minimum cutoff distance [A^2]
+ d_charges    vector of doubles with all partial charges of the atoms in a
+ molecule [e] d_pos        vector of doubles with all positions of the atoms in
+ a molecule [A] d_prefactor  prefactor taking into account the geometric
+ factors,natural constants and conversion of units
+ */
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT Coulomb {
+ public:
+  Coulomb() = default;
+  ~Coulomb() = default;
+
+  //! \brief constructs Coulomb object from vectors of charges and positions
+  /*!
+   \param charges     vector of charges [e]
+   \param pos         vector of postions [A]
+   \param probeCharge charge of probe [e] (default: 1.0)
+   \param absVal      if true, negative (favored) values of interactions are
+   calculated (default: false)
+   \param alpha       softcore interaction parameter [A^2]; if zero,
+   a minimum cutoff distance is used (default: 0.0)
+   \param cutoff      minimum cutoff distance [A] (default: 1.0)
+
+   */
+  Coulomb(const std::vector<double> &charges,
+          const std::vector<RDGeom::Point3D> &positions,
+          double probeCharge = 1.0, bool absVal = false, double alpha = 0.0,
+          double cutoff = 1.0);
+
+  //! \brief constructs Coulomb object from a molecule object
+  /*!
+   The molecule \c mol needs to have partial charges set as property of atoms.
+
+   \param mol         molecule object
+   \param confId      conformation id which is used to get positions of atoms
+   (default=-1)
+   \param absVal      if true, absolute values of interactions are
+   calculated (default: false)
+   \param probeCharge charge of probe [e] (default: 1.0)
+   \param prop	     property key for retrieving partial charges
+   of atoms (default: "_GasteigerCharge")
+   \param alpha       softcore interaction parameter [A^2], if zero, a minimum
+   cutoff distance is used (default: 0.0)
+   \param cutoff      minimum cutoff distance [A] (default: 1.0)
+
+   */
+  Coulomb(const RDKit::ROMol &mol, int confId = -1, double probeCharge = 1.0,
+          bool absVal = false, const std::string &prop = "_GasteigerCharge",
+          double alpha = 0.0, double cutoff = 1.0);
+
+  //! \brief calculated the electrostatic interaction at point \c pt in the
+  //! molecules field in [kJ mol^-1]
+  /*!
+   \param x, y, z     coordinates at which the interaction is calculated
+   \param thres       squared max distance until interactions are calc.
+
+   \return electrostatic interaction energy in [kJ mol^-1]
+   */
+  double operator()(double x, double y, double z, double thres) const;
+
+ private:
+  unsigned int d_nAtoms = 0;
+  bool d_softcore = false;
+  bool d_absVal = false;
+  double d_cutoff = 0.00001;
+  double d_probe = 1.0;
+  double d_alpha = 0.0;
+  std::vector<double> d_charges;
+  std::vector<double> d_pos;
+};
+
+//! \brief class for calculation of electrostatic interaction (Coulomb energy)
+//! between probe and molecule by taking a distance-dependent dielectric into
+//! account:
+/*!
+ Same energy term as used in GRID
+ References:
+ J. Med. Chem. 1985, 28, 849.
+ J. Comp. Chem. 1983, 4, 187.
+ */
+/*
+ member variables:
+ d_nAtoms:     No of atoms in molecule
+ d_softcore:   true if softcore interaction is used, else minimum cutoff
+ distance is used
+ d_dielectric: factor of dielectric constants
+ d_probe:      charge of probe [e]
+ d_absVal:     if true, negative (favored) values of interactions
+ are calculated (default: false)
+ d_epsilon:    relative permittivity of solvent
+ d_xi:         relative permittivity of solute
+ d_alpha:      softcore interaction parameter [A^2]
+ d_cutoff:     squared minimum cutoff distance [A^2]
+ d_charges:    vector of doubles with all partial charges of the atoms
+ in a molecule [e]
+ d_pos:        vector of Point3Ds with all positions of the atoms
+ in a molecule [A]
+ d_prefactor:  prefactor taking into account the geometric factors,
+ natural constants and conversion of units
+ */
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT CoulombDielectric {
+ public:
+  CoulombDielectric()
+      : d_nAtoms(0),
+        d_softcore(false),
+        d_absVal(false),
+        d_cutoff(0.001),
+        d_probe(1),
+        d_epsilon(1.0),
+        d_xi(1.0),
+        d_alpha(0.0) {
+    d_dielectric = (d_xi - d_epsilon) / (d_xi + d_epsilon);
+  }
+
+  //! \brief constructs CoulombDielectric object from vectors of charges and
+  //! positions
+  /*!
+   \param charges     vector of charges [e]
+   \param pos         vector of postions [A]
+   \param probeCharge charge of probe [e] (default: 1.0)
+   \param absVal      if true, negative (favored) values of interactions are
+   calculated (default: false)
+   \param alpha       softcore interaction parameter [A^2]; if zero,
+   a minimum cutoff distance is used (default: 0.0)
+   \param cutoff      minimum cutoff distance [A] (default: 1.0)
+   \param epsilon     relative permittivity of solvent (default: 80.0)
+   \param xi          relative permittivity of solute (default: 4.0)
+
+   */
+  CoulombDielectric(const std::vector<double> &charges,
+                    const std::vector<RDGeom::Point3D> &positions,
+                    double probeCharge = 1.0, bool absVal = false,
+                    double alpha = 0.0, double cutoff = 1.0,
+                    double epsilon = 80.0, double xi = 4.0);
+
+  //! \brief constructs Coulomb object from a molecule object
+  /*!
+   The molecule \c mol needs to have partial charges set as property of atoms.
+
+   \param mol         molecule object
+   \param confId      conformation id which is used to get positions of atoms
+   (default=-1)
+   \param probeCharge charge of probe [e] (default: 1.0)
+   \param absVal      if true, negative (favored) values of interactions are
+   calculated (default: false)
+   \param prop	       property key for retrieving partial charges of atoms
+   (default: "_GasteigerCharge")
+   \param alpha       softcore interaction parameter [A^2], if zero, a minimum
+   cutoff distance is used (default: 0.0)
+   \param cutoff      minimum cutoff distance [A] (default: 1.0)
+   \param epsilon     relative permittivity of solvent (default: 80.0)
+   \param xi          relative permittivity of solute (default: 4.0)
+   */
+  CoulombDielectric(const RDKit::ROMol &mol, int confId = -1,
+                    double probeCharge = 1.0, bool absVal = false,
+                    const std::string &prop = "_GasteigerCharge",
+                    double alpha = 0.0, double cutoff = 1.0,
+                    double epsilon = 80.0, double xi = 4.0);
+
+  //! \brief returns the electrostatic interaction at point \c pt in the
+  //! molecule's field in [kJ mol^-1]
+  /*!
+   \param x, y, z     coordinates at which the interaction is calculated
+   \param thres       squared threshold distance
+
+   \return electrostatic interaction energy in [kJ mol^-1]
+   */
+  double operator()(double x, double y, double z, double thres) const;
+
+ private:
+  unsigned int d_nAtoms;
+  bool d_softcore, d_absVal;
+  double d_cutoff, d_probe;
+  double d_epsilon, d_xi, d_alpha;
+  double d_dielectric;
+  std::vector<double> d_charges, d_sp;
+  std::vector<double> d_pos;
+  // used as a cache
+  mutable std::vector<double> d_dists;
+};
+
+//! \brief Abstract class for calculation of Van der Waals interaction between
+//! probe and molecule at gridpoint \c pt
+/*
+ * Either the MMFF94 or the UFF VdW term can be calculated with a defined probe
+ * atom and molecule.
+ * d_cutoff: minimum cutoff distance [A]
+ * d_nAtoms: No of atoms in molecule
+ * d_R_star_ij: vector of Lennard Jones parameters for every
+ * interaction (vdW-radii)
+ * d_wellDepth: vector of Lennard Jones parameters for
+ * every interaction (well depths)
+ * d_pos: vector of positions of atoms in  molecule
+ */
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT VdWaals {
+ public:
+  VdWaals() = default;
+  VdWaals(const RDKit::ROMol &mol, int confId = -1, double cutoff = 1.0);
+  VdWaals(const VdWaals &other) = delete;
+  VdWaals &operator=(const VdWaals &other) = delete;
+  VdWaals(VdWaals &&other) = default;
+  VdWaals &operator=(VdWaals &&other) = default;
+  virtual ~VdWaals() = default;
+
+  //! \brief returns the VdW interaction at point \c pt in the molecules field
+  //! in [kJ mol^-1]
+  /*!
+   \param x, y, z     coordinates at which the interaction is calculated
+   \param thres       squared max distance until interactions are calc.
+
+   \return vdW interaction energy in [kJ mol^-1]
+   */
+  double operator()(double x, double y, double z, double thres) const;
+
+ protected:
+  void fillVectors();
+  virtual double calcEnergy(double, double, double) const = 0;
+  virtual void fillVdwParamVectors(unsigned int atomIdx) = 0;
+  double d_cutoff = 1.0;
+  unsigned int d_nAtoms = 0;
+  std::vector<double> d_pos;
+  std::vector<double> d_R_star_ij;
+  std::vector<double> d_wellDepth;
+  std::unique_ptr<RDKit::ROMol> d_mol;
+};
+
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT MMFFVdWaals : public VdWaals {
+ public:
+  //! \brief constructs VdWaals object which uses MMFF94 from a molecule object
+  /*!
+  \param mol           molecule object
+  \param confId        conformation id which is used to get positions of atoms
+  (default: -1)
+  \param probeAtomType MMFF94 atom type for the probe atom
+  (default: 6, sp3 oxygen)
+  \param cutoff        minimum cutoff distance [A] (default: 1.0)
+  \param scaling       scaling of VdW parameters to take hydrogen bonds into
+  account (default: false)
+  */
+  MMFFVdWaals(const RDKit::ROMol &mol, int confId = -1,
+              unsigned int probeAtomType = 6, bool scaling = false,
+              double cutoff = 1.0);
+  MMFFVdWaals(const MMFFVdWaals &other) = delete;
+  MMFFVdWaals &operator=(const MMFFVdWaals &other) = delete;
+  MMFFVdWaals(MMFFVdWaals &&other) = default;
+  MMFFVdWaals &operator=(MMFFVdWaals &&other) = default;
+  ~MMFFVdWaals() = default;
+
+ private:
+  double calcEnergy(double, double, double) const;  // MMFF energy function
+  void fillVdwParamVectors(unsigned int atomIdx);
+  bool d_scaling;
+  std::unique_ptr<RDKit::MMFF::MMFFMolProperties> d_props;
+  const ForceFields::MMFF::MMFFVdWCollection *d_mmffVdW;
+  const ForceFields::MMFF::MMFFVdW *d_probeParams;
+};
+
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT UFFVdWaals : public VdWaals {
+ public:
+  //! \brief constructs VdWaals object which uses UFF from a molecule object
+  /*!
+  \param mol           molecule object
+  \param confId        conformation id which is used to get positions of atoms
+  (default: -1)
+  \param probeAtomType UFF atom type for the probe atom (e.g. "O_3", i.e. a
+  tetrahedral oxygen)
+  \param cutoff        minimum cutoff distance [A] (default: 1.0)
+  */
+  UFFVdWaals(const RDKit::ROMol &mol, int confId = -1,
+             const std::string &probeAtomType = "O_3", double cutoff = 1.0);
+  UFFVdWaals(const UFFVdWaals &other) = delete;
+  UFFVdWaals &operator=(const UFFVdWaals &other) = delete;
+  UFFVdWaals(UFFVdWaals &&other) = default;
+  UFFVdWaals &operator=(UFFVdWaals &&other) = default;
+  ~UFFVdWaals() = default;
+
+ private:
+  double calcEnergy(double, double, double) const;  // UFF energy function
+  void fillVdwParamVectors(unsigned int atomIdx);
+  const ForceFields::UFF::ParamCollection *d_uffParamColl;
+  const ForceFields::UFF::AtomicParams *d_probeParams;
+  RDKit::UFF::AtomicParamVect d_params;
+};
+
+//! \brief class for calculation of hydrogen bond potential between probe and
+//! molecule
+/*
+ implementation of GRID molecular H-Bond descriptor:
+ J.Med.Chem. 1989, 32, 1083.
+ J.Med.Chem. 1993, 36, 140.
+ J.Med.Chem. 1993, 36, 148.
+ member variables:
+ d_cutoff:       squared minimum cutoff distance [A^2]
+ d_nInteract:    No of interactions between probe and molecule
+ d_DAprop:		  either 'A' or 'D', defines whether target atoms (in
+ molecule) have to be acceptors or donors for interacting with probe
+ d_probetype:    defines type of probe atom, either N or O
+ d_targettypes:  vectors of types of target atoms (in molecule), either N or O
+ d_pos:          vector of positions of target atoms in molecule
+ d_direction:    if target accepting, vector of directions of lone pairs on
+ target atoms (if there are two lone pairs, the resulting direction is used) if
+ target donating, the X-H bond direction is saved
+ d_plane:        vector of plane vectors in which the lone pairs are
+ */
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT HBond {
+ public:
+  HBond() : d_cutoff(1.0), d_probetype(O), d_nInteract(0) {};
+
+  //! \brief constructs HBond object from a molecule object
+  /*!
+   The molecule \c mol needs to have partial charges set as property of atoms.
+   \param mol           molecule object
+   \param confId        conformation id which is used to get positions of atoms
+   (default=-1)
+   \param probeAtomType atom type for the probe atom ("OH", "O", "NH", "N")
+   \param fixed         for some groups, two different angle
+   dependencies are defined in GRID: one which takes some flexibility of groups
+   (rotation/swapping of lone pairs and hydrogen) into account and one for
+   strictly fixed conformations if true, strictly fixed conformations
+   (default: fixed)
+   \param cutoff        minimum cutoff distance [A] (default: 1.0)
+   */
+  HBond(const RDKit::ROMol &mol, int confId = -1,
+        const std::string &probeAtomType = "OH", bool fixed = true,
+        double cutoff = 1.0);
+  ~HBond() = default;
+
+  //! \brief returns the hydrogen bonding interaction at point \c pt in the
+  //! molecules field in [kJ mol^-1]
+  /*!
+   \param x, y, z     coordinates at which the interaction is calculated
+   \param thres       squared max distance until interactions are calc.
+
+   \return hydrogen bonding interaction energy in [kJ mol^-1]
+   */
+  double operator()(double x, double y, double z, double thres) const;
+
+  unsigned int getNumInteractions() const { return d_nInteract; }
+
+ private:
+  boost::uint8_t d_DAprop;
+
+  enum atomtype {
+    N,
+    O
+  };
+  double d_cutoff;
+  atomtype d_probetype;
+  unsigned int d_nInteract;  // number of HBond interactions
+
+  std::vector<atomtype> d_targettypes;
+  std::vector<double> d_pos;        // vector of positions of target atoms
+  std::vector<double> d_direction;  // vector of sum of lone pair vectors
+  std::vector<double> d_lengths;    // vector of lengths of direction vectors
+  std::vector<double> d_plane;      // vector of lone pair plane vectors
+  // these two are used as a computational cache during operator()
+  mutable std::vector<double>
+      d_eneContrib;  // energy contributions of all interactions
+  mutable std::vector<double>
+      d_vectTargetProbe;  // hydrogen bond direction of probe
+
+  // constructing functions
+  unsigned int findSpecials(const RDKit::ROMol &mol, int confId, bool fixed,
+                            std::vector<unsigned int> &specials);
+  unsigned int findAcceptors(const RDKit::ROMol &mol, int confId,
+                             const std::vector<unsigned int> &specials);
+  unsigned int findDonors(const RDKit::ROMol &mol, int confId,
+                          const std::vector<unsigned int> &specials);
+  unsigned int findAcceptorsUnfixed(const RDKit::ROMol &mol, int confId,
+                                    const std::vector<unsigned int> &specials);
+  unsigned int findDonorsUnfixed(const RDKit::ROMol &mol, int confId,
+                                 const std::vector<unsigned int> &specials);
+
+  void addVectElements(atomtype type, double (*funct)(double, double, double),
+                       const RDGeom::Point3D &pos, const RDGeom::Point3D &dir,
+                       const RDGeom::Point3D &plane = RDGeom::Point3D(0.0, 0.0,
+                                                                      0.0));
+
+  void normalize(double &x, double &y, double &z) const;
+  double angle(double x1, double y1, double z1, double x2, double y2,
+               double z2) const;
+
+  std::vector<double (*)(double, double, double)> d_function;
+};
+
+//! \brief class for calculation of hydrophilic field of molecule
+/*
+ interaction energy of hydrogen and oxygen of water with the molecule is
+ calculated at each point as a hydrogen bond interaction (either OH or O probe).
+ the favored interaction is returned. member variables: d_hbondO
+ HBond descriptor with "O" probe d_hbondOH				   HBond
+ descriptor with "OH" probe
+ */
+class RDKIT_MOLINTERACTIONFIELDS_EXPORT Hydrophilic {
+ public:
+  //! \brief constructs Hydrophilic object from a molecule object
+  /*!
+   params:
+   \param mol           		 molecule object
+   \param confId        		 conformation id which is used to get
+   positions of atoms (default: -1)
+   \param fixed         		 for some groups, two different angle
+   dependencies are defined in GRID: one which takes some flexibility of groups
+   (rotation/swapping of lone pairs and hydrogen) into account and one for
+   strictly fixed conformations if true, strictly fixed conformations
+   (default: fixed)
+   \param cutoff  minimum cutoff distance [A] (default: 1.0)
+   */
+  Hydrophilic(const RDKit::ROMol &mol, int confId = -1, bool fixed = true,
+              double cutoff = 1.0);
+  ~Hydrophilic() {};
+
+  double operator()(double x, double y, double z, double thres) const;
+
+ private:
+  HBond d_hbondOH, d_hbondO;
+};
+
+//! \brief writes the contents of the MIF to a stream
+/*!
+ The Grid \c grd is written in Gaussian Cube format
+ A molecule \c mol and a \c confId can optionally be provided
+ */
+RDKIT_MOLINTERACTIONFIELDS_EXPORT void writeToCubeStream(
+    const RDGeom::UniformRealValueGrid3D &grd, std::ostream &outStrm,
+    const RDKit::ROMol *mol = nullptr, int confid = -1);
+
+//! \brief writes the contents of the MIF to a file
+/*!
+ The Grid \c grd is written in Gaussian Cube format
+ A molecule \c mol and a \c confId can optionally be provided
+ */
+RDKIT_MOLINTERACTIONFIELDS_EXPORT void writeToCubeFile(
+    const RDGeom::UniformRealValueGrid3D &grd, const std::string &filename,
+    const RDKit::ROMol *mol = nullptr, int confid = -1);
+
+//! \brief reads the contents of the MIF from a stream in Gaussian cube format
+/*!
+ The Grid \c grd is modified according to input values
+ If a molecule was associated to the grid on write, a non-null pointer to a
+ molecule is returned with atoms and a conformer, but NO bond information.
+ If there is no atom information in the cube file, a null pointer is returned.
+ */
+RDKIT_MOLINTERACTIONFIELDS_EXPORT std::unique_ptr<RDKit::RWMol>
+readFromCubeStream(RDGeom::UniformRealValueGrid3D &grd, std::istream &inStrm);
+
+//! \brief reads the contents of the MIF from a file in Gaussian cube format
+/*!
+ The Grid \c grd is modified according to input values
+ If a molecule was associated to the grid on write, a non-null pointer to a
+ molecule is returned with atoms and a conformer, but NO bond information.
+ If there is no atom information in the cube file, a null pointer is returned.
+ */
+RDKIT_MOLINTERACTIONFIELDS_EXPORT std::unique_ptr<RDKit::RWMol>
+readFromCubeFile(RDGeom::UniformRealValueGrid3D &grd,
+                 const std::string &filename);
+
+}  // end of namespace RDMIF
+
+#endif /* RDMIF_DESCRIPTORS_H */

Code/GraphMol/MolInteractionFields/Wrap/CMakeLists.txt

@@ -0,0 +1,8 @@
+rdkit_python_extension(rdMIF 
+                       rdMIF.cpp
+                       DEST Chem
+                       LINK_LIBRARIES MolInteractionFields ForceFieldHelpers
+                       ForceField PartialCharges SubstructMatch SmilesParse GraphMol RDGeometryLib
+                       RDGeneral DataStructs RDBoost)
+
+add_pytest(pyMIF  ${CMAKE_CURRENT_SOURCE_DIR}/testMIF.py)

Code/GraphMol/MolInteractionFields/Wrap/rdMIF.cpp

@@ -0,0 +1,422 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <RDBoost/Wrap.h>
+#include <RDBoost/PySequenceHolder.h>
+#include <boost/python.hpp>
+#include <ForceField/MMFF/Nonbonded.h>
+#include <GraphMol/MolInteractionFields/MIFDescriptors.h>
+#include <Geometry/UniformRealValueGrid3D.h>
+#include <Geometry/point.h>
+
+#include <RDBoost/boost_numpy.h>
+
+namespace python = boost::python;
+using namespace RDMIF;
+
+void wrap_mif();
+
+BOOST_PYTHON_MODULE(rdMIF) {
+  python::scope().attr("__doc__") =
+      "Module containing functions for calculating molecular interaction fields (MIFs)\n\
+  NOTE: This functionality is experimental and the API and/or results may change in future releases.";
+  python::register_exception_translator<IndexErrorException>(
+      &translate_index_error);
+  python::register_exception_translator<ValueErrorException>(
+      &translate_value_error);
+
+  wrap_mif();
+}
+
+namespace RDMIF {
+
+RDGeom::UniformRealValueGrid3D *constructGridHelper(const RDKit::ROMol &mol,
+                                                    int confId, double margin,
+                                                    double spacing) {
+  return constructGrid(mol, confId, margin, spacing).release();
+}
+
+std::pair<std::vector<double>, std::vector<RDGeom::Point3D>>
+extractChargesAndPositions(const python::object &charges,
+                           const python::object &positions) {
+  const auto pyPos = positions.ptr();
+  const auto pyCharges = charges.ptr();
+  if (!pyPos || !PySequence_Check(pyPos)) {
+    throw_value_error("positions argument must be a sequence");
+  }
+  if (!pyCharges || !PySequence_Check(pyCharges)) {
+    throw_value_error("charges argument must be a sequence");
+  }
+  auto nrows = PySequence_Size(pyPos);
+  if (nrows != PySequence_Size(pyCharges)) {
+    throw_value_error("positions and charges must have the same length");
+  }
+  std::vector<RDGeom::Point3D> pos(nrows);
+  std::vector<double> ch(nrows);
+  for (unsigned int i = 0; i < nrows; ++i) {
+    const auto pyXyz = PySequence_GetItem(pyPos, i);
+    if (!pyXyz || !PySequence_Check(pyXyz) || PySequence_Size(pyXyz) != 3) {
+      throw_value_error(
+          "all elements in positions argument must be x,y,z sequences");
+    }
+    pos[i].x = python::extract<double>(PySequence_GetItem(pyXyz, 0));
+    pos[i].y = python::extract<double>(PySequence_GetItem(pyXyz, 1));
+    pos[i].z = python::extract<double>(PySequence_GetItem(pyXyz, 2));
+    ch[i] = python::extract<double>(PySequence_GetItem(pyCharges, i));
+  }
+  return std::make_pair(std::move(ch), std::move(pos));
+}
+
+Coulomb *makeAltCoulomb(const python::object &charges,
+                        const python::object &positions, double probecharge,
+                        bool absVal, double alpha, double cutoff) {
+  const auto [ch, pos] = extractChargesAndPositions(charges, positions);
+  return new Coulomb(ch, pos, probecharge, absVal, alpha, cutoff);
+}
+
+CoulombDielectric *makeAltCoulombDielectric(const python::object &charges,
+                                            const python::object &positions,
+                                            double probecharge, bool absVal,
+                                            double alpha, double cutoff,
+                                            double epsilon, double xi) {
+  const auto [ch, pos] = extractChargesAndPositions(charges, positions);
+  return new CoulombDielectric(ch, pos, probecharge, absVal, alpha, cutoff,
+                               epsilon, xi);
+}
+
+python::tuple readCubeFile(const std::string &filename) {
+  std::unique_ptr<RDGeom::UniformRealValueGrid3D> grd(
+      new RDGeom::UniformRealValueGrid3D());
+  auto res = readFromCubeFile(*grd, filename);
+  boost::python::manage_new_object::apply<
+      RDGeom::UniformRealValueGrid3D *>::type grdConverter;
+  boost::python::manage_new_object::apply<RDKit::ROMol *>::type molConverter;
+  return python::make_tuple(python::handle<>(grdConverter(grd.release())),
+                            python::handle<>(molConverter(
+                                static_cast<RDKit::ROMol *>(res.release()))));
+}
+
+struct mif_wrapper {
+  static void wrap() {
+    std::string docStringClass =
+        "Class for calculation of electrostatic interaction (Coulomb energy) between probe and molecule in\n\
+        vacuum (no dielectric).\n\n";
+    std::string docStringConst =
+        "Constructor for Coulomb class.\n\n\
+        ARGUMENTS:\n\
+        - mol:           the molecule of interest\n\
+        - confId:        the ID of the conformer to be used (defaults to -1)\n\
+        - probeCharge    charge of probe [e] (defaults to 1.0 e)\n\
+        - absVal:        if True, absolute values of interactions are calculated (defaults to False)\n\
+        - chargeKey      property key for retrieving partial charges of atoms from molecule (defaults to '_GasteigerCharge')\n\
+        - softcoreParam  softcore interaction parameter [A^2], if zero, a minimum cutoff distance is used (defaults to 0.0)\n\
+        - cutoff         minimum cutoff distance [A] (defaults to 1.0)\n";
+    std::string docStringConstAlt =
+        "Alternative constructor for Coulomb class.\n\n\
+        ARGUMENTS:\n\
+        - charges:       array of partial charges of a molecule's atoms\n\
+        - positions:     array of positions of a molecule's atoms\n\
+        - probeCharge    charge of probe [e] (defaults to 1.0 e)\n\
+        - absVal:        if True, absolute values of interactions are calculated (defaults to False)\n\
+        - softcoreParam  softcore interaction parameter [A^2], if zero, a minimum cutoff distance is used (defaults to 0.0)\n\
+        - cutoff         minimum cutoff distance [A] (defaults to 1.0)\n";
+    std::string docString =
+        "Calculates the electrostatic interaction (Coulomb energy) between probe and molecule in\n\
+        vacuum (no dielectric).\n\n\
+        ARGUMENTS:\n\
+        - x, y, z:   coordinates of probe position for energy calculation\n\
+        - threshold: maximal distance until which interactions are calculated\n\
+        RETURNS:\n\
+        - electrostatic potential in [kJ mol^-1]\n";
+    python::class_<Coulomb>(
+        "Coulomb", docStringClass.c_str(),
+        python::init<const RDKit::ROMol &, int, double, bool,
+                     const std::string &, double, double>(
+            (python::arg("mol"), python::arg("confId") = -1,
+             python::arg("probeCharge") = 1.0, python::arg("absVal") = false,
+             python::arg("chargeKey") = "_GasteigerCharge",
+             python::arg("softcoreParam") = 0.0, python::arg("cutoff") = 1.0),
+            docStringConst.c_str()))
+        .def("__init__",
+             python::make_constructor(
+                 makeAltCoulomb, python::default_call_policies(),
+                 (python::arg("charges"), python::arg("positions"),
+                  python::arg("probeCharge") = 1.0,
+                  python::arg("absVal") = false,
+                  python::arg("softcoreParam") = 0.0,
+                  python::arg("cutoff") = 1.0)),
+             docStringConstAlt.c_str())
+        .def("__call__", &Coulomb::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docStringClass =
+        "Class for calculation of electrostatic interaction (Coulomb energy) between probe and molecule in\n\
+        by taking a distance-dependent dielectric into account.\n\
+        Same energy term as used in GRID MIFs.\n\
+        References:\n\
+        - J. Med. Chem. 1985, 28, 849.\n\
+        - J. Comp. Chem. 1983, 4, 187.\n\n";
+    docStringConst =
+        "Constructor for CoulombDielectric class.\n\n\
+        ARGUMENTS:\n\
+        - mol:           the molecule of interest\n\
+        - confId:        the ID of the conformer to be used (defaults to -1)\n\
+        - probeCharge    charge of probe [e] (defaults to 1.0 e)\n\
+        - absVal:        if True, absolute values of interactions are calculated (defaults to False)\n\
+        - chargeKey       property key for retrieving partial charges of atoms from molecule (defaults to '_GasteigerCharge')\n\
+        - softcoreParam  softcore interaction parameter [A^2], if zero, a minimum cutoff distance is used (defaults to 0.0)\n\
+        - cutoff         minimum cutoff distance [A] (defaults to 1.0)\n\
+        - epsilon        relative permittivity of solvent (defaults to 80.0)\n\
+        - xi             relative permittivity of solute (defaults to 4.0)\n";
+    docStringConstAlt =
+        "Alternative constructor for CoulombDielectric class.\n\n\
+        ARGUMENTS:\n\
+      - charges:       array of partial charges of a molecule's atoms\n\
+      - positions:     array of positions of a molecule's atoms\n\
+      - probeCharge    charge of probe [e] (defaults to 1.0 e)\n\
+      - absVal:        if True, absolute values of interactions are calculated (defaults to False)\n\
+      - softcoreParam  softcore interaction parameter [A^2], if zero, a minimum cutoff distance is used (defaults to 0.0)\n\
+      - cutoff         minimum cutoff distance [A] (defaults to 1.0)\n\
+      - epsilon        relative permittivity of solvent (defaults to 80.0)\n\
+      - xi             relative permittivity of solute (defaults to 4.0)\n";
+    docString =
+        "Calculates the electrostatic interaction (Coulomb energy) between probe and molecule in\n\
+        by taking a distance-dependent dielectric into account.\n\n\
+        ARGUMENTS:\n\
+        - x, y, z:   coordinates of probe position for energy calculation\n\
+        - threshold: maximal distance until which interactions are calculated\n\
+        RETURNS:\n\
+        - electrostatic potential in [kJ mol^-1]\n";
+    python::class_<CoulombDielectric>(
+        "CoulombDielectric", docStringClass.c_str(),
+        python::init<const RDKit::ROMol &, int, double, bool,
+                     const std::string &, double, double, double, double>(
+            (python::arg("mol"), python::arg("confId") = -1,
+             python::arg("probeCharge") = 1.0, python::arg("absVal") = false,
+             python::arg("chargeKey") = "_GasteigerCharge",
+             python::arg("softcoreParam") = 0.0, python::arg("cutoff") = 1.0,
+             python::arg("epsilon") = 80.0, python::arg("xi") = 4.0),
+            docStringConst.c_str()))
+        .def("__init__",
+             python::make_constructor(
+                 makeAltCoulombDielectric, python::default_call_policies(),
+                 (python::arg("charges"), python::arg("positions"),
+                  python::arg("probeCharge") = 1.0,
+                  python::arg("absVal") = false,
+                  python::arg("softcoreParam") = 0.0,
+                  python::arg("cutoff") = 1.0, python::arg("epsilon") = 80.0,
+                  python::arg("xi") = 4.0)),
+             docStringConstAlt.c_str())
+        .def(python::init<const std::string &>(
+            python::args("self", "pklString")))
+        .def("__call__", &CoulombDielectric::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docStringClass =
+        "Class for calculating van der Waals interactions between molecule and a probe at a gridpoint\
+        based on the MMFF forcefield.\n";
+    docStringConst =
+        "ARGUMENTS:\n\
+        - mol           molecule object\n\
+        - confId        conformation id which is used to get positions of atoms (default=-1)\n\
+        - probeAtomType MMFF94 atom type for the probe atom (default=6, sp3 oxygen)\n\
+        - cutoff        minimum cutoff distance [A] (default:1.0)\n\
+        - scaling       scaling of VdW parameters to take hydrogen bonds into account (default=False)\n";
+    docString =
+        "Calculates the van der Waals interaction between molecule and a probe at a gridpoint.\n\n\
+        ARGUMENTS:\n\
+        - x, y, z:   coordinates of probe position for energy calculation\n\
+        - threshold: maximal distance until which interactions are calculated\n\
+        RETURNS:\n\
+        - van der Waals potential in [kJ mol^-1]\n";
+    python::class_<MMFFVdWaals, MMFFVdWaals *, boost::noncopyable>(
+        "MMFFVdWaals", docStringClass.c_str(),
+        python::init<const RDKit::ROMol &, int, unsigned int, bool, double>(
+            (python::arg("self"), python::arg("mol"),
+             python::arg("confId") = -1, python::arg("probeAtomType") = 6,
+             python::arg("scaling") = false, python::arg("cutoff") = 1.0),
+            docStringConst.c_str()))
+        .def("__call__", &MMFFVdWaals::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docStringClass =
+        "Class for calculating van der Waals interactions between molecule and a probe at a gridpoint\
+        based on the UFF forcefield.\n";
+    docStringConst =
+        "ARGUMENTS:\n\
+        - mol           molecule object\n\
+        - confId        conformation id which is used to get positions of atoms (default=-1)\n\
+        - probeAtomType UFF atom type for the probe atom (default='O_3', sp3 oxygen)\n\
+        - cutoff        minimum cutoff distance [A] (default:1.0)\n";
+    python::class_<UFFVdWaals, UFFVdWaals *, boost::noncopyable>(
+        "UFFVdWaals", docStringClass.c_str(),
+        python::init<const RDKit::ROMol &, int, const std::string &, double>(
+            (python::arg("self"), python::arg("mol"),
+             python::arg("confId") = -1, python::arg("probeAtomType") = "O_3",
+             python::arg("cutoff") = 1.0),
+            docStringConst.c_str()))
+        .def("__call__", &UFFVdWaals::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docStringClass =
+        "Class for calculation of hydrogen bonding energy between a probe and a molecule.\n\n\
+        Similar to GRID hydrogen bonding descriptors.\n\
+        References:\n\
+        - J.Med.Chem. 1989, 32, 1083.\n\
+        - J.Med.Chem. 1993, 36, 140.\n\
+        - J.Med.Chem. 1993, 36, 148.\n";
+    docStringConst =
+        "Constructor for HBond class.\n\n\
+        ARGUMENTS:\n\
+        - mol:           the molecule of interest\n\
+        - confId:        the ID of the conformer to be used (defaults to -1)\n\
+        - probeAtomType: atom type for the probe atom (either 'OH', 'O', 'NH' or 'N') (defaults to 'OH')\n\
+        - fixed:         for some groups, two different angle dependencies are defined:\n\
+                         one which takes some flexibility of groups (rotation/swapping of lone pairs and hydrogen)\n\
+                         into account and one for strictly fixed conformations\n\
+                         if True, strictly fixed conformations (defaults to True)\n\
+        - cutoff         minimum cutoff distance [A] (defaults to 1.0)\n";
+    docString =
+        "Calculates the hydrogen bonding energy between probe and molecule in\n\n\
+        ARGUMENTS:\n\
+        - x, y, z:   coordinates of probe position for energy calculation\n\
+        - threshold: maximal distance until which interactions are calculated\n\
+        RETURNS:\n\
+        hydrogen bonding energy in [kJ mol^-1]\n";
+    python::class_<HBond>(
+        "HBond", docStringClass.c_str(),
+        python::init<RDKit::ROMol &, int, const std::string &, bool, double>(
+            (python::arg("mol"), python::arg("confId") = -1,
+             python::arg("probeAtomType") = "OH", python::arg("fixed") = true,
+             python::arg("cutoff") = 1.0),
+            docStringConst.c_str()))
+        .def("__call__", &HBond::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docStringClass =
+        "Class for calculation of a hydrophilic potential of a molecule at a point.\n\n\
+        The interaction energy of hydrogen and oxygen of water is calculated at each point as a \n\
+        hydrogen bond interaction (either OH or O probe). The favored interaction is returned.\n";
+    docStringConst =
+        "Constructor for Hydrophilic class.\n\n\
+        ARGUMENTS:\n\
+        - mol:         the molecule of interest\n\
+        - confId:      the ID of the conformer to be used (defaults to -1)\n\
+        - fixed:       for some groups, two different angle dependencies are defined:\n\
+                       one which takes some flexibility of groups (rotation/swapping of lone pairs and hydrogen)\n\
+                       into account and one for strictly fixed conformations\n\
+                       if True, strictly fixed conformations (defaults to True)\n\
+        - cutoff       minimum cutoff distance [A] (default:1.0)\n";
+    docString =
+        "Calculates the hydrophilic field energy at a point.\n\n\
+        ARGUMENTS:\n\
+        - x, y, z:   coordinates of probe position for energy calculation\n\
+        - threshold: maximal distance until which interactions are calculated\n\
+        RETURNS:\n\
+        hydrophilic field energy in [kJ mol^-1]\n";
+    python::class_<Hydrophilic>(
+        "Hydrophilic", docStringClass.c_str(),
+        python::init<RDKit::ROMol &, int, bool, double>(
+            (python::arg("mol"), python::arg("confId") = -1,
+             python::arg("fixed") = true, python::arg("cutoff") = 1.0),
+            docStringConst.c_str()))
+        .def("__call__", &Hydrophilic::operator(),
+             (python::arg("x"), python::arg("y"), python::arg("z"),
+              python::arg("threshold")),
+             docString.c_str());
+
+    docString =
+        "Constructs a UniformRealValueGrid3D (3D grid with real values at gridpoints) fitting to a molecule.\n\n\
+        ARGUMENTS:\n\
+        - mol:     molecule of interest\n\
+        - confId:  the ID of the conformer to be used (defaults to -1)\n\
+        - margin:  minimum distance of molecule to surface of grid [A] (defaults to 5.0 A)\n\
+        - spacing: grid spacing [A] (defaults to 0.5 A)\n";
+    python::def("ConstructGrid", constructGridHelper,
+                (python::arg("mol"), python::arg("confId") = -1,
+                 python::arg("margin") = 5.0, python::arg("spacing") = 0.5),
+                docString.c_str(),
+                python::return_value_policy<python::manage_new_object>());
+
+    docString =
+        "Calculates descriptors (to be specified as parameter) of a molecule at every gridpoint of a grid.\n\n\
+        ARGUMENTS:\n\
+        - grid:      UniformRealValueGrid3D which get the MIF values\n\
+        - descriptor:  Descriptor class which is used to calculate values\n";
+    python::def("CalculateDescriptors", calculateDescriptors<Coulomb>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    python::def("CalculateDescriptors", calculateDescriptors<CoulombDielectric>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    python::def("CalculateDescriptors", calculateDescriptors<MMFFVdWaals>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    python::def("CalculateDescriptors", calculateDescriptors<UFFVdWaals>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    python::def("CalculateDescriptors", calculateDescriptors<HBond>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    python::def("CalculateDescriptors", calculateDescriptors<Hydrophilic>,
+                (python::arg("grid"), python::arg("descriptor"),
+                 python::arg("threshold") = -1.0),
+                docString.c_str());
+
+    docString =
+        "Writes Grid to a file in Gaussian CUBE format.\n\n\
+        ARGUMENTS:\n\
+        - grid:      UniformRealValueGrid3D to be stored\n\
+        - filename:  filename of file to be written\n\
+        - mol:       associated molecule (defaults to None)\n\
+        - confId:    the ID of the conformer to be used (defaults to -1)\n";
+    python::def(
+        "WriteToCubeFile", writeToCubeFile,
+        (python::arg("grid"), python::arg("filename"),
+         python::arg("mol") = python::object(), python::arg("confId") = -1),
+        docString.c_str());
+
+    docString =
+        "Reads Grid from a file in Gaussian CUBE format.\n\n\
+        ARGUMENTS:\n\
+        - filename:  filename of file to be read\n\
+        RETURNS:\n\
+        a tuple where the first element is the grid and\n\
+        the second element is the molecule object associated to the grid\n\
+        (only atoms and coordinates, no bonds;\n\
+        None if no molecule was associated to the grid)\n";
+    python::def("ReadFromCubeFile", readCubeFile, (python::arg("filename")),
+                docString.c_str());
+  }
+};
+}  // namespace RDMIF
+
+void wrap_mif() { mif_wrapper::wrap(); }

Code/GraphMol/MolInteractionFields/Wrap/testData/HCN.mol

@@ -0,0 +1,10 @@
+
+     RDKit          3D
+
+  3  2  0  0  0  0  0  0  0  0999 V2000
+    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1560    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0680    0.0000    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  3  0
+  1  3  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/HCl.mol

@@ -0,0 +1,8 @@
+
+     RDKit          3D
+
+  2  1  0  0  0  0  0  0  0  0999 V2000
+    0.6828    0.0000    0.0000 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6828    0.0000    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/aceticacid.mol

@@ -0,0 +1,20 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.9507   -0.1350   -0.0788 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4609    0.2321    0.2395 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3451   -0.3758   -0.5728 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8335    0.9884    1.1196 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.6234    0.3853    0.6088 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0902   -1.2118    0.0442 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1914    0.1699   -1.1002 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.2163   -0.0531   -0.2603 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  2  4  2  0
+  1  5  1  0
+  1  6  1  0
+  1  7  1  0
+  3  8  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/ethane.mol

@@ -0,0 +1,20 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.7559   -0.0040   -0.0117 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7559    0.0040    0.0117 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1297   -0.9664   -0.3736 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1568    0.1658    0.9921 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1343    0.7826   -0.6713 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1567   -0.1658   -0.9921 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1343   -0.7826    0.6713 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1297    0.9664    0.3736 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+  2  6  1  0
+  2  7  1  0
+  2  8  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/glucose.mol

@@ -0,0 +1,53 @@
+
+     RDKit          3D
+
+ 24 24  0  0  0  0  0  0  0  0999 V2000
+    3.1483    0.7685   -0.6067 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4736   -0.1497    0.2499 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0190   -0.3549   -0.1970 C   0  0  1  0  0  0  0  0  0  0  0  0
+    0.4035   -1.2513    0.7304 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9486   -1.5599    0.4098 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4301   -2.4068    1.4455 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8297   -0.2988    0.3349 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -3.1475   -0.6106   -0.1483 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2030    0.7203   -0.6126 C   0  0  1  0  0  0  0  0  0  0  0  0
+   -1.9497    1.9478   -0.5561 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2612    0.9823   -0.2638 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.7874    1.8664   -1.2622 O   0  0  0  0  0  0  0  0  0  0  0  0
+    4.0821    0.7758   -0.3230 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.5130    0.2327    1.2760 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.0077   -1.1063    0.2364 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0364   -0.8127   -1.1954 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.9954   -2.1202   -0.5323 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0743   -3.2877    1.2508 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.9558    0.1352    1.3345 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.5906    0.2559   -0.2419 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2910    0.3739   -1.6505 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.4004    2.5955   -1.0440 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.3171    1.5077    0.6975 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.7670    1.7970   -1.2063 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  3  2  1  1
+  3  4  1  0
+  4  5  1  0
+  5  6  1  0
+  5  7  1  0
+  7  8  1  6
+  7  9  1  0
+  9 10  1  6
+  9 11  1  0
+ 11 12  1  6
+ 11  3  1  0
+  1 13  1  0
+  2 14  1  0
+  2 15  1  0
+  3 16  1  0
+  5 17  1  0
+  6 18  1  0
+  7 19  1  0
+  8 20  1  0
+  9 21  1  0
+ 10 22  1  0
+ 11 23  1  0
+ 12 24  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/h2o.mol

@@ -0,0 +1,10 @@
+
+     RDKit          3D
+
+  3  2  0  0  0  0  0  0  0  0999 V2000
+    0.0075    0.3977    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7671   -0.1844    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.7596   -0.2134    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+M  END

Code/GraphMol/MolInteractionFields/Wrap/testData/test3.cube

@@ -0,0 +1,38 @@
+Gaussian cube format generated by RDKit
+*************************
+2                   0                   0                   0                   
+5                   1.88973             0                   0                   
+5                   0                   1.88973             0                   
+5                   0                   0                   1.88973             
+17                  17                  1.2903              0                   0                   
+1                   1                   -1.2903             0                   0                   
+0                   2.5                 5                   7.5                 10                  
+0.5                 3                   5.5                 8                   10.5                
+1                   3.5                 6                   8.5                 11                  
+1.5                 4                   6.5                 9                   11.5                
+2                   4.5                 7                   9.5                 12                  
+
+0.1                 2.6                 5.1                 7.6                 10.1                
+0.6                 3.1                 5.6                 8.1                 10.6                
+1.1                 3.6                 6.1                 8.6                 11.1                
+1.6                 4.1                 6.6                 9.1                 11.6                
+2.1                 4.6                 7.1                 9.6                 12.1                
+
+0.2                 2.7                 5.2                 7.7                 10.2                
+0.7                 3.2                 5.7                 8.2                 10.7                
+1.2                 3.7                 6.2                 8.7                 11.2                
+1.7                 4.2                 6.7                 9.2                 11.7                
+2.2                 4.7                 7.2                 9.7                 12.2                
+
+0.3                 2.8                 5.3                 7.8                 10.3                
+0.8                 3.3                 5.8                 8.3                 10.8                
+1.3                 3.8                 6.3                 8.8                 11.3                
+1.8                 4.3                 6.8                 9.3                 11.8                
+2.3                 4.8                 7.3                 9.8                 12.3                
+
+0.4                 2.9                 5.4                 7.9                 10.4                
+0.9                 3.4                 5.9                 8.4                 10.9                
+1.4                 3.9                 6.4                 8.9                 11.4                
+1.9                 4.4                 6.9                 9.4                 11.9                
+2.4                 4.9                 7.4                 9.9                 12.4                
+

Code/GraphMol/MolInteractionFields/Wrap/testMIF.py

@@ -0,0 +1,354 @@
+from __future__ import print_function
+import os, sys
+import unittest
+import copy
+from rdkit import DataStructs, RDConfig
+from rdkit.Geometry import rdGeometry as geom
+from rdkit.Chem import rdMIF, AllChem
+import math
+import numpy as np
+import pickle
+
+
+def feq(v1, v2, tol=1.0e-4):
+  return abs(v1 - v2) < tol
+
+
+class testfunctor():
+
+  def __call__(self, x, y, z):
+    return 1.0
+
+  def __call__(self, pt):
+    return 1.0
+
+
+class TestCase(unittest.TestCase):
+
+  def SetUp(self):
+    pass
+
+  def test1ConstructGrid(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/HCl.mol'),
+      removeHs=False)
+
+    grd = rdMIF.ConstructGrid(mol, margin=5.0, spacing=0.5)
+
+    bond = mol.GetConformer().GetAtomPosition(1) - mol.GetConformer().GetAtomPosition(0)
+    self.assertTrue(feq(grd.GetSpacing(), 0.5))
+    self.assertTrue(grd.GetNumX() == int((abs(bond.x) + 10.0) / 0.5 + 0.5))
+    self.assertTrue(grd.GetNumY() == int((abs(bond.y) + 10.0) / 0.5 + 0.5))
+    self.assertTrue(grd.GetNumZ() == int((abs(bond.z) + 10.0) / 0.5 + 0.5))
+
+  def test2CubeFiles(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/HCl.mol'),
+      removeHs=False)
+
+    grd = geom.UniformRealValueGrid3D(5.0, 5.0, 5.0, 1.0, geom.Point3D(0.0, 0.0, 0.0))
+    for i in range(grd.GetSize()):
+      grd.SetVal(i, float(i / 10.0))
+
+    rdMIF.WriteToCubeFile(
+      grd,
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/test3.cube'),
+      mol)
+
+    grd2, mol2 = rdMIF.ReadFromCubeFile(
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/test3.cube'))
+
+    self.assertTrue(grd.GetSize() == grd2.GetSize())
+
+    for i in range(grd.GetSize()):
+      self.assertTrue(feq(grd2.GetVal(i), float(i / 10.0)))
+
+    self.assertTrue(grd.CompareGrids(grd2))
+    self.assertTrue(mol.GetNumAtoms() == mol2.GetNumAtoms())
+
+    for i in range(mol.GetNumAtoms()):
+      self.assertTrue(mol.GetAtomWithIdx(i).GetAtomicNum() == mol2.GetAtomWithIdx(i).GetAtomicNum())
+      self.assertTrue(
+        feq(mol.GetConformer().GetAtomPosition(i).x,
+            mol2.GetConformer().GetAtomPosition(i).x))
+      self.assertTrue(
+        feq(mol.GetConformer().GetAtomPosition(i).y,
+            mol2.GetConformer().GetAtomPosition(i).y))
+      self.assertTrue(
+        feq(mol.GetConformer().GetAtomPosition(i).z,
+            mol2.GetConformer().GetAtomPosition(i).z))
+
+    rdMIF.WriteToCubeFile(
+      grd,
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/test4.cube'))
+
+    grd3, mol3 = rdMIF.ReadFromCubeFile(
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/test4.cube'))
+
+    self.assertTrue(grd.GetSize() == grd3.GetSize())
+
+    for i in range(grd.GetSize()):
+      self.assertTrue(feq(grd3.GetVal(i), float(i / 10.0)))
+
+    self.assertTrue(grd.CompareGrids(grd3))
+    self.assertIsNone(mol3)
+
+  def test3Coulomb(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/HCl.mol'),
+      removeHs=False)
+    AllChem.ComputeGasteigerCharges(mol)
+
+    conf = mol.GetConformer(0)
+    charges = [a.GetDoubleProp("_GasteigerCharge") for a in mol.GetAtoms()]
+    pos = conf.GetPositions()
+
+    grd = rdMIF.ConstructGrid(mol)
+
+    grd2 = copy.deepcopy(grd)
+
+    coul = rdMIF.Coulomb(mol)
+
+    rdMIF.CalculateDescriptors(grd, coul)
+    coul1 = rdMIF.Coulomb(charges, pos)
+    chargesMismatchingLen = charges[:-1]
+    self.assertRaises(ValueError, rdMIF.Coulomb, chargesMismatchingLen, pos)
+    rdMIF.CalculateDescriptors(grd2, coul1)
+    self.assertTrue(grd.CompareParams(grd2))
+    self.assertTrue(grd.CompareVectors(grd2))
+    self.assertTrue(grd.CompareGrids(grd2))
+    self.assertTrue(feq(coul(0.0, 0.0, 0.0, 1000), 0.0))
+    self.assertTrue(coul(2.0, 0.0, 0.0, 1000) < 0)
+    self.assertTrue(coul(-2.0, 0.0, 0.0, 1000) > 0)
+
+    rdMIF.CalculateDescriptors(grd, rdMIF.Coulomb(mol, absVal=True))
+
+    for i in range(grd.GetSize()):
+      self.assertTrue(grd.GetVal(i) <= 0.0)
+
+    coul1 = rdMIF.Coulomb(mol, probeCharge=-1.0)
+    self.assertTrue(coul1(-2.0, 0.0, 0.0, 1000) < 0)
+    self.assertTrue(coul1(2.0, 0.0, 0.0, 1000) > 0)
+
+    coul2 = rdMIF.Coulomb(mol, probeCharge=-.5)
+    self.assertTrue(coul1(-2.0, 0.0, 0.0, 1000) < coul2(-2.0, 0.0, 0.0, 1000))
+
+    coul3 = rdMIF.Coulomb(mol, confId=0, probeCharge=1.0, absVal=False,
+                          chargeKey="_GasteigerCharge", softcoreParam=0.01, cutoff=1.0)
+    self.assertTrue(coul3(0.0, 0.0, 0.0, 1000) > coul3(0.1, 0.0, 0.0, 1000))
+    self.assertTrue(coul3(0.66, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000))
+    self.assertTrue(coul3(0.70, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000))
+
+  def test4CoulombDielectric(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/HCl.mol'),
+      removeHs=False)
+    AllChem.ComputeGasteigerCharges(mol)
+
+    conf = mol.GetConformer(0)
+    charges = [a.GetDoubleProp("_GasteigerCharge") for a in mol.GetAtoms()]
+    pos = conf.GetPositions()
+
+    grd = rdMIF.ConstructGrid(mol, confId=0)
+    grd2 = rdMIF.ConstructGrid(mol, confId=0)
+
+    couldiele = rdMIF.CoulombDielectric(mol, confId=0)
+    couldiele1 = rdMIF.CoulombDielectric(charges, pos)
+    chargesMismatchingLen = charges[:-1]
+    self.assertRaises(ValueError, rdMIF.CoulombDielectric, chargesMismatchingLen, pos)
+
+    rdMIF.CalculateDescriptors(grd, couldiele)
+    rdMIF.CalculateDescriptors(grd2, couldiele1)
+
+    self.assertTrue(grd.CompareGrids(grd2))
+    self.assertTrue(feq(couldiele(0.0, 0.0, 0.0, 1000), 0.0))
+    self.assertTrue(couldiele(2.0, 0.0, 0.0, 1000) < 0)
+    self.assertTrue(couldiele(-2.0, 0.0, 0.0, 1000) > 0)
+
+    rdMIF.CalculateDescriptors(grd, rdMIF.CoulombDielectric(mol, absVal=True))
+
+    for i in range(grd.GetSize()):
+      self.assertTrue(grd.GetVal(i) <= 0.0)
+
+    couldiele1 = rdMIF.CoulombDielectric(mol, probeCharge=-1.0)
+    self.assertTrue(couldiele1(-2.0, 0.0, 0.0, 1000) < 0)
+    self.assertTrue(couldiele1(2.0, 0.0, 0.0, 1000) > 0)
+
+    couldiele2 = rdMIF.CoulombDielectric(mol, probeCharge=-.5)
+    self.assertTrue(couldiele1(-2.0, 0.0, 0.0, 1000) < couldiele2(-2.0, 0.0, 0.0, 1000))
+
+    couldiele3 = rdMIF.CoulombDielectric(mol, confId=0, probeCharge=1.0, absVal=False,
+                                         chargeKey="_GasteigerCharge", softcoreParam=0.01,
+                                         cutoff=1.0)
+    self.assertTrue(couldiele3(0.0, 0.0, 0.0, 1000) > couldiele3(0.1, 0.0, 0.0, 1000))
+    self.assertTrue(couldiele3(0.66, 0.0, 0.0, 1000) > couldiele3(0.68, 0.0, 0.0, 1000))
+    self.assertTrue(couldiele3(0.70, 0.0, 0.0, 1000) > couldiele3(0.68, 0.0, 0.0, 1000))
+
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/glucose.mol'), removeHs=False)
+    AllChem.ComputeGasteigerCharges(mol)
+
+    couldiele4 = rdMIF.CoulombDielectric(mol, confId=0, probeCharge=1.0, absVal=False,
+                                         chargeKey="_GasteigerCharge", softcoreParam=0.01,
+                                         cutoff=1.0, epsilon=80.0, xi=4.0)
+    couldiele5 = rdMIF.CoulombDielectric(mol, confId=0, probeCharge=1.0, absVal=False,
+                                         chargeKey="_GasteigerCharge", softcoreParam=0.01,
+                                         cutoff=1.0, epsilon=200.0, xi=4.0)
+    couldiele6 = rdMIF.CoulombDielectric(mol, confId=0, probeCharge=1.0, absVal=False,
+                                         chargeKey="_GasteigerCharge", softcoreParam=0.01,
+                                         cutoff=1.0, epsilon=80.0, xi=10.0)
+
+    self.assertTrue(couldiele5(-1.0, 0.0, 0.0, 1000) < couldiele4(-1.0, 0.0, 0.0, 1000))
+    self.assertTrue(couldiele6(-1.0, 0.0, 0.0, 1000) < couldiele4(-1.0, 0.0, 0.0, 1000))
+
+  def test5VdWaals(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/HCN.mol'),
+      removeHs=False)
+    vdw = rdMIF.MMFFVdWaals(mol, confId=0, probeAtomType=6, scaling=False, cutoff=1.0)
+
+    self.assertTrue(vdw(-5.0, 0, 0, 1000) < 0)
+    self.assertTrue(vdw(-1.68, 0, 0, 1000) > vdw(-5.0, 0, 0, 1000))
+    self.assertTrue(vdw(-5.0, 0, 0, 1000) < vdw(-10.0, 0, 0, 1000))
+
+    mol2 = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/h2o.mol'),
+      removeHs=False)
+    vdw = rdMIF.MMFFVdWaals(mol2, scaling=False)
+    vdw2 = rdMIF.MMFFVdWaals(mol2, scaling=True)
+
+    self.assertTrue(abs(vdw2(-3.0, 0, 0, 1000) - vdw(-3.0, 0, 0, 1000)) > 0.0001)
+
+    vdw3 = rdMIF.UFFVdWaals(mol, confId=0, probeAtomType="O_3", cutoff=1.0)
+
+    self.assertTrue(vdw3(-5.0, 0, 0, 1000) < 0)
+    self.assertTrue(vdw3(-1.68, 0, 0, 1000) > vdw3(-5.0, 0, 0, 1000))
+    self.assertTrue(vdw3(-5.0, 0, 0, 1000) < vdw3(-10.0, 0, 0, 1000))
+
+  def test6HBond(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/ethane.mol'), removeHs=False)
+    grd = rdMIF.ConstructGrid(mol, margin=5.0, spacing=0.5)
+    for i in range(grd.GetSize()):
+      grd.SetVal(i, 1.0)
+    grd1 = copy.deepcopy(grd)
+
+    hbonddes = rdMIF.HBond(mol)
+    rdMIF.CalculateDescriptors(grd, hbonddes)
+    self.assertTrue(not grd.CompareGrids(grd1))
+    self.assertTrue(
+      abs(int((grd.GetOccupancyVect() - grd1.GetOccupancyVect()).GetTotalVal())) == grd.GetSize())
+
+    hbonddes = rdMIF.HBond(mol, probeAtomType='O')
+    for i in range(grd.GetSize()):
+      grd.SetVal(i, 1.0)
+    rdMIF.CalculateDescriptors(grd, hbonddes)
+    self.assertTrue(not grd.CompareGrids(grd1))
+    self.assertTrue(
+      abs(int((grd.GetOccupancyVect() - grd1.GetOccupancyVect()).GetTotalVal())) == grd.GetSize())
+
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir,
+                   'Code/GraphMol/MolInteractionFields/Wrap/testData/aceticacid.mol'),
+      removeHs=False)
+    grd = rdMIF.ConstructGrid(mol, margin=5.0, spacing=0.5)
+
+    hbonddes = rdMIF.HBond(mol, probeAtomType="OH")
+    rdMIF.CalculateDescriptors(grd, hbonddes)
+
+    hbonddes1 = rdMIF.HBond(mol, probeAtomType="O", fixed=True)
+    rdMIF.CalculateDescriptors(grd, hbonddes1)
+
+    hbonddes2 = rdMIF.HBond(mol, probeAtomType="O", fixed=False)
+    self.assertTrue(hbonddes1(4.0, 0.0, 1.0, 1000) > hbonddes2(4.0, 0.0, 1.0, 1000))
+
+    hbonddes3 = rdMIF.HBond(mol, probeAtomType="NH")
+    self.assertTrue(hbonddes(2.0, 2.0, 1.0, 1000) < hbonddes3(2.0, 2.0, 1.0, 1000))
+
+    hbonddes4 = rdMIF.HBond(mol, probeAtomType="N")
+    self.assertTrue(hbonddes1(3.0, 0.0, 0.0, 1000) < hbonddes4(3.0, 0.0, 0.0, 1000))
+
+  def test7Hydrophilic(self):
+    mol = AllChem.MolFromMolFile(
+      os.path.join(RDConfig.RDBaseDir, 'Code/GraphMol/MolInteractionFields/Wrap/testData/h2o.mol'),
+      removeHs=False)
+    hydro = rdMIF.Hydrophilic(mol)
+    hbondOH = rdMIF.HBond(mol, probeAtomType="OH")
+    hbondO = rdMIF.HBond(mol, probeAtomType="O")
+
+    pt = geom.Point3D(0.0, 0.0, 0.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(1.0, 1.5, 2.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(2.0, 1.5, -3.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(-2.5, 0.5, 3.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(10.0, 1.5, 1.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(6.0, -5.0, 0.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(-3.0, -3.0, 7.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(1.0, 0.0, 0.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(0.0, 2.0, 2.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(2.0, -2.0, 0.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+    pt = geom.Point3D(2.0, -2.0, -3.0)
+    hyd = hydro(pt.x, pt.y, pt.z, 1000)
+    hOH = hbondOH(pt.x, pt.y, pt.z, 1000)
+    hO = hbondO(pt.x, pt.y, pt.z, 1000)
+    self.assertTrue(feq(min(hOH, hO), hyd))
+
+
+if __name__ == '__main__':
+  print("Testing MIF wrapper")
+  unittest.main()

Code/GraphMol/MolInteractionFields/testMIF.cpp

@@ -0,0 +1,569 @@
+//
+//  Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research and
+//  other RDKit contributors
+//
+//   @@ 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 <catch2/catch_all.hpp>
+
+#include <Geometry/UniformRealValueGrid3D.h>
+#include <Geometry/point.h>
+#include <RDGeneral/utils.h>
+#include <RDGeneral/RDLog.h>
+#include <GraphMol/MolInteractionFields/MIFDescriptors.h>
+#include <GraphMol/RDKitBase.h>
+#include <GraphMol/SmilesParse/SmilesParse.h>
+#include <GraphMol/FileParsers/FileParsers.h>
+#include <GraphMol/DistGeomHelpers/Embedder.h>
+#include <GraphMol/PartialCharges/GasteigerCharges.h>
+#include <GraphMol/PartialCharges/GasteigerParams.h>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+
+using namespace RDGeom;
+using namespace RDKit;
+using namespace RDMIF;
+
+namespace RDMIF {
+
+class testfunctor {
+ public:
+  testfunctor() {}
+  double operator()(const double &, const double &, const double &,
+                    double) const {
+    return 1.0;
+  }
+};
+}  // namespace RDMIF
+
+TEST_CASE("constructGrid") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+  //	//Generate Molecule with 3D Coordinates and Partial Charges
+  //	RWMol mol=*SmilesToMol("Cl");
+  //	MolOps::addHs(*mol);
+  //	DGeomHelpers::EmbedMolecule(*mol);
+
+  //	MolToMolFile(*mol, path + "HCl.mol");
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "HCl.mol", fopts);
+  REQUIRE(mol);
+  auto grd = *constructGrid(*mol, 0, 5.0, 0.5);
+
+  Point3D bond =
+      mol->getConformer().getAtomPos(1) - mol->getConformer().getAtomPos(0);
+  CHECK(feq(grd.getSpacing(), 0.5));
+  CHECK(feq(grd.getNumX(), std::floor((fabs(bond.x) + 10.0) / 0.5 + 0.5)));
+  CHECK(feq(grd.getNumY(), std::floor((fabs(bond.y) + 10.0) / 0.5 + 0.5)));
+  CHECK(feq(grd.getNumZ(), std::floor((fabs(bond.z) + 10.0) / 0.5 + 0.5)));
+}
+
+TEST_CASE("CalculateDescriptors") {
+  RealValueVect data(125, 0.0);
+  Point3D o(0.0, 0.0, 0.0);
+  UniformRealValueGrid3D grd(5.0, 5.0, 5.0, 1.0, &o, &data);
+
+  calculateDescriptors(grd, testfunctor());
+  CHECK(grd.getOccupancyVect()->getTotalVal() == grd.getSize());
+}
+
+TEST_CASE("CubeFiles") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "HCl.mol", fopts);
+  REQUIRE(mol);
+  RealValueVect data(125, 0.0);
+  Point3D o(0.0, 0.0, 0.0);
+  UniformRealValueGrid3D grd(5.0, 5.0, 5.0, 1.0, &o, &data);
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    grd.setVal(i, double(i / 10.0));
+  }
+  writeToCubeFile(grd, path + "test3.cube", mol.get());
+  UniformRealValueGrid3D grd2;
+  auto mol2 = readFromCubeFile(grd2, path + "test3.cube");
+
+  CHECK(grd.getSize() == grd2.getSize());
+  for (unsigned int i = 0; i < grd2.getSize(); i++) {
+    CHECK(feq(grd2.getVal(i), double(i / 10.0)));
+  }
+  CHECK(grd.getNumX() == grd2.getNumX());
+  CHECK(grd.getNumY() == grd2.getNumY());
+  CHECK(grd.getNumZ() == grd2.getNumZ());
+  CHECK(feq(grd.getOffset().x, grd2.getOffset().x));
+  CHECK(feq(grd.getSpacing(), grd2.getSpacing()));
+  CHECK(grd.compareVectors(grd2));
+  CHECK(grd.compareParams(grd2));
+  CHECK(grd.compareGrids(grd2));
+
+  CHECK(mol->getNumAtoms() == mol2->getNumAtoms());
+  for (unsigned int i = 0; i < mol->getNumAtoms(); i++) {
+    CHECK(mol->getAtomWithIdx(i)->getAtomicNum() ==
+          mol2->getAtomWithIdx(i)->getAtomicNum());
+    CHECK(feq(mol->getConformer().getAtomPos(i).x,
+              mol2->getConformer().getAtomPos(i).x));
+    CHECK(feq(mol->getConformer().getAtomPos(i).y,
+              mol2->getConformer().getAtomPos(i).y));
+    CHECK(feq(mol->getConformer().getAtomPos(i).z,
+              mol2->getConformer().getAtomPos(i).z));
+  }
+}
+
+TEST_CASE("Coulomb") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "HCl.mol", fopts);
+  REQUIRE(mol);
+
+  computeGasteigerCharges(*mol);
+
+  std::vector<double> charges;
+  std::vector<Point3D> pos;
+  Conformer conf = mol->getConformer(0);
+  for (auto i = 0u; i < mol->getNumAtoms(); ++i) {
+    charges.push_back(
+        mol->getAtomWithIdx(i)->getProp<double>("_GasteigerCharge"));
+    pos.push_back(conf.getAtomPos(i));
+  }
+
+  auto grd = *constructGrid(*mol);
+  auto grd2 = *constructGrid(*mol);
+
+  Coulomb coul(*mol);
+
+  calculateDescriptors<Coulomb>(grd, coul);
+  calculateDescriptors<Coulomb>(grd2, Coulomb(charges, pos));
+
+  CHECK(grd.compareGrids(grd2));
+  CHECK(feq(coul(0.0, 0.0, 0.0, 1000), 0.0));
+  CHECK(coul(2.0, 0.0, 0.0, 1000) < 0);
+  CHECK(coul(-2.0, 0.0, 0.0, 1000) > 0);
+  CHECK(feq(coul(0.0, 0.0, 0.0, 0.1), 0.0));
+
+  calculateDescriptors<Coulomb>(grd, Coulomb(*mol, 0, 1.0, true));
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    CHECK(grd.getVal(i) <= 0.0);
+  }
+
+  Coulomb coul1(*mol, 0, -1.0, false, "_GasteigerCharge", 0.0, 0.01);
+  CHECK(coul1(-2.0, 0.0, 0.0, 1000) < 0);
+  CHECK(coul1(2.0, 0.0, 0.0, 1000) > 0);
+
+  Coulomb coul2 = Coulomb(*mol, 0, -.5, false, "_GasteigerCharge", 0.0, 0.01);
+  CHECK(coul1(-2.0, 0.0, 0.0, 1000) < coul2(-2.0, 0.0, 0.0, 1000));
+
+  Coulomb coul3(*mol, 0, 1.0, false, "_GasteigerCharge", 0.01, 1.0);
+  CHECK(coul3(0.0, 0.0, 0.0, 1000) > coul3(0.1, 0.0, 0.0, 1000));
+  CHECK(coul3(0.66, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000));
+  CHECK(coul3(0.70, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000));
+  CHECK(feq(coul3(0.0, 0.0, 0.0, 0.1), 0.0));
+}
+
+TEST_CASE("CoulombDielectric") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "HCl.mol", fopts);
+  REQUIRE(mol);
+
+  computeGasteigerCharges(*mol);
+
+  std::vector<double> charges;
+  std::vector<Point3D> pos;
+  Conformer conf = mol->getConformer(0);
+  for (auto i = 0u; i < mol->getNumAtoms(); ++i) {
+    charges.push_back(
+        mol->getAtomWithIdx(i)->getProp<double>("_GasteigerCharge"));
+    pos.push_back(conf.getAtomPos(i));
+  }
+
+  auto grd = *constructGrid(*mol);
+  auto grd2 = *constructGrid(*mol);
+
+  CoulombDielectric couldiele(*mol);
+
+  calculateDescriptors<CoulombDielectric>(grd, couldiele);
+  calculateDescriptors<CoulombDielectric>(grd2,
+                                          CoulombDielectric(charges, pos));
+
+  CHECK(grd.compareGrids(grd2));
+  CHECK(feq(couldiele(0.0, 0.0, 0.0, 1000), 0.0));
+  CHECK(couldiele(2.0, 0.0, 0.0, 1000) < 0);
+  CHECK(couldiele(-2.0, 0.0, 0.0, 1000) > 0);
+
+  calculateDescriptors<CoulombDielectric>(
+      grd, CoulombDielectric(*mol, 0, 1.0, true));
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    CHECK(grd.getVal(i) <= 0.0);
+  }
+
+  CHECK(feq(couldiele(0.0, 0.0, 0.0, 1000), 0.0));
+  CHECK(couldiele(2.0, 0.0, 0.0, 1000) < 0);
+  CHECK(couldiele(-2.0, 0.0, 0.0, 1000) > 0);
+
+  calculateDescriptors<CoulombDielectric>(
+      grd, CoulombDielectric(*mol, 0, 1.0, true));
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    CHECK(grd.getVal(i) <= 0.0);
+  }
+
+  CoulombDielectric couldiele1(*mol, 0, -1.0);
+  CHECK(couldiele1(-2.0, 0.0, 0.0, 1000) < 0);
+  CHECK(couldiele1(2.0, 0.0, 0.0, 1000) > 0);
+
+  CoulombDielectric couldiele2(*mol, 0, -.5);
+
+  CHECK(couldiele1(-2.0, 0.0, 0.0, 1000) < couldiele2(-2.0, 0.0, 0.0, 1000));
+
+  CoulombDielectric couldiele3(*mol, 0, 1.0, false, "_GasteigerCharge", 0.01,
+                               1.0);
+  CHECK(couldiele3(0.0, 0.0, 0.0, 1000) > couldiele3(0.1, 0.0, 0.0, 1000));
+  CHECK(couldiele3(0.66, 0.0, 0.0, 1000) > couldiele3(0.68, 0.0, 0.0, 1000));
+  CHECK(couldiele3(0.70, 0.0, 0.0, 1000) > couldiele3(0.68, 0.0, 0.0, 1000));
+
+  mol = v2::FileParsers::MolFromMolFile(path + "glucose.mol", fopts);
+  REQUIRE(mol);
+  computeGasteigerCharges(*mol);
+
+  CoulombDielectric couldiele4(*mol, 0, 1.0, false, "_GasteigerCharge", 0.01,
+                               1.0, 80.0, 4.0);
+  CoulombDielectric couldiele5(*mol, 0, 1.0, false, "_GasteigerCharge", 0.01,
+                               1.0, 200.0, 4.0);
+  CoulombDielectric couldiele6(*mol, 0, 1.0, false, "_GasteigerCharge", 0.01,
+                               1.0, 80.0, 10.0);
+  CHECK(couldiele5(-1.0, 0.0, 0.0, 1000) < couldiele4(-1.0, 0.0, 0.0, 1000));
+
+  CHECK(couldiele6(-1.0, 0.0, 0.0, 1000) < couldiele4(-1.0, 0.0, 0.0, 1000));
+}
+
+TEST_CASE("VdWaals") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "HCl.mol", fopts);
+  REQUIRE(mol);
+  try {
+    MMFFVdWaals vdw(*mol, 0, 6, false, 1.0);
+  } catch (ValueErrorException &dexp) {
+    BOOST_LOG(rdInfoLog) << "Expected failure: " << dexp.what() << "\n";
+  }
+
+  mol = v2::FileParsers::MolFromMolFile(path + "HCN.mol", fopts);
+  REQUIRE(mol);
+  {
+    MMFFVdWaals vdw(*mol, 0, 6, false, 1.0);
+
+    CHECK(vdw(-5.0, 0, 0, 1000) < 0);
+    CHECK(vdw(-1.68, 0, 0, 1000) > vdw(-5.0, 0, 0, 1000));
+    CHECK(vdw(-5.0, 0, 0, 1000) < vdw(-10.0, 0, 0, 1000));
+  }
+
+  auto mol2 = v2::FileParsers::MolFromMolFile(path + "h2o.mol", fopts);
+  REQUIRE(mol2);
+  MMFFVdWaals vdw(*mol2, 0, 6, false, 1.0);
+  MMFFVdWaals vdw2(*mol2, 0, 6, true, 1.0);
+  CHECK(fabs(vdw2(-3.0, 0, 0, 1000) - vdw(-3.0, 0, 0, 1000)) > 0.0001);
+
+  UFFVdWaals vdw3(*mol, 0, "O_3", 1.0);
+  CHECK(vdw3(-5.0, 0, 0, 1000) < 0);
+  CHECK(vdw3(-1.68, 0, 0, 1000) > vdw3(-5.0, 0, 0, 1000));
+  CHECK(vdw3(-5.0, 0, 0, 1000) < vdw3(-10.0, 0, 0, 1000));
+
+  std::vector<std::string> names{
+      "acetone",       "aceticacid",    "phenol",       "phenolate",
+      "serine",        "threonine",     "ethanol",      "diethylether",
+      "h2o",           "ammonia",       "ethylamine",   "imine",
+      "acetonitrile",  "histidine",     "phenylamine",  "methylammonium",
+      "fluoromethane", "chloromethane", "bromomethane", "glycine",
+      "glyphe",        "glysergly",     "glythrgly",    "glucose"};
+  for (const auto &name : names) {
+    mol = v2::FileParsers::MolFromMolFile(path + name + ".mol", fopts);
+    REQUIRE(mol);
+    MMFFVdWaals mmffVdw(*mol);
+    CHECK(mmffVdw(0.0, 0.0, 0.0, 1000));
+    UFFVdWaals uffVdw(*mol);
+    CHECK(uffVdw(0.0, 0.0, 0.0, 1000));
+  }
+}
+
+TEST_CASE("HBond") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  // Generate Molecule with 3D Coordinates and Partial Charges
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol =
+      v2::FileParsers::MolFromMolFile(path + "ethane.mol", fopts);  // Ethane
+  REQUIRE(mol);
+  auto grd = *constructGrid(*mol, 0, 5.0, 1);
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    grd.setVal(i, 1.0);
+  }
+
+  UniformRealValueGrid3D grd1(grd);
+
+  HBond hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  TEST_ASSERT(!grd.compareGrids(grd1));
+  const RealValueVect *vect = grd.getOccupancyVect();
+  const RealValueVect *vect1 = grd1.getOccupancyVect();
+  CHECK((unsigned int)fabs(((*vect1 - *vect).getTotalVal())) == grd.getSize());
+
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+
+  for (unsigned int i = 0; i < grd.getSize(); i++) {
+    grd.setVal(i, 1.0);
+  }
+  calculateDescriptors<HBond>(grd, hbonddes);
+  TEST_ASSERT(!grd.compareGrids(grd1));
+  CHECK((unsigned int)fabs(((*vect1 - *vect).getTotalVal())));
+
+  mol = v2::FileParsers::MolFromMolFile(path + "aceticacid.mol",
+                                        fopts);  // Acetic Acid
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 2);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  HBond hbonddes1(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes1.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes1);
+  HBond hbonddes2(*mol, 0, "O", false, 0.001);
+  CHECK(hbonddes1(4.0, 0.0, 1.0, 1000) > hbonddes2(4.0, 0.0, 1.0, 1000));
+  HBond hbonddes3(*mol, 0, "NH", true, 0.001);
+  CHECK(hbonddes3.getNumInteractions() == 2);
+  CHECK(hbonddes(2.0, 2.0, 2.0, 1000) < hbonddes3(2.0, 2.0, 2.0, 1000));
+  HBond hbonddes4(*mol, 0, "N", true, 0.001);
+  CHECK(hbonddes4.getNumInteractions() == 1);
+  CHECK(hbonddes1(3.0, 0.0, 0.0, 1000) < hbonddes4(3.0, 0.0, 0.0, 1000));
+
+  mol =
+      v2::FileParsers::MolFromMolFile(path + "acetone.mol", fopts);  // Acetone
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "diethylether.mol",
+                                        fopts);  // Et2O
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "h2o.mol", fopts);  // H2O
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 2);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol =
+      v2::FileParsers::MolFromMolFile(path + "ammonia.mol", fopts);  // ammonia
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 3);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "imine.mol", fopts);  // imine
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "methylammonium.mol",
+                                        fopts);  // methylammonium
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 3);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "chloromethane.mol",
+                                        fopts);  // Chloromethane
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "phosphonate.mol",
+                                        fopts);  // Phosphonate
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 3);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "phosphatediester.mol",
+                                        fopts);  // Phosphatediester
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 4);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "hydrogenphosphatediester.mol",
+                                        fopts);  // Hydrogenphosphatediester
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 4);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "mustardgas.mol",
+                                        fopts);  // mustard gas
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 2);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "alicin.mol", fopts);  // Alicin
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 1);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 0);
+  calculateDescriptors<HBond>(grd, hbonddes);
+
+  mol = v2::FileParsers::MolFromMolFile(path + "sulfanilamide.mol",
+                                        fopts);  // Sulfanilamide
+  REQUIRE(mol);
+  grd = *constructGrid(*mol, 0, 5.0, 1);
+  hbonddes = HBond(*mol, 0, "OH", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 3);
+  calculateDescriptors<HBond>(grd, hbonddes);
+  hbonddes = HBond(*mol, 0, "O", true, 0.001);
+  CHECK(hbonddes.getNumInteractions() == 4);
+  calculateDescriptors<HBond>(grd, hbonddes);
+}
+
+TEST_CASE("Hydrophilic") {
+  std::string path = getenv("RDBASE");
+  path += "/Code/GraphMol/MolInteractionFields/test_data/";
+
+  auto fopts = v2::FileParsers::MolFileParserParams();
+  fopts.removeHs = false;
+  auto mol = v2::FileParsers::MolFromMolFile(path + "h2o.mol", fopts);
+  REQUIRE(mol);
+
+  Hydrophilic hydro(*mol);
+  HBond hbondOH(*mol, 0, "OH");
+  HBond hbondO(*mol, 0, "O");
+
+  double hyd = hydro(0.0, 0.0, 0.0, 1000), hOH = hbondOH(0.0, 0.0, 0.0, 1000),
+         hO = hbondO(0.0, 0.0, 0.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(1.0, 1.5, 2.0, 1000);
+  hOH = hbondOH(1.0, 1.5, 2.0, 1000);
+  hO = hbondO(1.0, 1.5, 2.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(2.0, 1.5, -3.0, 1000);
+  hOH = hbondOH(2.0, 1.5, -3.0, 1000);
+  hO = hbondO(2.0, 1.5, -3.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(-2.5, 0.5, 3.0, 1000);
+  hOH = hbondOH(-2.5, 0.5, 3.0, 1000);
+  hO = hbondO(-2.5, 0.5, 3.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(10.0, 1.5, 1.0, 1000);
+  hOH = hbondOH(10.0, 1.5, 1.0, 1000);
+  hO = hbondO(10.0, 1.5, 1.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(6.0, -5.0, 0.0, 1000);
+  hOH = hbondOH(6.0, -5.0, 0.0, 1000);
+  hO = hbondO(6.0, -5.0, 0.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(-3.0, -3.0, 7.0, 1000);
+  hOH = hbondOH(-3.0, -3.0, 7.0, 1000);
+  hO = hbondO(-3.0, -3.0, 7.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(1.0, 0.0, 0.0, 1000);
+  hOH = hbondOH(1.0, 0.0, 0.0, 1000);
+  hO = hbondO(1.0, 0.0, 0.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(0.0, 2.0, 2.0, 1000);
+  hOH = hbondOH(0.0, 2.0, 2.0, 1000);
+  hO = hbondO(0.0, 2.0, 2.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(2.0, -2.0, 0.0, 1000);
+  hOH = hbondOH(2.0, -2.0, 0.0, 1000);
+  hO = hbondO(2.0, -2.0, 0.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+
+  hyd = hydro(2.0, -2.0, -3.0, 1000);
+  hOH = hbondOH(2.0, -2.0, -3.0, 1000);
+  hO = hbondO(2.0, -2.0, -3.0, 1000);
+  CHECK(feq(std::min(hOH, hO), hyd));
+}

Code/GraphMol/MolInteractionFields/test_data/HCN.mol

@@ -0,0 +1,10 @@
+
+     RDKit          3D
+
+  3  2  0  0  0  0  0  0  0  0999 V2000
+    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1560    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0680    0.0000    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  3  0
+  1  3  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/HCl.mol

@@ -0,0 +1,8 @@
+
+     RDKit          3D
+
+  2  1  0  0  0  0  0  0  0  0999 V2000
+    0.6828    0.0000    0.0000 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6828    0.0000    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/aceticacid.mol

@@ -0,0 +1,20 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.9507   -0.1350   -0.0788 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4609    0.2321    0.2395 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3451   -0.3758   -0.5728 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8335    0.9884    1.1196 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.6234    0.3853    0.6088 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0902   -1.2118    0.0442 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1914    0.1699   -1.1002 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.2163   -0.0531   -0.2603 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  2  4  2  0
+  1  5  1  0
+  1  6  1  0
+  1  7  1  0
+  3  8  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/acetone.mol

@@ -0,0 +1,12 @@
+
+     RDKit          3D
+
+  4  3  0  0  0  0  0  0  0  0999 V2000
+   -0.0121    0.0021    0.0006 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1926   -0.2084   -0.0633 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7501   -0.8149    0.0408 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4303    1.0212    0.0218 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  1  3  1  0
+  1  4  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/acetonitrile.mol

@@ -0,0 +1,16 @@
+
+     RDKit          3D
+
+  6  5  0  0  0  0  0  0  0  0999 V2000
+   -0.4860   -0.0121    0.0118 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.9750    0.0242   -0.0236 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.1346    0.0530   -0.0516 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8807   -0.5958   -0.8251 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.8412   -0.4672    0.9411 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.9016    0.9979   -0.0525 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  3  0
+  1  4  1  0
+  1  5  1  0
+  1  6  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/alicin.mol

@@ -0,0 +1,42 @@
+
+     RDKit          3D
+
+ 19 18  0  0  0  0  0  0  0  0999 V2000
+    0.6440   -1.2602   -0.8499 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8855   -0.5513    0.4503 S   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9358    0.3975    0.6790 S   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.0943   -0.8949    0.1500 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.4951   -0.4536    0.4583 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.4404   -0.2335   -0.4658 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9854    0.8708    0.1332 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.3865    0.4018   -0.1274 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.4551    0.7529    0.6004 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9676   -1.0842   -0.9209 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.8859   -1.8242    0.6901 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.7501   -0.3259    1.5085 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.2407   -0.3523   -1.5266 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -5.4393    0.0779   -0.1757 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.9393    1.5323    1.0044 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.6101    1.4133   -0.7404 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.5232   -0.2608   -0.9815 H   0  0  0  0  0  1  0  0  0  0  0  0
+    4.3769    1.4163    1.4564 H   0  0  0  0  0  1  0  0  0  0  0  0
+    5.4431    0.3781    0.3480 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  2  3  1  0
+  3  4  1  0
+  4  5  1  0
+  5  6  2  3
+  2  7  1  0
+  7  8  1  0
+  8  9  2  3
+  4 10  1  0
+  4 11  1  0
+  5 12  1  0
+  6 13  1  0
+  6 14  1  0
+  7 15  1  0
+  7 16  1  0
+  8 17  1  0
+  9 18  1  0
+  9 19  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/ammonia.mol

@@ -0,0 +1,12 @@
+
+     RDKit          3D
+
+  4  3  0  0  0  0  0  0  0  0999 V2000
+    0.0064   -0.0023    0.2955 N   0  0  0  0  0  0  0  0  0  0  0  0
+    0.9211   -0.1732   -0.1197 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.6144   -0.7120   -0.0909 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.3130    0.8874   -0.0850 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/bromomethane.mol

@@ -0,0 +1,14 @@
+
+     RDKit          3D
+
+  5  4  0  0  0  0  0  0  0  0999 V2000
+   -0.1838   -0.0007    0.0043 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.7506    0.0063   -0.0407 Br  0  0  0  0  0  0  0  0  0  0  0  0
+   -0.5410   -0.5342   -0.8778 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.5263    1.0351   -0.0037 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4993   -0.5066    0.9179 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/chloromethane.mol

@@ -0,0 +1,14 @@
+
+     RDKit          3D
+
+  5  4  0  0  0  0  0  0  0  0999 V2000
+   -0.1389    0.0007    0.0088 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.6245   -0.0088   -0.1026 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+   -0.5553   -0.4796   -0.8793 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4431   -0.5465    0.9037 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4872    1.0341    0.0694 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/diethylether.mol

@@ -0,0 +1,22 @@
+
+     RDKit          3D
+
+  9  8  0  0  0  0  0  0  0  0999 V2000
+    1.1671   -0.0430    0.1359 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0465    0.8308    0.1212 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1625    0.1252   -0.1239 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2774   -0.5409   -0.8323 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0628   -0.7863    0.9323 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.0648    0.5508    0.3277 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.9868    0.8434   -0.1243 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.3419   -0.6127    0.6641 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1274   -0.3672   -1.1006 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  1  4  1  0
+  1  5  1  0
+  1  6  1  0
+  3  7  1  0
+  3  8  1  0
+  3  9  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/ethane.mol

@@ -0,0 +1,20 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.7559   -0.0040   -0.0117 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7559    0.0040    0.0117 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1297   -0.9664   -0.3736 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1568    0.1658    0.9921 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1343    0.7826   -0.6713 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1567   -0.1658   -0.9921 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1343   -0.7826    0.6713 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1297    0.9664    0.3736 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+  2  6  1  0
+  2  7  1  0
+  2  8  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/ethanol.mol

@@ -0,0 +1,22 @@
+
+     RDKit          3D
+
+  9  8  0  0  0  0  0  0  0  0999 V2000
+   -0.9024    0.0481    0.0321 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5429   -0.3442    0.2584 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3428    0.1039   -0.8261 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0006    1.1358   -0.0457 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2733   -0.3731   -0.9081 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.5350   -0.3046    0.8515 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.6383   -1.4318    0.3261 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.9284    0.0943    1.1838 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.2590    1.0716   -0.8719 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  1  4  1  0
+  1  5  1  0
+  1  6  1  0
+  2  7  1  0
+  2  8  1  0
+  3  9  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/ethylamine.mol

@@ -0,0 +1,24 @@
+
+     RDKit          3D
+
+ 10  9  0  0  0  0  0  0  0  0999 V2000
+   -0.9814    0.1961    0.2708 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1904   -0.6089   -0.2649 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4241    0.1665   -0.2251 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8979   -0.4014    0.2332 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1447    1.1018   -0.3229 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.8152    0.4936    1.3116 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.3182   -1.5200    0.3285 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.0094   -0.9153   -1.2968 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.6090    0.4651    0.7319 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.3070    1.0225   -0.7663 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  1  4  1  0
+  1  5  1  0
+  1  6  1  0
+  2  7  1  0
+  2  8  1  0
+  3  9  1  0
+  3 10  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/fluoromethane.mol

@@ -0,0 +1,14 @@
+
+     RDKit          3D
+
+  5  4  0  0  0  0  0  0  0  0999 V2000
+   -0.0603    0.0040    0.0067 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2853   -0.0850   -0.1424 F   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3051   -0.1174    1.0638 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.3884    0.9839   -0.3459 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.5315   -0.7855   -0.5822 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/glucose.mol

@@ -0,0 +1,53 @@
+
+     RDKit          3D
+
+ 24 24  0  0  0  0  0  0  0  0999 V2000
+    3.1483    0.7685   -0.6067 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4736   -0.1497    0.2499 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0190   -0.3549   -0.1970 C   0  0  1  0  0  0  0  0  0  0  0  0
+    0.4035   -1.2513    0.7304 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9486   -1.5599    0.4098 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4301   -2.4068    1.4455 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8297   -0.2988    0.3349 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -3.1475   -0.6106   -0.1483 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2030    0.7203   -0.6126 C   0  0  1  0  0  0  0  0  0  0  0  0
+   -1.9497    1.9478   -0.5561 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2612    0.9823   -0.2638 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.7874    1.8664   -1.2622 O   0  0  0  0  0  0  0  0  0  0  0  0
+    4.0821    0.7758   -0.3230 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.5130    0.2327    1.2760 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.0077   -1.1063    0.2364 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0364   -0.8127   -1.1954 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.9954   -2.1202   -0.5323 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0743   -3.2877    1.2508 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.9558    0.1352    1.3345 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.5906    0.2559   -0.2419 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2910    0.3739   -1.6505 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.4004    2.5955   -1.0440 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.3171    1.5077    0.6975 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.7670    1.7970   -1.2063 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  3  2  1  1
+  3  4  1  0
+  4  5  1  0
+  5  6  1  0
+  5  7  1  0
+  7  8  1  6
+  7  9  1  0
+  9 10  1  6
+  9 11  1  0
+ 11 12  1  6
+ 11  3  1  0
+  1 13  1  0
+  2 14  1  0
+  2 15  1  0
+  3 16  1  0
+  5 17  1  0
+  6 18  1  0
+  7 19  1  0
+  8 20  1  0
+  9 21  1  0
+ 10 22  1  0
+ 11 23  1  0
+ 12 24  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/glycine.mol

@@ -0,0 +1,24 @@
+
+     RDKit          3D
+
+ 10  9  0  0  0  0  0  0  0  0999 V2000
+    0.2930   -0.4542   -0.2562 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7942    0.3664    0.4374 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6862    1.5456    0.7433 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9158   -0.3486    0.6887 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4946    0.3556   -0.5124 N   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5618   -1.2981    0.3856 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.0975   -0.8242   -1.2086 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.8682   -1.2777    0.3928 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.2109    1.1799   -1.0495 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.8017    0.7553    0.3789 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  2  0
+  2  4  1  0
+  1  5  1  0
+  1  6  1  0
+  1  7  1  0
+  4  8  1  0
+  5  9  1  0
+  5 10  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/glyphe.mol

@@ -0,0 +1,65 @@
+
+     RDKit          3D
+
+ 30 30  0  0  0  0  0  0  0  0999 V2000
+    3.8045   -0.7137    0.8027 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.6586    0.1298   -0.1314 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.3179    0.6232   -1.1964 O   0  0  0  0  0  0  0  0  0  0  0  0
+    5.9273    0.2643    0.3073 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4098   -0.6463    0.4120 N   0  0  0  0  0  0  0  0  0  0  0  0
+    1.6109    0.3949    0.8392 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9158    1.1601    1.7480 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2696    0.4779    0.0885 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -0.7605   -0.3932    0.8247 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.1703   -0.2019    0.3140 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7418   -1.1324   -0.5642 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.0366   -0.9436   -1.0493 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.7710    0.1757   -0.6622 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.2121    1.1069    0.2112 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.9171    0.9206    0.6971 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4743    0.0323   -1.3205 N   0  0  0  0  0  0  0  0  0  0  0  0
+    3.9201   -0.3529    1.8288 H   0  0  0  0  0  1  0  0  0  0  0  0
+    4.1203   -1.7572    0.7254 H   0  0  0  0  0  1  0  0  0  0  0  0
+    6.3550    0.8275   -0.3715 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.1495   -0.9848   -0.5142 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.0368    1.5304    0.0717 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4817   -1.4540    0.7608 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.7480   -0.1619    1.8985 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.1798   -2.0103   -0.8761 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.4735   -1.6718   -1.7278 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -5.7801    0.3208   -1.0386 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.7858    1.9781    0.5161 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.4960    1.6550    1.3805 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0826    0.7164   -1.7788 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4252    0.1101   -1.8013 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  2  0
+  2  4  1  0
+  1  5  1  0
+  5  6  1  0
+  6  7  2  0
+  8  6  1  1
+  8  9  1  0
+  9 10  1  0
+ 10 11  2  0
+ 11 12  1  0
+ 12 13  2  0
+ 13 14  1  0
+ 14 15  2  0
+  8 16  1  0
+ 15 10  1  0
+  1 17  1  0
+  1 18  1  0
+  4 19  1  0
+  5 20  1  0
+  8 21  1  0
+  9 22  1  0
+  9 23  1  0
+ 11 24  1  0
+ 12 25  1  0
+ 13 26  1  0
+ 14 27  1  0
+ 15 28  1  0
+ 16 29  1  0
+ 16 30  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/glysergly.mol

@@ -0,0 +1,60 @@
+
+     RDKit          3D
+
+ 28 27  0  0  0  0  0  0  0  0999 V2000
+    3.6000    2.8791   -0.8706 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.8997    1.9897   -1.6008 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.0832    1.8793   -2.8023 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.8868    1.2239   -0.7685 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3215    0.1220   -1.5333 N   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3173   -0.7228   -1.1460 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.2733   -1.3623   -2.0164 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0167   -0.8395    0.3572 C   0  0  2  0  0  0  0  0  0  0  0  0
+    1.0408   -1.6961    1.0721 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8787   -3.0781    0.7545 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3472   -1.4221    0.5447 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.5533   -0.7883    0.4747 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.5925   -1.4488    0.4569 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.5589    0.7451    0.4581 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9347    1.2645    0.3976 N   0  0  0  0  0  0  0  0  0  0  0  0
+    4.1947    3.3133   -1.5201 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.3756    0.8168    0.1196 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0750    1.8976   -0.4801 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.4385    0.1605   -2.5448 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.0181    0.1662    0.7869 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.0598   -1.4023    0.8038 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.9232   -1.5952    2.1564 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.5228   -3.5813    1.2831 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.4040   -2.4247    0.3841 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.0787    1.1177    1.3670 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.0077    1.1014   -0.4172 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.4560    0.8336    1.1668 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.3765    0.8507   -0.4290 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  2  0
+  2  4  1  0
+  4  5  1  0
+  5  6  1  0
+  6  7  2  0
+  8  6  1  6
+  8  9  1  0
+  9 10  1  0
+  8 11  1  0
+ 11 12  1  0
+ 12 13  2  0
+ 12 14  1  0
+ 14 15  1  0
+  1 16  1  0
+  4 17  1  0
+  4 18  1  0
+  5 19  1  0
+  8 20  1  0
+  9 21  1  0
+  9 22  1  0
+ 10 23  1  0
+ 11 24  1  0
+ 14 25  1  0
+ 14 26  1  0
+ 15 27  1  0
+ 15 28  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/glythrgly.mol

@@ -0,0 +1,68 @@
+
+     RDKit          3D
+
+ 32 31  0  0  0  0  0  0  0  0999 V2000
+    5.0198    2.4382    0.8430 O   0  0  0  0  0  0  0  0  0  0  0  0
+    3.7744    2.5576    1.3496 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.5085    3.4675    2.1199 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.8175    1.5160    0.7945 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.6824    1.3419    1.7099 N   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6485    0.6698    0.9205 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1727   -0.4587    1.5532 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5281   -0.8798    2.6519 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8740   -1.1768    0.6878 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -0.3502   -2.5992    0.4023 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.6037   -2.5404   -0.6729 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4746   -3.5426   -0.0123 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1941   -0.5029   -0.5777 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.0738    0.5494   -0.6373 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.6893    0.9536    0.3456 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.2044    1.1717   -2.0321 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.2221    2.2321   -2.0381 N   0  0  0  0  0  0  0  0  0  0  0  0
+    5.5191    3.1601    1.2806 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.3453    0.5623    0.6958 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.4885    1.8618   -0.1914 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.2549    2.2509    1.9043 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.7864   -1.2121    1.2958 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.1589   -3.0177    1.2779 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.9905   -3.4319   -0.7577 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0767   -4.5433   -0.2129 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.2340   -3.6228    0.7715 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.9553   -3.2078   -0.9379 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.5666   -0.6952   -1.3530 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2369    1.5925   -2.3231 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.4848    0.3925   -2.7478 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.9877    2.8919   -1.2908 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -4.1020    1.8213   -1.7155 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  2  0
+  2  4  1  0
+  4  5  1  0
+  5  6  1  0
+  6  7  1  0
+  7  8  2  0
+  9  7  1  1
+  9 10  1  0
+ 10 11  1  6
+ 10 12  1  0
+  9 13  1  0
+ 13 14  1  0
+ 14 15  2  0
+ 14 16  1  0
+ 16 17  1  0
+  1 18  1  0
+  4 19  1  0
+  4 20  1  0
+  5 21  1  0
+  9 22  1  0
+ 10 23  1  0
+ 11 24  1  0
+ 12 25  1  0
+ 12 26  1  0
+ 12 27  1  0
+ 13 28  1  0
+ 16 29  1  0
+ 16 30  1  0
+ 17 31  1  0
+ 17 32  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/h2o.mol

@@ -0,0 +1,10 @@
+
+     RDKit          3D
+
+  3  2  0  0  0  0  0  0  0  0999 V2000
+    0.0075    0.3977    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7671   -0.1844    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.7596   -0.2134    0.0000 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/histidine.mol

@@ -0,0 +1,45 @@
+
+     RDKit          3D
+
+ 20 20  0  0  0  0  0  0  0  0999 V2000
+    2.1562   -1.6665   -0.0773 N   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1257   -0.6303    0.2028 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1156   -0.8754   -0.6728 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3492   -0.2257   -0.1368 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.4628   -0.7372    0.4940 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.3487    0.2604    0.8082 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7936    1.3684    0.3752 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5854    1.1264   -0.2109 N   0  0  0  0  0  0  0  0  0  0  0  0
+    1.7179    0.7695   -0.0507 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1262    1.7660   -0.4383 O   0  0  0  0  0  0  0  0  0  0  0  0
+    3.0400    0.8383    0.2350 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4249   -1.6132   -1.0600 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.9998   -1.3950    0.4373 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.8691   -0.6882    1.2670 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.0608   -0.5319   -1.7005 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.3120   -1.9545   -0.7364 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.6731   -1.7705    0.7401 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.2104    2.3629    0.4591 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.9507    1.8229   -0.5874 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.2807    1.7735    0.0651 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  3  4  1  0
+  4  5  2  0
+  5  6  1  0
+  6  7  2  0
+  7  8  1  0
+  2  9  1  0
+  9 10  2  0
+  9 11  1  0
+  8  4  1  0
+  1 12  1  0
+  1 13  1  0
+  2 14  1  0
+  3 15  1  0
+  3 16  1  0
+  5 17  1  0
+  7 18  1  0
+  8 19  1  0
+ 11 20  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/hydrogenphosphatediester.mol

@@ -0,0 +1,38 @@
+
+     RDKit          3D
+
+ 17 16  0  0  0  0  0  0  0  0999 V2000
+   -2.7923   -0.7921   -0.4606 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3878   -0.7696   -0.3155 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7787    0.3738    0.6309 P   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4381    0.5401    1.9596 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9117    1.7041   -0.2552 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8046    0.1224    0.6694 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4770   -0.1306   -0.5512 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.9680   -0.1634   -0.2844 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.2649   -1.0223    0.4980 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.0585   -1.5697   -1.1810 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.1529    0.1712   -0.8314 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.9020    2.4848    0.3287 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1495   -1.0941   -0.9533 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.2453    0.6548   -1.2770 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.5250   -0.3614   -1.2044 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.3056    0.7903    0.1348 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.2117   -0.9382    0.4500 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  3  4  2  0
+  3  5  1  0
+  3  6  1  0
+  6  7  1  0
+  7  8  1  0
+  1  9  1  0
+  1 10  1  0
+  1 11  1  0
+  5 12  1  0
+  7 13  1  0
+  7 14  1  0
+  8 15  1  0
+  8 16  1  0
+  8 17  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/imine.mol

@@ -0,0 +1,44 @@
+
+     RDKit          3D
+
+ 20 19  0  0  0  0  0  0  0  0999 V2000
+    1.8530    0.6711    1.0248 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3387    0.5832    0.9585 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1045   -0.0974   -0.3254 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4269    0.5039   -1.4202 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3803    1.9550   -1.4857 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1747   -1.6153   -0.3181 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5363   -2.0814    0.1634 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.2558    1.2416    0.1821 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.1652    1.1624    1.9515 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.3088   -0.3240    0.9971 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.0326    0.0015    1.8092 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.1031    1.5796    1.0428 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.5996    2.3504   -1.2021 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.5664    2.2691   -2.5177 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.1598    2.4053   -0.8631 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.0095   -1.9874   -1.3321 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.6131   -2.0128    0.3294 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.5898   -3.1743    0.1557 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.3340   -1.6942   -0.4792 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.7353   -1.7363    1.1832 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  3  4  2  3
+  4  5  1  0
+  3  6  1  0
+  6  7  1  0
+  1  8  1  0
+  1  9  1  0
+  1 10  1  0
+  2 11  1  0
+  2 12  1  0
+  5 13  1  0
+  5 14  1  0
+  5 15  1  0
+  6 16  1  0
+  6 17  1  0
+  7 18  1  0
+  7 19  1  0
+  7 20  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/methylammonium.mol

@@ -0,0 +1,21 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.7491    0.0327   -0.0707 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7187   -0.0313    0.0678 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1887   -0.3663    0.8454 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0255    1.0787   -0.2156 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0295   -0.5710   -0.9359 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0400   -0.9983    0.2076 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.0437    0.5220    0.8715 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.1904    0.3335   -0.7702 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  4  1  0
+  1  5  1  0
+  2  6  1  0
+  2  7  1  0
+  2  8  1  0
+M  CHG  1   2   1
+M  END

Code/GraphMol/MolInteractionFields/test_data/mustardgas.mol

@@ -0,0 +1,34 @@
+
+     RDKit          3D
+
+ 15 14  0  0  0  0  0  0  0  0999 V2000
+    4.0857   -0.2378   -0.1321 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+    2.6553    0.5677    0.5658 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3898   -0.1730    0.1579 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0755    0.6378    0.8801 S   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3630   -0.4651    0.1975 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7741   -0.0026    0.5483 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.2448    1.4640   -0.3487 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+    2.7868    0.5755    1.6524 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.6545    1.5999    0.2019 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.4242   -1.2096    0.5094 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.2870   -0.1773   -0.9324 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2016   -1.4615    0.6218 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.2429   -0.5352   -0.8888 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.4940   -0.7796    0.2732 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.8874    0.1967    1.6188 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  3  4  1  0
+  4  5  1  0
+  5  6  1  0
+  6  7  1  0
+  2  8  1  0
+  2  9  1  0
+  3 10  1  0
+  3 11  1  0
+  5 12  1  0
+  5 13  1  0
+  6 14  1  0
+  6 15  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/phenol.mol

@@ -0,0 +1,31 @@
+
+     RDKit          3D
+
+ 13 13  0  0  0  0  0  0  0  0999 V2000
+    0.1997    1.2435   -0.4769 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1699    0.9817   -0.4461 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.6331   -0.2518    0.0130 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7267   -1.2233    0.4413 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6445   -0.9633    0.4113 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0960    0.2707   -0.0486 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4261    0.5677   -0.0952 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5667    2.2016   -0.8331 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.8741    1.7390   -0.7799 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.7005   -0.4558    0.0369 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0904   -2.1833    0.7984 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.3376   -1.7276    0.7472 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.9241   -0.1992    0.2315 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  2  3  1  0
+  3  4  2  0
+  4  5  1  0
+  5  6  2  0
+  6  7  1  0
+  6  1  1  0
+  1  8  1  0
+  2  9  1  0
+  3 10  1  0
+  4 11  1  0
+  5 12  1  0
+  7 13  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/phenolate.mol

@@ -0,0 +1,30 @@
+
+     RDKit          3D
+
+ 12 12  0  0  0  0  0  0  0  0999 V2000
+    0.2382    1.3412   -0.0165 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0854    0.8829    0.0615 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3477   -0.4803    0.0794 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.2845   -1.3711    0.0192 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0306   -0.8890   -0.0583 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3508    0.4814   -0.0796 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5500    0.9088   -0.1503 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4199    2.4152   -0.0292 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.8952    1.6020    0.1073 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.3666   -0.8435    0.1395 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.4587   -2.4410    0.0315 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.8488   -1.6066   -0.1046 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  2  3  1  0
+  3  4  2  0
+  4  5  1  0
+  5  6  2  0
+  6  7  1  0
+  6  1  1  0
+  1  8  1  0
+  2  9  1  0
+  3 10  1  0
+  4 11  1  0
+  5 12  1  0
+M  CHG  1   7  -1
+M  END

Code/GraphMol/MolInteractionFields/test_data/phenylamine.mol

@@ -0,0 +1,33 @@
+
+     RDKit          3D
+
+ 14 14  0  0  0  0  0  0  0  0999 V2000
+    0.4170    1.1370   -0.1530 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9590    1.3799   -0.1766 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8555    0.3211   -0.0540 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3746   -0.9799    0.0729 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0024   -1.2169    0.0959 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.9105   -0.1566    0.0352 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.2846   -0.4086   -0.0596 N   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1021    1.9702   -0.2828 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.3281    2.3937   -0.3053 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.9258    0.5070   -0.0783 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.0689   -1.8130    0.1395 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.3606   -2.2405    0.1624 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.8627    0.3792    0.2145 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.5718   -1.2725    0.3891 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  2  3  1  0
+  3  4  2  0
+  4  5  1  0
+  5  6  2  0
+  6  7  1  0
+  6  1  1  0
+  1  8  1  0
+  2  9  1  0
+  3 10  1  0
+  4 11  1  0
+  5 12  1  0
+  7 13  1  0
+  7 14  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/phosphatediester.mol

@@ -0,0 +1,37 @@
+
+     RDKit          3D
+
+ 16 15  0  0  0  0  0  0  0  0999 V2000
+    3.1778    0.2302   -0.5824 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.8230   -0.0792   -0.8192 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8209    0.2782    0.4082 P   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7714    1.7772    0.5666 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1447   -0.6257    1.5712 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.5741   -0.2178   -0.2597 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.7148   -0.0527    0.5590 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.9294   -0.5561   -0.1937 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.7600   -0.0481   -1.4646 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.3000    1.3026   -0.4072 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.5541   -0.3324    0.2763 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.5892   -0.6207    1.4863 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.8420    1.0061    0.8062 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.8086   -1.6116   -0.4592 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.0582   -0.0047   -1.1310 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.8355   -0.4452    0.4084 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  3  4  2  0
+  3  5  1  0
+  3  6  1  0
+  6  7  1  0
+  7  8  1  0
+  1  9  1  0
+  1 10  1  0
+  1 11  1  0
+  7 12  1  0
+  7 13  1  0
+  8 14  1  0
+  8 15  1  0
+  8 16  1  0
+M  CHG  1   5  -1
+M  END

Code/GraphMol/MolInteractionFields/test_data/phosphonate.mol

@@ -0,0 +1,21 @@
+
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -0.9071    0.0653    0.0728 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.9339   -0.0672   -0.0750 P   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1001   -1.2018   -1.0783 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3539   -0.3874    1.3543 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3214    1.3176   -0.5790 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1576    0.8611    0.7789 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.3335    0.2968   -0.9066 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.3110   -0.8844    0.4329 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  1  0
+  2  3  2  0
+  2  4  1  0
+  2  5  1  0
+  1  6  1  0
+  1  7  1  0
+  1  8  1  0
+M  CHG  2   4  -1   5  -1
+M  END

Code/GraphMol/MolInteractionFields/test_data/serine.mol

@@ -0,0 +1,32 @@
+
+     RDKit          3D
+
+ 14 13  0  0  0  0  0  0  0  0999 V2000
+   -1.2618    0.3249    0.0306 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1892    0.0066    0.4133 C   0  0  1  0  0  0  0  0  0  0  0  0
+    0.7103   -1.1225   -0.4851 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0622   -2.0998   -0.8279 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0018   -0.9106   -0.8490 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1034    1.1813    0.2661 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.1525   -0.6759    0.5149 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3931    0.4136   -1.0536 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.5744    1.2704    0.4862 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.2544   -0.3589    1.4443 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.2606   -0.0451   -0.4425 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.8371    1.6943   -0.5818 H   0  0  0  0  0  1  0  0  0  0  0  0
+    0.8864    1.8335    1.0270 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.9236   -1.5118    0.0576 H   0  0  0  0  0  1  0  0  0  0  0  0
+  2  1  1  6
+  2  3  1  0
+  3  4  2  0
+  3  5  1  0
+  2  6  1  0
+  1  7  1  0
+  1  8  1  0
+  1  9  1  0
+  2 10  1  0
+  5 11  1  0
+  6 12  1  0
+  6 13  1  0
+  7 14  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/sulfanilamide.mol

@@ -0,0 +1,43 @@
+
+     RDKit          3D
+
+ 19 19  0  0  0  0  0  0  0  0999 V2000
+   -2.6013    0.8376    1.6879 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.3803   -0.0811    0.5954 S   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.8652   -1.4415    0.6185 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6647   -0.0806    0.1709 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1369    1.0048    0.5282 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4973    0.9891    0.2085 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0730   -0.1213   -0.4179 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.4060   -0.0774   -0.8404 N   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2462   -1.1797   -0.8091 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1158   -1.1772   -0.4956 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.0744    0.6332   -0.7556 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.2863    1.8661    1.0396 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.1081    1.8533    0.4575 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.9975    0.5260   -0.2774 H   0  0  0  0  0  1  0  0  0  0  0  0
+    3.8217   -0.9916   -0.9895 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.6595   -2.0207   -1.3603 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.7365   -2.0219   -0.7848 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.5217    1.5114   -0.4924 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -3.7001   -0.0287   -1.2151 H   0  0  0  0  0  1  0  0  0  0  0  0
+  1  2  2  0
+  2  3  2  0
+  2  4  1  0
+  4  5  2  0
+  5  6  1  0
+  6  7  2  0
+  7  8  1  0
+  7  9  1  0
+  9 10  2  0
+  2 11  1  0
+ 10  4  1  0
+  5 12  1  0
+  6 13  1  0
+  8 14  1  0
+  8 15  1  0
+  9 16  1  0
+ 10 17  1  0
+ 11 18  1  0
+ 11 19  1  0
+M  END

Code/GraphMol/MolInteractionFields/test_data/threonine.mol

@@ -0,0 +1,38 @@
+
+     RDKit          3D
+
+ 17 16  0  0  0  0  0  0  0  0999 V2000
+   -1.7910   -0.4061    0.3630 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7464   -0.4347   -0.7493 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.7131   -0.2400   -0.2787 C   0  0  1  0  0  0  0  0  0  0  0  0
+    0.9090    1.1255    0.3918 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2291    2.1267    0.2202 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0055    1.1764    1.1847 O   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1527   -1.3286    0.6270 N   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0931    0.5317   -1.7469 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8332    0.5706    0.8552 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.5991   -1.1724    1.1195 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -2.7884   -0.5803   -0.0561 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -0.8148   -1.4099   -1.2471 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.3719   -0.2276   -1.1564 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.0036    2.0922    1.5338 H   0  0  0  0  0  1  0  0  0  0  0  0
+    2.0630   -1.0596    1.0108 H   0  0  0  0  0  1  0  0  0  0  0  0
+    1.3091   -2.1730    0.0786 H   0  0  0  0  0  1  0  0  0  0  0  0
+   -1.0911    1.4091   -1.3123 H   0  0  0  0  0  1  0  0  0  0  0  0
+  2  1  1  1
+  2  3  1  0
+  3  4  1  1
+  4  5  2  0
+  4  6  1  0
+  3  7  1  0
+  2  8  1  0
+  1  9  1  0
+  1 10  1  0
+  1 11  1  0
+  2 12  1  0
+  3 13  1  0
+  6 14  1  0
+  7 15  1  0
+  7 16  1  0
+  8 17  1  0
+M  END