Region Potentials (#1658)

src/kernels/potentials/CMakeLists.txt

@@ -3,12 +3,14 @@ add_library(
   base.cpp
   directional.cpp
   producer.cpp
+  regional.cpp
   restraint.cpp
   simple.cpp
   types.cpp
   base.h
   directional.h
   producer.h
+  regional.h
   restraint.h
   simple.h
   types.h

src/kernels/potentials/producer.cpp

@@ -3,6 +3,7 @@
 
 #include "kernels/potentials/producer.h"
 #include "kernels/potentials/directional.h"
+#include "kernels/potentials/regional.h"
 #include "kernels/potentials/simple.h"
 
 // External Potential Producer
@@ -17,6 +18,8 @@ std::unique_ptr<ExternalPotential> create(ExternalPotentialTypes::ExternalPotent
             return std::make_unique<SimplePotential>();
         case (ExternalPotentialTypes::ExternalPotentialType::Directional):
             return std::make_unique<DirectionalPotential>();
+        case (ExternalPotentialTypes::ExternalPotentialType::Regional):
+            return std::make_unique<RegionalPotential>();
         default:
             throw(std::runtime_error(fmt::format("Creation of external potential type '{}' not implemented.\n",
                                                  ExternalPotentialTypes::keyword(type))));

src/kernels/potentials/regional.cpp

@@ -0,0 +1,144 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2023 Team Dissolve and contributors
+
+#include "kernels/potentials/regional.h"
+#include "classes/atom.h"
+#include "classes/box.h"
+#include "expression/variable.h"
+#include "kernels/potentials/types.h"
+#include "keywords/double.h"
+#include "keywords/nodeValue.h"
+#include "templates/algorithms.h"
+
+/*
+ * Regional Potential Voxel Kernel
+ */
+
+RegionalPotentialVoxelKernel::RegionalPotentialVoxelKernel(std::string_view expressionString,
+                                                           std::vector<std::shared_ptr<ExpressionVariable>> parameters,
+                                                           double minimumValue, double maximumValue, double valueOffset,
+                                                           double penaltyPower)
+{
+    // Create our local variables
+    x_ = expression_.addLocalVariable("x");
+    y_ = expression_.addLocalVariable("y");
+    z_ = expression_.addLocalVariable("z");
+    xFrac_ = expression_.addLocalVariable("xFrac");
+    yFrac_ = expression_.addLocalVariable("yFrac");
+    zFrac_ = expression_.addLocalVariable("zFrac");
+
+    // Set the expression
+    expression_.set(expressionString, std::move(parameters));
+
+    // Set limits
+    minimumValue_ = minimumValue;
+    maximumValue_ = maximumValue;
+    valueOffset_ = valueOffset;
+    penaltyPower_ = penaltyPower;
+}
+
+// Set voxel position variables
+void RegionalPotentialVoxelKernel::setVoxelPosition(const Box *box, Vec3<double> r) const
+{
+    x_->setValue(r.x);
+    y_->setValue(r.y);
+    z_->setValue(r.z);
+    box->toFractional(r);
+    xFrac_->setValue(r.x);
+    yFrac_->setValue(r.y);
+    zFrac_->setValue(r.z);
+}
+
+// Return current value of function, applying any threshold penalty
+double RegionalPotentialVoxelKernel::functionValue() const
+{
+    auto x = expression_.asDouble() + valueOffset_;
+    if (x < minimumValue_ || x > maximumValue_)
+        x = pow(x, penaltyPower_);
+    return x - valueOffset_;
+}
+
+// Calculate and store energy and force for the specified voxel centre
+void RegionalPotentialVoxelKernel::energyAndForce(const Box *box, const Vec3<double> &r, double &energy,
+                                                  Vec3<double> &force) const
+{
+    // Energy at the centre of the voxel
+    setVoxelPosition(box, r);
+    energy = functionValue();
+
+    // Force at the centre of the voxel
+    const auto rDelta = 0.01;
+    for (auto a = 0; a < 3; ++a)
+    {
+        // Evaluate r(i) + delta
+        setVoxelPosition(box, r.adjusted(a, rDelta));
+        force[a] = functionValue();
+
+        // Evaluate r(i) - delta
+        setVoxelPosition(box, r.adjusted(a, -rDelta));
+        force[a] -= functionValue();
+
+        // Finalise force
+        force[a] /= -rDelta * 2;
+    }
+}
+
+/*
+ * Regional Potential
+ */
+
+RegionalPotential::RegionalPotential() : ExternalPotential(ExternalPotentialTypes::ExternalPotentialType::Regional) {}
+
+/*
+ * Definition
+ */
+
+// Return functional form of the potential, as a string
+const std::string RegionalPotential::formString() const { return "Expression"; }
+
+// Return parameters of the potential, as a string
+const std::string RegionalPotential::formParametersString() const { return "N/A"; }
+
+// Set up potential for supplied box
+bool RegionalPotential::setUp(const Box *box, double voxelSize,
+                              const std::function<std::shared_ptr<RegionalPotentialVoxelKernel>(void)> &kernelGenerator)
+{
+    // Set fractional voxel sizes
+    Vec3<int> nVoxels;
+    for (auto n = 0; n < 3; ++n)
+        nVoxels.set(n, std::max(int(box->axisLength(n) / voxelSize), 1));
+    voxelSizeFrac_.set(1.0 / nVoxels.x, 1.0 / nVoxels.y, 1.0 / nVoxels.z);
+
+    // Initialise 3D maps
+    energyVoxels_.initialise(nVoxels.x, nVoxels.y, nVoxels.z);
+    forceVoxels_.initialise(nVoxels.x, nVoxels.y, nVoxels.z);
+
+    // Create a voxel combinable and check function
+    auto voxelCombinable = createCombinableVoxelKernel(kernelGenerator);
+    auto voxelFunction = [&](auto triplet, auto x, auto y, auto z)
+    {
+        auto r = box->getReal({(x + 0.5) * voxelSizeFrac_.x, (y + 0.5) * voxelSizeFrac_.y, (z + 0.5) * voxelSizeFrac_.z});
+        voxelCombinable.local()->energyAndForce(box, r, energyVoxels_[triplet], forceVoxels_[triplet]);
+    };
+
+    // Iterate over voxel indices
+    dissolve::for_each_triplet(ParallelPolicies::par, energyVoxels_.beginIndices(), energyVoxels_.endIndices(), voxelFunction);
+
+    return true;
+}
+
+/*
+ * Potential Calculation
+ */
+
+std::tuple<int, int, int> RegionalPotential::voxelIndices(const Atom &i, const Box *box) const
+{
+    auto r = box->getFractional(i.r());
+    return {r.x / voxelSizeFrac_.x, r.y / voxelSizeFrac_.y, r.z / voxelSizeFrac_.z};
+}
+
+// Calculate energy on specified atom
+double RegionalPotential::energy(const Atom &i, const Box *box) const { return energyVoxels_[voxelIndices(i, box)]; }
+
+// Calculate force on specified atom, summing in to supplied vector
+void RegionalPotential::force(const Atom &i, const Box *box, Vec3<double> &f) const { f = forceVoxels_[voxelIndices(i, box)]; }

src/kernels/potentials/regional.h

@@ -0,0 +1,87 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2023 Team Dissolve and contributors
+
+#pragma once
+
+#include "classes/interactionPotential.h"
+#include "classes/region.h"
+#include "kernels/potentials/base.h"
+
+// Regional Potential Voxel Kernel
+class RegionalPotentialVoxelKernel
+{
+    public:
+    explicit RegionalPotentialVoxelKernel(std::string_view expressionString = "",
+                                          std::vector<std::shared_ptr<ExpressionVariable>> = {}, double minimumValue = 0.0,
+                                          double maximumValue = 1.0, double valueOffset = 0.0, double penaltyPower = 1.0);
+
+    protected:
+    // Local variables, set when checking voxels
+    std::shared_ptr<ExpressionVariable> x_, y_, z_, xFrac_, yFrac_, zFrac_;
+    // Expression describing region
+    NodeValue expression_;
+    // Minimum threshold value for function
+    double minimumValue_{0.0};
+    // Maximum threshold value for function
+    double maximumValue_{1.0};
+    // Value offset to use when assessing threshold
+    double valueOffset_{0.0};
+    // Power law to apply when function value is outside of threshold limits
+    double penaltyPower_{1.0};
+
+    private:
+    // Set voxel position variables
+    void setVoxelPosition(const Box *box, Vec3<double> r) const;
+    // Return current value of function, applying any threshold penalty
+    double functionValue() const;
+
+    public:
+    // Calculate and store energy and force for the specified voxel centre
+    void energyAndForce(const Box *box, const Vec3<double> &r, double &energy, Vec3<double> &force) const;
+};
+
+// Regional Potential
+class RegionalPotential : public ExternalPotential
+{
+    public:
+    RegionalPotential();
+    ~RegionalPotential() = default;
+
+    /*
+     * Definition
+     */
+    private:
+    // Fractional voxel size
+    Vec3<double> voxelSizeFrac_;
+    // Vector fields for energy and derived force
+    Array3D<double> energyVoxels_;
+    Array3D<Vec3<double>> forceVoxels_;
+
+    private:
+    // Generate voxel combinable
+    static dissolve::CombinableFunctor<std::shared_ptr<RegionalPotentialVoxelKernel>>
+    createCombinableVoxelKernel(std::function<std::shared_ptr<RegionalPotentialVoxelKernel>(void)> kernelGenerator)
+    {
+        return kernelGenerator;
+    }
+    // Return voxel coordinates of supplied atom
+    std::tuple<int, int, int> voxelIndices(const Atom &i, const Box *box) const;
+
+    public:
+    // Set up potential for supplied box
+    bool setUp(const Box *box, double voxelSize,
+               const std::function<std::shared_ptr<RegionalPotentialVoxelKernel>(void)> &kernelGenerator);
+    // Return functional form of the potential, as a string
+    const std::string formString() const override;
+    // Return parameters of the potential, as a string
+    const std::string formParametersString() const override;
+
+    /*
+     * Potential Calculation
+     */
+    public:
+    // Calculate energy on specified atom
+    double energy(const Atom &i, const Box *box) const override;
+    // Calculate force on specified atom, summing in to supplied vector
+    void force(const Atom &i, const Box *box, Vec3<double> &f) const override;
+};

src/kernels/potentials/types.cpp

@@ -8,7 +8,8 @@ namespace ExternalPotentialTypes
 {
 // External Potential Types
 EnumOptions<ExternalPotentialType> types_("ExternalPotential", {{ExternalPotentialType::Simple, "Simple"},
-                                                                {ExternalPotentialType::Directional, "Directional"}});
+                                                                {ExternalPotentialType::Directional, "Directional"},
+                                                                {ExternalPotentialType::Regional, "Regional"}});
 
 // Return whether the supplied external potential type is valid
 std::optional<ExternalPotentialType> isType(std::string_view name)

src/kernels/potentials/types.h

@@ -13,7 +13,8 @@ namespace ExternalPotentialTypes
 enum class ExternalPotentialType
 {
     Simple,
-    Directional
+    Directional,
+    Regional
 };
 // Return whether the supplied external potential type is valid
 std::optional<ExternalPotentialType> isType(std::string_view name);

src/procedure/nodes/CMakeLists.txt

@@ -45,6 +45,7 @@ add_library(
   process2D.cpp
   process3D.cpp
   regionBase.cpp
+  regionalGlobalPotential.cpp
   registry.cpp
   remove.cpp
   restraintPotential.cpp
@@ -102,6 +103,7 @@ add_library(
   process2D.h
   process3D.h
   regionBase.h
+  regionalGlobalPotential.h
   registry.h
   remove.h
   restraintPotential.h

src/procedure/nodes/node.cpp

@@ -65,6 +65,7 @@ EnumOptions<ProcedureNode::NodeType> ProcedureNode::nodeTypes()
                      {ProcedureNode::NodeType::Process1D, "Process1D"},
                      {ProcedureNode::NodeType::Process2D, "Process2D"},
                      {ProcedureNode::NodeType::Process3D, "Process3D"},
+                     {ProcedureNode::NodeType::RegionalGlobalPotential, "RegionalGlobalPotential"},
                      {ProcedureNode::NodeType::Remove, "Remove"},
                      {ProcedureNode::NodeType::RestraintPotential, "RestraintPotential"},
                      {ProcedureNode::NodeType::RotateFragment, "RotateFragment"},

src/procedure/nodes/node.h

@@ -77,6 +77,7 @@ class ProcedureNode : public std::enable_shared_from_this<ProcedureNode>, public
         Process1D,
         Process2D,
         Process3D,
+        RegionalGlobalPotential,
         Remove,
         RestraintPotential,
         RotateFragment,

src/procedure/nodes/regionalGlobalPotential.cpp

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2023 Team Dissolve and contributors
+
+#include "procedure/nodes/regionalGlobalPotential.h"
+#include "classes/configuration.h"
+#include "kernels/potentials/regional.h"
+#include "keywords/double.h"
+#include "keywords/nodeValue.h"
+
+RegionalGlobalPotentialProcedureNode::RegionalGlobalPotentialProcedureNode()
+    : ProcedureNode(ProcedureNode::NodeType::RegionalGlobalPotential, {ProcedureNode::GenerationContext})
+{
+    keywords_.setOrganisation("Options", "Definition");
+    keywords_.add<NodeValueKeyword>("Expression", "Expression describing region", expression_, this);
+    keywords_.add<DoubleKeyword>("Minimum", "Minimum value for descriptive function defining region", minimumValue_);
+    keywords_.add<DoubleKeyword>("Maximum", "Maximum value for descriptive function defining region", maximumValue_);
+    keywords_.add<DoubleKeyword>("Offset", "Offset to apply to calculated value when assessing threshold and penalty",
+                                 valueOffset_);
+    keywords_.add<DoubleKeyword>("PenaltyPower", "Power to apply to the (offset) value if out of threshold range",
+                                 penaltyPower_);
+
+    keywords_.setOrganisation("Options", "Grid");
+    keywords_.add<DoubleKeyword>("VoxelSize", "Voxel size (length) guiding the coarseness / detail of the region", voxelSize_,
+                                 0.1);
+}
+
+/*
+ * Execute
+ */
+
+// Execute node
+bool RegionalGlobalPotentialProcedureNode::execute(const ProcedureContext &procedureContext)
+{
+    auto *cfg = procedureContext.configuration();
+
+    auto pot = std::make_unique<RegionalPotential>();
+    pot->setUp(cfg->box(), voxelSize_,
+               [&]()
+               {
+                   return std::make_shared<RegionalPotentialVoxelKernel>(expression_.asString(), getParameters(), minimumValue_,
+                                                                         maximumValue_, valueOffset_, penaltyPower_);
+               });
+
+    cfg->addGlobalPotential(std::unique_ptr<ExternalPotential>(std::move(pot)));
+
+    return true;
+}

src/procedure/nodes/regionalGlobalPotential.h

@@ -0,0 +1,26 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2023 Team Dissolve and contributors
+
+#pragma once
+
+#include "kernels/potentials/regional.h"
+#include "procedure/nodes/node.h"
+
+// Regional Global Potential Procedure Node
+class RegionalGlobalPotentialProcedureNode : public RegionalPotentialVoxelKernel, public ProcedureNode
+{
+    public:
+    RegionalGlobalPotentialProcedureNode();
+    ~RegionalGlobalPotentialProcedureNode() override = default;
+
+    private:
+    // Guide voxel size (Angstroms)
+    double voxelSize_{0.5};
+
+    /*
+     * Execute
+     */
+    public:
+    // Execute node
+    bool execute(const ProcedureContext &procedureContext) override;
+};