Improve doxygen in some nblib listed forces files
[alexxy/gromacs.git] / api / nblib / listed_forces / dataflow.hpp
index d078e176afea2fa665bc8ed58162efb92915bfd2..74b3bcca53c23f9efaff3e3a3f0297a407f08381 100644 (file)
@@ -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))
 namespace nblib
 {
 
-template<class TwoCenterType, class BasicVector>
+//! \brief returns the address of an element in the shiftForces buffer
+template<class ShiftForce>
+inline ShiftForce* accessShiftForces(int shiftIndex, gmx::ArrayRef<ShiftForce> shiftForces)
+{
+    return shiftForces.data() + shiftIndex;
+}
+
+//! \brief dummy op in case shift forces are not computed (will be removed by the compiler)
+inline std::nullptr_t accessShiftForces(int /* shiftIndex */,
+                                        gmx::ArrayRef<std::nullptr_t> /* shiftForces */)
+{
+    return nullptr;
+}
+
+template<class TwoCenterType, class BasicVector, class ShiftForce>
 inline NBLIB_ALWAYS_INLINE
-auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj)
+auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj,
+                      ShiftForce* sh_f, ShiftForce* sh_fc)
 {
     using ValueType = BasicVectorValueType_t<BasicVector>;
 
@@ -79,7 +93,7 @@ auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, Ba
     if (dr2 != 0.0)
     {
         force /= dr;
-        spreadTwoCenterForces(force, dx, fi, fj);
+        spreadTwoCenterForces(force, dx, fi, fj, sh_f, sh_fc);
     }
 
     return energy;
@@ -98,81 +112,166 @@ auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, Ba
  * \param[in] pbc Object used for computing distances accounting for PBC's
  * \return Computed kernel energies
  */
-template <class Buffer, class TwoCenterType, class BasicVector, class Pbc,
+template <class Buffer, class TwoCenterType, class BasicVector, class ShiftForce, 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,
+                         gmx::ArrayRef<ShiftForce> shiftForces,
                          const Pbc& pbc)
 {
     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 
     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   xi   = x[i];
+    BasicVector   xj   = x[j];
+    TwoCenterType bond = bondInstances[std::get<2>(index)];
 
-    gmx::RVec dx;
-    // calculate x1 - x2 modulo pbc
-    pbc.dxAiuc(x1, x2, dx);
+    BasicVector dx;
+    // calculate xi - xj modulo pbc
+    int sIdx = pbc.dxAiuc(xi, xj, dx);
 
-    energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j]);
+    ShiftForce* sh_f  = accessShiftForces(sIdx, shiftForces);
+    ShiftForce* sh_fc = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
+
+    energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j], sh_f, sh_fc);
     return energy;
 }
 
 
-template<class ThreeCenterType, class BasicVector>
+template<class ThreeCenterType, class BasicVector, class ShiftForce>
 inline NBLIB_ALWAYS_INLINE
 std::enable_if_t<HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
 addTwoCenterAggregate(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
-                      BasicVector* fi, BasicVector* fj, BasicVector* fk)
+                      BasicVector* fi, BasicVector* fj, BasicVector* fk,
+                      ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 {
-    if (parameters.manifest == ThreeCenterType::Cargo::ij)
+if (parameters.manifest == ThreeCenterType::Cargo::ij)
     {
         // i-j bond
-        return computeTwoCenter(parameters.twoCenter(), rij, fi, fj);
+        return computeTwoCenter(parameters.twoCenter(), rij, fi, fj, shf_ij, shf_c);
     }
     if (parameters.manifest == ThreeCenterType::Cargo::jk)
     {
         // j-k bond
-        return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj);
+        return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj, shf_kj, shf_c);
     }
 
     // aggregate is empty
     return 0.0;
 };
 
-template<class ThreeCenterType, class BasicVector>
+template<class ThreeCenterType, class BasicVector, class ShiftForce>
 inline NBLIB_ALWAYS_INLINE
 std::enable_if_t<!HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
-addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters,
-                      [[maybe_unused]] const BasicVector& rij,
-                      [[maybe_unused]] const BasicVector& rkj,
-                      [[maybe_unused]] BasicVector* fi,
-                      [[maybe_unused]] BasicVector* fj,
-                      [[maybe_unused]] BasicVector* fk)
+addTwoCenterAggregate(const ThreeCenterType& /* parameters */,
+                      const BasicVector& /* rij */,
+                      const BasicVector& /* rkj */,
+                      BasicVector* /* fi */,
+                      BasicVector* /* fj */,
+                      BasicVector* /* fk */,
+                      ShiftForce* /* shf_ij */,
+                      ShiftForce* /* shf_kj */,
+                      ShiftForce* /* shf_c */)
 {
     return 0.0;
 };
 
-template<class ThreeCenterType, class BasicVector>
+template<class ThreeCenterType, class BasicVector, class ShiftForce>
 inline NBLIB_ALWAYS_INLINE
 auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
-                        BasicVector* fi, BasicVector* fj, BasicVector* fk)
+                        const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+                        ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 {
     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.
+    // call type-specific angle kernel, e.g. harmonic, restricted, quartic, etc.
     auto [force, energy] = threeCenterKernel(theta, parameters);
 
-    spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk);
+    spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk, shf_ij, shf_kj, shf_c);
+
+    return energy;
+}
+
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const LinearAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
+                        const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+                        ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+    using ValueType = BasicVectorValueType_t<BasicVector>;
+
+    ValueType b  = parameters.equilConstant() - 1;
+    auto dr_vec  = b * rkj - parameters.equilConstant() * rij;
+    ValueType dr = norm(dr_vec);
+
+    auto [ci, ck, energy] = threeCenterKernel(dr, parameters);
+
+    BasicVector fi_loc = ci * dr_vec;
+    BasicVector fk_loc = ck * dr_vec;
+    BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+    *fi += fi_loc;
+    *fj += fj_loc;
+    *fk += fk_loc;
+
+    addShiftForce(fi_loc, shf_ij);
+    addShiftForce(fj_loc, shf_c);
+    addShiftForce(fk_loc, shf_kj);
+
+    return energy;
+}
+
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const CrossBondBond& parameters, const BasicVector& rij, const BasicVector& rkj,
+                        const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+                        ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+    using ValueType = BasicVectorValueType_t<BasicVector>;
+    // 28 flops from the norm() calls
+    auto [ci, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), parameters);
+
+    BasicVector fi_loc = ci * rij;
+    BasicVector fk_loc = ck * rkj;
+    BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+    *fi += fi_loc;
+    *fj += fj_loc;
+    *fk += fk_loc;
+
+    addShiftForce(fi_loc, shf_ij);
+    addShiftForce(fj_loc, shf_c);
+    addShiftForce(fk_loc, shf_kj);
+
+    return energy;
+}
+
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const CrossBondAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
+                        const BasicVector& rik, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+                        ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+    using ValueType = BasicVectorValueType_t<BasicVector>;
+    // 42 flops from the norm() calls
+    auto [ci, cj, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), norm(rik), parameters);
+
+    BasicVector fi_loc = ci * rij + ck * rik;
+    BasicVector fk_loc = cj * rkj - ck * rik;
+    BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+    *fi += fi_loc;
+    *fj += fj_loc;
+    *fk += fk_loc;
+
+    addShiftForce(fi_loc, shf_ij);
+    addShiftForce(fj_loc, shf_c);
+    addShiftForce(fk_loc, shf_kj);
 
     return energy;
 }
@@ -190,13 +289,14 @@ auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& ri
  * \param[in] PBC
  * \return Computed kernel energies
  */
-template <class Buffer, class ThreeCenterType, class BasicVector, class Pbc,
+template <class Buffer, class ThreeCenterType, class BasicVector, class ShiftForce, 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,
+                         gmx::ArrayRef<ShiftForce> shiftForces,
                          const Pbc& pbc)
 {
     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
@@ -205,19 +305,24 @@ 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;
+    int sIdx_ij = pbc.dxAiuc(xi, xj, rij); /* 3 */
+    int sIdx_kj = pbc.dxAiuc(xk, xj, rkj); /* 3 */
+    pbc.dxAiuc(xi, xk, rik); /* 3 */
 
-    energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
-    energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
+    ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
+    ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
+    ShiftForce* shf_c  = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
+
+    energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, rik, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
+    energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
 
     (*forces)[i] += fi;
     (*forces)[j] += fj;
@@ -251,14 +356,14 @@ addThreeCenterAggregate(const FourCenterType& parameters,
 template<class FourCenterType, class BasicVector>
 inline NBLIB_ALWAYS_INLINE
 std::enable_if_t<!HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
-addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
-                        [[maybe_unused]] const BasicVector& rij,
-                        [[maybe_unused]] const BasicVector& rkj,
-                        [[maybe_unused]] const BasicVector& rkl,
-                        [[maybe_unused]] BasicVector* fi,
-                        [[maybe_unused]] BasicVector* fj,
-                        [[maybe_unused]] BasicVector* fk,
-                        [[maybe_unused]] BasicVector* fl)
+addThreeCenterAggregate(const FourCenterType& /* parameters*/,
+                        const BasicVector& /* rij */,
+                        const BasicVector& /* rkj */,
+                        const BasicVector& /* rkl */,
+                        BasicVector* /* fi */,
+                        BasicVector* /* fj */,
+                        BasicVector* /* fk */,
+                        BasicVector* /* fl */)
 {
     return 0.0;
 };
@@ -276,45 +381,54 @@ addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
  * \param[in] pbc Object used for computing distances accounting for PBC's
  * \return Computed kernel energies
  */
-template <class Buffer, class FourCenterType, class BasicVector, class Pbc,
+template <class Buffer, class FourCenterType, class BasicVector, class ShiftForce, 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,
+                         gmx::ArrayRef<ShiftForce> shiftForces,
                          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;
-    pbc.dxAiuc(xi, xj, dxIJ);
-    pbc.dxAiuc(xk, xj, dxKJ);
+    BasicVector dxIJ, dxKJ, dxKL, dxLJ;
+    int sIdx_ij = pbc.dxAiuc(xi, xj, dxIJ);
+    int sIdx_kj = pbc.dxAiuc(xk, xj, dxKJ);
     pbc.dxAiuc(xk, xl, dxKL);
 
-    const FourCenterType& fourCenterTypeParams = parameters[std::get<4>(index)];
+    int sIdx_lj = pbc.dxAiuc(xl, xj, dxLJ);
+
+    ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
+    ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
+    ShiftForce* shf_lj = accessShiftForces(sIdx_lj, shiftForces);
+    ShiftForce* shf_c  = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
 
-    rvec m, n;
-    real phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
+    FourCenterType fourCenterTypeParams = parameters[std::get<4>(index)];
+
+    BasicVector m, n;
+    RealScalar phi = dihedralPhi(dxIJ, dxKJ, dxKL, &m, &n);
 
     auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
 
     energy.carrier()              = kernelEnergy;
     energy.threeCenterAggregate() = addThreeCenterAggregate(fourCenterTypeParams, dxIJ, dxKJ, dxKL, &fi, &fj, &fk, &fl);
 
-    spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl);
+    spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl, shf_ij, shf_kj, shf_lj, shf_c);
 
     (*forces)[i] += fi;
     (*forces)[j] += fj;
@@ -337,13 +451,14 @@ auto dispatchInteraction(const InteractionIndex<FourCenterType>& index,
  * \param[in] pbc Object used for computing distances accounting for PBC's
  * \return Computed kernel energies
  */
-template <class Buffer, class FiveCenterType, class BasicVector, class Pbc,
+template <class Buffer, class FiveCenterType, class BasicVector, class ShiftForce, 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,
+                         [[maybe_unused]] gmx::ArrayRef<ShiftForce> shiftForces,
                          const Pbc& pbc)
 {
     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
@@ -354,23 +469,24 @@ 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
-    // and we also don't want compiler warnings, so we cast to void
+    // we don't want to add [[maybe_unused]] in the signature since the params will
+    // be used once CMap is implemented, and we also don't want compiler warnings,
+    // so we cast to void.
     (void)fiveCenterTypeParams;
     (void)forces;
 
@@ -388,23 +504,35 @@ auto dispatchInteraction(const InteractionIndex<FiveCenterType>& index,
  * \param[in] pbc Object used for computing distances accounting for PBC's
  * \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,
+template <class Index, class InteractionType, class Buffer, class ShiftForce, class Pbc>
+auto computeForces(gmx::ArrayRef<const Index> indices,
+                   gmx::ArrayRef<const InteractionType> parameters,
                    gmx::ArrayRef<const Vec3> x,
                    Buffer* forces,
+                   gmx::ArrayRef<ShiftForce> shiftForces,
                    const Pbc& pbc)
 {
     KernelEnergy<BasicVectorValueType_t<Vec3>> energy;
 
     for (const auto& index : indices)
     {
-        energy += dispatchInteraction(index, interactionParameters, x, forces, pbc);
+        energy += dispatchInteraction(index, parameters, x, forces, shiftForces, pbc);
     }
 
     return energy;
 }
 
+//! \brief convenience overload without shift forces
+template <class Index, class InteractionType, class Buffer, class Pbc>
+auto computeForces(gmx::ArrayRef<const Index> indices,
+                   gmx::ArrayRef<const InteractionType> parameters,
+                   gmx::ArrayRef<const Vec3> x,
+                   Buffer* forces,
+                   const Pbc& pbc)
+{
+    return computeForces(indices, parameters, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
+}
+
 /*! \brief implement a loop over bond types and accumulate their force contributions
  *
  * \param[in] interactions interaction pairs and bond parameters
@@ -413,10 +541,11 @@ auto computeForces(const std::vector<Index>& indices,
  * \param[in] pbc Object used for computing distances accounting for PBC's
  * \return Computed kernel energies
  */
-template<class Buffer, class Pbc>
+template<class Buffer, class ShiftForce, class Pbc>
 auto reduceListedForces(const ListedInteractionData& interactions,
                         gmx::ArrayRef<const Vec3> x,
                         Buffer* forces,
+                        gmx::ArrayRef<ShiftForce> shiftForces,
                         const Pbc& pbc)
 {
     using ValueType = BasicVectorValueType_t<Vec3>;
@@ -424,10 +553,13 @@ auto reduceListedForces(const ListedInteractionData& interactions,
     energies.fill(0);
 
     // calculate one bond type
-    auto computeForceType = [forces, &x, &energies, &pbc](const auto& interactionElement) {
+    auto computeForceType = [forces, x, shiftForces, &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, shiftForces, pbc);
 
         energies[CarrierIndex<InteractionType>{}] += energy.carrier();
         energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();
@@ -441,6 +573,16 @@ auto reduceListedForces(const ListedInteractionData& interactions,
     return energies;
 }
 
+//! \brief convenience overload without shift forces
+template<class Buffer, class Pbc>
+auto reduceListedForces(const ListedInteractionData& interactions,
+                        gmx::ArrayRef<const Vec3> x,
+                        Buffer* forces,
+                        const Pbc& pbc)
+{
+    return reduceListedForces(interactions, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
+}
+
 } // namespace nblib
 
 #undef NBLIB_ALWAYS_INLINE