From 841fc6ffd1088a68aea51ba5898d00fbfebf8f08 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 10:32:09 +0100 Subject: [PATCH 1/9] Species populations should really be a std::map. --- src/classes/configuration.h | 4 ++-- src/classes/configuration_contents.cpp | 30 ++++++++++---------------- 2 files changed, 13 insertions(+), 21 deletions(-) diff --git a/src/classes/configuration.h b/src/classes/configuration.h index 35a1391a36..17f582934e 100644 --- a/src/classes/configuration.h +++ b/src/classes/configuration.h @@ -80,7 +80,7 @@ class Configuration : public Serialisable */ private: // Species populations present in the Configuration - std::vector> speciesPopulations_; + std::map speciesPopulations_; // AtomType populations in the configuration AtomTypeMix atomTypePopulations_; // Contents version, incremented whenever Configuration content or Atom positions change @@ -100,7 +100,7 @@ class Configuration : public Serialisable // Adjust population of specified Species in the Configuration void adjustSpeciesPopulation(const Species *sp, int delta); // Return Species populations within the Configuration - const std::vector> &speciesPopulations() const; + const std::map &speciesPopulations() const; // Return population of specified species within the Configuration int speciesPopulation(const Species *sp) const; // Return if the specified Species is present in the Configuration diff --git a/src/classes/configuration_contents.cpp b/src/classes/configuration_contents.cpp index b507d85ec6..77bd79cdcd 100644 --- a/src/classes/configuration_contents.cpp +++ b/src/classes/configuration_contents.cpp @@ -31,35 +31,27 @@ const AtomTypeMix &Configuration::atomTypePopulations() const { return atomTypeP // Adjust population of specified Species in the Configuration void Configuration::adjustSpeciesPopulation(const Species *sp, int delta) { - auto it = std::find_if(speciesPopulations_.begin(), speciesPopulations_.end(), - [sp](const auto &data) { return data.first == sp; }); - if (it == speciesPopulations_.end()) + if (speciesPopulations_.contains(sp)) + speciesPopulations_[sp] += delta; + else { if (delta < 0) Messenger::exception("Can't decrease population of Species '{}' as it is not in the list.\n", sp->name()); - speciesPopulations_.emplace_back(sp, delta); + speciesPopulations_[sp] = delta; } - else - it->second += delta; } // Return Species populations within the Configuration -const std::vector> &Configuration::speciesPopulations() const { return speciesPopulations_; } +const std::map &Configuration::speciesPopulations() const { return speciesPopulations_; } // Return population of specified species within the Configuration int Configuration::speciesPopulation(const Species *sp) const { - auto it = std::find_if(speciesPopulations_.begin(), speciesPopulations_.end(), - [sp](const auto &spInfo) { return spInfo.first == sp; }); - return it == speciesPopulations_.end() ? 0 : it->second; + return speciesPopulations_.contains(sp) ? speciesPopulations_.at(sp) : 0; } // Return if the specified Species is present in the Configuration -bool Configuration::containsSpecies(const Species *sp) -{ - return std::find_if(speciesPopulations_.begin(), speciesPopulations_.end(), - [sp](const auto &data) { return data.first == sp; }) != speciesPopulations_.end(); -} +bool Configuration::containsSpecies(const Species *sp) { return speciesPopulations_.contains(sp); } // Return the total charge of the Configuration double Configuration::totalCharge(bool ppIncludeCoulomb) const @@ -72,11 +64,11 @@ double Configuration::totalCharge(bool ppIncludeCoulomb) const // Return the total atomic mass present in the Configuration double Configuration::atomicMass() const { - double mass = 0.0; + auto mass = 0.0; // Get total molar mass in configuration - for (const auto &spPop : speciesPopulations_) - mass += spPop.first->mass() * spPop.second; + for (const auto &[sp, population] : speciesPopulations_) + mass += sp->mass() * population; // Convert to absolute mass return mass / DissolveMath::Avogadro; @@ -324,4 +316,4 @@ bool Configuration::energyIsStable() const { return energyIsStable_; } // Energy gradient void Configuration::setEnergyGradient(double grad) { energyGradient_ = grad; } -double Configuration::getEnergyGradient() const { return energyGradient_; } \ No newline at end of file +double Configuration::getEnergyGradient() const { return energyGradient_; } From 5914cd7fb04e886d2f0429b0eb6e7991b2edfd11 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 11:02:58 +0100 Subject: [PATCH 2/9] Clarify that PartialSet contains real (floating point) species populations by design. --- src/classes/partialSet.cpp | 7 +++++-- src/classes/partialSet.h | 10 ++++------ src/nodes/gr/gr.h | 2 +- src/nodes/gr/helpers.cpp | 20 +++++++++----------- src/nodes/gr/process.cpp | 7 ++++++- src/nodes/neutronSQ/helpers.cpp | 2 +- src/nodes/neutronSQ/process.cpp | 15 ++++++++++----- src/nodes/sq/process.cpp | 6 +++++- 8 files changed, 41 insertions(+), 28 deletions(-) diff --git a/src/classes/partialSet.cpp b/src/classes/partialSet.cpp index 53b45a617d..14ddf6f6bb 100644 --- a/src/classes/partialSet.cpp +++ b/src/classes/partialSet.cpp @@ -12,7 +12,10 @@ #include "math/mathFunc.h" #include "templates/algorithms.h" -PartialSet::PartialSet(const SpeciesPopulations &speciesPopulations) : speciesPopulations_(speciesPopulations) {} +PartialSet::PartialSet(const std::map &realSpeciesPopulations) + : realSpeciesPopulations_(realSpeciesPopulations) +{ +} PartialSet::~PartialSet() { @@ -268,7 +271,7 @@ Data1D &PartialSet::unboundTotal() { return unboundTotal_; } const Data1D &PartialSet::unboundTotal() const { return unboundTotal_; } // Species populations -SpeciesPopulations &PartialSet::speciesPopulations() { return speciesPopulations_; } +const std::map &PartialSet::realSpeciesPopulations() const { return realSpeciesPopulations_; } // Save all partials and total bool PartialSet::save(std::string_view prefix, std::string_view tag, std::string_view suffix, diff --git a/src/classes/partialSet.h b/src/classes/partialSet.h index 2705a132ef..99a8cfce60 100644 --- a/src/classes/partialSet.h +++ b/src/classes/partialSet.h @@ -13,14 +13,12 @@ class Configuration; class Interpolator; -using SpeciesPopulations = std::vector>; - // Set of Partials class PartialSet { public: PartialSet() = default; - PartialSet(const SpeciesPopulations &speciesPopulations); + PartialSet(const std::map &realSpeciesPopulations); ~PartialSet(); /* @@ -52,7 +50,7 @@ class PartialSet // Effective density double rho_; // Species populations - std::vector> speciesPopulations_; + std::map realSpeciesPopulations_; public: // Set up PartialSet, including initialising histograms for g(r) use @@ -104,8 +102,8 @@ class PartialSet // Return total unbound function Data1D &unboundTotal(); const Data1D &unboundTotal() const; - // Species populations - SpeciesPopulations &speciesPopulations(); + // Return real species populations + 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 diff --git a/src/nodes/gr/gr.h b/src/nodes/gr/gr.h index 38ad5faf69..3d4126e445 100644 --- a/src/nodes/gr/gr.h +++ b/src/nodes/gr/gr.h @@ -95,7 +95,7 @@ class GRNode : public Node // Calculate and return effective density based on target Configurations double effectiveDensity() const; // Calculate and return used species populations based on target Configurations - SpeciesPopulations speciesPopulations() const; + std::map realSpeciesPopulations() const; // (Re)calculate partial g(r) for the specified Configuration bool calculateGR(Configuration *cfg, PartialSet &originalgr, PartialsMethod method, const double rdfRange, const double rdfBinWidth, bool &alreadyUpToDate); diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp index 05c17bcd38..a3579f8141 100644 --- a/src/nodes/gr/helpers.cpp +++ b/src/nodes/gr/helpers.cpp @@ -261,7 +261,7 @@ bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const PartialSet &GRNode::originalGR(Configuration *cfg, const double rdfRange, const double rdfBinWidth) { if (!originalgr_) - originalgr_.emplace(speciesPopulations()); + originalgr_.emplace(realSpeciesPopulations()); originalgr_.value().setUp(cfg->atomTypePopulations(), rdfRange, rdfBinWidth); return originalgr_.value(); @@ -271,7 +271,7 @@ PartialSet &GRNode::originalGR(Configuration *cfg, const double rdfRange, const PartialSet &GRNode::unweightedGR() { if (!unweightedGR_) - unweightedGR_.emplace(speciesPopulations()); + unweightedGR_.emplace(realSpeciesPopulations()); return unweightedGR_.value(); } @@ -280,7 +280,7 @@ PartialSet &GRNode::unweightedGR() PartialSet &GRNode::summedUnweightedGR() { if (!summedUnweightedGR_) - summedUnweightedGR_.emplace(speciesPopulations()); + summedUnweightedGR_.emplace(realSpeciesPopulations()); return summedUnweightedGR_.value(); } @@ -312,23 +312,21 @@ double GRNode::effectiveDensity() const } // Calculate and return used species populations based on target Configurations -SpeciesPopulations GRNode::speciesPopulations() const +std::map GRNode::realSpeciesPopulations() const { - std::vector> populations; + std::map populations; for (auto *cfg : targetConfigurations_) { // TODO Get weight for configuration auto weight = 1.0; - for (const auto &spPop : cfg->speciesPopulations()) + for (const auto &[sp, population] : cfg->speciesPopulations()) { - auto it = std::find_if(populations.begin(), populations.end(), - [&spPop](auto &data) { return data.first == spPop.first; }); - if (it != populations.end()) - it->second += spPop.second * weight; + if (populations.contains(sp)) + populations[sp] += population * weight; else - populations.emplace_back(spPop.first, spPop.second * weight); + populations[sp] = population * weight; } } diff --git a/src/nodes/gr/process.cpp b/src/nodes/gr/process.cpp index 492a1c8584..b5eb0fbcb8 100644 --- a/src/nodes/gr/process.cpp +++ b/src/nodes/gr/process.cpp @@ -112,7 +112,12 @@ NodeConstants::ProcessResult GRNode::process() return NodeConstants::ProcessResult::Failed; unweightedGR().setEffectiveDensity(effectiveDensity()); - unweightedGR().speciesPopulations() = speciesPopulations(); + + // Set the real species populations + std::map realSpeciesPopulations; + // TODO + // std::transform(targetConfigurations_.) + // unweightedGR().realSpeciesPopulations() = speciesPopulations(); return NodeConstants::ProcessResult::Success; } diff --git a/src/nodes/neutronSQ/helpers.cpp b/src/nodes/neutronSQ/helpers.cpp index ef468a7242..f67c9751f3 100644 --- a/src/nodes/neutronSQ/helpers.cpp +++ b/src/nodes/neutronSQ/helpers.cpp @@ -101,7 +101,7 @@ void NeutronSQNode::calculateWeights(NeutronWeights &weights) const // Clear weights and get species speciesPopulations_ from GRModule weights.clear(); - for (auto &[sp, pop] : unweightedGR_->speciesPopulations()) + for (auto &[sp, pop] : unweightedGR_->realSpeciesPopulations()) { // Find the defined Isotopologue for this Species - if it doesn't exist, use the Natural one auto isoRef = isotopologueSet_.getIsotopologues(sp); diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp index 3c6a94b095..a3c9bc0c14 100644 --- a/src/nodes/neutronSQ/process.cpp +++ b/src/nodes/neutronSQ/process.cpp @@ -141,6 +141,9 @@ NodeConstants::ProcessResult NeutronSQNode::process() message("Isotopologue and isotope composition:\n\n"); weights_.print(); + // Get the real species populations from the input unweightedSQ + auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations(); + // Does a PartialSet for the weighted S(Q) already exist for this Configuration? /* auto [weightedSQ, wSQstatus] = dissolve.processingModuleData().realiseIf( @@ -149,14 +152,15 @@ NodeConstants::ProcessResult NeutronSQNode::process() weightedSQ.setUpPartials(unweightedSQ.atomTypeMix()); */ - if (!weightedGR_) + // Set up the weighted SQ storage if needed + if (!weightedSQ_) { - weightedSQ_.emplace(unweightedSQ_->speciesPopulations()); + weightedSQ_.emplace(realSpeciesPopulations); weightedSQ_->setUpPartials(unweightedSQ_->atomTypeMix()); } - auto &population = weightedGR_->speciesPopulations(); - for (const auto &[species, _] : population) + // Update the isotopologue set + for (const auto &[species, _] : realSpeciesPopulations) { for (const auto &isotopologue : species->isotopologues()) { @@ -202,9 +206,10 @@ NodeConstants::ProcessResult NeutronSQNode::process() weightedGR.setUpPartials(unweightedGR.atomTypeMix()); */ + // Set up weighted GR storage if we need it if (!weightedGR_) { - weightedGR_.emplace(unweightedGR_->speciesPopulations()); + weightedGR_.emplace(realSpeciesPopulations); weightedGR_->setUpPartials(unweightedGR_->atomTypeMix()); } diff --git a/src/nodes/sq/process.cpp b/src/nodes/sq/process.cpp index b1ce311e7a..e8e1a5226d 100644 --- a/src/nodes/sq/process.cpp +++ b/src/nodes/sq/process.cpp @@ -54,13 +54,17 @@ NodeConstants::ProcessResult SQNode::process() message("SQ: Save data is {}.\n", DissolveSys::onOff(save_)); message("\n"); + // Get the real species populations from the input unweightedSQ + auto &realSpeciesPopulations = unweightedGR_->realSpeciesPopulations(); + /* * Transform target UnweightedGR into the UnweightedSQ. */ + // Set up unweighted SQ storage if we need to if (!unweightedSQ_) { - unweightedSQ_.emplace(unweightedGR_->speciesPopulations()); + unweightedSQ_.emplace(realSpeciesPopulations); unweightedSQ_->setUpPartials(unweightedGR_->atomTypeMix()); } From e6265accab49b420cad8a172efa5e5dc33ae2f32 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 11:36:52 +0100 Subject: [PATCH 3/9] Simplify GR node making it act on a single configuration. --- src/modules/gr/process.cpp | 2 +- src/nodes/gr/gr.cpp | 6 +- src/nodes/gr/gr.h | 39 ++--- src/nodes/gr/helpers.cpp | 296 ++++++++++--------------------------- src/nodes/gr/process.cpp | 123 ++++++++------- 5 files changed, 150 insertions(+), 316 deletions(-) diff --git a/src/modules/gr/process.cpp b/src/modules/gr/process.cpp index fc49fbb29b..77db04e021 100644 --- a/src/modules/gr/process.cpp +++ b/src/modules/gr/process.cpp @@ -35,7 +35,7 @@ Module::ExecutionResult GRModule::process(Dissolve &dissolve) /* * Regardless of whether we are a main processing task (summing some combination of Configuration's partials) or - * multiple independent Configurations, we must loop over the specified targetConfigurations_ and calculate the partials + * multiple independent Configurations, we must loop over the specified targetConfiguration_ and calculate the partials * for each. */ diff --git a/src/nodes/gr/gr.cpp b/src/nodes/gr/gr.cpp index 210bb99b6e..d7fbb8aaa4 100644 --- a/src/nodes/gr/gr.cpp +++ b/src/nodes/gr/gr.cpp @@ -5,8 +5,7 @@ GRNode::GRNode(Graph *parentGraph) : Node(parentGraph) { - addInput>("Configurations", "Set target configuration(s) for the module", - targetConfigurations_) + addInput("Configurations", "Set target configuration(s) for the module", targetConfiguration_) ->setFlags({ParameterBase::Required, ParameterBase::ClearData}); addOption("BinWidth", "Bin width (spacing in r) to use", binWidth_); addOption>("Range", "Maximum r to calculate g(r) out to", requestedRange_); @@ -17,8 +16,7 @@ GRNode::GRNode(Graph *parentGraph) : Node(parentGraph) addOption("IntraBroadening", "Type of broadening to apply to intramolecular g(r)", intraBroadening_); addOption>("Smoothing", "Specifies the degree of smoothing to apply to calculated g(r)", nSmooths_); addOption("Save", "Whether to save partials and total functions to disk", save_); - addOption("SaveOriginal", "Whether to save original (unbroadened) partials and total functions to disk", - saveOriginal_); + addOption("SaveRaw", "Whether to save raw simulation partial and total functions to disk", saveRaw_); addOption( "InternalTest", "Perform internal check of calculated partials against a set calculated by a simple unoptimised double-loop", diff --git a/src/nodes/gr/gr.h b/src/nodes/gr/gr.h index 3d4126e445..d9e3b3ca38 100644 --- a/src/nodes/gr/gr.h +++ b/src/nodes/gr/gr.h @@ -44,13 +44,11 @@ class GRNode : public Node private: // Target configurations - std::vector targetConfigurations_; - // Original g(r) - std::optional originalgr_; + Configuration *targetConfiguration_{nullptr}; + // Raw simulation g(r) + std::optional rawGR_; // Unweighted g(r) std::optional unweightedGR_; - // Summed unweighted g(r) - std::optional summedUnweightedGR_; // Number of historical partial sets to combine into final partials std::optional averagingLength_{5}; // Weighting scheme to use when averaging partials @@ -71,40 +69,25 @@ class GRNode : public Node std::optional requestedRange_; // Whether to save partials and total functions to disk bool save_{false}; - // Whether to save original (unbroadened) partials and total functions to disk - bool saveOriginal_{false}; + // Whether to save raw partials and total functions to disk + bool saveRaw_{false}; /* * Functions */ private: // Calculate partial g(r) in serial with simple double-loop - bool calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet); + bool calculateGRTestSerial(); // Calculate partial g(r) with optimised double-loop - bool calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const double rdfRange); + bool calculateGRSimple(); // Calculate partial g(r) utilising Cell neighbour lists - bool calculateGRCells(Configuration *cfg, PartialSet &partialSet, const double binWidth); + bool calculateGRCells(double grRange); public: - // Get original g(r), constructing if empty - PartialSet &originalGR(Configuration *cfg, const double rdfRange, const double rdfBinWidth); - // Get unweighted g(r), constructing if empty - PartialSet &unweightedGR(); - // Get summed unweighted g(r), constructing if empty - PartialSet &summedUnweightedGR(); - // Calculate and return effective density based on target Configurations - double effectiveDensity() const; - // Calculate and return used species populations based on target Configurations - std::map realSpeciesPopulations() const; - // (Re)calculate partial g(r) for the specified Configuration - bool calculateGR(Configuration *cfg, PartialSet &originalgr, PartialsMethod method, const double rdfRange, - const double rdfBinWidth, bool &alreadyUpToDate); + // Calculate raw partials + bool calculateRawGR(const double grRange, bool &alreadyUpToDate); // Calculate smoothed/broadened partial g(r) from supplied partials - bool calculateUnweightedGR(Configuration *cfg, const PartialSet &originalgr, PartialSet &weightedgr, - const Function1DWrapper intraBroadening, int smoothing); - // Sum unweighted g(r) over the supplied Module's target Configurations - bool sumUnweightedGR(std::string_view targetPrefix, std::string_view parentPrefix, - const std::vector &parentCfgs, PartialSet &summedUnweightedGR); + bool calculateUnweightedGR(); // Test supplied PartialSets against each other bool testReferencePartials(PartialSet &setA, PartialSet &setB, double testThreshold); // Test calculated partial against supplied reference data diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp index a3579f8141..4d68be44f4 100644 --- a/src/nodes/gr/helpers.cpp +++ b/src/nodes/gr/helpers.cpp @@ -37,39 +37,39 @@ void addHistogramsToPartialSet(Array2D &histograms, PartialSet &tar */ // Calculate partial g(r) in serial with simple double-loop -bool GRNode::calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet) +bool GRNode::calculateGRTestSerial() { // Calculate radial distribution functions with a simple double loop, in serial - const auto *box = cfg->box(); + const auto *box = targetConfiguration_->box(); dissolve::for_each_pair( - ParallelPolicies::seq, cfg->atoms(), - [box, &partialSet](auto i, auto &ii, auto j, auto &jj) + ParallelPolicies::seq, targetConfiguration_->atoms(), + [&, 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())); + rawGR_->fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r())); }); return true; } // Calculate partial g(r) with optimised double-loop -bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const double binWidth) +bool GRNode::calculateGRSimple() { // Variables int n, m, nTypes, typeI, typeJ, i, j, nPoints; // Construct local arrays of atom type positions - nTypes = partialSet.nAtomTypes(); + nTypes = rawGR_->nAtomTypes(); message("Constructing local partial working arrays for {} types.\n", nTypes); - const auto *box = cfg->box(); + const auto *box = targetConfiguration_->box(); std::vector r(nTypes); std::vector maxr(nTypes), nr(nTypes); std::vector binss(nTypes); int *bins; n = 0; - for (auto &atd : cfg->atomTypePopulations()) + for (auto &atd : targetConfiguration_->atomTypePopulations()) { maxr[n] = atd.population(); nr[n] = 0; @@ -79,7 +79,7 @@ bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const } // Loop over Atoms and construct arrays - for (auto &atom : cfg->atoms()) + for (auto &atom : targetConfiguration_->atoms()) { m = atom.localTypeIndex(); if (m == AtomType::Ignore) @@ -91,7 +91,7 @@ bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const // Loop over assigned Atoms Vector3 centre, *ri, *rj, mim; - double rbin = 1.0 / binWidth; + double rbin = 1.0 / binWidth_.asDouble(); message("Self terms..\n"); @@ -99,9 +99,9 @@ bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const for (typeI = 0; typeI < nTypes; ++typeI) { ri = r[typeI]; - auto &histogram = partialSet.fullHistogram(typeI, typeI).bins(); + auto &histogram = rawGR_->fullHistogram(typeI, typeI).bins(); bins = binss[typeI]; - nPoints = partialSet.fullHistogram(typeI, typeI).nBins(); + nPoints = rawGR_->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 +135,9 @@ bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const continue; rj = r[typeJ]; - auto &histogram = partialSet.fullHistogram(typeI, typeJ).bins(); + auto &histogram = rawGR_->fullHistogram(typeI, typeJ).bins(); bins = binss[typeJ]; - nPoints = partialSet.fullHistogram(typeI, typeJ).nBins(); + nPoints = rawGR_->fullHistogram(typeI, typeJ).nBins(); for (i = 0; i < maxr[typeI]; ++i) { centre = ri[i]; @@ -160,35 +160,36 @@ bool GRNode::calculateGRSimple(Configuration *cfg, PartialSet &partialSet, const return true; } -bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const double rdfRange) +// Calculate partial g(r) utilising Cell neighbour lists +bool GRNode::calculateGRCells(double grRange) { - auto &cellArray = cfg->cells(); + auto &cellArray = targetConfiguration_->cells(); // Loop context is to use all processes in Pool as one group 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(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); return histograms; }); - auto unaryOp = [&combinableHistograms, cfg, &comb, rdfRange](const auto idx) + auto unaryOp = [&, grRange](const auto idx) { // auto &histograms = combinableHistograms.local().histograms_; auto &histograms = combinableHistograms.local(); - const auto *box = cfg->box(); - auto &cellArray = cfg->cells(); + const auto *box = targetConfiguration_->box(); + auto &cellArray = targetConfiguration_->cells(); auto [n, m] = comb.nthCombination(idx); auto *cellI = cellArray.cell(n); auto *cellJ = cellArray.cell(m); - if (!cellArray.withinMinimumImageRange(cellI, cellJ, rdfRange)) + if (!cellArray.withinMinimumImageRange(cellI, cellJ, grRange)) return; // Add contributions between atoms in cellI and cellJ @@ -222,7 +223,7 @@ bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const dissolve::for_each(ParallelPolicies::par, dissolve::counting_iterator(0), dissolve::counting_iterator(comb.getNumCombinations()), unaryOp); auto histograms = combinableHistograms.finalize(); - addHistogramsToPartialSet(histograms, partialSet); + addHistogramsToPartialSet(histograms, *rawGR_); // Atoms within the same cell for (int n = 0; n < cellArray.nCells(); ++n) @@ -234,7 +235,7 @@ bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const PairIterator pairs(atomsI.size()); std::for_each( pairs.begin(), pairs.end(), - [&atomsI, &partialSet](auto it) + [&, atomsI](auto it) { auto [idx, jdx] = it; if (idx == jdx) @@ -246,7 +247,7 @@ bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const 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()); + rawGR_->fullHistogram(i->localTypeIndex(), j->localTypeIndex()).bin((i->r() - j->r()).magnitude()); } }); } @@ -257,120 +258,43 @@ bool GRNode::calculateGRCells(Configuration *cfg, PartialSet &partialSet, const * Public Functions */ -// Get original g(r), constructing if empty -PartialSet &GRNode::originalGR(Configuration *cfg, const double rdfRange, const double rdfBinWidth) -{ - if (!originalgr_) - originalgr_.emplace(realSpeciesPopulations()); - originalgr_.value().setUp(cfg->atomTypePopulations(), rdfRange, rdfBinWidth); - - return originalgr_.value(); -} - -// Get original g(r), constructing if empty -PartialSet &GRNode::unweightedGR() -{ - if (!unweightedGR_) - unweightedGR_.emplace(realSpeciesPopulations()); - - return unweightedGR_.value(); -} - -// Get summed unweighted g(r), constructing if empty -PartialSet &GRNode::summedUnweightedGR() -{ - if (!summedUnweightedGR_) - summedUnweightedGR_.emplace(realSpeciesPopulations()); - - return summedUnweightedGR_.value(); -} - -// Calculate and return effective density based on target Configurations -double GRNode::effectiveDensity() const -{ - double rho0 = 0; - auto totalWeight = 0.0; - for (auto *cfg : targetConfigurations_) - { - auto cfgRho = cfg->atomicDensity(); - if (!cfgRho) - continue; - - // TODO Get weight for configuration - auto weight = 1.0; - - totalWeight += weight; - - // Add to sum - if (rho0) - rho0 += weight / *cfg->atomicDensity(); - else - rho0 = weight / *cfg->atomicDensity(); - } - - return 1.0 / (rho0 / totalWeight); -} - -// Calculate and return used species populations based on target Configurations -std::map GRNode::realSpeciesPopulations() const -{ - std::map populations; - - for (auto *cfg : targetConfigurations_) - { - // TODO Get weight for configuration - auto weight = 1.0; - - for (const auto &[sp, population] : cfg->speciesPopulations()) - { - if (populations.contains(sp)) - populations[sp] += population * weight; - else - populations[sp] = population * weight; - } - } - - return populations; -} - -// Calculate unweighted partials for the specified Configuration -bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::PartialsMethod method, const double rdfRange, - const double rdfBinWidth, bool &alreadyUpToDate) +// Calculate raw partials +bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) { // 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 alreadyUpToDate = false; - if (DissolveSys::sameString(originalgr_.value().fingerprint(), std::format("{}", cfg->contentsVersion())) && - (method != PartialsMethod::TestMethod)) + if (DissolveSys::sameString(rawGR_->fingerprint(), std::format("{}", targetConfiguration_->contentsVersion())) && + (partialsMethod_ != PartialsMethod::TestMethod)) { - message("Partial g(r) are up-to-date for Configuration '{}'.\n", cfg->name()); + message("Partial g(r) are up-to-date for Configuration '{}'.\n", targetConfiguration_->name()); alreadyUpToDate = true; return true; } - message("Calculating partial g(r) for Configuration '{}'...\n", cfg->name()); + message("Calculating partial g(r) for Configuration '{}'...\n", targetConfiguration_->name()); /* * Make sure histograms are set up, and reset any existing data */ - originalgr.setUpHistograms(rdfRange, rdfBinWidth); - originalgr.reset(); + rawGR_->setUpHistograms(grRange, binWidth_.asDouble()); + rawGR_->reset(); /* * Calculate full (intra+inter) partials */ Timer timer; - if (method == PartialsMethod::TestMethod) - calculateGRTestSerial(cfg, originalgr_.value()); - else if (method == PartialsMethod::SimpleMethod) - calculateGRSimple(cfg, originalgr, rdfBinWidth); - else if (method == PartialsMethod::CellsMethod) - calculateGRCells(cfg, originalgr, rdfRange); - else if (method == PartialsMethod::AutoMethod) + if (partialsMethod_ == PartialsMethod::TestMethod) + calculateGRTestSerial(); + else if (partialsMethod_ == PartialsMethod::SimpleMethod) + calculateGRSimple(); + else if (partialsMethod_ == PartialsMethod::CellsMethod) + calculateGRCells(grRange); + else if (partialsMethod_ == PartialsMethod::AutoMethod) { - cfg->nAtoms() > 10000 ? calculateGRCells(cfg, originalgr, rdfRange) : calculateGRSimple(cfg, originalgr, rdfBinWidth); + targetConfiguration_->nAtoms() > 10000 ? calculateGRCells(grRange) : calculateGRSimple(); } timer.stop(); message("Finished calculation of partials ({} elapsed).\n", timer.totalTimeString()); @@ -379,18 +303,18 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par * Calculate intramolecular partials */ - const auto *box = cfg->box(); - const auto &cells = cfg->cells(); + const auto *box = targetConfiguration_->box(); + const auto &cells = targetConfiguration_->cells(); timer.start(); // Loop over molecules - for (auto &mol : cfg->molecules()) + for (auto &mol : targetConfiguration_->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) @@ -404,7 +328,7 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par if (typeJ == AtomType::Ignore) return; - originalgr.boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); + rawGR_->boundHistogram(typeI, typeJ).bin(box->minimumDistance(i->r(), j->r())); }); } @@ -419,12 +343,12 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par timer.start(); auto success = - for_each_pair_early(originalgr.nAtomTypes(), - [&originalgr](auto typeI, auto typeJ) -> EarlyReturn + for_each_pair_early(rawGR_->nAtomTypes(), + [&](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); + rawGR_->unboundHistogram(typeI, typeJ) = rawGR_->fullHistogram(typeI, typeJ); + rawGR_->unboundHistogram(typeI, typeJ).add(rawGR_->boundHistogram(typeI, typeJ), -1.0); return EarlyReturn::Continue; }); @@ -432,10 +356,10 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par return false; // Transform histogram data into radial distribution functions - originalgr.formPartials(box->volume()); + rawGR_->formPartials(box->volume()); // Sum total functions - originalgr.formTotals(true); + rawGR_->formTotals(true); timer.stop(); message("Finished summation and normalisation of partial g(r) data ({}).\n", timer.totalTimeString()); @@ -447,123 +371,59 @@ bool GRNode::calculateGR(Configuration *cfg, PartialSet &originalgr, GRNode::Par } // Calculate smoothed/broadened partial g(r) from supplied partials -bool GRNode::calculateUnweightedGR(Configuration *cfg, const PartialSet &originalgr, PartialSet &unweightedgr, - const Function1DWrapper intraBroadening, int smoothing) +bool GRNode::calculateUnweightedGR() { - // If the unweightedgr is not yet initialised, copy the originalgr. Otherwise, just copy the values (in order to + // If the unweightedGR_ is not yet initialised, copy the rawGR_-> Otherwise, just copy the values (in order to // maintain the incremental versioning of the data) - if (unweightedgr.nAtomTypes() == 0) - unweightedgr = originalgr; + if (unweightedGR_->nAtomTypes() == 0) + *unweightedGR_ = *rawGR_; else { - for (auto i = 0; i < unweightedgr.nAtomTypes(); ++i) + for (auto i = 0; i < unweightedGR_->nAtomTypes(); ++i) { - for (auto j = i; j < unweightedgr.nAtomTypes(); ++j) + for (auto j = i; j < unweightedGR_->nAtomTypes(); ++j) { - unweightedgr.boundPartial(i, j).copyArrays(originalgr.boundPartial(i, j)); - unweightedgr.unboundPartial(i, j).copyArrays(originalgr.unboundPartial(i, j)); - unweightedgr.partial(i, j).copyArrays(originalgr.partial(i, j)); + unweightedGR_->boundPartial(i, j).copyArrays(rawGR_->boundPartial(i, j)); + unweightedGR_->unboundPartial(i, j).copyArrays(rawGR_->unboundPartial(i, j)); + unweightedGR_->partial(i, j).copyArrays(rawGR_->partial(i, j)); } } - unweightedgr.total().copyArrays(originalgr.total()); + unweightedGR_->total().copyArrays(rawGR_->total()); } // Remove bound partial from full partial - for (auto i = 0; i < unweightedgr.nAtomTypes(); ++i) + for (auto i = 0; i < unweightedGR_->nAtomTypes(); ++i) { - for (auto j = i; j < unweightedgr.nAtomTypes(); ++j) - unweightedgr.partial(i, j) -= originalgr.boundPartial(i, j); + for (auto j = i; j < unweightedGR_->nAtomTypes(); ++j) + unweightedGR_->partial(i, j) -= rawGR_->boundPartial(i, j); } // Broaden the bound partials according to the supplied PairBroadeningFunction - auto &types = unweightedgr.atomTypeMix(); + auto &types = unweightedGR_->atomTypeMix(); dissolve::for_each_pair(ParallelPolicies::seq, types, [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ) - { Filters::convolve(unweightedgr.boundPartial(i, j), intraBroadening, true, true); }); + { Filters::convolve(unweightedGR_->boundPartial(i, j), intraBroadening_, true, true); }); // Add broadened bound partials back in to full partials dissolve::for_each_pair(ParallelPolicies::seq, types, [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ) - { unweightedgr.partial(i, j) += unweightedgr.boundPartial(i, j); }); + { unweightedGR_->partial(i, j) += unweightedGR_->boundPartial(i, j); }); // Apply smoothing if requested + auto smoothing = nSmooths_.value_or(0).asInteger(); if (smoothing > 0) { dissolve::for_each_pair(ParallelPolicies::seq, types, [&](int i, const AtomTypeData &typeI, int j, const AtomTypeData &typeJ) { - Filters::movingAverage(unweightedgr.partial(i, j), smoothing); - Filters::movingAverage(unweightedgr.boundPartial(i, j), smoothing); - Filters::movingAverage(unweightedgr.unboundPartial(i, j), smoothing); + Filters::movingAverage(unweightedGR_->partial(i, j), smoothing); + Filters::movingAverage(unweightedGR_->boundPartial(i, j), smoothing); + Filters::movingAverage(unweightedGR_->unboundPartial(i, j), smoothing); }); } // Calculate total - unweightedgr.formTotals(true); - - return true; -} - -// Sum unweighted g(r) over the supplied Module's target Configurations -bool GRNode::sumUnweightedGR(std::string_view targetPrefix, std::string_view parentPrefix, - const std::vector &parentCfgs, PartialSet &summedUnweightedGR) -{ - combinedAtomTypes_.clear(); - for (Configuration *cfg : parentCfgs) - combinedAtomTypes_.add(cfg->atomTypePopulations()); - - // Finalise and save the combined AtomTypes matrix - combinedAtomTypes_.finalise(); - - // Set up PartialSet container - summedUnweightedGR.setUpPartials(combinedAtomTypes_); - - // Determine total weighting factors and combined density over all Configurations, and set up a Configuration/weight - // Vector for simplicity - std::vector> configWeights; - double totalWeight = 0.0; - for (Configuration *cfg : parentCfgs) - { - // Confirm atomic density is available (for the subsequent accumulator) - if (!cfg->atomicDensity()) - { - error("No density available for target configuration '{}'\n", cfg->name()); - return false; - } - - // TODO Assume weight of 1.0 - auto weight = 1.0; - - // Add our Configuration target - configWeights.emplace_back(cfg, weight); - totalWeight += weight; - } - - // Calculate overall density of combined system - double rho0 = std::accumulate(configWeights.begin(), configWeights.end(), 0.0, - [totalWeight](double acc, auto pair) - { return acc + pair.second / totalWeight / pair.first->atomicDensity().value(); }); - rho0 = 1.0 / rho0; - - // Sum Configurations into the PartialSet - std::string fingerprint; - for (auto [cfg, cfgWeight] : configWeights) - { - if (!cfg->atomicDensity()) - { - error("No density available for target configuration '{}'\n", cfg->name()); - return false; - } - - // Update fingerprint - fingerprint += - fingerprint.empty() ? std::format("{}", cfg->contentsVersion()) : std::format("_{}", cfg->contentsVersion()); - - // Calculate weighting factor - double weight = ((cfgWeight / totalWeight) * *cfg->atomicDensity()) / rho0; - - summedUnweightedGR.addPartials(unweightedGR(), weight); - } + unweightedGR_->formTotals(true); return true; } diff --git a/src/nodes/gr/process.cpp b/src/nodes/gr/process.cpp index b5eb0fbcb8..a15326c7e1 100644 --- a/src/nodes/gr/process.cpp +++ b/src/nodes/gr/process.cpp @@ -36,88 +36,81 @@ NodeConstants::ProcessResult GRNode::process() Functions1D::forms().keyword(intraBroadening_.form()), intraBroadening_.parameterSummary()); message("Calculation method is '{}'.\n", partialsMethods().keyword(partialsMethod_)); message("Save data is {}.\n", DissolveSys::onOff(save_)); - message("Save original (unbroadened) g(r) is {}.\n", DissolveSys::onOff(saveOriginal_)); + message("Save raw simulation g(r) is {}.\n", DissolveSys::onOff(saveRaw_)); if (nSmooths_) message("Degree of smoothing to apply to calculated partial g(r) is {}.\n", nSmooths_.value().asInteger()); message("\n"); - /* - * Regardless of whether we are a main processing task (summing some combination of Configuration's partials) or - * multiple independent Configurations, we must loop over the specified targetConfigurations_ and calculate the partials - * for each. - */ - - for (auto *cfg : targetConfigurations_) + // Check range + auto grRange = targetConfiguration_->box()->inscribedSphereRadius(); + if (!requestedRange_) + message("Maximal cutoff used for Configuration '{}' ({} Angstroms).\n", targetConfiguration_->niceName(), grRange); + else { - // Check RDF range - double rdfRange = cfg->box()->inscribedSphereRadius(); - if (!requestedRange_) - message("Maximal cutoff used for Configuration '{}' ({} Angstroms).\n", cfg->niceName(), rdfRange); - else + if (requestedRange_.value_or(Number(0.0)) > grRange) { - - if (requestedRange_.value_or(Number(0.0)) > rdfRange) - { - error("Specified RDF range of {} Angstroms is out of range for Configuration " - "'{}' (max = {} Angstroms).\n", - requestedRange_.value().asDouble(), cfg->niceName(), rdfRange); - return NodeConstants::ProcessResult::Failed; - } - - rdfRange = requestedRange_.value().asDouble(); - message("Cutoff for Configuration '{}' is {} Angstroms.\n", cfg->niceName(), rdfRange); + error("Specified RDF range of {} Angstroms is out of range for Configuration " + "'{}' (max = {} Angstroms).\n", + requestedRange_.value().asDouble(), targetConfiguration_->niceName(), grRange); + return NodeConstants::ProcessResult::Failed; } - // 'Snap' rdfRange_ to nearest bin width... - rdfRange = int(rdfRange / binWidth_.asDouble()) * binWidth_.asDouble(); - message("Cutoff (snapped to bin width) is {} Angstroms.\n", rdfRange); + grRange = requestedRange_.value().asDouble(); + message("Cutoff for Configuration '{}' is {} Angstroms.\n", targetConfiguration_->niceName(), grRange); + } - // Calculate unweighted partials for this Configuration - bool alreadyUpToDate; - calculateGR(cfg, originalGR(cfg, rdfRange, binWidth_.asDouble()), partialsMethod_, rdfRange, binWidth_.asDouble(), - alreadyUpToDate); + // 'Snap' grRange to nearest bin width... + grRange = int(grRange / binWidth_.asDouble()) * binWidth_.asDouble(); + message("Cutoff (snapped to bin width) is {} Angstroms.\n", grRange); - // Perform averagingLength_ of unweighted partials if requested, and if we're not already up-to-date - /* - if ((averagingLength_.value_or(1) > 1) && (!alreadyUpToDate)) - { - // Store the current fingerprint, since we must ensure we retain it in the averaged T. - std::string currentFingerprint{originalgr_.fingerprint()}; + // Convert configuration species populations into real species populations + std::map realSpeciesPopulations; + for (auto &[sp, iPop] : targetConfiguration_->speciesPopulations()) + realSpeciesPopulations[sp] = iPop; - Averaging::average(dissolve().processingModuleData(), std::format("{}//OriginalGR", cfg->niceName()), - name(), averagingLength_.value().asDouble(), averagingScheme_); - } - */ + // Create original GR storage if we need it + if (!rawGR_) + { + rawGR_.emplace(realSpeciesPopulations); + rawGR_->setUp(targetConfiguration_->atomTypePopulations(), grRange, binWidth_.asDouble()); + unweightedGR_->setEffectiveDensity(targetConfiguration_->atomicDensity().value_or(0.0)); + } - /* - // Perform internal test of original g(r)? - if (internalTest_) - { - // Copy the already-calculated g(r), then calculate a new set using the Test method - PartialSet referencePartials = originalgr; - calculateGR(dissolve.processingModuleData(), moduleContext.processPool(), cfg, GRModule::TestMethod, - rdfRange, binWidth_, alreadyUpToDate); - if (!testReferencePartials(referencePartials, originalgr, 1.0e-6)) - return ExecutionResult::Failed; - } - */ + // Calculate unweighted partials for this Configuration + bool alreadyUpToDate; + calculateRawGR(grRange, alreadyUpToDate); - // Form unweighted g(r) from original g(r), applying any requested nSmooths_.asInteger() / intramolecular broadening - calculateUnweightedGR(cfg, originalGR(cfg, rdfRange, binWidth_.asDouble()), unweightedGR(), intraBroadening_, - nSmooths_.value_or(0).asInteger()); + // Perform averagingLength_ of unweighted partials if requested, and if we're not already up-to-date + /* + if ((averagingLength_.value_or(1) > 1) && (!alreadyUpToDate)) + { + // Store the current fingerprint, since we must ensure we retain it in the averaged T. + std::string currentFingerprint{rawGR_.fingerprint()}; + + Averaging::average(dissolve().processingModuleData(), std::format("{}//OriginalGR", + targetConfiguration_->niceName()), name(), averagingLength_.value().asDouble(), averagingScheme_); } + */ - // Sum the partials from the associated Configurations - if (!sumUnweightedGR(name(), name(), targetConfigurations_, summedUnweightedGR())) - return NodeConstants::ProcessResult::Failed; + /* + // Perform internal test of original g(r)? + if (internalTest_) + { + // Copy the already-calculated g(r), then calculate a new set using the Test method + PartialSet referencePartials = originalgr; + calculateGR(dissolve.processingModuleData(), moduleContext.processPool(), cfg, GRModule::TestMethod, + grRange, binWidth_, alreadyUpToDate); + if (!testReferencePartials(referencePartials, originalgr, 1.0e-6)) + return ExecutionResult::Failed; + } + */ - unweightedGR().setEffectiveDensity(effectiveDensity()); + // Create unweighted GR storage if we need it + if (!unweightedGR_) + unweightedGR_.emplace(); - // Set the real species populations - std::map realSpeciesPopulations; - // TODO - // std::transform(targetConfigurations_.) - // unweightedGR().realSpeciesPopulations() = speciesPopulations(); + // Form unweighted g(r) from original g(r), applying any requested smoothing and/or intramolecular broadening + calculateUnweightedGR(); return NodeConstants::ProcessResult::Success; } From 060c93a9c01c331f27eba79f343a69f732c3d570 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 11:46:39 +0100 Subject: [PATCH 4/9] Tidy up SQ in a similar way. --- src/nodes/sq/helpers.cpp | 53 ++++++++++++++++++++++------------------ src/nodes/sq/process.cpp | 3 +-- src/nodes/sq/sq.h | 4 +-- 3 files changed, 32 insertions(+), 28 deletions(-) diff --git a/src/nodes/sq/helpers.cpp b/src/nodes/sq/helpers.cpp index 6f95527926..49c7333d9e 100644 --- a/src/nodes/sq/helpers.cpp +++ b/src/nodes/sq/helpers.cpp @@ -16,40 +16,45 @@ */ // Generate S(Q) from supplied g(r) -bool SQNode::calculateUnweightedSQ(const PartialSet &unweightedgr, PartialSet &unweightedsq, double qMin, double qDelta, - double qMax, double rho, const WindowFunction &windowFunction, Function1DWrapper broadening) +bool SQNode::calculateUnweightedSQ() { // Copy partial g(r) into our new S(Q) object - it should have been initialised already, so we will just check its size - if (unweightedgr.nAtomTypes() != unweightedsq.nAtomTypes()) + if (unweightedGR_->nAtomTypes() != unweightedSQ_->nAtomTypes()) + { error("SQNode::calculateUnweightedSQ - sizes of supplied partial sets are different.\n"); - return false; + return false; + } + + auto rho = unweightedGR_->effectiveDensity(); + auto qMin = qMin_.asDouble(), qDelta = qDelta_.asDouble(), qMax = qMax_.asDouble(); // Subtract 1.0 from the full and unbound partials so as to give (g(r)-1) and FT into S(Q) // Don't subtract 1.0 from the bound partials Timer timer; timer.start(); - dissolve::for_each_pair( - ParallelPolicies::par, unweightedgr.nAtomTypes(), - [&](int n, int m) - { - // Total partial - unweightedsq.partial(n, m).copyArrays(unweightedgr.partial(n, m)); - unweightedsq.partial(n, m) -= 1.0; - Fourier::sineFT(unweightedsq.partial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, windowFunction, broadening); - - // Bound partial - unweightedsq.boundPartial(n, m).copyArrays(unweightedgr.boundPartial(n, m)); - Fourier::sineFT(unweightedsq.boundPartial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, windowFunction, broadening); - - // Unbound partial - unweightedsq.unboundPartial(n, m).copyArrays(unweightedgr.unboundPartial(n, m)); - unweightedsq.unboundPartial(n, m) -= 1.0; - Fourier::sineFT(unweightedsq.unboundPartial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, windowFunction, - broadening); - }); + dissolve::for_each_pair(ParallelPolicies::par, unweightedGR_->nAtomTypes(), + [&](int n, int m) + { + // Total partial + unweightedSQ_->partial(n, m).copyArrays(unweightedGR_->partial(n, m)); + unweightedSQ_->partial(n, m) -= 1.0; + Fourier::sineFT(unweightedSQ_->partial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, + windowFunction_, qBroadening_); + + // Bound partial + unweightedSQ_->boundPartial(n, m).copyArrays(unweightedGR_->boundPartial(n, m)); + Fourier::sineFT(unweightedSQ_->boundPartial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, + windowFunction_, qBroadening_); + + // Unbound partial + unweightedSQ_->unboundPartial(n, m).copyArrays(unweightedGR_->unboundPartial(n, m)); + unweightedSQ_->unboundPartial(n, m) -= 1.0; + Fourier::sineFT(unweightedSQ_->unboundPartial(n, m), 4.0 * M_PI * rho, qMin, qDelta, qMax, + windowFunction_, qBroadening_); + }); // Sum into total - unweightedsq.formTotals(true); + unweightedSQ_->formTotals(true); timer.stop(); message("Finished Fourier transform and summation of partial g(r) into partial S(Q) ({} elapsed).\n", diff --git a/src/nodes/sq/process.cpp b/src/nodes/sq/process.cpp index e8e1a5226d..e5fee28ed5 100644 --- a/src/nodes/sq/process.cpp +++ b/src/nodes/sq/process.cpp @@ -78,8 +78,7 @@ NodeConstants::ProcessResult SQNode::process() */ // Transform g(r) into S(Q) - if (!calculateUnweightedSQ(*unweightedGR_, *unweightedSQ_, qMin, qDelta, qMax, unweightedGR_->effectiveDensity(), - WindowFunction(windowFunction_), qBroadening_)) + if (!calculateUnweightedSQ()) return NodeConstants::ProcessResult::Failed; /* diff --git a/src/nodes/sq/sq.h b/src/nodes/sq/sq.h index 55d20e2982..c95e6f883c 100644 --- a/src/nodes/sq/sq.h +++ b/src/nodes/sq/sq.h @@ -58,8 +58,8 @@ class SQNode : public Node */ public: // Calculate unweighted S(Q) from unweighted g(r) - bool calculateUnweightedSQ(const PartialSet &unweightedgr, PartialSet &unweightedsq, double qMin, double qDelta, - double qMax, double rho, const WindowFunction &windowFunction, Function1DWrapper broadening); + bool calculateUnweightedSQ(); + /* * Processing */ From 8c7529dcaec24bbac8bc30132c56ae8f532f1825 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 14:22:31 +0100 Subject: [PATCH 5/9] And now NeutronSQ. --- src/classes/neutronWeights.cpp | 24 +++++++ src/classes/neutronWeights.h | 6 ++ src/nodes/neutronSQ/helpers.cpp | 110 ++++++++++++++---------------- src/nodes/neutronSQ/neutronSQ.cpp | 6 +- src/nodes/neutronSQ/neutronSQ.h | 16 ++--- src/nodes/neutronSQ/process.cpp | 51 ++++---------- 6 files changed, 102 insertions(+), 111 deletions(-) diff --git a/src/classes/neutronWeights.cpp b/src/classes/neutronWeights.cpp index adb65c1886..f87c9ccea3 100644 --- a/src/classes/neutronWeights.cpp +++ b/src/classes/neutronWeights.cpp @@ -4,6 +4,7 @@ #include "classes/neutronWeights.h" #include "base/lineParser.h" #include "classes/atomType.h" +#include "classes/isotopologueSet.h" #include "classes/species.h" #include "data/isotopes.h" #include "items/deserialisers.h" @@ -215,6 +216,29 @@ void NeutronWeights::calculateWeightingMatrices() }); } +// Create from species populations and isotopologues +void NeutronWeights::create(const std::map &populations, const IsotopologueSet &isotopologues, + const std::vector> &exchangeableTypes) +{ + clear(); + + for (auto &[sp, pop] : populations) + { + // Find the defined Isotopologue for this Species - if it doesn't exist, use the Natural one + auto isoRef = isotopologues.getIsotopologues(sp); + if (isoRef) + { + const Isotopologues &topes = *isoRef; + for (const auto &isoWeight : topes.mix()) + addIsotopologue(sp, pop, isoWeight.isotopologue(), isoWeight.weight()); + } + else + addIsotopologue(sp, pop, sp->naturalIsotopologue(), 1.0); + } + + createFromIsotopologues(exchangeableTypes); +} + // Create AtomType list and matrices based on stored Isotopologues information void NeutronWeights::createFromIsotopologues(const std::vector> &exchangeableTypes) { diff --git a/src/classes/neutronWeights.h b/src/classes/neutronWeights.h index b915ee9cdf..02201e38b3 100644 --- a/src/classes/neutronWeights.h +++ b/src/classes/neutronWeights.h @@ -8,6 +8,9 @@ #include "templates/array2D.h" #include +// Forward Declarations +class IsotopologueSet; + // Neutron Weights Container class NeutronWeights { @@ -60,6 +63,9 @@ class NeutronWeights void calculateWeightingMatrices(); public: + // Create from species populations and isotopologues + void create(const std::map &populations, const IsotopologueSet &isotopologues, + const std::vector> &exchangeableTypes); // Create AtomType list and matrices based on stored Isotopologues information void createFromIsotopologues(const std::vector> &exchangeableTypes); // Reduce data to be naturally-weighted diff --git a/src/nodes/neutronSQ/helpers.cpp b/src/nodes/neutronSQ/helpers.cpp index f67c9751f3..8701b4bb7d 100644 --- a/src/nodes/neutronSQ/helpers.cpp +++ b/src/nodes/neutronSQ/helpers.cpp @@ -1,119 +1,113 @@ // SPDX-License-Identifier: GPL-3.0-or-later // Copyright (c) 2025 Team Dissolve and contributors -#include "classes/configuration.h" #include "classes/isotopologueSet.h" #include "classes/species.h" -#include "nodes/gr/gr.h" #include "nodes/neutronSQ/neutronSQ.h" -// Calculate weighted g(r) from supplied unweighted g(r) and neutron weights -bool NeutronSQNode::calculateWeightedGR(const PartialSet &unweightedgr, PartialSet &weightedgr, NeutronWeights &weights, - StructureFactors::NormalisationType normalisation) +// Calculate weighted g(r) +bool NeutronSQNode::calculateWeightedGR() { int typeI, typeJ; - for (typeI = 0; typeI < unweightedgr.nAtomTypes(); ++typeI) + for (typeI = 0; typeI < unweightedGR_->nAtomTypes(); ++typeI) { - for (typeJ = typeI; typeJ < unweightedgr.nAtomTypes(); ++typeJ) + for (typeJ = typeI; typeJ < unweightedGR_->nAtomTypes(); ++typeJ) { - double weight = weights.weight(typeI, typeJ); - double intraWeight = weights.intramolecularWeight(typeI, typeJ); + double weight = weights_.weight(typeI, typeJ); + double intraWeight = weights_.intramolecularWeight(typeI, typeJ); // Bound (intramolecular) partial (multiplied by the bound term weight) - weightedgr.boundPartial(typeI, typeJ).copyArrays(unweightedgr.boundPartial(typeI, typeJ)); - weightedgr.boundPartial(typeI, typeJ) *= intraWeight; + weightedGR_->boundPartial(typeI, typeJ).copyArrays(unweightedGR_->boundPartial(typeI, typeJ)); + weightedGR_->boundPartial(typeI, typeJ) *= intraWeight; // Unbound partial (multiplied by the full weight) - weightedgr.unboundPartial(typeI, typeJ).copyArrays(unweightedgr.unboundPartial(typeI, typeJ)); - weightedgr.unboundPartial(typeI, typeJ) -= 1.0; - weightedgr.unboundPartial(typeI, typeJ) *= weight; + weightedGR_->unboundPartial(typeI, typeJ).copyArrays(unweightedGR_->unboundPartial(typeI, typeJ)); + weightedGR_->unboundPartial(typeI, typeJ) -= 1.0; + weightedGR_->unboundPartial(typeI, typeJ) *= weight; // Full partial, summing bound and unbound terms - weightedgr.partial(typeI, typeJ).copyArrays(weightedgr.unboundPartial(typeI, typeJ)); - weightedgr.partial(typeI, typeJ) += weightedgr.boundPartial(typeI, typeJ); + weightedGR_->partial(typeI, typeJ).copyArrays(weightedGR_->unboundPartial(typeI, typeJ)); + weightedGR_->partial(typeI, typeJ) += weightedGR_->boundPartial(typeI, typeJ); } } // Calculate and normalise total to form factor if requested - weightedgr.formTotals(false); + weightedGR_->formTotals(false); // Normalise to Q=0.0 form factor if requested - if (normalisation != StructureFactors::NoNormalisation) + if (normaliseTo_ != StructureFactors::NoNormalisation) { - auto norm = normalisation == StructureFactors::AverageOfSquaresNormalisation ? weights.boundCoherentAverageOfSquares() - : weights.boundCoherentSquareOfAverage(); + auto norm = normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation ? weights_.boundCoherentAverageOfSquares() + : weights_.boundCoherentSquareOfAverage(); - weightedgr.total() /= norm; - weightedgr.boundTotal() /= norm; - weightedgr.unboundTotal() /= norm; + weightedGR_->total() /= norm; + weightedGR_->boundTotal() /= norm; + weightedGR_->unboundTotal() /= norm; } return true; } -// Calculate weighted S(Q) from supplied unweighted S(Q) and neutron weights -bool NeutronSQNode::calculateWeightedSQ(const PartialSet &unweightedsq, PartialSet &weightedsq, NeutronWeights &weights, - StructureFactors::NormalisationType normalisation) +// Calculate weighted S(Q) +bool NeutronSQNode::calculateWeightedSQ() { int typeI, typeJ; - for (typeI = 0; typeI < unweightedsq.nAtomTypes(); ++typeI) + for (typeI = 0; typeI < unweightedSQ_->nAtomTypes(); ++typeI) { - for (typeJ = typeI; typeJ < unweightedsq.nAtomTypes(); ++typeJ) + for (typeJ = typeI; typeJ < unweightedSQ_->nAtomTypes(); ++typeJ) { // Weight bound and unbound S(Q) and sum into full partial - double weight = weights.weight(typeI, typeJ); - double boundWeight = weights.intramolecularWeight(typeI, typeJ); + double weight = weights_.weight(typeI, typeJ); + double boundWeight = weights_.intramolecularWeight(typeI, typeJ); // Bound (intramolecular) partial (multiplied by the bound term weight) - weightedsq.boundPartial(typeI, typeJ).copyArrays(unweightedsq.boundPartial(typeI, typeJ)); - weightedsq.boundPartial(typeI, typeJ) *= boundWeight; + weightedSQ_->boundPartial(typeI, typeJ).copyArrays(unweightedSQ_->boundPartial(typeI, typeJ)); + weightedSQ_->boundPartial(typeI, typeJ) *= boundWeight; // Unbound partial (multiplied by the full weight) - weightedsq.unboundPartial(typeI, typeJ).copyArrays(unweightedsq.unboundPartial(typeI, typeJ)); - weightedsq.unboundPartial(typeI, typeJ) *= weight; + weightedSQ_->unboundPartial(typeI, typeJ).copyArrays(unweightedSQ_->unboundPartial(typeI, typeJ)); + weightedSQ_->unboundPartial(typeI, typeJ) *= weight; // Full partial (sum of bound and unbound terms) - weightedsq.partial(typeI, typeJ).copyArrays(weightedsq.unboundPartial(typeI, typeJ)); - weightedsq.partial(typeI, typeJ) += weightedsq.boundPartial(typeI, typeJ); + weightedSQ_->partial(typeI, typeJ).copyArrays(weightedSQ_->unboundPartial(typeI, typeJ)); + weightedSQ_->partial(typeI, typeJ) += weightedSQ_->boundPartial(typeI, typeJ); } } // Form total structure factor - weightedsq.formTotals(false); + weightedSQ_->formTotals(false); // Apply normalisation to all totals - if (normalisation != StructureFactors::NoNormalisation) + if (normaliseTo_ != StructureFactors::NoNormalisation) { - auto norm = normalisation == StructureFactors::AverageOfSquaresNormalisation ? weights.boundCoherentAverageOfSquares() - : weights.boundCoherentSquareOfAverage(); + auto norm = normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation ? weights_.boundCoherentAverageOfSquares() + : weights_.boundCoherentSquareOfAverage(); - weightedsq.total() /= norm; - weightedsq.boundTotal() /= norm; - weightedsq.unboundTotal() /= norm; + weightedSQ_->total() /= norm; + weightedSQ_->boundTotal() /= norm; + weightedSQ_->unboundTotal() /= norm; } return true; } -// Calculate neutron weights for relevant Configuration targets -void NeutronSQNode::calculateWeights(NeutronWeights &weights) const +// Calculate neutron weights matrix +void NeutronSQNode::calculateWeights(const std::map &realSpeciesPopulations) { - // Clear weights and get species speciesPopulations_ from GRModule - weights.clear(); - - for (auto &[sp, pop] : unweightedGR_->realSpeciesPopulations()) + // Create a set of named Isotopologues to use + IsotopologueSet topes; + for (const auto &[species, _] : realSpeciesPopulations) { - // Find the defined Isotopologue for this Species - if it doesn't exist, use the Natural one - auto isoRef = isotopologueSet_.getIsotopologues(sp); - if (isoRef) + for (const auto &isotopologue : species->isotopologues()) { - const Isotopologues &topes = *isoRef; - for (const auto &isoWeight : topes.mix()) - weights.addIsotopologue(sp, pop, isoWeight.isotopologue(), isoWeight.weight()); + auto iso = isotopologue.get(); + auto it = namedWeights_.find(iso->name()); + if (it != namedWeights_.end()) + topes.add(iso, it->second); } - else - weights.addIsotopologue(sp, pop, sp->naturalIsotopologue(), 1.0); } - weights.createFromIsotopologues(exchangeable_); + weights_.clear(); + + weights_.create(realSpeciesPopulations, topes, exchangeable_); } diff --git a/src/nodes/neutronSQ/neutronSQ.cpp b/src/nodes/neutronSQ/neutronSQ.cpp index 9973cc08df..8a86ab6d2e 100644 --- a/src/nodes/neutronSQ/neutronSQ.cpp +++ b/src/nodes/neutronSQ/neutronSQ.cpp @@ -13,10 +13,8 @@ NeutronSQNode::NeutronSQNode(Graph *parentGraph) : Node(parentGraph) { - addInput("UnweightedSQ", "Unweighted partials for target configuration", unweightedSQ_); - addInput("UnweightedGR", "Unweighted partials for target configuration", unweightedGR_); - addInput("Isotopologue", "Set/add an isotopologue and its population for a particular species", - isotopologueSet_); + addInput("UnweightedSQ", "Unweighted partial S(Q)", unweightedSQ_); + addInput("UnweightedGR", "Unweighted partials g(r)", unweightedGR_); addOption("NormaliseTo", "Normalisation to apply to total weighted F(Q)", normaliseTo_); addOption( diff --git a/src/nodes/neutronSQ/neutronSQ.h b/src/nodes/neutronSQ/neutronSQ.h index e2365bc2a6..8cf340cb59 100644 --- a/src/nodes/neutronSQ/neutronSQ.h +++ b/src/nodes/neutronSQ/neutronSQ.h @@ -47,8 +47,6 @@ class NeutronSQNode : public Node std::map namedWeights_{{"Ar36", 36}}; // Exchangeable atom types std::vector> exchangeable_; - // Isotopologues to use in weighting - IsotopologueSet isotopologueSet_; // Normalisation to apply to calculated total F(Q) StructureFactors::NormalisationType normaliseTo_{StructureFactors::NoNormalisation}; // Reference F(Q) file and format @@ -80,14 +78,12 @@ class NeutronSQNode : public Node * Functions */ public: - // Calculate weighted g(r) from supplied unweighted g(r) and neutron weights - bool calculateWeightedGR(const PartialSet &unweightedgr, PartialSet &weightedgr, NeutronWeights &weights, - StructureFactors::NormalisationType normalisation); - // Calculate weighted S(Q) from supplied unweighted S(Q) and neutron weights - bool calculateWeightedSQ(const PartialSet &unweightedsq, PartialSet &weightedsq, NeutronWeights &weights, - StructureFactors::NormalisationType normalisation); - // Calculate neutron weights for relevant Configuration targets - void calculateWeights(NeutronWeights &weights) const; + // Calculate weighted g(r) + bool calculateWeightedGR(); + // Calculate weighted S(Q) + bool calculateWeightedSQ(); + // Calculate neutron weights matrix + void calculateWeights(const std::map &realSpeciesPopulations); /* * Processing diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp index a3c9bc0c14..983aa898c4 100644 --- a/src/nodes/neutronSQ/process.cpp +++ b/src/nodes/neutronSQ/process.cpp @@ -132,18 +132,20 @@ NodeConstants::ProcessResult NeutronSQNode::process() /* * Transform UnweightedSQ from provided SQ data into WeightedSQ. */ + + // Get the real species populations from the input unweightedSQ + auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations(); + // Calculate and store weights_ /* auto& weights_ = dissolve.processingModuleData().realise("FullWeights", name(), GenericItem::InRestartFileFlag); */ - calculateWeights(weights_); + calculateWeights(realSpeciesPopulations); + message("Isotopologue and isotope composition:\n\n"); weights_.print(); - // Get the real species populations from the input unweightedSQ - auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations(); - // Does a PartialSet for the weighted S(Q) already exist for this Configuration? /* auto [weightedSQ, wSQstatus] = dissolve.processingModuleData().realiseIf( @@ -159,45 +161,17 @@ NodeConstants::ProcessResult NeutronSQNode::process() weightedSQ_->setUpPartials(unweightedSQ_->atomTypeMix()); } - // Update the isotopologue set - for (const auto &[species, _] : realSpeciesPopulations) - { - for (const auto &isotopologue : species->isotopologues()) - { - auto iso = isotopologue.get(); - auto it = namedWeights_.find(iso->name()); - if (it != namedWeights_.end()) - isotopologueSet_.add(iso, it->second); - } - } - // Calculate weighted S(Q) - calculateWeightedSQ(*unweightedSQ_, *weightedSQ_, weights_, normaliseTo_); + calculateWeightedSQ(); // Save data if requested - /* - if (saveSQ_ && (!MPIRunMaster(processPool(), weightedSQ.save(name(), "WeightedSQ", "sq", "Q, 1/Angstroms")))) + if (saveSQ_ && !weightedSQ_->save(name(), "WeightedSQ", "sq", "Q, 1/Angstroms")) return NodeConstants::ProcessResult::Failed; - */ /* - * Transform UnweightedGR from underlying RDF data into WeightedGR. + * Transform UnweightedGR from into WeightedGR. */ - // Get summed unweighted g(r) from the RDFMOdule - /* - if (!dissolve.processingModuleData().contains("UnweightedGR", rdfModule->name())) - { - error("Couldn't locate summed unweighted g(r) data.\n"); - return NodeConstants::ProcessResult::Failed; - } - */ - - /* - const auto& unweightedGR = - dissolve.processingModuleData().value("UnweightedGR", rdfModule->name()); - */ - // Create/retrieve PartialSet for summed weighted g(r) /* auto [weightedGR, wGRstatus] = dissolve.processingModuleData().realiseIf( @@ -214,13 +188,12 @@ NodeConstants::ProcessResult NeutronSQNode::process() } // Calculate weighted g(r) - calculateWeightedGR(*unweightedGR_, *weightedGR_, weights_, normaliseTo_); + calculateWeightedGR(); // Save data if requested - /* - if (saveGR_ && (!MPIRunMaster(processPool(), weightedGR.save(name(), "WeightedGR", "gr", "r, Angstroms")))) + if (saveGR_ && !weightedGR_->save(name(), "WeightedGR", "gr", "r, Angstroms")) return NodeConstants::ProcessResult::Failed; - */ + // Calculate representative total g(r) from FT of calculated F(Q) /* auto& repGR = dissolve.processingModuleData().realise("RepresentativeTotalGR", name(), From d3b4d6c230f257ded0a19deccf4f22b7a076cbb3 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 22 Jul 2025 14:31:29 +0100 Subject: [PATCH 6/9] Move code from (now redundant) set-up function into process() - had to do this anyway as we need NeutronWeights. --- src/nodes/neutronSQ/neutronSQ.h | 4 -- src/nodes/neutronSQ/process.cpp | 81 +++++++++++++++------------------ 2 files changed, 36 insertions(+), 49 deletions(-) diff --git a/src/nodes/neutronSQ/neutronSQ.h b/src/nodes/neutronSQ/neutronSQ.h index 8cf340cb59..7567441e2e 100644 --- a/src/nodes/neutronSQ/neutronSQ.h +++ b/src/nodes/neutronSQ/neutronSQ.h @@ -91,8 +91,4 @@ class NeutronSQNode : public Node private: // Run main processing NodeConstants::ProcessResult process() override; - - public: - // Run set-up stage - bool setUp(Flags actionSignals); }; diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp index 983aa898c4..6dd4c88a9c 100644 --- a/src/nodes/neutronSQ/process.cpp +++ b/src/nodes/neutronSQ/process.cpp @@ -15,27 +15,57 @@ #include "nodes/neutronSQ/neutronSQ.h" #include "nodes/sq/sq.h" -// Run set-up stage -bool NeutronSQNode::setUp(Flags actionSignals) +// Run main processing +NodeConstants::ProcessResult NeutronSQNode::process() { + // Print argument/parameter summary + if (referenceWindowFunction_ == WindowFunction::Form::None) + message("No window function will be applied when calculating representative g(r) from S(Q)."); + else + message("Window function to be applied when calculating representative g(r) from S(Q) is {}.", + WindowFunction::forms().keyword(referenceWindowFunction_)); + if (normaliseTo_ == StructureFactors::NoNormalisation) + message("NeutronSQ: No normalisation will be applied to total F(Q).\n"); + else if (normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation) + message("NeutronSQ: Total F(Q) will be normalised to "); + else if (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation) + message("NeutronSQ: Total F(Q) will be normalised to **2"); + if (saveSQ_) + message("NeutronSQ: Weighted partial S(Q) and total F(Q) will be saved.\n"); + if (saveGR_) + message("NeutronSQ: Weighted partial g(r) and total G(r) will be saved.\n"); + if (saveRepresentativeGR_) + message("NeutronSQ: Representative G(r) will be saved.\n"); + message("\n"); + + // Get the real species populations from the input unweightedSQ + auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations(); + + // Calculate and store weights_ + /* + auto& weights_ = dissolve.processingModuleData().realise("FullWeights", name(), + GenericItem::InRestartFileFlag); + */ + calculateWeights(realSpeciesPopulations); + message("Isotopologue and isotope composition:\n\n"); + weights_.print(); + /* * Load and set up reference data (if a file/format was given) */ - if (referenceFQ_.hasFilename() && actionSignals.isSetOrNone(KeywordBase::ReloadExternalData)) + if (referenceFQ_.hasFilename()) { // Load the data Data1D referenceData; if (!referenceFQ_.importData(referenceData)) { error("[SETUP {}] Failed to load reference data '{}'.\n", name(), referenceFQ_.filename()); - return false; + return NodeConstants::ProcessResult::Failed; } // Normalise reference data to be consistent with the calculated data if (referenceNormalisedTo_ != normaliseTo_) { - // We need the neutron weights_ in order to do the normalisation - calculateWeights(weights_); auto factor = 1.0; // Set up the multiplication factors @@ -103,49 +133,10 @@ bool NeutronSQNode::setUp(Flags actionSignals) } } - return true; -} - -// Run main processing -NodeConstants::ProcessResult NeutronSQNode::process() -{ - // Print argument/parameter summary - if (referenceWindowFunction_ == WindowFunction::Form::None) - message("No window function will be applied when calculating representative g(r) from S(Q)."); - else - message("Window function to be applied when calculating representative g(r) from S(Q) is {}.", - WindowFunction::forms().keyword(referenceWindowFunction_)); - if (normaliseTo_ == StructureFactors::NoNormalisation) - message("NeutronSQ: No normalisation will be applied to total F(Q).\n"); - else if (normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation) - message("NeutronSQ: Total F(Q) will be normalised to "); - else if (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation) - message("NeutronSQ: Total F(Q) will be normalised to **2"); - if (saveSQ_) - message("NeutronSQ: Weighted partial S(Q) and total F(Q) will be saved.\n"); - if (saveGR_) - message("NeutronSQ: Weighted partial g(r) and total G(r) will be saved.\n"); - if (saveRepresentativeGR_) - message("NeutronSQ: Representative G(r) will be saved.\n"); - message("\n"); - /* * Transform UnweightedSQ from provided SQ data into WeightedSQ. */ - // Get the real species populations from the input unweightedSQ - auto &realSpeciesPopulations = unweightedSQ_->realSpeciesPopulations(); - - // Calculate and store weights_ - /* - auto& weights_ = dissolve.processingModuleData().realise("FullWeights", name(), - GenericItem::InRestartFileFlag); - */ - calculateWeights(realSpeciesPopulations); - - message("Isotopologue and isotope composition:\n\n"); - weights_.print(); - // Does a PartialSet for the weighted S(Q) already exist for this Configuration? /* auto [weightedSQ, wSQstatus] = dissolve.processingModuleData().realiseIf( From 3beb008eae7f9fdb0dfca4721a813595940c83f3 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 23 Jul 2025 08:03:07 +0100 Subject: [PATCH 7/9] Correct return values. --- src/nodes/neutronSQ/process.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nodes/neutronSQ/process.cpp b/src/nodes/neutronSQ/process.cpp index 6dd4c88a9c..ad5f20b71f 100644 --- a/src/nodes/neutronSQ/process.cpp +++ b/src/nodes/neutronSQ/process.cpp @@ -126,10 +126,10 @@ NodeConstants::ProcessResult NeutronSQNode::process() { Data1DExportFileFormat exportFormat(std::format("{}-ReferenceData.q", name())); if (!exportFormat.exportData(storedData)) - return false; + return NodeConstants::ProcessResult::Failed; Data1DExportFileFormat exportFormatFT(std::format("{}-ReferenceData.r", name())); if (!exportFormatFT.exportData(storedDataFT)) - return false; + return NodeConstants::ProcessResult::Failed; } } From 7a1733b776836adbab2ef3a67a73cd160f3c3e46 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 23 Jul 2025 09:45:36 +0100 Subject: [PATCH 8/9] Format. --- src/nodes/gr/helpers.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nodes/gr/helpers.cpp b/src/nodes/gr/helpers.cpp index 4d68be44f4..13f5fbfe62 100644 --- a/src/nodes/gr/helpers.cpp +++ b/src/nodes/gr/helpers.cpp @@ -314,7 +314,7 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate) const auto &atoms = mol->atoms(); dissolve::for_each_pair(ParallelPolicies::seq, atoms, - [&, box](int index, auto &i, int jndex, auto &j) + [&, box](int index, auto &i, int jndex, auto &j) { // Ignore atom on itself if (index == jndex) From 05be970b4cc42bbcee52c30ae3f25146ed05a0e9 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Thu, 24 Jul 2025 09:23:22 +0100 Subject: [PATCH 9/9] Update src/nodes/gr/gr.cpp Co-authored-by: RobBuchanan <106311829+RobBuchananCompPhys@users.noreply.github.com> --- src/nodes/gr/gr.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nodes/gr/gr.cpp b/src/nodes/gr/gr.cpp index d7fbb8aaa4..2bf6698ace 100644 --- a/src/nodes/gr/gr.cpp +++ b/src/nodes/gr/gr.cpp @@ -5,7 +5,7 @@ GRNode::GRNode(Graph *parentGraph) : Node(parentGraph) { - addInput("Configurations", "Set target configuration(s) for the module", targetConfiguration_) + addInput("Configuration", "Set target configuration for the module", targetConfiguration_) ->setFlags({ParameterBase::Required, ParameterBase::ClearData}); addOption("BinWidth", "Bin width (spacing in r) to use", binWidth_); addOption>("Range", "Maximum r to calculate g(r) out to", requestedRange_);