diff --git a/src/classes/partialSet.cpp b/src/classes/partialSet.cpp index 5d614a19e3..4a0b8cd072 100644 --- a/src/classes/partialSet.cpp +++ b/src/classes/partialSet.cpp @@ -214,7 +214,7 @@ void PartialSet::formTotals(bool applyConcentrationWeights) } // Sum partials into total for TR -void PartialSet::formTRTotals(bool applyWeights, NeutronWeights weights) +void PartialSet::formTRTotals(NeutronWeights weights) { auto nTypes = atomTypeMix_.nItems(); if (nTypes == 0) @@ -238,7 +238,7 @@ void PartialSet::formTRTotals(bool applyWeights, NeutronWeights weights) [&](int typeI, const AtomTypeData &at1, int typeJ, const AtomTypeData &at2) { // Set weighting factor if requested - auto factor = applyWeights ? at1.fraction() * weights.boundCoherentProduct(typeI, typeJ) : 1.0; + auto factor = at1.fraction() * weights.boundCoherentProduct(typeI, typeJ); // Sum bound term std::transform(boundTotal_.values().begin(), boundTotal_.values().end(), diff --git a/src/classes/partialSet.h b/src/classes/partialSet.h index 460215e12c..99b517423e 100644 --- a/src/classes/partialSet.h +++ b/src/classes/partialSet.h @@ -84,7 +84,7 @@ class PartialSet // Sum partials into totals void formTotals(bool applyConcentrationWeights); // Sum partials into totals for TR - void formTRTotals(bool applyWeights, NeutronWeights weights); + void formTRTotals(NeutronWeights weights); // Return total function Data1D &total(); const Data1D &total() const; diff --git a/src/modules/tr/process.cpp b/src/modules/tr/process.cpp index 9b6b74822a..7382c0e3f4 100644 --- a/src/modules/tr/process.cpp +++ b/src/modules/tr/process.cpp @@ -52,10 +52,7 @@ Module::ExecutionResult TRModule::process(ModuleContext &moduleContext) if (wGRstatus == GenericItem::ItemStatus::Created) weightedTR.setUpPartials(unweightedGR.atomTypeMix(), false); - // Retrieve weights - const auto &weights = moduleData.value("FullWeights", sourceNeutronSQ_->name()); - - /* dissolve::for_each_pair( + dissolve::for_each_pair( ParallelPolicies::par, 0, weightedSQ.nAtomTypes(), [&](int n, int m) { @@ -65,34 +62,18 @@ Module::ExecutionResult TRModule::process(ModuleContext &moduleContext) Fourier::sineFT(weightedGR.boundPartial(n, m), 4.0 * PI * rho.value(), qMin_, qDelta_, qMax_, windowFunction_, qBroadening_); - weightedGR.boundPartial(n, m); + weightedGR.boundPartial(n, m) += 1; Fourier::sineFT(weightedGR.unboundPartial(n, m), 4.0 * PI * rho.value(), qMin_, qDelta_, qMax_, windowFunction_, qBroadening_); weightedGR.unboundPartial(n, m) += 1; Fourier::sineFT(weightedGR.partial(n, m), 4.0 * PI * rho.value(), qMin_, qDelta_, qMax_, windowFunction_, qBroadening_); weightedGR.partial(n, m) += 1; - - // Weight the GR - double weight = weights.weight(n, m); - double intraWeight = weights.intramolecularWeight(n, m); - - // Bound (intramolecular) partial (multiplied by the bound term weight) - weightedGR.boundPartial(n, m) *= intraWeight; - - // Unbound partial (multiplied by the full weight) - weightedGR.unboundPartial(n, m) *= weight; - - // Full partial, summing bound and unbound terms - weightedgr.partial(typeI, typeJ).copyArrays(weightedgr.unboundPartial(typeI, typeJ)); - weightedGR.partial(n, m) += weightedGR.boundPartial(n, m); }, - false); */ - - // Calculate and normalise total to form factor if requested - // weightedGR.formTotals(false); + false); + * / - auto [broadenedTR, bGRstatus] = moduleContext.dissolve().processingModuleData().realiseIf( + auto [broadenedTR, bGRstatus] = moduleContext.dissolve().processingModuleData().realiseIf( "BroadenedTR", name_, GenericItem::InRestartFileFlag); if (bGRstatus == GenericItem::ItemStatus::Created) broadenedTR.setUpPartials(weightedSQ.atomTypeMix(), false); @@ -133,6 +114,9 @@ Module::ExecutionResult TRModule::process(ModuleContext &moduleContext) }, false); + // Retrieve weights + const auto &weights = moduleData.value("FullWeights", sourceNeutronSQ_->name()); + dissolve::for_each_pair( ParallelPolicies::seq, 0, unweightedGR.nAtomTypes(), [&weights, &rho, &unweightedGR, &weightedTR](const auto typeI, const auto typeJ) @@ -166,45 +150,57 @@ Module::ExecutionResult TRModule::process(ModuleContext &moduleContext) }, false); - /* dissolve::for_each_pair( - ParallelPolicies::seq, 0, weightedGR.nAtomTypes(), - [&weights, &rho, &weightedGR, &broadenedTR](const auto typeI, const auto typeJ) + dissolve::for_each_pair( + ParallelPolicies::seq, 0, weightedGR.nAtomTypes(), + [&weights, &rho, &weightedGR, &broadenedTR](const auto typeI, const auto typeJ) + { + double intraWeight = weights.intramolecularWeight(typeI, typeJ); + auto cj = weights.atomTypes()[typeJ].fraction(); + auto factor = 4.0 * PI * rho.value() * cj; + + // Bound (intramolecular) partial (multiplied by the bound term weight) + broadenedTR.boundPartial(typeI, typeJ).copyArrays(weightedGR.boundPartial(typeI, typeJ)); + + for (auto &&[x, y] : + zip(broadenedTR.boundPartial(typeI, typeJ).xAxis(), broadenedTR.boundPartial(typeI, typeJ).values())) + { + double intraWeight = weights.intramolecularWeight(typeI, typeJ); + auto cj = weights.atomTypes()[typeJ].fraction(); + auto factor = 4.0 * PI * rho.value(); + + // Bound (intramolecular) partial (multiplied by the bound term weight) + broadenedTR.boundPartial(typeI, typeJ).copyArrays(weightedGR.boundPartial(typeI, typeJ)); + + for (auto &&[x, y] : + zip(broadenedTR.boundPartial(typeI, typeJ).xAxis(), broadenedTR.boundPartial(typeI, typeJ).values())) { - double intraWeight = weights.intramolecularWeight(typeI, typeJ); - auto cj = weights.atomTypes()[typeJ].fraction(); - auto factor = 4.0 * PI * rho.value(); - - // Bound (intramolecular) partial (multiplied by the bound term weight) - broadenedTR.boundPartial(typeI, typeJ).copyArrays(weightedGR.boundPartial(typeI, typeJ)); - - for (auto &&[x, y] : - zip(broadenedTR.boundPartial(typeI, typeJ).xAxis(), broadenedTR.boundPartial(typeI, typeJ).values())) - { - y *= x * factor; - } - // Unbound partial (multiplied by the full weight) - broadenedTR.unboundPartial(typeI, typeJ).copyArrays(weightedGR.unboundPartial(typeI, typeJ)); - for (auto &&[x, y] : - zip(broadenedTR.unboundPartial(typeI, typeJ).xAxis(), broadenedTR.unboundPartial(typeI, typeJ).values())) - { - y *= x * factor; - } - // Full partial, summing bound and unbound terms - broadenedTR.partial(typeI, typeJ).copyArrays(weightedGR.partial(typeI, typeJ)); - - for (auto &&[x, y] : zip(broadenedTR.partial(typeI, typeJ).xAxis(), broadenedTR.partial(typeI, typeJ).values())) - { - y *= x * factor; - } - }, - false); */ - - // Sum into total - weightedTR.formTRTotals(true, weights); - broadenedTR.formTRTotals(false, weights); - // Save data if requested - if (saveTR_ && (!MPIRunMaster(moduleContext.processPool(), weightedTR.save(name_, "WeightedTR", "tr", "Q, 1/Angstroms")))) - return ExecutionResult::Failed; + y *= x * factor; + } + // Unbound partial (multiplied by the full weight) + broadenedTR.unboundPartial(typeI, typeJ).copyArrays(weightedGR.unboundPartial(typeI, typeJ)); + for (auto &&[x, y] : + zip(broadenedTR.unboundPartial(typeI, typeJ).xAxis(), broadenedTR.unboundPartial(typeI, typeJ).values())) + { + y *= x * factor; + } + // Full partial, summing bound and unbound terms + broadenedTR.partial(typeI, typeJ).copyArrays(weightedGR.partial(typeI, typeJ)); - return ExecutionResult::Success; + for (auto &&[x, y] : zip(broadenedTR.partial(typeI, typeJ).xAxis(), broadenedTR.partial(typeI, typeJ).values())) + { + y *= x * factor; + } + }, + false); + * / + + // Sum into total + weightedTR.formTRTotals(true, weights); + broadenedTR.formTRTotals(false, weights); + // Save data if requested + if (saveTR_ && + (!MPIRunMaster(moduleContext.processPool(), weightedTR.save(name_, "WeightedTR", "tr", "Q, 1/Angstroms")))) + return ExecutionResult::Failed; + + return ExecutionResult::Success; } diff --git a/src/modules/tr/tr.h b/src/modules/tr/tr.h index ff15d7b170..827ce5d9af 100644 --- a/src/modules/tr/tr.h +++ b/src/modules/tr/tr.h @@ -4,7 +4,6 @@ #pragma once #include "classes/partialSet.h" -#include "data/structureFactors.h" #include "math/windowFunction.h" #include "module/module.h" @@ -28,13 +27,11 @@ class TRModule : public Module // Step size in Q for S(Q) calculation double qDelta_{0.05}; // Maximum Q for calculated S(Q) - double qMax_{50.0}; + double qMax_{30.0}; // Minimum Q for calculated S(Q) double qMin_{0.01}; // Window function to use when Fourier-transforming reference S(Q) to g(r)) - WindowFunction::Form windowFunction_{WindowFunction::Form::Lorch0}; - // Normalisation to apply to calculated total F(Q) - StructureFactors::NormalisationType normaliseTo_{StructureFactors::NoNormalisation}; + WindowFunction::Form windowFunction_{WindowFunction::Form::None}; // Broadening function to apply to S(Q) Function1DWrapper qBroadening_; // Source module for calculation