perf: Start using cells for bonding calculation

src/classes/atom.cpp

@@ -150,3 +150,9 @@ void AtomBase::deserialise(const SerialisedValue &node)
     set(Deserialisable::deser<Elements::Element>(node.at("z")), Deserialisable::deser<Vector3>(node.at("r")),
         Deserialisable::deser_or<double>(node, "q", 0));
 }
+
+// Set cell in which the atom exists
+void AtomBase::setCell(Cell *cell) { cell_ = cell; }
+
+// Return cell in which the atom exists
+Cell *AtomBase::cell() const { return cell_; }

src/classes/atom.h

@@ -9,6 +9,8 @@
 #include "data/elements.h"
 #include "math/vector3.h"
 
+class Cell;
+
 // Atom Base
 class AtomBase
 {
@@ -26,6 +28,8 @@ class AtomBase
     int index_{-1};
     // Atom type index in parent object
     int atomTypeIndex_{AtomConstants::TypeIndex::Ignore};
+    // Cell in which the atom exists
+    Cell *cell_{nullptr};
 
     public:
     // Set basic properties
@@ -52,6 +56,10 @@ class AtomBase
     void setAtomTypeIndex(int id);
     // Return associated atom type index
     int atomTypeIndex() const;
+    // Set cell in which the atom exists
+    void setCell(Cell *cell);
+    // Return cell in which the atom exists
+    Cell *cell() const;
 
     /*
      * General Connectivity

src/classes/cell.cpp

@@ -35,14 +35,14 @@ const Vector3 &Cell::centre() const { return centre_; }
  */
 
 // Return vector of contained Atoms
-std::vector<ConfigurationAtom *> &Cell::atoms() { return atoms_; }
-const std::vector<ConfigurationAtom *> &Cell::atoms() const { return atoms_; }
+std::vector<AtomBase *> &Cell::atoms() { return atoms_; }
+const std::vector<AtomBase *> &Cell::atoms() const { return atoms_; }
 
 // Return number of Atoms in list
 int Cell::nAtoms() const { return atoms_.size(); }
 
 // Add atom to Cell
-void Cell::addAtom(ConfigurationAtom *atom)
+void Cell::addAtom(AtomBase *atom)
 {
     assert(atom);
     assert(!atom->cell());
@@ -52,7 +52,7 @@ void Cell::addAtom(ConfigurationAtom *atom)
 }
 
 // Remove Atom from Cell
-void Cell::removeAtom(ConfigurationAtom *atom)
+void Cell::removeAtom(AtomBase *atom)
 {
     auto it = std::find(atoms_.begin(), atoms_.end(), atom);
     assert(it != atoms_.end());

src/classes/cell.h

@@ -47,18 +47,18 @@ class Cell
      */
     private:
     // Vector of Atoms contained in this Cell
-    std::vector<ConfigurationAtom *> atoms_;
+    std::vector<AtomBase *> atoms_;
 
     public:
     // Return array of contained Atoms
-    std::vector<ConfigurationAtom *> &atoms();
-    const std::vector<ConfigurationAtom *> &atoms() const;
+    std::vector<AtomBase *> &atoms();
+    const std::vector<AtomBase *> &atoms() const;
     // Return number of Atoms in array
     int nAtoms() const;
     // Add atom to Cell
-    void addAtom(ConfigurationAtom *atom);
+    void addAtom(AtomBase *atom);
     // Remove Atom from Cell
-    void removeAtom(ConfigurationAtom *atom);
+    void removeAtom(AtomBase *atom);
     // Clear all atoms from cell
     void clearAtoms();
 };

src/classes/cellArray.cpp

@@ -165,6 +165,16 @@ const std::vector<CellNeighbour> &CellArray::neighbours(const Cell &cell) const
 // Return vector of all unique cell neighbour pairs
 const std::vector<CellNeighbourPair> &CellArray::getCellNeighbourPairs() const { return neighbourPairs_; }
 
+// Return vector of all unique cell neighbour pairs with self cell pairs
+const std::vector<CellNeighbourPair> CellArray::getCellNeighbourPairsWithSelf() const
+{
+    auto result = neighbourPairs_;
+    for (auto &cell : cells_)
+        result.emplace_back(cell, cell, false);
+
+    return result;
+}
+
 // Return whether minimum image calculation is required between the supplied cells
 bool CellArray::minimumImageRequired(const Cell &a, const Cell &b) const
 {

src/classes/cellArray.h

@@ -110,6 +110,8 @@ class CellArray
     const std::vector<CellNeighbour> &neighbours(const Cell &cell) const;
     // Return vector of all unique cell neighbour pairs
     const std::vector<CellNeighbourPair> &getCellNeighbourPairs() const;
+    // Return vector of all unique cell neighbour pairs with self cell pairs
+    const std::vector<CellNeighbourPair> getCellNeighbourPairsWithSelf() const;
     // Return whether minimum image calculation is required between the supplied cells
     bool minimumImageRequired(const Cell &a, const Cell &b) const;
 

src/classes/configurationAtom.cpp

@@ -23,12 +23,6 @@ void ConfigurationAtom::setMolecule(std::shared_ptr<Molecule> mol) { molecule_ =
 // Return Molecule in which this Atom exists
 const std::shared_ptr<Molecule> &ConfigurationAtom::molecule() const { return molecule_; }
 
-// Set cell in which the atom exists
-void ConfigurationAtom::setCell(Cell *cell) { cell_ = cell; }
-
-// Return cell in which the atom exists
-Cell *ConfigurationAtom::cell() const { return cell_; }
-
 /*
  * Intramolecular Information
  */

src/classes/configurationAtom.h

@@ -24,8 +24,6 @@ class ConfigurationAtom : public Atom<SpeciesBond>
     const SpeciesAtom *speciesAtom_{nullptr};
     // Molecule in which this Atom exists
     std::shared_ptr<Molecule> molecule_{nullptr};
-    // Cell in which the atom exists
-    Cell *cell_{nullptr};
 
     public:
     // Set SpeciesAtom that this Atom represents
@@ -36,10 +34,6 @@ class ConfigurationAtom : public Atom<SpeciesBond>
     void setMolecule(std::shared_ptr<Molecule> mol);
     // Return Molecule in which this Atom exists
     const std::shared_ptr<Molecule> &molecule() const;
-    // Set cell in which the atom exists
-    void setCell(Cell *cell);
-    // Return cell in which the atom exists
-    Cell *cell() const;
 
     /*
      * Intramolecular Information

src/kernels/energy.cpp

@@ -5,6 +5,7 @@
 #include "classes/box.h"
 #include "classes/cell.h"
 #include "classes/configuration.h"
+#include "classes/configurationAtom.h"
 #include "classes/molecule.h"
 #include "classes/potentialMap.h"
 #include "classes/species.h"
@@ -46,15 +47,15 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellEnergy(const Cell &cell, bool
 
     for (auto i = 0; i < atoms.size(); ++i)
     {
-        auto &ii = atoms[i];
+        auto ii = static_cast<ConfigurationAtom *>(atoms[i]);
         auto molI = ii->molecule();
         auto &rI = ii->r();
 
         // Straight loop over other cell atoms
         for (auto j = i + 1; j < atoms.size(); ++j)
         {
             // Calculate rSquared distance between atoms, and check it against the stored cutoff distance
-            auto &jj = atoms[j];
+            auto jj = static_cast<ConfigurationAtom *>(atoms[j]);
 
             auto rSq = (rI - jj->r()).magnitudeSq();
             if (rSq > cutoffDistanceSquared_)
@@ -90,13 +91,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell &cent
     // Loop over central cell atoms
     if (applyMim)
     {
-        for (auto &ii : centralAtoms)
+        for (auto ii : centralAtoms | castView<const ConfigurationAtom *>())
         {
             auto molI = ii->molecule();
             auto &rI = ii->r();
 
             // Straight loop over other cell atoms
-            for (const auto &jj : otherAtoms)
+            for (const auto jj : otherAtoms | castView<const ConfigurationAtom *>())
             {
                 // Calculate rSquared distance between atoms, and check it against the stored cutoff distance
                 auto rSq = box_.minimumDistanceSquared(rI, jj->r());
@@ -122,13 +123,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell &cent
     }
     else
     {
-        for (auto &ii : centralAtoms)
+        for (auto ii : centralAtoms | castView<const ConfigurationAtom *>())
         {
             auto &molI = ii->molecule();
             auto &rI = ii->r();
 
             // Straight loop over other cell atoms
-            for (const auto &jj : otherAtoms)
+            for (const auto jj : otherAtoms | castView<const ConfigurationAtom *>())
             {
                 // Calculate rSquared distance between atoms, and check it against the stored cutoff distance
                 auto rSq = (rI - jj->r()).magnitudeSq();
@@ -172,7 +173,7 @@ double EnergyKernel::pairPotentialEnergy(const ConfigurationAtom &i) const
                                           return std::accumulate(nbrCellAtoms.begin(), nbrCellAtoms.end(), 0.0,
                                                                  [&i, mimRequired, this](const auto innerAcc, const auto *j)
                                                                  {
-                                                                     auto &jj = *j;
+                                                                     auto &jj = *static_cast<const ConfigurationAtom *>(j);
 
                                                                      // Calculate rSquared distance between atoms, and check it
                                                                      // against the stored cutoff distance
@@ -230,7 +231,7 @@ Kernel::PairPotentialEnergyValue EnergyKernel::pairPotentialEnergy(const Molecul
                                 return acc + std::accumulate(nbrCellAtoms.begin(), nbrCellAtoms.end(), 0.0,
                                                              [&ii, mimRequired, this](const auto innerAcc, const auto *j)
                                                              {
-                                                                 auto &jj = *j;
+                                                                 auto &jj = *static_cast<const ConfigurationAtom *>(j);
 
                                                                  // Don't consider atoms within the target molecule here - add
                                                                  // it on afterwards

src/kernels/force.cpp

@@ -96,24 +96,24 @@ void ForceKernel::cellToCellPairPotentialForces(const Cell *centralCell, const C
     // Loop over all atom pairs excluding any within the same molecule
     if (applyMim)
     {
-        for (const auto &i : centralAtoms)
+        for (const auto &i : centralAtoms | castView<ConfigurationAtom *>())
         {
             molI = i->molecule();
             auto indexI = molI->globalAtomIndex(i);
 
-            for (auto *j : otherAtoms)
+            for (auto *j : otherAtoms | castView<ConfigurationAtom *>())
                 if (molI != j->molecule())
                     forcesWithMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f);
         }
     }
     else
     {
-        for (const auto &i : centralAtoms)
+        for (const auto &i : centralAtoms | castView<ConfigurationAtom *>())
         {
             molI = i->molecule();
             auto indexI = molI->globalAtomIndex(i);
 
-            for (auto *j : otherAtoms)
+            for (auto *j : otherAtoms | castView<ConfigurationAtom *>())
                 if (molI != j->molecule())
                     forcesWithoutMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f);
         }
@@ -160,7 +160,7 @@ void ForceKernel::totalForces(std::vector<Vector3> &ppForceVector, std::vector<V
             auto &fLocal = combinablePP.local();
 
             // Interatomic interactions between atoms in this cell, excluding those within the same molecule
-            dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms(),
+            dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms() | castView<ConfigurationAtom *>(),
                                     [&](int indexI, const auto &i, int indexJ, const auto &j)
                                     {
                                         if (indexI == indexJ)

src/nodes/calculateBonding.cpp

@@ -2,6 +2,7 @@
 // Copyright (c) 2026 Team Dissolve and contributors
 
 #include "nodes/calculateBonding.h"
+#include "classes/cellArray.h"
 #include "data/atomicRadii.h"
 #include "main/dissolve.h"
 #include "templates/parallelDefs.h"
@@ -52,57 +53,65 @@ void CalculateBondingNode::calculate(Structure &structure, double tolerance, boo
         structure.clearBonds();
 
     auto box = structure.box();
-    auto nAtoms = PairIterator(structure.nAtoms());
 
-    PairIterator pairs(structure.nAtoms());
+    CellArray cells;
+    cells.generate(box, 6.0);
+
+    // Populate cells
+    for (auto &atom : structure.atoms())
+        cells.cell(atom->r())->addAtom(atom.get());
 
     // Look at two indices and check to see if it would be a vaild
     // bond.  If so, return a list with that bond, otherwise an empty
     // list.  From an information theory standpoint, a list of at most
     // one element is identical to a std::optional, but the lists can
     // be trivially combined during the reduce part of
     // transform_reduce
-    auto validBond = [&structure, &box, tolerance, preventMetallic,
-                      clearBefore](std::tuple<int, int> pair) -> std::vector<std::tuple<StructureAtom *, StructureAtom *>>
+    auto validBond = [&box, tolerance, preventMetallic](const CellNeighbourPair idx) -> std::set<std::pair<int, int>>
     {
-        auto [indexI, indexJ] = pair;
-        if (indexI == indexJ)
-            return {};
-        auto i = structure.atom(indexI);
-        // Get StructureAtom 'i' and its radius
-        auto radiusI = AtomicRadii::radius(i->Z());
-        // Get StructureAtom 'j'
-        auto j = structure.atom(indexJ);
-
-        // If the two atoms are both metal ions and prevent metallic bonds = true, continue
-        if (preventMetallic && Elements::isMetallic(i->Z()) && Elements::isMetallic(j->Z()))
-            return {};
-
-        // Calculate distance between atoms
-        auto r = box.minimumDistance(j->r(), i->r());
-
-        // Compare distance to sum of atomic radii (multiplied by tolerance factor)
-        if (r > (radiusI + AtomicRadii::radius(j->Z())) * tolerance)
-            return {};
-
-        if (structure.getBond(i, j))
-            return {};
-
-        return {{i, j}};
+        auto cellI = idx.cell;
+        auto cellJ = idx.neighbour;
+
+        std::set<std::pair<int, int>> result;
+
+        for (auto i : cellI.atoms() | castView<StructureAtom *>())
+            for (auto j : cellJ.atoms() | castView<StructureAtom *>())
+            {
+                // Don't bond atoms to themselves
+                if (i == j)
+                    continue;
+
+                auto radiusI = AtomicRadii::radius(i->Z());
+
+                // If the two atoms are both metal ions and prevent metallic bonds = true, continue
+                if (preventMetallic && Elements::isMetallic(i->Z()) && Elements::isMetallic(j->Z()))
+                    continue;
+
+                // Calculate distance between atoms
+                auto r = box.minimumDistance(j->r(), i->r());
+
+                // Compare distance to sum of atomic radii (multiplied by tolerance factor)
+                if (r > (radiusI + AtomicRadii::radius(j->Z())) * tolerance)
+                    continue;
+
+                result.insert({std::min(i->index(), j->index()), std::max(i->index(), j->index())});
+            }
+
+        return result;
     };
 
     // Combine two lists of bonds into a single list
     auto joinBonds = [](auto a, auto b)
     {
         auto ab = a;
-        ab.insert(ab.end(), b.begin(), b.end());
+        ab.merge(b);
         return ab;
     };
 
-    // Create an empty vector of the correct shape
-    std::vector<std::tuple<StructureAtom *, StructureAtom *>> empty;
     // In parallel, construct the list of the bonds that need to be added
-    auto results = std::transform_reduce(ParallelPolicies::par_unseq, pairs.begin(), pairs.end(), empty, joinBonds, validBond);
+    auto pairs = cells.getCellNeighbourPairsWithSelf();
+    auto results = std::transform_reduce(ParallelPolicies::par_unseq, pairs.begin(), pairs.end(),
+                                         std::set<std::pair<int, int>>(), joinBonds, validBond);
 
     // Add the bonds serially
     for (auto [i, j] : results)

src/templates/algorithms.h

@@ -304,3 +304,18 @@ template <class T, template <class...> class U> inline constexpr bool is_instanc
 template <template <class...> class U, class... Vs> inline constexpr bool is_instance_of_v<U<Vs...>, U> = std::true_type{};
 template <typename> constexpr bool is_optional = false;
 template <typename T> constexpr bool is_optional<std::optional<T>> = true;
+
+// View that statically casts the elements in a range into a new type
+// given by the template.  For example,
+//
+// for (auto *atom : atoms | castView<ConfigurationAtom *>())
+//    foo(atom)
+//
+// is equivalent to
+//
+// for (auto atom : atoms)
+//     for(static_cast<ConfigurationAtom *>(atom))
+template <typename T> auto castView()
+{
+    return std::ranges::views::transform([](auto x) { return static_cast<T>(x); });
+}