feat: more control of averaging potentials

src/classes/CMakeLists.txt

@@ -45,6 +45,7 @@ add_library(
   pairPotential.cpp
   pairPotentialOverride.cpp
   potentialMap.cpp
+  potentialSet.cpp
   pairIterator.cpp
   partialSet.cpp
   partialSetAccumulator.cpp
@@ -109,6 +110,7 @@ add_library(
   potentialMap.h
   partialSet.h
   partialSetAccumulator.h
+  potentialSet.h
   region.h
   regionalDistributor.h
   scatteringMatrix.h

src/classes/potentialSet.cpp

@@ -0,0 +1,104 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "classes/potentialSet.h"
+#include "base/lineParser.h"
+#include "classes/atomType.h"
+#include "classes/box.h"
+#include "classes/configuration.h"
+#include "io/export/data1D.h"
+#include "items/deserialisers.h"
+#include "items/serialisers.h"
+#include "templates/algorithms.h"
+
+PotentialSet::PotentialSet() { fingerprint_ = "NO_FINGERPRINT"; }
+
+PotentialSet::~PotentialSet() { potentials_.clear(); }
+
+// Reset partial arrays
+void PotentialSet::reset()
+{
+    potentials_.clear();
+    fingerprint_ = "NO_FINGERPRINT";
+}
+
+// Set new fingerprint
+void PotentialSet::setFingerprint(std::string_view fingerprint) { fingerprint_ = fingerprint; }
+
+// Return full map of potentials specified
+std::map<std::string, PotentialSet::EPData> &PotentialSet::potentialMap() { return potentials_; }
+const std::map<std::string, PotentialSet::EPData> &PotentialSet::potentialMap() const { return potentials_; }
+
+/*
+ * Operators
+ */
+
+PotentialSet &PotentialSet::operator+=(const double delta)
+{
+    for (auto &[key, potential] : potentials_)
+        potential.ep += delta;
+    return (*this);
+}
+
+PotentialSet &PotentialSet::operator+=(const PotentialSet &source)
+{
+    for (auto &[key, potential] : source.potentialMap())
+        potentials_[key].ep += potential.ep;
+    return (*this);
+}
+
+PotentialSet &PotentialSet::operator*=(const double factor)
+{
+    for (auto &[key, potential] : potentials_)
+        potential.ep *= factor;
+    return (*this);
+}
+
+/*
+ * Serialisation
+ */
+
+// Read data through specified LineParser
+bool PotentialSet::deserialise(LineParser &parser, const CoreData &coreData)
+{
+    if (parser.readNextLine(LineParser::Defaults, fingerprint_) != LineParser::Success)
+        return false;
+
+    if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success)
+        return false;
+    auto size = parser.argli(0);
+    for (auto n = 0; n < size; ++n)
+    {
+        if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success)
+            return false;
+        EPData value;
+        auto key = parser.args(0);
+        value.count = parser.argi(1);
+        value.at1 = coreData.findAtomType(parser.args(2));
+        value.at2 = coreData.findAtomType(parser.args(3));
+
+        if (!value.ep.deserialise(parser))
+            return false;
+
+        potentials_[key] = value;
+    }
+
+    return true;
+}
+
+// Write data through specified LineParser
+bool PotentialSet::serialise(LineParser &parser) const
+{
+    if (!parser.writeLineF("{}\n", fingerprint_))
+        return false;
+    if (!parser.writeLineF("{}\n", potentials_.size()))
+        return false;
+    for (auto &[key, value] : potentials_)
+    {
+        if (!parser.writeLineF("{} {} {} {}\n", key, value.count, value.at1->name(), value.at2->name()))
+            return false;
+        if (!value.ep.serialise(parser))
+            return false;
+    }
+    return true;
+}

src/classes/potentialSet.h

@@ -0,0 +1,67 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#pragma once
+
+#include "classes/atomTypeMix.h"
+#include "classes/neutronWeights.h"
+#include "math/data1D.h"
+#include "math/histogram1D.h"
+#include "modules/epsr/epsr.h"
+#include "modules/epsrManager/epsrManager.h"
+#include "templates/array2D.h"
+
+// Forward Declarations
+class Configuration;
+class Interpolator;
+
+// Set of Partials
+class PotentialSet
+{
+    public:
+    PotentialSet();
+    ~PotentialSet();
+
+    /*
+     * Potentials Data
+     */
+    private:
+    // Fingerprint for these partials (e.g. reflecting Configuration indices at which they were calculated)
+    std::string fingerprint_;
+    struct EPData
+    {
+        Data1D ep;
+        double count{0};
+        std::shared_ptr<AtomType> at1, at2;
+    };
+    // Pair matrix, containing full atom-atom partial
+    std::map<std::string, EPData> potentials_;
+
+    public:
+    // Reset partial arrays
+    void reset();
+    // Set new fingerprint
+    void setFingerprint(std::string_view fingerprint);
+    // Return fingerprint of partials
+    std::string_view fingerprint() const;
+    // Return full atom-atom partial specified
+    std::map<std::string, EPData> &potentialMap();
+    const std::map<std::string, EPData> &potentialMap() const;
+
+    /*
+     * Operators
+     */
+    public:
+    PotentialSet &operator+=(const double delta);
+    PotentialSet &operator+=(const PotentialSet &source);
+    PotentialSet &operator*=(const double factor);
+
+    /*
+     * Serialisation
+     */
+    public:
+    // Read data through specified LineParser
+    bool deserialise(LineParser &parser, const CoreData &coreData);
+    // Write data through specified LineParser
+    bool serialise(LineParser &parser) const;
+};

src/items/deserialisers.cpp

@@ -8,6 +8,7 @@
 #include "classes/neutronWeights.h"
 #include "classes/partialSet.h"
 #include "classes/partialSetAccumulator.h"
+#include "classes/potentialSet.h"
 #include "classes/xRayWeights.h"
 #include "items/legacy.h"
 #include "math/data1D.h"
@@ -168,6 +169,7 @@ GenericItemDeserialiser::GenericItemDeserialiser()
     registerDeserialiser<NeutronWeights>(simpleDeserialiseCore<NeutronWeights>);
     registerDeserialiser<PartialSet>(simpleDeserialiseCore<PartialSet>);
     registerDeserialiser<PartialSetAccumulator>(simpleDeserialise<PartialSetAccumulator>);
+    registerDeserialiser<PotentialSet>(simpleDeserialiseCore<PotentialSet>);
     registerDeserialiser<SampledData1D>(simpleDeserialise<SampledData1D>);
     registerDeserialiser<SampledDouble>(simpleDeserialise<SampledDouble>);
     registerDeserialiser<SampledVector>(simpleDeserialise<SampledVector>);

src/items/producers.cpp

@@ -7,6 +7,7 @@
 #include "classes/neutronWeights.h"
 #include "classes/partialSet.h"
 #include "classes/partialSetAccumulator.h"
+#include "classes/potentialSet.h"
 #include "classes/xRayWeights.h"
 #include "items/legacy.h"
 #include "math/data1D.h"
@@ -46,6 +47,7 @@ GenericItemProducer::GenericItemProducer()
     registerProducer<NeutronWeights>("NeutronWeights");
     registerProducer<PartialSet>("PartialSet");
     registerProducer<PartialSetAccumulator>("PartialSetAccumulator");
+    registerProducer<PotentialSet>("PotentialSet");
     registerProducer<SampledData1D>("SampledData1D");
     registerProducer<SampledDouble>("SampledDouble");
     registerProducer<SampledVector>("SampledVector");

src/items/serialisers.cpp

@@ -8,6 +8,7 @@
 #include "classes/neutronWeights.h"
 #include "classes/partialSet.h"
 #include "classes/partialSetAccumulator.h"
+#include "classes/potentialSet.h"
 #include "classes/xRayWeights.h"
 #include "math/data1D.h"
 #include "math/data2D.h"
@@ -111,6 +112,7 @@ GenericItemSerialiser::GenericItemSerialiser()
     registerSerialiser<NeutronWeights>(simpleSerialise<NeutronWeights>);
     registerSerialiser<PartialSet>(simpleSerialise<PartialSet>);
     registerSerialiser<PartialSetAccumulator>(simpleSerialise<PartialSetAccumulator>);
+    registerSerialiser<PotentialSet>(simpleSerialise<PotentialSet>);
     registerSerialiser<SampledData1D>(simpleSerialise<SampledData1D>);
     registerSerialiser<SampledDouble>(simpleSerialise<SampledDouble>);
     registerSerialiser<SampledVector>(simpleSerialise<SampledVector>);

src/modules/epsrManager/epsrManager.cpp

@@ -20,4 +20,10 @@ EPSRManagerModule::EPSRManagerModule() : Module(ModuleTypes::EPSRManager)
     keywords_.setOrganisation("Options", "Potential Scaling");
     keywords_.add<StringKeyword>("PotScale", "Comma-separated list of atom type pairs and scaling factors in the form A-B=1.0",
                                  potentialScalings_);
+    keywords_.setOrganisation("Options", "Averaging");
+    keywords_.add<OptionalIntegerKeyword>("Averaging", "Number of historical potential sets to combine into final potentials",
+                                          averagingLength_, 1, std::nullopt, 1, "Off");
+    keywords_.add<EnumOptionsKeyword<Averaging::AveragingScheme>>("AveragingScheme",
+                                                                  "Weighting scheme to use when averaging potentials",
+                                                                  averagingScheme_, Averaging::averagingSchemes());
 }

src/modules/epsrManager/epsrManager.h

@@ -3,8 +3,10 @@
 
 #pragma once
 
+#include "classes/potentialSet.h"
 #include "classes/scatteringMatrix.h"
 #include "generator/generator.h"
+#include "math/averaging.h"
 #include "module/groups.h"
 #include "module/module.h"
 #include <tuple>
@@ -27,14 +29,12 @@ class EPSRManagerModule : public Module
     std::optional<int> modifyPotential_{1};
     // Vector storing atom pairs and associated potentials
     std::vector<std::tuple<std::shared_ptr<AtomType>, std::shared_ptr<AtomType>, Data1D>> potentials_;
-    struct EPData
-    {
-        Data1D ep;
-        double count{0};
-        std::shared_ptr<AtomType> at1, at2;
-    };
     // Potential scalings
     std::string potentialScalings_;
+    // Number of historical partial sets to combine into final partials
+    std::optional<int> averagingLength_;
+    // Weighting scheme to use when averaging partials
+    Averaging::AveragingScheme averagingScheme_{Averaging::LinearAveraging};
 
     /*
      * Functions

src/modules/epsrManager/process.cpp

@@ -25,9 +25,19 @@ bool EPSRManagerModule::setUp(ModuleContext &moduleContext, Flags<KeywordBase::K
 // Run main processing
 Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
 {
+    if (averagingLength_)
+        Messenger::print("Potentials will be averaged over {} sets (scheme = {}).\n", averagingLength_.value(),
+                         Averaging::averagingSchemes().keyword(averagingScheme_));
+    else
+        Messenger::print("Potentials: No averaging of potentials will be performed.\n");
+
     auto &moduleData = moduleContext.dissolve().processingModuleData();
 
-    std::map<std::string, EPData> potentials;
+    // Does a PotentialSet already exist for this Configuration?
+    auto originalPotentialsObject =
+        moduleData.realiseIf<PotentialSet>(fmt::format("AveragingPotentials"), name_, GenericItem::InRestartFileFlag);
+
+    PotentialSet potentials = originalPotentialsObject.first;
 
     // Loop over target data
     for (auto *module : target_)
@@ -38,20 +48,45 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
         for (auto &&[at1, at2, ep] : eps)
         {
             auto key = EPSRManagerModule::pairKey(at1, at2);
-            auto keyIt = potentials.find(key);
-            if (keyIt == potentials.end())
-                potentials[key] = {ep, 1, at1, at2};
+            auto keyIt = potentials.potentialMap().find(key);
+            if (keyIt == potentials.potentialMap().end())
+                potentials.potentialMap()[key] = {ep, 1, at1, at2};
             else
             {
-                Interpolator::addInterpolated(ep, potentials[key].ep, 1.0);
-                ++potentials[key].count;
+                Interpolator::addInterpolated(ep, potentials.potentialMap()[key].ep, 1.0);
+                ++potentials.potentialMap()[key].count;
             }
         }
     }
 
-    // Form averages
-    for (auto &&[key, epData] : potentials)
-        epData.ep /= epData.count;
+    // Perform averaging of potentials data if requested
+    if (averagingLength_)
+        Averaging::average<PotentialSet>(moduleContext.dissolve().processingModuleData(), "AveragingPotentials", name(),
+                                         averagingLength_.value(), averagingScheme_);
+
+    /*    // Form averages
+       for (auto &&[key, epData] : potentials)
+           epData.ep /= epData.count;
+
+       std::map<std::string, EPData> averagedPotentials = potentials;
+
+       averagedPotentialsStore.emplace_back(potentials);
+       // Check if ran the right amount of iterations before averaging
+       if (averagedPotentialsStore.size() > averagingLength_)
+       {
+           averagedPotentialsStore.pop_back();
+       }
+
+       // Average the potentials and replace the map with the new averaged
+       for (const auto &pots : averagedPotentialsStore)
+       {
+           for (auto &&[key, epData] : pots)
+           {
+               averagedPotentials[key].ep += epData.ep;
+               averagedPotentials[key].ep /= averagingLength_.value();
+           }
+       }
+       potentials = averagedPotentials; */
 
     // Apply potential scalings
     auto scalings = DissolveSys::splitString(potentialScalings_, ",");
@@ -71,7 +106,7 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
 
         Messenger::print("Apply scaling factor of {} to potential(s) {}-{}...\n", scaleFactor, typeA, typeB);
         auto count = 0;
-        for (auto &&[key, epData] : potentials)
+        for (auto &&[key, epData] : potentials.potentialMap())
         {
             // Is this potential a match
             if ((DissolveSys::sameWildString(typeA, epData.at1->name()) &&
@@ -88,7 +123,7 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
     }
 
     // Adjust global potentials
-    for (auto &&[key, epData] : potentials)
+    for (auto &&[key, epData] : potentials.potentialMap())
     {
         // Grab pointer to the relevant pair potential (if it exists)
         auto *pp = moduleContext.dissolve().pairPotential(epData.at1, epData.at2);