initialise potential map from atom type pairs

Author: RobBuchananCompPhys

src/classes/potentialMap.cpp

@@ -103,6 +103,49 @@ bool PotentialMap::initialise(const std::vector<const AtomType*>& masterAtomType
     return true;
 }
 
+bool PotentialMap::initialise(const std::vector<const AtomType *> &masterAtomTypes,
+                              const DoubleKeyedMap<PairPotential*> &pairPotentials, double pairPotentialRange)
+{
+    // Clear old data first
+    clear();
+
+    // Create PairPotential matrix
+    nTypes_ = masterAtomTypes.size();
+    potentialMatrix_.initialise(nTypes_, nTypes_);
+
+    dissolve::for_each_pair(
+        ParallelPolicies::seq, masterAtomTypes,
+        [&](int i, const auto &atI, int j, const auto &atJ)
+        {
+            if (i == -1)
+                return Messenger::error("Couldn't find AtomType '{}' in typeIndex.\n", atI->name());
+            if (j == -1)
+                return Messenger::error("Couldn't find AtomType '{}' in typeIndex.\n", atJ->name());
+
+            auto pp = pairPotentials.get({atI->name(), atJ->name()});
+
+            // Store PairPotential pointer
+            if (i == j)
+            {
+                Messenger::print("Linking self-interaction PairPotential for '{}' (index {},{} in matrix).\n", atI->name(), i,
+                                 j);
+                potentialMatrix_[{i, j}] = pp;
+            }
+            else
+            {
+                Messenger::print("Linking PairPotential between '{}' and '{}' (indices {},{} and {},{} in matrix).\n",
+                                 atI->name(), atJ->name(), i, j, j, i);
+                potentialMatrix_[{i, j}] = pp;
+                potentialMatrix_[{j, i}] = pp;
+            }
+        });
+
+    // Store potential range
+    range_ = pairPotentialRange;
+
+    return true;
+}
+
 // Return PairPotential range
 double PotentialMap::range() const { return range_; }
 

src/classes/potentialMap.h

@@ -5,6 +5,7 @@
 
 #include "classes/pairPotential.h"
 #include "templates/array2D.h"
+#include "templates/doubleKeyedMap.h"
 
 // Forward Declarations
 class Atom;
@@ -37,6 +38,8 @@ class PotentialMap
                     const std::vector<PairPotential::Definition> &pairPotentials, double pairPotentialRange);
     bool initialise(const std::vector<const AtomType *> &atomTypes,
                     const std::vector<PairPotential::Definition> &pairPotentials, double pairPotentialRange);
+    bool initialise(const std::vector<const AtomType *> &atomTypes,
+                    const DoubleKeyedMap<PairPotential*> &pairPotentials, double pairPotentialRange);
     // Return PairPotential range
     double range() const;
 

src/nodes/dissolve.cpp

@@ -22,19 +22,15 @@ std::string_view DissolveGraph::summary() const { return "Parent node of all sim
 // Return dissolve
 Dissolve &DissolveGraph::dissolve() const { return dissolve_; }
 
+// Return pair potential store
+const DoubleKeyedMap<PairPotential*> &DissolveGraph::pairPotentialStore() { return pairPotentialStore_; }
+
 /*
  * Functions
  */
 
 // Return maximum distance for tabulated PairPotentials
-double DissolveGraph::pairPotentialRange() { return pairPotentialRange_; }
-
-// Return first PairPotential in list
-const std::vector<PairPotential::Definition> &DissolveGraph::pairPotentials()
-{
-    auto values = std::views::values(pairPotentialStore_.map());
-    return {values.begin(), values.end()};
-}
+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)
@@ -44,9 +40,7 @@ std::unique_ptr<EnergyKernel> DissolveGraph::prepareEnergyCalculation(Configurat
 
 	// Generate configuration potential map
     PotentialMap potentialMap;
-    auto atomTypeKeys = std::views::keys(cfg->atomTypeIndexMap());
-    std::vector<const AtomType *> atomTypes{atomTypeKeys.begin(), atomTypeKeys.end()};
-    potentialMap.initialise(atomTypes, pairPotentials(), pairPotentialRange());
+    potentialMap.initialise(cfg->atomTypeVector(), pairPotentialStore(), pairPotentialRange());
 
 	// Regenerate cells
     cfg->cells().generate(cfg->box(), cfg->requestedCellDivisionLength(), potentialMap.range());
@@ -61,7 +55,9 @@ void DissolveGraph::updatePairPotentials(const AtomType &i, const AtomType &j)
     if (pairPotentialStore_.contains(nameI, nameJ))
         return;
 
+    auto interactionPotential = ShortRangeFunctions::combine(i.interactionPotential(), j.interactionPotential());
     pairPotentialStore_.set(
         nameI, nameJ,
-        {std::make_shared<AtomType>(i), std::make_shared<AtomType>(j), std::make_unique<PairPotential>(nameI, nameJ)});
+        interactionPotential.has_value() ? new PairPotential(nameI, nameJ, *interactionPotential) : new PairPotential(nameI, nameJ)
+    );
 }

src/nodes/dissolve.h

@@ -36,22 +36,22 @@ class DissolveGraph : public Graph
     // Dissolve reference
     Dissolve &dissolve_;
     // Pair potential store
-    DoubleKeyedMap<PairPotential::Definition> pairPotentialStore_{true};
+    DoubleKeyedMap<PairPotential*> pairPotentialStore_{true};
     // Pair potential range
     double pairPotentialRange_{0.0};
 
     public:
     // Return dissolve
-    Dissolve &dissolve() const override;
+    Dissolve &dissolve() const;
+    // Return pair potential store
+    const DoubleKeyedMap<PairPotential*> &pairPotentialStore();
 
     /*
      * Functions
      */
     public:
     // Return maximum distance for tabulated PairPotentials
-    double pairPotentialRange();
-    // Return first PairPotential in list
-    const std::vector<PairPotential::Definition> &pairPotentials();
+    const double pairPotentialRange() const;
     // Return energy kernel containing potential map
     std::unique_ptr<EnergyKernel> prepareEnergyCalculation(Configuration *cfg, std::optional<double> energyCutoff = {});
     // Update pair potential store