refactor: Update ForceKernel and add water forces unit tests

benchmark/modules/forces.cpp

@@ -61,7 +61,7 @@ static void BM_CalculateForces_TotalIntraMolecular(benchmark::State &state)
 {
     Problem<speciesType, population> problemDef;
     auto *cfg = problemDef.configuration();
-    std::vector<Vector3> forces(cfg->nAtoms());
+    std::vector<Vector3> forces;
     auto forceKernel = createForceKernel(problemDef);
     for (auto _ : state)
         forceKernel->totalForces(forces, forces,
@@ -84,7 +84,7 @@ static void BM_CalculateForces_TotalInterAtomic(benchmark::State &state)
 {
     Problem<speciesType, population> problemDef;
     auto *cfg = problemDef.configuration();
-    std::vector<Vector3> forces(cfg->nAtoms());
+    std::vector<Vector3> forces;
     auto forceKernel = createForceKernel(problemDef);
     for (auto _ : state)
         forceKernel->totalForces(forces, forces, Kernel::ExcludeGeometric);
@@ -95,10 +95,10 @@ static void BM_CalculateForces_TotalForces(benchmark::State &state)
 {
     Problem<speciesType, population> problemDef;
     auto *cfg = problemDef.configuration();
-    std::vector<Vector3> forces(cfg->nAtoms());
-    const PotentialMap &potentialMap = problemDef.dissolve().potentialMap();
+    std::vector<Vector3> forces;
+    auto forceKernel = KernelProducer::forceKernel(cfg, problemDef.dissolve().potentialMap());
     for (auto _ : state)
-        ForcesModule::totalForces(cfg, potentialMap, ForcesModule::ForceCalculationType::Full, forces, forces);
+        forceKernel->totalForces(forces, forces);
 }
 // small molecule benchmarks
 BENCHMARK_TEMPLATE(BM_CalculateForces_SpeciesBond, SpeciesType::SmallMolecule, SpeciesPopulation::Small);

src/classes/pairPotential.cpp

@@ -10,6 +10,7 @@
 VersionCounter PairPotential::coreDefinitionsVersion_ = ++VersionCounter();
 double PairPotential::range_ = 15.0;
 double PairPotential::delta_ = 0.005;
+int PairPotential::nPoints_ = int(15.0 / 0.005);
 double PairPotential::rDelta_ = 1.0 / 0.005;
 PairPotential::CoulombTruncationScheme PairPotential::coulombTruncationScheme_ = PairPotential::ShiftedCoulombTruncation;
 PairPotential::ShortRangeTruncationScheme PairPotential::shortRangeTruncationScheme_ =
@@ -68,10 +69,13 @@ void PairPotential::setRange(double range, std::optional<double> delta)
 {
     range_ = range;
     if (delta)
-    {
         delta_ = *delta;
-        rDelta_ = 1.0 / delta_;
-    }
+
+    // Determine actual delta to cover requested pair potential range
+    nPoints_ = int(range_ / delta_);
+    delta_ += (range_ - (nPoints_ * delta_)) / nPoints_;
+    rDelta_ = 1.0 / delta_;
+
     ++coreDefinitionsVersion_;
 }
 
@@ -261,19 +265,15 @@ void PairPotential::tabulate()
 {
     // Are we already up do date?
     if (version_ == coreDefinitionsVersion_)
-    {
-        printf("SAME\n");
         return;
-    }
 
     // Precalculate some quantities
     shortRangeEnergyAtCutoff_ = analyticShortRangeEnergy(range_, PairPotential::NoShortRangeTruncation);
     shortRangeForceAtCutoff_ = analyticShortRangeForce(range_, PairPotential::NoShortRangeTruncation);
 
     // Set up containers
-    const auto nPoints = int(range_ / delta_);
-    referenceShortRangePotential_.initialise(nPoints);
-    for (auto n = 0; n < nPoints; ++n)
+    referenceShortRangePotential_.initialise(nPoints_);
+    for (auto n = 0; n < nPoints_; ++n)
         referenceShortRangePotential_.xAxis()[n] = n * delta_;
     coulombPotential_ = referenceShortRangePotential_;
     additionalShortRangePotential_ = referenceShortRangePotential_;

src/classes/pairPotential.h

@@ -54,6 +54,8 @@ class PairPotential : Serialisable<>
     static VersionCounter coreDefinitionsVersion_;
     // Pair potential range
     static double range_;
+    // Number of points in each pair potential
+    static int nPoints_;
     // Pair potential delta
     static double delta_, rDelta_;
     // Truncation scheme to apply to short-range part of potential

src/kernels/force.cpp

@@ -133,14 +133,20 @@ void ForceKernel::extendedForces(const Molecule &mol, std::vector<Vector3> &f) c
  */
 
 // Calculate total forces in the world
-void ForceKernel::totalForces(std::vector<Vector3> &fUnbound, std::vector<Vector3> &fBound,
+void ForceKernel::totalForces(std::vector<Vector3> &ppForceVector, std::vector<Vector3> &geometryForceVector,
                               Flags<Kernel::CalculationFlags> flags) const
 {
+    // Resize and zero force arrays
+    ppForceVector.resize(configuration_->nAtoms());
+    geometryForceVector.resize(configuration_->nAtoms());
+    std::fill(ppForceVector.begin(), ppForceVector.end(), Vector3());
+    std::fill(geometryForceVector.begin(), geometryForceVector.end(), Vector3());
+
     auto &cellArray = configuration_->cells();
     auto &molecules = configuration_->molecules();
 
-    auto combinableUnbound = Kernel::createCombinableVector3(fUnbound);
-    auto combinableBound = Kernel::createCombinableVector3(fBound);
+    auto combinablePP = Kernel::createCombinableVector3(ppForceVector);
+    auto combinableGeometric = Kernel::createCombinableVector3(geometryForceVector);
 
     // Pair potential forces between different molecules
     if (flags.isNotSet(Kernel::ExcludeInterMolecularPairPotential))
@@ -149,7 +155,7 @@ void ForceKernel::totalForces(std::vector<Vector3> &fUnbound, std::vector<Vector
         auto unaryOp = [&](const int id)
         {
             auto *cellI = cellArray.cell(id);
-            auto &fLocal = combinableUnbound.local();
+            auto &fLocal = combinablePP.local();
 
             // Interatomic interactions between atoms in this cell, excluding those within the same molecule
             dissolve::for_each_pair(ParallelPolicies::seq, cellI->atoms(),
@@ -179,8 +185,8 @@ void ForceKernel::totalForces(std::vector<Vector3> &fUnbound, std::vector<Vector
     // Other molecule forces
     auto moleculeForceOperator = [&](const auto &mol)
     {
-        auto &fLocalUnbound = combinableUnbound.local();
-        auto &fLocalBound = combinableBound.local();
+        auto &fLocalUnbound = combinablePP.local();
+        auto &fLocalBound = combinableGeometric.local();
 
         auto offset = mol->globalAtomOffset();
 
@@ -210,6 +216,129 @@ void ForceKernel::totalForces(std::vector<Vector3> &fUnbound, std::vector<Vector
 
     dissolve::for_each(ParallelPolicies::par, molecules.begin(), molecules.end(), moleculeForceOperator);
 
-    combinableUnbound.finalize();
-    combinableBound.finalize();
+    combinablePP.finalize();
+    combinableGeometric.finalize();
+
+    // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units)
+    std::transform(ppForceVector.begin(), ppForceVector.end(), ppForceVector.begin(), [](auto f) { return f * 100.0; });
+    std::transform(geometryForceVector.begin(), geometryForceVector.end(), geometryForceVector.begin(),
+                   [](auto f) { return f * 100.0; });
 }
+
+// Calculate total forces with simple loops for testing
+void ForceKernel::totalForcesSimple(std::vector<Vector3> &ppForceVector, std::vector<Vector3> &geometryForceVector,
+                                    Flags<Kernel::CalculationFlags> flags) const
+{
+    // Resize and zero force arrays
+    ppForceVector.resize(configuration_->nAtoms());
+    geometryForceVector.resize(configuration_->nAtoms());
+    std::fill(ppForceVector.begin(), ppForceVector.end(), Vector3());
+    std::fill(geometryForceVector.begin(), geometryForceVector.end(), Vector3());
+
+    const auto &molecules = configuration_->molecules();
+    std::shared_ptr<Molecule> molN, molM;
+
+    // Calculate interatomic and intramolecular energy in a loop over defined Molecules
+    for (auto n = 0; n < configuration_->nMolecules(); ++n)
+    {
+        molN = molecules[n];
+        auto offsetN = molN->globalAtomOffset();
+
+        // Intramolecular forces (excluding bound terms) in molecule N
+        if (flags.isNotSet(Kernel::CalculationFlags::ExcludeIntraMolecularPairPotential))
+            for (auto ii = 0; ii < molN->nAtoms() - 1; ++ii)
+            {
+                auto *i = molN->atom(ii);
+
+                for (auto jj = ii + 1; jj < molN->nAtoms(); ++jj)
+                {
+                    auto *j = molN->atom(jj);
+
+                    // Get intramolecular scaling of atom pair
+                    auto &&[scalingType, elec14, vdw14] = i->scaling(j);
+
+                    if (scalingType == SpeciesAtom::ScaledInteraction::Excluded)
+                        continue;
+
+                    // Determine final forces
+                    auto vij = box_->minimumVector(i->r(), j->r());
+                    auto magjisq = vij.magnitudeSq();
+                    if (magjisq > cutoffDistanceSquared_)
+                        continue;
+                    auto r = sqrt(magjisq);
+                    vij /= r;
+
+                    if (scalingType == SpeciesAtom::ScaledInteraction::NotScaled)
+                        vij *= potentialMap_.analyticForce(*molN->atom(ii), *molN->atom(jj), r);
+                    else if (scalingType == SpeciesAtom::ScaledInteraction::Scaled)
+                        vij *= potentialMap_.analyticForce(*molN->atom(ii), *molN->atom(jj), r, elec14, vdw14);
+
+                    ppForceVector[offsetN + ii] -= vij;
+                    ppForceVector[offsetN + jj] += vij;
+                }
+            }
+
+        // Forces between molecule N and molecule M
+        if (flags.isNotSet(Kernel::CalculationFlags::ExcludeInterMolecularPairPotential))
+            for (auto m = n + 1; m < configuration_->nMolecules(); ++m)
+            {
+                molM = molecules[m];
+                auto offsetM = molM->globalAtomOffset();
+
+                // Double loop over atoms
+                for (auto ii = 0; ii < molN->nAtoms(); ++ii)
+                {
+                    auto *i = molN->atom(ii);
+
+                    for (auto jj = 0; jj < molM->nAtoms(); ++jj)
+                    {
+                        auto *j = molM->atom(jj);
+
+                        // Determine final forces
+                        auto vij = box_->minimumVector(i->r(), j->r());
+                        auto magjisq = vij.magnitudeSq();
+                        if (magjisq > cutoffDistanceSquared_)
+                            continue;
+                        auto r = sqrt(magjisq);
+                        vij /= r;
+
+                        vij *= potentialMap_.analyticForce(*i, *j, r);
+
+                        ppForceVector[offsetN + ii] -= vij;
+                        ppForceVector[offsetM + jj] += vij;
+                    }
+                }
+            }
+
+        if (flags.isNotSet(Kernel::CalculationFlags::ExcludeGeometric))
+        {
+            // Bond forces
+            for (const auto &bond : molN->species()->bonds())
+                bondForces(bond, *molN->atom(bond.indexI()), offsetN + bond.indexI(), *molN->atom(bond.indexJ()),
+                           offsetN + bond.indexJ(), geometryForceVector);
+
+            // Angle forces
+            for (const auto &angle : molN->species()->angles())
+                angleForces(angle, *molN->atom(angle.indexI()), offsetN + angle.indexI(), *molN->atom(angle.indexJ()),
+                            offsetN + angle.indexJ(), *molN->atom(angle.indexK()), offsetN + angle.indexK(),
+                            geometryForceVector);
+
+            // Torsion forces
+            for (const auto &torsion : molN->species()->torsions())
+                torsionForces(torsion, *molN->atom(torsion.indexI()), offsetN + torsion.indexI(), *molN->atom(torsion.indexJ()),
+                              offsetN + torsion.indexJ(), *molN->atom(torsion.indexK()), offsetN + torsion.indexK(),
+                              *molN->atom(torsion.indexL()), offsetN + torsion.indexL(), geometryForceVector);
+
+            // Improper forces
+            for (const auto &imp : molN->species()->impropers())
+                improperForces(imp, *molN->atom(imp.indexI()), offsetN + imp.indexI(), *molN->atom(imp.indexJ()),
+                               offsetN + imp.indexJ(), *molN->atom(imp.indexK()), offsetN + imp.indexK(),
+                               *molN->atom(imp.indexL()), offsetN + imp.indexL(), geometryForceVector);
+        }
+    }
+
+    // Convert forces to 10J/mol
+    std::transform(ppForceVector.begin(), ppForceVector.end(), ppForceVector.begin(), [](auto &f) { return f * 100.0; });
+    std::transform(geometryForceVector.begin(), geometryForceVector.end(), geometryForceVector.begin(),
+                   [](auto &f) { return f * 100.0; });
+}

src/kernels/force.h

@@ -59,7 +59,10 @@ class ForceKernel : public GeometryKernel
      * Totals
      */
     public:
-    // Calculate total forces in the world
-    void totalForces(std::vector<Vector3> &fUnbound, std::vector<Vector3> &fBound,
+    // Calculate total forces
+    void totalForces(std::vector<Vector3> &ppForceVector, std::vector<Vector3> &geometryForceVector,
                      Flags<Kernel::CalculationFlags> flags = {}) const;
+    // Calculate total forces with simple loops for testing
+    void totalForcesSimple(std::vector<Vector3> &ppForceVector, std::vector<Vector3> &geometryForceVector,
+                           Flags<Kernel::CalculationFlags> flags = {}) const;
 };

src/math/interpolator.cpp

@@ -284,7 +284,11 @@ void Interpolator::interpolate()
         interpolateLinear();
     else if (scheme_ == Interpolator::ThreePointInterpolation)
     {
-        // No setup required
+        // Store delta y for each interval in a_
+        a_.clear();
+        a_.resize(x_.size() - 1);
+        for (auto i = 0; i < x_.size() - 1; ++i)
+            a_[i] = y_[i + 1] - y_[i];
     }
 }
 
@@ -322,36 +326,36 @@ double Interpolator::y(double x, int interval) const
         if (x >= x_.back())
             return y_.back();
 
-        double h = x - x_[interval];
-        double hh = h * h;
+        auto h = x - x_[interval];
+        auto hh = h * h;
         return a_[interval] + b_[interval] * h + c_[interval] * hh + d_[interval] * hh * h;
     }
     //	else if (scheme_ == Interpolator::ConstrainedSplineInterpolation)
     //	{
-    //		double h = x;
-    //		double hh = h*h;
+    //		auto h = x;
+    //		auto hh = h*h;
     //		return a_[interval] + b_[interval]*h + c_[interval]*hh + d_[interval]*hh*h;
     //	}
     else if (scheme_ == Interpolator::LinearInterpolation)
     {
         if (interval >= (x_.size() - 1))
             return y_.back();
 
-        double delta = (x - x_[interval]) / h_[interval];
+        auto delta = (x - x_[interval]) / h_[interval];
         return y_[interval] + delta * a_[interval];
     }
     else if (scheme_ == Interpolator::ThreePointInterpolation)
     {
         if (interval >= (x_.size() - 3))
-            return y_.back();
+            return y_.back() + (x - x_.back()) * a_.back() / h_.back();
 
-        double ppp = (x - x_[interval]) / h_[interval];
+        auto ppp = (x - x_[interval]) / h_[interval];
 
-        double vk0 = y_[interval];
-        double vk1 = y_[interval + 1];
-        double vk2 = y_[interval + 2];
-        double t1 = vk0 + (vk1 - vk0) * ppp;
-        double t2 = vk1 + (vk2 - vk1) * (ppp - 1.0);
+        auto vk0 = y_[interval];
+        auto vk1 = y_[interval + 1];
+        auto vk2 = y_[interval + 2];
+        auto t1 = vk0 + (vk1 - vk0) * ppp;
+        auto t2 = vk1 + (vk2 - vk1) * (ppp - 1.0);
         return t1 + (t2 - t1) * ppp * 0.5;
     }
 
@@ -408,14 +412,14 @@ double Interpolator::approximate(const Data1D &data, double x)
             left = i;
     }
 
-    double ppp = (x - xData[left]) / (xData[right] - xData[left]);
+    auto ppp = (x - xData[left]) / (xData[right] - xData[left]);
 
-    double vk0 = yData[left];
-    double vk1 = yData[left + 1];
-    double vk2 = yData[left + 2];
+    auto vk0 = yData[left];
+    auto vk1 = yData[left + 1];
+    auto vk2 = yData[left + 2];
 
-    double t1 = vk0 + (vk1 - vk0) * ppp;
-    double t2 = vk1 + (vk2 - vk1) * (ppp - 1.0);
+    auto t1 = vk0 + (vk1 - vk0) * ppp;
+    auto t2 = vk1 + (vk2 - vk1) * (ppp - 1.0);
 
     return t1 + (t2 - t1) * ppp * 0.5;
 }

src/modules/forces/forces.h

@@ -5,8 +5,8 @@
 
 #include "io/export/forces.h"
 #include "io/import/forces.h"
+#include "kernels/common.h"
 #include "module/module.h"
-#include <memory>
 
 // Forward Declarations
 class Molecule;
@@ -45,24 +45,4 @@ class ForcesModule : public Module
     public:
     // Run set-up stage
     bool setUp(Dissolve &dissolve, Flags<KeywordBase::KeywordSignal> actionSignals) override;
-
-    /*
-     * Functions
-     */
-    public:
-    // Force calculation type
-    enum class ForceCalculationType
-    {
-        Full,
-        PairPotentialOnly,
-        IntraMolecularFull,
-        IntraMolecularGeometry
-    };
-    // Calculate total forces within the specified Configuration
-    static void totalForces(Configuration *cfg, const PotentialMap &potentialMap, ForceCalculationType calculationType,
-                            std::vector<Vector3> &fUnbound, std::vector<Vector3> &fBound);
-    // Calculate forces acting on specific Molecules within the specified Configuration (arising from all atoms)
-    static void totalForces(Configuration *cfg, const std::vector<const Molecule *> &targetMolecules,
-                            const PotentialMap &potentialMap, ForceCalculationType calculationType,
-                            std::vector<Vector3> &fUnbound, std::vector<Vector3> &fBound);
 };