Create getExperimentalTorsions() function (#5969)

* a bit of minor refactoring

* create API function to get experimental torsion bonds

* put the torsion index in the dictionary
in this case it's silly to return tuples to python

* allow passing an EmbedParameters object from python

* capture the torsion index too

* add an atomIndices component
we should typedef the torsionBonds type

* test torsion index

* changers in response to review

* changes in response to review

* fix merge mistake

Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.cpp

@@ -1,5 +1,5 @@
 //
-//  Copyright (C) 2017 Sereina Riniker
+//  Copyright (C) 2017-2023 Sereina Riniker and other RDKit contributors
 //
 //   @@ All Rights Reserved @@
 //  This file is part of the RDKit.
@@ -50,15 +50,6 @@ const unsigned int MIN_MACROCYCLE_SIZE = 9;
 #include "torsionPreferences_smallrings.in"
 #include "torsionPreferences_macrocycles.in"
 
-//! A structure used to the experimental torsion patterns
-struct ExpTorsionAngle {
-  std::string smarts;
-  std::vector<double> V;
-  std::vector<int> signs;
-  boost::shared_ptr<const ROMol> dp_pattern;
-  unsigned int idx[4];
-};
-
 // class to store the experimental torsion angles
 class ExpTorsionAngleCollection {
  public:
@@ -83,7 +74,7 @@ const ExpTorsionAngleCollection *ExpTorsionAngleCollection::getParams(
     unsigned int version, bool useSmallRingTorsions, bool useMacrocycleTorsions,
     const std::string &paramData) {
   std::string params;
-  if (paramData == "") {
+  if (paramData.empty()) {
     switch (version) {
       case 1:
         params = torsionPreferencesV1;
@@ -106,8 +97,7 @@ const ExpTorsionAngleCollection *ExpTorsionAngleCollection::getParams(
     params += torsionPreferencesMacrocycles;
   }
 
-  const ExpTorsionAngleCollection *res = &(param_flyweight(params).get());
-  return res;
+  return &(param_flyweight(params).get());
 }
 
 ExpTorsionAngleCollection::ExpTorsionAngleCollection(
@@ -116,12 +106,14 @@ ExpTorsionAngleCollection::ExpTorsionAngleCollection(
   std::istringstream inStream(paramData);
 
   std::string inLine = RDKit::getLine(inStream);
+  unsigned int torsionIdx=0;
   while (!inStream.eof()) {
     if (inLine[0] != '#') {
       ExpTorsionAngle angle;
       tokenizer tokens(inLine, tabSep);
       tokenizer::iterator token = tokens.begin();
       angle.smarts = *token;
+      angle.torsionIdx = torsionIdx++;
       ++token;
       for (unsigned int i = 0; i < 12; i += 2) {
         angle.signs.push_back(boost::lexical_cast<int>(*token));
@@ -129,8 +121,7 @@ ExpTorsionAngleCollection::ExpTorsionAngleCollection(
         angle.V.push_back(boost::lexical_cast<double>(*token));
         ++token;
       }
-      angle.dp_pattern =
-          boost::shared_ptr<const ROMol>(SmartsToMol(angle.smarts));
+      angle.dp_pattern.reset(SmartsToMol(angle.smarts));
       // get the atom indices for atom 1, 2, 3, 4 in the pattern
       for (unsigned int i = 0; i < (angle.dp_pattern.get())->getNumAtoms();
            ++i) {
@@ -142,7 +133,7 @@ ExpTorsionAngleCollection::ExpTorsionAngleCollection(
           }
         }
       }
-      d_params.push_back(angle);
+      d_params.push_back(std::move(angle));
     }
     inLine = RDKit::getLine(inStream);
   }  // while loop
@@ -150,10 +141,12 @@ ExpTorsionAngleCollection::ExpTorsionAngleCollection(
   //    << d_params[d_params.size()-1].smarts << std::endl;
 }
 
-void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
-                             bool useExpTorsions, bool useSmallRingTorsions,
-                             bool useMacrocycleTorsions, bool useBasicKnowledge,
-                             unsigned int version, bool verbose) {
+void getExperimentalTorsions(
+    const RDKit::ROMol &mol, CrystalFFDetails &details,
+    std::vector<std::tuple<unsigned int, std::vector<unsigned int>, const ExpTorsionAngle *>> &torsionBonds,
+    bool useExpTorsions, bool useSmallRingTorsions, bool useMacrocycleTorsions,
+    bool useBasicKnowledge, unsigned int version, bool verbose) {
+  torsionBonds.clear();
   unsigned int nb = mol.getNumBonds();
   unsigned int na = mol.getNumAtoms();
   if (!na) {
@@ -174,13 +167,12 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
   boost::dynamic_bitset<> excludedBonds(nb);
   const RingInfo *rinfo = mol.getRingInfo();
   const VECT_INT_VECT &bondRings = rinfo->bondRings();
-  VECT_INT_VECT_CI rii, rjj;
-  for (rii = bondRings.begin(); rii != bondRings.end(); ++rii) {
+  for (auto rii = bondRings.begin(); rii != bondRings.end(); ++rii) {
     boost::dynamic_bitset<> rs1(nb);  // bitset for ring 1
     for (auto riiv : *rii) {
       rs1[riiv] = 1;
     }
-    for (rjj = rii + 1; rjj != bondRings.end(); ++rjj) {
+    for (auto rjj = rii + 1; rjj != bondRings.end(); ++rjj) {
       // we don't worry about the overlap if both rings are macrocycles:
       if (rii->size() >= MIN_MACROCYCLE_SIZE &&
           rjj->size() >= MIN_MACROCYCLE_SIZE) {
@@ -189,8 +181,7 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
       unsigned int nInCommon = 0;
       for (auto rjj_i : *rjj) {
         if (rs1[rjj_i]) {
-          ++nInCommon;
-          if (nInCommon > 1) {
+          if (++nInCommon > 1) {
             break;
           }
         }
@@ -215,29 +206,31 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
 
   if (useExpTorsions) {
     // we set the torsion angles with experimental data
-    const ExpTorsionAngleCollection *params =
-        ExpTorsionAngleCollection::getParams(version, useSmallRingTorsions,
-                                             useMacrocycleTorsions);
-
+    const auto *params = ExpTorsionAngleCollection::getParams(
+        version, useSmallRingTorsions, useMacrocycleTorsions);
+    CHECK_INVARIANT(params, "no parameters available");
     // loop over patterns
     for (const auto &param : *params) {
       std::vector<MatchVectType> matches;
       SubstructMatch(mol, *(param.dp_pattern.get()), matches, false, true);
       // loop over matches
-      for (std::vector<MatchVectType>::const_iterator matchIt = matches.begin();
-           matchIt != matches.end(); ++matchIt) {
+      for (const auto &match : matches) {
         // get bond indices
-        aid1 = (*matchIt)[param.idx[0]].second;
-        aid2 = (*matchIt)[param.idx[1]].second;
-        aid3 = (*matchIt)[param.idx[2]].second;
-        aid4 = (*matchIt)[param.idx[3]].second;
-        // FIX: check if bond is NULL
-        bid2 = mol.getBondBetweenAtoms(aid2, aid3)->getIdx();
+        aid1 = match[param.idx[0]].second;
+        aid2 = match[param.idx[1]].second;
+        aid3 = match[param.idx[2]].second;
+        aid4 = match[param.idx[3]].second;
+        const auto bnd = mol.getBondBetweenAtoms(aid2, aid3);
+        CHECK_INVARIANT(bnd, "bond between central atoms not found")
+        bid2 = bnd->getIdx();
+
         // check that a bond is part of maximum one ring
         if (excludedBonds[bid2] || mol.getRingInfo()->numBondRings(bid2) > 3) {
           doneBonds[bid2] = 1;
         }
         if (!doneBonds[bid2]) {
+          std::vector<unsigned int> aids{aid1,aid2,aid3,aid4};
+          torsionBonds.emplace_back(bid2, aids, &param);
           doneBonds[bid2] = 1;
           std::vector<int> atoms(4);
           atoms[0] = aid1;
@@ -251,11 +244,11 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
             // extra formatting provided by the logger after every entry;
             std::stringstream sstr;
             sstr << param.smarts << ": " << aid1 << " " << aid2 << " " << aid3
-                 << " " << aid4 << ", (";
+                 << " " << aid4 << ", [";
             for (unsigned int i = 0; i < param.V.size() - 1; ++i) {
-              sstr << param.V[i] << ", ";
+              sstr << "(" << param.signs[i] << " " << param.V[i] << "), ";
             }
-            sstr << param.V[param.V.size() - 1] << ") ";
+            sstr << "(" << param.signs.back() << " " << param.V.back() << ")] ";
             BOOST_LOG(rdInfoLog) << sstr.str() << std::endl;
           }
         }  // if not donePaths
@@ -268,8 +261,6 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
   // straight triple bonds, etc.
   if (useBasicKnowledge) {
     boost::dynamic_bitset<> doneAtoms(na);
-    ROMol::ADJ_ITER nbrIdx;
-    ROMol::ADJ_ITER endNbrs;
 
     // inversion terms (improper torsions / out-of-plane bends / inversion)
     // loop over atoms
@@ -280,23 +271,17 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
         const Atom *atom2 = mol.getAtomWithIdx(atoms[1]);
         int at2AtomicNum = atom2->getAtomicNum();
 
-        // if atom is a N,O or C and SP2-hybridized
+        // if atom is a N,O or C, SP2-hybridized, and has three neighbors
         if (((at2AtomicNum == 6) || (at2AtomicNum == 7) ||
              (at2AtomicNum == 8)) &&
-            (atom2->getHybridization() == Atom::SP2)) {
-          // get neighbors
-          boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(atom2);
-          // check if enough neighbours
-          if (mol.getAtomDegree(atom2) != 3) {
-            continue;
-          }
+            (atom2->getHybridization() == Atom::SP2) &&
+            mol.getAtomDegree(atom2) == 3) {
           unsigned int i = 0;
           unsigned int isBoundToSP2O = 0;  // false
-          for (; nbrIdx != endNbrs; ++nbrIdx) {
-            const Atom *atomX = mol[*nbrIdx];
+          for (const auto atomX : mol.atomNeighbors(atom2)) {
             atoms[i] = atomX->getIdx();
-            // if the central atom is sp2 carbon and is bound to sp2 oxygen, set
-            // a flag
+            // if the central atom is sp2 carbon and is bound to sp2 oxygen,
+            // set a flag
             if (!isBoundToSP2O) {
               isBoundToSP2O =
                   ((at2AtomicNum == 6) && (atomX->getAtomicNum() == 8) &&
@@ -311,8 +296,8 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
           atoms.push_back(isBoundToSP2O);
           details.improperAtoms.push_back(atoms);
           /*if (verbose) {
-            std::cout << "out-of-plane bend: " << atoms[0] << " " << atoms[1] <<
-          " "
+            std::cout << "out-of-plane bend: " << atoms[0] << " " << atoms[1]
+          << " "
                 << atoms[2] << " " << atoms[3] << std::endl;
           }*/
         }
@@ -320,11 +305,10 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
     }
 
     // torsions for flat rings
-    const RingInfo *rinfo =
-        mol.getRingInfo();  // FIX: make sure we have ring info
-    CHECK_INVARIANT(rinfo, "");
-    const VECT_INT_VECT &atomRings = rinfo->atomRings();
-    for (const auto &atomRing : atomRings) {
+    const RingInfo *rinfo = mol.getRingInfo();
+    CHECK_INVARIANT(rinfo, "no ring info");
+    CHECK_INVARIANT(rinfo->isInitialized(), "ring info not initialized");
+    for (const auto &atomRing : rinfo->atomRings()) {
       unsigned int rSize = atomRing.size();
       // we don't need to deal with 3 membered rings
       // and we do not treat rings greater than 6
@@ -352,14 +336,15 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
           atoms[2] = aid3;
           atoms[3] = aid4;
           details.expTorsionAtoms.push_back(atoms);
+
           std::vector<int> signs(6, 1);
           signs[1] = -1;  // MMFF sign for m = 2
           std::vector<double> fconsts(6, 0.0);
           fconsts[1] = 100.0;  // 7.0 is MMFF force constants for aromatic rings
           details.expTorsionAngles.emplace_back(signs, fconsts);
           /*if (verbose) {
-            std::cout << "SP2 ring: " << aid1 << " " << aid2 << " " << aid3 << "
-          " << aid4 << std::endl;
+            std::cout << "SP2 ring: " << aid1 << " " << aid2 << " " << aid3 <<
+          " " << aid4 << std::endl;
           }*/
         }
 
@@ -369,5 +354,15 @@ void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
 
 }  // end function
 
+void getExperimentalTorsions(const RDKit::ROMol &mol, CrystalFFDetails &details,
+                             bool useExpTorsions, bool useSmallRingTorsions,
+                             bool useMacrocycleTorsions, bool useBasicKnowledge,
+                             unsigned int version, bool verbose) {
+  std::vector<std::tuple<unsigned int, std::vector<unsigned int>, const ExpTorsionAngle *>> torsionBonds;
+  getExperimentalTorsions(mol, details, torsionBonds, useExpTorsions,
+                          useSmallRingTorsions, useMacrocycleTorsions,
+                          useBasicKnowledge, version, verbose);
+}
+
 }  // namespace CrystalFF
 }  // namespace ForceFields

Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.h

@@ -1,5 +1,5 @@
 //
-//  Copyright (C) 2017 Sereina Riniker
+//  Copyright (C) 2017-2023 Sereina Riniker and other RDKit contributors
 //
 //   @@ All Rights Reserved @@
 //  This file is part of the RDKit.
@@ -11,6 +11,8 @@
 #ifndef _RD_TORSIONPREFERENCES_H_
 #define _RD_TORSIONPREFERENCES_H_
 #include <vector>
+#include <string>
+#include <memory>
 
 namespace RDKit {
 class ROMol;
@@ -18,6 +20,17 @@ class ROMol;
 
 namespace ForceFields {
 namespace CrystalFF {
+
+//! A structure used to the experimental torsion patterns
+struct RDKIT_FORCEFIELDHELPERS_EXPORT ExpTorsionAngle {
+  unsigned int torsionIdx;
+  std::string smarts;
+  std::vector<double> V;
+  std::vector<int> signs;
+  std::unique_ptr<const RDKit::ROMol> dp_pattern;
+  unsigned int idx[4];
+};
+
 struct CrystalFFDetails {
   std::vector<std::vector<int>> expTorsionAtoms;
   std::vector<std::pair<std::vector<int>, std::vector<double>>>
@@ -35,6 +48,15 @@ RDKIT_FORCEFIELDHELPERS_EXPORT void getExperimentalTorsions(
     bool useExpTorsions = false, bool useSmallRingTorsions = false,
     bool useMacrocycleTorsions = false, bool useBasicKnowledge = false,
     unsigned int version = 1, bool verbose = false);
+
+//! \overload
+RDKIT_FORCEFIELDHELPERS_EXPORT void getExperimentalTorsions(
+    const RDKit::ROMol &mol, CrystalFFDetails &details,
+    std::vector<std::tuple<unsigned int, std::vector<unsigned int>, const ExpTorsionAngle *>> &torsionBonds,
+    bool useExpTorsions = false, bool useSmallRingTorsions = false,
+    bool useMacrocycleTorsions = false, bool useBasicKnowledge = false,
+    unsigned int version = 1, bool verbose = false);
+
 }  // namespace CrystalFF
 }  // namespace ForceFields
 

Code/GraphMol/ForceFieldHelpers/CrystalFF/testCrystalFF.cpp

@@ -132,14 +132,28 @@ void testTorsionPrefs() {
   TEST_ASSERT(details.expTorsionAngles[0].first.size() == 6);
   TEST_ASSERT(details.expTorsionAngles[0].second.size() == 6);
 
+  std::vector<
+      std::tuple<unsigned int, std::vector<unsigned int>, const ForceFields::CrystalFF::ExpTorsionAngle *>>
+      torsionBonds;
+  ForceFields::CrystalFF::getExperimentalTorsions(
+      *mol, details, torsionBonds, true, false, false, false, 2, false);
+  TEST_ASSERT(torsionBonds.size() == 1);
+  TEST_ASSERT(std::get<0>(torsionBonds[0]) == 1);
+  TEST_ASSERT(std::get<2>(torsionBonds[0])->smarts ==
+              "[!#1:1][CX4H2:2]!@;-[CX4H2:3][!#1:4]");
+  TEST_ASSERT(std::get<2>(torsionBonds[0])->torsionIdx == 229);
+
   delete mol;
   mol = SmilesToMol("CCCCC");
   TEST_ASSERT(mol);
 
-  ForceFields::CrystalFF::getExperimentalTorsions(*mol, details, true, false,
-                                                  false, false, 1, false);
+  ForceFields::CrystalFF::getExperimentalTorsions(
+      *mol, details, torsionBonds, true, false, false, false, 2, false);
   TEST_ASSERT(details.expTorsionAtoms.size() == 2);
   TEST_ASSERT(details.expTorsionAngles.size() == 2);
+  TEST_ASSERT(torsionBonds.size() == 2);
+  TEST_ASSERT(std::get<0>(torsionBonds[0]) == 1);
+  TEST_ASSERT(std::get<0>(torsionBonds[1]) == 2);
   delete mol;
 }