From b23fbb89b421d8fdceba666a3a909834e8d87d1c Mon Sep 17 00:00:00 2001 From: Benjamin Speake Date: Wed, 15 Oct 2025 16:36:57 +0100 Subject: [PATCH 1/2] Create CGNeutronSQ module and add functions to calculate Dodecane form factor and single bead F(Q) --- src/module/module.cpp | 1 + src/module/types.h | 1 + src/modules/cgNeutronSQ/CMakeLists.txt | 1 + src/modules/cgNeutronSQ/cgNeutronSQ.cpp | 79 +++++ src/modules/cgNeutronSQ/cgNeutronSQ.h | 93 ++++++ src/modules/cgNeutronSQ/functions.cpp | 202 ++++++++++++ src/modules/cgNeutronSQ/gui/CMakeLists.txt | 1 + .../cgNeutronSQ/gui/cgNeutronSQWidget.h | 58 ++++ .../cgNeutronSQ/gui/cgNeutronSQWidget.ui | 238 ++++++++++++++ .../gui/cgNeutronSQWidgetFuncs.cpp | 214 +++++++++++++ src/modules/cgNeutronSQ/process.cpp | 290 ++++++++++++++++++ src/modules/registry.cpp | 3 + src/modules/widgetProducer.cpp | 3 + 13 files changed, 1184 insertions(+) create mode 100644 src/modules/cgNeutronSQ/CMakeLists.txt create mode 100644 src/modules/cgNeutronSQ/cgNeutronSQ.cpp create mode 100644 src/modules/cgNeutronSQ/cgNeutronSQ.h create mode 100644 src/modules/cgNeutronSQ/functions.cpp create mode 100644 src/modules/cgNeutronSQ/gui/CMakeLists.txt create mode 100644 src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.h create mode 100644 src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.ui create mode 100644 src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp create mode 100644 src/modules/cgNeutronSQ/process.cpp diff --git a/src/module/module.cpp b/src/module/module.cpp index 8a1729e477..fbc4a50f6b 100644 --- a/src/module/module.cpp +++ b/src/module/module.cpp @@ -24,6 +24,7 @@ EnumOptions moduleTypes_("ModuleType", {{ModuleTypes::A {ModuleTypes::Checks, "Checks"}, {ModuleTypes::CheckSpecies, "CheckSpecies"}, {ModuleTypes::Clustering, "Clustering"}, + {ModuleTypes::CGNeutronSQ, "CGNeutronSQ"}, {ModuleTypes::Compare, "Compare"}, {ModuleTypes::DAngle, "DAngle"}, {ModuleTypes::DataTest, "DataTest"}, diff --git a/src/module/types.h b/src/module/types.h index 0da5470c00..a6fca7de31 100644 --- a/src/module/types.h +++ b/src/module/types.h @@ -18,6 +18,7 @@ enum ModuleType Checks, CheckSpecies, Clustering, + CGNeutronSQ, Compare, DAngle, DataTest, diff --git a/src/modules/cgNeutronSQ/CMakeLists.txt b/src/modules/cgNeutronSQ/CMakeLists.txt new file mode 100644 index 0000000000..e0b17a75cc --- /dev/null +++ b/src/modules/cgNeutronSQ/CMakeLists.txt @@ -0,0 +1 @@ +dissolve_add_module(cgNeutronSQ.h cgNeutronSQ) diff --git a/src/modules/cgNeutronSQ/cgNeutronSQ.cpp b/src/modules/cgNeutronSQ/cgNeutronSQ.cpp new file mode 100644 index 0000000000..f4249fe67d --- /dev/null +++ b/src/modules/cgNeutronSQ/cgNeutronSQ.cpp @@ -0,0 +1,79 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "keywords/atomTypeVector.h" +#include "keywords/bool.h" +#include "keywords/double.h" +#include "keywords/fileAndFormat.h" +#include "keywords/isotopologueSet.h" +#include "keywords/module.h" +#include "keywords/optionalDouble.h" +#include "modules/sq/sq.h" + +CGNeutronSQModule::CGNeutronSQModule() : Module(ModuleTypes::NeutronSQ) +{ + keywords_.addTarget>( + "SourceSQs", "Source unweighted S(Q) to transform into neutron-weighted S(Q)", sourceSQ_, ModuleTypes::SQ); + + keywords_.setOrganisation("Options", "Isotopes & Normalisation", + "Specify isotopologues to use for specific species, and which atoms are exchangeable."); + keywords_.add( + "Exchangeable", "A set of atom types in the system that are exchangeable with each other", exchangeable_); + keywords_.add("Isotopologue", "Set/add an isotopologue and its population for a particular species", + isotopologueSet_); + keywords_ + .add>("NormaliseTo", + "Normalisation to apply to total weighted F(Q)", + normaliseTo_, StructureFactors::normalisationTypes()) + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + + keywords_.setOrganisation("Options", "Reference Data", + "Reference (typically experimental) data set to display. The experimental data may be used by " + "other modules if present. The normalisation already applied to the reference data should be " + "specified here, and will be removed internally."); + keywords_.add("Reference", "F(Q) reference data", referenceFQ_, "EndReference") + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + keywords_ + .add>( + "ReferenceNormalisedTo", "Normalisation that has been applied to the reference data", referenceNormalisedTo_, + StructureFactors::normalisationTypes()) + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + keywords_ + .add("ReferenceFTQMin", + "Minimum Q value to use when Fourier-transforming reference data (0.0 for no minimum)", + referenceFTQMin_, 0.0, std::nullopt, 0.1, "No Minimum Limit") + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + keywords_ + .add("ReferenceFTQMax", + "Maximum Q value to use when Fourier-transforming reference data (0.0 for no maximum)", + referenceFTQMax_, 0.0, std::nullopt, 0.1, "No Maximum Limit") + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + keywords_ + .add("ReferenceFTDeltaR", "Spacing in r to use when generating the Fourier-transformed data", + referenceFTDeltaR_, 1.0e-4, 1.0) + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + keywords_ + .add>( + "ReferenceWindowFunction", "Window function to apply when Fourier-transforming reference S(Q) to g(r)", + referenceWindowFunction_, WindowFunction::forms()) + ->setEditSignals({KeywordBase::ReloadExternalData, KeywordBase::RecreateRenderables}); + + keywords_.setOrganisation("Export"); + keywords_.add("SaveGR", "Save weighted g(r) and G(r)", saveGR_); + keywords_.add("SaveReference", "Save the reference data and its Fourier transform", saveReference_); + keywords_.add("SaveRepresentativeGR", + "Save representative G(r), obtained from Fourier transform of the calculated F(Q)", + saveRepresentativeGR_); + keywords_.add("SaveSQ", "Save weighted partial and total structure factors", saveSQ_); + + // Deprecated keywords + keywords_.addDeprecated>( + "Normalisation", "Normalisation to apply to total weighted F(Q)", normaliseTo_, StructureFactors::normalisationTypes()); + keywords_.addDeprecated>( + "ReferenceNormalisation", "Normalisation to remove from reference data before use", referenceNormalisedTo_, + StructureFactors::normalisationTypes()); +} + +// Return file and format for reference total F(Q) +const Data1DImportFileFormat &CGNeutronSQModule::referenceFQFileAndFormat() { return referenceFQ_; } diff --git a/src/modules/cgNeutronSQ/cgNeutronSQ.h b/src/modules/cgNeutronSQ/cgNeutronSQ.h new file mode 100644 index 0000000000..8b6d3afb59 --- /dev/null +++ b/src/modules/cgNeutronSQ/cgNeutronSQ.h @@ -0,0 +1,93 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#pragma once + +#include "classes/isotopologueSet.h" +#include "classes/partialSet.h" +#include "data/structureFactors.h" +#include "io/import/data1D.h" +#include "math/windowFunction.h" +#include "module/module.h" + +// Forward Declarations +class PartialSet; +class GRModule; +class SQModule; + +// Neutron SQ Module +class CGNeutronSQModule : public Module +{ + public: + CGNeutronSQModule(); + ~CGNeutronSQModule() override = default; + + /* + * Definition + */ + private: + // 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 + Data1DImportFileFormat referenceFQ_{"", Data1DImportFileFormat::Data1DImportFormat::GudrunMint}; + // Minimum Q value to use when Fourier-transforming the data + std::optional referenceFTQMin_{0.3}; + // Maximum Q value to use when Fourier-transforming the data + std::optional referenceFTQMax_{30.0}; + // Spacing in r to use when generating the Fourier-transformed data + double referenceFTDeltaR_{0.05}; + // Normalisation to remove from reference total F(Q) + StructureFactors::NormalisationType referenceNormalisedTo_{StructureFactors::NoNormalisation}; + // Window function to use when Fourier transforming reference total F(Q) into g(r) + WindowFunction::Form referenceWindowFunction_{WindowFunction::Form::Lorch0}; + // Save weighted g(r) and G(r) + bool saveGR_{false}; + // Save the reference data and its Fourier transform + bool saveReference_{false}; + // Save representative G(r), obtained from Fourier transform of the calculated F(Q) + bool saveRepresentativeGR_{false}; + // Save weighted partial and total structure factors + bool saveSQ_{false}; + // Source module for calculation + const SQModule *sourceSQ_{nullptr}; + + public: + // Return file and format for reference total F(Q) + const Data1DImportFileFormat &referenceFQFileAndFormat(); + + /* + * 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, + const std::vector &ff, const std::vector &singleBead, + StructureFactors::NormalisationType normalisation); + // Calculate neutron weights for relevant Configuration targets + void calculateWeights(const GRModule *rdfModule, NeutronWeights &weights) const; + // + bool calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, const NeutronWeights &weights) const; + // + bool calculateSingleBead(std::vector &singleBead, const std::vector &ff, const NeutronWeights &weights) const; + + /* + * Processing + */ + private: + // Run main processing + Module::ExecutionResult process(ModuleContext &moduleContext) override; + + public: + // Set target data + void setTargets(const std::vector> &configurations, + const std::map> &moduleMap) override; + // Run set-up stage + bool setUp(ModuleContext &moduleContext, Flags actionSignals) override; +}; diff --git a/src/modules/cgNeutronSQ/functions.cpp b/src/modules/cgNeutronSQ/functions.cpp new file mode 100644 index 0000000000..a554678b5f --- /dev/null +++ b/src/modules/cgNeutronSQ/functions.cpp @@ -0,0 +1,202 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#include + +#include "classes/configuration.h" +#include "classes/isotopologueSet.h" +#include "classes/species.h" +#include "modules/gr/gr.h" +#include "modules/cgNeutronSQ/cgNeutronSQ.h" + +// Calculate weighted g(r) from supplied unweighted g(r) and neutron weights +bool CGNeutronSQModule::calculateWeightedGR(const PartialSet &unweightedgr, PartialSet &weightedgr, NeutronWeights &weights, + StructureFactors::NormalisationType normalisation) +{ + int typeI, typeJ; + for (typeI = 0; typeI < unweightedgr.nAtomTypes(); ++typeI) + { + for (typeJ = typeI; typeJ < unweightedgr.nAtomTypes(); ++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; + + // 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; + + // Full partial, summing bound and unbound terms + 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); + + // Normalise to Q=0.0 form factor if requested + if (normalisation != StructureFactors::NoNormalisation) + { + auto norm = normalisation == StructureFactors::AverageOfSquaresNormalisation ? weights.boundCoherentAverageOfSquares() + : weights.boundCoherentSquareOfAverage(); + + weightedgr.total() /= norm; + weightedgr.boundTotal() /= norm; + weightedgr.unboundTotal() /= norm; + } + + return true; +} + +// Calculate weighted S(Q) from supplied unweighted S(Q) and neutron weights +bool CGNeutronSQModule::calculateWeightedSQ(const PartialSet &unweightedsq, PartialSet &weightedsq, NeutronWeights &weights, + const std::vector &ff, const std::vector &singleBead, + StructureFactors::NormalisationType normalisation) +{ + int typeI, typeJ; + for (typeI = 0; typeI < unweightedsq.nAtomTypes(); ++typeI) + { + 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); + + // 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) *= ff[typeI].values(); + weightedsq.boundPartial(typeI, typeJ) *= ff[typeJ].values(); + + // 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) *= ff[typeI].values(); + weightedsq.unboundPartial(typeI, typeJ) *= ff[typeJ].values(); + + // 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); + } + } + + // Form total structure factor + weightedsq.formTotals(false); + for (typeI = 0; typeI < unweightedsq.nAtomTypes(); ++typeI) { + weightedsq.total() += singleBead[typeI]; + } + + // Apply normalisation to all totals + if (normalisation != StructureFactors::NoNormalisation) + { + auto norm = normalisation == StructureFactors::AverageOfSquaresNormalisation ? weights.boundCoherentAverageOfSquares() + : weights.boundCoherentSquareOfAverage(); + + weightedsq.total() /= norm; + weightedsq.boundTotal() /= norm; + weightedsq.unboundTotal() /= norm; + } + + return true; +} + +// Calculate neutron weights for relevant Configuration targets +void CGNeutronSQModule::calculateWeights(const GRModule *rdfModule, NeutronWeights &weights) const +{ + // Clear weights and get species populations from GRModule + weights.clear(); + auto populations = rdfModule->speciesPopulations(); + + for (auto &[sp, pop] : populations) + { + // Find the defined Isotopologue for this Species - if it doesn't exist, use the Natural one + auto isoRef = isotopologueSet_.getIsotopologues(sp); + if (isoRef) + { + const Isotopologues &topes = *isoRef; + for (const auto &isoWeight : topes.mix()) + weights.addIsotopologue(sp, pop, isoWeight.isotopologue(), isoWeight.weight()); + } + else + weights.addIsotopologue(sp, pop, sp->naturalIsotopologue(), 1.0); + } + + weights.createFromIsotopologues(exchangeable_); +} + +// +bool CGNeutronSQModule::calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, const NeutronWeights &weights) const +{ + int typeI; + int dens{0}; + //std::array sigma{1.63, 1.57}; + std::array sigma{3.0376754356895721, 3.7573162947024161}; + double f; + ff.clear(); + ff.reserve(qvals.size()); + + for (typeI = 0; typeI < weights.atomTypes().nItems(); ++typeI) + { + ff.emplace_back(Data1D()); + for (double q : qvals) + { + double qs = q * sigma[typeI]; + if (dens == 0) // gaussian + { + f = -0.5 * std::pow((0.51 * qs), 2); + f = std::exp(f); + ff[typeI].addPoint(q, f); + } + else if (dens == 1) // uniform + { + f = 3.0 / std::pow((qs), 3); + f *= (std::sin(qs) - qs * std::cos(qs)); + ff[typeI].addPoint(q, f); + } + else // no form factor + { + ff[typeI].addPoint(q, 1.0); + } + } + } + return true; +} + +// +bool CGNeutronSQModule::calculateSingleBead(std::vector &singleBead, const std::vector &ff, const NeutronWeights &weights) const { + + singleBead.clear(); + singleBead.reserve(weights.atomTypes().nItems()); + std::array b{0.6646, -0.3739}; + std::array, 2> n; + n[0][0] = 3; + n[0][1] = 7; + n[1][0] = 3; + n[1][1] = 6; + int atmI, atmJ; + for (auto bead = weights.atomTypes().begin(); bead != weights.atomTypes().end(); ++bead) + { + int s = bead - weights.atomTypes().begin(); + double innerSum{0.0}; + singleBead.emplace_back(Data1D()); + for (atmI = 0; atmI < 2; ++atmI) + { + for (atmJ = atmI; atmJ < 2; ++atmJ) + { + innerSum += n[s][atmI] * n[s][atmJ] * b[atmI] * b[atmJ] * (atmI == atmJ ? 1 : 2); + } + } + innerSum -= ((n[s][0] * b[0] * b[0]) + (n[s][1] * b[1] * b[1])); + innerSum *= 0.5; // (*bead).fraction(); + singleBead[s].copyArrays(ff[s]); + singleBead[s] *= ff[s].values(); + singleBead[s] *= innerSum * (1/9.5); + } + + return true; +} \ No newline at end of file diff --git a/src/modules/cgNeutronSQ/gui/CMakeLists.txt b/src/modules/cgNeutronSQ/gui/CMakeLists.txt new file mode 100644 index 0000000000..75ee6a560c --- /dev/null +++ b/src/modules/cgNeutronSQ/gui/CMakeLists.txt @@ -0,0 +1 @@ +dissolve_add_module_gui(cgNeutronSQ) diff --git a/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.h b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.h new file mode 100644 index 0000000000..e321ca12b5 --- /dev/null +++ b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.h @@ -0,0 +1,58 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#pragma once + +#include "modules/neutronSQ/gui/ui_cgNeutronSQWidget.h" +#include "modules/widget.h" + +// Forward Declarations +class Dissolve; +class CGNeutronSQModule; +class PartialSet; +class DataViewer; + +// Module Widget +class CGNeutronSQModuleWidget : public ModuleWidget +{ + // All Qt declarations derived from QObject must include this macro + Q_OBJECT + + public: + CGNeutronSQModuleWidget(QWidget *parent, CGNeutronSQModule *module, Dissolve &dissolve); + ~CGNeutronSQModuleWidget() override = default; + + private: + // Associated Module + CGNeutronSQModule *module_; + // Target partial data being displayed (if any) + OptionalReferenceWrapper targetPartials_; + + /* + * UI + */ + private: + // Main form declaration + Ui::CGNeutronSQModuleWidget ui_; + // DataViewers contained within this widget + DataViewer *graph_; + + private: + // Create renderables for current target PartialSet + void createPartialSetRenderables(std::string_view targetPrefix); + + public: + // Update controls within widget + void updateControls(const Flags &updateFlags = {}) override; + + /* + * Widgets / Functions + */ + private Q_SLOTS: + void on_TotalFQButton_clicked(bool checked); + void on_PartialSQButton_clicked(bool checked); + void on_TotalGRButton_clicked(bool checked); + void on_PartialGRButton_clicked(bool checked); + void on_FilterEdit_textChanged(QString text); + void on_ClearFilterButton_clicked(bool checked); +}; diff --git a/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.ui b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.ui new file mode 100644 index 0000000000..dfdb40da96 --- /dev/null +++ b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidget.ui @@ -0,0 +1,238 @@ + + + CGNeutronSQModuleWidget + + + + 0 + 0 + 846 + 631 + + + + + 0 + 0 + + + + + 8 + + + + NeutronSQ Controls + + + + 4 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + 0 + + + + + Show the neutron-weighted total F(Q) formed from the weighted partial S(Q) + + + Total F(Q) + + + true + + + true + + + true + + + + + + + Show the neutron-weighted partial S(Q) + + + Partial S(Q) + + + true + + + false + + + true + + + + + + + Show the total neutron-weighted G(r) formed from summation of the weighted partial g(r) + + + Total G(r) + + + true + + + false + + + true + + + + + + + Show the neutron-weighted partial g(r) + + + Partial g(r) + + + true + + + false + + + true + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + + 0 + 0 + + + + + 16 + 16 + + + + + + + :/general/icons/filter.svg + + + true + + + + + + + Filter displayed partials by text / atomtype + + + + + + + Clear partial filter + + + + + + + :/general/icons/cross.svg:/general/icons/cross.svg + + + + + + + + + + 0 + 0 + + + + true + + + QFrame::StyledPanel + + + QFrame::Sunken + + + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + + + + + + DataWidget + QWidget +
gui/dataWidget.h
+ 1 + + + + + + + diff --git a/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp new file mode 100644 index 0000000000..1016fd4a44 --- /dev/null +++ b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp @@ -0,0 +1,214 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#include "classes/atomType.h" +#include "classes/isotopeData.h" +#include "gui/dataViewer.h" +#include "gui/render/renderableData1D.h" +#include "main/dissolve.h" +#include "modules/cgNeutronSQ/gui/cgNeutronSQWidget.h" +#include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "templates/algorithms.h" + +CGNeutronSQModuleWidget::CGNeutronSQModuleWidget(QWidget *parent, CGNeutronSQModule *module, Dissolve &dissolve) + : ModuleWidget(parent, dissolve), module_(module) +{ + // Set up user interface + ui_.setupUi(this); + + // Set up graph (defaulting to total F(Q)) + graph_ = ui_.PlotWidget->dataViewer(); + // -- Set view + graph_->view().setViewType(View::FlatXYView); + graph_->view().axes().setTitle(0, "\\it{Q}, \\sym{angstrom}\\sup{-1}"); + graph_->view().axes().setMax(0, 10.0); + graph_->view().axes().setTitle(1, "F(Q)"); + graph_->view().axes().setMin(1, -1.0); + graph_->view().axes().setMax(1, 1.0); + graph_->view().setAutoFollowType(View::AllAutoFollow); + // -- Set group styling + graph_->groupManager().setGroupColouring("Full", RenderableGroup::AutomaticIndividualColouring); + graph_->groupManager().setGroupVerticalShifting("Full", RenderableGroup::IndividualVerticalShifting); + graph_->groupManager().setGroupColouring("Bound", RenderableGroup::AutomaticIndividualColouring); + graph_->groupManager().setGroupVerticalShifting("Bound", RenderableGroup::IndividualVerticalShifting); + graph_->groupManager().setGroupStipple("Bound", LineStipple::HalfDashStipple); + graph_->groupManager().setGroupColouring("Unbound", RenderableGroup::AutomaticIndividualColouring); + graph_->groupManager().setGroupVerticalShifting("Unbound", RenderableGroup::IndividualVerticalShifting); + graph_->groupManager().setGroupStipple("Unbound", LineStipple::DotStipple); + + refreshing_ = false; +} + +/* + * UI + */ + +// Create renderables for current target PartialSet +void CGNeutronSQModuleWidget::createPartialSetRenderables(std::string_view targetPrefix) +{ + if (!targetPartials_) + return; + + const PartialSet &ps = *targetPartials_; + + // Get the filter text (if there is any) + std::optional filterText; + if (!ui_.FilterEdit->text().isEmpty()) + filterText = ui_.FilterEdit->text().toStdString(); + + PairIterator pairs(ps.atomTypeMix().nItems()); + for (auto [first, second] : pairs) + { + auto &at1 = ps.atomTypeMix()[first]; + auto &at2 = ps.atomTypeMix()[second]; + const std::string id = std::format("{}-{}", at1.atomTypeName(), at2.atomTypeName()); + + // Filtering - does this 'id' match our filter? + if (filterText && id.find(filterText.value()) == std::string::npos) + continue; + + // Full partial + graph_->createRenderable(std::format("{}//{}//{}//Full", module_->name(), targetPrefix, id), + std::format("{} (Full)", id), "Full"); + + // Bound partial + graph_->createRenderable(std::format("{}//{}//{}//Bound", module_->name(), targetPrefix, id), + std::format("{} (Bound)", id), "Bound"); + + // Unbound partial + graph_->createRenderable(std::format("{}//{}//{}//Unbound", module_->name(), targetPrefix, id), + std::format("{} (Unbound)", id), "Unbound"); + } +} + +// Update controls within widget +void CGNeutronSQModuleWidget::updateControls(const Flags &updateFlags) +{ + refreshing_ = true; + + // Need to recreate renderables if requested as the updateType, or if we previously had no target PartialSet and have just + // located it + if (updateFlags.isSet(ModuleWidget::RecreateRenderablesFlag) || + (!ui_.TotalFQButton->isChecked() && !ui_.TotalGRButton->isChecked() && !targetPartials_)) + { + ui_.PlotWidget->clearRenderableData(); + + // Grab reference data file and format + const Data1DImportFileFormat &referenceFileAndFormat = module_->referenceFQFileAndFormat(); + + if (ui_.TotalFQButton->isChecked()) + { + graph_->createRenderable(std::format("{}//WeightedSQ//Total", module_->name()), "Total F(Q)", + "Calculated"); + auto boundTotal = graph_->createRenderable( + std::format("{}//WeightedSQ//BoundTotal", module_->name()), "Bound F(Q)", "Calculated"); + boundTotal->setColour(StockColours::GreenStockColour); + boundTotal->lineStyle().setStipple(LineStipple::DotStipple); + auto unboundTotal = graph_->createRenderable( + std::format("{}//WeightedSQ//UnboundTotal", module_->name()), "Unbound F(Q)", "Calculated"); + unboundTotal->setColour(StockColours::GreenStockColour); + unboundTotal->lineStyle().setStipple(LineStipple::HalfDashStipple); + + // Add on reference F(Q) data if present + if (referenceFileAndFormat.hasFilename()) + graph_ + ->createRenderable(std::format("{}//ReferenceData", module_->name()), "Reference F(Q)", + "Reference") + ->setColour(StockColours::RedStockColour); + } + else if (ui_.PartialSQButton->isChecked()) + { + targetPartials_ = dissolve_.processingModuleData().valueIf("WeightedSQ", module_->name()); + createPartialSetRenderables("WeightedSQ"); + } + else if (ui_.TotalGRButton->isChecked()) + { + graph_->createRenderable(std::format("{}//WeightedGR//Total", module_->name()), "Calculated", + "Calculated"); + auto repGR = graph_->createRenderable(std::format("{}//RepresentativeTotalGR", module_->name()), + "Via FT", "Calculated"); + repGR->lineStyle().setStipple(LineStipple::HalfDashStipple); + repGR->setColour(StockColours::GreenStockColour); + + // Add on reference G(r) (from FT of F(Q)) if present + if (referenceFileAndFormat.hasFilename()) + graph_ + ->createRenderable(std::format("{}//ReferenceDataFT", module_->name()), + "Reference G(r) (via FT)", "Reference") + ->setColour(StockColours::RedStockColour); + } + else if (ui_.PartialGRButton->isChecked()) + { + targetPartials_ = dissolve_.processingModuleData().valueIf("WeightedGR", module_->name()); + createPartialSetRenderables("WeightedGR"); + } + } + + // Validate renderables if they need it + graph_->validateRenderables(dissolve_.processingModuleData()); + + graph_->postRedisplay(); + ui_.PlotWidget->updateToolbar(); + + refreshing_ = false; +} + +/* + * Widgets / Functions + */ + +void CGNeutronSQModuleWidget::on_TotalFQButton_clicked(bool checked) +{ + if (!checked) + return; + + graph_->view().axes().setTitle(0, "\\it{Q}, \\sym{angstrom}\\sup{-1}"); + graph_->view().axes().setTitle(1, "F(Q)"); + graph_->groupManager().setVerticalShiftAmount(RenderableGroupManager::NoVerticalShift); + + updateControls(ModuleWidget::RecreateRenderablesFlag); +} + +void CGNeutronSQModuleWidget::on_PartialSQButton_clicked(bool checked) +{ + if (!checked) + return; + + graph_->view().axes().setTitle(0, "\\it{Q}, \\sym{angstrom}\\sup{-1}"); + graph_->view().axes().setTitle(1, "S(Q)"); + graph_->groupManager().setVerticalShiftAmount(RenderableGroupManager::TwoVerticalShift); + + updateControls(ModuleWidget::RecreateRenderablesFlag); +} + +void CGNeutronSQModuleWidget::on_TotalGRButton_clicked(bool checked) +{ + if (!checked) + return; + + graph_->view().axes().setTitle(0, "\\it{r}, \\sym{angstrom}"); + graph_->view().axes().setTitle(1, "G(r)"); + graph_->groupManager().setVerticalShiftAmount(RenderableGroupManager::NoVerticalShift); + + updateControls(ModuleWidget::RecreateRenderablesFlag); +} + +void CGNeutronSQModuleWidget::on_PartialGRButton_clicked(bool checked) +{ + if (!checked) + return; + + graph_->view().axes().setTitle(0, "\\it{r}, \\sym{angstrom}"); + graph_->view().axes().setTitle(1, "g(r)"); + graph_->groupManager().setVerticalShiftAmount(RenderableGroupManager::TwoVerticalShift); + + updateControls(ModuleWidget::RecreateRenderablesFlag); +} + +void CGNeutronSQModuleWidget::on_FilterEdit_textChanged(QString text) { updateControls(ModuleWidget::RecreateRenderablesFlag); } + +void CGNeutronSQModuleWidget::on_ClearFilterButton_clicked(bool checked) +{ + ui_.FilterEdit->setText(""); + updateControls(ModuleWidget::RecreateRenderablesFlag); +} diff --git a/src/modules/cgNeutronSQ/process.cpp b/src/modules/cgNeutronSQ/process.cpp new file mode 100644 index 0000000000..89b778b50f --- /dev/null +++ b/src/modules/cgNeutronSQ/process.cpp @@ -0,0 +1,290 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2025 Team Dissolve and contributors + +#include "classes/box.h" +#include "classes/configuration.h" +#include "classes/neutronWeights.h" +#include "classes/species.h" +#include "io/export/data1D.h" +#include "main/dissolve.h" +#include "math/filters.h" +#include "math/ft.h" +#include "module/context.h" +#include "modules/gr/gr.h" +#include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "modules/sq/sq.h" + +// Set target data +void CGNeutronSQModule::setTargets(const std::vector> &configurations, + const std::map> &moduleMap) +{ + auto sqIt = moduleMap.find(ModuleTypes::SQ); + if (sqIt != moduleMap.end()) + sourceSQ_ = dynamic_cast(sqIt->second.front()); +} + +// Run set-up stage +bool CGNeutronSQModule::setUp(ModuleContext &moduleContext, Flags actionSignals) +{ + /* + * Load and set up reference data (if a file/format was given) + */ + if (referenceFQ_.hasFilename() && actionSignals.isSetOrNone(KeywordBase::ReloadExternalData)) + { + // Load the data + Data1D referenceData; + if (!referenceFQ_.importData(referenceData, &moduleContext.processPool())) + return Messenger::error("[SETUP {}] Failed to load reference data '{}'.\n", name_, referenceFQ_.filename()); + + // Get dependent modules + if (!sourceSQ_) + return Messenger::error("[SETUP {}] A source SQ module must be provided.\n", name_); + auto *rdfModule = sourceSQ_->sourceGR(); + if (!rdfModule) + return Messenger::error("[SETUP {}] A source GR module (in the SQ module) must be provided.\n", name_); + + // Normalise reference data to be consistent with the calculated data + if (referenceNormalisedTo_ != normaliseTo_) + { + // We need the neutron weights in order to do the normalisation + NeutronWeights weights; + calculateWeights(rdfModule, weights); + auto factor = 1.0; + + // Set up the multiplication factors + switch (referenceNormalisedTo_) + { + case (StructureFactors::NoNormalisation): + factor = 1.0 / (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation + ? weights.boundCoherentSquareOfAverage() + : weights.boundCoherentAverageOfSquares()); + break; + case (StructureFactors::SquareOfAverageNormalisation): + factor = weights.boundCoherentSquareOfAverage(); + if (normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation) + factor /= weights.boundCoherentAverageOfSquares(); + break; + case (StructureFactors::AverageOfSquaresNormalisation): + factor = weights.boundCoherentAverageOfSquares(); + if (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation) + factor /= weights.boundCoherentSquareOfAverage(); + break; + default: + Messenger::exception("Unhandled StructureFactor::NormalisationType ({}).\n", + StructureFactors::normalisationTypes().keyword(referenceNormalisedTo_)); + } + + // Apply normalisation factors to the data + referenceData *= factor; + } + + // Get Q-range and window function to use for transformation of F(Q) to G(r) + auto ftQMin = referenceFTQMin_.value_or(0.0); + auto ftQMax = referenceFTQMax_.value_or(referenceData.xAxis().back() + 1.0); + if (referenceWindowFunction_ == WindowFunction::Form::None) + Messenger::print("[SETUP {}] No window function will be applied in Fourier transform of reference data to g(r).", + name_); + else + Messenger::print("[SETUP {}] Window function to be applied in Fourier transform of reference data is {}.", name_, + WindowFunction::forms().keyword(referenceWindowFunction_)); + + // Store the reference data in processing + referenceData.setTag(name()); + auto &storedData = moduleContext.dissolve().processingModuleData().realise("ReferenceData", name(), + GenericItem::ProtectedFlag); + storedData = referenceData; + + // Calculate and store the FT of the reference data in processing + referenceData.setTag(name()); + auto &storedDataFT = moduleContext.dissolve().processingModuleData().realise("ReferenceDataFT", name(), + GenericItem::ProtectedFlag); + storedDataFT = referenceData; + Filters::trim(storedDataFT, ftQMin, ftQMax); + auto rho = rdfModule->effectiveDensity(); + if (rho) + Messenger::print( + "[SETUP {}] Effective atomic density used in Fourier transform of reference data is {} atoms/Angstrom3.\n", + name_, rho.value()); + else + Messenger::warn("[SETUP {}] Effective atomic density used in Fourier transform of reference data not yet " + "available, so a default of 0.1 atoms/Angstrom3 used.\n", + name_); + Fourier::sineFT(storedDataFT, 1.0 / (2.0 * PI * PI * rho.value_or(0.1)), referenceFTDeltaR_, referenceFTDeltaR_, 30.0, + WindowFunction(referenceWindowFunction_)); + + // Save data? + if (saveReference_) + { + if (moduleContext.processPool().isMaster()) + { + Data1DExportFileFormat exportFormat(std::format("{}-ReferenceData.q", name())); + if (!exportFormat.exportData(storedData)) + return moduleContext.processPool().decideFalse(); + Data1DExportFileFormat exportFormatFT(std::format("{}-ReferenceData.r", name())); + if (!exportFormatFT.exportData(storedDataFT)) + return moduleContext.processPool().decideFalse(); + moduleContext.processPool().decideTrue(); + } + else if (!moduleContext.processPool().decision()) + return false; + } + } + + return true; +} + +// Run main processing +Module::ExecutionResult CGNeutronSQModule::process(ModuleContext &moduleContext) +{ + /* + * Calculate neutron structure factors from existing S(Q) data + * + * This is a serial routine, with each process constructing its own copy of the data. + * Partial calculation routines called by this routine are parallel. + */ + + if (!sourceSQ_) + { + Messenger::error("A source SQ module must be provided.\n"); + return ExecutionResult::Failed; + } + + const auto *rdfModule = sourceSQ_->sourceGR(); + if (!rdfModule) + { + Messenger::error("A source GR module (in the SQ module) must be provided.\n"); + return ExecutionResult::Failed; + } + + // Print argument/parameter summary + Messenger::print("NeutronSQ: Source unweighted S(Q) will be taken from module '{}'.\n", sourceSQ_->name()); + if (referenceWindowFunction_ == WindowFunction::Form::None) + Messenger::print("No window function will be applied when calculating representative g(r) from S(Q)."); + else + Messenger::print("Window function to be applied when calculating representative g(r) from S(Q) is {}.", + WindowFunction::forms().keyword(referenceWindowFunction_)); + if (normaliseTo_ == StructureFactors::NoNormalisation) + Messenger::print("NeutronSQ: No normalisation will be applied to total F(Q).\n"); + else if (normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation) + Messenger::print("NeutronSQ: Total F(Q) will be normalised to "); + else if (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation) + Messenger::print("NeutronSQ: Total F(Q) will be normalised to **2"); + if (saveSQ_) + Messenger::print("NeutronSQ: Weighted partial S(Q) and total F(Q) will be saved.\n"); + if (saveGR_) + Messenger::print("NeutronSQ: Weighted partial g(r) and total G(r) will be saved.\n"); + if (saveRepresentativeGR_) + Messenger::print("NeutronSQ: Representative G(r) will be saved.\n"); + Messenger::print("\n"); + + /* + * Transform UnweightedSQ from provided SQ data into WeightedSQ. + */ + + // Get unweighted S(Q) from the specified SQMOdule + if (!moduleContext.dissolve().processingModuleData().contains("UnweightedSQ", sourceSQ_->name())) + { + Messenger::error("Couldn't locate unweighted S(Q) data from the SQModule '{}'.\n", sourceSQ_->name()); + return ExecutionResult::Failed; + } + const auto &unweightedSQ = + moduleContext.dissolve().processingModuleData().value("UnweightedSQ", sourceSQ_->name()); + + // Calculate and store weights + auto &weights = moduleContext.dissolve().processingModuleData().realise("FullWeights", name_, + GenericItem::InRestartFileFlag); + calculateWeights(rdfModule, weights); + Messenger::print("Isotopologue and isotope composition:\n\n"); + weights.print(); + + // Does a PartialSet for the weighted S(Q) already exist for this Configuration? + auto [weightedSQ, wSQstatus] = moduleContext.dissolve().processingModuleData().realiseIf( + "WeightedSQ", name_, GenericItem::InRestartFileFlag); + if (wSQstatus == GenericItem::ItemStatus::Created) + weightedSQ.setUpPartials(unweightedSQ.atomTypeMix()); + + std::vector ff; + calculateBeadFormFactor(unweightedSQ.boundPartial(0, 0).xAxis(), ff, weights); + std::vector singleBead; + calculateSingleBead(singleBead, ff, weights); + + + // Calculate weighted S(Q) + calculateWeightedSQ(unweightedSQ, weightedSQ, weights, ff, singleBead, normaliseTo_); + + // Save data if requested + if (saveSQ_ && (!MPIRunMaster(moduleContext.processPool(), weightedSQ.save(name_, "WeightedSQ", "sq", "Q, 1/Angstroms")))) + return ExecutionResult::Failed; + + /* + * Transform UnweightedGR from underlying RDF data into WeightedGR. + */ + + // Get summed unweighted g(r) from the RDFMOdule + if (!moduleContext.dissolve().processingModuleData().contains("UnweightedGR", rdfModule->name())) + { + Messenger::error("Couldn't locate summed unweighted g(r) data.\n"); + return ExecutionResult::Failed; + } + + const auto &unweightedGR = + moduleContext.dissolve().processingModuleData().value("UnweightedGR", rdfModule->name()); + + // Create/retrieve PartialSet for summed weighted g(r) + auto [weightedGR, wGRstatus] = moduleContext.dissolve().processingModuleData().realiseIf( + "WeightedGR", name_, GenericItem::InRestartFileFlag); + if (wGRstatus == GenericItem::ItemStatus::Created) + weightedGR.setUpPartials(unweightedGR.atomTypeMix()); + + // Calculate weighted g(r) + calculateWeightedGR(unweightedGR, weightedGR, weights, normaliseTo_); + + // Save data if requested + if (saveGR_ && (!MPIRunMaster(moduleContext.processPool(), weightedGR.save(name_, "WeightedGR", "gr", "r, Angstroms")))) + return ExecutionResult::Failed; + + // Calculate representative total g(r) from FT of calculated F(Q) + auto &repGR = moduleContext.dissolve().processingModuleData().realise("RepresentativeTotalGR", name_, + GenericItem::InRestartFileFlag); + repGR = weightedSQ.total(); + auto ftQMax = 0.0; + if (referenceFTQMax_) + ftQMax = referenceFTQMax_.value(); + else if (referenceFQ_.hasFilename()) + { + // Take FT max Q limit from reference data + auto &referenceData = moduleContext.dissolve().processingModuleData().realise("ReferenceData", name(), + GenericItem::ProtectedFlag); + ftQMax = referenceData.xAxis().back(); + } + else + ftQMax = weightedSQ.total().xAxis().back(); + Filters::trim(repGR, referenceFTQMin_.value_or(0.0), ftQMax); + auto rMin = weightedGR.total().xAxis().front(); + auto rMax = weightedGR.total().xAxis().back(); + auto rho = rdfModule->effectiveDensity(); + if (!rho) + { + Messenger::error("No effective density available from RDF module '{}'\n", rdfModule->name()); + return ExecutionResult::Failed; + } + + Fourier::sineFT(repGR, 1.0 / (2.0 * PI * PI * *rho), rMin, 0.05, rMax, WindowFunction(referenceWindowFunction_)); + + // Save data if requested + if (saveRepresentativeGR_) + { + if (moduleContext.processPool().isMaster()) + { + Data1DExportFileFormat exportFormat(std::format("{}-weighted-total.gr.broad", name_)); + if (exportFormat.exportData(repGR)) + moduleContext.processPool().decideTrue(); + else + moduleContext.processPool().decideFalse(); + } + else if (!moduleContext.processPool().decision()) + return ExecutionResult::Failed; + } + + return ExecutionResult::Success; +} \ No newline at end of file diff --git a/src/modules/registry.cpp b/src/modules/registry.cpp index 9ec54fa8fe..db0f0d75f6 100644 --- a/src/modules/registry.cpp +++ b/src/modules/registry.cpp @@ -9,6 +9,7 @@ #include "modules/axisAngle/axisAngle.h" #include "modules/benchmark/benchmark.h" #include "modules/bragg/bragg.h" +#include "modules/cgNeutronSQ/cgNeutronSQ.h" #include "modules/clustering/clustering.h" #include "modules/compare/compare.h" #include "modules/dAngle/dAngle.h" @@ -55,6 +56,8 @@ ModuleRegistry::ModuleRegistry() "Checks & Tests"); registerProducer(ModuleTypes::Bragg, "Calculate Bragg intensities", "Correlation Functions"); registerProducer(ModuleTypes::Clustering, "Analyse clustering between sites", "Analysis"); + registerProducer(ModuleTypes::CGNeutronSQ, "Calculate CG neutron-weighted S(Q)", + "Correlation Functions"); registerProducer(ModuleTypes::Compare, "Compare data sets and calculate errors", "Checks & Tests"); registerProducer(ModuleTypes::DAngle, "Calculate distance/angle maps", "Analysis"); registerProducer(ModuleTypes::Energy, "Calculate the total energy of a Configuration", "Forcefield"); diff --git a/src/modules/widgetProducer.cpp b/src/modules/widgetProducer.cpp index b56f65b6d4..2bd255db90 100644 --- a/src/modules/widgetProducer.cpp +++ b/src/modules/widgetProducer.cpp @@ -16,6 +16,8 @@ #include "modules/benchmark/gui/benchmarkWidget.h" #include "modules/bragg/bragg.h" #include "modules/bragg/gui/braggWidget.h" +#include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "modules/cgNeutronSQ/gui/cgNeutronSQWidget.h" #include "modules/clustering/clustering.h" #include "modules/clustering/gui/clusteringWidget.h" #include "modules/compare/compare.h" @@ -75,6 +77,7 @@ ModuleWidgetProducer::ModuleWidgetProducer() registerProducer(); registerProducer(); registerProducer(); + registerProducer(); registerProducer(); registerProducer(); registerProducer(); From 149630bd0a1be18339f4bd94e70655f1475052a0 Mon Sep 17 00:00:00 2001 From: Benjamin Speake Date: Thu, 13 Nov 2025 16:07:07 +0000 Subject: [PATCH 2/2] Fix formatting --- src/modules/cgNeutronSQ/cgNeutronSQ.h | 12 ++++--- src/modules/cgNeutronSQ/functions.cpp | 32 +++++++++++-------- .../gui/cgNeutronSQWidgetFuncs.cpp | 2 +- src/modules/cgNeutronSQ/process.cpp | 5 ++- 4 files changed, 28 insertions(+), 23 deletions(-) diff --git a/src/modules/cgNeutronSQ/cgNeutronSQ.h b/src/modules/cgNeutronSQ/cgNeutronSQ.h index 8b6d3afb59..a4a5ebcfee 100644 --- a/src/modules/cgNeutronSQ/cgNeutronSQ.h +++ b/src/modules/cgNeutronSQ/cgNeutronSQ.h @@ -67,15 +67,17 @@ class CGNeutronSQModule : public Module 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, + bool calculateWeightedSQ(const PartialSet &unweightedsq, PartialSet &weightedsq, NeutronWeights &weights, const std::vector &ff, const std::vector &singleBead, StructureFactors::NormalisationType normalisation); // Calculate neutron weights for relevant Configuration targets void calculateWeights(const GRModule *rdfModule, NeutronWeights &weights) const; - // - bool calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, const NeutronWeights &weights) const; - // - bool calculateSingleBead(std::vector &singleBead, const std::vector &ff, const NeutronWeights &weights) const; + // Calculate the per bead form factor + bool calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, + const NeutronWeights &weights) const; + // Calculate the single bead internal scattering term + bool calculateSingleBead(std::vector &singleBead, const std::vector &ff, + const NeutronWeights &weights) const; /* * Processing diff --git a/src/modules/cgNeutronSQ/functions.cpp b/src/modules/cgNeutronSQ/functions.cpp index a554678b5f..b9417f66a5 100644 --- a/src/modules/cgNeutronSQ/functions.cpp +++ b/src/modules/cgNeutronSQ/functions.cpp @@ -6,12 +6,12 @@ #include "classes/configuration.h" #include "classes/isotopologueSet.h" #include "classes/species.h" -#include "modules/gr/gr.h" #include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "modules/gr/gr.h" // Calculate weighted g(r) from supplied unweighted g(r) and neutron weights bool CGNeutronSQModule::calculateWeightedGR(const PartialSet &unweightedgr, PartialSet &weightedgr, NeutronWeights &weights, - StructureFactors::NormalisationType normalisation) + StructureFactors::NormalisationType normalisation) { int typeI, typeJ; for (typeI = 0; typeI < unweightedgr.nAtomTypes(); ++typeI) @@ -55,7 +55,7 @@ bool CGNeutronSQModule::calculateWeightedGR(const PartialSet &unweightedgr, Part // Calculate weighted S(Q) from supplied unweighted S(Q) and neutron weights bool CGNeutronSQModule::calculateWeightedSQ(const PartialSet &unweightedsq, PartialSet &weightedsq, NeutronWeights &weights, - const std::vector &ff, const std::vector &singleBead, + const std::vector &ff, const std::vector &singleBead, StructureFactors::NormalisationType normalisation) { int typeI, typeJ; @@ -87,7 +87,8 @@ bool CGNeutronSQModule::calculateWeightedSQ(const PartialSet &unweightedsq, Part // Form total structure factor weightedsq.formTotals(false); - for (typeI = 0; typeI < unweightedsq.nAtomTypes(); ++typeI) { + for (typeI = 0; typeI < unweightedsq.nAtomTypes(); ++typeI) + { weightedsq.total() += singleBead[typeI]; } @@ -129,12 +130,13 @@ void CGNeutronSQModule::calculateWeights(const GRModule *rdfModule, NeutronWeigh weights.createFromIsotopologues(exchangeable_); } -// -bool CGNeutronSQModule::calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, const NeutronWeights &weights) const +// Calculate the per bead form factor +bool CGNeutronSQModule::calculateBeadFormFactor(const std::vector &qvals, std::vector &ff, + const NeutronWeights &weights) const { int typeI; int dens{0}; - //std::array sigma{1.63, 1.57}; + // std::array sigma{1.63, 1.57}; std::array sigma{3.0376754356895721, 3.7573162947024161}; double f; ff.clear(); @@ -167,9 +169,11 @@ bool CGNeutronSQModule::calculateBeadFormFactor(const std::vector &qvals return true; } -// -bool CGNeutronSQModule::calculateSingleBead(std::vector &singleBead, const std::vector &ff, const NeutronWeights &weights) const { - +// Calculate the single bead internal scattering term +bool CGNeutronSQModule::calculateSingleBead(std::vector &singleBead, const std::vector &ff, + const NeutronWeights &weights) const +{ + singleBead.clear(); singleBead.reserve(weights.atomTypes().nItems()); std::array b{0.6646, -0.3739}; @@ -191,12 +195,12 @@ bool CGNeutronSQModule::calculateSingleBead(std::vector &singleBead, con innerSum += n[s][atmI] * n[s][atmJ] * b[atmI] * b[atmJ] * (atmI == atmJ ? 1 : 2); } } - innerSum -= ((n[s][0] * b[0] * b[0]) + (n[s][1] * b[1] * b[1])); + innerSum -= ((n[s][0] * b[0] * b[0]) + (n[s][1] * b[1] * b[1])); innerSum *= 0.5; // (*bead).fraction(); singleBead[s].copyArrays(ff[s]); singleBead[s] *= ff[s].values(); - singleBead[s] *= innerSum * (1/9.5); + singleBead[s] *= innerSum * (1 / 9.5); } - - return true; + + return true; } \ No newline at end of file diff --git a/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp index 1016fd4a44..c1ae913651 100644 --- a/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp +++ b/src/modules/cgNeutronSQ/gui/cgNeutronSQWidgetFuncs.cpp @@ -6,8 +6,8 @@ #include "gui/dataViewer.h" #include "gui/render/renderableData1D.h" #include "main/dissolve.h" -#include "modules/cgNeutronSQ/gui/cgNeutronSQWidget.h" #include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "modules/cgNeutronSQ/gui/cgNeutronSQWidget.h" #include "templates/algorithms.h" CGNeutronSQModuleWidget::CGNeutronSQModuleWidget(QWidget *parent, CGNeutronSQModule *module, Dissolve &dissolve) diff --git a/src/modules/cgNeutronSQ/process.cpp b/src/modules/cgNeutronSQ/process.cpp index 89b778b50f..1048dff799 100644 --- a/src/modules/cgNeutronSQ/process.cpp +++ b/src/modules/cgNeutronSQ/process.cpp @@ -10,13 +10,13 @@ #include "math/filters.h" #include "math/ft.h" #include "module/context.h" -#include "modules/gr/gr.h" #include "modules/cgNeutronSQ/cgNeutronSQ.h" +#include "modules/gr/gr.h" #include "modules/sq/sq.h" // Set target data void CGNeutronSQModule::setTargets(const std::vector> &configurations, - const std::map> &moduleMap) + const std::map> &moduleMap) { auto sqIt = moduleMap.find(ModuleTypes::SQ); if (sqIt != moduleMap.end()) @@ -208,7 +208,6 @@ Module::ExecutionResult CGNeutronSQModule::process(ModuleContext &moduleContext) std::vector singleBead; calculateSingleBead(singleBead, ff, weights); - // Calculate weighted S(Q) calculateWeightedSQ(unweightedSQ, weightedSQ, weights, ff, singleBead, normaliseTo_);