Fixes #7181 (#7206)

* Fixes #7181

adds a new constrainedAtoms member to the etkdgdetails structure which is used to propagate info about the coordMap into the ETK minimizations

* Update Code/DistGeom/DistGeomUtils.cpp

Co-authored-by: Paolo Tosco <paolo.tosco.mail@gmail.com>

* Update Code/DistGeom/DistGeomUtils.cpp

Co-authored-by: Paolo Tosco <paolo.tosco.mail@gmail.com>

---------

Co-authored-by: Paolo Tosco <paolo.tosco.mail@gmail.com>

Code/DistGeom/DistGeomUtils.cpp

@@ -304,7 +304,8 @@ ForceFields::ForceField *construct3DForceField(
     }
   }
 
-  double fdist = 100.0;  // force constant
+  constexpr double knownDistanceConstraintForce = 100.0;
+  double fdist = knownDistanceConstraintForce;  // force constant
   // 1,2 distance constraints
   for (const auto &bnd : etkdgDetails.bonds) {
     unsigned int i = bnd.first;
@@ -347,13 +348,23 @@ ForceFields::ForceField *construct3DForceField(
     }
   }
 
-  // minimum distance for all other atom pairs
-  fdist = etkdgDetails.boundsMatForceScaling * 10.0;
+  // minimum distance for all other atom pairs that aren't constrained
   for (unsigned int i = 1; i < N; ++i) {
     for (unsigned int j = 0; j < i; ++j) {
       if (!atomPairs[j * N + i]) {
+        fdist = etkdgDetails.boundsMatForceScaling * 10.0;
         double l = mmat.getLowerBound(i, j);
         double u = mmat.getUpperBound(i, j);
+        if (!etkdgDetails.constrainedAtoms.empty() &&
+            etkdgDetails.constrainedAtoms[i] &&
+            etkdgDetails.constrainedAtoms[j]) {
+          // we're constrained, so use very tight bounds
+          l = u = ((*positions[i]) - (*positions[j])).length();
+          constexpr double INCR = 0.01;
+          l -= INCR;
+          u += INCR;
+          fdist = knownDistanceConstraintForce;
+        }
         auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
             field, i, j, l, u, fdist);
         field->contribs().push_back(ForceFields::ContribPtr(contrib));

Code/GraphMol/DistGeomHelpers/Embedder.cpp

@@ -1447,6 +1447,12 @@ void EmbedMultipleConfs(ROMol &mol, INT_VECT &res, unsigned int numConfs,
         << std::endl;
     coordMap = nullptr;
   }
+  boost::dynamic_bitset<> constrainedAtoms(mol.getNumAtoms());
+  if (coordMap) {
+    for (const auto &entry : *coordMap) {
+      constrainedAtoms.set(entry.first);
+    }
+  }
 
   if (molFrags.size() > 1 && params.boundsMat != nullptr) {
     BOOST_LOG(rdWarningLog)
@@ -1465,6 +1471,7 @@ void EmbedMultipleConfs(ROMol &mol, INT_VECT &res, unsigned int numConfs,
     unsigned int nAtoms = piece->getNumAtoms();
 
     ForceFields::CrystalFF::CrystalFFDetails etkdgDetails;
+    etkdgDetails.constrainedAtoms = constrainedAtoms;
     EmbeddingOps::initETKDG(piece.get(), params, etkdgDetails);
 
     DistGeom::BoundsMatPtr mmat;

Code/GraphMol/DistGeomHelpers/catch_tests.cpp

@@ -21,6 +21,7 @@
 #include <GraphMol/FileParsers/MolSupplier.h>
 #include <GraphMol/SmilesParse/SmilesParse.h>
 #include <GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.h>
+#include <GraphMol/MolAlign/AlignMolecules.h>
 #include "Embedder.h"
 #include "BoundsMatrixBuilder.h"
 #include <tuple>
@@ -930,4 +931,35 @@ TEST_CASE(
       CHECK(m->getAtomWithIdx(8)->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW);
     }
   }
+}
+
+TEST_CASE("Github #7181: ET terms applied to constrained atoms") {
+  SECTION("basics") {
+    auto templ =
+        "CNc1ccc(OC)cc1 |(-3.3363,0.129414,1.28582;-2.44714,-0.687978,0.507453;-1.11383,-0.29452,0.197587;-0.622766,0.911164,0.645083;0.652332,1.29026,0.350281;1.45603,0.462513,-0.400278;2.7718,0.891528,-0.684446;3.83908,0.0736652,-0.224516;0.984393,-0.734112,-0.850528;-0.300532,-1.12218,-0.556656)|"_smiles;
+    REQUIRE(templ);
+    auto mol = "COc1ccc(NC(C)C)cc1"_smiles;
+    REQUIRE(mol);
+    MolOps::addHs(*mol);
+    auto matches = SubstructMatch(*mol, *templ);
+    REQUIRE(matches.size() == 1);
+
+    auto tconf = templ->getConformer();
+    std::map<int, RDGeom::Point3D> cmap;
+    for (auto [ti, mi] : matches[0]) {
+      cmap[mi] = tconf.getAtomPos(ti);
+    }
+
+    DGeomHelpers::EmbedParameters ps = DGeomHelpers::ETKDGv3;
+    ps.randomSeed = 0xC0FFEE;
+    ps.coordMap = &cmap;
+    auto cid = DGeomHelpers::EmbedMolecule(*mol, ps);
+    CHECK(cid >= 0);
+    auto imatch = matches[0];
+    for (auto &[ti, mi] : imatch) {
+      std::swap(ti, mi);
+    }
+    auto rmsd = MolAlign::alignMol(*mol, *templ, cid, -1, &imatch);
+    CHECK(rmsd < 0.2);
+  }
 }

Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.cpp

@@ -230,6 +230,15 @@ void getExperimentalTorsions(
           doneBonds[bid2] = 1;
         }
         if (!doneBonds[bid2]) {
+          // do not add ET terms between constrained atoms
+          // REVIEW: do we really need to check all 4 atoms?
+          if (!details.constrainedAtoms.empty() &&
+              details.constrainedAtoms[aid1] &&
+              details.constrainedAtoms[aid2] &&
+              details.constrainedAtoms[aid3] &&
+              details.constrainedAtoms[aid4]) {
+            continue;
+          }
           std::vector<unsigned int> aids{aid1, aid2, aid3, aid4};
           torsionBonds.emplace_back(bid2, aids, &param);
           doneBonds[bid2] = 1;

Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.h

@@ -13,6 +13,7 @@
 #include <vector>
 #include <string>
 #include <memory>
+#include <boost/dynamic_bitset.hpp>
 
 namespace RDKit {
 class ROMol;
@@ -40,6 +41,7 @@ struct CrystalFFDetails {
   std::vector<std::vector<int>> angles;
   std::vector<int> atomNums;
   double boundsMatForceScaling;
+  boost::dynamic_bitset<> constrainedAtoms;
 };
 
 //! Get the experimental torsional angles in a molecule