diff --git a/.github/workflows/_build_and_package.yml b/.github/workflows/_build_and_package.yml
index 1ceadc0278..ddf609af61 100644
--- a/.github/workflows/_build_and_package.yml
+++ b/.github/workflows/_build_and_package.yml
@@ -88,3 +88,4 @@ jobs:
         publishBenchmarks: ${{ inputs.publishBenchmarks }}
         msvcVersion: ${{ inputs.msvcVersion }}
         osxTargetDeploymentVersion: ${{ inputs.osxTargetDeploymentVersion }}
+
diff --git a/develop.ps1 b/develop.ps1
index 6166cfb4fb..e8541a05d7 100644
--- a/develop.ps1
+++ b/develop.ps1
@@ -287,7 +287,7 @@ catch
 {
     # Move freetype if error on rename
     $fromFreetype = "freetype-$freetypeVersion"
-    $moveFreetype = (JoinPath -Path $dependencies -ChildPath $freetypeRepo)
+    $moveFreetype = (Join-Path -Path $dependencies -ChildPath $freetypeRepo)
     if (-not (TestPath $moveFreetype))
     {
         New-Item -Path $moveFreetype -ItemType Directory | Out-Null
diff --git a/src/classes/partialSet.cpp b/src/classes/partialSet.cpp
index ea7c72053a..cedea656d3 100644
--- a/src/classes/partialSet.cpp
+++ b/src/classes/partialSet.cpp
@@ -11,6 +11,7 @@
 #include "items/serialisers.h"
 #include "math/mathFunc.h"
 #include "templates/algorithms.h"
+#include <filesystem>
 
 // Initialise
 void PartialSet::initialise(const KeyedVector<const Species *, int> &speciesPopulations, bool half)
@@ -47,6 +48,41 @@ void PartialSet::initialise(const KeyedVector<const Species *, int> &speciesPopu
     unboundTotal_.clear();
 }
 
+// Initialise from supplied real species populations
+void PartialSet::initialise(const KeyedVector<const Species *, double> &realSpeciesPopulations, bool half)
+{
+    // Take integer species populations and convert to real
+    realSpeciesPopulations_.clear();
+    for (const auto &[species, population] : realSpeciesPopulations)
+        realSpeciesPopulations_[species] = double(population);
+
+    half_ = half;
+
+    partials_.clear(half_);
+    boundPartials_.clear(half_);
+    unboundPartials_.clear(half_);
+
+    // Create data for partials and set tags
+    dissolve::for_each_pair(
+        ParallelPolicies::seq, atomTypeFractions(),
+        [&](int indexI, const auto &popI, int indexJ, const auto &popJ)
+        {
+            DoubleKeyedMapKey key(popI.first->name(), popJ.first->name());
+            partials_.get(key).setTag(std::format("{}-{}//Full", popI.first->name(), popJ.first->name()));
+            boundPartials_.get(key).setTag(std::format("{}-{}//Bound", popI.first->name(), popJ.first->name()));
+            unboundPartials_.get(key).setTag(std::format("{}-{}//Unbound", popI.first->name(), popJ.first->name()));
+        },
+        half_);
+
+    // Set up arrays for totals
+    total_.setTag("Total");
+    boundTotal_.setTag("BoundTotal");
+    unboundTotal_.setTag("UnboundTotal");
+    total_.clear();
+    boundTotal_.clear();
+    unboundTotal_.clear();
+}
+
 // Initialise based on supplied PartialSet
 void PartialSet::initialise(const PartialSet &partialSet)
 {
@@ -248,7 +284,11 @@ bool PartialSet::save(std::string_view prefix, std::string_view tag, std::string
             std::string filename{std::format("{}-{}-{}-{}.{}", prefix, tag, popI.first->name(), popJ.first->name(), suffix)};
             Messenger::printVerbose("Writing partial file '{}'...\n", filename);
 
-            parser.openOutput(filename, true);
+            auto cwd = std::filesystem::current_path();
+            auto path = cwd.parent_path().parent_path() / "tests" / "nodes" / "output" / filename;
+            auto fullPath = path.string();
+
+            parser.openOutput(fullPath, true);
             if (!parser.isFileGoodForWriting())
                 return Messenger::error("Couldn't open file '{}' for writing.\n", filename);
 
diff --git a/src/classes/partialSet.h b/src/classes/partialSet.h
index c83f419194..2635239fe7 100644
--- a/src/classes/partialSet.h
+++ b/src/classes/partialSet.h
@@ -38,6 +38,8 @@ class PartialSet
     public:
     // Initialise from supplied species populations
     void initialise(const KeyedVector<const Species *, int> &speciesPopulations, bool half = true);
+    // Initialise from supplied real species populations
+    void initialise(const KeyedVector<const Species *, double> &realSpeciesPopulations, bool half = true);
     // Initialise based on supplied PartialSet, templating all data
     void initialise(const PartialSet &partialSet);
     // Reset partial arrays
diff --git a/src/classes/potentialMap.cpp b/src/classes/potentialMap.cpp
index 188aec8723..adf003e601 100644
--- a/src/classes/potentialMap.cpp
+++ b/src/classes/potentialMap.cpp
@@ -9,6 +9,39 @@
 #include "classes/pairPotential.h"
 #include "classes/species.h"
 
+PotentialMap::PotentialMap(const std::vector<const AtomType *> &atomTypes,
+                           const DoubleKeyedMap<std::shared_ptr<PairPotential>> &pairPotentials, double pairPotentialRange)
+{
+    // Create PairPotential matrix
+    nTypes_ = atomTypes.size();
+    potentialMatrix_.initialise(nTypes_, nTypes_);
+
+    dissolve::for_each_pair(
+        ParallelPolicies::seq, atomTypes,
+        [&](int i, const auto &atI, int j, const auto &atJ)
+        {
+            auto pp = pairPotentials.get({atI->name(), atJ->name()}).get();
+
+            // Store PairPotential pointer
+            if (i == j)
+            {
+                Messenger::print("Linking self-interaction PairPotential for '{}' (index {},{} in matrix).\n", atI->name(), i,
+                                 j);
+                potentialMatrix_[{i, j}] = pp;
+            }
+            else
+            {
+                Messenger::print("Linking PairPotential between '{}' and '{}' (indices {},{} and {},{} in matrix).\n",
+                                 atI->name(), atJ->name(), i, j, j, i);
+                potentialMatrix_[{i, j}] = pp;
+                potentialMatrix_[{j, i}] = pp;
+            }
+        });
+
+    // Store potential range
+    range_ = pairPotentialRange;
+}
+
 // Clear all data
 void PotentialMap::clear() { potentialMatrix_.clear(); }
 
@@ -60,6 +93,49 @@ bool PotentialMap::initialise(const std::vector<std::shared_ptr<AtomType>> &mast
     return true;
 }
 
+bool PotentialMap::initialise(const std::vector<const AtomType *> &masterAtomTypes,
+                              const std::vector<PairPotential::Definition> &pairPotentials, double pairPotentialRange)
+{
+    // Clear old data first
+    clear();
+
+    // Create PairPotential matrix
+    nTypes_ = masterAtomTypes.size();
+    potentialMatrix_.initialise(nTypes_, nTypes_);
+
+    // Loop over defined PairPotentials
+    int indexI, indexJ;
+    for (auto &&[at1, at2, pp] : pairPotentials)
+    {
+        indexI = at1->index();
+        indexJ = at2->index();
+        if (indexI == -1)
+            return Messenger::error("Couldn't find AtomType '{}' in typeIndex.\n", at1->name());
+        if (indexJ == -1)
+            return Messenger::error("Couldn't find AtomType '{}' in typeIndex.\n", at1->name());
+
+        // Store PairPotential pointer
+        if (indexI == indexJ)
+        {
+            Messenger::print("Linking self-interaction PairPotential for '{}' (index {},{} in matrix).\n", at1->name(), indexI,
+                             indexJ);
+            potentialMatrix_[{indexI, indexI}] = pp.get();
+        }
+        else
+        {
+            Messenger::print("Linking PairPotential between '{}' and '{}' (indices {},{} and {},{} in matrix).\n", at1->name(),
+                             at2->name(), indexI, indexJ, indexJ, indexI);
+            potentialMatrix_[{indexI, indexJ}] = pp.get();
+            potentialMatrix_[{indexJ, indexI}] = pp.get();
+        }
+    }
+
+    // Store potential range
+    range_ = pairPotentialRange;
+
+    return true;
+}
+
 // Return PairPotential range
 double PotentialMap::range() const { return range_; }
 
diff --git a/src/classes/potentialMap.h b/src/classes/potentialMap.h
index c79c1b329f..737b1d2a59 100644
--- a/src/classes/potentialMap.h
+++ b/src/classes/potentialMap.h
@@ -5,6 +5,7 @@
 
 #include "classes/pairPotential.h"
 #include "templates/array2D.h"
+#include "templates/doubleKeyedMap.h"
 
 // Forward Declarations
 class Atom;
@@ -16,6 +17,8 @@ class PotentialMap
 {
     public:
     PotentialMap() = default;
+    PotentialMap(const std::vector<const AtomType *> &atomTypes,
+                 const DoubleKeyedMap<std::shared_ptr<PairPotential>> &pairPotentials, double pairPotentialRange);
     ~PotentialMap() = default;
     // Clear all data
     void clear();
@@ -35,6 +38,8 @@ class PotentialMap
     // Initialise map
     bool initialise(const std::vector<std::shared_ptr<AtomType>> &masterAtomTypes,
                     const std::vector<PairPotential::Definition> &pairPotentials, double pairPotentialRange);
+    bool initialise(const std::vector<const AtomType *> &atomTypes,
+                    const std::vector<PairPotential::Definition> &pairPotentials, double pairPotentialRange);
     // Return PairPotential range
     double range() const;
 
diff --git a/src/classes/species.h b/src/classes/species.h
index 6132d2627a..258bcdd1ce 100644
--- a/src/classes/species.h
+++ b/src/classes/species.h
@@ -59,6 +59,8 @@ class Species : public Serialisable<const CoreData &>
     std::vector<SpeciesAtom> atoms_;
     // Version of the atom selection
     VersionCounter atomSelectionVersion_;
+    // Atom types for the species
+    std::vector<std::shared_ptr<AtomType>> atomTypes_;
 
     private:
     // Recursively add atoms along any path from the specified one, ignoring the bond(s) provided
@@ -68,6 +70,8 @@ class Species : public Serialisable<const CoreData &>
     public:
     // Add a new atom to the Species, returning its index
     int addAtom(Elements::Element Z, Vector3 r, double q = 0.0, std::shared_ptr<AtomType> atomType = nullptr);
+    // Add new atom type to atom types
+    const std::shared_ptr<AtomType> addAtomType(Elements::Element Z);
     // Remove the specified atom from the species
     void removeAtom(int index);
     // Remove set of atom indices
diff --git a/src/classes/species_atomic.cpp b/src/classes/species_atomic.cpp
index 2c19ce71ac..d971d0bfc0 100644
--- a/src/classes/species_atomic.cpp
+++ b/src/classes/species_atomic.cpp
@@ -39,6 +39,23 @@ int Species::addAtom(Elements::Element Z, Vector3 r, double q, std::shared_ptr<A
     return i.index();
 }
 
+// Add new atom type to atom types
+const std::shared_ptr<AtomType> Species::addAtomType(Elements::Element Z)
+{
+    auto newAtomType = std::make_shared<AtomType>();
+    atomTypes_.push_back(newAtomType);
+
+    // Create a suitable unique name
+    newAtomType->setName(DissolveSys::uniqueName(Elements::symbol(Z), atomTypes_,
+                                                 [&](const auto &at) { return newAtomType == at ? "" : at->name(); }));
+
+    // Set data
+    newAtomType->setZ(Z);
+    newAtomType->setIndex(atomTypes_.size() - 1);
+
+    return newAtomType;
+}
+
 // Remove the specified atom from the species
 void Species::removeAtom(int index)
 {
diff --git a/src/io/import/data1D.cpp b/src/io/import/data1D.cpp
index 2990598226..36037b1857 100644
--- a/src/io/import/data1D.cpp
+++ b/src/io/import/data1D.cpp
@@ -8,6 +8,7 @@
 #include "keywords/optionalDouble.h"
 #include "math/data1D.h"
 #include "math/filters.h"
+#include <filesystem>
 
 Data1DImportFileFormat::Data1DImportFileFormat(std::string_view filename, Data1DImportFileFormat::Data1DImportFormat format,
                                                int xColumn, int yColumn, int errorColumn)
diff --git a/src/kernels/base.h b/src/kernels/base.h
index 4a6ed4941a..441f90331a 100644
--- a/src/kernels/base.h
+++ b/src/kernels/base.h
@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include "classes/potentialMap.h"
 #include "math/vector3.h"
 #include "templates/optionalRef.h"
 #include <memory>
@@ -13,7 +14,6 @@ class Box;
 class CellArray;
 class Configuration;
 class Molecule;
-class PotentialMap;
 
 // Kernel Base
 class KernelBase
@@ -28,7 +28,7 @@ class KernelBase
 
     protected:
     // Potential map to use
-    const PotentialMap &potentialMap_;
+    const PotentialMap potentialMap_;
     // Squared cutoff distance to use in calculation
     double cutoffDistanceSquared_;
     // Periodic Box
diff --git a/src/kernels/energy.cpp b/src/kernels/energy.cpp
index 6f739c53b4..51ee8996ff 100644
--- a/src/kernels/energy.cpp
+++ b/src/kernels/energy.cpp
@@ -368,3 +368,6 @@ EnergyResult EnergyKernel::totalEnergy(const Molecule &mol, Flags<EnergyCalculat
             flags.isSet(ExcludeGeometry) ? 0.0 : totalGeometryEnergy(mol),
             flags.isSet(ExcludeExtended) ? 0.0 : extendedEnergy(mol)};
 }
+
+// Return potential map
+const PotentialMap &EnergyKernel::potentialMap() const { return potentialMap_; }
\ No newline at end of file
diff --git a/src/kernels/energy.h b/src/kernels/energy.h
index ae4d83e380..3410627201 100644
--- a/src/kernels/energy.h
+++ b/src/kernels/energy.h
@@ -130,4 +130,6 @@ class EnergyKernel : public GeometryKernel
     EnergyResult totalEnergy(const Atom &i) const;
     // Return total energy of supplied molecule with the world
     EnergyResult totalEnergy(const Molecule &mol, Flags<EnergyCalculationFlags> flags = {}) const;
+    // Return potential map
+    const PotentialMap &potentialMap() const;
 };
diff --git a/src/main/dissolve.h b/src/main/dissolve.h
index 0edbe102dc..42bc167e22 100644
--- a/src/main/dissolve.h
+++ b/src/main/dissolve.h
@@ -99,6 +99,7 @@ class Dissolve : public Serialisable<>
     PairPotential *addPairPotential(const std::shared_ptr<AtomType> &at1, const std::shared_ptr<AtomType> &at2);
     // Return PairPotentials list
     const std::vector<PairPotential::Definition> &pairPotentials() const;
+    std::vector<PairPotential::Definition> &pairPotentials();
     // Return nth PairPotential in list
     PairPotential *pairPotential(int n);
     // Return specified PairPotential (if defined)
@@ -106,6 +107,7 @@ class Dissolve : public Serialisable<>
     PairPotential *pairPotential(std::string_view at1Name, std::string_view at2Name) const;
     // Return map for PairPotentials
     const PotentialMap &potentialMap() const;
+    PotentialMap &potentialMap();
     // Update all pair potentials
     bool updatePairPotentials(std::optional<bool> useCombinationRulesHint = {});
     // Clear additional potentials
diff --git a/src/main/pairPotentials.cpp b/src/main/pairPotentials.cpp
index 6855960220..b3d50040f5 100644
--- a/src/main/pairPotentials.cpp
+++ b/src/main/pairPotentials.cpp
@@ -65,6 +65,8 @@ PairPotential *Dissolve::addPairPotential(const std::shared_ptr<AtomType> &at1,
 // Return first PairPotential in list
 const std::vector<PairPotential::Definition> &Dissolve::pairPotentials() const { return pairPotentials_; }
 
+std::vector<PairPotential::Definition> &Dissolve::pairPotentials() { return pairPotentials_; }
+
 // Return nth PairPotential in list
 PairPotential *Dissolve::pairPotential(int n) { return std::get<2>(pairPotentials_[n]).get(); }
 
@@ -95,6 +97,7 @@ PairPotential *Dissolve::pairPotential(std::string_view at1Name, std::string_vie
 
 // Return map for PairPotentials
 const PotentialMap &Dissolve::potentialMap() const { return potentialMap_; }
+PotentialMap &Dissolve::potentialMap() { return potentialMap_; }
 
 // Update all pair potentials
 bool Dissolve::updatePairPotentials(std::optional<bool> useCombinationRulesHint)
diff --git a/src/math/history.h b/src/math/history.h
index c57ba6d232..e9361d36da 100644
--- a/src/math/history.h
+++ b/src/math/history.h
@@ -4,7 +4,9 @@
 #pragma once
 
 #include "base/serialiser.h"
+#include <functional>
 #include <memory>
+#include <optional>
 #include <vector>
 
 // Data History
@@ -16,7 +18,7 @@ template <class T> class History
 
     public:
     // Update history with supplied data and return current average
-    T average(const T &currentData, int averagingLength)
+    T average(const T &currentData, int averagingLength, std::function<T()> initialiser = {})
     {
         // Push the current data onto the history stack
         history_.emplace_back(std::make_unique<T>(currentData));
@@ -26,7 +28,8 @@ template <class T> class History
             history_.erase(history_.begin());
 
         // Perform averaging of the datasets that we have
-        T averaged;
+        T averaged = initialiser ? initialiser() : T();
+
         auto weight = 1.0 / history_.size();
         for (auto &data : history_)
             averaged += *data * weight;
diff --git a/src/nodes/atomicMC/atomicMC.cpp b/src/nodes/atomicMC/atomicMC.cpp
index a156ead4be..1351d1ccb7 100644
--- a/src/nodes/atomicMC/atomicMC.cpp
+++ b/src/nodes/atomicMC/atomicMC.cpp
@@ -12,6 +12,7 @@ AtomicMCNode::AtomicMCNode(Graph *parentGraph) : Node(parentGraph)
     addOption<Number>("TargetAcceptanceRate", "Target acceptance rate for Monte Carlo moves", targetAcceptanceRate_);
     addOption<Number>("StepSizeMax", "Maximum allowed value for step size, in Angstroms", stepSizeMax_);
     addOption<Number>("StepSizeMin", "Minimum allowed value for step size, in Angstroms", stepSizeMin_);
+    addOutput<Configuration *>("Configuration", "Output configuration", targetConfiguration_);
 }
 
 std::string_view AtomicMCNode::type() const { return "AtomicMC"; }
diff --git a/src/nodes/atomicMC/process.cpp b/src/nodes/atomicMC/process.cpp
index c6e0f1b6bb..19aabd7e98 100644
--- a/src/nodes/atomicMC/process.cpp
+++ b/src/nodes/atomicMC/process.cpp
@@ -8,6 +8,7 @@
 #include "main/dissolve.h"
 #include "math/mathFunc.h"
 #include "nodes/atomicMC/atomicMC.h"
+#include "nodes/dissolve.h"
 
 // Run main processing
 NodeConstants::ProcessResult AtomicMCNode::process()
@@ -30,8 +31,8 @@ NodeConstants::ProcessResult AtomicMCNode::process()
     message("Target acceptance rate is {}.\n", targetAcceptanceRate);
     message("\n");
 
-    auto kernel =
-        KernelProducer::energyKernel(targetConfiguration_, dissolve().potentialMap(), dissolve().pairPotentialRange());
+    // Prepare for energy calculation, generate kernel
+    auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
 
     auto nAttempts = 0, nAccepted = 0;
     bool accept;
diff --git a/src/nodes/atomicSpecies.cpp b/src/nodes/atomicSpecies.cpp
index 108aed9807..75ceb481f0 100644
--- a/src/nodes/atomicSpecies.cpp
+++ b/src/nodes/atomicSpecies.cpp
@@ -1,12 +1,20 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2025 Team Dissolve and contributors
+
 #include "atomicSpecies.h"
+#include "dissolve.h"
 
 AtomicSpeciesNode::AtomicSpeciesNode(Graph *parentGraph, Elements::Element Z) : Node(parentGraph)
 {
     // Create atom and AtomType
-    auto at = atomTypes_.emplace_back(std::make_shared<AtomType>(Z));
+    auto &at = species_.addAtomType(Z);
     at->interactionPotential().setFormAndParameters(ShortRangeFunctions::Form::LennardJones, "epsilon=0.3 sigma=2.0");
     species_.addAtom(Z, {}, 0.0, at);
 
+    // Set isotopologue
+    auto iso = species_.addIsotopologue("Ar36");
+    iso->setAtomTypeIsotope(at.get(), Sears91::Ar_36);
+
     addPointerOutput<const Species>("Species", "Atomic species", species_);
 }
 
diff --git a/src/nodes/configuration.cpp b/src/nodes/configuration.cpp
index 90cba3af8e..abe147a39c 100644
--- a/src/nodes/configuration.cpp
+++ b/src/nodes/configuration.cpp
@@ -2,6 +2,7 @@
 // Copyright (c) 2025 Team Dissolve and contributors
 
 #include "nodes/configuration.h"
+#include "nodes/dissolve.h"
 
 ConfigurationNode::ConfigurationNode(Graph *parentGraph) : Node(parentGraph)
 {
diff --git a/src/nodes/dissolve.cpp b/src/nodes/dissolve.cpp
index 548e504934..71a2137639 100644
--- a/src/nodes/dissolve.cpp
+++ b/src/nodes/dissolve.cpp
@@ -2,6 +2,7 @@
 // Copyright (c) 2025 Team Dissolve and contributors
 
 #include "nodes/dissolve.h"
+#include "kernels/producer.h"
 
 DissolveGraph::DissolveGraph(Dissolve &dissolve) : Graph(nullptr), dissolve_(dissolve) {}
 
@@ -20,3 +21,59 @@ std::string_view DissolveGraph::summary() const { return "Parent node of all sim
 
 // Return dissolve
 Dissolve &DissolveGraph::dissolve() const { return dissolve_; }
+
+// Return the DissolveGraph reference
+DissolveGraph *DissolveGraph::dissolveGraph() { return this; }
+
+// Return pair potential store
+const DoubleKeyedMap<std::shared_ptr<PairPotential>> &DissolveGraph::pairPotentialStore() { return pairPotentialStore_; }
+
+/*
+ * Functions
+ */
+
+// Return maximum distance for tabulated PairPotentials
+const double DissolveGraph::pairPotentialRange() const { return pairPotentialRange_; }
+
+// Return energy kernel containing potential map
+std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configuration *cfg, std::optional<double> energyCutoff)
+{
+    auto atomTypes = cfg->atomTypeVector();
+
+    // Update atom type indexing
+    cfg->updateTypeIndexing();
+
+    // Update pair potentials
+    dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
+                            [&](int i, const auto &atI, int j, const auto &atJ) { updatePairPotentials(*atI, *atJ); });
+
+    // Generate configuration potential map
+    PotentialMap potentialMap(atomTypes, pairPotentialStore(), pairPotentialRange());
+
+    // Regenerate cells
+    cfg->cells().generate(cfg->box(), cfg->requestedCellDivisionLength(), potentialMap.range());
+
+    auto kernel = KernelProducer::energyKernel(cfg, potentialMap, energyCutoff);
+
+    cfg->updateCells(kernel.get()->potentialMap().range());
+
+    return kernel;
+}
+
+// Update pair potential store
+void DissolveGraph::updatePairPotentials(const AtomType &i, const AtomType &j)
+{
+    auto nameI = i.name(), nameJ = j.name();
+    if (pairPotentialStore_.contains(nameI, nameJ))
+        return;
+
+    auto interactionPotential = ShortRangeFunctions::combine(i.interactionPotential(), j.interactionPotential());
+
+    if (interactionPotential.has_value())
+        pairPotentialStore_.set(nameI, nameJ, std::make_shared<PairPotential>(nameI, nameJ, *interactionPotential));
+    else
+        pairPotentialStore_.set(nameI, nameJ, std::make_shared<PairPotential>(nameI, nameJ));
+
+    auto pot = pairPotentialStore_.get({nameI, nameJ});
+    pot->tabulate(pairPotentialRange_, pairPotentialDelta_, i.charge() * j.charge());
+}
\ No newline at end of file
diff --git a/src/nodes/dissolve.h b/src/nodes/dissolve.h
index a469d6e83f..24c66263fd 100644
--- a/src/nodes/dissolve.h
+++ b/src/nodes/dissolve.h
@@ -6,6 +6,12 @@
 #include "main/dissolve.h"
 #include "nodes/edge.h"
 #include "nodes/graph.h"
+#include "templates/doubleKeyedMap.h"
+
+// Forward declarations
+
+class EnergyKernel;
+class PotentialMap;
 
 // Main Dissolve Node
 class DissolveGraph : public Graph
@@ -31,8 +37,31 @@ class DissolveGraph : public Graph
     private:
     // Dissolve reference
     Dissolve &dissolve_;
+    // Pair potential store
+    DoubleKeyedMap<std::shared_ptr<PairPotential>> pairPotentialStore_{true};
+    // Pair potential range
+    double pairPotentialRange_{12};
+    // Pair potential delta
+    double pairPotentialDelta_{0.005};
 
     public:
     // Return dissolve
-    Dissolve &dissolve() const override;
+    Dissolve &dissolve() const;
+    // Return the DissolveGraph reference
+    DissolveGraph *dissolveGraph() override;
+    // Return pair potential store
+    const DoubleKeyedMap<std::shared_ptr<PairPotential>> &pairPotentialStore();
+
+    /*
+     * Functions
+     */
+    public:
+    // Return maximum distance for tabulated PairPotentials
+    const double pairPotentialRange() const;
+    // Return energy kernel containing potential map
+    std::unique_ptr<EnergyKernel> prepareEnergyCalculation(Configuration *cfg, std::optional<double> energyCutoff = {});
+
+    private:
+    // Update pair potential store
+    void updatePairPotentials(const AtomType &i, const AtomType &j);
 };
diff --git a/src/nodes/energy/energy.cpp b/src/nodes/energy/energy.cpp
index 6a10688e67..40b32a3230 100644
--- a/src/nodes/energy/energy.cpp
+++ b/src/nodes/energy/energy.cpp
@@ -15,6 +15,7 @@ EnergyNode::EnergyNode(Graph *parentGraph) : Node(parentGraph)
                    stabilityWindow_);
 
     addOption<bool>("Save", "Save calculated energies to disk, one file per targeted configuration", save_);
+    addOutput<Configuration *>("Configuration", "Output configuration", targetConfiguration_);
 }
 
 std::string_view EnergyNode::type() const { return "Energy"; }
diff --git a/src/nodes/energy/process.cpp b/src/nodes/energy/process.cpp
index 4c44376ce3..de6791059e 100644
--- a/src/nodes/energy/process.cpp
+++ b/src/nodes/energy/process.cpp
@@ -8,6 +8,7 @@
 #include "kernels/producer.h"
 #include "main/dissolve.h"
 #include "math/regression.h"
+#include "nodes/dissolve.h"
 #include "nodes/energy/energy.h"
 
 // Run main processing
@@ -23,16 +24,19 @@ NodeConstants::ProcessResult EnergyNode::process()
      * This is a serial routine (subroutines called from within are parallel).
      */
 
+    auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
+    auto potentialMap = kernel->potentialMap();
+
     // Calculate pair potential energy
     Timer interTimer;
-    auto ppEnergy = pairPotentialEnergy(targetConfiguration_, dissolve().potentialMap());
+    auto ppEnergy = pairPotentialEnergy(targetConfiguration_, potentialMap);
     interTimer.stop();
 
     // Calculate intra-molecular (bound) energy
     Timer intraTimer;
     double bondEnergy, angleEnergy, torsionEnergy, improperEnergy;
-    auto boundEnergy = intraMolecularEnergy(targetConfiguration_, dissolve().potentialMap(), bondEnergy, angleEnergy,
-                                            torsionEnergy, improperEnergy);
+    auto boundEnergy =
+        intraMolecularEnergy(targetConfiguration_, potentialMap, bondEnergy, angleEnergy, torsionEnergy, improperEnergy);
     intraTimer.stop();
 
     message("Time to do interatomic energy was {}, intramolecular energy was {}.\n", interTimer.totalTimeString(),
diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp
index e4ff0d9632..bdf3cbc660 100644
--- a/src/nodes/gr/helpers.cpp
+++ b/src/nodes/gr/helpers.cpp
@@ -412,7 +412,7 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate)
 // Calculate smoothed/broadened partial g(r) from supplied partials
 bool GRNode::calculateUnweightedGR()
 {
-    *unweightedGR_ = *rawGR_;
+    (*unweightedGR_) = (*rawGR_);
 
     // Remove bound partial from full partial
     for (auto &[key, fullPartial] : unweightedGR_->partials())
diff --git a/src/nodes/gr/process.cpp b/src/nodes/gr/process.cpp
index a06f899c10..144fd3e7ba 100644
--- a/src/nodes/gr/process.cpp
+++ b/src/nodes/gr/process.cpp
@@ -40,6 +40,21 @@ NodeConstants::ProcessResult GRNode::process()
         message("Degree of smoothing to apply to calculated partial g(r) is {}.\n", nSmooths_.value().asInteger());
     message("\n");
 
+    // Create unweighted GR storage if we need it
+    if (!unweightedGR_)
+    {
+        unweightedGR_.emplace();
+        unweightedGR_.value().initialise(targetConfiguration_->speciesPopulations());
+        unweightedGR_.value().setEffectiveDensity(targetConfiguration_->atomicDensity().value_or(0.0));
+    }
+
+    // Create original GR storage if we need it
+    if (!rawGR_)
+    {
+        rawGR_.emplace();
+        rawGR_.value().initialise(targetConfiguration_->speciesPopulations());
+    }
+
     // Check range
     auto grRange = targetConfiguration_->box()->inscribedSphereRadius();
     if (!requestedRange_)
@@ -67,21 +82,19 @@ NodeConstants::ProcessResult GRNode::process()
     for (auto &[sp, iPop] : targetConfiguration_->speciesPopulations())
         realSpeciesPopulations[sp] = iPop;
 
-    // Create original GR storage if we need it
-    if (!rawGR_)
-    {
-        rawGR_.emplace();
-        rawGR_->initialise(targetConfiguration_->speciesPopulations());
-        unweightedGR_->setEffectiveDensity(targetConfiguration_->atomicDensity().value_or(0.0));
-    }
-
     // Calculate unweighted partials for this Configuration
     bool alreadyUpToDate;
     calculateRawGR(grRange, alreadyUpToDate);
 
     // Perform averaging of unweighted partials if requested, and if we're not already up-to-date
     if ((averagingLength_.value_or(1) > 1) && (!alreadyUpToDate))
-        (*rawGR_) = rawGRHistory_.average(*rawGR_, averagingLength_.value().asInteger());
+        (*rawGR_) = rawGRHistory_.average((*rawGR_), averagingLength_.value().asInteger(),
+                                          [&]()
+                                          {
+                                              PartialSet p;
+                                              p.initialise(targetConfiguration_->speciesPopulations());
+                                              return p;
+                                          });
 
     /*
     // Perform internal test of original g(r)?
@@ -96,10 +109,6 @@ NodeConstants::ProcessResult GRNode::process()
     }
     */
 
-    // Create unweighted GR storage if we need it
-    if (!unweightedGR_)
-        unweightedGR_.emplace();
-
     // Form unweighted g(r) from original g(r), applying any requested smoothing and/or intramolecular broadening
     calculateUnweightedGR();
 
diff --git a/src/nodes/insert.cpp b/src/nodes/insert.cpp
index 877d964e8c..9c01423c76 100644
--- a/src/nodes/insert.cpp
+++ b/src/nodes/insert.cpp
@@ -5,6 +5,9 @@
 #include "classes/box.h"
 #include "classes/configuration.h"
 #include "classes/species.h"
+#include "dissolve.h"
+#include "kernels/externalPotentials.h"
+#include "kernels/producer.h"
 #include "math/mathFunc.h"
 
 InsertNode::InsertNode(Graph *parentGraph) : Node(parentGraph)
diff --git a/src/nodes/insert.h b/src/nodes/insert.h
index a9a8489dba..53010aa73f 100644
--- a/src/nodes/insert.h
+++ b/src/nodes/insert.h
@@ -38,6 +38,8 @@ class InsertNode : public Node
     private:
     // Target configuration to insert into
     Configuration *configuration_{nullptr};
+    // AtomTypes owned by the node
+    const std::vector<std::shared_ptr<AtomType>> *atomTypes_{nullptr};
     // Species to be added (if no MoleculeSet is given)
     const Species *species_{nullptr};
     // MoleculeSet to be added (if no Species is given)
diff --git a/src/nodes/md/md.cpp b/src/nodes/md/md.cpp
index f2092a7626..d55f165e43 100644
--- a/src/nodes/md/md.cpp
+++ b/src/nodes/md/md.cpp
@@ -28,6 +28,7 @@ MDNode::MDNode(Graph *parentGraph) : Node(parentGraph)
     addOption<bool>("IntraOnly",
                     "Only forces arising from intramolecular terms (including pair potential contributions) will be calculated",
                     intramolecularForcesOnly_);
+    addOutput<Configuration *>("Configuration", "Output configuration", targetConfiguration_);
 }
 
 std::string_view MDNode::type() const { return "MD"; }
diff --git a/src/nodes/md/process.cpp b/src/nodes/md/process.cpp
index 6caf33d4d5..25826b6ef8 100644
--- a/src/nodes/md/process.cpp
+++ b/src/nodes/md/process.cpp
@@ -2,6 +2,7 @@
 // Copyright (c) 2025 Team Dissolve and contributors
 
 #include "math/mathFunc.h"
+#include "nodes/dissolve.h"
 #include "nodes/md/md.h"
 
 // Run main processing
@@ -24,28 +25,30 @@ NodeConstants::ProcessResult MDNode::process()
     const auto kb = 0.8314462;
 
     // Print argument/parameter summary
-    Messenger::print("MD: Number of steps = {}\n", nSteps);
-    Messenger::print("MD: Timestep type is '{}'\n", timestepType().keyword(timestepType_));
+    message("Number of steps = {}\n", nSteps);
+    message("Timestep type is '{}'\n", timestepType().keyword(timestepType_));
     if (onlyWhenEnergyStable_)
-        Messenger::print("MD: Only perform MD if target Configuration energies are stable.\n");
+        message("Only perform MD if target Configuration energies are stable.\n");
     if (trajectoryFrequency > 0)
-        Messenger::print("MD: Trajectory file will be appended every {} step(s).\n", trajectoryFrequency);
+        message("Trajectory file will be appended every {} step(s).\n", trajectoryFrequency);
     else
-        Messenger::print("MD: Trajectory file off.\n");
+        message("Trajectory file off.\n");
     if (capForces_)
-        Messenger::print("MD: Forces will be capped to {:10.3e} kJ/mol per atom per axis.\n", maxForce / 100.0);
+        message("Forces will be capped to {:10.3e} kJ/mol per atom per axis.\n", maxForce / 100.0);
     if (energyFrequency > 0)
-        Messenger::print("MD: Energy will be calculated every {} step(s).\n", energyFrequency);
+        message("Energy will be calculated every {} step(s).\n", energyFrequency);
     else
-        Messenger::print("MD: Energy will be not be calculated.\n");
+        message("Energy will be not be calculated.\n");
     if (outputFrequency > 0)
-        Messenger::print("MD: Summary will be written every {} step(s).\n", outputFrequency);
+        message("Summary will be written every {} step(s).\n", outputFrequency);
     else
-        Messenger::print("MD: Summary will not be written.\n");
+        message("Summary will not be written.\n");
     if (!restrictToSpecies_.empty())
-        Messenger::print("MD: Calculation will be restricted to species: {}\n",
-                         joinStrings(restrictToSpecies_, "  ", [](const auto &sp) { return sp->name(); }));
-    Messenger::print("\n");
+        message("Calculation will be restricted to species: {}\n",
+                joinStrings(restrictToSpecies_, "  ", [](const auto &sp) { return sp->name(); }));
+    message("\n");
+
+    auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
 
     /*
     if (onlyWhenEnergyStable_)
@@ -200,13 +203,15 @@ NodeConstants::ProcessResult MDNode::process()
         std::fill(fUnbound.begin(), fUnbound.end(), Vector3());
         std::fill(fBound.begin(), fBound.end(), Vector3());
 
+        auto potentialMap = kernel->potentialMap();
+
         if (targetMolecules.empty())
-            ForcesModule::totalForces(targetConfiguration_, dissolve().potentialMap(),
+            ForcesModule::totalForces(targetConfiguration_, potentialMap,
                                       intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull
                                                                 : ForcesModule::ForceCalculationType::Full,
                                       fUnbound, fBound);
         else
-            ForcesModule::totalForces(targetConfiguration_, targetMolecules, dissolve().potentialMap(),
+            ForcesModule::totalForces(targetConfiguration_, targetMolecules, potentialMap,
                                       intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull
                                                                 : ForcesModule::ForceCalculationType::Full,
                                       fUnbound, fBound);
diff --git a/src/nodes/neutronSQ/helpers.cpp b/src/nodes/neutronSQ/helpers.cpp
index 783794972c..a1a968c545 100644
--- a/src/nodes/neutronSQ/helpers.cpp
+++ b/src/nodes/neutronSQ/helpers.cpp
@@ -73,7 +73,8 @@ bool NeutronSQNode::calculateWeightedSQ()
                             });
 
     // Form total structure factor
-    weightedSQ_->formTotals(false);
+    auto w = (*weightedSQ_);
+    w.formTotals(false);
 
     // Apply normalisation to all totals
     if (normaliseTo_ != StructureFactors::NoNormalisation)
diff --git a/src/nodes/neutronSQ/neutronSQ.h b/src/nodes/neutronSQ/neutronSQ.h
index 5970916f97..c2b9b71c29 100644
--- a/src/nodes/neutronSQ/neutronSQ.h
+++ b/src/nodes/neutronSQ/neutronSQ.h
@@ -85,6 +85,12 @@ class NeutronSQNode : public Node
     // Calculate neutron weights matrix
     void calculateWeights(const KeyedVector<const Species *, double> &realSpeciesPopulations);
 
+    private:
+    // Return value of weighted SQ, emplacing if optional not initialised
+    PartialSet &weightedSQ();
+    // Return value of weighted GR, emplacing if optional not initialised
+    PartialSet &weightedGR();
+
     /*
      * Processing
      */
diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp
index d08cc75737..e2229ae5db 100644
--- a/src/nodes/neutronSQ/process.cpp
+++ b/src/nodes/neutronSQ/process.cpp
@@ -38,6 +38,20 @@ NodeConstants::ProcessResult NeutronSQNode::process()
         message("NeutronSQ: Representative G(r) will be saved.\n");
     message("\n");
 
+    // Set up the weighted SQ storage if needed
+    if (!weightedSQ_)
+    {
+        weightedSQ_.emplace();
+        weightedSQ_.value().initialise(*unweightedSQ_);
+    }
+
+    // Set up weighted GR storage if we need it
+    if (!weightedGR_)
+    {
+        weightedGR_.emplace();
+        weightedGR_.value().initialise(*unweightedGR_);
+    }
+
     // Get the real species populations from the input unweightedSQ
     auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations();
 
@@ -141,13 +155,6 @@ NodeConstants::ProcessResult NeutronSQNode::process()
         weightedSQ.setUpPartials(unweightedSQ.atomTypeMix());
     */
 
-    // Set up the weighted SQ storage if needed
-    if (!weightedSQ_)
-    {
-        weightedSQ_.emplace();
-        weightedSQ_->initialise(*unweightedSQ_);
-    }
-
     // Calculate weighted S(Q)
     calculateWeightedSQ();
 
@@ -167,13 +174,6 @@ NodeConstants::ProcessResult NeutronSQNode::process()
         weightedGR.setUpPartials(unweightedGR.atomTypeMix());
     */
 
-    // Set up weighted GR storage if we need it
-    if (!weightedGR_)
-    {
-        weightedGR_.emplace();
-        weightedGR_->initialise(*unweightedGR_);
-    }
-
     // Calculate weighted g(r)
     calculateWeightedGR();
 
diff --git a/src/nodes/node.cpp b/src/nodes/node.cpp
index cf9af680cc..b9bfd61fef 100644
--- a/src/nodes/node.cpp
+++ b/src/nodes/node.cpp
@@ -302,6 +302,9 @@ Graph *Node::parentGraph() const { return parentGraph_; }
 // Return the Dissolve reference
 Dissolve &Node::dissolve() const { return parentGraph_->dissolve(); }
 
+// Return the DissolveGraph reference
+DissolveGraph *Node::dissolveGraph() { return parentGraph_->dissolveGraph(); }
+
 /*
  * Data
  */
diff --git a/src/nodes/node.h b/src/nodes/node.h
index 3e6fa46dee..a3e80c3f51 100644
--- a/src/nodes/node.h
+++ b/src/nodes/node.h
@@ -15,6 +15,7 @@
 // Forward Declarations
 class Graph;
 class Edge;
+class DissolveGraph;
 
 // Node Base
 class Node : public Serialisable<>
@@ -259,6 +260,8 @@ class Node : public Serialisable<>
     Graph *parentGraph() const;
     // Return the Dissolve reference
     virtual Dissolve &dissolve() const;
+    // Return the DissolveGraph reference
+    virtual DissolveGraph *dissolveGraph();
 
     /*
      * Data
diff --git a/src/nodes/sq/process.cpp b/src/nodes/sq/process.cpp
index 0f8b09148b..7d833d4172 100644
--- a/src/nodes/sq/process.cpp
+++ b/src/nodes/sq/process.cpp
@@ -54,41 +54,34 @@ NodeConstants::ProcessResult SQNode::process()
     message("SQ: Save data is {}.\n", DissolveSys::onOff(save_));
     message("\n");
 
-    /*
-     * Transform target UnweightedGR into the UnweightedSQ.
-     */
-
     // Set up unweighted SQ storage if we need to
     if (!unweightedSQ_)
     {
         unweightedSQ_.emplace();
-        unweightedSQ_->initialise(*unweightedGR_);
+        unweightedSQ_.value().initialise(*unweightedGR_);
     }
 
     /*
-    // Is the PartialSet already up-to-date?
-    if (DissolveSys::sameString(unweightedSQ_.fingerprint(), std::format("{}/{}", -1), -1))
-    {
-        message("SQ: Unweighted partial S(Q) are up-to-date.\n");
-        return NodeConstants::ProcessResult::Failed;
-    }
-    */
+     * Transform target UnweightedGR into the UnweightedSQ.
+     */
 
     // Transform g(r) into S(Q)
     if (!calculateUnweightedSQ())
         return NodeConstants::ProcessResult::Failed;
 
-    /*
     // Perform averaging of unweighted partials if requested, and if we're not already up-to-date
     if (averagingLength_)
-    {
-        // Store the current fingerprint, since we must ensure we retain it in the averaged data.
-        std::string currentFingerprint{unweightedSQ_.fingerprint()};
+        (*unweightedSQ_) = unweightedSQHistory_.average((*unweightedSQ_), averagingLength_.value().asInteger(),
+                                                        [&]()
+                                                        {
+                                                            PartialSet p;
+                                                            p.initialise(unweightedGR_->realSpeciesPopulations());
+                                                            return p;
+                                                        });
 
-        Averaging::average<PartialSet>(dissolve().processingModuleData(), "UnweightedSQ", name_, averagingLength_.value(),
-                                       averagingScheme_);
-    }
-    */
+    // Save data if requested
+    if (save_ && !unweightedSQ_->save(name(), "UnweightedSQ", "sq", "Q, 1/Angstroms"))
+        return NodeConstants::ProcessResult::Failed;
 
     return NodeConstants::ProcessResult::Success;
 }
diff --git a/src/nodes/sq/sq.cpp b/src/nodes/sq/sq.cpp
index 718422e2c6..a57a41f127 100644
--- a/src/nodes/sq/sq.cpp
+++ b/src/nodes/sq/sq.cpp
@@ -20,8 +20,9 @@ SQNode::SQNode(Graph *parentGraph) : Node(parentGraph)
                                      averagingLength_);
     addOption<Averaging::AveragingScheme>("AveragingScheme", "Weighting scheme to use when averaging partials",
                                           averagingScheme_);
-    addOptionalPointerOutput<PartialSet>("UnweightedSQ", "Unweighted partials for target configuration", unweightedSQ_);
 
+    addOptionalPointerOutput<PartialSet>("UnweightedSQ", "Unweighted partials for target configuration", unweightedSQ_);
+    addOutput<PartialSet *>("UnweightedGR", "Unweighted partials for target configuration", unweightedGR_);
     addOption<bool>("Save", "Whether to save partials to disk after calculation", save_);
 }
 
diff --git a/src/nodes/sq/sq.h b/src/nodes/sq/sq.h
index c95e6f883c..ae6f11ff5e 100644
--- a/src/nodes/sq/sq.h
+++ b/src/nodes/sq/sq.h
@@ -6,6 +6,7 @@
 #include "classes/partialSet.h"
 #include "math/averaging.h"
 #include "math/function1D.h"
+#include "math/history.h"
 #include "math/windowFunction.h"
 #include "module/module.h"
 #include "nodes/graph.h"
@@ -34,6 +35,8 @@ class SQNode : public Node
     PartialSet *unweightedGR_{nullptr};
     // Unweighted S(Q)
     std::optional<PartialSet> unweightedSQ_;
+    // Historical unweighted S(Q)
+    History<PartialSet> unweightedSQHistory_;
     // Number of historical partial sets to combine into final partials
     std::optional<Number> averagingLength_;
     // Weighting scheme to use when averaging partials
@@ -56,7 +59,7 @@ class SQNode : public Node
     /*
      * Functions
      */
-    public:
+    private:
     // Calculate unweighted S(Q) from unweighted g(r)
     bool calculateUnweightedSQ();
 
diff --git a/tests/nodes/graph_argon.cpp b/tests/nodes/graph_argon.cpp
index 85cdb2c141..1be14ea3c9 100644
--- a/tests/nodes/graph_argon.cpp
+++ b/tests/nodes/graph_argon.cpp
@@ -1,11 +1,17 @@
 // SPDX-License-Identifier: GPL-3.0-or-later
 // Copyright (c) 2025 Team Dissolve and contributors
 
+#include "nodes/atomicMC/atomicMC.h"
 #include "nodes/atomicSpecies.h"
 #include "nodes/configuration.h"
 #include "nodes/dissolve.h"
+#include "nodes/energy/energy.h"
+#include "nodes/gr/gr.h"
 #include "nodes/insert.h"
+#include "nodes/md/md.h"
+#include "nodes/neutronSQ/neutronSQ.h"
 #include "nodes/number.h"
+#include "nodes/sq/sq.h"
 #include "tests/testData.h"
 #include <gtest/gtest.h>
 
@@ -14,10 +20,10 @@ namespace UnitTest
 class GraphArgonTest : public ::testing::Test
 {
     public:
-    GraphArgonTest() : dissolve_(coreData_), root_(dissolve_) {}
+    GraphArgonTest() : dissolve_(coreData_), root_(dissolve_) { Node::echo_ = true; }
 
     // Create a graph for testing
-    void createGraph()
+    void createGraph(bool advanced = false)
     {
         /*
          *    Configuration (Bulk)
@@ -45,6 +51,31 @@ class GraphArgonTest : public ::testing::Test
 
         ASSERT_TRUE(root_.addEdge({"Ar", "Species", "Insert", "Species"}));
         ASSERT_TRUE(root_.addEdge({"Bulk", "Configuration", "Insert", "Configuration"}));
+
+        if (advanced)
+        {
+            atomicMCNode_ = dynamic_cast<AtomicMCNode *>(root_.createNode("AtomicMC", "AtomicMC"));
+            mdNode_ = dynamic_cast<MDNode *>(root_.createNode("MD", "MD"));
+            energyNode_ = dynamic_cast<EnergyNode *>(root_.createNode("Energy", "Energy"));
+            grNode_ = dynamic_cast<GRNode *>(root_.createNode("GR", "GR"));
+            sqNode_ = dynamic_cast<SQNode *>(root_.createNode("SQ", "SQ"));
+            neutronSQNode_ = dynamic_cast<NeutronSQNode *>(root_.createNode("NeutronSQ", "NeutronSQ"));
+
+            ASSERT_TRUE(atomicMCNode_);
+            ASSERT_TRUE(mdNode_);
+            ASSERT_TRUE(energyNode_);
+            ASSERT_TRUE(grNode_);
+            ASSERT_TRUE(sqNode_);
+            ASSERT_TRUE(neutronSQNode_);
+
+            ASSERT_TRUE(root_.addEdge({"Insert", "Configuration", "AtomicMC", "Configuration"}));
+            ASSERT_TRUE(root_.addEdge({"AtomicMC", "Configuration", "MD", "Configuration"}));
+            ASSERT_TRUE(root_.addEdge({"MD", "Configuration", "Energy", "Configuration"}));
+            ASSERT_TRUE(root_.addEdge({"Energy", "Configuration", "GR", "Configuration"}));
+            ASSERT_TRUE(root_.addEdge({"GR", "UnweightedGR", "SQ", "UnweightedGR"}));
+            ASSERT_TRUE(root_.addEdge({"SQ", "UnweightedGR", "NeutronSQ", "UnweightedGR"}));
+            ASSERT_TRUE(root_.addEdge({"SQ", "UnweightedSQ", "NeutronSQ", "UnweightedSQ"}));
+        }
     }
 
     protected:
@@ -55,9 +86,15 @@ class GraphArgonTest : public ::testing::Test
     AtomicSpeciesNode *arNode_{nullptr};
     ConfigurationNode *configurationNode_{nullptr};
     InsertNode *insertNode_{nullptr};
+    AtomicMCNode *atomicMCNode_{nullptr};
+    MDNode *mdNode_{nullptr};
+    EnergyNode *energyNode_{nullptr};
+    GRNode *grNode_{nullptr};
+    SQNode *sqNode_{nullptr};
+    NeutronSQNode *neutronSQNode_{nullptr};
 };
 
-TEST_F(GraphArgonTest, Simulation)
+TEST_F(GraphArgonTest, InitSimulation)
 {
     createGraph();
 
@@ -69,4 +106,10 @@ TEST_F(GraphArgonTest, Simulation)
     EXPECT_EQ(cfg->nMolecules(), insertNode_->getInputValue<Number>("Population").asInteger());
 };
 
+TEST_F(GraphArgonTest, AdvancedSimulation)
+{
+    createGraph(true);
+    ASSERT_EQ(neutronSQNode_->run(), NodeConstants::ProcessResult::Success);
+}
+
 } // namespace UnitTest