Skip to content

Commit 1dc24f2

Browse files
trisyoungsTristan Youngs
authored and
rprospero
committed
refactor: Separate histogram calculation from PartialSet (#2216)
Co-authored-by: Tristan Youngs <trisyoungs@googlemail.com>
1 parent d2e669f commit 1dc24f2

17 files changed

Lines changed: 310 additions & 253 deletions

File tree

‎src/classes/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ add_library(
3333
distributor.cpp
3434
empiricalFormula.cpp
3535
fullPairIterator.cpp
36+
histogramSet.cpp
3637
isotopeData.cpp
3738
isotopologue.cpp
3839
isotopologues.cpp
@@ -96,6 +97,7 @@ add_library(
9697
dataSource.h
9798
distributor.h
9899
empiricalFormula.h
100+
histogramSet.h
99101
interactionPotential.h
100102
isotopeData.h
101103
isotopologue.h

‎src/classes/histogramSet.cpp‎

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
// SPDX-License-Identifier: GPL-3.0-or-later
2+
// Copyright (c) 2025 Team Dissolve and contributors
3+
4+
#include "classes/histogramSet.h"
5+
#include "base/lineParser.h"
6+
#include "classes/atomType.h"
7+
#include "items/deserialisers.h"
8+
#include "math/mathFunc.h"
9+
#include "templates/algorithms.h"
10+
11+
HistogramSet::~HistogramSet()
12+
{
13+
fullHistograms_.clear();
14+
boundHistograms_.clear();
15+
unboundHistograms_.clear();
16+
}
17+
18+
/*
19+
* Data
20+
*/
21+
22+
// Initialise histograms
23+
void HistogramSet::initialise(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth)
24+
{
25+
atomTypeMix_ = atomTypeMix;
26+
27+
auto nTypes = atomTypeMix_.nItems();
28+
29+
fullHistograms_.initialise(nTypes, nTypes, half_);
30+
boundHistograms_.initialise(nTypes, nTypes, half_);
31+
unboundHistograms_.initialise(nTypes, nTypes, half_);
32+
33+
dissolve::for_each_pair(
34+
ParallelPolicies::par, nTypes,
35+
[&](int i, int j)
36+
{
37+
fullHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
38+
boundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
39+
unboundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
40+
},
41+
half_);
42+
}
43+
44+
// Clear all histogram data
45+
void HistogramSet::clear()
46+
{
47+
fullHistograms_.clear();
48+
boundHistograms_.clear();
49+
unboundHistograms_.clear();
50+
}
51+
52+
// Zero histogram bins
53+
void HistogramSet::zeroBins()
54+
{
55+
for (auto &histo : fullHistograms_.linearArray())
56+
histo.zeroBins();
57+
for (auto &histo : boundHistograms_.linearArray())
58+
histo.zeroBins();
59+
for (auto &histo : unboundHistograms_.linearArray())
60+
histo.zeroBins();
61+
}
62+
63+
// Return atom types list
64+
const AtomTypeMix &HistogramSet::atomTypeMix() const { return atomTypeMix_; }
65+
66+
// Set new fingerprint
67+
void HistogramSet::setFingerprint(std::string_view fingerprint) { fingerprint_ = fingerprint; }
68+
69+
// Return fingerprint of partials
70+
std::string_view HistogramSet::fingerprint() const { return fingerprint_; }
71+
72+
// Return full histogram specified
73+
Histogram1D &HistogramSet::fullHistogram(int i, int j) { return fullHistograms_[{i, j}]; }
74+
75+
// Return bound histogram specified
76+
Histogram1D &HistogramSet::boundHistogram(int i, int j) { return boundHistograms_[{i, j}]; }
77+
78+
// Return unbound histogram specified
79+
Histogram1D &HistogramSet::unboundHistogram(int i, int j) { return unboundHistograms_[{i, j}]; }
80+
81+
/*
82+
* Manipulation
83+
*/
84+
85+
// Create partials from stored Histogram data
86+
void HistogramSet::formPartials(PartialSet &partials, double boxVolume)
87+
{
88+
dissolve::for_each_pair(
89+
ParallelPolicies::seq, atomTypeMix_,
90+
[&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
91+
{
92+
// Calculate RDFs from histogram data
93+
calculateRDF(partials.partial(n, m), fullHistograms_[{n, m}], boxVolume, at1.population(), at2.population(),
94+
&at1 == &at2 ? 2.0 : 1.0);
95+
calculateRDF(partials.boundPartial(n, m), boundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(),
96+
&at1 == &at2 ? 2.0 : 1.0);
97+
calculateRDF(partials.unboundPartial(n, m), unboundHistograms_[{n, m}], boxVolume, at1.population(),
98+
at2.population(), &at1 == &at2 ? 2.0 : 1.0);
99+
100+
// Set flags for bound partials specifying if they are empty (i.e. there are no
101+
// contributions of that type)
102+
partials.emptyBoundPartial(n, m) = boundHistograms_[{n, m}].nBinned() == 0;
103+
},
104+
half_);
105+
}
106+
107+
// Calculate RDF from supplied Histogram and normalisation data
108+
void HistogramSet::calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres,
109+
int nSurrounding, double multiplier)
110+
{
111+
auto nBins = histogram.nBins();
112+
double delta = histogram.binWidth();
113+
const auto &bins = histogram.bins();
114+
115+
destination.clear();
116+
117+
double shellVolume, factor, r = 0.5 * delta, lowerShellLimit = 0.0, numberDensity = nSurrounding / boxVolume;
118+
for (auto n = 0; n < nBins; ++n)
119+
{
120+
shellVolume = (4.0 / 3.0) * M_PI * (pow(lowerShellLimit + delta, 3.0) - pow(lowerShellLimit, 3.0));
121+
factor = nCentres * (shellVolume * numberDensity);
122+
123+
destination.addPoint(r, bins[n] * (multiplier / factor));
124+
125+
r += delta;
126+
lowerShellLimit += delta;
127+
}
128+
}

‎src/classes/histogramSet.h‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
// SPDX-License-Identifier: GPL-3.0-or-later
2+
// Copyright (c) 2025 Team Dissolve and contributors
3+
4+
#pragma once
5+
6+
#include "classes/atomTypeMix.h"
7+
#include "classes/partialSet.h"
8+
#include "math/histogram1D.h"
9+
#include "templates/array2D.h"
10+
11+
// Set of Histograms
12+
class HistogramSet
13+
{
14+
public:
15+
HistogramSet() = default;
16+
~HistogramSet();
17+
18+
/*
19+
* Data
20+
*/
21+
private:
22+
// AtomTypeMix used to generate matrices
23+
AtomTypeMix atomTypeMix_;
24+
// Fingerprint for these partials (e.g. reflecting Configuration indices at which they were calculated)
25+
std::string fingerprint_;
26+
// Histograms used for calculating full atom-atom partials in r
27+
Array2D<Histogram1D> fullHistograms_;
28+
// Histograms used for calculating bound atom-atom partials in r
29+
Array2D<Histogram1D> boundHistograms_;
30+
// Histograms used for deriving unbound atom-atom partials in r
31+
Array2D<Histogram1D> unboundHistograms_;
32+
// Check for full or half matrix
33+
bool half_{true};
34+
35+
public:
36+
// Set up histograms
37+
void initialise(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth);
38+
// Clear all histogram data
39+
void clear();
40+
// Zero histogram bins
41+
void zeroBins();
42+
// Return atom types mis
43+
const AtomTypeMix &atomTypeMix() const;
44+
// Set new fingerprint
45+
void setFingerprint(std::string_view fingerprint);
46+
// Return fingerprint of partials
47+
std::string_view fingerprint() const;
48+
// Return full histogram specified
49+
Histogram1D &fullHistogram(int i, int j);
50+
// Return bound histogram specified
51+
Histogram1D &boundHistogram(int i, int j);
52+
// Return unbound histogram specified
53+
Histogram1D &unboundHistogram(int i, int j);
54+
55+
/*
56+
* Manipulation
57+
*/
58+
public:
59+
// Form partials from stored Histogram data
60+
void formPartials(PartialSet &partials, double boxVolume);
61+
// Calculate RDF from supplied Histogram and normalisation data
62+
static void calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres,
63+
int nSurrounding, double multiplier);
64+
};

‎src/classes/partialSet.cpp‎

Lines changed: 6 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -19,10 +19,6 @@ PartialSet::PartialSet(const std::map<const Species *, double> &realSpeciesPopul
1919

2020
PartialSet::~PartialSet()
2121
{
22-
fullHistograms_.clear();
23-
boundHistograms_.clear();
24-
unboundHistograms_.clear();
25-
2622
partials_.clear();
2723
boundPartials_.clear();
2824
emptyBoundPartials_.clear();
@@ -33,21 +29,8 @@ PartialSet::~PartialSet()
3329
* Set of Partials
3430
*/
3531

36-
// Set up PartialSet
37-
bool PartialSet::setUp(const AtomTypeMix &atomTypeMix, double rdfRange, double binWidth)
38-
{
39-
// Set up partial arrays
40-
if (!setUpPartials(atomTypeMix, half_))
41-
return false;
42-
43-
// Initialise histograms for g(r) calculation
44-
setUpHistograms(rdfRange, binWidth);
45-
46-
return true;
47-
}
48-
49-
// Set up PartialSet without initialising histogram arrays
50-
bool PartialSet::setUpPartials(const AtomTypeMix &atomTypeMix, bool half)
32+
// Initialise
33+
void PartialSet::initialise(const AtomTypeMix &atomTypeMix, bool half)
5134
{
5235
// Copy type array
5336
atomTypeMix_ = atomTypeMix;
@@ -78,42 +61,11 @@ bool PartialSet::setUpPartials(const AtomTypeMix &atomTypeMix, bool half)
7861
total_.clear();
7962
boundTotal_.clear();
8063
unboundTotal_.clear();
81-
82-
return true;
83-
}
84-
85-
// Set up histogram arrays for g(r) calculation
86-
void PartialSet::setUpHistograms(double rdfRange, double binWidth)
87-
{
88-
auto nTypes = atomTypeMix_.nItems();
89-
90-
fullHistograms_.initialise(nTypes, nTypes, half_);
91-
boundHistograms_.initialise(nTypes, nTypes, half_);
92-
unboundHistograms_.initialise(nTypes, nTypes, half_);
93-
94-
dissolve::for_each_pair(
95-
ParallelPolicies::par, nTypes,
96-
[&](int i, int j)
97-
{
98-
fullHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
99-
boundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
100-
unboundHistograms_[{i, j}].initialise(0.0, rdfRange, binWidth);
101-
},
102-
half_);
10364
}
10465

10566
// Reset partial arrays
10667
void PartialSet::reset()
10768
{
108-
// Zero histogram bins if present
109-
for (auto n = 0; n < fullHistograms_.nRows(); ++n)
110-
for (auto m = n; m < fullHistograms_.nColumns(); ++m)
111-
{
112-
fullHistograms_[{n, m}].zeroBins();
113-
boundHistograms_[{n, m}].zeroBins();
114-
unboundHistograms_[{n, m}].zeroBins();
115-
}
116-
11769
// Zero partials
11870
dissolve::for_each_pair(
11971
ParallelPolicies::par, atomTypeMix_.nItems(),
@@ -144,15 +96,6 @@ void PartialSet::setFingerprint(std::string_view fingerprint) { fingerprint_ = f
14496
// Return fingerprint of partials
14597
std::string_view PartialSet::fingerprint() const { return fingerprint_; }
14698

147-
// Return full histogram specified
148-
Histogram1D &PartialSet::fullHistogram(int i, int j) { return fullHistograms_[{i, j}]; }
149-
150-
// Return bound histogram specified
151-
Histogram1D &PartialSet::boundHistogram(int i, int j) { return boundHistograms_[{i, j}]; }
152-
153-
// Return unbound histogram specified
154-
Histogram1D &PartialSet::unboundHistogram(int i, int j) { return unboundHistograms_[{i, j}]; }
155-
15699
// Return full atom-atom partial specified
157100
Data1D &PartialSet::partial(int i, int j) { return partials_[{i, j}]; }
158101
const Data1D &PartialSet::partial(int i, int j) const { return partials_[{i, j}]; }
@@ -165,6 +108,10 @@ const Data1D &PartialSet::unboundPartial(int i, int j) const { return unboundPar
165108
Data1D &PartialSet::boundPartial(int i, int j) { return boundPartials_[{i, j}]; }
166109
const Data1D &PartialSet::boundPartial(int i, int j) const { return boundPartials_[{i, j}]; }
167110

111+
// Return emptyBound flag
112+
char &PartialSet::emptyBoundPartial(int i, int j) { return emptyBoundPartials_[{i, j}]; }
113+
const char &PartialSet::emptyBoundPartial(int i, int j) const { return emptyBoundPartials_[{i, j}]; }
114+
168115
// Return whether specified bound partial is empty
169116
bool PartialSet::isBoundPartialEmpty(int i, int j) const { return emptyBoundPartials_[{i, j}]; }
170117

@@ -344,28 +291,6 @@ void PartialSet::adjust(double delta)
344291
unboundTotal_ += delta;
345292
}
346293

347-
// Form partials from stored Histogram data
348-
void PartialSet::formPartials(double boxVolume)
349-
{
350-
dissolve::for_each_pair(
351-
ParallelPolicies::seq, atomTypeMix_,
352-
[&](int n, const AtomTypeData &at1, int m, const AtomTypeData &at2)
353-
{
354-
// Calculate RDFs from histogram data
355-
calculateRDF(partials_[{n, m}], fullHistograms_[{n, m}], boxVolume, at1.population(), at2.population(),
356-
&at1 == &at2 ? 2.0 : 1.0);
357-
calculateRDF(boundPartials_[{n, m}], boundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(),
358-
&at1 == &at2 ? 2.0 : 1.0);
359-
calculateRDF(unboundPartials_[{n, m}], unboundHistograms_[{n, m}], boxVolume, at1.population(), at2.population(),
360-
&at1 == &at2 ? 2.0 : 1.0);
361-
362-
// Set flags for bound partials specifying if they are empty (i.e. there are no
363-
// contributions of that type)
364-
emptyBoundPartials_[{n, m}] = boundHistograms_[{n, m}].nBinned() == 0;
365-
},
366-
half_);
367-
}
368-
369294
// Add in partials from source PartialSet to our own
370295
bool PartialSet::addPartials(PartialSet &source, double weighting)
371296
{
@@ -410,29 +335,6 @@ bool PartialSet::addPartials(PartialSet &source, double weighting)
410335
return true;
411336
}
412337

413-
// Calculate RDF from supplied Histogram and normalisation data
414-
void PartialSet::calculateRDF(Data1D &destination, const Histogram1D &histogram, double boxVolume, int nCentres,
415-
int nSurrounding, double multiplier)
416-
{
417-
auto nBins = histogram.nBins();
418-
double delta = histogram.binWidth();
419-
const auto &bins = histogram.bins();
420-
421-
destination.clear();
422-
423-
double shellVolume, factor, r = 0.5 * delta, lowerShellLimit = 0.0, numberDensity = nSurrounding / boxVolume;
424-
for (auto n = 0; n < nBins; ++n)
425-
{
426-
shellVolume = (4.0 / 3.0) * M_PI * (pow(lowerShellLimit + delta, 3.0) - pow(lowerShellLimit, 3.0));
427-
factor = nCentres * (shellVolume * numberDensity);
428-
429-
destination.addPoint(r, bins[n] * (multiplier / factor));
430-
431-
r += delta;
432-
lowerShellLimit += delta;
433-
}
434-
}
435-
436338
/*
437339
* Operators
438340
*/

0 commit comments

Comments
 (0)