Fixing pickle bug: Number of rings now maintained when unpickling into existing molecules

Code/DataStructs/ExplicitBitVect.h

@@ -28,7 +28,9 @@
  */
 class RDKIT_DATASTRUCTS_EXPORT ExplicitBitVect : public BitVect {
  public:
-  ExplicitBitVect() {}
+  ExplicitBitVect() : dp_bits(nullptr), d_size(0), d_numOnBits(0) {
+    _initForSize(0);
+  }
   //! initialize with a particular size;
   explicit ExplicitBitVect(unsigned int size)
       : dp_bits(nullptr), d_size(0), d_numOnBits(0) {

Code/GraphMol/Abbreviations/Abbreviations.cpp

@@ -25,6 +25,8 @@ void applyMatches(RWMol &mol, const std::vector<AbbreviationMatch> &matches) {
   std::vector<unsigned int> prevBondMapping;
   std::vector<unsigned int> addedBonds;
   addedBonds.reserve(mol.getNumBonds());
+  // Store bonds to add - pairs of (atom1, atom2, originalBondIdx)
+  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>> bondsToAdd;
   bool hasPrevMapping =
       mol.getPropIfPresent(common_properties::origAtomMapping,
                            prevAtomMapping) &&
@@ -74,10 +76,10 @@ void applyMatches(RWMol &mol, const std::vector<AbbreviationMatch> &matches) {
                            [nbrIdx](const std::pair<int, int> &tpr) {
                              return tpr.second == rdcast<int>(nbrIdx);
                            })) {
-            mol.addBond(nbrIdx, connectIdx, Bond::BondType::SINGLE);
-            addedBonds.push_back(hasPrevMapping
-                                     ? prevBondMapping.at(bond->getIdx())
-                                     : bond->getIdx());
+            // Defer bond addition to avoid iterator invalidation
+            bondsToAdd.emplace_back(nbrIdx, connectIdx,
+                                    hasPrevMapping ? prevBondMapping.at(bond->getIdx())
+                                                   : bond->getIdx());
           }
         }
       }
@@ -92,11 +94,19 @@ void applyMatches(RWMol &mol, const std::vector<AbbreviationMatch> &matches) {
         }
       }
       CHECK_INVARIANT(bondIdx != -1, "bondIdx not found");
-      mol.addBond(oaidx, connectIdx, Bond::BondType::SINGLE);
-      addedBonds.push_back(hasPrevMapping ? prevBondMapping.at(bondIdx)
-                                          : bondIdx);
+      // Defer bond addition to avoid iterator invalidation
+      bondsToAdd.emplace_back(oaidx, connectIdx,
+                              hasPrevMapping ? prevBondMapping.at(bondIdx)
+                                             : bondIdx);
     }
   }
+
+  // Now add all the deferred bonds, checking for duplicates
+  for (const auto &[atom1, atom2, origBondIdx] : bondsToAdd) {
+    mol.addBond(atom1, atom2, Bond::BondType::SINGLE);
+    addedBonds.push_back(origBondIdx);
+  }
+
   std::vector<unsigned int> atomMapping;
   atomMapping.reserve(mol.getNumAtoms());
   std::vector<unsigned int> bondMapping;

Code/GraphMol/AddHs.cpp

@@ -549,12 +549,27 @@ void addHs(RWMol &mol, const AddHsParameters &params,
   } else {
     onAtoms.set();
   }
+  std::vector<unsigned int> numExplicitHs(mol.getNumAtoms(), 0);
+  std::vector<unsigned int> numImplicitHs(mol.getNumAtoms(), 0);
+
+  // Ensure all atoms have implicit valence calculated (e.g., after unpickling)
+  // If noImplicit is false, the atom may need implicit valence calculated
+  // Skip query atoms (or atoms connected to query bonds) as they are handled separately
+  for (auto at : mol.atoms()) {
+    if (!at->getNoImplicit() && !isQueryAtom(mol, *at)) {
+      at->updatePropertyCache(false);
+    }
+  }
+
   for (auto at : mol.atoms()) {
+    numExplicitHs[at->getIdx()] = at->getNumExplicitHs();
+    numImplicitHs[at->getIdx()] = at->getNumImplicitHs();
     if (onAtoms[at->getIdx()]) {
       if (params.skipQueries && isQueryAtom(mol, *at)) {
         onAtoms.set(at->getIdx(), 0);
         continue;
       }
+
       numAddHyds += at->getNumExplicitHs();
       if (!params.explicitOnly) {
         numAddHyds += at->getNumImplicitHs();
@@ -586,7 +601,7 @@ void addHs(RWMol &mol, const AddHsParameters &params,
     unsigned int newIdx;
     newAt->clearComputedProps();
     // always convert explicit Hs
-    unsigned int onumexpl = newAt->getNumExplicitHs();
+    unsigned int onumexpl = numExplicitHs[aidx];
     for (unsigned int i = 0; i < onumexpl; i++) {
       newIdx = mol.addAtom(new Atom(1), false, true);
       mol.addBond(aidx, newIdx, Bond::SINGLE);
@@ -605,8 +620,7 @@ void addHs(RWMol &mol, const AddHsParameters &params,
 
     if (!params.explicitOnly) {
       // take care of implicits
-      for (unsigned int i = 0; i < mol.getAtomWithIdx(aidx)->getNumImplicitHs();
-           i++) {
+      for (unsigned int i = 0; i < mol.getAtomWithIdx(aidx)->getNumImplicitHs(); i++) {
         newIdx = mol.addAtom(new Atom(1), false, true);
         mol.addBond(aidx, newIdx, Bond::SINGLE);
         // set the isImplicit label so that we can strip these back

Code/GraphMol/AdjustQuery.cpp

@@ -349,16 +349,16 @@ void adjustQueryProperties(RWMol &mol, const AdjustQueryParameters *inParams) {
   if (inParams) {
     params = *inParams;
   }
-  const auto ringInfo = mol.getRingInfo();
-
   if (params.aromatizeIfPossible) {
     unsigned int failed;
     sanitizeMol(mol, failed, SANITIZE_SYMMRINGS | SANITIZE_SETAROMATICITY);
   } else {
+    const auto ringInfo = mol.getRingInfo();
     if (!ringInfo->isSymmSssr()) {
       MolOps::symmetrizeSSSR(mol);
     }
   }
+  const auto ringInfo = mol.getRingInfo();
   QueryAtom qaTmpl;
   QueryBond qbTmpl;
 

Code/GraphMol/Atom.cpp

@@ -286,14 +286,14 @@ void Atom::setAtomicNum(int newNum) {
 }
 std::string Atom::getSymbol() const {
   int atomicNum = getAtomicNum();
-  if (atomicNum == 0) {
-    // handle dummies differently:
-    std::string res;
-    dp_dataMol->getAtomPropIfPresent<std::string>(
-        common_properties::dummyLabelToken, d_index, res);
-    return res;
+  std::string res;
+  // handle dummies differently:
+  if (atomicNum != 0 ||
+      !dp_dataMol->getAtomPropIfPresent<std::string>(
+          common_properties::dummyLabelToken, d_index, res)) {
+    res = PeriodicTable::getTable()->getElementSymbol(atomicNum);
   }
-  return PeriodicTable::getTable()->getElementSymbol(atomicNum);
+  return res;
 }
 ROMol &Atom::getOwningMol() const {
   PRECONDITION(dp_owningMol != nullptr, "no owner");
@@ -412,13 +412,14 @@ bool Atom::Match(Atom const* what) const {
   //   [*] matches [*],[1*],[2*],etc.
   //   [1*] only matches [*] and [1*]
   if (res) {
-    if (hasOwningMol() && what->hasOwningMol() &&
-        this->getOwningMol().getRingInfo()->isInitialized() &&
-        what->getOwningMol().getRingInfo()->isInitialized() &&
-        this->getOwningMol().getRingInfo()->numAtomRings(d_index) >
-            what->getOwningMol().getRingInfo()->numAtomRings(what->d_index)) {
-      res = false;
-    } else if (!this->getAtomicNum()) {
+    // if (hasOwningMol() && what->hasOwningMol() &&
+        // this->getOwningMol().getRingInfo()->isInitialized() &&
+        // what->getOwningMol().getRingInfo()->isInitialized() &&
+        // this->getOwningMol().getRingInfo()->numAtomRings(d_index) >
+            // what->getOwningMol().getRingInfo()->numAtomRings(what->d_index)) {
+      // res = false;
+  // } else if (!this->getAtomicNum()) {
+  if (!this->getAtomicNum()) {
       // this is the new behavior, based on the isotopes:
       int tgt = this->getIsotope();
       int test = what->getIsotope();

Code/GraphMol/Bond.cpp

@@ -101,8 +101,6 @@ void Bond::setOwningMol(ROMol *other) {
 }
 
 void Bond::setOwningMol(RDMol* other) {
-  PRECONDITION(dp_owningMol == nullptr || dp_owningMol == other,
-               "setOwningMol called twice");
   dp_owningMol = other;
 }
 

Code/GraphMol/ChemReactions/ReactionRunner.cpp

@@ -1546,7 +1546,7 @@ generateOneProductSet(const ChemicalReaction &rxn,
     if (!(*pTemplIt)->getStereoGroups().empty()) {
       copyTemplateStereoGroupsToMol(**pTemplIt, product);
     }
-
+    product->updatePropertyCache(false);
     res[prodId] = product;
     ++prodId;
   }

Code/GraphMol/ChemTransforms/ChemTransforms.cpp

@@ -227,19 +227,25 @@ std::vector<ROMOL_SPTR> replaceSubstructs(
 
     // loop over the central atom's (the first atom in match) bonds
     // and duplicate any that connect to the remainder of the molecule:
+    // Collect neighbor indices first since addBond will invalidate iterators
     Atom *origAtom = newMol->getAtomWithIdx(match[0]);
+    std::vector<unsigned int> neighbors;
     ROMol::ADJ_ITER nbrIdx, endNbrs;
     boost::tie(nbrIdx, endNbrs) = newMol->getAtomNeighbors(origAtom);
     while (nbrIdx != endNbrs) {
+      neighbors.push_back(*nbrIdx);
+      nbrIdx++;
+    }
+
+    for (auto neighborIdx : neighbors) {
       // we don't want to duplicate any "intra-match" bonds:
       if (!std::binary_search(sortMatch.begin(), sortMatch.end(),
-                              int(*nbrIdx))) {
-        Bond *oBond = newMol->getBondBetweenAtoms(match[0], *nbrIdx);
+                              int(neighborIdx))) {
+        Bond *oBond = newMol->getBondBetweenAtoms(match[0], neighborIdx);
         CHECK_INVARIANT(oBond, "required bond not found");
-        newMol->addBond(numOrigAtoms + replacementConnectionPoint, *nbrIdx,
+        newMol->addBond(numOrigAtoms + replacementConnectionPoint, neighborIdx,
                         oBond->getBondType());
       }
-      nbrIdx++;
     }
 
     if (replaceAll) {

Code/GraphMol/FileParsers/MaeWriter.cpp

@@ -160,14 +160,28 @@ void copyProperties(
       continue;
     }
 
-    switch (prop.getRDValueTag()) {
+    // Get the RDValue tag based on scope
+    auto tag = (scope == RDMol::Scope::MOL) ? prop.getRDValueTag() : prop.getRDValueTag(idx);
+
+    // Helper lambda to get property value based on scope
+    auto getPropValue = [&](const auto& token) -> RDValue {
+      if (scope == RDMol::Scope::MOL) {
+        return mol.getMolProp<RDValue>(token);
+      } else if (scope == RDMol::Scope::ATOM) {
+        return mol.getAtomProp<RDValue>(token, idx);
+      } else { // BOND
+        return mol.getBondProp<RDValue>(token, idx);
+      }
+    };
+
+    switch (tag) {
       case RDTypeTag::BoolTag: {
         auto propName = prop.name().getString();
         if (!std::regex_match(prop.name().getString(), MMCT_PROP_REGEX)) {
           propName.insert(0, "b_rdkit_");
         }
 
-        bool v = mol.getAtomProp<bool>(prop.name(), idx);
+        bool v = rdvalue_cast<bool>(getPropValue(prop.name()));
         boolSetter(propName, idx, v);
         break;
       }
@@ -182,7 +196,7 @@ void copyProperties(
           propName.insert(0, "i_rdkit_");
         }
 
-        int v = mol.getAtomProp<int>(prop.name(), idx);
+        int v = rdvalue_cast<int>(getPropValue(prop.name()));
         intSetter(propName, idx, v);
         break;
       }
@@ -194,7 +208,7 @@ void copyProperties(
           propName.insert(0, "r_rdkit_");
         }
 
-        double v = mol.getAtomProp<double>(prop.name(), idx);
+        double v = rdvalue_cast<double>(getPropValue(prop.name()));
         realSetter(propName, idx, v);
         break;
       }
@@ -205,7 +219,7 @@ void copyProperties(
           propName.insert(0, "s_rdkit_");
         }
 
-        std::string v = mol.getAtomProp<std::string>(prop.name(), idx);
+        std::string v = rdvalue_cast<std::string>(getPropValue(prop.name()));
         stringSetter(propName, idx, v);
         break;
       }

Code/GraphMol/FileParsers/MolFileParser.cpp

@@ -1348,7 +1348,7 @@ void ParseAtomAlias(RWMol *mol, std::string text, const std::string &nextLine,
   }
   URANGE_CHECK(idx, mol->getNumAtoms());
   Atom *at = mol->getAtomWithIdx(idx);
-  at->setProp(common_properties::molFileAlias, nextLine);
+  setAtomAlias(at, nextLine);
 }
 
 void ParseAtomValue(RWMol *mol, std::string text, unsigned int line) {
@@ -1367,8 +1367,7 @@ void ParseAtomValue(RWMol *mol, std::string text, unsigned int line) {
   }
   URANGE_CHECK(idx, mol->getNumAtoms());
   Atom *at = mol->getAtomWithIdx(idx);
-  at->setProp(common_properties::molFileValue,
-              text.substr(7, text.length() - 7));
+  setAtomValue(at, text.substr(7, text.length() - 7));
 }
 
 namespace {