refactor: Rework Function1D, add derivatives (#1876)

disorderedmaterials · May 14, 2024 · d572e8c · d572e8c
1 parent 420d363
commit d572e8c
Show file tree

Hide file tree

Showing 38 changed files with 387 additions and 337 deletions.
diff --git a/src/classes/coreData.cpp b/src/classes/coreData.cpp
@@ -128,7 +128,7 @@ void CoreData::clearAtomTypes() { atomTypes_.clear(); }
 // Create new pair potential override
 PairPotentialOverride *CoreData::addPairPotentialOverride(std::string_view matchI, std::string_view matchJ,
                                                           PairPotentialOverride::PairPotentialOverrideType overrideType,
-                                                          const InteractionPotential<ShortRangeFunctions> &potential)
+                                                          const InteractionPotential<Functions1D> &potential)
 {
     auto &pp =
         pairPotentialOverrides_.emplace_back(std::make_unique<PairPotentialOverride>(matchI, matchJ, overrideType, potential));

diff --git a/src/classes/coreData.h b/src/classes/coreData.h
@@ -77,8 +77,7 @@ class CoreData
     PairPotentialOverride *addPairPotentialOverride(
         std::string_view matchI = "", std::string_view matchJ = "",
         PairPotentialOverride::PairPotentialOverrideType overrideType = PairPotentialOverride::PairPotentialOverrideType::Off,
-        const InteractionPotential<ShortRangeFunctions> &potential =
-            InteractionPotential<ShortRangeFunctions>(ShortRangeFunctions::Form::None));
+        const InteractionPotential<Functions1D> &potential = {});
     // Return defined overrides for PairPotentials
     std::vector<std::unique_ptr<PairPotentialOverride>> &pairPotentialOverrides();
     const std::vector<std::unique_ptr<PairPotentialOverride>> &pairPotentialOverrides() const;

diff --git a/src/classes/pairPotential.cpp b/src/classes/pairPotential.cpp
@@ -15,22 +15,22 @@ PairPotential::ShortRangeTruncationScheme PairPotential::shortRangeTruncationSch
     PairPotential::ShiftedShortRangeTruncation;
 
 PairPotential::PairPotential()
-    : interactionPotential_(ShortRangeFunctions::Form::None), uFullInterpolation_(uFull_), dUFullInterpolation_(dUFull_)
+    : interactionPotential_(Functions1D::Form::None), uFullInterpolation_(uFull_), dUFullInterpolation_(dUFull_)
 {
 }
 
 PairPotential::PairPotential(const std::shared_ptr<AtomType> &typeI, const std::shared_ptr<AtomType> &typeJ,
                              bool includeCharges)
-    : interactionPotential_(ShortRangeFunctions::Form::None), uFullInterpolation_(uFull_), dUFullInterpolation_(dUFull_)
+    : interactionPotential_(Functions1D::Form::None), uFullInterpolation_(uFull_), dUFullInterpolation_(dUFull_)
 {
     setUp(typeI, typeJ, includeCharges);
 }
 
-PairPotential::PairPotential(std::string_view nameI, std::string_view nameJ,
-                             const InteractionPotential<ShortRangeFunctions> &potential)
+PairPotential::PairPotential(std::string_view nameI, std::string_view nameJ, const InteractionPotential<Functions1D> &potential)
     : includeAtomTypeCharges_(false), nameI_(nameI), nameJ_(nameJ), interactionPotential_(potential),
       uFullInterpolation_(uFull_), dUFullInterpolation_(dUFull_)
 {
+    potentialFunction_.setFormAndParameters(interactionPotential_.form(), interactionPotential_.parameters());
 }
 
 // Return enum option info for CoulombTruncationScheme
@@ -91,7 +91,7 @@ bool PairPotential::setUp(const std::shared_ptr<AtomType> &typeI, const std::sha
         throw(std::runtime_error("Invalid AtomType pointer (typeJ) given to PairPotential::setUp().\n"));
 
     includeAtomTypeCharges_ = includeCharges;
-    interactionPotential_.setFormAndParameters(ShortRangeFunctions::Form::None, "");
+    interactionPotential_.setFormAndParameters(Functions1D::Form::None, "");
 
     nameI_ = typeI->name();
     nameJ_ = typeJ->name();
@@ -105,7 +105,10 @@ bool PairPotential::setUp(const std::shared_ptr<AtomType> &typeI, const std::sha
         return Messenger::error("Can't set parameters for PairPotential since atom type {} has no valid short range form.\n",
                                 typeJ->name());
 
+    // Combine the atom type parameters into potential function parameters
     interactionPotential_ = ShortRangeFunctions::combine(typeI->interactionPotential(), typeJ->interactionPotential());
+    potentialFunction_.setFormAndParameters(interactionPotential_.form(), interactionPotential_.parameters());
+
     if (!interactionPotential_.hasValidForm())
         return false;
 
@@ -135,8 +138,8 @@ std::string_view PairPotential::nameI() const { return nameI_; }
 std::string_view PairPotential::nameJ() const { return nameJ_; };
 
 // Return interaction potential
-InteractionPotential<ShortRangeFunctions> &PairPotential::interactionPotential() { return interactionPotential_; }
-const InteractionPotential<ShortRangeFunctions> &PairPotential::interactionPotential() const { return interactionPotential_; }
+InteractionPotential<Functions1D> &PairPotential::interactionPotential() { return interactionPotential_; }
+const InteractionPotential<Functions1D> &PairPotential::interactionPotential() const { return interactionPotential_; }
 
 // Set charge I
 void PairPotential::setChargeI(double value) { chargeI_ = value; }
@@ -157,32 +160,8 @@ double PairPotential::chargeJ() const { return chargeJ_; }
 // Return analytic short range potential energy
 double PairPotential::analyticShortRangeEnergy(double r, PairPotential::ShortRangeTruncationScheme truncation) const
 {
-    auto &params = interactionPotential_.parameters();
-
-    auto energy = 0.0;
-    switch (interactionPotential_.form())
-    {
-        case (ShortRangeFunctions::Form::None):
-            break;
-        case (ShortRangeFunctions::Form::LennardJones):
-        case (ShortRangeFunctions::Form::LennardJonesGeometric):
-        {
-            /*
-             * Standard Lennard-Jones potential
-             * Parameter 0 = Epsilon
-             * Parameter 1 = Sigma
-             */
-            auto sigmar = params[1] / r;
-            auto sigmar6 = pow(sigmar, 6.0);
-            auto sigmar12 = sigmar6 * sigmar6;
-            energy = 4.0 * params[0] * (sigmar12 - sigmar6);
-        }
-        break;
-        default:
-            throw(std::runtime_error(fmt::format(
-                "Short-range interaction type {} is not accounted for in PairPotential::analyticShortRangeEnergy().\n",
-                ShortRangeFunctions::forms().keyword(interactionPotential_.form()))));
-    }
+    // Assess stored potential function at specified r
+    auto energy = potentialFunction_.y(r);
 
     // Apply the selected truncation scheme
     if (truncation == PairPotential::ShiftedShortRangeTruncation)
@@ -196,40 +175,8 @@ double PairPotential::analyticShortRangeEnergy(double r, PairPotential::ShortRan
 // Return analytic short range force
 double PairPotential::analyticShortRangeForce(double r, PairPotential::ShortRangeTruncationScheme truncation) const
 {
-    auto &params = interactionPotential_.parameters();
-
-    auto force = 0.0;
-    switch (interactionPotential_.form())
-    {
-        case (ShortRangeFunctions::Form::None):
-            break;
-        case (ShortRangeFunctions::Form::LennardJones):
-        case (ShortRangeFunctions::Form::LennardJonesGeometric):
-        {
-            /*
-             * Standard Lennard-Jones potential
-             * Parameter 0 = Epsilon
-             * Parameter 1 = Sigma
-             */
-
-            /*
-             * dU/dr = -48*epsilon*((sigma**12/x**13)-0.5*(sigma**6/x**7))
-             *
-             *                        signa**6   (   sigma**6        )
-             *       = 48 * epsilon * -------- * ( - --------- + 0.5 ) * r**-1
-             *                          r**6     (     r**6          )
-             */
-
-            auto sigmar = params[1] / r;
-            auto sigmar6 = pow(sigmar, 6.0);
-            force = -48.0 * params[0] * sigmar6 * (-sigmar6 + 0.5) / r;
-        }
-        break;
-        default:
-            throw(std::runtime_error(fmt::format(
-                "Short-range interaction type {} is not accounted for in PairPotential::analyticShortRangeForce().\n",
-                ShortRangeFunctions::forms().keyword(interactionPotential_.form()))));
-    }
+    // Assess stored potential function derivative at specified r and negate to get force
+    auto force = -potentialFunction_.dYdX(r);
 
     // Apply the selected truncation scheme
     if (truncation == PairPotential::ShiftedShortRangeTruncation)
@@ -335,18 +282,14 @@ double PairPotential::delta() const { return delta_; }
 // (Re)generate potential from current parameters
 void PairPotential::calculateUOriginal(bool recalculateUFull)
 {
-    double r;
-
+    // Loop over points
     for (auto n = 1; n < nPoints_; ++n)
     {
-        r = n * delta_;
+        auto r = n * delta_;
         uOriginal_.xAxis(n) = r;
 
-        // Construct potential
-        uOriginal_.value(n) = 0.0;
-
-        // Short-range potential contribution
-        uOriginal_.value(n) += analyticShortRangeEnergy(r);
+        // Set short-range potential contribution
+        uOriginal_.value(n) = analyticShortRangeEnergy(r);
 
         // -- Add Coulomb contribution
         if (includeAtomTypeCharges_)

diff --git a/src/classes/pairPotential.h b/src/classes/pairPotential.h
@@ -5,6 +5,7 @@
 
 #include "data/ff/ff.h"
 #include "math/data1D.h"
+#include "math/function1D.h"
 #include "math/interpolator.h"
 #include <memory>
 
@@ -19,7 +20,7 @@ class PairPotential
     public:
     PairPotential();
     PairPotential(const std::shared_ptr<AtomType> &typeI, const std::shared_ptr<AtomType> &typeJ, bool includeCharges);
-    PairPotential(std::string_view nameI, std::string_view nameJ, const InteractionPotential<ShortRangeFunctions> &potential);
+    PairPotential(std::string_view nameI, std::string_view nameJ, const InteractionPotential<Functions1D> &potential);
     // Coulomb Truncation Scheme enum
     enum CoulombTruncationScheme
     {
@@ -78,8 +79,9 @@ class PairPotential
     private:
     // Names reflecting source parameters
     std::string nameI_, nameJ_;
-    // Interaction potential
-    InteractionPotential<ShortRangeFunctions> interactionPotential_;
+    // Interaction potential and function
+    InteractionPotential<Functions1D> interactionPotential_;
+    Function1DWrapper potentialFunction_;
     // Charge on I (taken from AtomType)
     double chargeI_{0.0};
     // Charge on J (taken from AtomType)
@@ -97,8 +99,8 @@ class PairPotential
     // Return name for second source parameters
     std::string_view nameJ() const;
     // Return interaction potential
-    InteractionPotential<ShortRangeFunctions> &interactionPotential();
-    const InteractionPotential<ShortRangeFunctions> &interactionPotential() const;
+    InteractionPotential<Functions1D> &interactionPotential();
+    const InteractionPotential<Functions1D> &interactionPotential() const;
     // Set charge I
     void setChargeI(double value);
     // Return charge I

diff --git a/src/classes/pairPotentialOverride.cpp b/src/classes/pairPotentialOverride.cpp
@@ -4,11 +4,11 @@
 #include "classes/pairPotentialOverride.h"
 #include "classes/atomType.h"
 
-PairPotentialOverride::PairPotentialOverride() : interactionPotential_(ShortRangeFunctions::Form::None) {}
+PairPotentialOverride::PairPotentialOverride() : interactionPotential_(Functions1D::Form::None) {}
 
 PairPotentialOverride::PairPotentialOverride(std::string_view matchI, std::string_view matchJ,
                                              PairPotentialOverride::PairPotentialOverrideType overrideType,
-                                             const InteractionPotential<ShortRangeFunctions> &potential)
+                                             const InteractionPotential<Functions1D> &potential)
     : matchI_(matchI), matchJ_(matchJ), type_(overrideType), interactionPotential_(potential)
 {
 }
@@ -41,11 +41,8 @@ void PairPotentialOverride::setType(PairPotentialOverrideType overrideType) { ty
 PairPotentialOverride::PairPotentialOverrideType PairPotentialOverride::type() const { return type_; }
 
 // Return interaction potential
-InteractionPotential<ShortRangeFunctions> &PairPotentialOverride::interactionPotential() { return interactionPotential_; }
-const InteractionPotential<ShortRangeFunctions> &PairPotentialOverride::interactionPotential() const
-{
-    return interactionPotential_;
-}
+InteractionPotential<Functions1D> &PairPotentialOverride::interactionPotential() { return interactionPotential_; }
+const InteractionPotential<Functions1D> &PairPotentialOverride::interactionPotential() const { return interactionPotential_; }
 
 /*
  * I/O
@@ -60,12 +57,12 @@ SerialisedValue PairPotentialOverride::serialise() const
     value["matchJ"] = matchJ_;
     value["type"] = pairPotentialOverrideTypes().keyword(type_);
 
-    value["form"] = ShortRangeFunctions::forms().keyword(interactionPotential_.form());
+    value["form"] = Functions1D::forms().keyword(interactionPotential_.form());
     auto &values = interactionPotential().parameters();
     if (!values.empty())
     {
         SerialisedValue interactionParameters;
-        auto &parameters = ShortRangeFunctions::parameters(interactionPotential_.form());
+        auto &parameters = Functions1D::parameters(interactionPotential_.form());
         for (auto &&[parameter, value] : zip(parameters, values))
             interactionParameters[parameter] = value;
         value["parameters"] = interactionParameters;
@@ -84,15 +81,15 @@ void PairPotentialOverride::deserialise(const SerialisedValue &node)
                              [this](const auto node)
                              { type_ = pairPotentialOverrideTypes().enumeration(std::string(node.as_string())); });
 
-    Serialisable::optionalOn(
-        node, "form",
-        [this](const auto node)
-        { interactionPotential_.setForm(ShortRangeFunctions::forms().enumeration(std::string(node.as_string()))); });
+    Serialisable::optionalOn(node, "form",
+                             [this](const auto node) {
+                                 interactionPotential_.setForm(Functions1D::forms().enumeration(std::string(node.as_string())));
+                             });
 
     Serialisable::optionalOn(node, "parameters",
                              [this](const auto node)
                              {
-                                 auto &parameters = ShortRangeFunctions::parameters(interactionPotential_.form());
+                                 auto &parameters = Functions1D::parameters(interactionPotential_.form());
                                  std::vector<double> values;
                                  std::transform(parameters.begin(), parameters.end(), std::back_inserter(values),
                                                 [&node](const auto parameter) { return node.at(parameter).as_floating(); });

diff --git a/src/classes/pairPotentialOverride.h b/src/classes/pairPotentialOverride.h
@@ -4,6 +4,7 @@
 #pragma once
 
 #include "classes/atomType.h"
+#include "math/function1D.h"
 #include <memory>
 #include <string>
 
@@ -23,15 +24,15 @@ class PairPotentialOverride : public Serialisable<>
     PairPotentialOverride();
     PairPotentialOverride(std::string_view matchI, std::string_view matchJ,
                           PairPotentialOverride::PairPotentialOverrideType overrideType,
-                          const InteractionPotential<ShortRangeFunctions> &potential);
+                          const InteractionPotential<Functions1D> &potential);
 
     private:
     // AtomType names to match
     std::string matchI_{"*"}, matchJ_{"*"};
     // Override type
     PairPotentialOverrideType type_{PairPotentialOverrideType::Add};
     // Interaction potential
-    InteractionPotential<ShortRangeFunctions> interactionPotential_;
+    InteractionPotential<Functions1D> interactionPotential_;
 
     public:
     // Set first AtomType name to match
@@ -47,8 +48,8 @@ class PairPotentialOverride : public Serialisable<>
     // Return override type
     PairPotentialOverrideType type() const;
     // Return interaction potential
-    InteractionPotential<ShortRangeFunctions> &interactionPotential();
-    const InteractionPotential<ShortRangeFunctions> &interactionPotential() const;
+    InteractionPotential<Functions1D> &interactionPotential();
+    const InteractionPotential<Functions1D> &interactionPotential() const;
 
     /*
      * I/O

diff --git a/src/classes/shortRangeFunctions.cpp b/src/classes/shortRangeFunctions.cpp
@@ -36,32 +36,32 @@ std::optional<int> ShortRangeFunctions::parameterIndex(Form form, std::string_vi
  */
 
 // Combine parameters for the two atom types using suitable rules
-InteractionPotential<ShortRangeFunctions> ShortRangeFunctions::combine(const InteractionPotential<ShortRangeFunctions> &srI,
-                                                                       const InteractionPotential<ShortRangeFunctions> &srJ)
+InteractionPotential<Functions1D> ShortRangeFunctions::combine(const InteractionPotential<ShortRangeFunctions> &srI,
+                                                               const InteractionPotential<ShortRangeFunctions> &srJ)
 {
     // Equivalent functional forms?
     if (srI.form() == srJ.form())
     {
         switch (srI.form())
         {
             case (Form::None):
-                break;
+                return {Functions1D::Form::None};
             case (Form::LennardJones):
                 /*
                  * Combine parameters (Lorentz-Berthelot):
                  * Parameter 0 = Epsilon
                  * Parameter 1 = Sigma
                  */
-                return {srI.form(), std::vector<double>{sqrt(srI.parameters()[0] * srJ.parameters()[0]),
-                                                        (srI.parameters()[1] + srJ.parameters()[1]) * 0.5}};
+                return {Functions1D::Form::LennardJones126,
+                        {sqrt(srI.parameters()[0] * srJ.parameters()[0]), (srI.parameters()[1] + srJ.parameters()[1]) * 0.5}};
             case (Form::LennardJonesGeometric):
                 /*
                  * Combine parameters (Geometric):
                  * Parameter 0 = Epsilon
                  * Parameter 1 = Sigma
                  */
-                return {srI.form(), std::vector<double>{sqrt(srI.parameters()[0] * srJ.parameters()[0]),
-                                                        sqrt(srI.parameters()[1] * srJ.parameters()[1])}};
+                return {Functions1D::Form::LennardJones126,
+                        {sqrt(srI.parameters()[0] * srJ.parameters()[0]), sqrt(srI.parameters()[1] * srJ.parameters()[1])}};
                 break;
             default:
                 throw(std::runtime_error(
@@ -85,8 +85,8 @@ InteractionPotential<ShortRangeFunctions> ShortRangeFunctions::combine(const Int
              * Parameter 0 = Epsilon
              * Parameter 1 = Sigma
              */
-            return {Form::LennardJones, std::vector<double>{sqrt(srI.parameters()[0] * srJ.parameters()[0]),
-                                                            (srI.parameters()[1] + srJ.parameters()[1]) * 0.5}};
+            return {Functions1D::Form::LennardJones126,
+                    {sqrt(srI.parameters()[0] * srJ.parameters()[0]), (srI.parameters()[1] + srJ.parameters()[1]) * 0.5}};
         }
         else
             Messenger::error("Can't combine parameters between dislike ShortRangeFunctions {} and {}.\n",

diff --git a/src/classes/shortRangeFunctions.h b/src/classes/shortRangeFunctions.h
@@ -5,6 +5,7 @@
 
 #include "base/enumOptions.h"
 #include "classes/interactionPotential.h"
+#include "math/function1D.h"
 #include <optional>
 
 // Short-range functional forms
@@ -31,6 +32,6 @@ class ShortRangeFunctions
      */
     public:
     // Combine parameters for the two atom types using suitable rules
-    static InteractionPotential<ShortRangeFunctions> combine(const InteractionPotential<ShortRangeFunctions> &srI,
-                                                             const InteractionPotential<ShortRangeFunctions> &srJ);
+    static InteractionPotential<Functions1D> combine(const InteractionPotential<ShortRangeFunctions> &srI,
+                                                     const InteractionPotential<ShortRangeFunctions> &srJ);
 };