test: Dissolve2 AtomShake unit test

src/classes/atom.cpp

@@ -31,12 +31,6 @@ void Atom::setConfigurationTypeIndex(int id) { configurationTypeIndex_ = id; }
 // Return AtomType index in parent Configuration
 int Atom::configurationTypeIndex() const { return configurationTypeIndex_; }
 
-// Set master AtomType index
-void Atom::setMasterTypeIndex(int id) { masterTypeIndex_ = id; }
-
-// Return master AtomType index
-int Atom::masterTypeIndex() const { return masterTypeIndex_; }
-
 // Return global index of the atom
 int Atom::globalIndex() const
 {

src/classes/atom.h

@@ -24,8 +24,6 @@ class Atom
     Vector3 r_;
     // Assigned AtomType index, local to Configuration (for partial indexing etc.)
     int configurationTypeIndex_{-1};
-    // Assigned master AtomType index (for pair potential indexing)
-    int masterTypeIndex_{-1};
 
     public:
     // Set coordinates
@@ -44,10 +42,6 @@ class Atom
     void setConfigurationTypeIndex(int id);
     // Return AtomType index in parent Configuration
     int configurationTypeIndex() const;
-    // Set master AtomType index
-    void setMasterTypeIndex(int id);
-    // Return master AtomType index
-    int masterTypeIndex() const;
     // Return global index of the atom
     int globalIndex() const;
 

src/classes/cellArray.cpp

@@ -190,11 +190,16 @@ bool CellArray::minimumImageRequired(const Cell &a, const Cell &b) const
 // Generate Cells for Box
 bool CellArray::generate(const Box *box, double cellSize, double pairPotentialRange)
 {
+    // We need to regenerate the cell array only if it is currently empty or the pair potential range has changed
+    if (!cells_.empty() && pairPotentialRangeCreatedAt_ && pairPotentialRangeCreatedAt_.value() == pairPotentialRange)
+        return true;
+
     clear();
 
     const auto minCellsPerSide = 3;
     const auto tolerance = 0.01;
 
+    pairPotentialRangeCreatedAt_ = pairPotentialRange;
     box_ = box;
 
     Messenger::print("Generating cells for box - minimum cells per side is {}, cell size is {}...\n", minCellsPerSide,

src/classes/cellArray.h

@@ -70,6 +70,8 @@ class CellArray
     Matrix3 axes_;
     // Cell array (one-dimensional)
     std::vector<Cell> cells_;
+    // Pair potential range at which the array was last created
+    std::optional<double> pairPotentialRangeCreatedAt_;
     // Box associated with this cell division scheme
     const Box *box_{nullptr};
 

src/classes/configuration_box.cpp

@@ -43,6 +43,9 @@ void Configuration::createBoxAndCells(const Vector3 lengths, const Vector3 angle
 // Create Box definition from axes matrix, and initialise cell array
 void Configuration::createBoxAndCells(const Matrix3 axes, double pairPotentialRange)
 {
+    // Forcibly clear the cell array so we ensure that it is regenerated following the box change
+    cells_.clear();
+
     createBox(axes);
     cells_.generate(box_.get(), requestedCellDivisionLength_, pairPotentialRange);
 }

src/classes/configuration_contents.cpp

@@ -27,8 +27,9 @@ Atom &Configuration::addAtom(const SpeciesAtom *sourceAtom, const std::shared_pt
     // Set the position
     newAtom.setCoordinates(r);
 
-    // Set master index for pair potential lookup
-    newAtom.setMasterTypeIndex(sourceAtom->atomType()->index());
+    // Set configuration type index for pair potential lookup
+    // TODO This can be removed once the Dissolve1 unit tests have been ported over to Dissolve2
+    newAtom.setConfigurationTypeIndex(sourceAtom->atomType()->index());
 
     return newAtom;
 }

src/classes/potentialMap.cpp

@@ -148,7 +148,7 @@ double PotentialMap::energy(const Atom &i, const Atom &j, double r) const
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return pp->energy(r) + (PairPotential::includeCoulombPotential()
                                 ? 0
                                 : pp->analyticCoulombEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r));
@@ -162,7 +162,7 @@ double PotentialMap::energy(const Atom &i, const Atom &j, double r, double elecS
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->energy(r, elecScale, srScale)
                : pp->energy(r) * srScale +
@@ -203,7 +203,7 @@ double PotentialMap::analyticEnergy(const Atom &i, const Atom &j, double r) cons
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being local to the atom
     // types
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->analyticEnergy(r, 1.0, 1.0)
                : pp->analyticEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, 1.0, 1.0);
@@ -216,7 +216,7 @@ double PotentialMap::analyticEnergy(const Atom &i, const Atom &j, double r, doub
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being local to the atom
     // types
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->analyticEnergy(r, elecScale, srScale)
                : pp->analyticEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, elecScale, srScale);
@@ -227,7 +227,7 @@ double PotentialMap::force(const Atom &i, const Atom &j, double r) const
 {
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->force(r)
                : pp->force(r) + pp->analyticCoulombForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r);
@@ -238,7 +238,7 @@ double PotentialMap::force(const Atom &i, const Atom &j, double r, double elecSc
 {
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->force(r, elecScale, srScale)
                : pp->force(r) * srScale +
@@ -280,7 +280,7 @@ double PotentialMap::analyticForce(const Atom &i, const Atom &j, double r) const
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->analyticForce(r, 1.0, 1.0)
                : pp->analyticForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, 1.0, 1.0);
@@ -294,7 +294,7 @@ double PotentialMap::analyticForce(const Atom &i, const Atom &j, double r, doubl
 
     // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the
     // interpolated potential
-    auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}];
+    auto *pp = potentialMatrix_[{i.configurationTypeIndex(), j.configurationTypeIndex()}];
     return PairPotential::includeCoulombPotential()
                ? pp->analyticForce(r, elecScale, srScale)
                : pp->analyticForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, elecScale, srScale);

src/classes/species_atomic.cpp

@@ -349,12 +349,6 @@ double Species::totalCharge(bool useAtomTypes) const
 // Apply random noise to atoms
 void Species::randomiseCoordinates(double maxDisplacement)
 {
-    Vector3 deltaR;
-
     for (auto &i : atoms_)
-    {
-        deltaR.randomUnit();
-        deltaR *= maxDisplacement;
-        i.translateCoordinates(deltaR);
-    }
+        i.translateCoordinates(Vector3::randomUnit() * maxDisplacement);
 }

src/math/vector3.cpp

@@ -405,14 +405,13 @@ void Vector3::orthogonalise(const Vector3 &reference1, const Vector3 &reference2
 // Prints the contents of the vector
 void Vector3::print() const { std::cout << std::format("{} {} {}", x, y, z) << std::endl; }
 
-// Generate random unit vector
-void Vector3::randomUnit()
-{
-    // Generates a random unit vector
-    x = DissolveMath::random() - 0.5;
-    y = DissolveMath::random() - 0.5;
-    z = DissolveMath::random() - 0.5;
-    normalise();
+// Generate random unit vector (on a unit sphere)
+Vector3 Vector3::randomUnit()
+{
+    auto y = 2.0 * (DissolveMath::random() - 0.5);
+    auto r = sqrt(1 - y * y);
+    auto lambda = 2.0 * M_PI * (DissolveMath::random() - 0.5);
+    return {r * sin(lambda), y, r * cos(lambda)};
 }
 
 // Convert spherical who,phi,theta coordinates into cartesian x,y,z

src/math/vector3.h

@@ -113,7 +113,7 @@ class Vector3 : public Serialisable<>
     // Prints the contents of the vector
     void print() const;
     // Generate random unit vector
-    void randomUnit();
+    static Vector3 randomUnit();
     // Convert spherical who,phi,theta coordinates into cartesian x,y,z
     void toCartesian();
     // Convert cartesian x,y,z coordinates into spherical (rho,phi/zenith,theta/azimuthal)

src/nodes/atomicMC/process.cpp

@@ -35,10 +35,7 @@ NodeConstants::ProcessResult AtomicMCNode::process()
     auto kernel = dissolveGraph()->prepareEnergyCalculation(targetConfiguration_);
 
     auto nAttempts = 0, nAccepted = 0;
-    bool accept;
-    double currentEnergy, currentIntraEnergy, newEnergy, newIntraEnergy, delta, totalDelta = 0.0;
-    Vector3 rDelta;
-    EnergyResult er;
+    auto totalDelta = 0.0;
 
     Timer timer;
     // Loop over target Molecules
@@ -52,36 +49,32 @@ NodeConstants::ProcessResult AtomicMCNode::process()
         for (const auto &i : mol->atoms())
         {
             // Calculate reference energies for the Atom
-            er = kernel->totalEnergy(*i);
-            currentEnergy = er.totalUnbound();
-            currentIntraEnergy = er.geometry() * termScale;
+            auto eCurrent = kernel->totalEnergy(*i);
 
             // Loop over number of shakes per Atom
             for (auto n = 0; n < nShakesPerAtom; ++n)
             {
                 auto moveInitialPos = i->r();
 
-                // Create a random translation vector
-                rDelta.set(DissolveMath::randomPlusMinusOne() * stepSize, DissolveMath::randomPlusMinusOne() * stepSize,
-                           DissolveMath::randomPlusMinusOne() * stepSize);
-
-                // Translate Atom and update its Cell position
-                i->translateCoordinates(rDelta);
+                // Translate Atom randomly according to the stepsize and update its Cell position
+                i->translateCoordinates(Vector3::randomUnit() * stepSize);
                 targetConfiguration_->updateAtomLocation(i);
 
                 // Calculate new energy
-                er = kernel->totalEnergy(*i);
-                newEnergy = er.totalUnbound();
-                newIntraEnergy = er.geometry() * termScale;
+                auto eNew = kernel->totalEnergy(*i);
 
                 // Trial the transformed Atom position
-                delta = (newEnergy + newIntraEnergy) - (currentEnergy + currentIntraEnergy);
-                accept = delta < 0 ? true : (DissolveMath::random() < exp(-delta * rRT));
+                auto delta = (eNew.totalUnbound() + eNew.geometry() * termScale) -
+                             (eCurrent.totalUnbound() + eCurrent.geometry() * termScale);
+                auto accept = delta < 0 ? true : (DissolveMath::random() < exp(-delta * rRT));
 
-                // Increase attempt counters
                 if (accept)
                 {
+                    // Store incremental total energy and new reference energy
                     totalDelta += delta;
+                    eCurrent = eNew;
+
+                    // Increase attempt counter
                     ++nAccepted;
                 }
                 else

src/nodes/configuration.cpp

@@ -5,11 +5,16 @@
 
 ConfigurationNode::ConfigurationNode(Graph *parentGraph) : Node(parentGraph)
 {
+    addOption<Number>("Temperature", "Configuration temperature (K)", temperature_);
     addPointerOutput<Configuration>("Configuration", "Configuration object", configuration_);
 }
 
 std::string_view ConfigurationNode::type() const { return "Configuration"; }
 
 std::string_view ConfigurationNode::summary() const { return "Produce an empty atomic configuration."; }
 
-NodeConstants::ProcessResult ConfigurationNode::process() { return NodeConstants::ProcessResult::Unchanged; }
+NodeConstants::ProcessResult ConfigurationNode::process()
+{
+    configuration_.setTemperature(temperature_.asInteger());
+    return NodeConstants::ProcessResult::Unchanged;
+}

src/nodes/configuration.h

@@ -26,6 +26,8 @@ class ConfigurationNode : public Node
     private:
     // Configuration object
     Configuration configuration_;
+    // Configuration temperature
+    Number temperature_;
 
     /*
      * Processing

src/nodes/dissolve.cpp

@@ -36,7 +36,7 @@ const DoubleKeyedMap<PairPotential> &DissolveGraph::pairPotentialStore() { retur
 const double DissolveGraph::pairPotentialRange() const { return pairPotentialRange_; }
 
 // Return energy kernel containing potential map
-std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configuration *cfg, std::optional<double> energyCutoff)
+std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configuration *cfg)
 {
     auto atomTypes = cfg->atomTypeVector();
 
@@ -50,12 +50,11 @@ std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configurat
     // Generate configuration potential map
     PotentialMap potentialMap(atomTypes, pairPotentialStore(), pairPotentialRange());
 
-    // Regenerate cells
-    cfg->cells().generate(cfg->box(), cfg->requestedCellDivisionLength(), potentialMap.range());
+    // Regenerate cells in the configuration if necessary
+    // TODO
+    cfg->updateCells(potentialMap.range());
 
-    auto kernel = KernelProducer::energyKernel(cfg, potentialMap, energyCutoff);
-
-    cfg->updateCells(kernel.get()->potentialMap().range());
+    auto kernel = KernelProducer::energyKernel(cfg, potentialMap);
 
     return kernel;
 }

src/nodes/dissolve.h

@@ -58,8 +58,8 @@ class DissolveGraph : public Graph
     public:
     // Return maximum distance for tabulated PairPotentials
     const double pairPotentialRange() const;
-    // Return energy kernel containing potential map
-    std::unique_ptr<EnergyKernel> prepareEnergyCalculation(Configuration *cfg, std::optional<double> energyCutoff = {});
+    // Return suitable energy kernel for the supplied Configuration
+    std::unique_ptr<EnergyKernel> prepareEnergyCalculation(Configuration *cfg);
 
     private:
     // Update pair potential store

src/nodes/loopGraph.cpp

@@ -5,7 +5,7 @@
 
 LoopGraph::LoopGraph(Graph *parentGraph) : Graph(parentGraph)
 {
-    addOption<Number>("NLoops", "Number of loops (iterations) to perform", nLoops_);
+    addOption<Number>("N", "Number of loops (iterations) to perform", nLoops_);
     loopBacks_ = dynamic_cast<LoopBacksNode *>(addNode(std::make_unique<LoopBacksNode>(this), "LoopBacks"));
 }
 
@@ -149,5 +149,5 @@ NodeConstants::ProcessResult LoopGraph::process()
 
     resetLoopCounter();
 
-    return Graph::process();
+    return NodeConstants::ProcessResult::Success;
 }

src/nodes/node.h

@@ -241,6 +241,15 @@ class Node : public Serialisable<>
     bool ownParameter(std::shared_ptr<ParameterBase> &parameter, bool isOutput = false);
     // Return named input parameter if it exists
     std::shared_ptr<ParameterBase> findInput(std::string_view inputName) const;
+    template <class T> bool setInput(std::string_view name, const T &value)
+    {
+        auto i = findInput(name);
+        if (!i)
+            return false;
+
+        i->set<T>(value);
+        return true;
+    }
     // Return input parameters
     NodeParameterMap &inputs();
     // Get input parameter value

src/nodes/temperature.cpp

@@ -0,0 +1,28 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2025 Team Dissolve and contributors
+
+#include "nodes/temperature.h"
+
+TemperatureNode::TemperatureNode(Graph *parentGraph) : Node(parentGraph)
+{
+    addOption<Number>("Temperature", "Temperature (K)", temperature_);
+    addInput<Configuration *>("Configuration", "Configuration to be modified", targetConfiguration_);
+    addOutput<Configuration *>("Configuration", "Modified configuration", targetConfiguration_);
+}
+
+/*
+ * Definition (Virtuals)
+ */
+
+// Return type of the node
+std::string_view TemperatureNode::type() const { return "Temperature"; }
+
+// Return short summary of the node's purpose
+std::string_view TemperatureNode::summary() const { return "Sets a new temperature value for a given target configuration"; }
+
+// Perform processing
+NodeConstants::ProcessResult TemperatureNode::process()
+{
+    targetConfiguration_->setTemperature(temperature_.asDouble());
+    return NodeConstants::ProcessResult::Success;
+}