From 3ef0fad3373fd0b6d289dc41a0598a43085b463e Mon Sep 17 00:00:00 2001 From: Adam Washington Date: Wed, 23 Sep 2026 12:39:30 +0100 Subject: [PATCH 1/4] perf: Start using cells for bonding calculation --- src/classes/atom.cpp | 6 +++ src/classes/atom.h | 8 ++++ src/classes/cell.cpp | 8 ++-- src/classes/cell.h | 10 ++-- src/classes/cellArray.cpp | 10 ++++ src/classes/cellArray.h | 2 + src/classes/configurationAtom.cpp | 6 --- src/classes/configurationAtom.h | 6 --- src/kernels/energy.cpp | 20 ++++---- src/kernels/force.cpp | 12 +++-- src/nodes/calculateBonding.cpp | 80 +++++++++++++++++++------------ 11 files changed, 104 insertions(+), 64 deletions(-) diff --git a/src/classes/atom.cpp b/src/classes/atom.cpp index 14ecb6644e..b0f0901091 100644 --- a/src/classes/atom.cpp +++ b/src/classes/atom.cpp @@ -150,3 +150,9 @@ void AtomBase::deserialise(const SerialisedValue &node) set(Deserialisable::deser(node.at("z")), Deserialisable::deser(node.at("r")), Deserialisable::deser_or(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_; } diff --git a/src/classes/atom.h b/src/classes/atom.h index 4357241f5e..4acca28732 100644 --- a/src/classes/atom.h +++ b/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 diff --git a/src/classes/cell.cpp b/src/classes/cell.cpp index 9b7a392252..91e1da615e 100644 --- a/src/classes/cell.cpp +++ b/src/classes/cell.cpp @@ -35,14 +35,14 @@ const Vector3 &Cell::centre() const { return centre_; } */ // Return vector of contained Atoms -std::vector &Cell::atoms() { return atoms_; } -const std::vector &Cell::atoms() const { return atoms_; } +std::vector &Cell::atoms() { return atoms_; } +const std::vector &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()); diff --git a/src/classes/cell.h b/src/classes/cell.h index a0bd8e78ef..dbd65dbdb0 100644 --- a/src/classes/cell.h +++ b/src/classes/cell.h @@ -47,18 +47,18 @@ class Cell */ private: // Vector of Atoms contained in this Cell - std::vector atoms_; + std::vector atoms_; public: // Return array of contained Atoms - std::vector &atoms(); - const std::vector &atoms() const; + std::vector &atoms(); + const std::vector &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(); }; diff --git a/src/classes/cellArray.cpp b/src/classes/cellArray.cpp index ad9f9fef7e..6d85850467 100644 --- a/src/classes/cellArray.cpp +++ b/src/classes/cellArray.cpp @@ -165,6 +165,16 @@ const std::vector &CellArray::neighbours(const Cell &cell) const // Return vector of all unique cell neighbour pairs const std::vector &CellArray::getCellNeighbourPairs() const { return neighbourPairs_; } +// Return vector of all unique cell neighbour pairs with self cell pairs +const std::vector 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 { diff --git a/src/classes/cellArray.h b/src/classes/cellArray.h index 9879dff8a9..40b8d0da93 100644 --- a/src/classes/cellArray.h +++ b/src/classes/cellArray.h @@ -110,6 +110,8 @@ class CellArray const std::vector &neighbours(const Cell &cell) const; // Return vector of all unique cell neighbour pairs const std::vector &getCellNeighbourPairs() const; + // Return vector of all unique cell neighbour pairs with self cell pairs + const std::vector getCellNeighbourPairsWithSelf() const; // Return whether minimum image calculation is required between the supplied cells bool minimumImageRequired(const Cell &a, const Cell &b) const; diff --git a/src/classes/configurationAtom.cpp b/src/classes/configurationAtom.cpp index af4af168ca..f5d050bb20 100644 --- a/src/classes/configurationAtom.cpp +++ b/src/classes/configurationAtom.cpp @@ -23,12 +23,6 @@ void ConfigurationAtom::setMolecule(std::shared_ptr mol) { molecule_ = // Return Molecule in which this Atom exists const std::shared_ptr &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 */ diff --git a/src/classes/configurationAtom.h b/src/classes/configurationAtom.h index 07598e34d1..6b4126548c 100644 --- a/src/classes/configurationAtom.h +++ b/src/classes/configurationAtom.h @@ -24,8 +24,6 @@ class ConfigurationAtom : public Atom const SpeciesAtom *speciesAtom_{nullptr}; // Molecule in which this Atom exists std::shared_ptr 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 void setMolecule(std::shared_ptr mol); // Return Molecule in which this Atom exists const std::shared_ptr &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 diff --git a/src/kernels/energy.cpp b/src/kernels/energy.cpp index b3fcd21a00..80e5be289a 100644 --- a/src/kernels/energy.cpp +++ b/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" @@ -12,6 +13,9 @@ #include #include +auto castToConstConfiguration = + std::ranges::views::transform([](const AtomBase *x) { return static_cast(x); }); + EnergyKernel::EnergyKernel(const Configuration *cfg, const PotentialMap &potentialMap) : GeometryKernel(cfg->box(), potentialMap), configuration_(cfg) { @@ -46,7 +50,7 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellEnergy(const Cell &cell, bool for (auto i = 0; i < atoms.size(); ++i) { - auto &ii = atoms[i]; + auto ii = static_cast(atoms[i]); auto molI = ii->molecule(); auto &rI = ii->r(); @@ -54,7 +58,7 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellEnergy(const Cell &cell, bool 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(atoms[j]); auto rSq = (rI - jj->r()).magnitudeSq(); if (rSq > cutoffDistanceSquared_) @@ -90,13 +94,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell ¢ // Loop over central cell atoms if (applyMim) { - for (auto &ii : centralAtoms) + for (auto ii : centralAtoms | castToConstConfiguration) { auto molI = ii->molecule(); auto &rI = ii->r(); // Straight loop over other cell atoms - for (const auto &jj : otherAtoms) + for (const auto jj : otherAtoms | castToConstConfiguration) { // Calculate rSquared distance between atoms, and check it against the stored cutoff distance auto rSq = box_.minimumDistanceSquared(rI, jj->r()); @@ -122,13 +126,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell ¢ } else { - for (auto &ii : centralAtoms) + for (auto ii : centralAtoms | castToConstConfiguration) { auto &molI = ii->molecule(); auto &rI = ii->r(); // Straight loop over other cell atoms - for (const auto &jj : otherAtoms) + for (const auto jj : otherAtoms | castToConstConfiguration) { // Calculate rSquared distance between atoms, and check it against the stored cutoff distance auto rSq = (rI - jj->r()).magnitudeSq(); @@ -172,7 +176,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(j); // Calculate rSquared distance between atoms, and check it // against the stored cutoff distance @@ -230,7 +234,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(j); // Don't consider atoms within the target molecule here - add // it on afterwards diff --git a/src/kernels/force.cpp b/src/kernels/force.cpp index d66f63a0b4..6640fdb28b 100644 --- a/src/kernels/force.cpp +++ b/src/kernels/force.cpp @@ -11,6 +11,8 @@ #include "templates/algorithms.h" #include +auto castToConfiguration = std::ranges::views::transform([](AtomBase *x) { return static_cast(x); }); + ForceKernel::ForceKernel(const Configuration *cfg, const PotentialMap &potentialMap) : GeometryKernel(cfg->box(), potentialMap), configuration_(cfg) { @@ -96,24 +98,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 | castToConfiguration) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms) + for (auto *j : otherAtoms | castToConfiguration) if (molI != j->molecule()) forcesWithMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f); } } else { - for (const auto &i : centralAtoms) + for (const auto &i : centralAtoms | castToConfiguration) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms) + for (auto *j : otherAtoms | castToConfiguration) if (molI != j->molecule()) forcesWithoutMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f); } @@ -160,7 +162,7 @@ void ForceKernel::totalForces(std::vector &ppForceVector, std::vectoratoms(), + dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms() | castToConfiguration, [&](int indexI, const auto &i, int indexJ, const auto &j) { if (indexI == indexJ) diff --git a/src/nodes/calculateBonding.cpp b/src/nodes/calculateBonding.cpp index 05aea2ddb1..451ccf8ecd 100644 --- a/src/nodes/calculateBonding.cpp +++ b/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,9 +53,20 @@ 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, box.inscribedSphereRadius() * 2.0 / 3.0); + + // Populate cells + // + // As an ugly hack for the moment, we'll shove the StructureAtom + // pointers into the ConfigurationAtom points that the cell + // expects. + for (auto &atom : structure.atoms()) + cells.cell(atom->r())->addAtom(atom.get()); + + // We'll eventually need a view to transform these pointers back + auto castToStructure = std::ranges::views::transform([](const auto x) { return (StructureAtom *)x; }); // 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 @@ -62,46 +74,54 @@ void CalculateBondingNode::calculate(Structure &structure, double tolerance, boo // 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 pair) -> std::vector> + auto validBond = [&castToStructure, &structure, &box, tolerance, preventMetallic, + clearBefore](const CellNeighbourPair idx) -> std::set> { - 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> result; + + bool sameCell = cellI.index() == cellJ.index(); + + for (auto i : cellI.atoms() | castToStructure) + for (auto j : cellJ.atoms() | castToStructure) + { + // Don't bond atoms to themselves + if (sameCell && i->index() == j->index()) + 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({i, j}); + } + + 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> empty; + std::set> empty; // In parallel, construct the list of the bonds that need to be added + auto pairs = cells.getCellNeighbourPairsWithSelf(); auto results = std::transform_reduce(ParallelPolicies::par_unseq, pairs.begin(), pairs.end(), empty, joinBonds, validBond); // Add the bonds serially From cd1777e68216f94b06343289f05e98a18745e428 Mon Sep 17 00:00:00 2001 From: Adam Washington Date: Thu, 24 Sep 2026 14:08:08 +0100 Subject: [PATCH 2/4] Templated casting view --- src/kernels/energy.cpp | 11 ++++------- src/kernels/force.cpp | 12 +++++------- src/nodes/calculateBonding.cpp | 13 +++---------- src/templates/algorithms.h | 15 +++++++++++++++ 4 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/kernels/energy.cpp b/src/kernels/energy.cpp index 80e5be289a..f3b3172df1 100644 --- a/src/kernels/energy.cpp +++ b/src/kernels/energy.cpp @@ -13,9 +13,6 @@ #include #include -auto castToConstConfiguration = - std::ranges::views::transform([](const AtomBase *x) { return static_cast(x); }); - EnergyKernel::EnergyKernel(const Configuration *cfg, const PotentialMap &potentialMap) : GeometryKernel(cfg->box(), potentialMap), configuration_(cfg) { @@ -94,13 +91,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell ¢ // Loop over central cell atoms if (applyMim) { - for (auto ii : centralAtoms | castToConstConfiguration) + for (auto ii : centralAtoms | castView()) { auto molI = ii->molecule(); auto &rI = ii->r(); // Straight loop over other cell atoms - for (const auto jj : otherAtoms | castToConstConfiguration) + for (const auto jj : otherAtoms | castView()) { // Calculate rSquared distance between atoms, and check it against the stored cutoff distance auto rSq = box_.minimumDistanceSquared(rI, jj->r()); @@ -126,13 +123,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell ¢ } else { - for (auto ii : centralAtoms | castToConstConfiguration) + for (auto ii : centralAtoms | castView()) { auto &molI = ii->molecule(); auto &rI = ii->r(); // Straight loop over other cell atoms - for (const auto jj : otherAtoms | castToConstConfiguration) + for (const auto jj : otherAtoms | castView()) { // Calculate rSquared distance between atoms, and check it against the stored cutoff distance auto rSq = (rI - jj->r()).magnitudeSq(); diff --git a/src/kernels/force.cpp b/src/kernels/force.cpp index 6640fdb28b..0217f6ee25 100644 --- a/src/kernels/force.cpp +++ b/src/kernels/force.cpp @@ -11,8 +11,6 @@ #include "templates/algorithms.h" #include -auto castToConfiguration = std::ranges::views::transform([](AtomBase *x) { return static_cast(x); }); - ForceKernel::ForceKernel(const Configuration *cfg, const PotentialMap &potentialMap) : GeometryKernel(cfg->box(), potentialMap), configuration_(cfg) { @@ -98,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 | castToConfiguration) + for (const auto &i : centralAtoms | castView()) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms | castToConfiguration) + for (auto *j : otherAtoms | castView()) if (molI != j->molecule()) forcesWithMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f); } } else { - for (const auto &i : centralAtoms | castToConfiguration) + for (const auto &i : centralAtoms | castView()) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms | castToConfiguration) + for (auto *j : otherAtoms | castView()) if (molI != j->molecule()) forcesWithoutMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f); } @@ -162,7 +160,7 @@ void ForceKernel::totalForces(std::vector &ppForceVector, std::vectoratoms() | castToConfiguration, + dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms() | castView(), [&](int indexI, const auto &i, int indexJ, const auto &j) { if (indexI == indexJ) diff --git a/src/nodes/calculateBonding.cpp b/src/nodes/calculateBonding.cpp index 451ccf8ecd..a1142f130b 100644 --- a/src/nodes/calculateBonding.cpp +++ b/src/nodes/calculateBonding.cpp @@ -58,23 +58,16 @@ void CalculateBondingNode::calculate(Structure &structure, double tolerance, boo cells.generate(box, box.inscribedSphereRadius() * 2.0 / 3.0); // Populate cells - // - // As an ugly hack for the moment, we'll shove the StructureAtom - // pointers into the ConfigurationAtom points that the cell - // expects. for (auto &atom : structure.atoms()) cells.cell(atom->r())->addAtom(atom.get()); - // We'll eventually need a view to transform these pointers back - auto castToStructure = std::ranges::views::transform([](const auto x) { return (StructureAtom *)x; }); - // 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 = [&castToStructure, &structure, &box, tolerance, preventMetallic, + auto validBond = [&structure, &box, tolerance, preventMetallic, clearBefore](const CellNeighbourPair idx) -> std::set> { auto cellI = idx.cell; @@ -84,8 +77,8 @@ void CalculateBondingNode::calculate(Structure &structure, double tolerance, boo bool sameCell = cellI.index() == cellJ.index(); - for (auto i : cellI.atoms() | castToStructure) - for (auto j : cellJ.atoms() | castToStructure) + for (auto i : cellI.atoms() | castView()) + for (auto j : cellJ.atoms() | castView()) { // Don't bond atoms to themselves if (sameCell && i->index() == j->index()) diff --git a/src/templates/algorithms.h b/src/templates/algorithms.h index c66e2492c6..7ad726ec62 100644 --- a/src/templates/algorithms.h +++ b/src/templates/algorithms.h @@ -304,3 +304,18 @@ template class U> inline constexpr bool is_instanc template