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..f3b3172df1 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" @@ -46,7 +47,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 +55,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 +91,13 @@ Kernel::PairPotentialEnergyValue EnergyKernel::cellToCellEnergy(const Cell ¢ // Loop over central cell atoms if (applyMim) { - for (auto &ii : centralAtoms) + for (auto ii : centralAtoms | castView()) { 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()) { // 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 ¢ } else { - for (auto &ii : centralAtoms) + for (auto ii : centralAtoms | castView()) { 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()) { // 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(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(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..0217f6ee25 100644 --- a/src/kernels/force.cpp +++ b/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()) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms) + for (auto *j : otherAtoms | castView()) 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()) { molI = i->molecule(); auto indexI = molI->globalAtomIndex(i); - for (auto *j : otherAtoms) + for (auto *j : otherAtoms | castView()) if (molI != j->molecule()) forcesWithoutMim(*i, indexI, *j, j->molecule()->globalAtomIndex(j), f); } @@ -160,7 +160,7 @@ void ForceKernel::totalForces(std::vector &ppForceVector, std::vectoratoms(), + 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 05aea2ddb1..393176010d 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,13 @@ 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 @@ -62,47 +67,51 @@ 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 = [&box, tolerance, preventMetallic](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; + + for (auto i : cellI.atoms() | castView()) + for (auto j : cellJ.atoms() | castView()) + { + // 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> 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>(), joinBonds, validBond); // Add the bonds serially for (auto [i, j] : results) 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