refactor: Argon nodes testing

.github/workflows/_build_and_package.yml

@@ -88,3 +88,4 @@ jobs:
         publishBenchmarks: ${{ inputs.publishBenchmarks }}
         msvcVersion: ${{ inputs.msvcVersion }}
         osxTargetDeploymentVersion: ${{ inputs.osxTargetDeploymentVersion }}
+

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

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);
 

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

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_; }
 

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;
 

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

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)
 {

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)

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

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_; }

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;
 };

src/main/dissolve.h

@@ -99,13 +99,15 @@ 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)
     PairPotential *pairPotential(const AtomType *at1, const AtomType *at2) const;
     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

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)

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;

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"; }

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;