feat: more control of averaging potentials (#1970)

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

src/modules/epsrManager/epsrManager.cpp

@@ -22,7 +22,7 @@ EPSRManagerModule::EPSRManagerModule() : Module(ModuleTypes::EPSRManager)
                                  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");
+                                          averagingLength_, 0, 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,6 +3,7 @@
 
 #pragma once
 
+#include "classes/potentialSet.h"
 #include "classes/scatteringMatrix.h"
 #include "generator/generator.h"
 #include "math/averaging.h"
@@ -28,20 +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 potential sets to combine into final potential
-    std::optional<int> averagingLength_;
-    // Weighting scheme to use when averaging potential
+    // Number of historical potentials sets to combine into final potentials
+    std::optional<int> averagingLength_{};
+    // Weighting scheme to use when averaging potentials
     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

@@ -25,53 +25,47 @@ bool EPSRManagerModule::setUp(ModuleContext &moduleContext, Flags<KeywordBase::K
 // Run main processing
 Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
 {
-    auto &moduleData = moduleContext.dissolve().processingModuleData();
+    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");
 
-    std::map<std::string, EPData> potentials;
+    auto &moduleData = moduleContext.dissolve().processingModuleData();
 
-    // Loop over target data
+    // Loop over target data and form summed / averaged potentials
+    PotentialSet newPotentials;
     for (auto *module : target_)
     {
         auto *epsrModule = dynamic_cast<EPSRModule *>(module);
         auto eps = epsrModule->empiricalPotentials();
 
-        for (auto &&[at1, at2, ep] : eps)
+        for (auto &&[at1, at2, potential] : eps)
         {
             auto key = EPSRManagerModule::pairKey(at1, at2);
-            auto keyIt = potentials.find(key);
-            if (keyIt == potentials.end())
-                potentials[key] = {ep, 1, at1, at2};
+            auto keyIt = newPotentials.potentialMap().find(key);
+            if (keyIt == newPotentials.potentialMap().end())
+                newPotentials.potentialMap()[key] = {potential, 1, at1, at2};
             else
             {
-                Interpolator::addInterpolated(ep, potentials[key].ep, 1.0);
-                ++potentials[key].count;
+                Interpolator::addInterpolated(potential, newPotentials.potentialMap()[key].potential, 1.0);
+                ++newPotentials.potentialMap()[key].count;
             }
         }
     }
+    for (auto &&[key, epData] : newPotentials.potentialMap())
+        epData.potential /= newPotentials.potentialMap()[key].count;
 
-    // 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;
+    // Does a PotentialSet already exist for this Configuration?
+    auto originalPotentialsObject = moduleData.realiseIf<PotentialSet>("PotentialSet", name_, GenericItem::InRestartFileFlag);
+    // Set restart equal to changes
+    originalPotentialsObject.first = newPotentials;
+    // Reference to the current potentials
+    auto &currentPotentials = moduleData.realise<PotentialSet>("PotentialSet", name_, GenericItem::InRestartFileFlag);
+    // Average the Potentials
+    if (averagingLength_)
+        Averaging::average<PotentialSet>(moduleContext.dissolve().processingModuleData(), "PotentialSet", name(),
+                                         averagingLength_.value(), averagingScheme_);
 
     // Apply potential scalings
     auto scalings = DissolveSys::splitString(potentialScalings_, ",");
@@ -91,7 +85,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] : newPotentials.potentialMap())
         {
             // Is this potential a match
             if ((DissolveSys::sameWildString(typeA, epData.at1->name()) &&
@@ -100,20 +94,20 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext)
                  DissolveSys::sameWildString(typeA, epData.at2->name())))
             {
                 Messenger::print(" ... matched and scaled potential {}-{}\n", epData.at1->name(), epData.at2->name());
-                epData.ep *= scaleFactor;
+                epData.potential *= scaleFactor;
                 ++count;
             }
         }
         Messenger::print(" ... matched {} potential(s) in total.\n", count);
     }
 
     // Adjust global potentials
-    for (auto &&[key, epData] : potentials)
+    for (auto &&[key, epData] : currentPotentials.potentialMap())
     {
         // Grab pointer to the relevant pair potential (if it exists)
         auto *pp = moduleContext.dissolve().pairPotential(epData.at1, epData.at2);
         if (pp)
-            pp->setAdditionalPotential(epData.ep);
+            pp->setAdditionalPotential(epData.potential);
     }
 
     return ExecutionResult::Success;

tests/classes/potentialSet.cpp

@@ -1,7 +1,8 @@
 // SPDX-License-Identifier: GPL-3.0-or-later
-// Copyright (c) 2024 Team Dissolve and contributors
+// Copyright (c) 2025 Team Dissolve and contributors
 
 #include "classes/potentialSet.h"
+#include "classes/atomType.h"
 #include "math/data1D.h"
 #include "tests/testData.h"
 #include <gtest/gtest.h>
@@ -63,4 +64,152 @@ TEST(PotentialSetTest, ComplexAddition)
     EXPECT_EQ(2.0, pots.potentialMap()["A-D"].potential.value(0));
 }
 
+TEST(PotentialSetTest, Averaging)
+{
+    GenericList moduleData;
+    auto originalPotentialsObject =
+        moduleData.realiseIf<PotentialSet>("PotentialSet", "module", GenericItem::InRestartFileFlag);
+
+    std::string filename{std::format("{}.restart.txt", ::testing::UnitTest::GetInstance()->current_test_info()->name())};
+    // Open the file
+    LineParser parser;
+    ASSERT_TRUE(parser.openOutput(filename));
+
+    Data1D x;
+    Data1D y;
+    const auto value = 2.0;
+    const auto value2 = 4.0;
+    const auto averagingLength = 10;
+
+    auto A = std::make_shared<AtomType>();
+    auto B = std::make_shared<AtomType>();
+    auto C = std::make_shared<AtomType>();
+    auto D = std::make_shared<AtomType>();
+
+    A->setName("A");
+    B->setName("B");
+    C->setName("C");
+    D->setName("D");
+
+    x.addPoint(1, value);
+    y.addPoint(1, value2);
+
+    for (auto n = 0; n <= 2 * averagingLength; n++)
+    {
+        PotentialSet pots;
+        pots.potentialMap()["A-A"].at1 = A;
+        pots.potentialMap()["A-A"].at2 = A;
+        pots.potentialMap()["A-A"].potential = x;
+
+        pots.potentialMap()["A-B"].at1 = A;
+        pots.potentialMap()["A-B"].at2 = B;
+        pots.potentialMap()["A-B"].potential = x;
+
+        pots.potentialMap()["A-C"].at1 = A;
+        pots.potentialMap()["A-C"].at2 = C;
+        pots.potentialMap()["A-C"].potential = y;
+
+        pots.potentialMap()["A-D"].at1 = A;
+        pots.potentialMap()["A-D"].at2 = D;
+        pots.potentialMap()["A-D"].potential = y;
+
+        originalPotentialsObject.first = pots;
+        Averaging::average<PotentialSet>(moduleData, "PotentialSet", "module", averagingLength,
+                                         Averaging::AveragingScheme::LinearAveraging);
+
+        EXPECT_EQ(A, pots.potentialMap()["A-A"].at1);
+        EXPECT_EQ(A, pots.potentialMap()["A-A"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-B"].at1);
+        EXPECT_EQ(B, pots.potentialMap()["A-B"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-C"].at1);
+        EXPECT_EQ(C, pots.potentialMap()["A-C"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-D"].at1);
+        EXPECT_EQ(D, pots.potentialMap()["A-D"].at2);
+
+        EXPECT_EQ(2.0, pots.potentialMap()["A-A"].potential.value(0));
+        EXPECT_EQ(2.0, pots.potentialMap()["A-A"].potential.value(0));
+        EXPECT_EQ(2.0, pots.potentialMap()["A-B"].potential.value(0));
+        EXPECT_EQ(4.0, pots.potentialMap()["A-C"].potential.value(0));
+        EXPECT_EQ(4.0, pots.potentialMap()["A-D"].potential.value(0));
+
+        pots.serialise(parser);
+    }
+}
+
+TEST(PotentialSetTest, Averaging2)
+{
+    GenericList moduleData;
+    auto originalPotentialsObject =
+        moduleData.realiseIf<PotentialSet>("PotentialSet", "module", GenericItem::InRestartFileFlag);
+
+    std::string filename{std::format("{}.restart.txt", ::testing::UnitTest::GetInstance()->current_test_info()->name())};
+
+    // Open the file
+    LineParser parser;
+    if (!parser.openOutput(filename))
+    {
+        parser.closeFiles();
+    }
+
+    Data1D x;
+    Data1D y;
+    const auto value = 2.0;
+    const auto value2 = 4.0;
+    const auto averagingLength = 10;
+
+    auto A = std::make_shared<AtomType>();
+    auto B = std::make_shared<AtomType>();
+    auto C = std::make_shared<AtomType>();
+    auto D = std::make_shared<AtomType>();
+
+    A->setName("A");
+    B->setName("B");
+    C->setName("C");
+    D->setName("D");
+
+    x.addPoint(1, value);
+    y.addPoint(1, value2);
+
+    for (auto n = 0; n <= 2 * averagingLength; n++)
+    {
+        PotentialSet pots;
+        pots.potentialMap()["A-A"].at1 = A;
+        pots.potentialMap()["A-A"].at2 = A;
+        pots.potentialMap()["A-A"].potential = x;
+
+        pots.potentialMap()["A-B"].at1 = A;
+        pots.potentialMap()["A-B"].at2 = B;
+        pots.potentialMap()["A-B"].potential = x;
+
+        pots.potentialMap()["A-C"].at1 = A;
+        pots.potentialMap()["A-C"].at2 = C;
+        pots.potentialMap()["A-C"].potential = y;
+
+        pots.potentialMap()["A-D"].at1 = A;
+        pots.potentialMap()["A-D"].at2 = D;
+        pots.potentialMap()["A-D"].potential = y;
+
+        originalPotentialsObject.first = pots;
+        Averaging::average<PotentialSet>(moduleData, "PotentialSet", "module", averagingLength,
+                                         Averaging::AveragingScheme::LinearAveraging);
+
+        EXPECT_EQ(A, pots.potentialMap()["A-A"].at1);
+        EXPECT_EQ(A, pots.potentialMap()["A-A"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-B"].at1);
+        EXPECT_EQ(B, pots.potentialMap()["A-B"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-C"].at1);
+        EXPECT_EQ(C, pots.potentialMap()["A-C"].at2);
+        EXPECT_EQ(A, pots.potentialMap()["A-D"].at1);
+        EXPECT_EQ(D, pots.potentialMap()["A-D"].at2);
+
+        EXPECT_EQ(2.0, pots.potentialMap()["A-A"].potential.value(0));
+        EXPECT_EQ(2.0, pots.potentialMap()["A-A"].potential.value(0));
+        EXPECT_EQ(2.0, pots.potentialMap()["A-B"].potential.value(0));
+        EXPECT_EQ(4.0, pots.potentialMap()["A-C"].potential.value(0));
+        EXPECT_EQ(4.0, pots.potentialMap()["A-D"].potential.value(0));
+
+        pots.serialise(parser);
+    }
+}
+
 } // namespace UnitTest

web/docs/userguide/modules/forcefield/epsrmanager/_index.md

@@ -0,0 +1,30 @@
+---
+title: EPSR Manager (Module)
+linkTitle: EPSR Manager
+description: Modify and average the empirical potential.
+---
+
+## Overview
+
+The `EPSR Manager` module allows for the modification, by a fixed multiplier, and averaging over a set number of iterations, of empirial potentials that is calculated via the {{< module "EPSR" >}} module. It can also be used to combine the potentials from several {{< module "EPSR" >}} modules.
+
+## Options
+
+### Targets
+
+|Keyword|Arguments|Default|Description|
+|:------|:--:|:-----:|-----------|
+|`Target`|`EPSR`|--|{{< required-label >}}Specifies the EPSR modules on which to operate.|
+
+### Potential Scaling
+
+|Keyword|Arguments|Default|Description|
+|:------|:--:|:-----:|-----------|
+|`PotScale`|`Potential Set`|--|{{< required-label >}}Specifies the Potential Pair and the modifier to use e.g. Si-O=2.0 would increase the Si-O potential by a factor of 2.0.|
+
+### Averaging
+
+|Keyword|Arguments|Default|Description|
+|:------|:--:|:-----:|-----------|
+|`Averaging`|`int`|`OFF`|Number of historical datasets to combine into final potentials|
+|`AveragingScheme`|[`AveragingScheme`]({{< ref "averagingscheme" >}})|`Linear`|Weighting scheme to use when averaging data|