1+ // SPDX-License-Identifier: GPL-3.0-or-later
2+ // Copyright (c) 2026 Team Dissolve and contributors
3+
4+ #define _USE_MATH_DEFINES
5+ #include " classes/configuration.h"
6+ #include " classes/species.h"
7+ #include " classes/xRayWeights.h"
8+ #include " io/export/data1D.h"
9+ #include " math/filters.h"
10+ #include " math/ft.h"
11+ #include " nodes/xRaySQ/xRaySQ.h"
12+
13+ bool XRaySQNode::setReferenceData ()
14+ {
15+ // Normalise reference data to be consistent with the calculated data
16+ if (referenceNormalisedTo_ != normaliseTo_)
17+ {
18+ auto bBarSquareOfAverage = weights_.boundCoherentSquareOfAverage (referenceFQ_->xAxis ());
19+ auto bBarAverageOfSquares = weights_.boundCoherentAverageOfSquares (referenceFQ_->xAxis ());
20+ std::vector<double > factors;
21+
22+ // Set up the multiplication factors
23+ switch (referenceNormalisedTo_)
24+ {
25+ case (StructureFactors::NoNormalisation):
26+ factors =
27+ normaliseTo_ == StructureFactors::SquareOfAverageNormalisation ? bBarSquareOfAverage : bBarAverageOfSquares;
28+ std::transform (factors.begin (), factors.end (), factors.begin (), [](const auto factor) { return 1.0 / factor; });
29+ break ;
30+ case (StructureFactors::SquareOfAverageNormalisation):
31+ factors = bBarSquareOfAverage;
32+ if (normaliseTo_ == StructureFactors::AverageOfSquaresNormalisation)
33+ std::transform (factors.begin (), factors.end (), bBarAverageOfSquares.begin (), factors.begin (),
34+ std::divides<>());
35+ break ;
36+ case (StructureFactors::AverageOfSquaresNormalisation):
37+ factors = bBarAverageOfSquares;
38+ if (normaliseTo_ == StructureFactors::SquareOfAverageNormalisation)
39+ std::transform (factors.begin (), factors.end (), bBarSquareOfAverage.begin (), factors.begin (),
40+ std::divides<>());
41+ break ;
42+ default :
43+ Messenger::exception (" Unhandled StructureFactor::NormalisationType ({}).\n "
44+ StructureFactors::normalisationTypes ().keyword (referenceNormalisedTo_));
45+ }
46+
47+ // Apply normalisation factors to the data
48+ std::transform (referenceFQ_->values ().begin (), referenceFQ_->values ().end (), factors.begin (),
49+ referenceFQ_->values ().begin (), std::multiplies<>());
50+ }
51+
52+ // Get Q-range and window function to use for transformation of F(Q) to G(r)
53+ auto ftQMin = referenceFTQMin_.value_or (0.0 );
54+ auto ftQMax = referenceFTQMax_.value_or (referenceFQ_->xAxis ().back () + 1.0 );
55+ if (referenceWindowFunction_ == WindowFunction::Form::None)
56+ Node::message (" [SETUP {}] No window function will be applied in Fourier transform of reference data to g(r)." name ());
57+ else
58+ Node::message (" [SETUP {}] Window function to be applied in Fourier transform of reference data is {}." name (),
59+ WindowFunction::forms ().keyword (referenceWindowFunction_));
60+
61+ referenceGR_ = *referenceFQ_;
62+ Filters::trim (referenceGR_, ftQMin, ftQMax);
63+ auto rho = unweightedGR_->effectiveDensity ();
64+ if (rho)
65+ Node::message (
66+ " [SETUP {}] Effective atomic density used in Fourier transform of reference data is {} atoms/Angstrom3.\n " name (),
67+ rho);
68+ else
69+ Node::message (" [SETUP {}] Effective atomic density used in Fourier transform of reference data not yet "
70+ " available, so a default of 0.1 atoms/Angstrom3 used.\n "
71+ name ());
72+ Fourier::sineFT (referenceGR_, 1.0 / (2.0 * M_PI * M_PI * rho), referenceFTDeltaR_, referenceFTDeltaR_, 30.0 ,
73+ WindowFunction (referenceWindowFunction_));
74+
75+ // Save data?
76+ if (saveReference_)
77+ {
78+ Data1DExportFileFormat exportFormat (std::format (" {}-ReferenceData.q" name ()));
79+ if (!exportFormat.exportData (*referenceFQ_))
80+ return false ;
81+ Data1DExportFileFormat exportFormatFT (std::format (" {}-ReferenceData.r" name ()));
82+ if (!exportFormatFT.exportData (referenceGR_))
83+ return false ;
84+ }
85+
86+ return true ;
87+ }
88+
89+ // Return xRay weights
90+ const XRayWeights &XRaySQNode::weights () const { return weights_; }
91+
92+ // Calculate weighted g(r) from supplied unweighted g(r) and Weights
93+ bool XRaySQNode::calculateWeightedGR (const PartialSet &unweightedgr, PartialSet &weightedgr, const XRayWeights &weights,
94+ StructureFactors::NormalisationType normalisation)
95+ {
96+ dissolve::for_each_pair (ParallelPolicies::seq, unweightedgr.atomTypeFractions (),
97+ [&](int indexI, const auto &popI, int indexJ, const auto &popJ)
98+ {
99+ auto key = DoubleKeyedMapKey{popI.first ->name (), popJ.first ->name ()};
100+
101+ auto weight = weights.weight (popI.first , popJ.first , 0.0 );
102+
103+ // Bound (intramolecular) partial (multiplied by the bound term weight)
104+ weightedgr.boundPartials ().get (key).copyArrays (unweightedgr.boundPartials ().get (key));
105+ weightedgr.boundPartials ().get (key) *= weight;
106+
107+ // Unbound partial (multiplied by the full weight)
108+ weightedgr.unboundPartials ().get (key).copyArrays (unweightedgr.unboundPartials ().get (key));
109+ weightedgr.unboundPartials ().get (key) -= 1.0 ;
110+ weightedgr.unboundPartials ().get (key) *= weight;
111+
112+ // Full partial, summing bound and unbound terms
113+ weightedgr.partials ().get (key).copyArrays (weightedgr.unboundPartials ().get (key));
114+ weightedgr.partials ().get (key) += weightedgr.boundPartials ().get (key);
115+ });
116+
117+ // Form total G(r)
118+ weightedgr.formTotals (false );
119+
120+ // Normalise to Q=0.0 form factor if requested
121+ if (normalisation != StructureFactors::NoNormalisation)
122+ {
123+ auto norm = normalisation == StructureFactors::AverageOfSquaresNormalisation
124+ ? weights.boundCoherentAverageOfSquares (0.0 )
125+ : weights.boundCoherentSquareOfAverage (0.0 );
126+
127+ weightedgr.total () /= norm;
128+ weightedgr.boundTotal () /= norm;
129+ weightedgr.unboundTotal () /= norm;
130+ }
131+
132+ return true ;
133+ }
134+
135+ // Calculate weighted S(Q) from supplied unweighted S(Q) and Weights
136+ bool XRaySQNode::calculateWeightedSQ (const PartialSet &unweightedsq, PartialSet &weightedsq, const XRayWeights &weights,
137+ StructureFactors::NormalisationType normalisation)
138+ {
139+ dissolve::for_each_pair (ParallelPolicies::seq, unweightedsq.atomTypeFractions (),
140+ [&](int indexI, const auto &popI, int indexJ, const auto &popJ)
141+ {
142+ auto key = DoubleKeyedMapKey{popI.first ->name (), popJ.first ->name ()};
143+
144+ // Weight bound and unbound S(Q) and sum into full partial
145+ auto qWeights =
146+ weights.weight (popI.first , popJ.first , unweightedsq.boundPartials ().get (key).xAxis ());
147+
148+ // Bound (intramolecular) and unbound partials
149+ weightedsq.boundPartials ().get (key).copyArrays (unweightedsq.boundPartials ().get (key));
150+ weightedsq.boundPartials ().get (key) *= qWeights;
151+ weightedsq.unboundPartials ().get (key).copyArrays (unweightedsq.unboundPartials ().get (key));
152+ weightedsq.unboundPartials ().get (key) *= qWeights;
153+
154+ // Full partial (sum of bound and unbound terms)
155+ weightedsq.partials ().get (key).copyArrays (weightedsq.unboundPartials ().get (key));
156+ weightedsq.partials ().get (key) += weightedsq.boundPartials ().get (key);
157+ });
158+
159+ // Form total structure factor
160+ weightedsq.formTotals (false );
161+
162+ // Apply normalisation to all totals
163+ if (normalisation != StructureFactors::NoNormalisation)
164+ {
165+ auto bbar = normalisation == StructureFactors::SquareOfAverageNormalisation
166+ ? weights.boundCoherentSquareOfAverage (weightedsq.total ().xAxis ())
167+ : weights.boundCoherentAverageOfSquares (weightedsq.total ().xAxis ());
168+
169+ std::transform (weightedsq.total ().values ().begin (), weightedsq.total ().values ().end (), bbar.begin (),
170+ weightedsq.total ().values ().begin (), std::divides<>());
171+ std::transform (weightedsq.boundTotal ().values ().begin (), weightedsq.boundTotal ().values ().end (), bbar.begin (),
172+ weightedsq.boundTotal ().values ().begin (), std::divides<>());
173+ std::transform (weightedsq.unboundTotal ().values ().begin (), weightedsq.unboundTotal ().values ().end (), bbar.begin (),
174+ weightedsq.unboundTotal ().values ().begin (), std::divides<>());
175+ }
176+
177+ return true ;
178+ }
0 commit comments