Cookbook entry

.gitignore

@@ -50,6 +50,8 @@ __pycache__/
 /Code/GraphMol/SmilesParse/smiles.tab.cpp
 /Code/GraphMol/SmilesParse/smiles.tab.hpp
 
+/Data/ANIParams
+
 /External/**/*.tar.gz
 /External/**/*.zip
 /External/catch/catch/

CMakeLists.txt

@@ -23,6 +23,7 @@ add_dependencies(rdkit_base catch)
 option(RDK_BUILD_SWIG_WRAPPERS "build the SWIG wrappers" OFF )
 option(RDK_BUILD_PYTHON_WRAPPERS "build the standard python wrappers" ON )
 option(RDK_BUILD_COMPRESSED_SUPPLIERS "build in support for compressed MolSuppliers" OFF )
+option(RDK_BUILD_ANI "build support for ANI style ForceField" OFF)
 option(RDK_BUILD_INCHI_SUPPORT "build the rdkit inchi wrapper" OFF )
 option(RDK_BUILD_AVALON_SUPPORT "install support for the avalon toolkit. Use the variable AVALONTOOLS_DIR to set the location of the source." OFF )
 option(RDK_BUILD_PGSQL "build the PostgreSQL cartridge" OFF )
@@ -380,6 +381,35 @@ if(RDK_BUILD_DESCRIPTORS3D)
   endif(EIGEN3_FOUND)
 endif(RDK_BUILD_DESCRIPTORS3D)
 
+if(RDK_BUILD_ANI)
+  find_package(Eigen3)
+  if(EIGEN3_FOUND)
+    set(RDK_HAS_EIGEN ON)
+    if(NOT WIN32)
+      target_compile_definitions(rdkit_base INTERFACE "-DRDK_HAS_EIGEN3" "-DRDK_BUILD_ANI")
+      target_link_libraries(rdkit_base INTERFACE Eigen3::Eigen)
+      set(ANIParamsDir "${CMAKE_CURRENT_SOURCE_DIR}/Data/ANIParams")
+      if(NOT EXISTS "${ANIParamsDir}/ANI-1x/selfEnergies")
+        set(RELEASE_NO "1.0")
+        downloadAndCheckMD5("https://github.com/rdkit/rdkit_ANI_weights/archive/v${RELEASE_NO}.tar.gz"
+                "${CMAKE_CURRENT_SOURCE_DIR}/ANIWeights-v${RELEASE_NO}.tar.gz" "")
+        execute_process(COMMAND ${CMAKE_COMMAND} -E tar zxf
+          ${CMAKE_CURRENT_SOURCE_DIR}/ANIWeights-v${RELEASE_NO}.tar.gz
+          WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+        file(RENAME "${CMAKE_CURRENT_SOURCE_DIR}/rdkit_ANI_weights-${RELEASE_NO}/Params" "${CMAKE_CURRENT_SOURCE_DIR}/Data/ANIParams")
+        file(REMOVE_RECURSE "${CMAKE_CURRENT_SOURCE_DIR}/rdkit_ANI_weights-${RELEASE_NO}")
+        file(REMOVE "${CMAKE_CURRENT_SOURCE_DIR}/ANIWeights-v${RELEASE_NO}.tar.gz")
+      endif()
+    endif()
+  else(EIGEN3_FOUND)
+    message("Eigen3 not found, disabling the ANI build.")
+    set(RDK_BUILD_ANI OFF)
+  endif(EIGEN3_FOUND)
+  if(WIN32)
+    message("ANI not supported on Windows")
+    set(RDK_BUILD_ANI OFF)
+  endif(WIN32)
+endif(RDK_BUILD_ANI)
 
 find_package (Threads)
 set(RDKit_THREAD_LIBS Threads::Threads)

Code/ForceField/ANI/AtomicContrib.cpp

@@ -61,7 +61,7 @@ ANIAtomContrib::ANIAtomContrib(ForceField *owner, int atomType,
   // Different values for means of the gaussian symmetry functions
   std::string path = getenv("RDBASE");
   std::string paramFilePath =
-      path + "/Code/ForceField/ANI/Params/" + modelType + "/AEVParams/";
+      path + "/Data/ANIParams/" + modelType + "/AEVParams/";
 
   // Weights for the radial symmetry functions
   ArrayXd ShfR;
@@ -121,8 +121,10 @@ double ANIAtomContrib::forwardProp(ArrayXXd &aev) const {
 }
 
 double ANIAtomContrib::getEnergy(double *pos) const {
-  auto aev = RDKit::Descriptors::ANI::AtomicEnvironmentVector(
-      pos, this->d_speciesVec, this->d_numAtoms, &(this->d_aevParams));
+  PRECONDITION(pos != nullptr, "Positions array is NULL");
+  ArrayXXd aev;
+  RDKit::Descriptors::ANI::AtomicEnvironmentVector(
+      aev, pos, this->d_speciesVec, this->d_numAtoms, &(this->d_aevParams));
   ArrayXXd row = aev.row(this->d_atomIdx);
   return this->ANIAtomContrib::forwardProp(row) + this->d_selfEnergy;
 }
@@ -132,7 +134,43 @@ double ANIAtomContrib::getEnergy(Eigen::ArrayXXd &aev) const {
   return this->ANIAtomContrib::forwardProp(row) + this->d_selfEnergy;
 }
 
-void ANIAtomContrib::getGrad(double *pos, double *grad) const {}
+void ANIAtomContrib::getGrad(double *pos, double *grad) const {
+  PRECONDITION(pos != nullptr, "Positions are null");
+  PRECONDITION(grad != nullptr, "Gradient array is null");
+  double displacement = 1e-5;
+
+  // + - x movement
+  pos[3 * this->d_atomIdx] += displacement;
+  auto posXEnergy = this->dp_forceField->calcEnergy(pos);
+
+  pos[3 * this->d_atomIdx] -= 2 * displacement;
+  auto negXEnergy = this->dp_forceField->calcEnergy(pos);
+
+  grad[3 * this->d_atomIdx] = (posXEnergy - negXEnergy) / (2 * displacement);
+  pos[3 * this->d_atomIdx] += displacement;
+
+  // + - Y movement
+  pos[3 * this->d_atomIdx + 1] += displacement;
+  auto posYEnergy = this->dp_forceField->calcEnergy(pos);
+
+  pos[3 * this->d_atomIdx + 1] -= 2 * displacement;
+  auto negYEnergy = this->dp_forceField->calcEnergy(pos);
+
+  grad[3 * this->d_atomIdx + 1] =
+      (posYEnergy - negYEnergy) / (2 * displacement);
+  pos[3 * this->d_atomIdx + 1] += displacement;
+
+  // + - Z movement
+  pos[3 * this->d_atomIdx + 2] += displacement;
+  auto posZEnergy = this->dp_forceField->calcEnergy(pos);
+
+  pos[3 * this->d_atomIdx + 2] -= 2 * displacement;
+  auto negZEnergy = this->dp_forceField->calcEnergy(pos);
+
+  grad[3 * this->d_atomIdx + 2] =
+      (posZEnergy - negZEnergy) / (2 * displacement);
+  pos[3 * this->d_atomIdx + 2] += displacement;
+}
 
 namespace Utils {
 
@@ -155,8 +193,9 @@ std::vector<std::string> tokenize(const std::string &s) {
 void loadFromBin(std::vector<MatrixXd> *weights, unsigned int model,
                  std::string weightType, unsigned int layer,
                  std::string atomType, std::string modelType) {
+  PRECONDITION(weights != nullptr, "Vector of weights is NULL");
   std::string path = getenv("RDBASE");
-  std::string paramFile = path + "/Code/ForceField/ANI/Params/" + modelType +
+  std::string paramFile = path + "/Data/ANIParams/" + modelType +
                           "/model" + std::to_string(model) + "/" + atomType +
                           "_" + std::to_string(layer) + "_" + weightType +
                           ".bin";
@@ -168,8 +207,10 @@ void loadFromBin(std::vector<MatrixXd> *weights, unsigned int model,
 void loadFromBin(std::vector<MatrixXd> *weights, std::vector<MatrixXd> *biases,
                  unsigned int model, std::string atomType,
                  std::string modelType) {
+  PRECONDITION(weights != nullptr, "Weights array is null");
+  PRECONDITION(biases != nullptr, "Biases array is null");
   std::string path = getenv("RDBASE");
-  std::string paramFile = path + "/Code/ForceField/ANI/Params/" + modelType +
+  std::string paramFile = path + "/Data/ANIParams/" + modelType +
                           "/model" + std::to_string(model) + ".bin";
   std::vector<MatrixXf> floatWeights, floatBiases;
   RDNumeric::EigenSerializer::deserializeAll(&floatWeights, &floatBiases,
@@ -183,8 +224,9 @@ void loadFromBin(std::vector<MatrixXd> *weights, std::vector<MatrixXd> *biases,
 void loadFromCSV(std::vector<MatrixXd> *weights, unsigned int model,
                  std::string weightType, unsigned int layer,
                  std::string atomType, std::string modelType) {
+  PRECONDITION(weights != nullptr, "Weights array is null");
   std::string path = getenv("RDBASE");
-  std::string paramFile = path + "/Code/ForceField/ANI/Params/" + modelType +
+  std::string paramFile = path + "/Data/ANIParams/" + modelType +
                           "/model" + std::to_string(model) + "/" + atomType +
                           "_" + std::to_string(layer) + "_" + weightType;
 
@@ -227,7 +269,7 @@ void loadSelfEnergy(double *energy, std::string atomType,
                     std::string modelType) {
   std::string path = getenv("RDBASE");
   std::string filePath =
-      path + "/Code/ForceField/ANI/Params/" + modelType + "/selfEnergies";
+      path + "/Data/ANIParams/" + modelType + "/selfEnergies";
 
   std::ifstream selfEnergyFile(filePath.c_str());
   if (!selfEnergyFile.good()) {

Code/ForceField/ANI/CMakeLists.txt

@@ -1,7 +1,7 @@
-if(RDK_HAS_EIGEN)
+if(RDK_BUILD_ANI)
 rdkit_catch_test(ANIForceFieldCatchTest ANIForceField_catch.cpp  catch_main.cpp
     LINK_LIBRARIES 
     DistGeomHelpers ForceFieldHelpers
     FileParsers MolTransforms SmilesParse 
     SubstructMatch MolAlign ANIDescriptors ANIForceField)
-endif(RDK_HAS_EIGEN)
+endif(RDK_BUILD_ANI)

Code/ForceField/CMakeLists.txt

@@ -14,11 +14,12 @@ rdkit_library(ForceField
               LINK_LIBRARIES Optimizer Trajectory)
 target_compile_definitions(ForceField PRIVATE RDKIT_FORCEFIELD_BUILD)
 
-if(RDK_HAS_EIGEN)
+if(RDK_BUILD_ANI)
 rdkit_library(ANIForceField
              ANI/AtomicContrib.cpp
-             LINK_LIBRARIES Optimizer Trajectory ANIDescriptors EigenSerializer)
-endif(RDK_HAS_EIGEN)
+             LINK_LIBRARIES Optimizer Trajectory ANIDescriptors EigenSerializer ForceField)
+target_compile_definitions(ANIForceField PRIVATE RDKIT_ANIFORCEFIELD_BUILD)
+endif(RDK_BUILD_ANI)
 
 rdkit_headers(Contrib.h
               ForceField.h DEST ForceField)
@@ -47,15 +48,15 @@ rdkit_headers(MMFF/AngleBend.h
               MMFF/TorsionConstraint.h
               MMFF/PositionConstraint.h
               MMFF/Params.h DEST ForceField/MMFF)
-if(RDK_HAS_EIGEN)
+if(RDK_BUILD_ANI)
 rdkit_headers(ANI/AtomicContrib.h DEST ForceField/ANI)
-endif(RDK_HAS_EIGEN)
+endif(RDK_BUILD_ANI)
 
 add_subdirectory(UFF)
 add_subdirectory(MMFF)
-if(RDK_HAS_EIGEN)
+if(RDK_BUILD_ANI)
 add_subdirectory(ANI)
-endif(RDK_HAS_EIGEN)
+endif(RDK_BUILD_ANI)
 if(RDK_BUILD_PYTHON_WRAPPERS)
 add_subdirectory(Wrap)
 endif()

Code/ForceField/Wrap/CMakeLists.txt

@@ -1,6 +1,13 @@
 remove_definitions(-DRDKIT_FORCEFIELD_BUILD)
+if(RDK_BUILD_ANI)
 rdkit_python_extension(rdForceField ForceField.cpp
                        DEST ForceField
-                       LINK_LIBRARIES ForceFieldHelpers ANIForceField)
+                       LINK_LIBRARIES ForceFieldHelpers ANIForceField ANIForceFieldHelpers)
+add_pytest(pyForceFieldANI ${CMAKE_CURRENT_SOURCE_DIR}/testANI.py)
+else()
+rdkit_python_extension(rdForceField ForceField.cpp
+                       DEST ForceField
+                       LINK_LIBRARIES ForceFieldHelpers)
+endif(RDK_BUILD_ANI)
 add_pytest(pyForceFieldConstraints ${CMAKE_CURRENT_SOURCE_DIR}/testConstraints.py)
 add_pytest(pyForceFieldANI ${CMAKE_CURRENT_SOURCE_DIR}/testANI.py)

Code/ForceField/Wrap/ForceField.cpp

@@ -22,13 +22,15 @@
 #include <ForceField/MMFF/AngleConstraint.h>
 #include <ForceField/MMFF/TorsionConstraint.h>
 #include <ForceField/MMFF/PositionConstraint.h>
+#ifdef RDK_BUILD_ANI
 #include <ForceField/ANI/AtomicContrib.h>
+#endif
 #include "PyForceField.h"
 
 using namespace ForceFields;
 namespace python = boost::python;
 
-#ifdef RDK_HAS_EIGEN3
+#ifdef RDK_BUILD_ANI
 void ANIAddAtomContrib(PyForceField *self, python::list speciesVec,
                        int atomType, unsigned int atomIdx,
                        unsigned int numAtoms, unsigned int numLayers,
@@ -428,7 +430,7 @@ BOOST_PYTHON_MODULE(rdForceField) {
            (python::arg("self"), python::arg("idx"), python::arg("maxDispl"),
             python::arg("forceConstant")),
            "Adds a position constraint to the MMFF force field.")
-#ifdef RDK_HAS_EIGEN3
+#ifdef RDK_BUILD_ANI
       .def("AddANIAtomContrib", ANIAddAtomContrib,
            (python::arg("self"), python::arg("speciesVec"),
             python::arg("atomType"), python::arg("atomIdx"),