diff --git a/src/classes/CMakeLists.txt b/src/classes/CMakeLists.txt index d41b8fa4c6..d73ffce036 100644 --- a/src/classes/CMakeLists.txt +++ b/src/classes/CMakeLists.txt @@ -33,6 +33,7 @@ add_library( distributor.cpp empiricalFormula.cpp fullPairIterator.cpp + histogramSet.cpp isotopeData.cpp isotopologue.cpp isotopologues.cpp @@ -96,6 +97,7 @@ add_library( dataSource.h distributor.h empiricalFormula.h + histogramSet.h interactionPotential.h isotopeData.h isotopologue.h diff --git a/src/classes/histogramSet.cpp b/src/classes/histogramSet.cpp new file mode 100644 index 0000000000..26bd223e77 --- /dev/null +++ b/src/classes/histogramSet.cpp @@ -0,0 +1,128 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#include "classes/histogramSet.h" +#include "base/lineParser.h" +#include "classes/atomType.h" +#include "items/deserialisers.h" +#include "math/mathFunc.h" +#include "templates/algorithms.h" + +HistogramSet::~HistogramSet() +{ + fullHistograms_.clear(); + boundHistograms_.clear(); + unboundHistograms_.clear(); +} + +/* + * Data + */ + +// Initialise histograms +void HistogramSet::initialise(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth) +{ + atomTypeMix_ = atomTypeMix; + + auto nTypes = atomTypeMix_.nItems(); + + fullHistograms_.initialise(nTypes, nTypes, half_); + boundHistograms_.initialise(nTypes, nTypes, half_); + unboundHistograms_.initialise(nTypes, nTypes, half_); + + dissolve::for_each_pair( + ParallelPolicies::par, nTypes, + [&](int i, int j) + { + fullHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); + boundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); + unboundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); + }, + half_); +} + +// Clear all histogram data +void HistogramSet::clear() +{ + fullHistograms_.clear(); + boundHistograms_.clear(); + unboundHistograms_.clear(); +} + +// Zero histogram bins +void HistogramSet::zeroBins() +{ + for (auto &histo : fullHistograms_.linearArray()) + histo.zeroBins(); + for (auto &histo : boundHistograms_.linearArray()) + histo.zeroBins(); + for (auto &histo : unboundHistograms_.linearArray()) + histo.zeroBins(); +} + +// Return atom types list +const AtomTypeMix &HistogramSet::atomTypeMix() const { return atomTypeMix_; } + +// Set new fingerprint +void HistogramSet::setFingerprint(std::string_view fingerprint) { fingerprint_ = fingerprint; } + +// Return fingerprint of partials +std::string_view HistogramSet::fingerprint() const { return fingerprint_; } + +// Return full histogram specified +Histogram1D &HistogramSet::fullHistogram(int i, int j) { return fullHistograms_[{i, j}]; } + +// Return bound histogram specified +Histogram1D &HistogramSet::boundHistogram(int i, int j) { return boundHistograms_[{i, j}]; } + +// Return unbound histogram specified +Histogram1D &HistogramSet::unboundHistogram(int i, int j) { return unboundHistograms_[{i, j}]; } + +/* + * Manipulation + */ + +// Create partials from stored Histogram data +void HistogramSet::formPartials(PartialSet &partials, double boxVolume) +{ + dissolve::for_each_pair( + ParallelPolicies::seq, atomTypeMix_, + [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2) + { + // Calculate RDFs from histogram data + calculateRDF(partials.partial(n, m), fullHistograms_[{n, m}], boxVolume, at1.population(), at2.population(), + &at1 == &at2 ? 2.0 : 1.0); + calculateRDF(partials.boundPartial(n, m), boundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(), + &at1 == &at2 ? 2.0 : 1.0); + calculateRDF(partials.unboundPartial(n, m), unboundHistograms_[{n, m}], boxVolume, at1.population(), + at2.population(), &at1 == &at2 ? 2.0 : 1.0); + + // Set flags for bound partials specifying if they are empty (i.e. there are no + // contributions of that type) + partials.emptyBoundPartial(n, m) = boundHistograms_[{n, m}].nBinned() == 0; + }, + half_); +} + +// Calculate RDF from supplied Histogram and normalisation data +void HistogramSet::calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres, + int nSurrounding, double multiplier) +{ + auto nBins = histogram.nBins(); + double delta = histogram.binWidth(); + const auto &bins = histogram.bins(); + + destination.clear(); + + double shellVolume, factor, r = 0.5 * delta, lowerShellLimit = 0.0, numberDensity = nSurrounding / boxVolume; + for (auto n = 0; n < nBins; ++n) + { + shellVolume = (4.0 / 3.0) * M_PI * (pow(lowerShellLimit + delta, 3.0) - pow(lowerShellLimit, 3.0)); + factor = nCentres * (shellVolume * numberDensity); + + destination.addPoint(r, bins[n] * (multiplier / factor)); + + r += delta; + lowerShellLimit += delta; + } +} diff --git a/src/classes/histogramSet.h b/src/classes/histogramSet.h new file mode 100644 index 0000000000..37526a130f --- /dev/null +++ b/src/classes/histogramSet.h @@ -0,0 +1,64 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#pragma once + +#include "classes/atomTypeMix.h" +#include "classes/partialSet.h" +#include "math/histogram1D.h" +#include "templates/array2D.h" + +// Set of Histograms +class HistogramSet +{ + public: + HistogramSet() = default; + ~HistogramSet(); + + /* + * Data + */ + private: + // AtomTypeMix used to generate matrices + AtomTypeMix atomTypeMix_; + // Fingerprint for these partials (e.g. reflecting Configuration indices at which they were calculated) + std::string fingerprint_; + // Histograms used for calculating full atom-atom partials in r + Array2D fullHistograms_; + // Histograms used for calculating bound atom-atom partials in r + Array2D boundHistograms_; + // Histograms used for deriving unbound atom-atom partials in r + Array2D unboundHistograms_; + // Check for full or half matrix + bool half_{true}; + + public: + // Set up histograms + void initialise(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth); + // Clear all histogram data + void clear(); + // Zero histogram bins + void zeroBins(); + // Return atom types mis + const AtomTypeMix &atomTypeMix() const; + // Set new fingerprint + void setFingerprint(std::string_view fingerprint); + // Return fingerprint of partials + std::string_view fingerprint() const; + // Return full histogram specified + Histogram1D &fullHistogram(int i, int j); + // Return bound histogram specified + Histogram1D &boundHistogram(int i, int j); + // Return unbound histogram specified + Histogram1D &unboundHistogram(int i, int j); + + /* + * Manipulation + */ + public: + // Form partials from stored Histogram data + void formPartials(PartialSet &partials, double boxVolume); + // Calculate RDF from supplied Histogram and normalisation data + static void calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres, + int nSurrounding, double multiplier); +}; diff --git a/src/classes/partialSet.cpp b/src/classes/partialSet.cpp index 476bbff605..2b7eb0f9c5 100644 --- a/src/classes/partialSet.cpp +++ b/src/classes/partialSet.cpp @@ -19,10 +19,6 @@ PartialSet::PartialSet(const std::map &realSpeciesPopul PartialSet::~PartialSet() { - fullHistograms_.clear(); - boundHistograms_.clear(); - unboundHistograms_.clear(); - partials_.clear(); boundPartials_.clear(); emptyBoundPartials_.clear(); @@ -33,21 +29,8 @@ PartialSet::~PartialSet() * Set of Partials */ -// Set up PartialSet -bool PartialSet::setUp(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth) -{ - // Set up partial arrays - if (!setUpPartials(atomTypeMix, half_)) - return false; - - // Initialise histograms for g(r) calculation - setUpHistograms(rdfRange, binWidth); - - return true; -} - -// Set up PartialSet without initialising histogram arrays -bool PartialSet::setUpPartials(const AtomTypeMix &atomTypeMix, bool half) +// Initialise +void PartialSet::initialise(const AtomTypeMix &atomTypeMix, bool half) { // Copy type array atomTypeMix_ = atomTypeMix; @@ -78,42 +61,11 @@ bool PartialSet::setUpPartials(const AtomTypeMix &atomTypeMix, bool half) total_.clear(); boundTotal_.clear(); unboundTotal_.clear(); - - return true; -} - -// Set up histogram arrays for g(r) calculation -void PartialSet::setUpHistograms(double rdfRange, double binWidth) -{ - auto nTypes = atomTypeMix_.nItems(); - - fullHistograms_.initialise(nTypes, nTypes, half_); - boundHistograms_.initialise(nTypes, nTypes, half_); - unboundHistograms_.initialise(nTypes, nTypes, half_); - - dissolve::for_each_pair( - ParallelPolicies::par, nTypes, - [&](int i, int j) - { - fullHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); - boundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); - unboundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth); - }, - half_); } // Reset partial arrays void PartialSet::reset() { - // Zero histogram bins if present - for (auto n = 0; n < fullHistograms_.nRows(); ++n) - for (auto m = n; m < fullHistograms_.nColumns(); ++m) - { - fullHistograms_[{n, m}].zeroBins(); - boundHistograms_[{n, m}].zeroBins(); - unboundHistograms_[{n, m}].zeroBins(); - } - // Zero partials dissolve::for_each_pair( ParallelPolicies::par, atomTypeMix_.nItems(), @@ -144,15 +96,6 @@ void PartialSet::setFingerprint(std::string_view fingerprint) { fingerprint_ = f // Return fingerprint of partials std::string_view PartialSet::fingerprint() const { return fingerprint_; } -// Return full histogram specified -Histogram1D &PartialSet::fullHistogram(int i, int j) { return fullHistograms_[{i, j}]; } - -// Return bound histogram specified -Histogram1D &PartialSet::boundHistogram(int i, int j) { return boundHistograms_[{i, j}]; } - -// Return unbound histogram specified -Histogram1D &PartialSet::unboundHistogram(int i, int j) { return unboundHistograms_[{i, j}]; } - // Return full atom-atom partial specified Data1D &PartialSet::partial(int i, int j) { return partials_[{i, j}]; } const Data1D &PartialSet::partial(int i, int j) const { return partials_[{i, j}]; } @@ -165,6 +108,10 @@ const Data1D &PartialSet::unboundPartial(int i, int j) const { return unboundPar Data1D &PartialSet::boundPartial(int i, int j) { return boundPartials_[{i, j}]; } const Data1D &PartialSet::boundPartial(int i, int j) const { return boundPartials_[{i, j}]; } +// Return emptyBound flag +char &PartialSet::emptyBoundPartial(int i, int j) { return emptyBoundPartials_[{i, j}]; } +const char &PartialSet::emptyBoundPartial(int i, int j) const { return emptyBoundPartials_[{i, j}]; } + // Return whether specified bound partial is empty bool PartialSet::isBoundPartialEmpty(int i, int j) const { return emptyBoundPartials_[{i, j}]; } @@ -344,28 +291,6 @@ void PartialSet::adjust(double delta) unboundTotal_ += delta; } -// Form partials from stored Histogram data -void PartialSet::formPartials(double boxVolume) -{ - dissolve::for_each_pair( - ParallelPolicies::seq, atomTypeMix_, - [&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2) - { - // Calculate RDFs from histogram data - calculateRDF(partials_[{n, m}], fullHistograms_[{n, m}], boxVolume, at1.population(), at2.population(), - &at1 == &at2 ? 2.0 : 1.0); - calculateRDF(boundPartials_[{n, m}], boundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(), - &at1 == &at2 ? 2.0 : 1.0); - calculateRDF(unboundPartials_[{n, m}], unboundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(), - &at1 == &at2 ? 2.0 : 1.0); - - // Set flags for bound partials specifying if they are empty (i.e. there are no - // contributions of that type) - emptyBoundPartials_[{n, m}] = boundHistograms_[{n, m}].nBinned() == 0; - }, - half_); -} - // Add in partials from source PartialSet to our own bool PartialSet::addPartials(PartialSet &source, double weighting) { @@ -410,29 +335,6 @@ bool PartialSet::addPartials(PartialSet &source, double weighting) return true; } -// Calculate RDF from supplied Histogram and normalisation data -void PartialSet::calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres, - int nSurrounding, double multiplier) -{ - auto nBins = histogram.nBins(); - double delta = histogram.binWidth(); - const auto &bins = histogram.bins(); - - destination.clear(); - - double shellVolume, factor, r = 0.5 * delta, lowerShellLimit = 0.0, numberDensity = nSurrounding / boxVolume; - for (auto n = 0; n < nBins; ++n) - { - shellVolume = (4.0 / 3.0) * M_PI * (pow(lowerShellLimit + delta, 3.0) - pow(lowerShellLimit, 3.0)); - factor = nCentres * (shellVolume * numberDensity); - - destination.addPoint(r, bins[n] * (multiplier / factor)); - - r += delta; - lowerShellLimit += delta; - } -} - /* * Operators */ diff --git a/src/classes/partialSet.h b/src/classes/partialSet.h index 0324abad89..96a5ceb253 100644 --- a/src/classes/partialSet.h +++ b/src/classes/partialSet.h @@ -6,7 +6,6 @@ #include "classes/atomTypeMix.h" #include "classes/neutronWeights.h" #include "math/data1D.h" -#include "math/histogram1D.h" #include "templates/array2D.h" // Forward Declarations @@ -29,12 +28,6 @@ class PartialSet AtomTypeMix atomTypeMix_; // Fingerprint for these partials (e.g. reflecting Configuration indices at which they were calculated) std::string fingerprint_; - // Histograms used for calculating full atom-atom partials in r - Array2D fullHistograms_; - // Histograms used for calculating bound atom-atom partials in r - Array2D boundHistograms_; - // Histograms used for deriving unbound atom-atom partials in r - Array2D unboundHistograms_; // Pair matrix, containing full atom-atom partial Array2D partials_; // Unbound matrix, containing atom-atom partial of unbound pairs @@ -53,12 +46,8 @@ class PartialSet std::map realSpeciesPopulations_; public: - // Set up PartialSet, including initialising histograms for g(r) use - bool setUp(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth); - // Set up PartialSet without initialising histogram arrays - bool setUpPartials(const AtomTypeMix &atomTypMix, bool half = true); - // Set up histogram arrays for g(r) calculation - void setUpHistograms(double rdfRange, double binWidth); + // Initialise + void initialise(const AtomTypeMix &atomTypMix, bool half = true); // Reset partial arrays void reset(); // Return number of AtomTypes used to generate matrices @@ -69,12 +58,6 @@ class PartialSet void setFingerprint(std::string_view fingerprint); // Return fingerprint of partials std::string_view fingerprint() const; - // Return full histogram specified - Histogram1D &fullHistogram(int i, int j); - // Return bound histogram specified - Histogram1D &boundHistogram(int i, int j); - // Return unbound histogram specified - Histogram1D &unboundHistogram(int i, int j); // Return full atom-atom partial specified Data1D &partial(int i, int j); const Data1D &partial(int i, int j) const; @@ -84,6 +67,9 @@ class PartialSet // Return atom-atom partial for bound pairs Data1D &boundPartial(int i, int j); const Data1D &boundPartial(int i, int j) const; + // Return empty bound partial flag + char &emptyBoundPartial(int i, int j); + const char &emptyBoundPartial(int i, int j) const; // Return whether specified bound partial is empty bool isBoundPartialEmpty(int i, int j) const; // Sum partials into totals @@ -106,10 +92,6 @@ class PartialSet const std::map &realSpeciesPopulations() const; // Save all partials and total bool save(std::string_view prefix, std::string_view tag, std::string_view suffix, std::string_view abscissaUnits) const; - // Name all object based on the supplied prefix - void setObjectTags(std::string_view prefix, std::string_view suffix = ""); - // Return prefix applied to object names - std::string_view objectNamePrefix() const; /* * Manipulation @@ -119,13 +101,8 @@ class PartialSet void adjust(double delta); public: - // Form partials from stored Histogram data - void formPartials(double boxVolume); // Add in partials from source PartialSet to our own, with specified weighting bool addPartials(PartialSet &source, double weighting); - // Calculate RDF from supplied Histogram and normalisation data - static void calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres, - int nSurrounding, double multiplier); /* * Operators diff --git a/src/items/searchers.cpp b/src/items/searchers.cpp index d4cccf0398..f63a6f7b31 100644 --- a/src/items/searchers.cpp +++ b/src/items/searchers.cpp @@ -9,6 +9,7 @@ #include "math/data2D.h" #include "math/data3D.h" #include "math/sampledData1D.h" +#include "math/sampledDouble.h" /* * Data1D diff --git a/src/modules/gr/functions.cpp b/src/modules/gr/functions.cpp index 29f9850b1a..346fc04f24 100644 --- a/src/modules/gr/functions.cpp +++ b/src/modules/gr/functions.cpp @@ -9,9 +9,6 @@ #include "classes/cell.h" #include "classes/configuration.h" #include "classes/species.h" -#include "classes/speciesAngle.h" -#include "classes/speciesBond.h" -#include "classes/speciesTorsion.h" #include "main/dissolve.h" #include "math/combinations.h" #include "math/error.h" @@ -20,51 +17,37 @@ #include "modules/gr/gr.h" #include "templates/algorithms.h" #include "templates/combinable.h" -#include #include -namespace -{ -void addHistogramsToPartialSet(Array2D &histograms, PartialSet &target) -{ - for (auto k = 0; k < target.nAtomTypes(); ++k) - for (auto j = 0; j < target.nAtomTypes(); ++j) - { - auto &histo = target.fullHistogram(k, j); - histo = std::move(histograms[{k, j}]); - } -} -} // namespace - /* * Private Functions */ // Calculate partial g(r) in serial with simple double-loop -bool GRModule::calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet) +bool GRModule::calculateGRTestSerial(Configuration *cfg) { // Calculate radial distribution functions with a simple double loop, in serial const auto *box = cfg->box(); dissolve::for_each_pair( ParallelPolicies::seq, cfg->atoms(), - [box, &partialSet](auto i, auto &ii, auto j, auto &jj) + [&, box](auto i, auto &ii, auto j, auto &jj) { if (&ii != &jj) - partialSet.fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r())); + histograms_->fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r())); }); return true; } // Calculate partial g(r) with optimised double-loop -bool GRModule::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const double binWidth) +bool GRModule::calculateGRSimple(Configuration *cfg, const double binWidth) { // Variables int n, m, nTypes, typeI, typeJ, i, j, nPoints; // Construct local arrays of atom type positions - nTypes = partialSet.nAtomTypes(); + nTypes = histograms_->atomTypeMix().nItems(); Messenger::printVerbose("Constructing local partial working arrays for {} types.\n", nTypes); const auto *box = cfg->box(); std::vector r(nTypes); @@ -103,9 +86,9 @@ bool GRModule::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, con for (typeI = 0; typeI < nTypes; ++typeI) { ri = r[typeI]; - auto &histogram = partialSet.fullHistogram(typeI, typeI).bins(); + auto &histogram = histograms_->fullHistogram(typeI, typeI).bins(); bins = binss[typeI]; - nPoints = partialSet.fullHistogram(typeI, typeI).nBins(); + nPoints = histograms_->fullHistogram(typeI, typeI).nBins(); PairIterator pairs(maxr[typeI]); std::for_each(pairs.begin(), pairs.end(), [box, bins, rbin, ri, nPoints, &histogram](auto it) @@ -139,9 +122,9 @@ bool GRModule::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, con continue; rj = r[typeJ]; - auto &histogram = partialSet.fullHistogram(typeI, typeJ).bins(); + auto &histogram = histograms_->fullHistogram(typeI, typeJ).bins(); bins = binss[typeJ]; - nPoints = partialSet.fullHistogram(typeI, typeJ).nBins(); + nPoints = histograms_->fullHistogram(typeI, typeJ).nBins(); for (i = 0; i < maxr[typeI]; ++i) { centre = ri[i]; @@ -164,19 +147,19 @@ bool GRModule::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, con return true; } -bool GRModule::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const double rdfRange) +bool GRModule::calculateGRCells(Configuration *cfg, const double rdfRange) { auto &cellArray = cfg->cells(); Combinations comb(cellArray.nCells()); auto combinableHistograms = dissolve::CombinableValue>( - [&partialSet]() + [&]() { Array2D histograms; - histograms.initialise(partialSet.nAtomTypes(), partialSet.nAtomTypes(), true); - for (auto i = 0; i < partialSet.nAtomTypes(); ++i) - for (auto j = i; j < partialSet.nAtomTypes(); ++j) - histograms[{i, j}] = partialSet.fullHistogram(i, j); + histograms.initialise(histograms_->atomTypeMix().nItems(), histograms_->atomTypeMix().nItems(), true); + for (auto i = 0; i < histograms_->atomTypeMix().nItems(); ++i) + for (auto j = i; j < histograms_->atomTypeMix().nItems(); ++j) + histograms[{i, j}] = histograms_->fullHistogram(i, j); return histograms; }); @@ -224,7 +207,11 @@ bool GRModule::calculateGRCells(Configuration *cfg, PartialSet &partialSet, cons dissolve::for_each(ParallelPolicies::par, dissolve::counting_iterator(0), dissolve::counting_iterator(comb.getNumCombinations()), unaryOp); auto histograms = combinableHistograms.finalize(); - addHistogramsToPartialSet(histograms, partialSet); + + // Copy the final calculated full histograms to the HistogramSet + for (auto k = 0; k < histograms_->atomTypeMix().nItems(); ++k) + for (auto j = 0; j < histograms_->atomTypeMix().nItems(); ++j) + histograms_->fullHistogram(k, j) = histograms[{k, j}]; // Atoms within the same cell for (int n = 0; n < cellArray.nCells(); ++n) @@ -236,7 +223,7 @@ bool GRModule::calculateGRCells(Configuration *cfg, PartialSet &partialSet, cons PairIterator pairs(atomsI.size()); std::for_each( pairs.begin(), pairs.end(), - [&atomsI, &partialSet](auto it) + [&](auto it) { auto [idx, jdx] = it; if (idx == jdx) @@ -248,7 +235,7 @@ bool GRModule::calculateGRCells(Configuration *cfg, PartialSet &partialSet, cons if (typeI != AtomType::Ignore && typeJ != AtomType::Ignore) { // No need to perform MIM since we're in the same cell - partialSet.fullHistogram(i->localTypeIndex(), j->localTypeIndex()).bin((i->r() - j->r()).magnitude()); + histograms_->fullHistogram(i->localTypeIndex(), j->localTypeIndex()).bin((i->r() - j->r()).magnitude()); } }); } @@ -321,7 +308,7 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo GenericItem::InRestartFileFlag); auto &originalgr = originalGRObject.first; if (originalGRObject.second == GenericItem::ItemStatus::Created) - originalgr.setUp(cfg->atomTypePopulations(), rdfRange, rdfBinWidth); + originalgr.initialise(cfg->atomTypePopulations()); // Is the PartialSet already up-to-date? // If so, can exit now, *unless* the Test method is requested, in which case we go ahead and calculate anyway @@ -339,24 +326,28 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo /* * Make sure histograms are set up, and reset any existing data */ - - originalgr.setUpHistograms(rdfRange, rdfBinWidth); - originalgr.reset(); + if (!histograms_) + { + histograms_.emplace(); + histograms_->initialise(cfg->atomTypePopulations(), rdfRange, rdfBinWidth); + } + histograms_->zeroBins(); /* * Calculate full (intra+inter) partials */ Timer timer; + originalgr.reset(); if (method == GRModule::TestMethod) - calculateGRTestSerial(cfg, originalgr); + calculateGRTestSerial(cfg); else if (method == GRModule::SimpleMethod) - calculateGRSimple(cfg, originalgr, rdfBinWidth); + calculateGRSimple(cfg, rdfBinWidth); else if (method == GRModule::CellsMethod) - calculateGRCells(cfg, originalgr, rdfRange); + calculateGRCells(cfg, rdfRange); else if (method == GRModule::AutoMethod) { - cfg->nAtoms() > 10000 ? calculateGRCells(cfg, originalgr, rdfRange) : calculateGRSimple(cfg, originalgr, rdfBinWidth); + cfg->nAtoms() > 10000 ? calculateGRCells(cfg, rdfRange) : calculateGRSimple(cfg, rdfBinWidth); } timer.stop(); Messenger::print("Finished calculation of partials ({} elapsed).\n", timer.totalTimeString()); @@ -366,21 +357,14 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo */ const auto *box = cfg->box(); - const auto &cells = cfg->cells(); - - // Set start/stride for parallel loop (pool solo) - auto offset = 0; - auto nChunks = 1; - timer.start(); - // Loop over molecules for (auto &mol : cfg->molecules()) { const auto &atoms = mol->atoms(); dissolve::for_each_pair(ParallelPolicies::seq, atoms, - [box, &originalgr](int index, auto &i, int jndex, auto &j) + [&, box](int index, auto &i, int jndex, auto &j) { // Ignore atom on itself if (index == jndex) @@ -394,7 +378,7 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo if (typeJ == AtomType::Ignore) return; - originalgr.boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); + histograms_->boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); }); } @@ -408,21 +392,21 @@ bool GRModule::calculateGR(GenericList &processingData, Configuration *cfg, GRMo */ timer.start(); - auto success = - for_each_pair_early(originalgr.nAtomTypes(), - [&originalgr](auto typeI, auto typeJ) -> EarlyReturn - { - // Create unbound histogram from total and bound data - originalgr.unboundHistogram(typeI, typeJ) = originalgr.fullHistogram(typeI, typeJ); - originalgr.unboundHistogram(typeI, typeJ).add(originalgr.boundHistogram(typeI, typeJ), -1.0); - - return EarlyReturn::Continue; - }); + auto success = for_each_pair_early( + originalgr.nAtomTypes(), + [&](auto typeI, auto typeJ) -> EarlyReturn + { + // Create unbound histogram from total and bound data + histograms_->unboundHistogram(typeI, typeJ) = histograms_->fullHistogram(typeI, typeJ); + histograms_->unboundHistogram(typeI, typeJ).add(histograms_->boundHistogram(typeI, typeJ), -1.0); + + return EarlyReturn::Continue; + }); if (success.has_value() && !success.value()) return false; - // Transform histogram data into radial distribution functions - originalgr.formPartials(box->volume()); + // Transform histogram data into partials and store + histograms_->formPartials(originalgr, box->volume()); // Sum total functions originalgr.formTotals(true); @@ -511,7 +495,7 @@ bool GRModule::sumUnweightedGR(GenericList &processingData, std::string_view tar combinedAtomTypes.finalise(); // Set up PartialSet container - summedUnweightedGR.setUpPartials(combinedAtomTypes); + summedUnweightedGR.initialise(combinedAtomTypes); // Determine total weighting factors and combined density over all Configurations, and set up a Configuration/weight // Vector for simplicity diff --git a/src/modules/gr/gr.h b/src/modules/gr/gr.h index 52a122f053..d283b96383 100644 --- a/src/modules/gr/gr.h +++ b/src/modules/gr/gr.h @@ -3,6 +3,7 @@ #pragma once +#include "classes/histogramSet.h" #include "classes/partialSet.h" #include "math/averaging.h" #include "math/function1D.h" @@ -59,17 +60,19 @@ class GRModule : public Module bool save_{false}; // Whether to save original (unbroadened) partials and total functions to disk bool saveOriginal_{false}; + // Histograms for RDF calculation + std::optional histograms_; /* * Functions */ private: // Calculate partial g(r) in serial with simple double-loop - bool calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet); + bool calculateGRTestSerial(Configuration *cfg); // Calculate partial g(r) with optimised double-loop - bool calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const double rdfRange); + bool calculateGRSimple(Configuration *cfg, const double rdfRange); // Calculate partial g(r) utilising Cell neighbour lists - bool calculateGRCells(Configuration *cfg, PartialSet &partialSet, const double binWidth); + bool calculateGRCells(Configuration *cfg, const double binWidth); public: // Calculate and return effective density based on target Configurations diff --git a/src/modules/neutronSQ/process.cpp b/src/modules/neutronSQ/process.cpp index 94d12e4304..e0314ca268 100644 --- a/src/modules/neutronSQ/process.cpp +++ b/src/modules/neutronSQ/process.cpp @@ -194,7 +194,7 @@ Module::ExecutionResult NeutronSQModule::process(Dissolve &dissolve) auto [weightedSQ, wSQstatus] = dissolve.processingModuleData().realiseIf("WeightedSQ", name_, GenericItem::InRestartFileFlag); if (wSQstatus == GenericItem::ItemStatus::Created) - weightedSQ.setUpPartials(unweightedSQ.atomTypeMix()); + weightedSQ.initialise(unweightedSQ.atomTypeMix()); // Calculate weighted S(Q) calculateWeightedSQ(unweightedSQ, weightedSQ, weights, normaliseTo_); @@ -220,7 +220,7 @@ Module::ExecutionResult NeutronSQModule::process(Dissolve &dissolve) auto [weightedGR, wGRstatus] = dissolve.processingModuleData().realiseIf("WeightedGR", name_, GenericItem::InRestartFileFlag); if (wGRstatus == GenericItem::ItemStatus::Created) - weightedGR.setUpPartials(unweightedGR.atomTypeMix()); + weightedGR.initialise(unweightedGR.atomTypeMix()); // Calculate weighted g(r) calculateWeightedGR(unweightedGR, weightedGR, weights, normaliseTo_); diff --git a/src/modules/sq/process.cpp b/src/modules/sq/process.cpp index 59948da382..8612703405 100644 --- a/src/modules/sq/process.cpp +++ b/src/modules/sq/process.cpp @@ -89,7 +89,7 @@ Module::ExecutionResult SQModule::process(Dissolve &dissolve) dissolve.processingModuleData().realiseIf("UnweightedSQ", name_, GenericItem::InRestartFileFlag); auto &unweightedsq = uSQObject.first; if (uSQObject.second == GenericItem::ItemStatus::Created) - unweightedsq.setUpPartials(unweightedgr.atomTypeMix()); + unweightedsq.initialise(unweightedgr.atomTypeMix()); // Is the PartialSet already up-to-date? if (DissolveSys::sameString( diff --git a/src/modules/tr/process.cpp b/src/modules/tr/process.cpp index aef7d455c0..862684def9 100644 --- a/src/modules/tr/process.cpp +++ b/src/modules/tr/process.cpp @@ -54,7 +54,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve) // Make weightedGR Partial set PartialSet representativeGR; - representativeGR.setUpPartials(unweightedSQ.atomTypeMix(), false); + representativeGR.initialise(unweightedSQ.atomTypeMix(), false); // Get effective atomic density of underlying g(r) const auto rho = grModule->effectiveDensity(); @@ -63,7 +63,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve) auto [weightedTR, wGRstatus] = dissolve.processingModuleData().realiseIf("WeightedTR", name_, GenericItem::InRestartFileFlag); if (wGRstatus == GenericItem::ItemStatus::Created) - weightedTR.setUpPartials(unweightedGR.atomTypeMix(), false); + weightedTR.initialise(unweightedGR.atomTypeMix(), false); // Get Q-range and window function to use for transformation of reference F(Q) to G(r) auto refftQMin = refQMin_.value_or(0.0); @@ -126,7 +126,7 @@ Module::ExecutionResult TRModule::process(Dissolve &dissolve) auto [representativeTR, rTRstatus] = dissolve.processingModuleData().realiseIf("RepresentativeTR", name_, GenericItem::InRestartFileFlag); if (rTRstatus == GenericItem::ItemStatus::Created) - representativeTR.setUpPartials(representativeGR.atomTypeMix(), false); + representativeTR.initialise(representativeGR.atomTypeMix(), false); dissolve::for_each_pair( ParallelPolicies::par, representativeGR.nAtomTypes(), diff --git a/src/modules/xRaySQ/process.cpp b/src/modules/xRaySQ/process.cpp index a87551c307..9017b2c3e4 100644 --- a/src/modules/xRaySQ/process.cpp +++ b/src/modules/xRaySQ/process.cpp @@ -202,7 +202,7 @@ Module::ExecutionResult XRaySQModule::process(Dissolve &dissolve) auto [weightedSQ, wSQtatus] = dissolve.processingModuleData().realiseIf("WeightedSQ", name_, GenericItem::InRestartFileFlag); if (wSQtatus == GenericItem::ItemStatus::Created) - weightedSQ.setUpPartials(unweightedSQ.atomTypeMix()); + weightedSQ.initialise(unweightedSQ.atomTypeMix()); // Calculate weighted S(Q) calculateWeightedSQ(unweightedSQ, weightedSQ, weights, normaliseTo_); @@ -258,7 +258,7 @@ Module::ExecutionResult XRaySQModule::process(Dissolve &dissolve) auto [weightedGR, wGRstatus] = dissolve.processingModuleData().realiseIf("WeightedGR", name_, GenericItem::InRestartFileFlag); if (wGRstatus == GenericItem::ItemStatus::Created) - weightedGR.setUpPartials(unweightedSQ.atomTypeMix()); + weightedGR.initialise(unweightedSQ.atomTypeMix()); // Calculate weighted g(r) calculateWeightedGR(unweightedGR, weightedGR, weights, normaliseTo_); diff --git a/src/nodes/gr/gr.h b/src/nodes/gr/gr.h index b2cc9da832..90624f4ae1 100644 --- a/src/nodes/gr/gr.h +++ b/src/nodes/gr/gr.h @@ -5,6 +5,7 @@ #include "base/enumOptions.h" #include "classes/configuration.h" +#include "classes/histogramSet.h" #include "classes/partialSet.h" #include "classes/species.h" #include "items/list.h" @@ -68,6 +69,8 @@ class GRNode : public Node bool save_{false}; // Whether to save raw partials and total functions to disk bool saveRaw_{false}; + // Histograms for RDF calculation + std::optional histograms_; /* * Functions diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp index 13f5fbfe62..e70e172a41 100644 --- a/src/nodes/gr/helpers.cpp +++ b/src/nodes/gr/helpers.cpp @@ -5,9 +5,6 @@ #include "classes/atomType.h" #include "classes/box.h" #include "classes/cell.h" -#include "classes/speciesAngle.h" -#include "classes/speciesBond.h" -#include "classes/speciesTorsion.h" #include "main/dissolve.h" #include "math/combinations.h" #include "math/error.h" @@ -16,22 +13,8 @@ #include "nodes/gr/gr.h" #include "templates/algorithms.h" #include "templates/combinable.h" -#include #include -namespace -{ -void addHistogramsToPartialSet(Array2D &histograms, PartialSet &target) -{ - for (auto k = 0; k < target.nAtomTypes(); ++k) - for (auto j = 0; j < target.nAtomTypes(); ++j) - { - auto &histo = target.fullHistogram(k, j); - histo = std::move(histograms[{k, j}]); - } -} -} // namespace - /* * Private Functions */ @@ -47,7 +30,7 @@ bool GRNode::calculateGRTestSerial() [&, box](auto i, auto &ii, auto j, auto &jj) { if (&ii != &jj) - rawGR_->fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r())); + histograms_->fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r())); }); return true; @@ -99,9 +82,9 @@ bool GRNode::calculateGRSimple() for (typeI = 0; typeI < nTypes; ++typeI) { ri = r[typeI]; - auto &histogram = rawGR_->fullHistogram(typeI, typeI).bins(); + auto &histogram = histograms_->fullHistogram(typeI, typeI).bins(); bins = binss[typeI]; - nPoints = rawGR_->fullHistogram(typeI, typeI).nBins(); + nPoints = histograms_->fullHistogram(typeI, typeI).nBins(); PairIterator pairs(maxr[typeI]); std::for_each(pairs.begin(), pairs.end(), [box, bins, rbin, ri, nPoints, &histogram](auto it) @@ -135,9 +118,9 @@ bool GRNode::calculateGRSimple() continue; rj = r[typeJ]; - auto &histogram = rawGR_->fullHistogram(typeI, typeJ).bins(); + auto &histogram = histograms_->fullHistogram(typeI, typeJ).bins(); bins = binss[typeJ]; - nPoints = rawGR_->fullHistogram(typeI, typeJ).nBins(); + nPoints = histograms_->fullHistogram(typeI, typeJ).nBins(); for (i = 0; i < maxr[typeI]; ++i) { centre = ri[i]; @@ -175,7 +158,7 @@ bool GRNode::calculateGRCells(double grRange) histograms.initialise(rawGR_->nAtomTypes(), rawGR_->nAtomTypes(), true); for (auto i = 0; i < rawGR_->nAtomTypes(); ++i) for (auto j = i; j < rawGR_->nAtomTypes(); ++j) - histograms[{i, j}] = rawGR_->fullHistogram(i, j); + histograms[{i, j}] = histograms_->fullHistogram(i, j); return histograms; }); @@ -223,7 +206,15 @@ bool GRNode::calculateGRCells(double grRange) dissolve::for_each(ParallelPolicies::par, dissolve::counting_iterator(0), dissolve::counting_iterator(comb.getNumCombinations()), unaryOp); auto histograms = combinableHistograms.finalize(); - addHistogramsToPartialSet(histograms, *rawGR_); + // Copy the final calculated full histograms to the HistogramSet + for (auto k = 0; k < histograms_->atomTypeMix().nItems(); ++k) + for (auto j = 0; j < histograms_->atomTypeMix().nItems(); ++j) + histograms_->fullHistogram(k, j) = histograms[{k, j}]; + + // Copy the final calculated full histograms to the HistogramSet + for (auto k = 0; k < histograms_->atomTypeMix().nItems(); ++k) + for (auto j = 0; j < histograms_->atomTypeMix().nItems(); ++j) + histograms_->fullHistogram(k, j) = histograms[{k, j}]; // Atoms within the same cell for (int n = 0; n < cellArray.nCells(); ++n) @@ -247,7 +238,7 @@ bool GRNode::calculateGRCells(double grRange) if (typeI != AtomType::Ignore && typeJ != AtomType::Ignore) { // No need to perform MIM since we're in the same cell - rawGR_->fullHistogram(i->localTypeIndex(), j->localTypeIndex()).bin((i->r() - j->r()).magnitude()); + histograms_->fullHistogram(i->localTypeIndex(), j->localTypeIndex()).bin((i->r() - j->r()).magnitude()); } }); } @@ -278,8 +269,12 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) * Make sure histograms are set up, and reset any existing data */ - rawGR_->setUpHistograms(grRange, binWidth_.asDouble()); - rawGR_->reset(); + if (!histograms_) + { + histograms_.emplace(); + histograms_->initialise(targetConfiguration_->atomTypePopulations(), grRange, binWidth_.asDouble()); + } + histograms_->zeroBins(); /* * Calculate full (intra+inter) partials @@ -304,8 +299,6 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) */ const auto *box = targetConfiguration_->box(); - const auto &cells = targetConfiguration_->cells(); - timer.start(); // Loop over molecules @@ -328,7 +321,7 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) if (typeJ == AtomType::Ignore) return; - rawGR_->boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); + histograms_->boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); }); } @@ -342,21 +335,21 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) */ timer.start(); - auto success = - for_each_pair_early(rawGR_->nAtomTypes(), - [&](auto typeI, auto typeJ) -> EarlyReturn - { - // Create unbound histogram from total and bound data - rawGR_->unboundHistogram(typeI, typeJ) = rawGR_->fullHistogram(typeI, typeJ); - rawGR_->unboundHistogram(typeI, typeJ).add(rawGR_->boundHistogram(typeI, typeJ), -1.0); - - return EarlyReturn::Continue; - }); + auto success = for_each_pair_early( + rawGR_->nAtomTypes(), + [&](auto typeI, auto typeJ) -> EarlyReturn + { + // Create unbound histogram from total and bound data + histograms_->unboundHistogram(typeI, typeJ) = histograms_->fullHistogram(typeI, typeJ); + histograms_->unboundHistogram(typeI, typeJ).add(histograms_->boundHistogram(typeI, typeJ), -1.0); + + return EarlyReturn::Continue; + }); if (success.has_value() && !success.value()) return false; // Transform histogram data into radial distribution functions - rawGR_->formPartials(box->volume()); + histograms_->formPartials(*rawGR_, box->volume()); // Sum total functions rawGR_->formTotals(true); diff --git a/src/nodes/gr/process.cpp b/src/nodes/gr/process.cpp index 95932f8830..70ab5e470e 100644 --- a/src/nodes/gr/process.cpp +++ b/src/nodes/gr/process.cpp @@ -71,7 +71,7 @@ NodeConstants::ProcessResult GRNode::process() if (!rawGR_) { rawGR_.emplace(realSpeciesPopulations); - rawGR_->setUp(targetConfiguration_->atomTypePopulations(), grRange, binWidth_.asDouble()); + rawGR_->initialise(targetConfiguration_->atomTypePopulations()); unweightedGR_->setEffectiveDensity(targetConfiguration_->atomicDensity().value_or(0.0)); } diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp index ad5f20b71f..e29c04cc57 100644 --- a/src/nodes/neutronSQ/process.cpp +++ b/src/nodes/neutronSQ/process.cpp @@ -149,7 +149,7 @@ NodeConstants::ProcessResult NeutronSQNode::process() if (!weightedSQ_) { weightedSQ_.emplace(realSpeciesPopulations); - weightedSQ_->setUpPartials(unweightedSQ_->atomTypeMix()); + weightedSQ_->initialise(unweightedSQ_->atomTypeMix()); } // Calculate weighted S(Q) @@ -175,7 +175,7 @@ NodeConstants::ProcessResult NeutronSQNode::process() if (!weightedGR_) { weightedGR_.emplace(realSpeciesPopulations); - weightedGR_->setUpPartials(unweightedGR_->atomTypeMix()); + weightedGR_->initialise(unweightedGR_->atomTypeMix()); } // Calculate weighted g(r) diff --git a/src/nodes/sq/process.cpp b/src/nodes/sq/process.cpp index e5fee28ed5..56f9aa74f5 100644 --- a/src/nodes/sq/process.cpp +++ b/src/nodes/sq/process.cpp @@ -65,7 +65,7 @@ NodeConstants::ProcessResult SQNode::process() if (!unweightedSQ_) { unweightedSQ_.emplace(realSpeciesPopulations); - unweightedSQ_->setUpPartials(unweightedGR_->atomTypeMix()); + unweightedSQ_->initialise(unweightedGR_->atomTypeMix()); } /*