feat: addition of potentialSet class (#1991)

Co-authored-by: Adam Washington <[email protected]>

src/classes/CMakeLists.txt

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

src/classes/potentialSet.cpp

@@ -0,0 +1,110 @@
+// 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 Potentials
+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::PotentialData> &PotentialSet::potentialMap() { return potentials_; }
+const std::map<std::string, PotentialSet::PotentialData> &PotentialSet::potentialMap() const { return potentials_; }
+
+/*
+ * Operators
+ */
+
+PotentialSet &PotentialSet::operator+=(const double delta)
+{
+    for (auto &[key, pot] : potentials_)
+        pot.potential += delta;
+    return (*this);
+}
+
+PotentialSet &PotentialSet::operator+=(const PotentialSet &source)
+{
+    for (auto &[key, pot] : source.potentialMap())
+    {
+        auto it = potentials_.find(key);
+        if (it != potentials_.end())
+            it->second.potential += pot.potential;
+        else
+            potentials_[key] = pot;
+    }
+    return (*this);
+}
+
+PotentialSet &PotentialSet::operator*=(const double factor)
+{
+    for (auto &[key, pot] : potentials_)
+        pot.potential *= 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;
+        PotentialData 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.potential.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.potential.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 Potentials
+class PotentialSet
+{
+    public:
+    PotentialSet();
+    ~PotentialSet();
+
+    /*
+     * Potentials Data
+     */
+    private:
+    // Fingerprint for these potentials
+    std::string fingerprint_;
+    struct PotentialData
+    {
+        Data1D potential;
+        double count{0};
+        std::shared_ptr<AtomType> at1, at2;
+    };
+    // Map of named potentials to data
+    std::map<std::string, PotentialData> potentials_;
+
+    public:
+    // Reset Potentials
+    void reset();
+    // Set new fingerprint
+    void setFingerprint(std::string_view fingerprint);
+    // Return fingerprint of potentials
+    std::string_view fingerprint() const;
+    // Return full map of potentials specified
+    std::map<std::string, PotentialData> &potentialMap();
+    const std::map<std::string, PotentialData> &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"
@@ -184,6 +185,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"
@@ -47,6 +48,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"
@@ -122,6 +123,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

@@ -5,6 +5,7 @@
 
 #include "classes/scatteringMatrix.h"
 #include "generator/generator.h"
+#include "math/averaging.h"
 #include "module/groups.h"
 #include "module/module.h"
 #include <tuple>
@@ -35,6 +36,12 @@ class EPSRManagerModule : public Module
     };
     // Potential scalings
     std::string potentialScalings_;
+    // Number of historical potential sets to combine into final potential
+    std::optional<int> averagingLength_;
+    // Weighting scheme to use when averaging potential
+    Averaging::AveragingScheme averagingScheme_{Averaging::LinearAveraging};
+    // Vector of averaged potentials over multiple iterations
+    std::vector<std::map<std::string, EPData>> averagedPotentialsStore;
 
     /*
      * Functions

src/modules/epsrManager/process.cpp

@@ -53,6 +53,26 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
     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_, ",");
     for (const auto &scaling : scalings)

tests/classes/CMakeLists.txt

@@ -8,6 +8,7 @@ dissolve_add_test(SRC flags.cpp)
 dissolve_add_test(SRC genericList.cpp)
 dissolve_add_test(SRC interactionPotential.cpp)
 dissolve_add_test(SRC neutronWeights.cpp)
+dissolve_add_test(SRC potentialSet.cpp)
 dissolve_add_test(SRC spaceGroups.cpp)
 dissolve_add_test(SRC species.cpp)
 dissolve_add_test(SRC speciesSite.cpp)

tests/classes/potentialSet.cpp

@@ -0,0 +1,66 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "classes/potentialSet.h"
+#include "math/data1D.h"
+#include "tests/testData.h"
+#include <gtest/gtest.h>
+
+namespace UnitTest
+{
+TEST(PotentialSetTest, SimpleAddition)
+{
+    PotentialSet pots;
+    Data1D x;
+    const auto value = 2.0;
+    x.addPoint(1, value);
+    pots.potentialMap()["A-A"].potential = x;
+    pots.potentialMap()["A-B"].potential = x;
+    pots.potentialMap()["A-C"].potential = x;
+
+    pots += pots;
+    EXPECT_EQ(4.0, pots.potentialMap()["A-A"].potential.value(0));
+    EXPECT_EQ(4.0, pots.potentialMap()["A-B"].potential.value(0));
+    EXPECT_EQ(4.0, pots.potentialMap()["A-C"].potential.value(0));
+}
+
+TEST(PotentialSetTest, Multiplication)
+{
+    PotentialSet pots;
+    Data1D x;
+    const auto value = 3.0;
+    x.addPoint(1, value);
+    pots.potentialMap()["A-A"].potential = x;
+    pots.potentialMap()["A-B"].potential = x;
+    pots.potentialMap()["A-C"].potential = x;
+
+    pots *= 2;
+    EXPECT_EQ(6.0, pots.potentialMap()["A-A"].potential.value(0));
+    EXPECT_EQ(6.0, pots.potentialMap()["A-B"].potential.value(0));
+    EXPECT_EQ(6.0, pots.potentialMap()["A-C"].potential.value(0));
+}
+
+TEST(PotentialSetTest, ComplexAddition)
+{
+    PotentialSet pots;
+    PotentialSet pots2;
+    Data1D x;
+    const auto value = 2.0;
+    x.addPoint(1, value);
+    pots.potentialMap()["A-A"].potential = x;
+    pots.potentialMap()["A-B"].potential = x;
+    pots.potentialMap()["A-C"].potential = x;
+
+    pots2.potentialMap()["A-A"].potential = x;
+    pots2.potentialMap()["A-B"].potential = x;
+    pots2.potentialMap()["A-C"].potential = x;
+    pots2.potentialMap()["A-D"].potential = x;
+
+    pots += pots2;
+    EXPECT_EQ(4.0, pots.potentialMap()["A-A"].potential.value(0));
+    EXPECT_EQ(4.0, pots.potentialMap()["A-B"].potential.value(0));
+    EXPECT_EQ(4.0, pots.potentialMap()["A-C"].potential.value(0));
+    EXPECT_EQ(2.0, pots.potentialMap()["A-D"].potential.value(0));
+}
+
+} // namespace UnitTest