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/add.cpp

@@ -5,8 +5,8 @@
 
 AddNode::AddNode(Graph *parentGraph) : Node(parentGraph)
 {
-    addSerialisableInput<Number>("A", "First operand to the addition", a_);
-    addSerialisableInput<Number>("B", "Second operand to the addition", b_);
+    addSerialisableInput<Number>("X", "First operand to the addition", x_);
+    addSerialisableInput<Number>("Y", "Second operand to the addition", y_);
     addOutput<Number>("Result", "The sum of the operands", result_);
 }
 
@@ -18,12 +18,12 @@ AddNode::AddNode(Graph *parentGraph) : Node(parentGraph)
 std::string_view AddNode::type() const { return "Add"; }
 
 // Return short summary of the node's purpose
-std::string_view AddNode::summary() const { return "Performs addition of operands A and B"; }
+std::string_view AddNode::summary() const { return "Performs addition of operands X and Y"; }
 
 // Perform processing
 NodeConstants::ProcessResult AddNode::process()
 {
-    result_ = a_ + b_;
+    result_ = x_ + y_;
 
     return NodeConstants::ProcessResult::Success;
 }

src/nodes/add.h

@@ -26,11 +26,11 @@ class AddNode : public Node
      * Processing & Validity
      */
     private:
-    // Operand A
-    Number a_;
-    // Operand B
-    Number b_;
-    // Result (sum of A and B)
+    // Operand X
+    Number x_;
+    // Operand Y
+    Number y_;
+    // Result (sum of X and Y)
     Number result_;
 
     public:

src/nodes/atomicMC/atomicMC.cpp

@@ -8,6 +8,8 @@ AtomicMCNode::AtomicMCNode(Graph *parentGraph) : Node(parentGraph)
     // Inputs
     addInput<Configuration *>("Configuration", "Set target configuration for the module", targetConfiguration_)
         ->setFlags({ParameterBase::Required, ParameterBase::ClearData});
+    addInput<Number>("Temperature", "Temperature (K)", temperature_)
+        ->setFlags({ParameterBase::Required, ParameterBase::ClearData});
 
     // Options
     addOption<Number>("ShakesPerAtom", "Number of shakes to attempt per atom", nShakesPerAtom_);

src/nodes/atomicMC/atomicMC.h

@@ -26,16 +26,18 @@ class AtomicMCNode : public Node
     private:
     // Target configurations
     Configuration *targetConfiguration_{nullptr};
+    // Temperature (K)
+    Number temperature_{300};
     // Interatomic cutoff distance to use for energy calculation
     std::optional<double> cutoffDistance_;
     // Number of shakes to attempt per atom
     Number nShakesPerAtom_{1};
     // Step size in Angstroms to use in Monte Carlo moves
-    Number stepSize_{0.05};
+    Number stepSize_{0.0001};
     // Maximum allowed value for step size, in Angstroms
     Number stepSizeMax_{1.0};
     // Minimum allowed value for step size, in Angstroms
-    Number stepSizeMin_{0.001};
+    Number stepSizeMin_{0.0001};
     // Target acceptance rate for Monte Carlo moves
     Number targetAcceptanceRate_{0.33};
 

src/nodes/atomicMC/process.cpp

@@ -22,7 +22,7 @@ NodeConstants::ProcessResult AtomicMCNode::process()
 
     // Retrieve control parameters from Configuration
     const auto termScale = 1.0;
-    const auto rRT = 1.0 / (.008314472 * targetConfiguration_->temperature());
+    const auto rRT = 1.0 / (.008314472 * temperature_.asDouble());
 
     // Print argument/parameter summary
     message("Performing {} shake(s) per Atom\n", nShakesPerAtom);
@@ -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/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,10 @@ 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
+    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/iterableGraph.cpp

@@ -5,7 +5,7 @@
 
 IterableGraph::IterableGraph(Graph *parentGraph) : Graph(parentGraph)
 {
-    addOption<Number>("NLoops", "Number of loops (iterations) to perform", nIterations_);
+    addOption<Number>("N", "Number of loops (iterations) to perform", nIterations_);
     loopBacks_ = dynamic_cast<LoopBacksNode *>(addNode(std::make_unique<LoopBacksNode>(this), "LoopBacks"));
 }