Reinventing SiteRDF (#1778)

Co-authored-by: Tristan Youngs <tristan.youngs@stfc.ac.uk>

src/analyser/CMakeLists.txt

@@ -1,4 +1,13 @@
-add_library(analyser siteFilter.cpp siteFilter.h siteSelector.cpp siteSelector.h typeDefs.h)
+add_library(
+  analyser
+  dataNormaliser1D.cpp
+  dataNormaliser1D.h
+  siteFilter.cpp
+  siteFilter.h
+  siteSelector.cpp
+  siteSelector.h
+  typeDefs.h
+)
 
 target_include_directories(analyser PRIVATE ${PROJECT_SOURCE_DIR}/src ${PROJECT_BINARY_DIR}/src)
 

src/analyser/dataNormaliser1D.cpp

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "analyser/dataNormaliser1D.h"
+#include "classes/configuration.h"
+#include "math/data1D.h"
+
+DataNormaliser1D::DataNormaliser1D(Data1D &targetData) : targetData_(targetData) {}
+
+void DataNormaliser1D::normaliseByNumberDensity(double population, Configuration *targetConfiguration)
+{
+    targetData_ /= (population / targetConfiguration->box()->volume());
+}
+
+void DataNormaliser1D::normaliseBySitePopulation(double population) { targetData_ /= population; }
+
+void DataNormaliser1D::normaliseBySphericalShell()
+{
+    // We expect x values to be centre-bin values, and regularly spaced
+    const auto &xAxis = targetData_.xAxis();
+    auto &values = targetData_.values();
+
+    if (xAxis.size() < 2)
+        return;
+
+    // Derive first left-bin boundary from the delta between points 0 and 1
+    auto leftBin = xAxis[0] - (xAxis[1] - xAxis[0]) * 0.5;
+    auto r1Cubed = pow(leftBin, 3);
+    for (auto n = 0; n < xAxis.size(); ++n)
+    {
+        // Get new right-bin from existing left bin boundary and current bin centre
+        auto rightBin = leftBin + 2 * (xAxis[n] - leftBin);
+        auto r2Cubed = pow(rightBin, 3);
+
+        // Calculate divisor for normalisation
+        auto divisor = (4.0 / 3.0) * PI * (r2Cubed - r1Cubed);
+
+        // Peform normalisation step
+        values[n] /= divisor;
+        if (targetData_.valuesHaveErrors())
+            targetData_.error(n) /= divisor;
+
+        // Overwrite old values for next iteration
+        r1Cubed = r2Cubed;
+        leftBin = rightBin;
+    }
+}

src/analyser/dataNormaliser1D.h

@@ -0,0 +1,30 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#pragma once
+
+#include "analyser/typeDefs.h"
+
+class Configuration;
+class Data1D;
+
+class DataNormaliser1D
+{
+    public:
+    DataNormaliser1D(Data1D &targetData);
+
+    /*
+     * Target
+     */
+    private:
+    // Target data to normalise
+    Data1D &targetData_;
+
+    /*
+     * Normalisation functions
+     */
+    public:
+    void normaliseByNumberDensity(double population, Configuration *targetConfiguration);
+    void normaliseBySitePopulation(double population);
+    void normaliseBySphericalShell();
+};

src/modules/siteRDF/functions.cpp

@@ -2,23 +2,6 @@
 // Copyright (c) 2024 Team Dissolve and contributors
 
 #include "modules/siteRDF/siteRDF.h"
-#include "procedure/nodes/process1D.h"
-#include "procedure/nodes/sum1D.h"
-
-// Return Collect1DNode for A-B RDF
-std::shared_ptr<Collect1DProcedureNode> SiteRDFModule::collectDistanceNode() const { return collectDistance_; }
-
-// Return SelectNode for site A
-std::shared_ptr<SelectProcedureNode> SiteRDFModule::selectANode() const { return selectA_; }
 
 // Return whether specified coordination number range is enabled
-bool SiteRDFModule::isRangeEnabled(int id) const { return sumCN_ && sumCN_->rangeEnabled(id); }
-
-// Return Process1DNode result (i.e. RDF)
-std::shared_ptr<Process1DProcedureNode> SiteRDFModule::rdfResult() const
-{
-    if ((!processDistance_) || (!processDistance_->hasProcessedData()))
-        return nullptr;
-
-    return processDistance_;
-}
+bool SiteRDFModule::isRangeEnabled(int id) const { return rangeEnabled_[id]; }

src/modules/siteRDF/gui/siteRDFWidgetFuncs.cpp

@@ -34,17 +34,14 @@ void SiteRDFModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFlags>
     // Update CN labels
     auto rangeAOn = module_->isRangeEnabled(0);
     ui_.RegionAResultFrame->setText(
-        rangeAOn ? dissolve_.processingModuleData().valueOr("Sum1D//CN//A", module_->name(), SampledDouble())
-                 : SampledDouble());
+        rangeAOn ? dissolve_.processingModuleData().valueOr("CN//A", module_->name(), SampledDouble()) : SampledDouble());
     auto rangeBOn = module_->isRangeEnabled(1);
     ui_.RegionBResultFrame->setText(
-        rangeBOn ? dissolve_.processingModuleData().valueOr("Sum1D//CN//B", module_->name(), SampledDouble())
-                 : SampledDouble());
+        rangeBOn ? dissolve_.processingModuleData().valueOr("CN//B", module_->name(), SampledDouble()) : SampledDouble());
     ui_.RegionBResultFrame->setEnabled(rangeBOn);
     auto rangeCOn = module_->isRangeEnabled(2);
     ui_.RegionCResultFrame->setText(
-        rangeCOn ? dissolve_.processingModuleData().valueOr("Sum1D//CN//C", module_->name(), SampledDouble())
-                 : SampledDouble());
+        rangeCOn ? dissolve_.processingModuleData().valueOr("CN//C", module_->name(), SampledDouble()) : SampledDouble());
     ui_.RegionCResultFrame->setEnabled(rangeCOn);
 
     if (updateFlags.isSet(ModuleWidget::RecreateRenderablesFlag))
@@ -55,7 +52,7 @@ void SiteRDFModuleWidget::updateControls(const Flags<ModuleWidget::UpdateFlags>
         auto *cfg = module_->keywords().getConfiguration("Configuration");
         if (cfg)
             rdfGraph_
-                ->createRenderable<RenderableData1D>(fmt::format("{}//Process1D//RDF", module_->name()),
+                ->createRenderable<RenderableData1D>(fmt::format("{}//RDF", module_->name()),
                                                      fmt::format("RDF//{}", cfg->niceName()), cfg->niceName())
                 ->setColour(StockColours::BlueStockColour);
     }

src/modules/siteRDF/process.cpp

@@ -1,16 +1,14 @@
 // SPDX-License-Identifier: GPL-3.0-or-later
 // Copyright (c) 2024 Team Dissolve and contributors
 
+#include "analyser/dataNormaliser1D.h"
 #include "base/sysFunc.h"
 #include "io/export/data1D.h"
 #include "main/dissolve.h"
+#include "math/integrator.h"
+#include "math/sampledDouble.h"
 #include "module/context.h"
 #include "modules/siteRDF/siteRDF.h"
-#include "procedure/nodes/collect1D.h"
-#include "procedure/nodes/operateSitePopulationNormalise.h"
-#include "procedure/nodes/select.h"
-#include "procedure/nodes/sequence.h"
-#include "procedure/nodes/sum1D.h"
 
 // Run main processing
 Module::ExecutionResult SiteRDFModule::process(ModuleContext &moduleContext)
@@ -22,68 +20,85 @@ Module::ExecutionResult SiteRDFModule::process(ModuleContext &moduleContext)
         return ExecutionResult::Failed;
     }
 
-    // Ensure any parameters in our nodes are set correctly
-    collectDistance_->keywords().set("RangeX", distanceRange_);
-    if (excludeSameMolecule_)
-        selectB_->setSameMoleculeExclusions({selectA_});
-    else
-        selectB_->setSameMoleculeExclusions({});
-    cnNormaliser_->keywords().set("Site", ConstNodeVector<SelectProcedureNode>{selectA_});
+    auto &processingData = moduleContext.dissolve().processingModuleData();
 
-    // Execute the analysis
-    if (!analyser_.execute({moduleContext.dissolve(), targetConfiguration_, name()}))
-    {
-        Messenger::error("CalculateRDF experienced problems with its analysis.\n");
-        return ExecutionResult::Failed;
-    }
+    // Select site A
+    SiteSelector a(targetConfiguration_, a_);
 
-    // Accumulate instantaneous coordination number data
-    if (instantaneous_)
+    // Select site B
+    SiteSelector b(targetConfiguration_, b_);
+
+    // Calculate rAB
+    auto [histAB, status] = processingData.realiseIf<Histogram1D>("Histo-AB", name(), GenericItem::InRestartFileFlag);
+    if (status == GenericItem::ItemStatus::Created)
+        histAB.initialise(distanceRange_.x, distanceRange_.y, distanceRange_.z);
+
+    histAB.zeroBins();
+    for (const auto &[siteA, indexA] : a.sites())
     {
-        if (isRangeEnabled(0))
+        for (const auto &[siteB, indexB] : b.sites())
         {
-            auto &sumA =
-                moduleContext.dissolve().processingModuleData().realise<Data1D>("SumA", name(), GenericItem::InRestartFileFlag);
-            sumA.addPoint(moduleContext.dissolve().iteration(), sumCN_->sum(0).value());
-            if (exportInstantaneous_)
-            {
-                Data1DExportFileFormat exportFormat(fmt::format("{}_SumA.txt", name()));
-                if (!exportFormat.exportData(sumA))
-                {
-                    Messenger::error("Failed to write instantaneous coordination number data for range A.\n");
-                    return ExecutionResult::Failed;
-                }
-            }
+            if (excludeSameMolecule_ && (siteB->molecule() == siteA->molecule()))
+                continue;
+            histAB.bin(targetConfiguration_->box()->minimumDistance(siteA->origin(), siteB->origin()));
         }
-        if (isRangeEnabled(1))
+    }
+    histAB.accumulate();
+
+    // RDF
+    auto &dataRDF = processingData.realise<Data1D>("RDF", name(), GenericItem::InRestartFileFlag);
+    dataRDF = histAB.accumulatedData();
+    DataNormaliser1D normaliserRDF(dataRDF);
+
+    // Normalise by A site population
+    normaliserRDF.normaliseBySitePopulation(double(a.sites().size()));
+
+    // Normalise by B site population density
+    normaliserRDF.normaliseByNumberDensity(double(b.sites().size()), targetConfiguration_);
+
+    // Normalise by spherical shell
+    normaliserRDF.normaliseBySphericalShell();
+
+    // CN
+    auto &dataCN = processingData.realise<Data1D>("HistogramNorm", name(), GenericItem::InRestartFileFlag);
+    dataCN = histAB.accumulatedData();
+    DataNormaliser1D normaliserCN(dataCN);
+    normaliserCN.normaliseBySitePopulation(double(a.sites().size()));
+
+    const std::vector<std::string> rangeNames = {"A", "B", "C"};
+    for (int i = 0; i < 3; ++i)
+        if (rangeEnabled_[i])
         {
-            auto &sumB =
-                moduleContext.dissolve().processingModuleData().realise<Data1D>("SumB", name(), GenericItem::InRestartFileFlag);
-            sumB.addPoint(moduleContext.dissolve().iteration(), sumCN_->sum(1).value());
-            if (exportInstantaneous_)
+            auto &sumN = processingData.realise<SampledDouble>(fmt::format("CN//{}", rangeNames[i]), name(),
+                                                               GenericItem::InRestartFileFlag);
+            sumN += Integrator::sum(dataCN, range_[i]);
+            if (instantaneous_)
             {
-                Data1DExportFileFormat exportFormat(fmt::format("{}_SumB.txt", name()));
-                if (!exportFormat.exportData(sumB))
+                auto &sumNInst = processingData.realise<Data1D>(fmt::format("CN//{}Inst", rangeNames[i]), name(),
+                                                                GenericItem::InRestartFileFlag);
+                sumNInst.addPoint(moduleContext.dissolve().iteration(), sumN.value());
+                if (exportInstantaneous_)
                 {
-                    Messenger::error("Failed to write instantaneous coordination number data for range B.\n");
-                    return ExecutionResult::Failed;
+                    Data1DExportFileFormat exportFormat(fmt::format("{}_Sum{}.txt", name(), rangeNames[i]));
+                    if (!exportFormat.exportData(sumNInst))
+                    {
+                        Messenger::error("Failed to write instantaneous coordination number data for range {}.\n",
+                                         rangeNames[i]);
+                        return ExecutionResult::Failed;
+                    }
                 }
             }
         }
-        if (isRangeEnabled(2))
+
+    // Save data?
+    if (exportFileAndFormat_.hasFilename())
+    {
+        if (moduleContext.processPool().isMaster())
         {
-            auto &sumC =
-                moduleContext.dissolve().processingModuleData().realise<Data1D>("SumC", name(), GenericItem::InRestartFileFlag);
-            sumC.addPoint(moduleContext.dissolve().iteration(), sumCN_->sum(2).value());
-            if (exportInstantaneous_)
-            {
-                Data1DExportFileFormat exportFormat(fmt::format("{}_SumC.txt", name()));
-                if (!exportFormat.exportData(sumC))
-                {
-                    Messenger::error("Failed to write instantaneous coordination number data for range C.\n");
-                    return ExecutionResult::Failed;
-                }
-            }
+            if (exportFileAndFormat_.exportData(dataCN))
+                moduleContext.processPool().decideTrue();
+            else
+                moduleContext.processPool().decideFalse();
         }
     }