From 7e01ece71f200740556f2e3bf6b0da925f5aa8fb Mon Sep 17 00:00:00 2001 From: Dan Bradley <115464393+Danbr4d@users.noreply.github.com> Date: Fri, 17 Jan 2025 14:56:34 +0000 Subject: [PATCH] feat: addition of potentialSet class (#1991) Co-authored-by: Adam Washington --- src/classes/CMakeLists.txt | 2 + src/classes/potentialSet.cpp | 110 ++++++++++++++++++++++++ src/classes/potentialSet.h | 67 +++++++++++++++ src/items/deserialisers.cpp | 2 + src/items/producers.cpp | 2 + src/items/serialisers.cpp | 2 + src/modules/epsrManager/epsrManager.cpp | 6 ++ src/modules/epsrManager/epsrManager.h | 7 ++ src/modules/epsrManager/process.cpp | 20 +++++ tests/classes/CMakeLists.txt | 1 + tests/classes/potentialSet.cpp | 66 ++++++++++++++ 11 files changed, 285 insertions(+) create mode 100644 src/classes/potentialSet.cpp create mode 100644 src/classes/potentialSet.h create mode 100644 tests/classes/potentialSet.cpp diff --git a/src/classes/CMakeLists.txt b/src/classes/CMakeLists.txt index 4e5754c81c..e1990bdc2b 100644 --- a/src/classes/CMakeLists.txt +++ b/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 diff --git a/src/classes/potentialSet.cpp b/src/classes/potentialSet.cpp new file mode 100644 index 0000000000..590e56ac8e --- /dev/null +++ b/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 &PotentialSet::potentialMap() { return potentials_; } +const std::map &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; +} diff --git a/src/classes/potentialSet.h b/src/classes/potentialSet.h new file mode 100644 index 0000000000..b907f72475 --- /dev/null +++ b/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 at1, at2; + }; + // Map of named potentials to data + std::map 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 &potentialMap(); + const std::map &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; +}; diff --git a/src/items/deserialisers.cpp b/src/items/deserialisers.cpp index f0264a5e38..ab31b6425b 100644 --- a/src/items/deserialisers.cpp +++ b/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(simpleDeserialiseCore); registerDeserialiser(simpleDeserialiseCore); registerDeserialiser(simpleDeserialise); + registerDeserialiser(simpleDeserialiseCore); registerDeserialiser(simpleDeserialise); registerDeserialiser(simpleDeserialise); registerDeserialiser(simpleDeserialise); diff --git a/src/items/producers.cpp b/src/items/producers.cpp index c576d221e6..e3b27902af 100644 --- a/src/items/producers.cpp +++ b/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"); registerProducer("PartialSet"); registerProducer("PartialSetAccumulator"); + registerProducer("PotentialSet"); registerProducer("SampledData1D"); registerProducer("SampledDouble"); registerProducer("SampledVector"); diff --git a/src/items/serialisers.cpp b/src/items/serialisers.cpp index 0244f32ec4..335f1db198 100644 --- a/src/items/serialisers.cpp +++ b/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(simpleSerialise); registerSerialiser(simpleSerialise); registerSerialiser(simpleSerialise); + registerSerialiser(simpleSerialise); registerSerialiser(simpleSerialise); registerSerialiser(simpleSerialise); registerSerialiser(simpleSerialise); diff --git a/src/modules/epsrManager/epsrManager.cpp b/src/modules/epsrManager/epsrManager.cpp index d840202426..23d7809ee2 100644 --- a/src/modules/epsrManager/epsrManager.cpp +++ b/src/modules/epsrManager/epsrManager.cpp @@ -20,4 +20,10 @@ EPSRManagerModule::EPSRManagerModule() : Module(ModuleTypes::EPSRManager) keywords_.setOrganisation("Options", "Potential Scaling"); keywords_.add("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("Averaging", "Number of historical potential sets to combine into final potentials", + averagingLength_, 1, 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 de6d2da606..cabe6d6423 100644 --- a/src/modules/epsrManager/epsrManager.h +++ b/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 @@ -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 averagingLength_; + // Weighting scheme to use when averaging potential + 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 e28310c4c7..0324daff68 100644 --- a/src/modules/epsrManager/process.cpp +++ b/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 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) diff --git a/tests/classes/CMakeLists.txt b/tests/classes/CMakeLists.txt index 6f1d0368e1..62678b4d24 100644 --- a/tests/classes/CMakeLists.txt +++ b/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) diff --git a/tests/classes/potentialSet.cpp b/tests/classes/potentialSet.cpp new file mode 100644 index 0000000000..4793b6a364 --- /dev/null +++ b/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 + +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