Skip to content

Commit

Permalink
Revert "wightedTR equation corrected"
Browse files Browse the repository at this point in the history 
This reverts commit 1192907.
  • Loading branch information
@Danbr4d
Danbr4d committed Jan 16, 2025
1 parent ed2c730 commit 1502107
Show file tree
Hide file tree
Showing 4 changed files with 64 additions and 71 deletions.
4 changes: 2 additions & 2 deletions src/classes/partialSet.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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(),
Expand Down
2 changes: 1 addition & 1 deletion src/classes/partialSet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
122 changes: 59 additions & 63 deletions src/modules/tr/process.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<NeutronWeights>("FullWeights", sourceNeutronSQ_->name());

/* dissolve::for_each_pair(
dissolve::for_each_pair(
ParallelPolicies::par, 0, weightedSQ.nAtomTypes(),
[&](int n, int m)
{
Expand All @@ -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<PartialSet>(
auto [broadenedTR, bGRstatus] = moduleContext.dissolve().processingModuleData().realiseIf<PartialSet>(
"BroadenedTR", name_, GenericItem::InRestartFileFlag);
if (bGRstatus == GenericItem::ItemStatus::Created)
broadenedTR.setUpPartials(weightedSQ.atomTypeMix(), false);
Expand Down Expand Up @@ -133,6 +114,9 @@ Module::ExecutionResult TRModule::process(ModuleContext &moduleContext)
},
false);

// Retrieve weights
const auto &weights = moduleData.value<NeutronWeights>("FullWeights", sourceNeutronSQ_->name());

dissolve::for_each_pair(
ParallelPolicies::seq, 0, unweightedGR.nAtomTypes(),
[&weights, &rho, &unweightedGR, &weightedTR](const auto typeI, const auto typeJ)
Expand Down Expand Up @@ -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;
}
7 changes: 2 additions & 5 deletions src/modules/tr/tr.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
#pragma once

#include "classes/partialSet.h"
#include "data/structureFactors.h"
#include "math/windowFunction.h"
#include "module/module.h"

Expand All @@ -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
Expand Down

0 comments on commit 1502107

Please sign in to comment.