diff --git a/src/classes/kVector.cpp b/src/classes/kVector.cpp
index 3fe2013604..b065e01523 100644
--- a/src/classes/kVector.cpp
+++ b/src/classes/kVector.cpp
@@ -72,7 +72,7 @@ void KVector::calculateIntensities(std::vector<BraggReflection> &reflections)
     auto &braggReflection = reflections[braggReflectionIndex_];
     braggReflection.addKVectors(halfSphereNorm);
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, cosTerms_.size(),
+        ParallelPolicies::par, cosTerms_.size(),
         [&](auto i, auto j)
         { braggReflection.addIntensity(i, j, (cosTerms_[i] * cosTerms_[j] + sinTerms_[i] * sinTerms_[j]) * halfSphereNorm); });
 }
diff --git a/src/classes/neutronWeights.cpp b/src/classes/neutronWeights.cpp
index b0e79b7d31..adb65c1886 100644
--- a/src/classes/neutronWeights.cpp
+++ b/src/classes/neutronWeights.cpp
@@ -118,7 +118,7 @@ void NeutronWeights::calculateWeightingMatrices()
     double ci, cj, bi, bj;
 
     // Determine atomic concentration products, bound coherent products, and full scattering weights
-    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_.begin(), atomTypes_.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_,
                             [&](int typeI, const AtomTypeData &atd1, int typeJ, const AtomTypeData &atd2)
                             {
                                 ci = atd1.fraction();
@@ -166,7 +166,7 @@ void NeutronWeights::calculateWeightingMatrices()
             auto weight = speciesWeight * isoWeight.weight();
             const auto *tope = isoWeight.isotopologue();
 
-            dissolve::for_each_pair(ParallelPolicies::seq, speciesAtomTypes.begin(), speciesAtomTypes.end(),
+            dissolve::for_each_pair(ParallelPolicies::seq, speciesAtomTypes,
                                     [&](int spTypeI, const AtomTypeData &atd1, int spTypeJ, const AtomTypeData &atd2)
                                     {
                                         // First, check that both of atom types used in the species are present in the weights
@@ -200,7 +200,7 @@ void NeutronWeights::calculateWeightingMatrices()
     }
 
     // Normalise the boundWeights_ array, and multiply by atomic concentrations and Kronecker delta
-    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_.begin(), atomTypes_.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_,
                             [&](int typeI, const AtomTypeData &atd1, int typeJ, const AtomTypeData &atd2)
                             {
                                 // Skip this pair if there are no such intramolecular interactions
diff --git a/src/classes/partialSet.cpp b/src/classes/partialSet.cpp
index c003321365..53b45a617d 100644
--- a/src/classes/partialSet.cpp
+++ b/src/classes/partialSet.cpp
@@ -59,7 +59,7 @@ bool PartialSet::setUpPartials(const AtomTypeMix &atomTypeMix, bool half)
 
     // Set up array matrices for partials
     dissolve::for_each_pair(
-        ParallelPolicies::par, atomTypeMix_.begin(), atomTypeMix_.end(),
+        ParallelPolicies::par, atomTypeMix_,
         [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
         {
             partials_[{n, m}].setTag(std::format("{}-{}//Full", at1.atomTypeName(), at2.atomTypeName()));
@@ -89,7 +89,7 @@ void PartialSet::setUpHistograms(double rdfRange, double binWidth)
     unboundHistograms_.initialise(nTypes, nTypes, half_);
 
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, nTypes,
+        ParallelPolicies::par, nTypes,
         [&](int i, int j)
         {
             fullHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
@@ -113,7 +113,7 @@ void PartialSet::reset()
 
     // Zero partials
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, atomTypeMix_.nItems(),
+        ParallelPolicies::par, atomTypeMix_.nItems(),
         [&](int i, int j)
         {
             std::ranges::fill(partials_[{i, j}].values(), 0.0);
@@ -186,7 +186,7 @@ void PartialSet::formTotals(bool applyConcentrationWeights)
     std::fill(total_.values().begin(), total_.values().end(), 0.0);
 
     dissolve::for_each_pair(
-        ParallelPolicies::seq, atomTypeMix_.begin(), atomTypeMix_.end(),
+        ParallelPolicies::seq, atomTypeMix_,
         [&](int typeI, const AtomTypeData &at1, int typeJ, const AtomTypeData &at2)
         {
             // Set weighting factor if requested
@@ -229,7 +229,7 @@ void PartialSet::formTRTotals(NeutronWeights weights)
     std::fill(total_.values().begin(), total_.values().end(), 0.0);
 
     dissolve::for_each_pair(
-        ParallelPolicies::seq, atomTypeMix_.begin(), atomTypeMix_.end(),
+        ParallelPolicies::seq, atomTypeMix_,
         [&](int typeI, const AtomTypeData &at1, int typeJ, const AtomTypeData &at2)
         {
             // Set weighting factor if requested
@@ -280,7 +280,7 @@ bool PartialSet::save(std::string_view prefix, std::string_view tag, std::string
 
     // Write partials
     for_each_pair_early(
-        atomTypeMix_.begin(), atomTypeMix_.end(),
+        atomTypeMix_,
         [&](int typeI, const AtomTypeData &at1, int typeJ, const AtomTypeData &at2) -> EarlyReturn<bool>
         {
             // Open file and check that we're OK to proceed writing to it
@@ -327,7 +327,7 @@ bool PartialSet::save(std::string_view prefix, std::string_view tag, std::string
 void PartialSet::adjust(double delta)
 {
     dissolve::for_each_pair(
-        ParallelPolicies::par, atomTypeMix_.begin(), atomTypeMix_.end(),
+        ParallelPolicies::par, atomTypeMix_,
         [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
         {
             partials_[{n, m}] += delta;
@@ -345,7 +345,7 @@ void PartialSet::adjust(double delta)
 void PartialSet::formPartials(double boxVolume)
 {
     dissolve::for_each_pair(
-        ParallelPolicies::seq, atomTypeMix_.begin(), atomTypeMix_.end(),
+        ParallelPolicies::seq, atomTypeMix_,
         [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
         {
             // Calculate RDFs from histogram data
@@ -450,7 +450,7 @@ void PartialSet::operator+=(const PartialSet &source)
     // Loop over partials in source set
     const auto &types = source.atomTypeMix();
     dissolve::for_each_pair(
-        ParallelPolicies::seq, types.begin(), types.end(),
+        ParallelPolicies::seq, types,
         [&](int typeI, const AtomTypeData &atd1, int typeJ, const AtomTypeData &atd2)
         {
             auto optPairIndex = atomTypeMix_.indexOf(atd1.atomType(), atd2.atomType());
@@ -482,7 +482,7 @@ void PartialSet::operator-=(const double delta) { adjust(-delta); }
 void PartialSet::operator*=(const double factor)
 {
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, atomTypeMix_.nItems(),
+        ParallelPolicies::par, atomTypeMix_.nItems(),
         [&](auto n, auto m)
         {
             partials_[{n, m}] *= factor;
@@ -557,7 +557,7 @@ bool PartialSet::deserialise(LineParser &parser, const CoreData &coreData)
     emptyBoundPartials_ = false;
 
     for_each_pair_early(
-        0, nTypes,
+        nTypes,
         [&](int typeI, int typeJ) -> EarlyReturn<bool>
         {
             auto &part = partials_[{typeI, typeJ}];
@@ -657,7 +657,7 @@ bool PartialSet::serialise(LineParser &parser) const
 
     // Write individual Data1D
     auto success = for_each_pair_early(
-        0, nTypes,
+        nTypes,
         [&](int typeI, int typeJ) -> EarlyReturn<bool>
         {
             const auto &part = partials_[{typeI, typeJ}];
diff --git a/src/classes/partialSetAccumulator.cpp b/src/classes/partialSetAccumulator.cpp
index 8c6902a4cb..6b41dfad02 100644
--- a/src/classes/partialSetAccumulator.cpp
+++ b/src/classes/partialSetAccumulator.cpp
@@ -23,7 +23,7 @@ void PartialSetAccumulator::operator+=(const PartialSet &source)
         total_.clear();
 
         // Copy tags (for retrieval and sanity checking purposes)
-        dissolve::for_each_pair(ParallelPolicies::par, 0, n,
+        dissolve::for_each_pair(ParallelPolicies::par, n,
                                 [&](auto i, auto j)
                                 {
                                     partials_[{i, j}].setTag(source.partial(i, j).tag());
@@ -37,7 +37,7 @@ void PartialSetAccumulator::operator+=(const PartialSet &source)
 
     // Accumulate the data, ensuring tags are identical - if not, we really don't want to be blindly accumulating
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, n,
+        ParallelPolicies::par, n,
         [&](auto i, auto j)
         {
             // Full partials
diff --git a/src/classes/scatteringMatrix.cpp b/src/classes/scatteringMatrix.cpp
index 6d97e2c2f1..6758071bf8 100644
--- a/src/classes/scatteringMatrix.cpp
+++ b/src/classes/scatteringMatrix.cpp
@@ -355,7 +355,7 @@ void ScatteringMatrix::initialise(const AtomTypeMix &typeMix, Array2D<Data1D> &e
     // Copy atom types and construct pairs
     atomTypes_.resize(typeMix.nItems());
     std::transform(typeMix.begin(), typeMix.end(), atomTypes_.begin(), [](const auto &atd) { return atd.atomType(); });
-    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_.begin(), atomTypes_.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes_,
                             [this](int i, auto &at1, int j, auto &at2) { typePairs_.emplace_back(at1, at2); });
 
     // Create partials array
diff --git a/src/classes/xRayWeights.cpp b/src/classes/xRayWeights.cpp
index a3590144e9..15f1f30784 100644
--- a/src/classes/xRayWeights.cpp
+++ b/src/classes/xRayWeights.cpp
@@ -135,7 +135,7 @@ void XRayWeights::setUpMatrices()
     preFactors_.initialise(atomTypeMix_.nItems(), atomTypeMix_.nItems(), true);
 
     // Determine atomic concentration products and full pre-factor
-    dissolve::for_each_pair(ParallelPolicies::seq, atomTypeMix_.begin(), atomTypeMix_.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypeMix_,
                             [&](int typeI, const AtomTypeData &atd1, int typeJ, const AtomTypeData &atd2)
                             {
                                 double ci = atd1.fraction();
diff --git a/src/kernels/energy.cpp b/src/kernels/energy.cpp
index 5ae171fda5..6f739c53b4 100644
--- a/src/kernels/energy.cpp
+++ b/src/kernels/energy.cpp
@@ -282,7 +282,7 @@ PairPotentialEnergyValue EnergyKernel::pairPotentialEnergy(const Molecule &mol,
     if (includeIntraMolecular)
     {
         auto intra = 0.0;
-        dissolve::for_each_pair(ParallelPolicies::seq, 0, mol.nAtoms(),
+        dissolve::for_each_pair(ParallelPolicies::seq, mol.nAtoms(),
                                 [&](int i, int j)
                                 {
                                     if (i == j)
diff --git a/src/kernels/force.cpp b/src/kernels/force.cpp
index aeff612f65..4db915055c 100644
--- a/src/kernels/force.cpp
+++ b/src/kernels/force.cpp
@@ -159,7 +159,7 @@ void ForceKernel::totalForces(ForceVector &fUnbound, ForceVector &fBound, Flags<
             auto &fLocal = combinableUnbound.local();
 
             // Interatomic interactions between atoms in this cell, excluding those within the same molecule
-            dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms().begin(), cellI->atoms().end(),
+            dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms(),
                                     [&](int indexI, const auto &i, int indexJ, const auto &j)
                                     {
                                         if (indexI == indexJ)
@@ -199,7 +199,7 @@ void ForceKernel::totalForces(ForceVector &fUnbound, ForceVector &fBound, Flags<
 
             // Pair potential interactions between atoms within the molecule
             if (!flags.isSet(ExcludeIntraMolecularPairPotential))
-                dissolve::for_each_pair(ParallelPolicies::seq, mol->atoms().begin(), mol->atoms().end(),
+                dissolve::for_each_pair(ParallelPolicies::seq, mol->atoms(),
                                         [&](int indexI, const auto &i, int indexJ, const auto &j)
                                         {
                                             if (indexI == indexJ)
diff --git a/src/main/pairPotentials.cpp b/src/main/pairPotentials.cpp
index 06bf397579..b51f2e21e6 100644
--- a/src/main/pairPotentials.cpp
+++ b/src/main/pairPotentials.cpp
@@ -119,7 +119,7 @@ bool Dissolve::updatePairPotentials(std::optional<bool> useCombinationRulesHint)
                           pairPotentials_.end());
 
     // Second step - add or update tabulated pair potentials defined by the parameters and form of the associated atom types
-    if (!for_each_pair_early(coreData_.atomTypes().begin(), coreData_.atomTypes().end(),
+    if (!for_each_pair_early(coreData_.atomTypes(),
                              [&](int typeI, const auto &at1, int typeJ, const auto &at2) -> EarlyReturn<bool>
                              {
                                  // Try to locate existing pair potential between these atom types
diff --git a/src/modules/bragg/functions.cpp b/src/modules/bragg/functions.cpp
index c985a9e68b..2f9542cc63 100644
--- a/src/modules/bragg/functions.cpp
+++ b/src/modules/bragg/functions.cpp
@@ -347,7 +347,7 @@ bool BraggModule::formReflectionFunctions(GenericList &moduleData, Configuration
     double qCentre;
     int bin;
     auto &types = cfg->atomTypePopulations();
-    dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, types,
                             [&](int typeI, auto &atd1, int typeJ, auto &atd2)
                             {
                                 // Retrieve partial container and make sure its tag is set
diff --git a/src/modules/bragg/process.cpp b/src/modules/bragg/process.cpp
index 42e815b9c9..3c0a527744 100644
--- a/src/modules/bragg/process.cpp
+++ b/src/modules/bragg/process.cpp
@@ -89,7 +89,7 @@ Module::ExecutionResult BraggModule::process(Dissolve &dissolve)
         // Save intensity data
         auto &types = targetConfiguration_->atomTypePopulations();
         auto success = for_each_pair_early(
-            types.begin(), types.end(),
+            types,
             [&](int i, const AtomTypeData &atd1, int j, const AtomTypeData &atd2) -> EarlyReturn<bool>
             {
                 LineParser intensityParser;
diff --git a/src/modules/epsr/functions.cpp b/src/modules/epsr/functions.cpp
index 0e858625f3..9f3a9b2135 100644
--- a/src/modules/epsr/functions.cpp
+++ b/src/modules/epsr/functions.cpp
@@ -74,7 +74,7 @@ bool EPSRModule::generateEmpiricalPotentials(Dissolve &dissolve, double averaged
     Array2D<std::vector<double>> &coefficients = potentialCoefficients(dissolve.processingModuleData(), nAtomTypes, ncoeffp);
 
     auto result = for_each_pair_early(
-        atomTypes.begin(), atomTypes.end(),
+        atomTypes,
         [&](int i, auto at1, int j, auto at2) -> EarlyReturn<bool>
         {
             auto &potCoeff = coefficients[{i, j}];
diff --git a/src/modules/epsr/gui/epsrWidgetFuncs.cpp b/src/modules/epsr/gui/epsrWidgetFuncs.cpp
index 1c2680cc0c..8459a44090 100644
--- a/src/modules/epsr/gui/epsrWidgetFuncs.cpp
+++ b/src/modules/epsr/gui/epsrWidgetFuncs.cpp
@@ -109,7 +109,7 @@ void EPSRModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFlags> &up
         {
             // Add experimentally-determined partial S(Q)
             dissolve::for_each_pair(
-                ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+                ParallelPolicies::seq, atomTypes,
                 [&](int typeI, const auto &at1, int typeJ, const auto &at2)
                 {
                     const std::string id = std::format("{}-{}", at1->name(), at2->name());
@@ -152,7 +152,7 @@ void EPSRModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFlags> &up
 
             // Add experimentally-determined partial g(r)
             dissolve::for_each_pair(
-                ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+                ParallelPolicies::seq, atomTypes,
                 [&](int typeI, const auto &at1, int typeJ, const auto &at2)
                 {
                     const std::string id = std::format("{}-{}", at1->name(), at2->name());
@@ -183,7 +183,7 @@ void EPSRModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFlags> &up
         else if (ui_.PotentialsButton->isChecked())
         {
             // Add on additional potentials
-            dissolve::for_each_pair(ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+            dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
                                     [&](int typeI, const auto &at1, int typeJ, const auto &at2)
                                     {
                                         const std::string id = std::format("{}-{}", at1->name(), at2->name());
diff --git a/src/modules/epsr/process.cpp b/src/modules/epsr/process.cpp
index 3439b318f4..a4f5dfa84b 100644
--- a/src/modules/epsr/process.cpp
+++ b/src/modules/epsr/process.cpp
@@ -217,7 +217,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
     // Create storage for our summed UnweightedSQ
     auto &calculatedUnweightedSQ = moduleData.realise<Array2D<Data1D>>("UnweightedSQ", name_, GenericItem::InRestartFileFlag);
     calculatedUnweightedSQ.initialise(nAtomTypes, nAtomTypes, true);
-    dissolve::for_each_pair(ParallelPolicies::par, atomTypes.begin(), atomTypes.end(),
+    dissolve::for_each_pair(ParallelPolicies::par, atomTypes,
                             [&](int i, auto at1, int j, auto at2) {
                                 calculatedUnweightedSQ[{i, j}].setTag(std::format("{}-{}", at1->name(), at2->name()));
                             });
@@ -498,7 +498,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
 
         // Add the unweighted from this target to our combined, unweighted S(Q) data
         auto &types = scatteringMatrix_.atomTypes();
-        dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, types,
                                 [&](int i, const auto &at1, int j, const auto &at2)
                                 {
                                     auto pairIndex = scatteringMatrix_.pairIndexOf(at1, at2);
@@ -549,7 +549,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
 
     // Add a contribution from each interatomic partial S(Q), weighted according to the feedback factor
     auto success = for_each_pair_early(
-        atomTypes.begin(), atomTypes.end(),
+        atomTypes,
         [&](int i, auto at1, int j, auto at2) -> EarlyReturn<bool>
         {
             // Copy and rename the data for clarity
@@ -608,7 +608,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
      */
     auto &estimatedGR = moduleData.realise<Array2D<Data1D>>("EstimatedGR", name_, GenericItem::InRestartFileFlag);
     estimatedGR.initialise(nAtomTypes, nAtomTypes, true);
-    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
                             [&](int i, auto at1, int j, auto at2)
                             {
                                 auto &expGR = estimatedGR[{i, j}];
@@ -635,7 +635,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
 
         // Loop over pair potentials and retrieve the inverse weight from the scattering matrix
         dissolve::for_each_pair(
-            ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+            ParallelPolicies::seq, atomTypes,
             [&](int i, auto at1, int j, auto at2)
             {
                 auto weight = scatteringMatrix_.qZeroMatrixInverse()[{scatteringMatrix_.columnIndex(at1, at2), dataIndex}];
@@ -673,7 +673,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
     {
         // Sum fluctuation coefficients in to the potential coefficients
         auto &coefficients = potentialCoefficients(moduleData, nAtomTypes, ncoeffp);
-        dissolve::for_each_pair(ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
                                 [&](int i, auto at1, int j, auto at2)
                                 {
                                     auto &potCoeff = coefficients[{i, j}];
@@ -741,7 +741,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
     // Save data?
     if (saveEmpiricalPotentials_)
     {
-        if (!for_each_pair_early(atomTypes.begin(), atomTypes.end(),
+        if (!for_each_pair_early(atomTypes,
                                  [&](int i, auto at1, int j, auto at2) -> EarlyReturn<bool>
                                  {
                                      // Grab pointer to the relevant pair potential
@@ -760,7 +760,7 @@ Module::ExecutionResult EPSRModule::process(Dissolve &dissolve)
     {
         auto &coefficients = potentialCoefficients(moduleData, nAtomTypes, ncoeffp);
 
-        if (!for_each_pair_early(atomTypes.begin(), atomTypes.end(),
+        if (!for_each_pair_early(atomTypes,
                                  [&](int i, auto at1, int j, auto at2) -> EarlyReturn<bool>
                                  {
                                      // Grab reference to coefficients
diff --git a/src/modules/epsrManager/gui/epsrManagerFuncs.cpp b/src/modules/epsrManager/gui/epsrManagerFuncs.cpp
index 80542c5525..a915308155 100644
--- a/src/modules/epsrManager/gui/epsrManagerFuncs.cpp
+++ b/src/modules/epsrManager/gui/epsrManagerFuncs.cpp
@@ -53,7 +53,7 @@ void EPSRManagerModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFla
         if (ui_.PotentialsButton->isChecked())
         {
             // Add on additional potentials
-            dissolve::for_each_pair(ParallelPolicies::seq, atomTypes.begin(), atomTypes.end(),
+            dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
                                     [&](int typeI, const auto &at1, int typeJ, const auto &at2)
                                     {
                                         const std::string id = std::format("{}-{}", at1->name(), at2->name());
diff --git a/src/modules/forces/functions.cpp b/src/modules/forces/functions.cpp
index f29a2de4c8..7aeb360810 100644
--- a/src/modules/forces/functions.cpp
+++ b/src/modules/forces/functions.cpp
@@ -115,9 +115,9 @@ void ForcesModule::totalForces(const Species *sp, const PotentialMap &potentialM
     if (calculationType == ForceCalculationType::Full || calculationType == ForceCalculationType::PairPotentialOnly)
     {
         if (sp->nAtoms() < 100)
-            dissolve::for_each_pair(ParallelPolicies::seq, sp->atoms().begin(), sp->atoms().end(), pairwiseForceOperator);
+            dissolve::for_each_pair(ParallelPolicies::seq, std::span(sp->atoms().begin(), sp->nAtoms()), pairwiseForceOperator);
         else
-            dissolve::for_each_pair(ParallelPolicies::par, sp->atoms().begin(), sp->atoms().end(), pairwiseForceOperator);
+            dissolve::for_each_pair(ParallelPolicies::par, std::span(sp->atoms().begin(), sp->nAtoms()), pairwiseForceOperator);
         combinableUnbound.finalize();
     }
 
diff --git a/src/modules/gr/functions.cpp b/src/modules/gr/functions.cpp
index 75e06091e7..29f9850b1a 100644
--- a/src/modules/gr/functions.cpp
+++ b/src/modules/gr/functions.cpp
@@ -1,6 +1,8 @@
 // SPDX-License-Identifier: GPL-3.0-or-later
 // Copyright (c) 2025 Team Dissolve and contributors
 
+#include <ranges>
+
 #include "classes/atom.h"
 #include "classes/atomType.h"
 #include "classes/box.h"
@@ -45,7 +47,7 @@ bool GRModule::calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet)
     const auto *box = cfg->box();
 
     dissolve::for_each_pair(
-        ParallelPolicies::seq, cfg->atoms().begin(), cfg->atoms().end(),
+        ParallelPolicies::seq, cfg->atoms(),
         [box, &partialSet](auto i, auto &ii, auto j, auto &jj)
         {
             if (&ii != &jj)
@@ -377,7 +379,7 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo
     {
         const auto &atoms = mol->atoms();
 
-        dissolve::for_each_pair(ParallelPolicies::seq, atoms.begin(), atoms.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, atoms,
                                 [box, &originalgr](int index, auto &i, int jndex, auto &j)
                                 {
                                     // Ignore atom on itself
@@ -407,7 +409,7 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo
 
     timer.start();
     auto success =
-        for_each_pair_early(0, originalgr.nAtomTypes(),
+        for_each_pair_early(originalgr.nAtomTypes(),
                             [&originalgr](auto typeI, auto typeJ) -> EarlyReturn<bool>
                             {
                                 // Create unbound histogram from total and bound data
@@ -467,19 +469,19 @@ bool GRModule::calculateUnweightedGR(Configuration *cfg, const PartialSet &origi
 
     // Broaden the bound partials according to the supplied PairBroadeningFunction
     auto &types = unweightedgr.atomTypeMix();
-    dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, types,
                             [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                             { Filters::convolve(unweightedgr.boundPartial(i, j), intraBroadening, true, true); });
 
     // Add broadened bound partials back in to full partials
-    dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, types,
                             [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                             { unweightedgr.partial(i, j) += unweightedgr.boundPartial(i, j); });
 
     // Apply smoothing if requested
     if (smoothing > 0)
     {
-        dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, types,
                                 [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                                 {
                                     Filters::movingAverage(unweightedgr.partial(i, j), smoothing);
@@ -568,7 +570,7 @@ bool GRModule::testReferencePartials(PartialSet &setA, PartialSet &setB, double
 
     for_each_pair_early(
 
-        atomTypes.begin(), atomTypes.end(),
+        atomTypes,
         [&](int n, const AtomTypeData &typeI, int m, const AtomTypeData &typeJ) -> EarlyReturn<bool>
         {
             // Full partial
diff --git a/src/modules/sq/functions.cpp b/src/modules/sq/functions.cpp
index 002ab3463e..e723e7ea00 100644
--- a/src/modules/sq/functions.cpp
+++ b/src/modules/sq/functions.cpp
@@ -10,6 +10,8 @@
 #include "modules/sq/sq.h"
 #include "templates/algorithms.h"
 
+#include <ranges>
+
 /*
  * Public Functions
  */
@@ -28,7 +30,7 @@ bool SQModule::calculateUnweightedSQ(const PartialSet &unweightedgr, PartialSet
     Timer timer;
     timer.start();
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, unweightedgr.nAtomTypes(),
+        ParallelPolicies::par, unweightedgr.nAtomTypes(),
         [&](int n, int m)
         {
             // Total partial
diff --git a/src/modules/sq/process.cpp b/src/modules/sq/process.cpp
index ada5f7bca9..59948da382 100644
--- a/src/modules/sq/process.cpp
+++ b/src/modules/sq/process.cpp
@@ -12,6 +12,8 @@
 #include "modules/sq/sq.h"
 #include "templates/algorithms.h"
 
+#include <ranges>
+
 // Set target data
 void SQModule::setTargets(const std::vector<std::unique_ptr<Configuration>> &configurations,
                           const std::map<ModuleTypes::ModuleType, std::vector<const Module *>> &moduleMap)
@@ -135,7 +137,7 @@ Module::ExecutionResult SQModule::process(Dissolve &dissolve)
 
         // For each partial in our S(Q) array, calculate the broadened Bragg function and blend it
         auto success = for_each_pair_early(
-            unweightedsq.atomTypeMix().begin(), unweightedsq.atomTypeMix().end(),
+            unweightedsq.atomTypeMix(),
             [&](auto i, auto &at1, auto j, auto &at2) -> EarlyReturn<bool>
             {
                 // Locate the corresponding Bragg intensities for this atom type pair
@@ -170,7 +172,7 @@ Module::ExecutionResult SQModule::process(Dissolve &dissolve)
                            [v0](auto &val) { return val * 2.0 * pow(M_PI, 2) / v0; });
 
         // Remove self-scattering level from partials between the same atom type and remove normalisation from atomic fractions
-        dissolve::for_each_pair(ParallelPolicies::par, unweightedsq.atomTypeMix().begin(), unweightedsq.atomTypeMix().end(),
+        dissolve::for_each_pair(ParallelPolicies::par, unweightedsq.atomTypeMix(),
                                 [&braggPartials](auto i, auto &atd1, auto j, auto &atd2)
                                 {
                                     // Subtract self-scattering level if types are equivalent
@@ -182,7 +184,7 @@ Module::ExecutionResult SQModule::process(Dissolve &dissolve)
                                 });
 
         // Blend the bound/unbound and Bragg partials at the higher Q limit
-        dissolve::for_each_pair(ParallelPolicies::par, 0, unweightedsq.nAtomTypes(),
+        dissolve::for_each_pair(ParallelPolicies::par, unweightedsq.nAtomTypes(),
                                 [&](const int i, const int j)
                                 {
                                     // Note: Intramolecular broadening will not be applied to bound terms within the
diff --git a/src/modules/tr/gui/trWidgetFuncs.cpp b/src/modules/tr/gui/trWidgetFuncs.cpp
index 552c109848..36f8950aed 100644
--- a/src/modules/tr/gui/trWidgetFuncs.cpp
+++ b/src/modules/tr/gui/trWidgetFuncs.cpp
@@ -61,7 +61,7 @@ void TRModuleWidget::createPartialSetRenderables(std::string_view targetPrefix)
 
     // Set up array matrices for partials
     dissolve::for_each_pair(
-        ParallelPolicies::seq, ps.atomTypeMix().begin(), ps.atomTypeMix().end(),
+        ParallelPolicies::seq, ps.atomTypeMix(),
         [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
         {
             const std::string id = std::format("{}-{}", at1.atomTypeName(), at2.atomTypeName());
diff --git a/src/modules/tr/process.cpp b/src/modules/tr/process.cpp
index aa32da257c..aef7d455c0 100644
--- a/src/modules/tr/process.cpp
+++ b/src/modules/tr/process.cpp
@@ -11,6 +11,8 @@
 #include "modules/sq/sq.h"
 #include "modules/tr/tr.h"
 
+#include <ranges>
+
 // Run main processing
 Module::ExecutionResult TRModule::process(Dissolve &dissolve)
 {
@@ -86,7 +88,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve)
                          WindowFunction::forms().keyword(repWindowFunction_));
     // FT unweightedSQ to unweightedGR to get better representation of calculations
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, unweightedSQ.nAtomTypes(),
+        ParallelPolicies::par, unweightedSQ.nAtomTypes(),
         [&](int n, int m)
         {
             // Total partial
@@ -127,7 +129,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve)
         representativeTR.setUpPartials(representativeGR.atomTypeMix(), false);
 
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, representativeGR.nAtomTypes(),
+        ParallelPolicies::par, representativeGR.nAtomTypes(),
         [&weights, &rho, &representativeGR, &representativeTR](const auto typeI, const auto typeJ)
         {
             double intraWeight = weights.intramolecularWeight(typeI, typeJ);
@@ -161,7 +163,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve)
 
     // Calculate weightedTR
     dissolve::for_each_pair(
-        ParallelPolicies::seq, 0, unweightedGR.nAtomTypes(),
+        ParallelPolicies::seq, unweightedGR.nAtomTypes(),
         [&weights, &rho, &unweightedGR, &weightedTR](const auto typeI, const auto typeJ)
         {
             double intraWeight = weights.intramolecularWeight(typeI, typeJ);
diff --git a/src/modules/xRaySQ/process.cpp b/src/modules/xRaySQ/process.cpp
index 93f5ca69ed..a87551c307 100644
--- a/src/modules/xRaySQ/process.cpp
+++ b/src/modules/xRaySQ/process.cpp
@@ -212,7 +212,7 @@ Module::ExecutionResult XRaySQModule::process(Dissolve &dissolve)
         return ExecutionResult::Failed;
     if (saveFormFactors_)
     {
-        auto result = for_each_pair_early(unweightedSQ.atomTypeMix().begin(), unweightedSQ.atomTypeMix().end(),
+        auto result = for_each_pair_early(unweightedSQ.atomTypeMix(),
                                           [&](int i, auto &at1, int j, auto &at2) -> EarlyReturn<bool>
                                           {
                                               if (i == j)
diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp
index 1d2d00af3d..05c17bcd38 100644
--- a/src/nodes/gr/helpers.cpp
+++ b/src/nodes/gr/helpers.cpp
@@ -43,7 +43,7 @@ bool GRNode::calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet)
     const auto *box = cfg->box();
 
     dissolve::for_each_pair(
-        ParallelPolicies::seq, cfg->atoms().begin(), cfg->atoms().end(),
+        ParallelPolicies::seq, cfg->atoms(),
         [box, &partialSet](auto i, auto &ii, auto j, auto &jj)
         {
             if (&ii != &jj)
@@ -391,7 +391,7 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par
     {
         const auto &atoms = mol->atoms();
 
-        dissolve::for_each_pair(ParallelPolicies::seq, atoms.begin(), atoms.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, atoms,
                                 [box, &originalgr](int index, auto &i, int jndex, auto &j)
                                 {
                                     // Ignore atom on itself
@@ -421,7 +421,7 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par
 
     timer.start();
     auto success =
-        for_each_pair_early(0, originalgr.nAtomTypes(),
+        for_each_pair_early(originalgr.nAtomTypes(),
                             [&originalgr](auto typeI, auto typeJ) -> EarlyReturn<bool>
                             {
                                 // Create unbound histogram from total and bound data
@@ -479,19 +479,19 @@ bool GRNode::calculateUnweightedGR(Configuration *cfg, const PartialSet &origina
 
     // Broaden the bound partials according to the supplied PairBroadeningFunction
     auto &types = unweightedgr.atomTypeMix();
-    dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, types,
                             [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                             { Filters::convolve(unweightedgr.boundPartial(i, j), intraBroadening, true, true); });
 
     // Add broadened bound partials back in to full partials
-    dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+    dissolve::for_each_pair(ParallelPolicies::seq, types,
                             [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                             { unweightedgr.partial(i, j) += unweightedgr.boundPartial(i, j); });
 
     // Apply smoothing if requested
     if (smoothing > 0)
     {
-        dissolve::for_each_pair(ParallelPolicies::seq, types.begin(), types.end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, types,
                                 [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ)
                                 {
                                     Filters::movingAverage(unweightedgr.partial(i, j), smoothing);
@@ -577,8 +577,7 @@ bool GRNode::testReferencePartials(PartialSet &setA, PartialSet &setB, double te
     auto atomTypes = setA.atomTypeMix();
 
     for_each_pair_early(
-
-        atomTypes.begin(), atomTypes.end(),
+        atomTypes,
         [&](int n, const AtomTypeData &typeI, int m, const AtomTypeData &typeJ) -> EarlyReturn<bool>
         {
             // Full partial
diff --git a/src/nodes/registry.cpp b/src/nodes/registry.cpp
index c7d40db36c..3791393f8a 100644
--- a/src/nodes/registry.cpp
+++ b/src/nodes/registry.cpp
@@ -44,7 +44,7 @@ void NodeRegistry::instantiateNodeProducers()
                   {"Configuration", makeDerivedNode<ConfigurationNode>()},
                   {"Derivative", makeDerivedNode<DerivativeNode>()},
                   {"DotProduct", makeDerivedNode<DotProductNode>()},
-                  {"Energy", makeDerivedNode<EnergyNode>()},
+                  // {"Energy", makeDerivedNode<EnergyNode>()},
                   {"GR", makeDerivedNode<GRNode>()},
                   {"Graph", makeDerivedNode<Graph>()},
                   {"GR", makeDerivedNode<GRNode>()},
@@ -52,7 +52,7 @@ void NodeRegistry::instantiateNodeProducers()
                   {"Integrator", makeDerivedNode<Integrator1DNode>()},
                   {"MD", makeDerivedNode<MDNode>()},
                   {"Multiply", makeDerivedNode<MultiplyNode>()},
-                  {"NeutronSQ", makeDerivedNode<NeutronSQNode>()},
+                  // {"NeutronSQ", makeDerivedNode<NeutronSQNode>()},
                   {"Number", makeDerivedNode<NumberNode>()},
                   {"SQ", makeDerivedNode<SQNode>()},
                   {"Subtract", makeDerivedNode<SubtractNode>()},
diff --git a/src/nodes/sq/helpers.cpp b/src/nodes/sq/helpers.cpp
index eaa98fa602..6f95527926 100644
--- a/src/nodes/sq/helpers.cpp
+++ b/src/nodes/sq/helpers.cpp
@@ -29,7 +29,7 @@ bool SQNode::calculateUnweightedSQ(const PartialSet &unweightedgr, PartialSet &u
     Timer timer;
     timer.start();
     dissolve::for_each_pair(
-        ParallelPolicies::par, 0, unweightedgr.nAtomTypes(),
+        ParallelPolicies::par, unweightedgr.nAtomTypes(),
         [&](int n, int m)
         {
             // Total partial
diff --git a/src/templates/algorithms.h b/src/templates/algorithms.h
index 92e82d3a34..42978fdbfc 100644
--- a/src/templates/algorithms.h
+++ b/src/templates/algorithms.h
@@ -9,6 +9,7 @@
 #include <format>
 #include <functional>
 #include <optional>
+#include <ranges>
 #include <sstream>
 #include <tuple>
 #include <utility>
@@ -35,21 +36,29 @@ template <typename T> class EarlyReturn
     std::optional<T> value() const { return value_; }
 };
 
+// A way to check if a lambda *only* takes two int parameters
+template <typename T>
+concept PairIndexLambda = requires(T lam, int x, int y) { lam(x, y); };
+
 // Perform an operation on every pair of elements in a container,
 // or the half-matrix only ([i,j] == [j,i])
 // Please note that this can *not* be transformed to use the
 // FullPairIterator, since it would prevent using `Break` to move to
 // the next loop iteration
-template <class Iter, class Lam>
-auto for_each_pair_early(Iter begin, Iter end, Lam lambda, bool half = true) -> decltype(lambda(0, *begin, 0, *end).value())
+template <std::ranges::range Range, class Lam>
+auto for_each_pair_early(Range range, Lam lambda, bool half = true) -> std::optional<bool>
 {
     int i = 0;
-    for (auto elem1 = begin; elem1 != end; ++elem1, ++i)
+    for (auto elem1 = range.begin(); elem1 != range.end(); ++elem1, ++i)
     {
         int j = half ? i : 0;
-        for (auto elem2 = half ? elem1 : begin; elem2 != end; ++elem2, ++j)
+        for (auto elem2 = half ? elem1 : range.begin(); elem2 != range.end(); ++elem2, ++j)
         {
-            auto result = lambda(i, *elem1, j, *elem2);
+            EarlyReturn<bool> result;
+            if constexpr (PairIndexLambda<Lam>)
+                result = lambda(i, j);
+            else
+                result = lambda(i, *elem1, j, *elem2);
             switch (result.type())
             {
                 case EarlyReturn<typename decltype(result)::inner>::Return:
@@ -64,29 +73,10 @@ auto for_each_pair_early(Iter begin, Iter end, Lam lambda, bool half = true) ->
     return std::nullopt;
 }
 
-// Perform an operation on every pair of elements in a range, or the half-matrix only ([i,j] == [j,i])
-// Please note that this can *not* be transformed to use the
-// FullPairIterator, since it would prevent using `Break` to move to
-// the next loop iteration
-template <class Lam>
-auto for_each_pair_early(int begin, int end, Lam lambda, bool half = true) -> decltype(lambda(0, 0).value())
+// Overload to avoid using iota everywhere
+template <class Lam> auto for_each_pair_early(int count, Lam lambda, bool half = true) -> std::optional<bool>
 {
-    for (auto i = begin; i < end; ++i)
-        for (auto j = half ? i : begin; j < end; ++j)
-        {
-            auto result = lambda(i, j);
-            switch (result.type())
-            {
-                case EarlyReturn<typename decltype(result)::inner>::Return:
-                    return result.value();
-                case EarlyReturn<typename decltype(result)::inner>::Break:
-                    break;
-                case EarlyReturn<typename decltype(result)::inner>::Continue:
-                    continue;
-            }
-        }
-
-    return std::nullopt;
+    return for_each_pair_early(std::views::iota(0, count), lambda, half);
 }
 
 template <typename... Args> class ZipIterator
@@ -183,8 +173,8 @@ T transform_reduce(ParallelPolicy policy, Iter begin, Iter end, T initialVal, Bi
     return std::transform_reduce(policy, begin, end, initialVal, binaryOp, unaryOp);
 }
 
-// Enabled if parallelpolicy is not a real execution policy, i.e. we haven't compiled with multithreading but attempted to set a
-// parallel policy
+// Enabled if parallelpolicy is not a real execution policy, i.e. we haven't compiled with multithreading but attempted to
+// set a parallel policy
 template <typename ParallelPolicy, class Iter, typename T, class UnaryOp, class BinaryOp,
           std::enable_if_t<std::is_same_v<ParallelPolicy, FakeParallelPolicy>, bool> = true>
 T transform_reduce(ParallelPolicy, Iter begin, Iter end, T initialVal, BinaryOp binaryOp, UnaryOp unaryOp)
@@ -207,8 +197,8 @@ void for_each(ParallelPolicy policy, Iter begin, Iter end, UnaryOp unaryOp)
     std::for_each(policy, begin, end, unaryOp);
 }
 
-// Enabled if parallelpolicy is not a real execution policy, i.e. we haven't compiled with multithreading but attempted to set a
-// parallel policy
+// Enabled if parallelpolicy is not a real execution policy, i.e. we haven't compiled with multithreading but attempted to
+// set a parallel policy
 template <typename ParallelPolicy, class Iter, class UnaryOp,
           std::enable_if_t<std::is_same_v<ParallelPolicy, FakeParallelPolicy>, bool> = true>
 void for_each(ParallelPolicy, Iter begin, Iter end, UnaryOp unaryOp)
@@ -217,26 +207,37 @@ void for_each(ParallelPolicy, Iter begin, Iter end, UnaryOp unaryOp)
 }
 
 // Perform an operation on every pair of elements in a contained, or the half-matrix only ([i,j] == [j,i])
-template <typename ParallelPolicy, class Iter, class Lam>
-void for_each_pair(ParallelPolicy policy, Iter begin, Iter end, Lam lambda, bool half = true)
+template <typename ParallelPolicy, std::ranges::range Range, class Lam>
+void for_each_pair(ParallelPolicy policy, Range range, Lam lambda, bool half = true)
 {
-    auto actions = [&lambda, &begin](const auto pair)
+    auto actions = [&lambda, &range](const auto pair)
     {
         auto &[i, j] = pair;
-        lambda(i, begin[i], j, begin[j]);
+        if constexpr (PairIndexLambda<Lam>)
+            lambda(i, j);
+        else
+            lambda(i, range.begin()[i], j, range.begin()[j]);
     };
     if (half)
     {
-        PairIterator start(end - begin), stop(end - begin, ((end - begin) * (end - begin + 1)) / 2);
+        PairIterator start(range.end() - range.begin()),
+            stop(range.end() - range.begin(), ((range.end() - range.begin()) * (range.end() - range.begin() + 1)) / 2);
         for_each(policy, start, stop, actions);
     }
     else
     {
-        FullPairIterator start(end - begin), stop(end - begin, (end - begin) * (end - begin));
+        FullPairIterator start(range.end() - range.begin()),
+            stop(range.end() - range.begin(), (range.end() - range.begin()) * (range.end() - range.begin()));
         for_each(policy, start, stop, actions);
     }
 }
 
+// Overload to avoid using iota everywhere
+template <typename ParallelPolicy, class Lam> void for_each_pair(ParallelPolicy policy, int count, Lam lambda, bool half = true)
+{
+    for_each_pair(policy, std::views::iota(0, count), lambda, half);
+}
+
 template <typename ParallelPolicy, class Iter, class Lam>
 void for_each_triplet(ParallelPolicy policy, Iter begin, Iter end, Lam lambda)
 {
@@ -247,27 +248,6 @@ void for_each_triplet(ParallelPolicy policy, Iter begin, Iter end, Lam lambda)
                  lambda(triplet, x, y, z);
              });
 }
-
-// Perform an operation on every pair of elements in a range, or the half-matrix only ([i,j] == [j,i])
-template <typename ParallelPolicy, class Lam>
-void for_each_pair(ParallelPolicy policy, int begin, int end, Lam lambda, bool half = true)
-{
-    auto actions = [&lambda](const auto pair)
-    {
-        auto [i, j] = pair;
-        lambda(i, j);
-    };
-    if (half)
-    {
-        PairIterator start(end), stop(end, end * (end + 1) / 2);
-        for_each(policy, start, stop, actions);
-    }
-    else
-    {
-        FullPairIterator start(end - begin), stop(end - begin, (end - begin) * (end - begin));
-        for_each(policy, start, stop, actions);
-    }
-}
 } // namespace dissolve
 
 // Join elements into a delimited string
diff --git a/tests/algorithms/pairIterator.cpp b/tests/algorithms/pairIterator.cpp
index 90533f19a0..0888395de4 100644
--- a/tests/algorithms/pairIterator.cpp
+++ b/tests/algorithms/pairIterator.cpp
@@ -32,7 +32,7 @@ void for_each_test(bool unordered)
 
     int sum = 0;
     dissolve::for_each_pair(
-        ParallelPolicies::seq, 0, size,
+        ParallelPolicies::seq, size,
         [&sum, &store](const auto i, const auto j) {
             sum += store[{i, j}];
         },
diff --git a/tests/classes/cells2.cpp b/tests/classes/cells2.cpp
index bea413e9bf..40f5e8aca8 100644
--- a/tests/classes/cells2.cpp
+++ b/tests/classes/cells2.cpp
@@ -94,7 +94,7 @@ class CellsPBCTest : public ::testing::Test
         {
             const auto &nbrs = cellArray.neighbours(*cellArray.cell(n));
 
-            dissolve::for_each_pair(ParallelPolicies::seq, nbrs.begin(), nbrs.end(),
+            dissolve::for_each_pair(ParallelPolicies::seq, nbrs,
                                     [&](auto i, auto &nbri, auto j, auto &nbrj)
                                     {
                                         if (i != j)
diff --git a/tests/classes/cells3.cpp b/tests/classes/cells3.cpp
index 3e6dd50602..2a205e71b6 100644
--- a/tests/classes/cells3.cpp
+++ b/tests/classes/cells3.cpp
@@ -59,7 +59,7 @@ class CellsEnergyTest : public ::testing::Test
         auto *box = cfg->box();
         auto *pp = dissolve_.pairPotential("Ar", "Ar");
         auto energy = 0.0;
-        dissolve::for_each_pair(ParallelPolicies::seq, cfg->molecules().begin(), cfg->molecules().end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, cfg->molecules(),
                                 [&](int i, const auto &molI, int j, const auto &molJ)
                                 {
                                     if (i == j)
@@ -79,7 +79,7 @@ class CellsEnergyTest : public ::testing::Test
         auto *box = cfg->box();
         auto *pp = dissolve_.pairPotential("Ar", "Ar");
         auto energy = 0.0;
-        dissolve::for_each_pair(ParallelPolicies::seq, cfg->molecules().begin(), cfg->molecules().end(),
+        dissolve::for_each_pair(ParallelPolicies::seq, cfg->molecules(),
                                 [&](int i, const auto &molI, int j, const auto &molJ)
                                 {
                                     if (i == j)
diff --git a/tests/data/Averaging.restart.txt b/tests/data/Averaging.restart.txt
new file mode 100644
index 0000000000..253acabf8e
--- /dev/null
+++ b/tests/data/Averaging.restart.txt
@@ -0,0 +1,378 @@
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
diff --git a/tests/data/Averaging2.restart.txt b/tests/data/Averaging2.restart.txt
new file mode 100644
index 0000000000..253acabf8e
--- /dev/null
+++ b/tests/data/Averaging2.restart.txt
@@ -0,0 +1,378 @@
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4
+NO_FINGERPRINT
+4
+A-A 0 A A
+
+1 False
+1  2
+A-B 0 A B
+
+1 False
+1  2
+A-C 0 A C
+
+1 False
+1  4
+A-D 0 A D
+
+1 False
+1  4