Extend nblib listed pbc holder and template kernels

diff --git a/api/nblib/listed_forces/calculator.cpp b/api/nblib/listed_forces/calculator.cpp
index 7122001a6fb7efb301a052bad47af96e1a8f5fd5..5ded958e743d7edd5bef5af185ee3c53fca0d0dc 100644 (file)
--- a/api/nblib/listed_forces/calculator.cpp
+++ b/api/nblib/listed_forces/calculator.cpp
@@ -60,7 +60,7 @@ ListedForceCalculator::ListedForceCalculator(const ListedInteractionData& intera
                                              const Box&                   box) :
     numThreads(nthr),
     masterForceBuffer_(bufferSize, Vec3{ 0, 0, 0 }),
-    pbcHolder_(std::make_unique<PbcHolder>(box))
+    pbcHolder_(std::make_unique<PbcHolder>(PbcType::Xyz, box))
 {
     // split up the work
     threadedInteractions_ = splitListedWork(interactions, bufferSize, numThreads);

diff --git a/api/nblib/listed_forces/dataflow.hpp b/api/nblib/listed_forces/dataflow.hpp
index d078e176afea2fa665bc8ed58162efb92915bfd2..1f39f84863e019a2de24c896b4682386d8c68343 100644 (file)
--- a/api/nblib/listed_forces/dataflow.hpp
+++ b/api/nblib/listed_forces/dataflow.hpp
@@ -57,7 +57,6 @@
 #include "nblib/util/util.hpp"
 #include "nblib/pbc.hpp"
 #include "nblib/vector.h"
-#include "gromacs/math/vec.h"
 #include "gromacs/utility/arrayref.h"
 
 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
@@ -101,8 +100,8 @@ auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, Ba
 template <class Buffer, class TwoCenterType, class BasicVector, class Pbc,
           std::enable_if_t<Contains<TwoCenterType, SupportedTwoCenterTypes>{}>* = nullptr>
 inline NBLIB_ALWAYS_INLINE
-auto dispatchInteraction(const InteractionIndex<TwoCenterType>& index,
-                         const std::vector<TwoCenterType>& bondInstances,
+auto dispatchInteraction(InteractionIndex<TwoCenterType> index,
+                         gmx::ArrayRef<const TwoCenterType> bondInstances,
                          gmx::ArrayRef<const BasicVector> x,
                          Buffer* forces,
                          const Pbc& pbc)
@@ -111,11 +110,11 @@ auto dispatchInteraction(const InteractionIndex<TwoCenterType>& index,
 
     int i = std::get<0>(index);
     int j = std::get<1>(index);
-    const gmx::RVec& x1 = x[i];
-    const gmx::RVec& x2 = x[j];
-    const TwoCenterType& bond = bondInstances[std::get<2>(index)];
+    BasicVector x1 = x[i];
+    BasicVector x2 = x[j];
+    TwoCenterType bond = bondInstances[std::get<2>(index)];
 
-    gmx::RVec dx;
+    BasicVector dx;
     // calculate x1 - x2 modulo pbc
     pbc.dxAiuc(x1, x2, dx);
 
@@ -161,13 +160,13 @@ addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters,
 template<class ThreeCenterType, class BasicVector>
 inline NBLIB_ALWAYS_INLINE
 auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
-                        BasicVector* fi, BasicVector* fj, BasicVector* fk)
+                        [[maybe_unused]] const BasicVector& rik, BasicVector* fi, BasicVector* fj, BasicVector* fk)
 {
     using ValueType = BasicVectorValueType_t<BasicVector>;
     // calculate 3-center common quantities: angle between x1-x2 and x2-x3
     // Todo: after sufficient evaluation, switch over to atan2 based algorithm
-    ValueType costh = cos_angle(rij, rkj); /* 25             */
-    ValueType theta = std::acos(costh);    /* 10             */
+    ValueType costh = basicVectorCosAngle(rij, rkj); /* 25 */
+    ValueType theta = std::acos(costh);    /* 10 */
 
     // call type-specific angle kernel, e.g. harmonic, linear, quartic,  etc.
     auto [force, energy] = threeCenterKernel(theta, parameters);
@@ -193,8 +192,8 @@ auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& ri
 template <class Buffer, class ThreeCenterType, class BasicVector, class Pbc,
           std::enable_if_t<Contains<ThreeCenterType, SupportedThreeCenterTypes>{}>* = nullptr>
 inline NBLIB_ALWAYS_INLINE
-auto dispatchInteraction(const InteractionIndex<ThreeCenterType>& index,
-                         const std::vector<ThreeCenterType>& parameters,
+auto dispatchInteraction(InteractionIndex<ThreeCenterType> index,
+                         gmx::ArrayRef<const ThreeCenterType> parameters,
                          gmx::ArrayRef<const BasicVector> x,
                          Buffer* forces,
                          const Pbc& pbc)
@@ -205,18 +204,19 @@ auto dispatchInteraction(const InteractionIndex<ThreeCenterType>& index,
     int i = std::get<0>(index);
     int j = std::get<1>(index);
     int k = std::get<2>(index);
-    const gmx::RVec& xi = x[i];
-    const gmx::RVec& xj = x[j];
-    const gmx::RVec& xk = x[k];
-    const ThreeCenterType& threeCenterParameters = parameters[std::get<3>(index)];
+    BasicVector xi = x[i];
+    BasicVector xj = x[j];
+    BasicVector xk = x[k];
+    ThreeCenterType threeCenterParameters = parameters[std::get<3>(index)];
 
-    gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
+    BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
 
-    gmx::RVec rij, rkj;
-    pbc.dxAiuc(xi, xj, rij); /*  3              */
-    pbc.dxAiuc(xk, xj, rkj); /*  3              */
+    BasicVector rij, rkj, rik;
+    pbc.dxAiuc(xi, xj, rij); /* 3 */
+    pbc.dxAiuc(xk, xj, rkj); /* 3 */
+    pbc.dxAiuc(xi, xk, rik); /* 3 */
 
-    energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
+    energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, rik, &fi, &fj, &fk);
     energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
 
     (*forces)[i] += fi;
@@ -279,35 +279,36 @@ addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
 template <class Buffer, class FourCenterType, class BasicVector, class Pbc,
           std::enable_if_t<Contains<FourCenterType, SupportedFourCenterTypes>{}>* = nullptr>
 inline NBLIB_ALWAYS_INLINE
-auto dispatchInteraction(const InteractionIndex<FourCenterType>& index,
-                         const std::vector<FourCenterType>& parameters,
+auto dispatchInteraction(InteractionIndex<FourCenterType> index,
+                         gmx::ArrayRef<const FourCenterType> parameters,
                          gmx::ArrayRef<const BasicVector> x,
                          Buffer* forces,
                          const Pbc& pbc)
 {
-    KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
+    using RealScalar = BasicVectorValueType_t<BasicVector>;
+    KernelEnergy<RealScalar> energy;
 
     int i = std::get<0>(index);
     int j = std::get<1>(index);
     int k = std::get<2>(index);
     int l = std::get<3>(index);
 
-    const gmx::RVec& xi = x[i];
-    const gmx::RVec& xj = x[j];
-    const gmx::RVec& xk = x[k];
-    const gmx::RVec& xl = x[l];
+    BasicVector xi = x[i];
+    BasicVector xj = x[j];
+    BasicVector xk = x[k];
+    BasicVector xl = x[l];
 
-    gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
+    BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
 
-    gmx::RVec dxIJ, dxKJ, dxKL;
+    BasicVector dxIJ, dxKJ, dxKL;
     pbc.dxAiuc(xi, xj, dxIJ);
     pbc.dxAiuc(xk, xj, dxKJ);
     pbc.dxAiuc(xk, xl, dxKL);
 
-    const FourCenterType& fourCenterTypeParams = parameters[std::get<4>(index)];
+    FourCenterType fourCenterTypeParams = parameters[std::get<4>(index)];
 
-    rvec m, n;
-    real phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
+    BasicVector m, n;
+    RealScalar phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
 
     auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
 
@@ -340,8 +341,8 @@ auto dispatchInteraction(const InteractionIndex<FourCenterType>& index,
 template <class Buffer, class FiveCenterType, class BasicVector, class Pbc,
           std::enable_if_t<Contains<FiveCenterType, SupportedFiveCenterTypes>{}>* = nullptr>
 inline NBLIB_ALWAYS_INLINE
-auto dispatchInteraction(const InteractionIndex<FiveCenterType>& index,
-                         const std::vector<FiveCenterType>& parameters,
+auto dispatchInteraction(InteractionIndex<FiveCenterType> index,
+                         gmx::ArrayRef<const FiveCenterType> parameters,
                          gmx::ArrayRef<const BasicVector> x,
                          Buffer* forces,
                          const Pbc& pbc)
@@ -354,19 +355,19 @@ auto dispatchInteraction(const InteractionIndex<FiveCenterType>& index,
     int l = std::get<3>(index);
     int m = std::get<4>(index);
 
-    const gmx::RVec& xi = x[i];
-    const gmx::RVec& xj = x[j];
-    const gmx::RVec& xk = x[k];
-    const gmx::RVec& xl = x[l];
-    const gmx::RVec& xm = x[m];
+    BasicVector xi = x[i];
+    BasicVector xj = x[j];
+    BasicVector xk = x[k];
+    BasicVector xl = x[l];
+    BasicVector xm = x[m];
 
-    gmx::RVec dxIJ, dxJK, dxKL, dxLM;
+    BasicVector dxIJ, dxJK, dxKL, dxLM;
     pbc.dxAiuc(xi, xj, dxIJ);
     pbc.dxAiuc(xj, xk, dxJK);
     pbc.dxAiuc(xk, xl, dxKL);
     pbc.dxAiuc(xl, xm, dxLM);
 
-    const FiveCenterType& fiveCenterTypeParams = parameters[std::get<5>(index)];
+    FiveCenterType fiveCenterTypeParams = parameters[std::get<5>(index)];
 
     // this dispatch function is not in use yet, because CMap is not yet implemented
     // we don't want to add [[maybe_unused]] in the signature
@@ -389,8 +390,8 @@ auto dispatchInteraction(const InteractionIndex<FiveCenterType>& index,
  * \return Computed kernel energies
  */
 template <class Index, class InteractionType, class Buffer, class Pbc>
-auto computeForces(const std::vector<Index>& indices,
-                   const std::vector<InteractionType>& interactionParameters,
+auto computeForces(gmx::ArrayRef<const Index> indices,
+                   gmx::ArrayRef<const InteractionType> parameters,
                    gmx::ArrayRef<const Vec3> x,
                    Buffer* forces,
                    const Pbc& pbc)
@@ -399,7 +400,7 @@ auto computeForces(const std::vector<Index>& indices,
 
     for (const auto& index : indices)
     {
-        energy += dispatchInteraction(index, interactionParameters, x, forces, pbc);
+        energy += dispatchInteraction(index, parameters, x, forces, pbc);
     }
 
     return energy;
@@ -427,7 +428,9 @@ auto reduceListedForces(const ListedInteractionData& interactions,
     auto computeForceType = [forces, &x, &energies, &pbc](const auto& interactionElement) {
         using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
 
-        KernelEnergy<ValueType> energy = computeForces(interactionElement.indices, interactionElement.parameters, x, forces, pbc);
+        gmx::ArrayRef<const InteractionIndex<InteractionType>> indices(interactionElement.indices);
+        gmx::ArrayRef<const InteractionType>                   parameters(interactionElement.parameters);
+        KernelEnergy<ValueType> energy = computeForces(indices, parameters, x, forces, pbc);
 
         energies[CarrierIndex<InteractionType>{}] += energy.carrier();
         energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();

diff --git a/api/nblib/listed_forces/kernels.hpp b/api/nblib/listed_forces/kernels.hpp
index 1f5b339ea77da58cffa0dc92b292955dde3f6c7c..5b98c1c32f842eabbfbbcb9d4d146e020a021b2f 100644 (file)
--- a/api/nblib/listed_forces/kernels.hpp
+++ b/api/nblib/listed_forces/kernels.hpp
@@ -426,6 +426,33 @@ static inline void makeAnglePeriodic(real& angle)
     }
 }
 
+/*! \brief calculate the cosine of the angle between a and b
+ *
+ * \tparam T  float or double
+ * \param a   a 3D vector
+ * \param b   another 3D vector
+ * \return    the cosine of the angle between a and b
+ *
+ *                  ax*bx + ay*by + az*bz
+ * cos-vec (a,b) =  ---------------------
+ *                      ||a|| * ||b||
+ */
+template<class T>
+inline T basicVectorCosAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
+{
+    gmx::BasicVector<double> a_double(a[0], a[1], a[2]);
+    gmx::BasicVector<double> b_double(b[0], b[1], b[2]);
+
+    double numerator     = dot(a_double, b_double);
+    double denominatorSq = dot(a_double, a_double) * dot(b_double, b_double);
+
+    T cosval = (denominatorSq > 0) ? static_cast<T>(numerator * gmx::invsqrt(denominatorSq)) : 1;
+    cosval   = std::min(cosval, T(1.0));
+
+    /* 25 TOTAL */
+    return std::max(cosval, T(-1.0));
+}
+
 //! \brief Computes and returns a dihedral phi angle
 static inline real dihedralPhi(rvec dxIJ, rvec dxKJ, rvec dxKL, rvec m, rvec n)
 {

diff --git a/api/nblib/listed_forces/tests/calculator.cpp b/api/nblib/listed_forces/tests/calculator.cpp
index 4093f00c2849a64ee1501dc74aa7e92904ebb211..4a56cca2731b5da01b69fc110e82d5d45ab979d3 100644 (file)
--- a/api/nblib/listed_forces/tests/calculator.cpp
+++ b/api/nblib/listed_forces/tests/calculator.cpp
@@ -119,7 +119,7 @@ protected:
         forces = std::vector<gmx::RVec>(3, gmx::RVec{ 0, 0, 0 });
 
         box.reset(new Box(3, 3, 3));
-        pbc.reset(new PbcHolder(*box));
+        pbc.reset(new PbcHolder(PbcType::Xyz, *box));
     }
 
     std::vector<gmx::RVec> x;
@@ -133,8 +133,9 @@ protected:
 
 TEST_F(ListedExampleData, ComputeHarmonicBondForces)
 {
-    auto indices = pickType<HarmonicBondType>(interactions).indices;
-    auto bonds   = pickType<HarmonicBondType>(interactions).parameters;
+    gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
+            pickType<HarmonicBondType>(interactions).indices;
+    gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
     computeForces(indices, bonds, x, &forces, *pbc);
 
     RefDataChecker vector3DTest(1e-3);
@@ -143,9 +144,10 @@ TEST_F(ListedExampleData, ComputeHarmonicBondForces)
 
 TEST_F(ListedExampleData, ComputeHarmonicBondEnergies)
 {
-    auto indices = pickType<HarmonicBondType>(interactions).indices;
-    auto bonds   = pickType<HarmonicBondType>(interactions).parameters;
-    real energy  = computeForces(indices, bonds, x, &forces, *pbc);
+    gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
+            pickType<HarmonicBondType>(interactions).indices;
+    gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
+    real                                  energy = computeForces(indices, bonds, x, &forces, *pbc);
 
     RefDataChecker vector3DTest(1e-4);
     vector3DTest.testReal(energy, "Bond energy");
@@ -153,8 +155,9 @@ TEST_F(ListedExampleData, ComputeHarmonicBondEnergies)
 
 TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
 {
-    auto indices = pickType<HarmonicAngle>(interactions).indices;
-    auto angles  = pickType<HarmonicAngle>(interactions).parameters;
+    gmx::ArrayRef<const InteractionIndex<HarmonicAngle>> indices =
+            pickType<HarmonicAngle>(interactions).indices;
+    gmx::ArrayRef<const HarmonicAngle> angles = pickType<HarmonicAngle>(interactions).parameters;
     computeForces(indices, angles, x, &forces, *pbc);
 
     RefDataChecker vector3DTest(1e-4);

diff --git a/api/nblib/listed_forces/tests/typetests.cpp b/api/nblib/listed_forces/tests/typetests.cpp
index 83d3610e75ff2610006a0d1e7164577f447ebe6b..15d7e35fdbe5386e55cdfcbb1f37806fb9917117 100644 (file)
--- a/api/nblib/listed_forces/tests/typetests.cpp
+++ b/api/nblib/listed_forces/tests/typetests.cpp
@@ -90,11 +90,15 @@ public:
         input_(std::move(input)),
         x_(std::move(coordinates)),
         indices_(std::move(indices)),
-        pbcHolder_(Box(1.5)),
+        pbcHolder_(PbcType::Xyz, Box(1.5)),
         checker_(refData_.rootChecker()),
         forces_(c_numAtoms, gmx::RVec{ 0, 0, 0 })
     {
-        energy_ = computeForces(indices_, input_, x_, &forces_, pbcHolder_);
+        energy_ = computeForces(gmx::ArrayRef<const InteractionIndex<Interaction>>(indices_),
+                                gmx::ArrayRef<const Interaction>(input_),
+                                x_,
+                                &forces_,
+                                pbcHolder_);
     }
 
     void checkForcesAndEnergies()

diff --git a/api/nblib/pbc.hpp b/api/nblib/pbc.hpp
index b52ab9eea3f0f08dce631eb387fe3effda2d3487..bc9523c14b0d018be87e93e01afb38b1c8da447c 100644 (file)
--- a/api/nblib/pbc.hpp
+++ b/api/nblib/pbc.hpp
@@ -34,8 +34,9 @@
  */
 /*! \internal \file
  * \brief
- * Implements a RAII encapsulation for PbcAiuc
- * note this is a header implementation to make inlining easier if needed
+ * Implements a RAII encapsulation objects for PbcAiuc and t_pbc that
+ * can be passed to listed forces.
+ * Note: this is a header implementation to make inlining easier if needed.
  *
  * \author Victor Holanda <victor.holanda@cscs.ch>
  * \author Joe Jordan <ejjordan@kth.se>
@@ -48,21 +49,23 @@
 
 #include "nblib/box.h"
 #include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pbcutil/pbc_aiuc.h"
 
 namespace nblib
 {
 
-/*! \brief life-time manager for PbcAiuc
+/*! \brief life-time manager for pbcAiuc
  *
  */
-class PbcHolder
+class PbcHolderAiuc
 {
 public:
     //! \brief ctor
-    explicit PbcHolder(const Box& box)
+    explicit PbcHolderAiuc(PbcType pbcType, const Box& box)
     {
-        setPbcAiuc(3, box.legacyMatrix(), &m_pbc);
+        setPbcAiuc(numPbcDimensions(pbcType), box.legacyMatrix(), &m_pbc);
     }
 
     //! \brief calculate pbc-aware r1-r2
@@ -75,6 +78,30 @@ private:
    PbcAiuc m_pbc;
 };
 
+/*! \brief life-time manager for t_pbc
+ *
+ */
+class PbcHolder
+{
+public:
+    //! \brief ctor
+    explicit PbcHolder(PbcType pbcType, const Box& box)
+    {
+        set_pbc(&m_pbc, pbcType, box.legacyMatrix());
+    }
+
+    //! \brief calculate pbc-aware r1-r2, including the shift index
+    template<class T>
+    inline int dxAiuc(const gmx::BasicVector<T>& r1, const gmx::BasicVector<T>& r2,
+                      gmx::BasicVector<T>& dr) const
+    {
+        return pbc_dx_aiuc(&m_pbc, r1.as_vec(), r2.as_vec(), dr.as_vec());
+    }
+
+private:
+    t_pbc m_pbc;
+};
+
 /*! \brief dummy class used to turn off Pbc in listed forces
  *
  */
@@ -83,9 +110,11 @@ class NoPbc
 public:
 
     //! \brief calculate r1-r2, ignore pbc
-    inline void dxAiuc(const rvec& r1, const rvec& r2, rvec dr) const
+    template<class T>
+    inline int dxAiuc(const gmx::BasicVector<T>& r1, const gmx::BasicVector<T>& r2, gmx::BasicVector<T>& dr) const
     {
-        rvec_sub(r1, r2, dr);
+        dr = r1 - r2;
+        return gmx::c_centralShiftIndex;
     }
 };
 

diff --git a/api/nblib/tests/pbcholder.cpp b/api/nblib/tests/pbcholder.cpp
index af53143aaa02188819530dd66ba7cfd2f39d16b4..db484bf94930b7a108aa42cd47c004f0fe6d41ac 100644 (file)
--- a/api/nblib/tests/pbcholder.cpp
+++ b/api/nblib/tests/pbcholder.cpp
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2020, by the GROMACS development team, led by
+ * Copyright (c) 2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -57,17 +57,17 @@ TEST(NBlibTest, PbcHolderWorks)
 {
     Box box(10, 10, 10);
 
-    PbcHolder pbcHolder(box);
+    PbcHolder pbcHolder(PbcType::Xyz, box);
 
-    rvec x1{ 1, 1, 1 }, x2{ 9, 9, 9 };
-    rvec dx;
+    gmx::RVec x1{ 1.0, 1.1, 0.9 }, x2{ 9, 8.9, 9.1 };
+    gmx::RVec dx;
 
     pbcHolder.dxAiuc(x1, x2, dx);
-    rvec ref{ 2, 2, 2 };
+    gmx::RVec ref{ 2, 2.2, 1.8 };
 
-    EXPECT_REAL_EQ_TOL(ref[0], dx[0], gmx::test::defaultRealTolerance());
-    EXPECT_REAL_EQ_TOL(ref[1], dx[1], gmx::test::defaultRealTolerance());
-    EXPECT_REAL_EQ_TOL(ref[2], dx[2], gmx::test::defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(ref[0], dx[0], gmx::test::relativeToleranceAsFloatingPoint(ref[0], 1e-6));
+    EXPECT_REAL_EQ_TOL(ref[1], dx[1], gmx::test::relativeToleranceAsFloatingPoint(ref[0], 1e-6));
+    EXPECT_REAL_EQ_TOL(ref[2], dx[2], gmx::test::relativeToleranceAsFloatingPoint(ref[0], 1e-6));
 }
 
 } // namespace nblib