api/nblib/listed_forces/dataflow.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements the data flow from ListedInteractionData and coordinates
  38  * down to the individual interaction type kernels
  39  *
  40  * Intended for internal use inside ListedCalculator only.
  41  *
  42  * \author Victor Holanda <victor.holanda@cscs.ch>
  43  * \author Joe Jordan <ejjordan@kth.se>
  44  * \author Prashanth Kanduri <kanduri@cscs.ch>
  45  * \author Sebastian Keller <keller@cscs.ch>
  46  * \author Artem Zhmurov <zhmurov@gmail.com>
  47  */
  48
  49 #ifndef NBLIB_LISTEDFORCES_DATAFLOW_HPP
  50 #define NBLIB_LISTEDFORCES_DATAFLOW_HPP
  51
  52 #include <tuple>
  53 #include <vector>
  54
  55 #include "nblib/listed_forces/traits.h"
  56 #include "nblib/listed_forces/kernels.hpp"
  57 #include "nblib/util/util.hpp"
  58 #include "nblib/pbc.hpp"
  59 #include "nblib/vector.h"
  60 #include "gromacs/utility/arrayref.h"
  61
  62 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
  63
  64 namespace nblib
  65 {
  66
  67 //! \brief returns the address of an element in the shiftForces buffer
  68 template<class ShiftForce>
  69 inline ShiftForce* accessShiftForces(int shiftIndex, gmx::ArrayRef<ShiftForce> shiftForces)
  70 {
  71     return shiftForces.data() + shiftIndex;
  72 }
  73
  74 //! \brief dummy op in case shift forces are not computed (will be removed by the compiler)
  75 inline std::nullptr_t accessShiftForces(int /* shiftIndex */,
  76                                         gmx::ArrayRef<std::nullptr_t> /* shiftForces */)
  77 {
  78     return nullptr;
  79 }
  80
  81 template<class TwoCenterType, class BasicVector, class ShiftForce>
  82 inline NBLIB_ALWAYS_INLINE
  83 auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj,
  84                       ShiftForce* sh_f, ShiftForce* sh_fc)
  85 {
  86     using ValueType = BasicVectorValueType_t<BasicVector>;
  87
  88     ValueType dr2 = dot(dx, dx);
  89     ValueType dr  = std::sqrt(dr2);
  90     auto [force, energy] = bondKernel(dr, parameters);
  91
  92     // avoid division by 0
  93     if (dr2 != 0.0)
  94     {
  95         force /= dr;
  96         spreadTwoCenterForces(force, dx, fi, fj, sh_f, sh_fc);
  97     }
  98
  99     return energy;
 100 }
 101
 102 /*! \brief calculate two-center interactions
 103  *
 104  * \tparam Force buffer type
 105  * \tparam TwoCenterType The bond type to compute; used for type deduction
 106  * \tparam Cartesian vector type
 107  * \tparam PBC type
 108  * \param[in] index The atom and parameter indices used for computing the interaction
 109  * \param[in] bondInstances The full type-specific interaction list
 110  * \param[in] x The coordinates
 111  * \param[in/out] forces The forces
 112  * \param[in] pbc Object used for computing distances accounting for PBC's
 113  * \return Computed kernel energies
 114  */
 115 template <class Buffer, class TwoCenterType, class BasicVector, class ShiftForce, class Pbc,
 116           std::enable_if_t<Contains<TwoCenterType, SupportedTwoCenterTypes>{}>* = nullptr>
 117 inline NBLIB_ALWAYS_INLINE
 118 auto dispatchInteraction(InteractionIndex<TwoCenterType> index,
 119                          gmx::ArrayRef<const TwoCenterType> bondInstances,
 120                          gmx::ArrayRef<const BasicVector> x,
 121                          Buffer* forces,
 122                          gmx::ArrayRef<ShiftForce> shiftForces,
 123                          const Pbc& pbc)
 124 {
 125     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 126
 127     int i = std::get<0>(index);
 128     int j = std::get<1>(index);
 129     BasicVector   xi   = x[i];
 130     BasicVector   xj   = x[j];
 131     TwoCenterType bond = bondInstances[std::get<2>(index)];
 132
 133     BasicVector dx;
 134     // calculate xi - xj modulo pbc
 135     int sIdx = pbc.dxAiuc(xi, xj, dx);
 136
 137     ShiftForce* sh_f  = accessShiftForces(sIdx, shiftForces);
 138     ShiftForce* sh_fc = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
 139
 140     energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j], sh_f, sh_fc);
 141     return energy;
 142 }
 143
 144
 145 template<class ThreeCenterType, class BasicVector, class ShiftForce>
 146 inline NBLIB_ALWAYS_INLINE
 147 std::enable_if_t<HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
 148 addTwoCenterAggregate(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
 149                       BasicVector* fi, BasicVector* fj, BasicVector* fk,
 150                       ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 151 {
 152 if (parameters.manifest == ThreeCenterType::Cargo::ij)
 153     {
 154         // i-j bond
 155         return computeTwoCenter(parameters.twoCenter(), rij, fi, fj, shf_ij, shf_c);
 156     }
 157     if (parameters.manifest == ThreeCenterType::Cargo::jk)
 158     {
 159         // j-k bond
 160         return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj, shf_kj, shf_c);
 161     }
 162
 163     // aggregate is empty
 164     return 0.0;
 165 };
 166
 167 template<class ThreeCenterType, class BasicVector, class ShiftForce>
 168 inline NBLIB_ALWAYS_INLINE
 169 std::enable_if_t<!HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
 170 addTwoCenterAggregate(const ThreeCenterType& /* parameters */,
 171                       const BasicVector& /* rij */,
 172                       const BasicVector& /* rkj */,
 173                       BasicVector* /* fi */,
 174                       BasicVector* /* fj */,
 175                       BasicVector* /* fk */,
 176                       ShiftForce* /* shf_ij */,
 177                       ShiftForce* /* shf_kj */,
 178                       ShiftForce* /* shf_c */)
 179 {
 180     return 0.0;
 181 };
 182
 183 template<class ThreeCenterType, class BasicVector, class ShiftForce>
 184 inline NBLIB_ALWAYS_INLINE
 185 auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
 186                         const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
 187                         ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 188 {
 189     using ValueType = BasicVectorValueType_t<BasicVector>;
 190     // calculate 3-center common quantities: angle between x1-x2 and x2-x3
 191     // Todo: after sufficient evaluation, switch over to atan2 based algorithm
 192     ValueType costh = basicVectorCosAngle(rij, rkj); /* 25 */
 193     ValueType theta = std::acos(costh);    /* 10 */
 194
 195     // call type-specific angle kernel, e.g. harmonic, restricted, quartic, etc.
 196     auto [force, energy] = threeCenterKernel(theta, parameters);
 197
 198     spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk, shf_ij, shf_kj, shf_c);
 199
 200     return energy;
 201 }
 202
 203 template<class BasicVector, class ShiftForce>
 204 inline NBLIB_ALWAYS_INLINE
 205 auto computeThreeCenter(const LinearAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
 206                         const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
 207                         ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 208 {
 209     using ValueType = BasicVectorValueType_t<BasicVector>;
 210
 211     ValueType b  = parameters.equilConstant() - 1;
 212     auto dr_vec  = b * rkj - parameters.equilConstant() * rij;
 213     ValueType dr = norm(dr_vec);
 214
 215     auto [ci, ck, energy] = threeCenterKernel(dr, parameters);
 216
 217     BasicVector fi_loc = ci * dr_vec;
 218     BasicVector fk_loc = ck * dr_vec;
 219     BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
 220     *fi += fi_loc;
 221     *fj += fj_loc;
 222     *fk += fk_loc;
 223
 224     addShiftForce(fi_loc, shf_ij);
 225     addShiftForce(fj_loc, shf_c);
 226     addShiftForce(fk_loc, shf_kj);
 227
 228     return energy;
 229 }
 230
 231 template<class BasicVector, class ShiftForce>
 232 inline NBLIB_ALWAYS_INLINE
 233 auto computeThreeCenter(const CrossBondBond& parameters, const BasicVector& rij, const BasicVector& rkj,
 234                         const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
 235                         ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 236 {
 237     using ValueType = BasicVectorValueType_t<BasicVector>;
 238     // 28 flops from the norm() calls
 239     auto [ci, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), parameters);
 240
 241     BasicVector fi_loc = ci * rij;
 242     BasicVector fk_loc = ck * rkj;
 243     BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
 244     *fi += fi_loc;
 245     *fj += fj_loc;
 246     *fk += fk_loc;
 247
 248     addShiftForce(fi_loc, shf_ij);
 249     addShiftForce(fj_loc, shf_c);
 250     addShiftForce(fk_loc, shf_kj);
 251
 252     return energy;
 253 }
 254
 255 template<class BasicVector, class ShiftForce>
 256 inline NBLIB_ALWAYS_INLINE
 257 auto computeThreeCenter(const CrossBondAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
 258                         const BasicVector& rik, BasicVector* fi, BasicVector* fj, BasicVector* fk,
 259                         ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
 260 {
 261     using ValueType = BasicVectorValueType_t<BasicVector>;
 262     // 42 flops from the norm() calls
 263     auto [ci, cj, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), norm(rik), parameters);
 264
 265     BasicVector fi_loc = ci * rij + ck * rik;
 266     BasicVector fk_loc = cj * rkj - ck * rik;
 267     BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
 268     *fi += fi_loc;
 269     *fj += fj_loc;
 270     *fk += fk_loc;
 271
 272     addShiftForce(fi_loc, shf_ij);
 273     addShiftForce(fj_loc, shf_c);
 274     addShiftForce(fk_loc, shf_kj);
 275
 276     return energy;
 277 }
 278
 279 /*! \brief Calculate three-center interactions
 280  *
 281  * \tparam Force buffer type
 282  * \tparam Three centre interaction parameters
 283  * \tparam Cartesian vector type
 284  * \tparam PBC type
 285  * \param[in] index
 286  * \param[in] Bond parameters
 287  * \param[in] x coordinate array
 288  * \param[in/out] Force buffer
 289  * \param[in] PBC
 290  * \return Computed kernel energies
 291  */
 292 template <class Buffer, class ThreeCenterType, class BasicVector, class ShiftForce, class Pbc,
 293           std::enable_if_t<Contains<ThreeCenterType, SupportedThreeCenterTypes>{}>* = nullptr>
 294 inline NBLIB_ALWAYS_INLINE
 295 auto dispatchInteraction(InteractionIndex<ThreeCenterType> index,
 296                          gmx::ArrayRef<const ThreeCenterType> parameters,
 297                          gmx::ArrayRef<const BasicVector> x,
 298                          Buffer* forces,
 299                          gmx::ArrayRef<ShiftForce> shiftForces,
 300                          const Pbc& pbc)
 301 {
 302     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 303
 304     //! fetch input data: position vectors x1-x3 and interaction parameters
 305     int i = std::get<0>(index);
 306     int j = std::get<1>(index);
 307     int k = std::get<2>(index);
 308     BasicVector xi = x[i];
 309     BasicVector xj = x[j];
 310     BasicVector xk = x[k];
 311     ThreeCenterType threeCenterParameters = parameters[std::get<3>(index)];
 312
 313     BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
 314
 315     BasicVector rij, rkj, rik;
 316     int sIdx_ij = pbc.dxAiuc(xi, xj, rij); /* 3 */
 317     int sIdx_kj = pbc.dxAiuc(xk, xj, rkj); /* 3 */
 318     pbc.dxAiuc(xi, xk, rik); /* 3 */
 319
 320     ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
 321     ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
 322     ShiftForce* shf_c  = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
 323
 324     energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, rik, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
 325     energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
 326
 327     (*forces)[i] += fi;
 328     (*forces)[j] += fj;
 329     (*forces)[k] += fk;
 330
 331     return energy;
 332 }
 333
 334 template<class FourCenterType, class BasicVector>
 335 inline NBLIB_ALWAYS_INLINE
 336 std::enable_if_t<HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
 337 addThreeCenterAggregate(const FourCenterType& parameters,
 338                         const BasicVector& rij, const BasicVector& rkj, const BasicVector& rkl,
 339                         BasicVector* fi, BasicVector* fj, BasicVector* fk, BasicVector* fl)
 340 {
 341     using ValueType = BasicVectorValueType_t<BasicVector>;
 342     if (parameters.manifest == FourCenterType::Cargo::j)
 343     {
 344         return computeThreeCenter(parameters.threeCenter(), rij, rkj, fi, fj, fk);
 345     }
 346     if (parameters.manifest == FourCenterType::Cargo::k)
 347     {
 348         return computeThreeCenter(parameters.threeCenter(), ValueType(-1.0)*rkj, ValueType(-1.0)*rkl, fj, fk, fl);
 349         //return computeThreeCenter(parameters.threeCenter(), rkj, rkl, fj, fk, fl);
 350     }
 351
 352     // aggregate is empty
 353     return 0.0;
 354 };
 355
 356 template<class FourCenterType, class BasicVector>
 357 inline NBLIB_ALWAYS_INLINE
 358 std::enable_if_t<!HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
 359 addThreeCenterAggregate(const FourCenterType& /* parameters*/,
 360                         const BasicVector& /* rij */,
 361                         const BasicVector& /* rkj */,
 362                         const BasicVector& /* rkl */,
 363                         BasicVector* /* fi */,
 364                         BasicVector* /* fj */,
 365                         BasicVector* /* fk */,
 366                         BasicVector* /* fl */)
 367 {
 368     return 0.0;
 369 };
 370
 371 /*! \brief Calculate four-center interactions
 372  *
 373  * \tparam Force buffer type
 374  * \tparam FourCenterType The bond type to compute; used for type deduction
 375  * \tparam Cartesian vector type
 376  * \tparam PBC type
 377  * \param[in] index The atom and parameter indices used for computing the interaction
 378  * \param[in] parameters The full type-specific interaction list
 379  * \param[in] x The coordinates
 380  * \param[in/out] forces The forces
 381  * \param[in] pbc Object used for computing distances accounting for PBC's
 382  * \return Computed kernel energies
 383  */
 384 template <class Buffer, class FourCenterType, class BasicVector, class ShiftForce, class Pbc,
 385           std::enable_if_t<Contains<FourCenterType, SupportedFourCenterTypes>{}>* = nullptr>
 386 inline NBLIB_ALWAYS_INLINE
 387 auto dispatchInteraction(InteractionIndex<FourCenterType> index,
 388                          gmx::ArrayRef<const FourCenterType> parameters,
 389                          gmx::ArrayRef<const BasicVector> x,
 390                          Buffer* forces,
 391                          gmx::ArrayRef<ShiftForce> shiftForces,
 392                          const Pbc& pbc)
 393 {
 394     using RealScalar = BasicVectorValueType_t<BasicVector>;
 395     KernelEnergy<RealScalar> energy;
 396
 397     int i = std::get<0>(index);
 398     int j = std::get<1>(index);
 399     int k = std::get<2>(index);
 400     int l = std::get<3>(index);
 401
 402     BasicVector xi = x[i];
 403     BasicVector xj = x[j];
 404     BasicVector xk = x[k];
 405     BasicVector xl = x[l];
 406
 407     BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
 408
 409     BasicVector dxIJ, dxKJ, dxKL, dxLJ;
 410     int sIdx_ij = pbc.dxAiuc(xi, xj, dxIJ);
 411     int sIdx_kj = pbc.dxAiuc(xk, xj, dxKJ);
 412     pbc.dxAiuc(xk, xl, dxKL);
 413
 414     int sIdx_lj = pbc.dxAiuc(xl, xj, dxLJ);
 415
 416     ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
 417     ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
 418     ShiftForce* shf_lj = accessShiftForces(sIdx_lj, shiftForces);
 419     ShiftForce* shf_c  = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
 420
 421     FourCenterType fourCenterTypeParams = parameters[std::get<4>(index)];
 422
 423     BasicVector m, n;
 424     RealScalar phi = dihedralPhi(dxIJ, dxKJ, dxKL, &m, &n);
 425
 426     auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
 427
 428     energy.carrier()              = kernelEnergy;
 429     energy.threeCenterAggregate() = addThreeCenterAggregate(fourCenterTypeParams, dxIJ, dxKJ, dxKL, &fi, &fj, &fk, &fl);
 430
 431     spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl, shf_ij, shf_kj, shf_lj, shf_c);
 432
 433     (*forces)[i] += fi;
 434     (*forces)[j] += fj;
 435     (*forces)[k] += fk;
 436     (*forces)[l] += fl;
 437
 438     return energy;
 439 }
 440
 441 /*! \brief Calculate five-center interactions
 442  *
 443  * \tparam Force buffer type
 444  * \tparam FiveCenterType The bond type to compute; used for type deduction
 445  * \tparam Cartesian vector type
 446  * \tparam PBC type
 447  * \param[in] index The atom and parameter indices used for computing the interaction
 448  * \param[in] parameters The full type-specific interaction list
 449  * \param[in] x The coordinates
 450  * \param[in/out] forces The forces
 451  * \param[in] pbc Object used for computing distances accounting for PBC's
 452  * \return Computed kernel energies
 453  */
 454 template <class Buffer, class FiveCenterType, class BasicVector, class ShiftForce, class Pbc,
 455           std::enable_if_t<Contains<FiveCenterType, SupportedFiveCenterTypes>{}>* = nullptr>
 456 inline NBLIB_ALWAYS_INLINE
 457 auto dispatchInteraction(InteractionIndex<FiveCenterType> index,
 458                          gmx::ArrayRef<const FiveCenterType> parameters,
 459                          gmx::ArrayRef<const BasicVector> x,
 460                          Buffer* forces,
 461                          [[maybe_unused]] gmx::ArrayRef<ShiftForce> shiftForces,
 462                          const Pbc& pbc)
 463 {
 464     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 465
 466     int i = std::get<0>(index);
 467     int j = std::get<1>(index);
 468     int k = std::get<2>(index);
 469     int l = std::get<3>(index);
 470     int m = std::get<4>(index);
 471
 472     BasicVector xi = x[i];
 473     BasicVector xj = x[j];
 474     BasicVector xk = x[k];
 475     BasicVector xl = x[l];
 476     BasicVector xm = x[m];
 477
 478     BasicVector dxIJ, dxJK, dxKL, dxLM;
 479     pbc.dxAiuc(xi, xj, dxIJ);
 480     pbc.dxAiuc(xj, xk, dxJK);
 481     pbc.dxAiuc(xk, xl, dxKL);
 482     pbc.dxAiuc(xl, xm, dxLM);
 483
 484     FiveCenterType fiveCenterTypeParams = parameters[std::get<5>(index)];
 485
 486     // this dispatch function is not in use yet, because CMap is not yet implemented
 487     // we don't want to add [[maybe_unused]] in the signature since the params will
 488     // be used once CMap is implemented, and we also don't want compiler warnings,
 489     // so we cast to void.
 490     (void)fiveCenterTypeParams;
 491     (void)forces;
 492
 493     return energy;
 494 }
 495
 496
 497 /*! \brief implement a loop over bonds for a given BondType and Kernel
 498  *  corresponds to e.g. the "bonds" function at Gromacs:bonded.cpp@450
 499  *
 500  * \param[in] indices interaction atom pair indices + bond parameter index
 501  * \param[in] interactionParameters bond/interaction parameters
 502  * \param[in] x coordinate input
 503  * \param[in/out] forces The forces
 504  * \param[in] pbc Object used for computing distances accounting for PBC's
 505  * \return Computed kernel energies
 506  */
 507 template <class Index, class InteractionType, class Buffer, class ShiftForce, class Pbc>
 508 auto computeForces(gmx::ArrayRef<const Index> indices,
 509                    gmx::ArrayRef<const InteractionType> parameters,
 510                    gmx::ArrayRef<const Vec3> x,
 511                    Buffer* forces,
 512                    gmx::ArrayRef<ShiftForce> shiftForces,
 513                    const Pbc& pbc)
 514 {
 515     KernelEnergy<BasicVectorValueType_t<Vec3>> energy;
 516
 517     for (const auto& index : indices)
 518     {
 519         energy += dispatchInteraction(index, parameters, x, forces, shiftForces, pbc);
 520     }
 521
 522     return energy;
 523 }
 524
 525 //! \brief convenience overload without shift forces
 526 template <class Index, class InteractionType, class Buffer, class Pbc>
 527 auto computeForces(gmx::ArrayRef<const Index> indices,
 528                    gmx::ArrayRef<const InteractionType> parameters,
 529                    gmx::ArrayRef<const Vec3> x,
 530                    Buffer* forces,
 531                    const Pbc& pbc)
 532 {
 533     return computeForces(indices, parameters, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
 534 }
 535
 536 /*! \brief implement a loop over bond types and accumulate their force contributions
 537  *
 538  * \param[in] interactions interaction pairs and bond parameters
 539  * \param[in] x coordinate input
 540  * \param[in/out] forces output force buffer
 541  * \param[in] pbc Object used for computing distances accounting for PBC's
 542  * \return Computed kernel energies
 543  */
 544 template<class Buffer, class ShiftForce, class Pbc>
 545 auto reduceListedForces(const ListedInteractionData& interactions,
 546                         gmx::ArrayRef<const Vec3> x,
 547                         Buffer* forces,
 548                         gmx::ArrayRef<ShiftForce> shiftForces,
 549                         const Pbc& pbc)
 550 {
 551     using ValueType = BasicVectorValueType_t<Vec3>;
 552     std::array<ValueType, std::tuple_size<ListedInteractionData>::value> energies{0};
 553     energies.fill(0);
 554
 555     // calculate one bond type
 556     auto computeForceType = [forces, x, shiftForces, &energies, &pbc](const auto& interactionElement) {
 557         using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
 558
 559         gmx::ArrayRef<const InteractionIndex<InteractionType>> indices(interactionElement.indices);
 560         gmx::ArrayRef<const InteractionType>                   parameters(interactionElement.parameters);
 561
 562         KernelEnergy<ValueType> energy = computeForces(indices, parameters, x, forces, shiftForces, pbc);
 563
 564         energies[CarrierIndex<InteractionType>{}] += energy.carrier();
 565         energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();
 566         energies[ThreeCenterAggregateIndex<InteractionType>{}] += energy.threeCenterAggregate();
 567         // energies[FepIndex{}] += energy.freeEnergyDerivative();
 568     };
 569
 570     // calculate all bond types, returns a tuple with the energies for each type
 571     for_each_tuple(computeForceType, interactions);
 572
 573     return energies;
 574 }
 575
 576 //! \brief convenience overload without shift forces
 577 template<class Buffer, class Pbc>
 578 auto reduceListedForces(const ListedInteractionData& interactions,
 579                         gmx::ArrayRef<const Vec3> x,
 580                         Buffer* forces,
 581                         const Pbc& pbc)
 582 {
 583     return reduceListedForces(interactions, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
 584 }
 585
 586 } // namespace nblib
 587
 588 #undef NBLIB_ALWAYS_INLINE
 589
 590 #endif // NBLIB_LISTEDFORCES_DATAFLOW_HPP