diff --git a/Code/GraphMol/ChemReactions/MDLParser.cpp b/Code/GraphMol/ChemReactions/MDLParser.cpp
index 3d9fa7026..8975fcc62 100644
--- a/Code/GraphMol/ChemReactions/MDLParser.cpp
+++ b/Code/GraphMol/ChemReactions/MDLParser.cpp
@@ -1,5 +1,5 @@
 //
-//  Copyright (c) 2007-2024, Novartis Institutes for BioMedical Research Inc.
+//  Copyright (c) 2007-2026, Novartis Institutes for BioMedical Research Inc.
 //  and other RDKit contributors
 //
 //  All rights reserved.
@@ -350,17 +350,15 @@ std::unique_ptr<ChemicalReaction> ReactionFromRxnDataStream(
   for (MOL_SPTR_VECT::const_iterator iter = res->beginReactantTemplates();
        iter != res->endReactantTemplates(); ++iter) {
     // to write the mol block, we need ring information:
-    for (ROMol::AtomIterator atomIt = (*iter)->beginAtoms();
-         atomIt != (*iter)->endAtoms(); ++atomIt) {
-      QueryOps::replaceAtomWithQueryAtom((RWMol *)iter->get(), (*atomIt));
+    for (auto atom : (*iter)->atoms()) {
+      QueryOps::replaceAtomWithQueryAtom((RWMol *)iter->get(), atom);
     }
   }
   for (MOL_SPTR_VECT::const_iterator iter = res->beginProductTemplates();
        iter != res->endProductTemplates(); ++iter) {
     // to write the mol block, we need ring information:
-    for (ROMol::AtomIterator atomIt = (*iter)->beginAtoms();
-         atomIt != (*iter)->endAtoms(); ++atomIt) {
-      QueryOps::replaceAtomWithQueryAtom((RWMol *)iter->get(), (*atomIt));
+    for (auto atom : (*iter)->atoms()) {
+      QueryOps::replaceAtomWithQueryAtom((RWMol *)iter->get(), atom);
     }
   }
   updateProductsStereochem(res.get());
diff --git a/Code/GraphMol/ChemReactions/Reaction.cpp b/Code/GraphMol/ChemReactions/Reaction.cpp
index ec1cdef45..6021dbb2d 100644
--- a/Code/GraphMol/ChemReactions/Reaction.cpp
+++ b/Code/GraphMol/ChemReactions/Reaction.cpp
@@ -107,11 +107,9 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
   for (auto molIter = this->beginReactantTemplates();
        molIter != this->endReactantTemplates(); ++molIter) {
     bool thisMolMapped = false;
-    for (ROMol::AtomIterator atomIt = (*molIter)->beginAtoms();
-         atomIt != (*molIter)->endAtoms(); ++atomIt) {
+    for (const auto atom : (*molIter)->atoms()) {
       int mapNum;
-      if ((*atomIt)->getPropIfPresent(common_properties::molAtomMapNumber,
-                                      mapNum)) {
+      if (atom->getPropIfPresent(common_properties::molAtomMapNumber, mapNum)) {
         thisMolMapped = true;
         if (std::find(mapNumbersSeen.begin(), mapNumbersSeen.end(), mapNum) !=
             mapNumbersSeen.end()) {
@@ -123,7 +121,7 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
           res = false;
         } else {
           mapNumbersSeen.push_back(mapNum);
-          reactingAtoms[mapNum] = *atomIt;
+          reactingAtoms[mapNum] = atom;
         }
       }
     }
@@ -143,27 +141,24 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
        molIter != this->endProductTemplates(); ++molIter) {
     // clear out some possible cached properties to prevent
     // misleading warnings
-    for (ROMol::AtomIterator atomIt = (*molIter)->beginAtoms();
-         atomIt != (*molIter)->endAtoms(); ++atomIt) {
-      if ((*atomIt)->hasProp(common_properties::_QueryFormalCharge)) {
-        (*atomIt)->clearProp(common_properties::_QueryFormalCharge);
+    for (auto atom : (*molIter)->atoms()) {
+      if (atom->hasProp(common_properties::_QueryFormalCharge)) {
+        atom->clearProp(common_properties::_QueryFormalCharge);
       }
-      if ((*atomIt)->hasProp(common_properties::_QueryHCount)) {
-        (*atomIt)->clearProp(common_properties::_QueryHCount);
+      if (atom->hasProp(common_properties::_QueryHCount)) {
+        atom->clearProp(common_properties::_QueryHCount);
       }
-      if ((*atomIt)->hasProp(common_properties::_QueryMass)) {
-        (*atomIt)->clearProp(common_properties::_QueryMass);
+      if (atom->hasProp(common_properties::_QueryMass)) {
+        atom->clearProp(common_properties::_QueryMass);
       }
-      if ((*atomIt)->hasProp(common_properties::_QueryIsotope)) {
-        (*atomIt)->clearProp(common_properties::_QueryIsotope);
+      if (atom->hasProp(common_properties::_QueryIsotope)) {
+        atom->clearProp(common_properties::_QueryIsotope);
       }
     }
     bool thisMolMapped = false;
-    for (ROMol::AtomIterator atomIt = (*molIter)->beginAtoms();
-         atomIt != (*molIter)->endAtoms(); ++atomIt) {
+    for (auto atom : (*molIter)->atoms()) {
       int mapNum;
-      if ((*atomIt)->getPropIfPresent(common_properties::molAtomMapNumber,
-                                      mapNum)) {
+      if (atom->getPropIfPresent(common_properties::molAtomMapNumber, mapNum)) {
         thisMolMapped = true;
         bool seenAlready =
             std::find(productNumbersSeen.begin(), productNumbersSeen.end(),
@@ -180,8 +175,8 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
           // ------------
           const Atom *rAtom = reactingAtoms[mapNum];
           CHECK_INVARIANT(rAtom, "missing atom");
-          if (rAtom->getDegree() != (*atomIt)->getDegree()) {
-            (*atomIt)->setProp(common_properties::_ReactionDegreeChanged, 1);
+          if (rAtom->getDegree() != atom->getDegree()) {
+            atom->setProp(common_properties::_ReactionDegreeChanged, 1);
           }
 
         } else {
@@ -205,8 +200,8 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
           // ------------
           const Atom *rAtom = reactingAtoms[mapNum];
           CHECK_INVARIANT(rAtom, "missing atom");
-          if (rAtom->getDegree() != (*atomIt)->getDegree()) {
-            (*atomIt)->setProp(common_properties::_ReactionDegreeChanged, 1);
+          if (rAtom->getDegree() != atom->getDegree()) {
+            atom->setProp(common_properties::_ReactionDegreeChanged, 1);
           }
         }
       }
@@ -214,9 +209,9 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
       // ------------
       //    Deal with queries
       // ------------
-      if ((*atomIt)->hasQuery()) {
+      if (atom->hasQuery()) {
         std::list<const Atom::QUERYATOM_QUERY *> queries;
-        queries.push_back((*atomIt)->getQuery());
+        queries.push_back(atom->getQuery());
         while (!queries.empty()) {
           const Atom::QUERYATOM_QUERY *query = queries.front();
           queries.pop_front();
@@ -229,74 +224,72 @@ bool ChemicalReaction::validate(unsigned int &numWarnings,
             int qval;
             int neg =
                 query->getDescription() == "AtomNegativeFormalCharge" ? -1 : 1;
-            if ((*atomIt)->getPropIfPresent(
-                    common_properties::_QueryFormalCharge, qval) &&
+            if (atom->getPropIfPresent(common_properties::_QueryFormalCharge,
+                                       qval) &&
                 (neg * qval) !=
                     static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal()) {
               if (!silent) {
                 BOOST_LOG(rdWarningLog)
-                    << "atom " << (*atomIt)->getIdx() << " in product "
-                    << molIdx << " has multiple charge specifications.\n";
+                    << "atom " << atom->getIdx() << " in product " << molIdx
+                    << " has multiple charge specifications.\n";
               }
               numWarnings++;
             } else {
               int neg = query->getDescription() == "AtomNegativeFormalCharge"
                             ? -1
                             : 1;
-              (*atomIt)->setProp(
+              atom->setProp(
                   common_properties::_QueryFormalCharge,
                   neg *
                       static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal());
             }
           } else if (query->getDescription() == "AtomHCount") {
             int qval;
-            if ((*atomIt)->getPropIfPresent(common_properties::_QueryHCount,
-                                            qval) &&
+            if (atom->getPropIfPresent(common_properties::_QueryHCount, qval) &&
                 qval !=
                     static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal()) {
               if (!silent) {
                 BOOST_LOG(rdWarningLog)
-                    << "atom " << (*atomIt)->getIdx() << " in product "
-                    << molIdx << " has multiple H count specifications.\n";
+                    << "atom " << atom->getIdx() << " in product " << molIdx
+                    << " has multiple H count specifications.\n";
               }
               numWarnings++;
             } else {
-              (*atomIt)->setProp(
+              atom->setProp(
                   common_properties::_QueryHCount,
                   static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal());
             }
           } else if (query->getDescription() == "AtomMass") {
             int qval;
-            if ((*atomIt)->getPropIfPresent(common_properties::_QueryMass,
-                                            qval) &&
+            if (atom->getPropIfPresent(common_properties::_QueryMass, qval) &&
                 qval !=
                     static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal()) {
               if (!silent) {
                 BOOST_LOG(rdWarningLog)
-                    << "atom " << (*atomIt)->getIdx() << " in product "
-                    << molIdx << " has multiple mass specifications.\n";
+                    << "atom " << atom->getIdx() << " in product " << molIdx
+                    << " has multiple mass specifications.\n";
               }
               numWarnings++;
             } else {
-              (*atomIt)->setProp(
+              atom->setProp(
                   common_properties::_QueryMass,
                   static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal() /
                       massIntegerConversionFactor);
             }
           } else if (query->getDescription() == "AtomIsotope") {
             int qval;
-            if ((*atomIt)->getPropIfPresent(common_properties::_QueryIsotope,
-                                            qval) &&
+            if (atom->getPropIfPresent(common_properties::_QueryIsotope,
+                                       qval) &&
                 qval !=
                     static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal()) {
               if (!silent) {
                 BOOST_LOG(rdWarningLog)
-                    << "atom " << (*atomIt)->getIdx() << " in product "
-                    << molIdx << " has multiple isotope specifications.\n";
+                    << "atom " << atom->getIdx() << " in product " << molIdx
+                    << " has multiple isotope specifications.\n";
               }
               numWarnings++;
             } else {
-              (*atomIt)->setProp(
+              atom->setProp(
                   common_properties::_QueryIsotope,
                   static_cast<const ATOM_EQUALS_QUERY *>(query)->getVal());
             }
diff --git a/Code/GraphMol/ChemReactions/ReactionUtils.cpp b/Code/GraphMol/ChemReactions/ReactionUtils.cpp
index 69b50d57d..57380d5c2 100644
--- a/Code/GraphMol/ChemReactions/ReactionUtils.cpp
+++ b/Code/GraphMol/ChemReactions/ReactionUtils.cpp
@@ -340,10 +340,9 @@ void removeMappingNumbersFromReactionMoleculeTemplate(
     const MOL_SPTR_VECT &molVec) {
   for (const auto &begin : molVec) {
     ROMol &mol = *begin.get();
-    for (ROMol::AtomIterator atomIt = mol.beginAtoms();
-         atomIt != mol.endAtoms(); ++atomIt) {
-      if ((*atomIt)->hasProp(common_properties::molAtomMapNumber)) {
-        (*atomIt)->clearProp(common_properties::molAtomMapNumber);
+    for (auto atom : mol.atoms()) {
+      if (atom->hasProp(common_properties::molAtomMapNumber)) {
+        atom->clearProp(common_properties::molAtomMapNumber);
       }
     }
   }
diff --git a/Code/GraphMol/Chirality.cpp b/Code/GraphMol/Chirality.cpp
index 517663afe..275583150 100644
--- a/Code/GraphMol/Chirality.cpp
+++ b/Code/GraphMol/Chirality.cpp
@@ -2924,12 +2924,10 @@ void findPotentialStereoBonds(ROMol &mol, bool cleanIt) {
     ranks.resize(mol.getNumAtoms());
     bool cipDone = false;
 
-    ROMol::BondIterator bondIt;
-    for (bondIt = mol.beginBonds(); bondIt != mol.endBonds(); ++bondIt) {
-      if ((*bondIt)->getBondType() == Bond::DOUBLE &&
-          !(mol.getRingInfo()->numBondRings((*bondIt)->getIdx()))) {
+    for (auto dblBond : mol.bonds()) {
+      if (dblBond->getBondType() == Bond::DOUBLE &&
+          !(mol.getRingInfo()->numBondRings(dblBond->getIdx()))) {
         // we are ignoring ring bonds here - read the FIX above
-        Bond *dblBond = *bondIt;
         // proceed only if we either want to clean the stereocode on this bond,
         // if none is set on it yet, or it is STEREOANY and we need to find
         // stereoatoms
@@ -3024,10 +3022,10 @@ void findPotentialStereoBonds(ROMol &mol, bool cleanIt) {
               }
             }  // end of check that beg and end atoms have at least 1
                // neighbor:
-          }    // end of 2 and 3 coordinated atoms only
-        }      // end of we want it or CIP code is not set
-      }        // end of double bond
-    }          // end of for loop over all bonds
+          }  // end of 2 and 3 coordinated atoms only
+        }  // end of we want it or CIP code is not set
+      }  // end of double bond
+    }  // end of for loop over all bonds
     mol.setProp(common_properties::_BondsPotentialStereo, 1, true);
   }
 }
@@ -3661,7 +3659,7 @@ void setDoubleBondNeighborDirections(ROMol &mol, const Conformer *conf) {
 
   // oof, now loop over the double bonds in that order and
   // update their neighbor directionalities:
-  for (const auto& pairIter : orderedBondsInPlay) {
+  for (const auto &pairIter : orderedBondsInPlay) {
     updateDoubleBondNeighbors(mol, pairIter.second, conf, needsDir,
                               singleBondCounts, singleBondNbrs);
   }
diff --git a/Code/GraphMol/Descriptors/Lipinski.cpp b/Code/GraphMol/Descriptors/Lipinski.cpp
index fc8d48dee..a64fd400d 100644
--- a/Code/GraphMol/Descriptors/Lipinski.cpp
+++ b/Code/GraphMol/Descriptors/Lipinski.cpp
@@ -76,9 +76,8 @@ namespace Descriptors {
 
 unsigned int calcLipinskiHBA(const ROMol &mol) {
   unsigned int res = 0;
-  for (ROMol::ConstAtomIterator iter = mol.beginAtoms(); iter != mol.endAtoms();
-       ++iter) {
-    if ((*iter)->getAtomicNum() == 7 || (*iter)->getAtomicNum() == 8) {
+  for (const auto atom : mol.atoms()) {
+    if (atom->getAtomicNum() == 7 || atom->getAtomicNum() == 8) {
       ++res;
     }
   }
@@ -87,10 +86,9 @@ unsigned int calcLipinskiHBA(const ROMol &mol) {
 
 unsigned int calcLipinskiHBD(const ROMol &mol) {
   unsigned int res = 0;
-  for (ROMol::ConstAtomIterator iter = mol.beginAtoms(); iter != mol.endAtoms();
-       ++iter) {
-    if (((*iter)->getAtomicNum() == 7 || (*iter)->getAtomicNum() == 8)) {
-      res += (*iter)->getTotalNumHs(true);
+  for (const auto atom : mol.atoms()) {
+    if (atom->getAtomicNum() == 7 || atom->getAtomicNum() == 8) {
+      res += atom->getTotalNumHs(true);
     }
   }
   return res;
diff --git a/Code/GraphMol/Descriptors/MolSurf.cpp b/Code/GraphMol/Descriptors/MolSurf.cpp
index 29bcf0300..03c712bce 100644
--- a/Code/GraphMol/Descriptors/MolSurf.cpp
+++ b/Code/GraphMol/Descriptors/MolSurf.cpp
@@ -40,22 +40,21 @@ double getLabuteAtomContribs(const ROMol &mol, std::vector<double> &Vi,
     Vi[i] = 0.0;
   }
 
-  for (ROMol::ConstBondIterator bondIt = mol.beginBonds();
-       bondIt != mol.endBonds(); ++bondIt) {
+  for (const auto bond : mol.bonds()) {
     const double bondScaleFacts[4] = {.1, 0, .2, .3};
-    double Ri = rads[(*bondIt)->getBeginAtomIdx()];
-    double Rj = rads[(*bondIt)->getEndAtomIdx()];
+    double Ri = rads[bond->getBeginAtomIdx()];
+    double Rj = rads[bond->getEndAtomIdx()];
     double bij = Ri + Rj;
-    if (!(*bondIt)->getIsAromatic()) {
-      if ((*bondIt)->getBondType() < 4) {
-        bij -= bondScaleFacts[(*bondIt)->getBondType()];
+    if (!bond->getIsAromatic()) {
+      if (bond->getBondType() < 4) {
+        bij -= bondScaleFacts[bond->getBondType()];
       }
     } else {
       bij -= bondScaleFacts[0];
     }
     double dij = std::min(std::max(fabs(Ri - Rj), bij), Ri + Rj);
-    Vi[(*bondIt)->getBeginAtomIdx()] += Rj * Rj - (Ri - dij) * (Ri - dij) / dij;
-    Vi[(*bondIt)->getEndAtomIdx()] += Ri * Ri - (Rj - dij) * (Rj - dij) / dij;
+    Vi[bond->getBeginAtomIdx()] += Rj * Rj - (Ri - dij) * (Ri - dij) / dij;
+    Vi[bond->getEndAtomIdx()] += Ri * Ri - (Rj - dij) * (Rj - dij) / dij;
   }
   hContrib = 0.0;
   if (includeHs) {
@@ -118,31 +117,29 @@ double getTPSAAtomContribs(const ROMol &mol, std::vector<double> &Vi,
   unsigned int nAtoms = mol.getNumAtoms();
   std::vector<int> nNbrs(nAtoms, 0), nSing(nAtoms, 0), nDoub(nAtoms, 0),
       nTrip(nAtoms, 0), nArom(nAtoms, 0), nHs(nAtoms, 0);
-  for (ROMol::ConstBondIterator bIt = mol.beginBonds(); bIt != mol.endBonds();
-       ++bIt) {
-    const Bond *bnd = (*bIt);
-    if (bnd->getBeginAtom()->getAtomicNum() == 1) {
-      nNbrs[bnd->getEndAtomIdx()] -= 1;
-      nHs[bnd->getEndAtomIdx()] += 1;
-    } else if (bnd->getEndAtom()->getAtomicNum() == 1) {
-      nNbrs[bnd->getBeginAtomIdx()] -= 1;
-      nHs[bnd->getBeginAtomIdx()] += 1;
-    } else if (bnd->getIsAromatic()) {
-      nArom[bnd->getBeginAtomIdx()] += 1;
-      nArom[bnd->getEndAtomIdx()] += 1;
+  for (const auto bond : mol.bonds()) {
+    if (bond->getBeginAtom()->getAtomicNum() == 1) {
+      nNbrs[bond->getEndAtomIdx()] -= 1;
+      nHs[bond->getEndAtomIdx()] += 1;
+    } else if (bond->getEndAtom()->getAtomicNum() == 1) {
+      nNbrs[bond->getBeginAtomIdx()] -= 1;
+      nHs[bond->getBeginAtomIdx()] += 1;
+    } else if (bond->getIsAromatic()) {
+      nArom[bond->getBeginAtomIdx()] += 1;
+      nArom[bond->getEndAtomIdx()] += 1;
     } else {
-      switch (bnd->getBondType()) {
+      switch (bond->getBondType()) {
         case Bond::SINGLE:
-          nSing[bnd->getBeginAtomIdx()] += 1;
-          nSing[bnd->getEndAtomIdx()] += 1;
+          nSing[bond->getBeginAtomIdx()] += 1;
+          nSing[bond->getEndAtomIdx()] += 1;
           break;
         case Bond::DOUBLE:
-          nDoub[bnd->getBeginAtomIdx()] += 1;
-          nDoub[bnd->getEndAtomIdx()] += 1;
+          nDoub[bond->getBeginAtomIdx()] += 1;
+          nDoub[bond->getEndAtomIdx()] += 1;
           break;
         case Bond::TRIPLE:
-          nTrip[bnd->getBeginAtomIdx()] += 1;
-          nTrip[bnd->getEndAtomIdx()] += 1;
+          nTrip[bond->getBeginAtomIdx()] += 1;
+          nTrip[bond->getEndAtomIdx()] += 1;
           break;
         default:
           break;
diff --git a/Code/GraphMol/Descriptors/PMI.cpp b/Code/GraphMol/Descriptors/PMI.cpp
index cc4c85512..c2628682a 100644
--- a/Code/GraphMol/Descriptors/PMI.cpp
+++ b/Code/GraphMol/Descriptors/PMI.cpp
@@ -42,9 +42,8 @@ bool getMoments(const ROMol &mol, int confId, bool useAtomicMasses, double &pm1,
   if (useAtomicMasses) {
     std::vector<double> weights;
     weights.resize(mol.getNumAtoms());
-    for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
-         ++cai) {
-      weights[(*cai)->getIdx()] = (*cai)->getMass();
+    for (const auto atom : mol.atoms()) {
+      weights[atom->getIdx()] = atom->getMass();
     }
     res = MolTransforms::computePrincipalAxesAndMoments(
         conf, axes, moments, ignoreHs, force, &weights);
@@ -85,9 +84,8 @@ bool getMomentsFromGyration(const ROMol &mol, int confId, bool useAtomicMasses,
   if (useAtomicMasses) {
     std::vector<double> weights;
     weights.resize(mol.getNumAtoms());
-    for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
-         ++cai) {
-      weights[(*cai)->getIdx()] = (*cai)->getMass();
+    for (const auto atom : mol.atoms()) {
+      weights[atom->getIdx()] = atom->getMass();
     }
     res = MolTransforms::computePrincipalAxesAndMomentsFromGyrationMatrix(
         conf, axes, moments, ignoreHs, force, &weights);
diff --git a/Code/GraphMol/FileParsers/CMLWriter.cpp b/Code/GraphMol/FileParsers/CMLWriter.cpp
index 12d32273d..5c0b26a9b 100644
--- a/Code/GraphMol/FileParsers/CMLWriter.cpp
+++ b/Code/GraphMol/FileParsers/CMLWriter.cpp
@@ -156,14 +156,10 @@ boost::property_tree::ptree molToPTree(const ROMol &mol, int confId,
   unsigned bond_id = 0u;
 
   auto &bondArray = molecule.add("bondArray", "");
-  for (auto atom_itr = rwmol.beginAtoms(), atom_itr_end = rwmol.endAtoms();
-       atom_itr != atom_itr_end; ++atom_itr) {
-    const auto &atom = *atom_itr;
+  for (auto atom : rwmol.atoms()) {
     PRECONDITION(atom, "bad atom");
     const auto src = atom->getIdx();
-    for (auto bond_itrs = rwmol.getAtomBonds(atom);
-         bond_itrs.first != bond_itrs.second; ++bond_itrs.first) {
-      auto *bptr = rwmol[*bond_itrs.first];
+    for (auto bptr : rwmol.atomBonds(atom)) {
       auto *nptr = bptr->getOtherAtom(atom);
       const auto dst = nptr->getIdx();
       if (dst < src) {
diff --git a/Code/GraphMol/FileParsers/Mol2FileParser.cpp b/Code/GraphMol/FileParsers/Mol2FileParser.cpp
index 912e1879e..8f10cdf00 100644
--- a/Code/GraphMol/FileParsers/Mol2FileParser.cpp
+++ b/Code/GraphMol/FileParsers/Mol2FileParser.cpp
@@ -138,9 +138,7 @@ void readFormalChargesFromAttr(std::istream *inStream, RWMol *res) {
 
 void guessFormalCharges(RWMol *res) {
   // FIX: this whole thing has problems with positively charged pyridines et al.
-  for (RWMol::AtomIterator atomIt = res->beginAtoms();
-       atomIt != res->endAtoms(); ++atomIt) {
-    Atom *at = (*atomIt);
+  for (auto at : res->atoms()) {
     // assign only if no formal charge set on that atom and atom is not carbon
     // (the latter
     // might needs changing later on - let's see) and not for query atoms (dummy
@@ -951,10 +949,11 @@ std::unique_ptr<RWMol> MolFromMol2DataStream(std::istream &inStream,
     MolOps::cleanUp(*res);
 
     try {
-      // when we sanitize for mol2, we skip the cleanup organometallic step since it's
-      // not really compatible with the semantics of mol2 files
-      constexpr auto sanitizeFlags = MolOps::SanitizeFlags::SANITIZE_ALL ^
-                            MolOps::SanitizeFlags::SANITIZE_CLEANUP_ORGANOMETALLICS;
+      // when we sanitize for mol2, we skip the cleanup organometallic step
+      // since it's not really compatible with the semantics of mol2 files
+      constexpr auto sanitizeFlags =
+          MolOps::SanitizeFlags::SANITIZE_ALL ^
+          MolOps::SanitizeFlags::SANITIZE_CLEANUP_ORGANOMETALLICS;
       if (params.removeHs) {
         // Bond stereo detection must happen before H removal, or
         // else we might be removing stereogenic H atoms in double
@@ -965,7 +964,8 @@ std::unique_ptr<RWMol> MolFromMol2DataStream(std::istream &inStream,
         // rings in bond stereo detection, and another in
         // sanitization's SSSR symmetrization).
         unsigned int failedOp = 0;
-        MolOps::sanitizeMol(*res, failedOp, MolOps::SanitizeFlags::SANITIZE_CLEANUP);
+        MolOps::sanitizeMol(*res, failedOp,
+                            MolOps::SanitizeFlags::SANITIZE_CLEANUP);
         MolOps::detectBondStereochemistry(*res);
         MolOps::RemoveHsParameters rhp;
         bool sanitize = false;
diff --git a/Code/GraphMol/FileParsers/MolFileParser.cpp b/Code/GraphMol/FileParsers/MolFileParser.cpp
index 9b7ca1bc5..d5642c578 100644
--- a/Code/GraphMol/FileParsers/MolFileParser.cpp
+++ b/Code/GraphMol/FileParsers/MolFileParser.cpp
@@ -367,9 +367,8 @@ void ParseChargeLine(RWMol *mol, const std::string &text, bool firstCall,
   // if this line is specified all the atom other than those specified
   // here should carry a charge of 0; but we should only do this once:
   if (firstCall) {
-    for (ROMol::AtomIterator ai = mol->beginAtoms(); ai != mol->endAtoms();
-         ++ai) {
-      (*ai)->setFormalCharge(0);
+    for (auto at : mol->atoms()) {
+      at->setFormalCharge(0);
     }
   }
 
@@ -408,9 +407,8 @@ void ParseRadicalLine(RWMol *mol, const std::string &text, bool firstCall,
   // if this line is specified all the atom other than those specified
   // here should carry a charge of 0; but we should only do this once:
   if (firstCall) {
-    for (ROMol::AtomIterator ai = mol->beginAtoms(); ai != mol->endAtoms();
-         ++ai) {
-      (*ai)->setFormalCharge(0);
+    for (auto at : mol->atoms()) {
+      at->setFormalCharge(0);
     }
   }
 
diff --git a/Code/GraphMol/FileParsers/MolFileWriter.cpp b/Code/GraphMol/FileParsers/MolFileWriter.cpp
index 90cbd14f1..b1a75ee97 100644
--- a/Code/GraphMol/FileParsers/MolFileWriter.cpp
+++ b/Code/GraphMol/FileParsers/MolFileWriter.cpp
@@ -156,9 +156,7 @@ const std::string GetMolFileChargeInfo(const RWMol &mol) {
   unsigned int nChgs = 0;
   unsigned int nRads = 0;
   unsigned int nMassDiffs = 0;
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    const Atom *atom = *atomIt;
+  for (const auto atom : mol.atoms()) {
     if (atom->getFormalCharge() != 0) {
       ++nChgs;
       chgss << boost::format(" %3d %3d") % (atom->getIdx() + 1) %
@@ -285,12 +283,11 @@ const std::string GetMolFileQueryInfo(
 const std::string GetMolFileRGroupInfo(const RWMol &mol) {
   std::stringstream ss;
   unsigned int nEntries = 0;
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
+  for (const auto atom : mol.atoms()) {
     unsigned int lbl;
-    if ((*atomIt)->getPropIfPresent(common_properties::_MolFileRLabel, lbl)) {
-      ss << " " << std::setw(3) << (*atomIt)->getIdx() + 1 << " "
-         << std::setw(3) << lbl;
+    if (atom->getPropIfPresent(common_properties::_MolFileRLabel, lbl)) {
+      ss << " " << std::setw(3) << atom->getIdx() + 1 << " " << std::setw(3)
+         << lbl;
       ++nEntries;
     }
   }
@@ -303,12 +300,11 @@ const std::string GetMolFileRGroupInfo(const RWMol &mol) {
 
 const std::string GetMolFileAliasInfo(const RWMol &mol) {
   std::stringstream ss;
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
+  for (const auto atom : mol.atoms()) {
     std::string lbl;
-    if ((*atomIt)->getPropIfPresent(common_properties::molFileAlias, lbl)) {
+    if (atom->getPropIfPresent(common_properties::molFileAlias, lbl)) {
       if (!lbl.empty()) {
-        ss << "A  " << std::setw(3) << (*atomIt)->getIdx() + 1 << "\n"
+        ss << "A  " << std::setw(3) << atom->getIdx() + 1 << "\n"
            << lbl << "\n";
       }
     }
@@ -332,19 +328,18 @@ const std::string GetMolFileZBOInfo(const RWMol &mol) {
   std::stringstream ss;
   unsigned int nEntries = 0;
   boost::dynamic_bitset<> atomsAffected(mol.getNumAtoms(), 0);
-  for (ROMol::ConstBondIterator bondIt = mol.beginBonds();
-       bondIt != mol.endBonds(); ++bondIt) {
-    if ((*bondIt)->getBondType() == Bond::ZERO) {
+  for (const auto bond : mol.bonds()) {
+    if (bond->getBondType() == Bond::ZERO) {
       ++nEntries;
-      ss << " " << std::setw(3) << (*bondIt)->getIdx() + 1 << " "
-         << std::setw(3) << 0;
+      ss << " " << std::setw(3) << bond->getIdx() + 1 << " " << std::setw(3)
+         << 0;
       if (nEntries == 8) {
         res << "M  ZBO" << std::setw(3) << nEntries << ss.str() << "\n";
         nEntries = 0;
         ss.str("");
       }
-      atomsAffected[(*bondIt)->getBeginAtomIdx()] = 1;
-      atomsAffected[(*bondIt)->getEndAtomIdx()] = 1;
+      atomsAffected[bond->getBeginAtomIdx()] = 1;
+      atomsAffected[bond->getEndAtomIdx()] = 1;
     }
   }
   if (nEntries) {
@@ -1190,9 +1185,8 @@ std::string getV3000CTAB(const ROMol &tmol,
 
   boost::dynamic_bitset<> queryListAtoms(tmol.getNumAtoms());
   res += "M  V30 BEGIN ATOM\n";
-  for (ROMol::ConstAtomIterator atomIt = tmol.beginAtoms();
-       atomIt != tmol.endAtoms(); ++atomIt) {
-    res += GetV3000MolFileAtomLine(*atomIt, conf, queryListAtoms, precision);
+  for (const auto atom : tmol.atoms()) {
+    res += GetV3000MolFileAtomLine(atom, conf, queryListAtoms, precision);
     res += "\n";
   }
   res += "M  V30 END ATOM\n";
@@ -1378,9 +1372,8 @@ std::string outputMolToMolBlock(const RWMol &tmol, int confId,
   boost::dynamic_bitset<> queryListAtoms(tmol.getNumAtoms());
   if (!isV3000) {
     // V2000 output.
-    for (ROMol::ConstAtomIterator atomIt = tmol.beginAtoms();
-         atomIt != tmol.endAtoms(); ++atomIt) {
-      res += GetMolFileAtomLine(*atomIt, conf, queryListAtoms);
+    for (const auto atom : tmol.atoms()) {
+      res += GetMolFileAtomLine(atom, conf, queryListAtoms);
       res += "\n";
     }
 
diff --git a/Code/GraphMol/FileParsers/PDBParser.cpp b/Code/GraphMol/FileParsers/PDBParser.cpp
index 8cdf570ef..e8a5dedcb 100644
--- a/Code/GraphMol/FileParsers/PDBParser.cpp
+++ b/Code/GraphMol/FileParsers/PDBParser.cpp
@@ -510,9 +510,7 @@ bool StandardPDBChiralAtom(const char *resnam, const char *atmnam) {
 }
 
 void StandardPDBResidueChirality(RWMol *mol) {
-  for (ROMol::AtomIterator atomIt = mol->beginAtoms();
-       atomIt != mol->endAtoms(); ++atomIt) {
-    Atom *atom = *atomIt;
+  for (auto atom : mol->atoms()) {
     if (atom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
       auto *info = (AtomPDBResidueInfo *)atom->getMonomerInfo();
       if (info && info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE &&
@@ -529,10 +527,7 @@ void StandardPDBResidueChirality(RWMol *mol) {
 }
 
 void BasicPDBCleanup(RWMol &mol) {
-  ROMol::VERTEX_ITER atBegin, atEnd;
-  boost::tie(atBegin, atEnd) = mol.getVertices();
-  while (atBegin != atEnd) {
-    Atom *atom = mol[*atBegin];
+  for (auto atom : mol.atoms()) {
     atom->calcExplicitValence(false);
 
     // correct four-valent neutral N -> N+
@@ -541,7 +536,6 @@ void BasicPDBCleanup(RWMol &mol) {
         atom->getValence(Atom::ValenceType::EXPLICIT) == 4) {
       atom->setFormalCharge(1);
     }
-    ++atBegin;
   }
 }
 
diff --git a/Code/GraphMol/FileParsers/PDBWriter.cpp b/Code/GraphMol/FileParsers/PDBWriter.cpp
index 18cc70bef..b32ed2014 100644
--- a/Code/GraphMol/FileParsers/PDBWriter.cpp
+++ b/Code/GraphMol/FileParsers/PDBWriter.cpp
@@ -231,9 +231,8 @@ std::string MolToPDBBody(const ROMol &mol, const Conformer *conf,
   std::string res;
   std::string last;
   std::map<unsigned int, unsigned int> elem;
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    last = GetPDBAtomLine(*atomIt, conf, elem);
+  for (auto atom : mol.atoms()) {
+    last = GetPDBAtomLine(atom, conf, elem);
     res += last;
     res += '\n';
     atm_count++;
@@ -256,9 +255,8 @@ std::string MolToPDBBody(const ROMol &mol, const Conformer *conf,
   bool both = (flavor & 4) != 0;
   bool mult = (flavor & 8) == 0;
   if (all || mult) {
-    for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-         atomIt != mol.endAtoms(); ++atomIt) {
-      res += GetPDBBondLines(*atomIt, all, both, mult, conect_count);
+    for (const auto atom : mol.atoms()) {
+      res += GetPDBBondLines(atom, all, both, mult, conect_count);
     }
   }
   return res;
diff --git a/Code/GraphMol/FileParsers/ProximityBonds.cpp b/Code/GraphMol/FileParsers/ProximityBonds.cpp
index 30b8ca378..dd6a0e423 100644
--- a/Code/GraphMol/FileParsers/ProximityBonds.cpp
+++ b/Code/GraphMol/FileParsers/ProximityBonds.cpp
@@ -534,9 +534,7 @@ static bool StandardPDBDoubleBond(RWMol *mol, Atom *beg, Atom *end) {
 }
 
 void StandardPDBResidueBondOrders(RWMol *mol) {
-  RWMol::BondIterator bondIt;
-  for (bondIt = mol->beginBonds(); bondIt != mol->endBonds(); ++bondIt) {
-    Bond *bond = *bondIt;
+  for (const auto bond : mol->bonds()) {
     if (bond->getBondType() == Bond::SINGLE) {
       Atom *beg = bond->getBeginAtom();
       Atom *end = bond->getEndAtom();
diff --git a/Code/GraphMol/FileParsers/SequenceWriters.cpp b/Code/GraphMol/FileParsers/SequenceWriters.cpp
index 09f8d6952..15ca9e93e 100644
--- a/Code/GraphMol/FileParsers/SequenceWriters.cpp
+++ b/Code/GraphMol/FileParsers/SequenceWriters.cpp
@@ -185,9 +185,7 @@ std::string MolToSequence(const ROMol &mol) {
   std::string result;
   std::string chain;
 
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    const Atom *atom = *atomIt;
+  for (const auto atom : mol.atoms()) {
     auto *info = (AtomPDBResidueInfo *)(atom->getMonomerInfo());
     if (info && info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE) {
       if (info->getName() == " CA ") {
@@ -550,9 +548,7 @@ std::string MolToHELM(const ROMol &mol) {
   int id = 1;
 
   /* First pass: Monomers */
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    const Atom *atom = *atomIt;
+  for (const auto atom : mol.atoms()) {
     auto *info = (AtomPDBResidueInfo *)(atom->getMonomerInfo());
     // We can only write HELM if all atoms have PDB residue information
     if (!info || info->getMonomerType() != AtomMonomerInfo::PDBRESIDUE) {
@@ -706,9 +702,7 @@ std::string MolToHELM(const ROMol &mol) {
   result += "}$";
 
   first = true;
-  for (ROMol::ConstBondIterator bondIt = mol.beginBonds();
-       bondIt != mol.endBonds(); ++bondIt) {
-    const Bond *bond = *bondIt;
+  for (const auto bond : mol.bonds()) {
     Atom *beg = bond->getBeginAtom();
     Atom *end = bond->getEndAtom();
     if (!beg || !end) {
diff --git a/Code/GraphMol/Fingerprints/AtomPairGenerator.cpp b/Code/GraphMol/Fingerprints/AtomPairGenerator.cpp
index 2951782cc..50ad6835a 100644
--- a/Code/GraphMol/Fingerprints/AtomPairGenerator.cpp
+++ b/Code/GraphMol/Fingerprints/AtomPairGenerator.cpp
@@ -32,10 +32,9 @@ std::vector<std::uint32_t> *AtomPairAtomInvGenerator::getAtomInvariants(
     const ROMol &mol) const {
   auto *atomInvariants = new std::vector<std::uint32_t>(mol.getNumAtoms());
 
-  for (ROMol::ConstAtomIterator atomItI = mol.beginAtoms();
-       atomItI != mol.endAtoms(); ++atomItI) {
-    (*atomInvariants)[(*atomItI)->getIdx()] =
-        getAtomCode(*atomItI, 0, df_includeChirality) -
+  for (const auto atom : mol.atoms()) {
+    (*atomInvariants)[atom->getIdx()] =
+        getAtomCode(atom, 0, df_includeChirality) -
         (df_topologicalTorsionCorrection ? 2 : 0);
   }
 
@@ -200,17 +199,13 @@ AtomPairEnvGenerator<OutputType>::getEnvironments(
     distanceMatrix = MolOps::get3DDistanceMat(mol, confId);
   }
 
-  for (ROMol::ConstAtomIterator atomItI = mol.beginAtoms();
-       atomItI != mol.endAtoms(); ++atomItI) {
-    unsigned int i = (*atomItI)->getIdx();
+  for (unsigned int i = 0; i < atomCount; ++i) {
     if (ignoreAtoms && std::find(ignoreAtoms->begin(), ignoreAtoms->end(), i) !=
                            ignoreAtoms->end()) {
       continue;
     }
 
-    for (ROMol::ConstAtomIterator atomItJ = atomItI + 1;
-         atomItJ != mol.endAtoms(); ++atomItJ) {
-      unsigned int j = (*atomItJ)->getIdx();
+    for (unsigned int j = i + 1; j < atomCount; ++j) {
       if (ignoreAtoms && std::find(ignoreAtoms->begin(), ignoreAtoms->end(),
                                    j) != ignoreAtoms->end()) {
         continue;
diff --git a/Code/GraphMol/Fingerprints/AtomPairs.cpp b/Code/GraphMol/Fingerprints/AtomPairs.cpp
index fa453c51b..6c0ea4415 100644
--- a/Code/GraphMol/Fingerprints/AtomPairs.cpp
+++ b/Code/GraphMol/Fingerprints/AtomPairs.cpp
@@ -183,15 +183,14 @@ SparseIntVect<boost::int64_t> *getTopologicalTorsionFingerprint(
 
   std::vector<std::uint32_t> atomCodes;
   atomCodes.reserve(lmol->getNumAtoms());
-  for (ROMol::ConstAtomIterator atomItI = lmol->beginAtoms();
-       atomItI != lmol->endAtoms(); ++atomItI) {
+  for (const auto atom : lmol->atoms()) {
     if (!atomInvariants) {
-      atomCodes.push_back(getAtomCode(*atomItI, 0, includeChirality));
+      atomCodes.push_back(getAtomCode(atom, 0, includeChirality));
     } else {
       // need to add to the atomCode here because we subtract off up to 2 below
       // as part of the branch correction
       atomCodes.push_back(
-          (*atomInvariants)[(*atomItI)->getIdx()] % ((1 << codeSize) - 1) + 2);
+          (*atomInvariants)[atom->getIdx()] % ((1 << codeSize) - 1) + 2);
     }
   }
 
@@ -204,7 +203,7 @@ SparseIntVect<boost::int64_t> *getTopologicalTorsionFingerprint(
   }
   boost::dynamic_bitset<> *ignoreAtomsBV = nullptr;
   if (ignoreAtoms) {
-    ignoreAtomsBV = new boost::dynamic_bitset<>(mol.getNumAtoms());
+    ignoreAtomsBV = new boost::dynamic_bitset<>(lmol->getNumAtoms());
     for (auto fAt : *ignoreAtoms) {
       ignoreAtomsBV->set(fAt);
     }
diff --git a/Code/GraphMol/Fingerprints/FingerprintUtil.cpp b/Code/GraphMol/Fingerprints/FingerprintUtil.cpp
index 24e3b3c7b..6370f1f57 100644
--- a/Code/GraphMol/Fingerprints/FingerprintUtil.cpp
+++ b/Code/GraphMol/Fingerprints/FingerprintUtil.cpp
@@ -272,10 +272,9 @@ void buildDefaultRDKitFingerprintAtomInvariants(
     const ROMol &mol, std::vector<std::uint32_t> &lAtomInvariants) {
   lAtomInvariants.clear();
   lAtomInvariants.reserve(mol.getNumAtoms());
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    unsigned int aHash = ((*atomIt)->getAtomicNum() % 128) << 1 |
-                         static_cast<unsigned int>((*atomIt)->getIsAromatic());
+  for (const auto atom : mol.atoms()) {
+    unsigned int aHash = (atom->getAtomicNum() % 128) << 1 |
+                         static_cast<unsigned int>(atom->getIsAromatic());
     lAtomInvariants.push_back(aHash);
   }
 }
diff --git a/Code/GraphMol/Fingerprints/MACCS.cpp b/Code/GraphMol/Fingerprints/MACCS.cpp
index 43cdc75b0..389afdf65 100644
--- a/Code/GraphMol/Fingerprints/MACCS.cpp
+++ b/Code/GraphMol/Fingerprints/MACCS.cpp
@@ -327,8 +327,8 @@ void GenerateFP(const RDKit::ROMol &mol, ExplicitBitVect &fp) {
   RDKit::MatchVectType match;
   unsigned int count;
 
-  for (atom = mol.beginAtoms(); atom != mol.endAtoms(); ++atom) {
-    switch ((*atom)->getAtomicNum()) {
+  for (const auto atom : mol.atoms()) {
+    switch (atom->getAtomicNum()) {
       case 3:
       case 11:
       case 19:
diff --git a/Code/GraphMol/Fingerprints/RDKitFPGenerator.cpp b/Code/GraphMol/Fingerprints/RDKitFPGenerator.cpp
index 0fd65d0a6..8d0a27ee3 100644
--- a/Code/GraphMol/Fingerprints/RDKitFPGenerator.cpp
+++ b/Code/GraphMol/Fingerprints/RDKitFPGenerator.cpp
@@ -45,10 +45,9 @@ std::vector<std::uint32_t> *RDKitFPAtomInvGenerator::getAtomInvariants(
     const ROMol &mol) const {
   auto *result = new std::vector<std::uint32_t>();
   result->reserve(mol.getNumAtoms());
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    unsigned int aHash = ((*atomIt)->getAtomicNum() % 128) << 1 |
-                         static_cast<unsigned int>((*atomIt)->getIsAromatic());
+  for (const auto atom : mol.atoms()) {
+    unsigned int aHash = (atom->getAtomicNum() % 128) << 1 |
+                         static_cast<unsigned int>(atom->getIsAromatic());
     result->push_back(aHash);
   }
   return result;
diff --git a/Code/GraphMol/ForceFieldHelpers/MMFF/Builder.cpp b/Code/GraphMol/ForceFieldHelpers/MMFF/Builder.cpp
index 4ac87b8fd..9b1af24a3 100644
--- a/Code/GraphMol/ForceFieldHelpers/MMFF/Builder.cpp
+++ b/Code/GraphMol/ForceFieldHelpers/MMFF/Builder.cpp
@@ -52,10 +52,9 @@ void addBonds(const ROMol &mol, MMFFMolProperties *mmffMolProperties,
 
   auto contrib = std::make_unique<BondStretchContrib>(field);
   bool hasContrib = false;
-  for (ROMol::ConstBondIterator bi = mol.beginBonds(); bi != mol.endBonds();
-       ++bi) {
-    unsigned int idx1 = (*bi)->getBeginAtomIdx();
-    unsigned int idx2 = (*bi)->getEndAtomIdx();
+  for (const auto bond : mol.bonds()) {
+    unsigned int idx1 = bond->getBeginAtomIdx();
+    unsigned int idx2 = bond->getEndAtomIdx();
     unsigned int bondType;
     MMFFBond mmffBondParams;
     if (mmffMolProperties->getMMFFBondStretchParams(mol, idx1, idx2, bondType,
@@ -65,8 +64,8 @@ void addBonds(const ROMol &mol, MMFFMolProperties *mmffMolProperties,
       if (mmffMolProperties->getMMFFVerbosity()) {
         unsigned int iAtomType = mmffMolProperties->getMMFFAtomType(idx1);
         unsigned int jAtomType = mmffMolProperties->getMMFFAtomType(idx2);
-        const Atom *iAtom = (*bi)->getBeginAtom();
-        const Atom *jAtom = (*bi)->getEndAtom();
+        const Atom *iAtom = bond->getBeginAtom();
+        const Atom *jAtom = bond->getEndAtom();
         const double dist = field->distance(idx1, idx2);
         const double bondStretchEnergy = MMFF::Utils::calcBondStretchEnergy(
             mmffBondParams.r0, mmffBondParams.kb, dist);
diff --git a/Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp b/Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp
index f79c9fcab..1750059b4 100644
--- a/Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp
+++ b/Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp
@@ -34,18 +34,17 @@ void addBonds(const ROMol &mol, const AtomicParamVect &params,
   PRECONDITION(mol.getNumAtoms() == params.size(), "bad parameters");
   PRECONDITION(field, "bad forcefield");
 
-  for (ROMol::ConstBondIterator bi = mol.beginBonds(); bi != mol.endBonds();
-       bi++) {
-    int idx1 = (*bi)->getBeginAtomIdx();
-    int idx2 = (*bi)->getEndAtomIdx();
+  for (const auto bond : mol.bonds()) {
+    int idx1 = bond->getBeginAtomIdx();
+    int idx2 = bond->getEndAtomIdx();
 
     // FIX: recognize amide bonds here.
 
     if (params[idx1] && params[idx2]) {
       BondStretchContrib *contrib;
-      contrib = new BondStretchContrib(field, idx1, idx2,
-                                       (*bi)->getBondTypeAsDouble(),
-                                       params[idx1], params[idx2]);
+      contrib =
+          new BondStretchContrib(field, idx1, idx2, bond->getBondTypeAsDouble(),
+                                 params[idx1], params[idx2]);
       field->contribs().push_back(ForceFields::ContribPtr(contrib));
     }
   }
@@ -96,19 +95,21 @@ boost::shared_array<std::uint8_t> buildNeighborMatrix(const ROMol &mol) {
   unsigned nTwoBitCells = (nAtoms * (nAtoms + 1) - 1) / 8 + 1;
   boost::shared_array<std::uint8_t> res(new std::uint8_t[nTwoBitCells]);
   std::memset(res.get(), RELATION_1_X_INIT, nTwoBitCells);
-  for (ROMol::ConstBondIterator bondi = mol.beginBonds();
-       bondi != mol.endBonds(); ++bondi) {
-    setTwoBitCell(res,
-                  twoBitCellPos(nAtoms, (*bondi)->getBeginAtomIdx(),
-                                (*bondi)->getEndAtomIdx()),
-                  RELATION_1_2);
-    unsigned int bondiBeginAtomIdx = (*bondi)->getBeginAtomIdx();
-    unsigned int bondiEndAtomIdx = (*bondi)->getEndAtomIdx();
-    for (ROMol::ConstBondIterator bondj = bondi; ++bondj != mol.endBonds();) {
+  for (const auto bondi : mol.bonds()) {
+    setTwoBitCell(
+        res,
+        twoBitCellPos(nAtoms, bondi->getBeginAtomIdx(), bondi->getEndAtomIdx()),
+        RELATION_1_2);
+    unsigned int bondiBeginAtomIdx = bondi->getBeginAtomIdx();
+    unsigned int bondiEndAtomIdx = bondi->getEndAtomIdx();
+    for (const auto bondj : mol.bonds()) {
+      if (bondj == bondi) {
+        continue;
+      }
       int idx1 = -1;
       int idx3 = -1;
-      unsigned int bondjBeginAtomIdx = (*bondj)->getBeginAtomIdx();
-      unsigned int bondjEndAtomIdx = (*bondj)->getEndAtomIdx();
+      unsigned int bondjBeginAtomIdx = bondj->getBeginAtomIdx();
+      unsigned int bondjEndAtomIdx = bondj->getEndAtomIdx();
       if (bondiBeginAtomIdx == bondjBeginAtomIdx) {
         idx1 = bondiEndAtomIdx;
         idx3 = bondjEndAtomIdx;
diff --git a/Code/GraphMol/FragCatalog/FragCatalogEntry.cpp b/Code/GraphMol/FragCatalog/FragCatalogEntry.cpp
index ba814926f..b1e881275 100644
--- a/Code/GraphMol/FragCatalog/FragCatalogEntry.cpp
+++ b/Code/GraphMol/FragCatalog/FragCatalogEntry.cpp
@@ -154,9 +154,7 @@ Subgraphs::DiscrimTuple FragCatalogEntry::getDiscrims() const {
     // the atoms
     std::vector<std::uint32_t> funcGpInvars;
     gboost::hash<INT_VECT> vectHasher;
-    for (ROMol::AtomIterator atomIt = dp_mol->beginAtoms();
-         atomIt != dp_mol->endAtoms(); ++atomIt) {
-      unsigned int aid = (*atomIt)->getIdx();
+    for (unsigned int aid = 0; aid < dp_mol->getNumAtoms(); ++aid) {
       std::uint32_t invar = 0;
       auto mapPos = d_aToFmap.find(aid);
       if (mapPos != d_aToFmap.end()) {
diff --git a/Code/GraphMol/MarvinParse/MarvinParser.cpp b/Code/GraphMol/MarvinParse/MarvinParser.cpp
index 9767fbae2..19c3882ae 100644
--- a/Code/GraphMol/MarvinParse/MarvinParser.cpp
+++ b/Code/GraphMol/MarvinParse/MarvinParser.cpp
@@ -600,10 +600,9 @@ class MarvinCMLReader {
       mol->clearAllBondBookmarks();
 
       // calculate explicit valence on each atom:
-      for (RWMol::AtomIterator atomIt = mol->beginAtoms();
-           atomIt != mol->endAtoms(); ++atomIt) {
-        (*atomIt)->calcExplicitValence(false);
-        (*atomIt)->calcImplicitValence(false);
+      for (auto atom : mol->atoms()) {
+        atom->calcExplicitValence(false);
+        atom->calcImplicitValence(false);
       }
 
       // update the chirality and stereo-chemistry
diff --git a/Code/GraphMol/Matrices.cpp b/Code/GraphMol/Matrices.cpp
index 72c0e3e4f..00d22094d 100644
--- a/Code/GraphMol/Matrices.cpp
+++ b/Code/GraphMol/Matrices.cpp
@@ -309,23 +309,22 @@ double *getAdjacencyMatrix(const ROMol &mol, bool useBO, int emptyVal,
   auto *res = new double[nAts * nAts];
   memset(static_cast<void *>(res), emptyVal, nAts * nAts * sizeof(double));
 
-  for (ROMol::ConstBondIterator bondIt = mol.beginBonds();
-       bondIt != mol.endBonds(); bondIt++) {
-    if (bondsToUse && !(*bondsToUse)[(*bondIt)->getIdx()]) {
+  for (const auto bond : mol.bonds()) {
+    if (bondsToUse && !(*bondsToUse)[bond->getIdx()]) {
       continue;
     }
     if (!useBO) {
-      int beg = (*bondIt)->getBeginAtomIdx();
-      int end = (*bondIt)->getEndAtomIdx();
+      int beg = bond->getBeginAtomIdx();
+      int end = bond->getEndAtomIdx();
       res[beg * nAts + end] = 1;
       res[end * nAts + beg] = 1;
     } else {
-      int begIdx = (*bondIt)->getBeginAtomIdx();
-      int endIdx = (*bondIt)->getEndAtomIdx();
+      int begIdx = bond->getBeginAtomIdx();
+      int endIdx = bond->getEndAtomIdx();
       Atom const *beg = mol.getAtomWithIdx(begIdx);
       Atom const *end = mol.getAtomWithIdx(endIdx);
-      res[begIdx * nAts + endIdx] = (*bondIt)->getValenceContrib(beg);
-      res[endIdx * nAts + begIdx] = (*bondIt)->getValenceContrib(end);
+      res[begIdx * nAts + endIdx] = bond->getValenceContrib(beg);
+      res[endIdx * nAts + begIdx] = bond->getValenceContrib(end);
     }
   }
   sptr.reset(res);
diff --git a/Code/GraphMol/MolOps.cpp b/Code/GraphMol/MolOps.cpp
index 5836d0e3c..1f8f28926 100644
--- a/Code/GraphMol/MolOps.cpp
+++ b/Code/GraphMol/MolOps.cpp
@@ -1011,12 +1011,11 @@ template RDKIT_GRAPHMOL_EXPORT unsigned int getMolFragsWithQuery(
     bool sanitizeFrags, const std::vector<unsigned int> *, bool);
 
 int getFormalCharge(const ROMol &mol) {
-  int accum = 0;
-  for (ROMol::ConstAtomIterator atomIt = mol.beginAtoms();
-       atomIt != mol.endAtoms(); ++atomIt) {
-    accum += (*atomIt)->getFormalCharge();
-  }
-  return accum;
+  auto res = std::accumulate(mol.atoms().begin(), mol.atoms().end(), 0,
+                             [](int accum, const auto atom) {
+                               return accum + atom->getFormalCharge();
+                             });
+  return res;
 };
 
 unsigned getNumAtomsWithDistinctProperty(const ROMol &mol,
diff --git a/Code/GraphMol/MolPickler.cpp b/Code/GraphMol/MolPickler.cpp
index 00602474d..efcbcd731 100644
--- a/Code/GraphMol/MolPickler.cpp
+++ b/Code/GraphMol/MolPickler.cpp
@@ -843,16 +843,16 @@ void pickleAtomPDBResidueInfo(std::ostream &ss,
   }
 }
 
-AtomMonomerInfo* unpickleAtomPDBResidueInfo(std::istream &ss, int version) {
+AtomMonomerInfo *unpickleAtomPDBResidueInfo(std::istream &ss, int version) {
   std::string nm;
   streamRead(ss, nm, version);
   unsigned int typ;
   streamRead(ss, typ, version);
 
-  // As of version 16.3, some member fields from AtomPDBResidueInfo were moved to
-  // the AtomMonomerInfo base class. Pickles made from 16.3+ will use the new tags
-  // and unpickleAtomMonomerInfo to unpickle, but we need to continue to support old
-  // pickles that don't use the new tags.
+  // As of version 16.3, some member fields from AtomPDBResidueInfo were moved
+  // to the AtomMonomerInfo base class. Pickles made from 16.3+ will use the new
+  // tags and unpickleAtomMonomerInfo to unpickle, but we need to continue to
+  // support old pickles that don't use the new tags.
   auto type = static_cast<AtomMonomerInfo::AtomMonomerType>(typ);
   if (type != AtomMonomerInfo::AtomMonomerType::PDBRESIDUE) {
     auto info = new AtomMonomerInfo(type, nm);
@@ -945,8 +945,7 @@ void pickleAtomMonomerInfo(std::ostream &ss, const AtomMonomerInfo *info) {
                 info->getResidueNumber());
   }
   if (!info->getChainId().empty()) {
-    streamWrite(ss, MolPickler::ATOM_MONOMER_INFO_CHAINID,
-                info->getChainId());
+    streamWrite(ss, MolPickler::ATOM_MONOMER_INFO_CHAINID, info->getChainId());
   }
   if (!info->getMonomerClass().empty()) {
     streamWrite(ss, MolPickler::ATOM_MONOMER_INFO_MONOMERCLASS,
@@ -959,7 +958,8 @@ AtomMonomerInfo *unpickleAtomMonomerInfo(std::istream &ss, int version) {
   std::uint8_t typ;
   streamRead(ss, nm, version);
   streamRead(ss, typ, version);
-  auto info = new AtomMonomerInfo(static_cast<RDKit::AtomMonomerInfo::AtomMonomerType>(typ), nm);
+  auto info = new AtomMonomerInfo(
+      static_cast<RDKit::AtomMonomerInfo::AtomMonomerType>(typ), nm);
 
   std::string residueName = "";
   int residueNumber = 0;
@@ -990,8 +990,8 @@ AtomMonomerInfo *unpickleAtomMonomerInfo(std::istream &ss, int version) {
       default:
         // None of the ATOM_PDB_RESIDUE_XXX tags should appear here
         throw MolPicklerException(
-            "unrecognized tag while parsing atom monomer info " + std::to_string(
-                static_cast<int>(tag)));
+            "unrecognized tag while parsing atom monomer info " +
+            std::to_string(static_cast<int>(tag)));
     }
   }
   return info;
@@ -1205,9 +1205,9 @@ void MolPickler::_pickle(const ROMol *mol, std::ostream &ss,
   streamWrite(ss, BEGINATOM);
   ROMol::ConstAtomIterator atIt;
   int nWritten = 0;
-  for (atIt = mol->beginAtoms(); atIt != mol->endAtoms(); ++atIt) {
-    _pickleAtom<T>(ss, *atIt);
-    atomIdxMap[(*atIt)->getIdx()] = nWritten;
+  for (auto atom : mol->atoms()) {
+    _pickleAtom<T>(ss, atom);
+    atomIdxMap[atom->getIdx()] = nWritten;
     nWritten++;
   }
 
@@ -1217,10 +1217,9 @@ void MolPickler::_pickle(const ROMol *mol, std::ostream &ss,
   //
   // -------------------
   streamWrite(ss, BEGINBOND);
-  for (unsigned int i = 0; i < mol->getNumBonds(); i++) {
-    auto bond = mol->getBondWithIdx(i);
+  for (auto bond : mol->bonds()) {
     _pickleBond<T>(ss, bond, atomIdxMap);
-    bondIdxMap[bond->getIdx()] = i;
+    bondIdxMap[bond->getIdx()] = bond->getIdx();
   }
 
   // -------------------
@@ -1786,9 +1785,11 @@ void MolPickler::_pickleAtom(std::ostream &ss, const Atom *atom) {
                 atom->getProp<std::string>(common_properties::dummyLabel));
   }
   if (atom->getMonomerInfo()) {
-    if (atom->getMonomerInfo()->getMonomerType() == AtomMonomerInfo::PDBRESIDUE) {
+    if (atom->getMonomerInfo()->getMonomerType() ==
+        AtomMonomerInfo::PDBRESIDUE) {
       streamWrite(ss, BEGIN_PDB_RESIDUE);
-      pickleAtomPDBResidueInfo(ss, static_cast<const AtomPDBResidueInfo*>(atom->getMonomerInfo()));
+      pickleAtomPDBResidueInfo(
+          ss, static_cast<const AtomPDBResidueInfo *>(atom->getMonomerInfo()));
       streamWrite(ss, END_PDB_RESIDUE);
     } else {
       streamWrite(ss, BEGIN_ATOM_MONOMER_INFO);
@@ -2600,9 +2601,7 @@ void MolPickler::_pickleV1(const ROMol *mol, std::ostream &ss) {
   if (mol->getNumConformers() > 0) {
     conf = &(mol->getConformer());
   }
-  for (atIt = mol->beginAtoms(); atIt != mol->endAtoms(); ++atIt) {
-    const Atom *atom = *atIt;
-
+  for (const auto atom : mol->atoms()) {
     streamWrite(ss, BEGINATOM);
     streamWrite(ss, ATOM_NUMBER, atom->getAtomicNum());
 
@@ -2633,9 +2632,7 @@ void MolPickler::_pickleV1(const ROMol *mol, std::ostream &ss) {
     streamWrite(ss, ENDATOM);
   }
 
-  ROMol::ConstBondIterator bondIt;
-  for (bondIt = mol->beginBonds(); bondIt != mol->endBonds(); ++bondIt) {
-    const Bond *bond = *bondIt;
+  for (const auto bond : mol->bonds()) {
     streamWrite(ss, BEGINBOND);
     streamWrite(ss, BOND_INDEX, bond->getIdx());
     streamWrite(ss, BOND_BEGATOMIDX, bond->getBeginAtomIdx());
diff --git a/Code/GraphMol/MolTransforms/MolTransforms.cpp b/Code/GraphMol/MolTransforms/MolTransforms.cpp
index 18ddb12de..296874765 100644
--- a/Code/GraphMol/MolTransforms/MolTransforms.cpp
+++ b/Code/GraphMol/MolTransforms/MolTransforms.cpp
@@ -37,9 +37,8 @@ void transformAtom(Atom *atom, RDGeom::Transform3D &tform) {
 void transformMolsAtoms(ROMol *mol, RDGeom::Transform3D &tform) {
   PRECONDITION(mol, "no molecule");
 
-  ROMol::AtomIterator atomIt;
-  for (atomIt = mol->beginAtoms(); atomIt != mol->endAtoms(); atomIt++) {
-    transformAtom(*atomIt, tform);
+  for (auto atom : mol->atoms()) {
+    transformAtom(atom, tform);
   }
 }
 
@@ -74,16 +73,15 @@ void computeCovarianceTerms(const Conformer &conf,
   xx = xy = xz = yy = yz = zz = 0.0;
   const ROMol &mol = conf.getOwningMol();
   double wSum = 0.0;
-  for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
-       cai++) {
-    if (((*cai)->getAtomicNum() == 1) && (ignoreHs)) {
+  for (const auto atom : mol.atoms()) {
+    if ((atom->getAtomicNum() == 1) && (ignoreHs)) {
       continue;
     }
-    RDGeom::Point3D loc = conf.getAtomPos((*cai)->getIdx());
+    RDGeom::Point3D loc = conf.getAtomPos(atom->getIdx());
     loc -= center;
     double w = 1.0;
     if (weights) {
-      w = (*weights)[(*cai)->getIdx()];
+      w = (*weights)[atom->getIdx()];
     }
     wSum += w;
     xx += w * loc.x * loc.x;
@@ -130,16 +128,15 @@ void computeInertiaTerms(const Conformer &conf, const RDGeom::Point3D &center,
 
   xx = xy = xz = yy = yz = zz = 0.0;
   const ROMol &mol = conf.getOwningMol();
-  for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
-       cai++) {
-    if (((*cai)->getAtomicNum() == 1) && (ignoreHs)) {
+  for (const auto atom : mol.atoms()) {
+    if ((atom->getAtomicNum() == 1) && (ignoreHs)) {
       continue;
     }
-    RDGeom::Point3D loc = conf.getAtomPos((*cai)->getIdx());
+    RDGeom::Point3D loc = conf.getAtomPos(atom->getIdx());
     loc -= center;
     double w = 1.0;
     if (weights) {
-      w = (*weights)[(*cai)->getIdx()];
+      w = (*weights)[atom->getIdx()];
     }
     xx += w * (loc.y * loc.y + loc.z * loc.z);
     yy += w * (loc.x * loc.x + loc.z * loc.z);
diff --git a/Code/GraphMol/PartialCharges/GasteigerCharges.cpp b/Code/GraphMol/PartialCharges/GasteigerCharges.cpp
index 5e4bc7201..786c3df83 100644
--- a/Code/GraphMol/PartialCharges/GasteigerCharges.cpp
+++ b/Code/GraphMol/PartialCharges/GasteigerCharges.cpp
@@ -18,66 +18,48 @@
 namespace Gasteiger {
 using namespace RDKit;
 /*! \brief split the formal charge across atoms of same type if we have a
- *conjugated system
+ * conjugated system
  *
- *  This function is called before the charge equivalization iteration is
- *started for the
- *  Gasteiger charges. If any of the atom involved in conjugated system have
- *formal charges
- *  set on them, this charges is equally distributed across atoms of the same
- *type in that
- *  conjugated system. So for example the two nitrogens in the benzamidine
- *system start the iteration
- * with equal charges of 0.5
+ * This function is called before the charge equivalization iteration is
+ * started for the Gasteiger charges. If any of the atom involved in conjugated
+ * system have formal charges set on them, this charges is equally distributed
+ * across atoms of the same type in that conjugated system. So for example the
+ * two nitrogens in the benzamidine system start the iteration with equal
+ * charges of 0.5
  */
 void splitChargeConjugated(const ROMol &mol, DOUBLE_VECT &charges) {
-  int aix;
-  int natms = mol.getNumAtoms();
   INT_VECT marker;
-  INT_VECT_CI mci;
-  int aax, yax;
-  double formal;
-  const Atom *at, *aat, *yat;
-  for (aix = 0; aix < natms; aix++) {
-    at = mol.getAtomWithIdx(aix);
-    formal = at->getFormalCharge();
+  for (const auto at : mol.atoms()) {
+    auto aix = at->getIdx();
+    double formal = at->getFormalCharge();
     // std::cout << aix << " formal charges:" << formal << "\n";
     marker.resize(0);
     if ((fabs(formal) > EPS_DOUBLE) && (fabs(charges[aix]) < EPS_DOUBLE)) {
       marker.push_back(aix);
-      ROMol::OEDGE_ITER bnd1, end1, bnd2, end2;
-      boost::tie(bnd1, end1) = mol.getAtomBonds(at);
-      while (bnd1 != end1) {
-        if (mol[*bnd1]->getIsConjugated()) {
-          aax = mol[*bnd1]->getOtherAtomIdx(aix);
-          aat = mol.getAtomWithIdx(aax);
-          boost::tie(bnd2, end2) = mol.getAtomBonds(aat);
-          while (bnd2 != end2) {
-            if ((*bnd1) != (*bnd2)) {
-              if (mol[*bnd2]->getIsConjugated()) {
-                yax = mol[*bnd2]->getOtherAtomIdx(aax);
-                yat = mol.getAtomWithIdx(yax);
+      for (const auto bnd1 : mol.atomBonds(at)) {
+        if (bnd1->getIsConjugated()) {
+          auto aax = bnd1->getOtherAtomIdx(aix);
+          auto aat = mol.getAtomWithIdx(aax);
+          for (const auto bnd2 : mol.atomBonds(aat)) {
+            if (bnd1 != bnd2) {
+              if (bnd2->getIsConjugated()) {
+                auto yax = bnd2->getOtherAtomIdx(aax);
+                auto yat = mol.getAtomWithIdx(yax);
                 if (at->getAtomicNum() == yat->getAtomicNum()) {
                   formal += yat->getFormalCharge();
                   marker.push_back(yax);
                 }
               }
             }
-            bnd2++;
           }
         }
-        bnd1++;
       }
 
-      for (mci = marker.begin(); mci != marker.end(); mci++) {
-        charges[*mci] = (formal / marker.size());
+      for (const auto mci : marker) {
+        charges[mci] = (formal / marker.size());
       }
     }
   }
-  /*
-  for (aix = 0; aix < natms; aix++) {
-    std::cout << "In splitter: " << " charges:" << charges[aix] << "\n";
-    }*/
 }
 }  // end of namespace Gasteiger
 
@@ -104,42 +86,32 @@ void computeGasteigerCharges(const ROMol &mol, std::vector<double> &charges,
                              int nIter, bool throwOnParamFailure) {
   PRECONDITION(charges.size() >= mol.getNumAtoms(), "bad array size");
 
-  PeriodicTable *table = PeriodicTable::getTable();
+  const PeriodicTable *table = PeriodicTable::getTable();
   const GasteigerParams *params = GasteigerParams::getParams();
 
-  double damp = DAMP;
-  int natms = mol.getNumAtoms();
+  const auto natms = mol.getNumAtoms();
   // space for parameters for each atom in the molecule
   std::vector<DOUBLE_VECT> atmPs;
   atmPs.reserve(natms);
 
   std::fill(charges.begin(), charges.end(), 0.0);
 
-  DOUBLE_VECT
-  hChrg;  // total charge on the implicit hydrogen on each heavy atom
-  hChrg.resize(natms, 0.0);
+  DOUBLE_VECT hChrg(
+      natms,
+      0.0);  // total charge for the implicit hydrogens on each heavy atom
 
-  DOUBLE_VECT ionX;
-  ionX.resize(natms, 0.0);
-
-  DOUBLE_VECT energ;
-  energ.resize(natms, 0.0);
-
-  ROMol::ADJ_ITER nbrIdx, endIdx;
+  DOUBLE_VECT ionX(natms, 0.0);
+  DOUBLE_VECT energ(natms, 0.0);
 
   // deal with the conjugated system - distribute the formal charges on atoms of
-  // same type in each
-  // conjugated system
+  // same type in each conjugated system
   Gasteiger::splitChargeConjugated(mol, charges);
 
-  // now read in the parameters
-  ROMol::ConstAtomIterator ai;
-
-  for (ai = mol.beginAtoms(); ai != mol.endAtoms(); ai++) {
-    std::string elem = table->getElementSymbol((*ai)->getAtomicNum());
+  for (const auto atom : mol.atoms()) {
+    const std::string elem = table->getElementSymbol(atom->getAtomicNum());
     std::string mode;
 
-    switch ((*ai)->getHybridization()) {
+    switch (atom->getHybridization()) {
       case Atom::SP3:
         mode = "sp3";
         break;
@@ -150,19 +122,17 @@ void computeGasteigerCharges(const ROMol &mol, std::vector<double> &charges,
         mode = "sp";
         break;
       default:
-        if ((*ai)->getAtomicNum() == 1) {
+        if (atom->getAtomicNum() == 1) {
           // if it is hydrogen
           mode = "*";
-        } else if ((*ai)->getAtomicNum() == 16) {
+        } else if (atom->getAtomicNum() == 16) {
           // we have a sulfur atom with no hybridization information
           // check how many oxygens we have on the sulfur
-          boost::tie(nbrIdx, endIdx) = mol.getAtomNeighbors(*ai);
           int no = 0;
-          while (nbrIdx != endIdx) {
-            if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() == 8) {
+          for (const auto nbrAt : mol.atomNeighbors(atom)) {
+            if (nbrAt->getAtomicNum() == 8) {
               no++;
             }
-            nbrIdx++;
           }
           if (no == 2) {
             mode = "so2";
@@ -175,77 +145,66 @@ void computeGasteigerCharges(const ROMol &mol, std::vector<double> &charges,
         }
     }
 
-    // if we get a unknown mode or element type the
-    // following will will throw an exception
+    // if we get a unknown mode or element type the following will will throw an
+    // exception
     atmPs.push_back(params->getParams(elem, mode, throwOnParamFailure));
 
     // set ionX parameters
     // if Hydrogen treat differently
-    int idx = (*ai)->getIdx();
-    if ((*ai)->getAtomicNum() == 1) {
+    const auto idx = atom->getIdx();
+    if (atom->getAtomicNum() == 1) {
       ionX[idx] = IONXH;
     } else {
       ionX[idx] = atmPs[idx][0] + atmPs[idx][1] + atmPs[idx][2];
     }
   }
 
-  // do the iteration here
-  int itx, aix, sgn, niHs;
-  double enr, dq, dx, qHs, dqH;
-  // parameters for hydrogen atoms (for case where the hydrogen are not in the
+  // parameters for hydrogen atoms for case where the hydrogen are not in the
   // graph (implicit hydrogens)
-  DOUBLE_VECT hParams;
-  hParams = params->getParams("H", "*", throwOnParamFailure);
+  DOUBLE_VECT hParams = params->getParams("H", "*", throwOnParamFailure);
 
-  /*
-  int itmp;
-  for (itmp = 0; itmp < 5; itmp++) {
-    std::cout << " aq:" << charges[itmp] << "\n";
-    }*/
-
-  for (itx = 0; itx < nIter; itx++) {
-    for (aix = 0; aix < natms; aix++) {
-      // enr = p0 + charge*(p1 + p2*charge)
-      enr = atmPs[aix][0] +
-            charges[aix] * (atmPs[aix][1] + atmPs[aix][2] * charges[aix]);
+  double damp = DAMP;
+  for (auto itx = 0u; itx < static_cast<unsigned int>(nIter); itx++) {
+    for (auto aix = 0u; aix < natms; aix++) {
+      auto enr = atmPs[aix][0] +
+                 charges[aix] * (atmPs[aix][1] + atmPs[aix][2] * charges[aix]);
       energ[aix] = enr;
     }
 
-    for (aix = 0; aix < natms; aix++) {
-      dq = 0.0;
-      boost::tie(nbrIdx, endIdx) =
-          mol.getAtomNeighbors(mol.getAtomWithIdx(aix));
-      while (nbrIdx != endIdx) {
-        dx = energ[*nbrIdx] - energ[aix];
+    for (const auto atom : mol.atoms()) {
+      const auto aix = atom->getIdx();
+      double dq = 0.0;
+      for (const auto nbrAt : mol.atomNeighbors(atom)) {
+        const auto nbr = nbrAt->getIdx();
+        const auto dx = energ[nbr] - energ[aix];
+        int sgn;
         if (dx < 0.0) {
           sgn = 0;
         } else {
           sgn = 1;
         }
-        dq += dx / ((sgn * (ionX[aix] - ionX[*nbrIdx])) + ionX[*nbrIdx]);
-        nbrIdx++;
+        dq += dx / ((sgn * (ionX[aix] - ionX[nbr])) + ionX[nbr]);
       }
       // now loop over the implicit hydrogens and get their contributions
       // since hydrogens don't connect to anything else, update their charges at
       // the same time
-      niHs = mol.getAtomWithIdx(aix)->getTotalNumHs();
+      const auto niHs = atom->getTotalNumHs();
       if (niHs > 0) {
-        qHs = hChrg[aix] / niHs;
-        enr = hParams[0] + qHs * (hParams[1] + hParams[2] * qHs);
-        dx = enr - energ[aix];
+        const auto qHs = hChrg[aix] / niHs;
+        const auto enr = hParams[0] + qHs * (hParams[1] + hParams[2] * qHs);
+        const auto dx = enr - energ[aix];
+        int sgn;
         if (dx < 0.0) {
           sgn = 0;
         } else {
           sgn = 1;
         }
 
-        dqH = dx / ((sgn * (ionX[aix] - IONXH)) + IONXH);
-
+        const auto dqH = dx / ((sgn * (ionX[aix] - IONXH)) + IONXH);
         dq += (niHs * dqH);
 
         // adjust the charges on the hydrogens simultaneously (possible because
-        // each of the
-        // hydrogens have no other neighbors)
+        // each of the hydrogens have no other neighbors)
         hChrg[aix] -= (niHs * dqH * damp);
       }
       charges[aix] += (damp * dq);
@@ -254,12 +213,12 @@ void computeGasteigerCharges(const ROMol &mol, std::vector<double> &charges,
     damp *= DAMP_SCALE;
   }
 
-  for (aix = 0; aix < natms; aix++) {
-    mol.getAtomWithIdx(aix)->setProp(common_properties::_GasteigerCharge,
-                                     charges[aix], true);
+  for (auto atom : mol.atoms()) {
+    const auto aix = atom->getIdx();
+    // set the charge on the heavy atom
+    atom->setProp(common_properties::_GasteigerCharge, charges[aix], true);
     // set the implicit hydrogen charges
-    mol.getAtomWithIdx(aix)->setProp(common_properties::_GasteigerHCharge,
-                                     hChrg[aix], true);
+    atom->setProp(common_properties::_GasteigerHCharge, hChrg[aix], true);
   }
 }
 }  // namespace RDKit
diff --git a/Code/GraphMol/PartialCharges/GasteigerParams.h b/Code/GraphMol/PartialCharges/GasteigerParams.h
index 42b6cac2c..a3bf1988f 100644
--- a/Code/GraphMol/PartialCharges/GasteigerParams.h
+++ b/Code/GraphMol/PartialCharges/GasteigerParams.h
@@ -22,10 +22,10 @@ extern std::string additionalParamData;
 
 // this is a constant used during the iteration procedure for the hydrogen atoms
 // for the remaining atoms it is computed on the fly
-const double IONXH = 20.02;
+constexpr double IONXH = 20.02;
 
-const double DAMP_SCALE = 0.5;
-const double DAMP = 0.5;
+constexpr double DAMP_SCALE = 0.5;
+constexpr double DAMP = 0.5;
 
 class RDKIT_PARTIALCHARGES_EXPORT GasteigerParams {
   /* \brief Container for all the partial charge parameters
diff --git a/Code/GraphMol/PartialCharges/Wrap/testPartialCharges.py b/Code/GraphMol/PartialCharges/Wrap/testPartialCharges.py
index f3894fcee..f4b5316e6 100644
--- a/Code/GraphMol/PartialCharges/Wrap/testPartialCharges.py
+++ b/Code/GraphMol/PartialCharges/Wrap/testPartialCharges.py
@@ -44,7 +44,7 @@ class TestCase(unittest.TestCase):
 
     i = 0
     for line in lines:
-      self.assertTrue(line.strip() == olst[i])
+      self.assertEqual(line.strip(), olst[i])
       i += 1
 
   def test1PPDataset(self):
@@ -92,7 +92,7 @@ class TestCase(unittest.TestCase):
     for i in range(m1.GetNumAtoms()):
       c1 = float(m1.GetAtomWithIdx(i).GetProp('_GasteigerCharge'))
       c2 = float(m2.GetAtomWithIdx(i).GetProp('_GasteigerCharge'))
-      self.assertTrue(feq(c1, c2, 1e-4))
+      self.assertAlmostEqual(c1, c2, 4)
 
   def test3Params(self):
     """ tests handling of Issue187 """
diff --git a/Code/GraphMol/RGroupDecomposition/RGroupDecompData.cpp b/Code/GraphMol/RGroupDecomposition/RGroupDecompData.cpp
index 89d14cfef..07cb27db8 100644
--- a/Code/GraphMol/RGroupDecomposition/RGroupDecompData.cpp
+++ b/Code/GraphMol/RGroupDecomposition/RGroupDecompData.cpp
@@ -460,9 +460,7 @@ void RGroupDecompData::relabelRGroup(RGroupData &rgroup,
   std::vector<std::pair<Atom *, Atom *>> atomsToAdd;  // adds -R if necessary
   std::map<int, int> rLabelCoreIndexToAtomicWt;
 
-  for (RWMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
-       ++atIt) {
-    Atom *atom = *atIt;
+  for (auto atom : mol.atoms()) {
     if (atom->hasProp(SIDECHAIN_RLABELS)) {
       atom->setIsotope(0);
       const std::vector<int> &rlabels =
diff --git a/Code/GraphMol/RGroupDecomposition/RGroupUtils.cpp b/Code/GraphMol/RGroupDecomposition/RGroupUtils.cpp
index 90a820938..99d4e081c 100644
--- a/Code/GraphMol/RGroupDecomposition/RGroupUtils.cpp
+++ b/Code/GraphMol/RGroupDecomposition/RGroupUtils.cpp
@@ -127,9 +127,8 @@ bool isAtomWithMultipleNeighborsOrNotDummyRGroupAttachment(const Atom &atom) {
 }
 
 bool hasDummy(const RWMol &core) {
-  for (RWMol::ConstAtomIterator atIt = core.beginAtoms();
-       atIt != core.endAtoms(); ++atIt) {
-    if ((*atIt)->getAtomicNum() == 0) {
+  for (auto atom : core.atoms()) {
+    if (atom->getAtomicNum() == 0) {
       return true;
     }
   }
diff --git a/Code/GraphMol/ReducedGraphs/ReducedGraphs.cpp b/Code/GraphMol/ReducedGraphs/ReducedGraphs.cpp
index c8e10b9c5..d6d09d0fd 100644
--- a/Code/GraphMol/ReducedGraphs/ReducedGraphs.cpp
+++ b/Code/GraphMol/ReducedGraphs/ReducedGraphs.cpp
@@ -146,10 +146,9 @@ RDNumeric::DoubleVector *generateErGFingerprintForReducedGraph(
   // cache the atom type vectors:
   std::vector<std::vector<int>> tvs;
   tvs.reserve(mol.getNumAtoms());
-  for (ROMol::ConstAtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
-       ++atIt) {
+  for (const auto atom : mol.atoms()) {
     const std::vector<int> &tv =
-        (*atIt)->getProp<std::vector<int>>("_ErGAtomTypes");
+        atom->getProp<std::vector<int>>("_ErGAtomTypes");
     tvs.push_back(tv);
   }
 
@@ -212,16 +211,15 @@ ROMol *generateMolExtendedReducedGraph(
   const int aliphaticFlag = atomTypes->size() - 1;  // the last type
   const int aromaticFlag = atomTypes->size();
 
-  for (ROMol::AtomIterator atIt = res->beginAtoms(); atIt != res->endAtoms();
-       ++atIt) {
+  for (const auto atom : res->atoms()) {
     std::vector<int> tv;
     tv.clear();
     for (unsigned int i = 0; i < atomTypes->size(); ++i) {
-      if ((*atomTypes)[i][(*atIt)->getIdx()]) {
+      if ((*atomTypes)[i][atom->getIdx()]) {
         tv.push_back(i);
       }
     }
-    (*atIt)->setProp("_ErGAtomTypes", tv);
+    atom->setProp("_ErGAtomTypes", tv);
   }
 
   // start by adding dummies at the ring centroids
diff --git a/Code/GraphMol/Renumber.cpp b/Code/GraphMol/Renumber.cpp
index 2b94fc3a0..62f38a14e 100644
--- a/Code/GraphMol/Renumber.cpp
+++ b/Code/GraphMol/Renumber.cpp
@@ -65,9 +65,7 @@ ROMol *renumberAtoms(const ROMol &mol,
   }
 
   // now the bonds:
-  for (ROMol::ConstBondIterator bi = mol.beginBonds(); bi != mol.endBonds();
-       ++bi) {
-    const Bond *oBond = (*bi);
+  for (const auto oBond : mol.bonds()) {
     Bond *nBond = oBond->copy();
     nBond->setBeginAtomIdx(revOrder[oBond->getBeginAtomIdx()]);
     nBond->setEndAtomIdx(revOrder[oBond->getEndAtomIdx()]);
diff --git a/Code/GraphMol/Resonance.cpp b/Code/GraphMol/Resonance.cpp
index 119b56e76..d3da2d34d 100644
--- a/Code/GraphMol/Resonance.cpp
+++ b/Code/GraphMol/Resonance.cpp
@@ -85,7 +85,10 @@ class ConjElectrons {
     HAVE_CATION = (1 << 1),
     HAVE_ANION = (1 << 2)
   } ConjElectronsFlags;
-  typedef enum { FP_BONDS = (1 << 0), FP_ATOMS = (1 << 1) } FPFlags;
+  typedef enum {
+    FP_BONDS = (1 << 0),
+    FP_ATOMS = (1 << 1)
+  } FPFlags;
   ConjElectrons(ResonanceMolSupplier *parent, unsigned int groupIdx);
   ConjElectrons(const ConjElectrons &ce);
   ~ConjElectrons();
@@ -179,7 +182,9 @@ class AtomElectrons {
     STACKED = (1 << 2),
     HAS_MULTIPLE_BOND = (1 << 3),
   } AtomElectronsFlags;
-  typedef enum { NEED_CHARGE_BIT = 1 } AllowedBondFlag;
+  typedef enum {
+    NEED_CHARGE_BIT = 1
+  } AllowedBondFlag;
   AtomElectrons(ConjElectrons *parent, const Atom *a);
   AtomElectrons(ConjElectrons *parent, const AtomElectrons &ae);
   ~AtomElectrons() = default;
@@ -227,7 +232,9 @@ class AtomElectrons {
 
 class BondElectrons {
  public:
-  typedef enum { DEFINITIVE = (1 << 0) } BondElectronsFlags;
+  typedef enum {
+    DEFINITIVE = (1 << 0)
+  } BondElectronsFlags;
   BondElectrons(ConjElectrons *parent, const Bond *b);
   BondElectrons(ConjElectrons *parent, const BondElectrons &be);
   ~BondElectrons() = default;
@@ -285,11 +292,10 @@ void sanitizeMol(RWMol &mol) {
 // resonance structure which has been assembled
 void fixExplicitImplicitHs(ROMol &mol) {
   mol.clearComputedProps(false);
-  for (ROMol::AtomIterator ai = mol.beginAtoms(); ai != mol.endAtoms(); ++ai) {
-    (*ai)->clearComputedProps();
-    (*ai)->setNumExplicitHs((*ai)->getNumImplicitHs() +
-                            (*ai)->getNumExplicitHs());
-    (*ai)->updatePropertyCache();
+  for (const auto atom : mol.atoms()) {
+    atom->clearComputedProps();
+    atom->setNumExplicitHs(atom->getNumImplicitHs() + atom->getNumExplicitHs());
+    atom->updatePropertyCache();
   }
 }
 }  // end of namespace ResonanceUtils
diff --git a/Code/GraphMol/SLNParse/SLNAttribs.cpp b/Code/GraphMol/SLNParse/SLNAttribs.cpp
index 01f6f12fe..407c20f34 100644
--- a/Code/GraphMol/SLNParse/SLNAttribs.cpp
+++ b/Code/GraphMol/SLNParse/SLNAttribs.cpp
@@ -470,10 +470,9 @@ void parseMolAttribs(ROMol *mol, AttribListType attribs) {
 }
 
 void adjustAtomChiralities(RWMol *mol) {
-  for (RWMol::AtomIterator atomIt = mol->beginAtoms();
-       atomIt != mol->endAtoms(); atomIt++) {
+  for (auto atom : mol->atoms()) {
     std::string attribVal;
-    if ((*atomIt)->getPropIfPresent(common_properties::_SLN_s, attribVal)) {
+    if (atom->getPropIfPresent(common_properties::_SLN_s, attribVal)) {
       // the atom is marked as chiral, translate the sln chirality into
       // RDKit chirality
 
@@ -483,17 +482,14 @@ void adjustAtomChiralities(RWMol *mol) {
       //      ClC[s=n]H(F)Br  <->  Cl[C@H](F)Br (CHI_TETRAHEDRAL_CCW)
       //      FC[1:s=n](Cl)OCH2@1   <->  F[C@@]1(Cl)OC1 (CHI_TETRAHEDRAL_CW)
       if (attribVal[0] == 'n') {
-        (*atomIt)->setChiralTag(Atom::CHI_TETRAHEDRAL_CW);
+        atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CW);
       } else if (attribVal[0] == 'i') {
-        (*atomIt)->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
+        atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
       }
       std::list<std::pair<int, int>> neighbors;
-      RWMol::ADJ_ITER nbrIdx, endNbrs;
-      boost::tie(nbrIdx, endNbrs) = mol->getAtomNeighbors(*atomIt);
-      while (nbrIdx != endNbrs) {
-        Bond *nbrBond = mol->getBondBetweenAtoms((*atomIt)->getIdx(), *nbrIdx);
-        neighbors.emplace_back(*nbrIdx, nbrBond->getIdx());
-        ++nbrIdx;
+      for (auto nbrBond : mol->atomBonds(atom)) {
+        neighbors.emplace_back(nbrBond->getOtherAtomIdx(atom->getIdx()),
+                               nbrBond->getIdx());
       }
 
       // std::cerr << "CHIRAL " << (*atomIt)->getIdx();
@@ -512,13 +508,12 @@ void adjustAtomChiralities(RWMol *mol) {
       // ok, we now have the ordering of the bonds (used for RDKit chirality),
       // figure out the permutation order relative to the atom numbering
       // (sln chirality):
-      int nSwaps = (*atomIt)->getPerturbationOrder(bondOrdering);
+      int nSwaps = atom->getPerturbationOrder(bondOrdering);
 
       if (nSwaps % 2) {
-        (*atomIt)->setChiralTag((*atomIt)->getChiralTag() ==
-                                        Atom::CHI_TETRAHEDRAL_CW
-                                    ? Atom::CHI_TETRAHEDRAL_CCW
-                                    : Atom::CHI_TETRAHEDRAL_CW);
+        atom->setChiralTag(atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW
+                               ? Atom::CHI_TETRAHEDRAL_CCW
+                               : Atom::CHI_TETRAHEDRAL_CW);
       }
     }
   }
diff --git a/Code/GraphMol/SLNParse/SLNParse.cpp b/Code/GraphMol/SLNParse/SLNParse.cpp
index cccf28b1b..0ad8eed12 100644
--- a/Code/GraphMol/SLNParse/SLNParse.cpp
+++ b/Code/GraphMol/SLNParse/SLNParse.cpp
@@ -81,12 +81,10 @@ void finalizeQueryMol(ROMol *mol, bool mergeHs) {
 
   // do we need to remove the Hs from the molecule?
   if (mergeHs) {
-    for (ROMol::AtomIterator atomIt = mol->beginAtoms();
-         atomIt != mol->endAtoms(); ++atomIt) {
+    for (auto atom : mol->atoms()) {
       // set a query for the H count:
-      if ((*atomIt)->getNumExplicitHs()) {
-        (*atomIt)->expandQuery(
-            makeAtomHCountQuery((*atomIt)->getNumExplicitHs()));
+      if (atom->getNumExplicitHs()) {
+        atom->expandQuery(makeAtomHCountQuery(atom->getNumExplicitHs()));
       }
     }
   }
@@ -96,16 +94,15 @@ void finalizeQueryMol(ROMol *mol, bool mergeHs) {
   VECT_INT_VECT sssr;
   MolOps::symmetrizeSSSR(*mol, sssr);
   int rootIdx = -1;
-  for (ROMol::AtomIterator atomIt = mol->beginAtoms();
-       atomIt != mol->endAtoms(); ++atomIt) {
-    SLNParse::parseFinalAtomAttribs(*atomIt, true);
-    if ((*atomIt)->hasProp(common_properties::_starred)) {
+  for (auto atom : mol->atoms()) {
+    SLNParse::parseFinalAtomAttribs(atom, true);
+    if (atom->hasProp(common_properties::_starred)) {
       if (rootIdx > -1) {
         BOOST_LOG(rdErrorLog) << "SLN Error: multiple starred atoms in a "
                                  "recursive query. Extra stars ignored"
                               << std::endl;
       } else {
-        rootIdx = (*atomIt)->getIdx();
+        rootIdx = atom->getIdx();
       }
     }
   }
@@ -162,13 +159,10 @@ RWMol *toMol(std::string inp, bool doQueries, int debugParse) {
     if (molVect.size() > 0) {
       res = molVect[0];
 
-      for (ROMol::BOND_BOOKMARK_MAP::const_iterator bmIt =
-               res->getBondBookmarks()->begin();
-           bmIt != res->getBondBookmarks()->end(); ++bmIt) {
-        if (bmIt->first > 0 &&
-            bmIt->first < static_cast<int>(res->getNumAtoms())) {
+      for (auto [first, _] : *res->getBondBookmarks()) {
+        if (first > 0 && first < static_cast<int>(res->getNumAtoms())) {
           std::stringstream err;
-          err << "SLN Parser error: Ring closure " << bmIt->first
+          err << "SLN Parser error: Ring closure " << first
               << " does not have a corresponding opener.";
           throw SLNParseException(err.str());
         }
@@ -210,9 +204,8 @@ RWMol *SLNToMol(const std::string &sln, bool sanitize, int debugParse) {
 
   RWMol *res = SLNParse::toMol(sln, false, debugParse);
   if (res) {
-    for (ROMol::AtomIterator atomIt = res->beginAtoms();
-         atomIt != res->endAtoms(); ++atomIt) {
-      SLNParse::parseFinalAtomAttribs(*atomIt, false);
+    for (auto atom : res->atoms()) {
+      SLNParse::parseFinalAtomAttribs(atom, false);
     }
     if (sanitize) {
       // we're going to remove explicit Hs from the graph,
diff --git a/Code/GraphMol/StructChecker/StripSmallFragments.cpp b/Code/GraphMol/StructChecker/StripSmallFragments.cpp
index 34c18068a..756025aad 100644
--- a/Code/GraphMol/StructChecker/StripSmallFragments.cpp
+++ b/Code/GraphMol/StructChecker/StripSmallFragments.cpp
@@ -14,7 +14,6 @@
 #include "../SmilesParse/SmilesWrite.h"
 #include "../FileParsers/MolFileStereochem.h"
 
-
 // define snprintf for msvc
 #if _MSC_VER
 #if _MSC_VER < 1900
@@ -92,9 +91,8 @@ bool StripSmallFragments(RWMol &mol, bool verbose) {
       MolOps::clearSingleBondDirFlags(copy);
       MolOps::detectBondStereochemistry(copy);
       MolOps::assignStereochemistry(copy, true, true, true);
-      for (ROMol::AtomIterator atIt = copy.beginAtoms();
-           atIt != copy.endAtoms(); ++atIt) {
-        if ((*atIt)->hasProp(common_properties::_ChiralityPossible)) {
+      for (const auto atom : copy.atoms()) {
+        if (atom->hasProp(common_properties::_ChiralityPossible)) {
           ischiral = true;
           checkChiral = false;
           break;
@@ -105,19 +103,17 @@ bool StripSmallFragments(RWMol &mol, bool verbose) {
 
     // are chiral tags set
     if (checkChiral) {
-      for (ROMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
-           ++atIt) {
-        if ((*atIt)->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
-            (*atIt)->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) {
+      for (const auto atom : mol.atoms()) {
+        if (atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
+            atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) {
           ischiral = true;
           break;
         }
       }
 
-      for (ROMol::BondIterator bondIt = mol.beginBonds();
-           bondIt != mol.endBonds(); ++bondIt) {
-        if ((*bondIt)->getBondDir() == Bond::BEGINDASH ||
-            (*bondIt)->getBondDir() == Bond::BEGINWEDGE) {
+      for (const auto bond : mol.bonds()) {
+        if (bond->getBondDir() == Bond::BEGINDASH ||
+            bond->getBondDir() == Bond::BEGINWEDGE) {
           ischiral = true;
           break;
         }
diff --git a/Code/GraphMol/StructChecker/Tautomer.cpp b/Code/GraphMol/StructChecker/Tautomer.cpp
index 2744296fb..e483e1194 100644
--- a/Code/GraphMol/StructChecker/Tautomer.cpp
+++ b/Code/GraphMol/StructChecker/Tautomer.cpp
@@ -57,16 +57,15 @@ bool StructCheckTautomer::applyTautomer(unsigned it) {
   }
   // scan for completely mapped bonds and replace bond order with mapped bond
   // from to_tautomer
-  for (RDKit::BondIterator_ bond = Mol.beginBonds(); bond != Mol.endBonds();
-       ++bond) {
-    unsigned ti = atomIdxMap[(*bond)->getBeginAtomIdx()];
-    unsigned tj = atomIdxMap[(*bond)->getEndAtomIdx()];
+  for (auto bond : Mol.bonds()) {
+    unsigned ti = atomIdxMap[bond->getBeginAtomIdx()];
+    unsigned tj = atomIdxMap[bond->getEndAtomIdx()];
     if (invalid_idx == ti || invalid_idx == tj) {
       continue;
     }
     const Bond *tb = toTautomer.getBondBetweenAtoms(ti, tj);
-    if (tb && (*bond)->getBondType() != tb->getBondType()) {
-      (*bond)->setBondType(tb->getBondType());
+    if (tb && bond->getBondType() != tb->getBondType()) {
+      bond->setBondType(tb->getBondType());
     }
   }
   // apply charge/radical fixes if any
diff --git a/Code/GraphMol/Subgraphs/Subgraphs.cpp b/Code/GraphMol/Subgraphs/Subgraphs.cpp
index d94b13db3..caac0012b 100644
--- a/Code/GraphMol/Subgraphs/Subgraphs.cpp
+++ b/Code/GraphMol/Subgraphs/Subgraphs.cpp
@@ -463,10 +463,9 @@ findAllPathsOfLengthsMtoN(const ROMol &mol, unsigned int lowerLen,
 
   if (!distMat) {
     // generate the adjacency matrix by hand by looping over the bonds
-    ROMol::ConstBondIterator bondIt;
-    for (bondIt = mol.beginBonds(); bondIt != mol.endBonds(); bondIt++) {
-      Atom *beg = (*bondIt)->getBeginAtom();
-      Atom *end = (*bondIt)->getEndAtom();
+    for (const auto bond : mol.bonds()) {
+      Atom *beg = bond->getBeginAtom();
+      Atom *end = bond->getEndAtom();
       // check for H, which we might be skipping
       if (useHs || (beg->getAtomicNum() != 1 && end->getAtomicNum() != 1)) {
         adjMat[beg->getIdx() * dim + end->getIdx()] = 1;
diff --git a/Code/GraphMol/resMolSupplierTest.cpp b/Code/GraphMol/resMolSupplierTest.cpp
index a35ce7dd5..ad65a65cd 100644
--- a/Code/GraphMol/resMolSupplierTest.cpp
+++ b/Code/GraphMol/resMolSupplierTest.cpp
@@ -627,9 +627,8 @@ void getBtVectVect(ResonanceMolSupplier *resMolSuppl,
     ROMol *resMol = resMolSuppl->next();
     std::vector<unsigned int> bt;
     bt.reserve(resMol->getNumBonds());
-    for (ROMol::BondIterator bi = resMol->beginBonds();
-         bi != resMol->endBonds(); ++bi) {
-      bt.push_back(static_cast<unsigned int>((*bi)->getBondTypeAsDouble()));
+    for (const auto bond : resMol->bonds()) {
+      bt.push_back(static_cast<unsigned int>(bond->getBondTypeAsDouble()));
     }
     btVect2.push_back(bt);
     delete resMol;
diff --git a/External/FreeSASA/RDFreeSASA.cpp b/External/FreeSASA/RDFreeSASA.cpp
index 055134a3d..b654d3729 100644
--- a/External/FreeSASA/RDFreeSASA.cpp
+++ b/External/FreeSASA/RDFreeSASA.cpp
@@ -89,8 +89,7 @@ bool classifyAtoms(ROMol &mol, std::vector<double> &radii,
   }
 
   bool success = true;
-  for (ROMol::AtomIterator at = mol.beginAtoms(); at != mol.endAtoms(); ++at) {
-    Atom *atom = *at;
+  for (const auto atom : mol.atoms()) {
     freesasa_atom_class cls = FREESASA_ATOM_UNKNOWN;
     std::string classification = "Unclassified";
     double radius = 0.0;
@@ -168,9 +167,9 @@ double internalCalcSASA(const ROMol &mol, const std::vector<double> &radii,
   double sasa = res->total;
   mol.setProp(common_properties::Molecule::SASA, sasa);
   size_t i = 0;
-  for (ROMol::ConstAtomIterator at = mol.beginAtoms(); at != mol.endAtoms();
-       ++at, ++i) {
-    (*at)->setProp(common_properties::Atom::SASA, res->sasa[i]);
+  for (const auto atom : mol.atoms()) {
+    atom->setProp(common_properties::Atom::SASA, res->sasa[i]);
+    ++i;
   }
 
   freesasa_result_free(res);
@@ -184,10 +183,10 @@ double calcSASA(const RDKit::ROMol &mol, const std::vector<double> &radii,
   double result = internalCalcSASA(mol, radii, confIdx, opts);
   if (query) {
     result = 0.0f;
-    for (ROMol::ConstQueryAtomIterator at = mol.beginQueryAtoms(query);
-         at != mol.endQueryAtoms(); ++at) {
-      const Atom *atom = *at;
-      result += atom->getProp<double>("SASA");
+    for (const auto atom : mol.atoms()) {
+      if (query->Match(atom)) {
+        result += atom->getProp<double>("SASA");
+      }
     }
   }
   return result;
diff --git a/External/INCHI-API/inchi.cpp b/External/INCHI-API/inchi.cpp
index 63a710aca..dc90ed83f 100644
--- a/External/INCHI-API/inchi.cpp
+++ b/External/INCHI-API/inchi.cpp
@@ -1149,97 +1149,96 @@ bool _Valence3ClCleanUp1(RWMol &mol, Atom *atom) {
 }
 
 void cleanUp(RWMol &mol) {
-  ROMol::AtomIterator ai;
   bool aromHolder;
-  for (ai = mol.beginAtoms(); ai != mol.endAtoms(); ++ai) {
-    switch ((*ai)->getAtomicNum()) {
+  for (const auto atom : mol.atoms()) {
+    switch (atom->getAtomicNum()) {
       case 7:
-        if ((*ai)->calcExplicitValence(false) == 4) {
-          if (_Valence4NCleanUp1(mol, *ai)) {
+        if (atom->calcExplicitValence(false) == 4) {
+          if (_Valence4NCleanUp1(mol, atom)) {
             continue;
           }
-          if ((*ai)->getFormalCharge() == -1) {
-            if (_Valence4NCleanUp2(mol, *ai)) {
+          if (atom->getFormalCharge() == -1) {
+            if (_Valence4NCleanUp2(mol, atom)) {
               continue;
             }
           }
           continue;
         }
 
-        if ((*ai)->getFormalCharge()) {
+        if (atom->getFormalCharge()) {
           continue;
         }
-        aromHolder = (*ai)->getIsAromatic();
-        (*ai)->setIsAromatic(0);
+        aromHolder = atom->getIsAromatic();
+        atom->setIsAromatic(0);
 
-        if ((*ai)->calcExplicitValence(false) == 5) {
+        if (atom->calcExplicitValence(false) == 5) {
           // rings CN1=CCN=CC=1, CN1=NCOCC=1, [N]=C1N=CN=N1, [N]=C1C=CN=N1
-          (_Valence5NCleanUp6(mol, *ai)) || (_Valence5NCleanUp7(mol, *ai)) ||
-              (_Valence5NCleanUp8(mol, *ai)) ||
-              (_Valence5NCleanUp9(mol, *ai)) ||
-              (_Valence5NCleanUpA(mol, *ai)) ||
+          (_Valence5NCleanUp6(mol, atom)) || (_Valence5NCleanUp7(mol, atom)) ||
+              (_Valence5NCleanUp8(mol, atom)) ||
+              (_Valence5NCleanUp9(mol, atom)) ||
+              (_Valence5NCleanUpA(mol, atom)) ||
               // try search for valence-5 N connected to a N+
-              (_Valence5NCleanUp1(mol, *ai)) ||
+              (_Valence5NCleanUp1(mol, atom)) ||
               // connected to N- through a tiple then single bond
-              (_Valence5NCleanUp2(mol, *ai)) ||
+              (_Valence5NCleanUp2(mol, atom)) ||
               // directly to a N
-              (_Valence5NCleanUp3(mol, *ai)) ||
+              (_Valence5NCleanUp3(mol, atom)) ||
               // to two Si- via double bonds
-              (_Valence5NCleanUp4(mol, *ai)) ||
+              (_Valence5NCleanUp4(mol, atom)) ||
               // alternating bonds to O
-              (_Valence5NCleanUp5(mol, *ai, 8)) ||
+              (_Valence5NCleanUp5(mol, atom, 8)) ||
               // alternating bonds to S
-              (_Valence5NCleanUp5(mol, *ai, 16)) ||
+              (_Valence5NCleanUp5(mol, atom, 16)) ||
               // alternating bonds to S
-              (_Valence5NCleanUp5(mol, *ai, 9)) ||
+              (_Valence5NCleanUp5(mol, atom, 9)) ||
               // alternating bonds to S
-              (_Valence5NCleanUp5(mol, *ai, 17)) ||
+              (_Valence5NCleanUp5(mol, atom, 17)) ||
               // last resort
-              (_Valence5NCleanUpB(mol, *ai));
+              (_Valence5NCleanUpB(mol, atom));
         }
         if (aromHolder) {
-          (*ai)->setIsAromatic(1);
+          atom->setIsAromatic(1);
         }
         break;
       case 17:
-        if ((*ai)->calcExplicitValence(false) == 8 &&
-            _Valence8ClCleanUp1(mol, *ai)) {
+        if (atom->calcExplicitValence(false) == 8 &&
+            _Valence8ClCleanUp1(mol, atom)) {
           continue;
         }
-        if ((*ai)->calcExplicitValence(false) == 5 &&
-            _Valence5ClCleanUp1(mol, *ai)) {
+        if (atom->calcExplicitValence(false) == 5 &&
+            _Valence5ClCleanUp1(mol, atom)) {
           continue;
         }
-        if ((*ai)->calcExplicitValence(false) == 3 &&
-            _Valence3ClCleanUp1(mol, *ai)) {
+        if (atom->calcExplicitValence(false) == 3 &&
+            _Valence3ClCleanUp1(mol, atom)) {
           continue;
         }
         break;
       case 16:
-        if ((*ai)->calcExplicitValence(false) == 7) {
-          if (_Valence7SCleanUp1(mol, *ai)) {
+        if (atom->calcExplicitValence(false) == 7) {
+          if (_Valence7SCleanUp1(mol, atom)) {
             continue;
           }
-          if (_Valence7SCleanUp2(mol, *ai)) {
+          if (_Valence7SCleanUp2(mol, atom)) {
             continue;
           }
-          if (_Valence7SCleanUp3(mol, *ai)) {
+          if (_Valence7SCleanUp3(mol, atom)) {
             continue;
           }
-          _Valence8SCleanUp1(mol, *ai);
-        } else if ((*ai)->calcExplicitValence(false) == 8) {
-          _Valence8SCleanUp1(mol, *ai);
+          _Valence8SCleanUp1(mol, atom);
+        } else if (atom->calcExplicitValence(false) == 8) {
+          _Valence8SCleanUp1(mol, atom);
         }
         break;
       case 35:
-        if ((*ai)->calcExplicitValence(false) == 3 &&
-            (*ai)->getFormalCharge() == 0) {
+        if (atom->calcExplicitValence(false) == 3 &&
+            atom->getFormalCharge() == 0) {
           // connected to Se. Example: PubChem 10787526
-          if ((*ai)->getDegree() == 1) {
+          if (atom->getDegree() == 1) {
             RWMol::ADJ_ITER nid, end;
-            boost::tie(nid, end) = mol.getAtomNeighbors(*ai);
+            boost::tie(nid, end) = mol.getAtomNeighbors(atom);
             if (mol.getAtomWithIdx(*nid)->getAtomicNum() == 34) {
-              mol.getBondBetweenAtoms((*ai)->getIdx(), *nid)
+              mol.getBondBetweenAtoms(atom->getIdx(), *nid)
                   ->setBondType(Bond::SINGLE);
             }
           }