feat: Reimplement GR unit test

src/classes/partialSet.cpp

@@ -91,6 +91,10 @@ void PartialSet::initialise(const PartialSet &partialSet)
     triangular_ = partialSet.triangular_;
     rho_ = partialSet.rho_;
 
+    partials_.clear(triangular_);
+    boundPartials_.clear(triangular_);
+    unboundPartials_.clear(triangular_);
+
     // Template data from source PartialSet and set tags
     dissolve::for_each_pair(
         ParallelPolicies::seq, atomTypeFractions(),
@@ -289,11 +293,7 @@ bool PartialSet::save(std::string_view prefix, std::string_view tag, std::string
             std::string filename{std::format("{}-{}-{}-{}.{}", prefix, tag, popI.first->name(), popJ.first->name(), suffix)};
             Messenger::printVerbose("Writing partial file '{}'...\n", filename);
 
-            auto cwd = std::filesystem::current_path();
-            auto path = cwd.parent_path().parent_path() / "tests" / "nodes" / "output" / filename;
-            auto fullPath = path.string();
-
-            parser.openOutput(fullPath, true);
+            parser.openOutput(filename, true);
             if (!parser.isFileGoodForWriting())
                 return Messenger::error("Couldn't open file '{}' for writing.\n", filename);
 

src/classes/species.h

@@ -71,7 +71,7 @@ class Species : public Serialisable<const CoreData &>
     // Add a new atom to the Species, returning its index
     int addAtom(Elements::Element Z, Vector3 r, double q = 0.0, std::shared_ptr<AtomType> atomType = nullptr);
     // Add new atom type to atom types
-    const std::shared_ptr<AtomType> addAtomType(Elements::Element Z);
+    const std::shared_ptr<AtomType> addAtomType(Elements::Element Z, std::string_view name = "");
     // Remove the specified atom from the species
     void removeAtom(int index);
     // Remove set of atom indices

src/classes/species_atomic.cpp

@@ -40,17 +40,15 @@ int Species::addAtom(Elements::Element Z, Vector3 r, double q, std::shared_ptr<A
 }
 
 // Add new atom type to atom types
-const std::shared_ptr<AtomType> Species::addAtomType(Elements::Element Z)
+const std::shared_ptr<AtomType> Species::addAtomType(Elements::Element Z, std::string_view name)
 {
-    auto newAtomType = std::make_shared<AtomType>();
-    atomTypes_.push_back(newAtomType);
-
     // Create a suitable unique name
-    newAtomType->setName(DissolveSys::uniqueName(Elements::symbol(Z), atomTypes_,
-                                                 [&](const auto &at) { return newAtomType == at ? "" : at->name(); }));
+    auto uniqueName = DissolveSys::uniqueName(name == "" ? Elements::symbol(Z) : name, atomTypes_,
+                                              [&](const auto &at) { return at->name(); });
 
-    // Set data
-    newAtomType->setZ(Z);
+    // Create atom type and set data
+    auto newAtomType = std::make_shared<AtomType>(Z, uniqueName);
+    atomTypes_.push_back(newAtomType);
     newAtomType->setIndex(atomTypes_.size() - 1);
 
     return newAtomType;

src/io/import/coordinates.h

@@ -6,8 +6,6 @@
 #include "base/enumOptions.h"
 #include "io/fileAndFormat.h"
 
-// Forward Declarations
-
 // Coordinate Import Formats
 class CoordinateImportFileFormat : public FileAndFormat
 {

src/math/mathFunc.cpp

@@ -78,6 +78,16 @@ double sgn(double a, double signOf) { return signOf >= 0.0 ? fabs(a) : -fabs(a);
 // Return the cyclic permutation of the integer 'i', span 3
 int cp3(int i) { return (i % 3); }
 
+// Return the integer wrapped into the specified range
+int wrap(int i, int lower, int upper)
+{
+    const auto range = upper - lower + 1;
+
+    if (i < 0)
+        i += range * ((lower - i) / range + 1);
+    return lower + (i - lower) % range;
+}
+
 /*
  * Conversion
  */

src/math/mathFunc.h

@@ -39,6 +39,8 @@ int sgn(double x);
 double sgn(double a, double signOf);
 // Return the cyclic permutation of the integer 'i', span 3
 int cp3(int i);
+// Return the integer wrapped into the specified range
+int wrap(int i, int lower, int upper);
 
 /*
  * Conversion

src/nodes/atomicMC/process.cpp

@@ -32,7 +32,7 @@ NodeConstants::ProcessResult AtomicMCNode::process()
     message("\n");
 
     // Prepare for energy calculation, generate kernel
-    auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
+    auto kernel = dissolveGraph()->createEnergyCalculation(targetConfiguration_);
 
     auto nAttempts = 0, nAccepted = 0;
     auto totalDelta = 0.0;

src/nodes/dissolve.cpp

@@ -32,30 +32,32 @@ const DoubleKeyedMap<PairPotential> &DissolveGraph::pairPotentialStore() { retur
  * Functions
  */
 
-// Return maximum distance for tabulated PairPotentials
-const double DissolveGraph::pairPotentialRange() const { return pairPotentialRange_; }
+// Return range of tabulated pari potentials
+double DissolveGraph::pairPotentialRange() const { return pairPotentialRange_; }
 
-// Return energy kernel containing potential map
-std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configuration *cfg)
+// Ensure that the specified Configuration has updated type indexing, cells etc.
+void DissolveGraph::updateIndexingAndCells(Configuration *cfg) const
 {
-    auto atomTypes = cfg->atomTypeVector();
-
     // Update atom type indexing
     cfg->updateTypeIndexing();
 
+    // Regenerate cells in the configuration if necessary
+    cfg->updateCells(pairPotentialRange());
+}
+
+// Create an energy kernel suitable for the supplied Configuration
+std::unique_ptr<EnergyKernel> DissolveGraph::createEnergyCalculation(Configuration *cfg)
+{
+    // Update types and cells in Configuration
+    updateIndexingAndCells(cfg);
+    auto atomTypes = cfg->atomTypeVector();
+
     // Update pair potentials
     dissolve::for_each_pair(ParallelPolicies::seq, atomTypes,
                             [&](int i, const auto &atI, int j, const auto &atJ) { updatePairPotentials(*atI, *atJ); });
 
-    // Generate configuration potential map
-    PotentialMap potentialMap(atomTypes, pairPotentialStore(), pairPotentialRange());
-
-    // Regenerate cells in the configuration if necessary
-    cfg->updateCells(potentialMap.range());
-
-    auto kernel = KernelProducer::energyKernel(cfg, potentialMap);
-
-    return kernel;
+    // Generate and return kernel
+    return KernelProducer::energyKernel(cfg, PotentialMap(atomTypes, pairPotentialStore(), pairPotentialRange()));
 }
 
 // Update pair potential store

src/nodes/dissolve.h

@@ -56,10 +56,12 @@ class DissolveGraph : public Graph
      * Functions
      */
     public:
-    // Return maximum distance for tabulated PairPotentials
-    const double pairPotentialRange() const;
-    // Return suitable energy kernel for the supplied Configuration
-    std::unique_ptr<EnergyKernel> prepareEnergyCalculation(Configuration *cfg);
+    // Return range of tabulated pari potentials
+    double pairPotentialRange() const;
+    // Ensure that the specified Configuration has updated type indexing, cells etc.
+    void updateIndexingAndCells(Configuration *cfg) const;
+    // Create an energy kernel suitable for the supplied Configuration
+    std::unique_ptr<EnergyKernel> createEnergyCalculation(Configuration *cfg);
 
     private:
     // Update pair potential store

src/nodes/edge.cpp

@@ -32,7 +32,7 @@ class EdgeConstructor : public Edge
 std::unique_ptr<Edge> Edge::create(Graph *parent, const EdgeDefinition &definition)
 {
     // Get source node and output
-    auto sourceNode = parent->node(definition.sourceNode);
+    auto sourceNode = parent->findNode(definition.sourceNode);
     if (!sourceNode)
     {
         Messenger::error("Source node '{}' does not exist in the graph.\n", definition.sourceNode);
@@ -54,7 +54,7 @@ std::unique_ptr<Edge> Edge::create(Graph *parent, const EdgeDefinition &definiti
     }
 
     // Get target node and input
-    auto targetNode = parent->node(definition.targetNode);
+    auto targetNode = parent->findNode(definition.targetNode);
     if (!targetNode)
     {
         Messenger::error("Target node '{}' does not exist in the graph.\n", definition.targetNode);

src/nodes/energy/process.cpp

@@ -24,8 +24,8 @@ NodeConstants::ProcessResult EnergyNode::process()
      * This is a serial routine (subroutines called from within are parallel).
      */
 
-    auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
-    auto potentialMap = kernel->potentialMap();
+    auto kernel = dissolveGraph()->createEnergyCalculation(targetConfiguration_);
+    auto &potentialMap = kernel->potentialMap();
 
     // Calculate pair potential energy
     Timer interTimer;

src/nodes/gr/gr.cpp

@@ -29,10 +29,6 @@ GRNode::GRNode(Graph *parentGraph)
     addOption<std::optional<Number>>("Smoothing", "Specifies the degree of smoothing to apply to calculated g(r)", nSmooths_);
     addOption<bool>("Save", "Whether to save partials and total functions to disk", save_);
     addOption<bool>("SaveRaw", "Whether to save raw simulation partial and total functions to disk", saveRaw_);
-    addOption<bool>(
-        "InternalTest",
-        "Perform internal check of calculated partials against a set calculated by a simple unoptimised double-loop",
-        internalTest_);
     addOption<GRNode::PartialsMethod>("Method", "Calculation method for partial radial distribution functions",
                                       partialsMethod_);
 

src/nodes/gr/gr.h

@@ -53,8 +53,6 @@ class GRNode : public Node
     std::optional<Number> averagingLength_{5};
     // Bin width (spacing in r) to use
     Number binWidth_{0.001};
-    // Perform internal check of calculated partials against a set calculated by a simple unoptimised double-loop
-    bool internalTest_{false};
     // Type of broadening to apply to intramolecular g(r)
     Function1DWrapper intraBroadening_{Functions1D::Form::Gaussian, {0.18}};
     // Degree of smoothing to apply
@@ -89,12 +87,6 @@ class GRNode : public Node
     bool calculateRawGR(const double grRange, bool &alreadyUpToDate);
     // Calculate smoothed/broadened partial g(r) from supplied partials
     bool calculateUnweightedGR();
-    // Test supplied PartialSets against each other
-    bool testReferencePartials(const std::vector<const AtomType *> &types, PartialSet &setA, PartialSet &setB,
-                               double testThreshold);
-    // Test calculated partial against supplied reference data
-    bool testReferencePartial(const PartialSet &partials, double testThreshold, const Data1D &testData,
-                              std::string_view typeIorTotal, std::string_view typeJ = "", std::string_view target = "");
 
     /*
      * Processing

src/nodes/gr/helpers.cpp

@@ -8,9 +8,9 @@
 #include "classes/cell.h"
 #include "main/dissolve.h"
 #include "math/combinations.h"
-#include "math/error.h"
 #include "math/filters.h"
 #include "module/group.h"
+#include "nodes/dissolve.h"
 #include "nodes/gr/gr.h"
 #include "templates/algorithms.h"
 #include "templates/combinable.h"
@@ -241,6 +241,7 @@ bool GRNode::calculateGRCells(double grRange, const Array2D<typename std::map<st
                           }
                       });
     }
+
     return true;
 }
 
@@ -319,10 +320,19 @@ bool GRNode::calculateRawGR(const double grRange, bool &alreadyUpToDate)
     else if (partialsMethod_ == PartialsMethod::SimpleMethod)
         calculateGRSimple(fullLUT);
     else if (partialsMethod_ == PartialsMethod::CellsMethod)
+    {
+        dissolveGraph()->updateIndexingAndCells(targetConfiguration_);
         calculateGRCells(grRange, fullLUT);
+    }
     else if (partialsMethod_ == PartialsMethod::AutoMethod)
     {
-        targetConfiguration_->nAtoms() > 10000 ? calculateGRCells(grRange, fullLUT) : calculateGRSimple(fullLUT);
+        if (targetConfiguration_->nAtoms() > 10000)
+        {
+            dissolveGraph()->updateIndexingAndCells(targetConfiguration_);
+            calculateGRCells(grRange, fullLUT);
+        }
+        else
+            calculateGRSimple(fullLUT);
     }
     timer.stop();
     message("Finished calculation of partials ({} elapsed).\n", timer.totalTimeString());
@@ -444,116 +454,3 @@ bool GRNode::calculateUnweightedGR()
 
     return true;
 }
-
-// Test supplied PartialSets against each other
-bool GRNode::testReferencePartials(const std::vector<const AtomType *> &types, PartialSet &setA, PartialSet &setB,
-                                   double testThreshold)
-{
-    for_each_pair_early(types,
-                        [&](int n, const auto *typeI, int m, const auto *typeJ) -> EarlyReturn<bool>
-                        {
-                            DoubleKeyedMapKey key{typeI->name(), typeJ->name()};
-
-                            // Full partial
-                            auto errorReport = Error::percent(setA.partials().get(key), setB.partials().get(key));
-                            message("{}", Error::errorReportString(errorReport));
-                            message("Test reference full partial '{}-{}' has {} error of {:7.3f}{} with calculated data and is "
-                                    "{} (threshold is {:6.3f}%)\n\n",
-                                    typeI->name(), typeJ->name(), Error::errorTypes().keyword(errorReport.errorType),
-                                    errorReport.error, errorReport.errorType == Error::ErrorType::PercentError ? "%" : "",
-                                    errorReport.error <= testThreshold ? "OK" : "NOT OK", testThreshold);
-                            if (errorReport.error > testThreshold)
-                                return false;
-
-                            // Bound partial
-                            errorReport = Error::percent(setA.boundPartials().get(key), setB.boundPartials().get(key));
-                            message("{}", Error::errorReportString(errorReport));
-                            message("Test reference bound partial '{}-{}' has {} error of {:7.3f}{} with calculated data and "
-                                    "is {} (threshold is {:6.3f}%)\n\n",
-                                    typeI->name(), typeJ->name(), Error::errorTypes().keyword(errorReport.errorType),
-                                    errorReport.error, errorReport.errorType == Error::ErrorType::PercentError ? "%" : "",
-                                    errorReport.error <= testThreshold ? "OK" : "NOT OK", testThreshold);
-                            if (errorReport.error > testThreshold)
-                                return false;
-
-                            // Unbound reference
-                            errorReport = Error::percent(setA.unboundPartials().get(key), setB.unboundPartials().get(key));
-                            message("{}", Error::errorReportString(errorReport));
-                            message("Test reference unbound partial '{}-{}' has {} error of {:7.3f}{} with calculated data and "
-                                    "is {} (threshold is {:6.3f}%)\n\n",
-                                    typeI->name(), typeJ->name(), Error::errorTypes().keyword(errorReport.errorType),
-                                    errorReport.error, errorReport.errorType == Error::ErrorType::PercentError ? "%" : "",
-                                    errorReport.error <= testThreshold ? "OK" : "NOT OK", testThreshold);
-                            if (errorReport.error > testThreshold)
-                                return false;
-
-                            return EarlyReturn<bool>::Continue;
-                        });
-
-    // Total reference data supplied?
-    auto errorReport = Error::percent(setA.total(), setB.total());
-    message("{}", Error::errorReportString(errorReport));
-    message("Test reference total has {} error of {:7.3f}{} with calculated data and is {} (threshold is {:6.3f}%)\n\n",
-            Error::errorTypes().keyword(errorReport.errorType), errorReport.error,
-            errorReport.errorType == Error::ErrorType::PercentError ? "%" : "",
-            errorReport.error <= testThreshold ? "OK" : "NOT OK", testThreshold);
-    if (errorReport.error > testThreshold)
-        return false;
-
-    return true;
-}
-
-// Test calculated partial against supplied reference data
-bool GRNode::testReferencePartial(const PartialSet &partials, double testThreshold, const Data1D &testData,
-                                  std::string_view typeIorTotal, std::string_view typeJ, std::string_view target)
-{
-    // We either expect two AtomType names and a target next, or the target 'total'
-    auto testResult = false;
-    if (DissolveSys::sameString(typeIorTotal, "total") && typeJ.empty() && target.empty())
-    {
-        auto errorReport = Error::percent(partials.total(), testData);
-        message("{}", Error::errorReportString(errorReport));
-        testResult = (errorReport.error <= testThreshold);
-        message("Test reference data '{}' has {} error of {:7.3f}{} with calculated data and is {} (threshold is "
-                "{:6.3f}%)\n\n",
-                testData.tag(), Error::errorTypes().keyword(errorReport.errorType), errorReport.error,
-                errorReport.errorType == Error::ErrorType::PercentError ? "%" : "", testResult ? "OK" : "NOT OK",
-                testThreshold);
-    }
-    else
-    {
-        DoubleKeyedMapKey key{typeIorTotal, typeJ};
-        Error::ErrorReport errorReport;
-        if (DissolveSys::sameString(target, "bound"))
-        {
-            errorReport = Error::percent(partials.boundPartials().get(key), testData);
-            message("{}", Error::errorReportString(errorReport));
-        }
-
-        else if (DissolveSys::sameString(target, "unbound"))
-        {
-            errorReport = Error::percent(partials.unboundPartials().get(key), testData);
-            message("{}", Error::errorReportString(errorReport));
-        }
-        else if (DissolveSys::sameString(target, "full"))
-        {
-            errorReport = Error::percent(partials.partials().get(key), testData);
-            message("{}", Error::errorReportString(errorReport));
-        }
-
-        else
-        {
-            error("Unrecognised test data name '{}'.\n", testData.tag());
-            return false;
-        }
-
-        testResult = (errorReport.error <= testThreshold);
-        message("Test reference data '{}' has {} error of {:7.3f}{} with calculated data and is {} (threshold is "
-                "{:6.3f}%)\n\n",
-                testData.tag(), Error::errorTypes().keyword(errorReport.errorType), errorReport.error,
-                errorReport.errorType == Error::ErrorType::PercentError ? "%" : "", testResult ? "OK" : "NOT OK",
-                testThreshold);
-    }
-
-    return testResult;
-}