diff --git a/src/classes/potentialSet.cpp b/src/classes/potentialSet.cpp index 590e56ac8e..d86bcbb284 100644 --- a/src/classes/potentialSet.cpp +++ b/src/classes/potentialSet.cpp @@ -1,5 +1,5 @@ // 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 "base/lineParser.h" diff --git a/src/classes/potentialSet.h b/src/classes/potentialSet.h index b907f72475..665c1b3109 100644 --- a/src/classes/potentialSet.h +++ b/src/classes/potentialSet.h @@ -1,5 +1,5 @@ // SPDX-License-Identifier: GPL-3.0-or-later -// Copyright (c) 2024 Team Dissolve and contributors +// Copyright (c) 2025 Team Dissolve and contributors #pragma once diff --git a/src/modules/epsrManager/epsrManager.cpp b/src/modules/epsrManager/epsrManager.cpp index 23d7809ee2..7af22f6356 100644 --- a/src/modules/epsrManager/epsrManager.cpp +++ b/src/modules/epsrManager/epsrManager.cpp @@ -22,7 +22,7 @@ EPSRManagerModule::EPSRManagerModule() : Module(ModuleTypes::EPSRManager) potentialScalings_); keywords_.setOrganisation("Options", "Averaging"); keywords_.add("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>("AveragingScheme", "Weighting scheme to use when averaging potentials", averagingScheme_, Averaging::averagingSchemes()); diff --git a/src/modules/epsrManager/epsrManager.h b/src/modules/epsrManager/epsrManager.h index cabe6d6423..000f72113c 100644 --- a/src/modules/epsrManager/epsrManager.h +++ b/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 modifyPotential_{1}; // Vector storing atom pairs and associated potentials std::vector, std::shared_ptr, Data1D>> potentials_; - struct EPData - { - Data1D ep; - double count{0}; - std::shared_ptr at1, at2; - }; // Potential scalings std::string potentialScalings_; - // Number of historical potential sets to combine into final potential - std::optional averagingLength_; - // Weighting scheme to use when averaging potential + // Number of historical partial sets to combine into final partials + std::optional averagingLength_{}; + // Weighting scheme to use when averaging partials Averaging::AveragingScheme averagingScheme_{Averaging::LinearAveraging}; - // Vector of averaged potentials over multiple iterations - std::vector> averagedPotentialsStore; /* * Functions diff --git a/src/modules/epsrManager/process.cpp b/src/modules/epsrManager/process.cpp index 0324daff68..4dfb35ad3d 100644 --- a/src/modules/epsrManager/process.cpp +++ b/src/modules/epsrManager/process.cpp @@ -25,53 +25,49 @@ bool EPSRManagerModule::setUp(ModuleContext &moduleContext, Flags 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(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 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(fmt::format("PotentialSet"), name_, GenericItem::InRestartFileFlag); + // Set restart equal to changes + originalPotentialsObject.first = newPotentials; + // Reference to the current potentials + auto ¤tPotentials = + moduleData.realise(fmt::format("PotentialSet"), name_, GenericItem::InRestartFileFlag); + // Average the Potentials + if (averagingLength_) + Averaging::average(moduleContext.dissolve().processingModuleData(), "PotentialSet", name(), + averagingLength_.value(), averagingScheme_); // Apply potential scalings auto scalings = DissolveSys::splitString(potentialScalings_, ","); @@ -91,7 +87,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,7 +96,7 @@ 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; } } @@ -108,12 +104,12 @@ Module::ExecutionResult EPSRManagerModule::process(ModuleContext &moduleContext) } // 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; diff --git a/tests/classes/potentialSet.cpp b/tests/classes/potentialSet.cpp index 4793b6a364..32c028fabb 100644 --- a/tests/classes/potentialSet.cpp +++ b/tests/classes/potentialSet.cpp @@ -2,6 +2,7 @@ // Copyright (c) 2024 Team Dissolve and contributors #include "classes/potentialSet.h" +#include "classes/atomType.h" #include "math/data1D.h" #include "tests/testData.h" #include @@ -63,4 +64,156 @@ TEST(PotentialSetTest, ComplexAddition) EXPECT_EQ(2.0, pots.potentialMap()["A-D"].potential.value(0)); } +TEST(PotentialSetTest, Averaging) +{ + GenericList moduleData; + auto originalPotentialsObject = + moduleData.realiseIf(fmt::format("PotentialSet"), "module", GenericItem::InRestartFileFlag); + + std::string filename_{"test_restart.txt"}; + // 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(); + auto B = std::make_shared(); + auto C = std::make_shared(); + auto D = std::make_shared(); + + 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; + fmt::print("Run {}\n", n); + Averaging::average(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(fmt::format("PotentialSet"), "module", GenericItem::InRestartFileFlag); + + std::string filename_{"test_restart.txt"}; + // 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(); + auto B = std::make_shared(); + auto C = std::make_shared(); + auto D = std::make_shared(); + + 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; + fmt::print("Run {}\n", n); + Averaging::average(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 diff --git a/web/docs/userguide/modules/forcefield/epsrmanager/_index.md b/web/docs/userguide/modules/forcefield/epsrmanager/_index.md new file mode 100644 index 0000000000..3cef199378 --- /dev/null +++ b/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 the empirial potential that is calculated via the `EPSR` module. + +## 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|