Fixes github 7331 (#7352)

* Fixes #7331 Essentially this modifies the code to use the isolobal analogy Also improves error reporting a bit. * generalize that a bit * more testing
rdkit · Apr 19, 2024 · 4f5d641 · 4f5d641
1 parent 848b57f
commit 4f5d641
Show file tree

Hide file tree

Showing 2 changed files with 54 additions and 9 deletions.
diff --git a/Code/GraphMol/DetermineBonds/DetermineBonds.cpp b/Code/GraphMol/DetermineBonds/DetermineBonds.cpp
@@ -75,18 +75,22 @@ std::vector<T> LazyCartesianProduct<T>::entryAt(uint1024_t pos) const {
 std::vector<unsigned int> possibleValences(
     const RDKit::Atom *atom,
     const std::unordered_map<int, std::vector<unsigned int>> &atomicValence) {
-  auto atomNum = atom->getAtomicNum();
+  auto atomNum = atom->getAtomicNum() - atom->getFormalCharge();
   auto numBonds = atom->getDegree();
 
+  std::vector<unsigned int> avalences;
   auto valences = atomicValence.find(atomNum);
-  if (valences == atomicValence.end()) {
-    std::stringstream ss;
-    ss << "determineBondOrdering() does not work with element "
-       << RDKit::PeriodicTable::getTable()->getElementSymbol(atomNum);
-    throw ValueErrorException(ss.str());
+  if (valences != atomicValence.end()) {
+    avalences = valences->second;
+  } else {
+    for (auto v : RDKit::PeriodicTable::getTable()->getValenceList(atomNum)) {
+      if (v >= 0) {
+        avalences.push_back(v);
+      }
+    }
   }
   std::vector<unsigned int> possible;
-  for (const auto &valence : valences->second) {
+  for (const auto &valence : avalences) {
     if (numBonds <= valence) {
       possible.push_back(valence);
     }
@@ -264,7 +268,8 @@ void setAtomicCharges(RWMol &mol, const std::vector<unsigned int> &valency,
   int molCharge = 0;
   for (unsigned int i = 0; i < mol.getNumAtoms(); i++) {
     auto atom = mol.getAtomWithIdx(i);
-    int atomCharge = getAtomicCharge(atom->getAtomicNum(), valency[i]);
+    int atomCharge = getAtomicCharge(
+        atom->getAtomicNum() - atom->getFormalCharge(), valency[i]);
     molCharge += atomCharge;
     if (atom->getAtomicNum() == 6) {
       if (atom->getDegree() == 2 && valency[i] == 2) {
@@ -342,8 +347,11 @@ void addBondOrdering(RWMol &mol,
   }
 
   if (MolOps::getFormalCharge(mol) != charge) {
+    mol.debugMol(std::cerr);
     std::stringstream ss;
-    ss << "Final molecular charge does not match input; could not find valid bond ordering";
+    ss << "Final molecular charge (" << charge << ") does not match input ("
+       << MolOps::getFormalCharge(mol)
+       << "); could not find valid bond ordering";
     throw ValueErrorException(ss.str());
   }
 

diff --git a/Code/GraphMol/DetermineBonds/catch_tests.cpp b/Code/GraphMol/DetermineBonds/catch_tests.cpp
@@ -439,4 +439,41 @@ TEST_CASE(
       }
     }
   }
+}
+
+TEST_CASE(
+    "github #7331: DetermineBondOrders() makes incorrect assumptions about valence") {
+  SECTION("as reported") {
+    // do not anything here that needs implicit Hs
+    std::vector<std::string> smiles = {
+        "O=NO[Cl+][O-]",
+        "[O-][I+3]([O-])([O-])[O-]",
+        "[O-][I+2]([O-])[O-]",
+        "F[P-](F)(F)(F)(F)F",
+        "F[C+](F)F",
+        "F[C-](F)F",
+        "F[N+](F)(F)F",
+        "F[N-]F",
+        "F[Cl+]F",
+        "F[Br+]F",
+        "O=[Cl+]",
+    };
+    for (const auto &smi : smiles) {
+      INFO(smi);
+      auto m = v2::SmilesParse::MolFromSmiles(smi);
+      REQUIRE(m);
+      int charge = 0;
+      for (auto atom : m->atoms()) {
+        charge += atom->getFormalCharge();
+      }
+      bool allowChargedFragments = true;
+      bool embedChiral = false;
+      RWMol m2(*m);
+      determineBondOrders(m2, charge, allowChargedFragments, embedChiral);
+      for (auto bnd : m2.bonds()) {
+        CHECK(bnd->getBondType() ==
+              m->getBondWithIdx(bnd->getIdx())->getBondType());
+      }
+    }
+  }
 }