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/math/vec.h"
  61 #include "gromacs/utility/arrayref.h"
  62
  63 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
  64
  65 namespace nblib
  66 {
  67
  68 template<class TwoCenterType, class BasicVector>
  69 inline NBLIB_ALWAYS_INLINE
  70 auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj)
  71 {
  72     using ValueType = BasicVectorValueType_t<BasicVector>;
  73
  74     ValueType dr2 = dot(dx, dx);
  75     ValueType dr  = std::sqrt(dr2);
  76     auto [force, energy] = bondKernel(dr, parameters);
  77
  78     // avoid division by 0
  79     if (dr2 != 0.0)
  80     {
  81         force /= dr;
  82         spreadTwoCenterForces(force, dx, fi, fj);
  83     }
  84
  85     return energy;
  86 }
  87
  88 /*! \brief calculate two-center interactions
  89  *
  90  * \tparam Force buffer type
  91  * \tparam TwoCenterType The bond type to compute; used for type deduction
  92  * \tparam Cartesian vector type
  93  * \tparam PBC type
  94  * \param[in] index The atom and parameter indices used for computing the interaction
  95  * \param[in] bondInstances The full type-specific interaction list
  96  * \param[in] x The coordinates
  97  * \param[in/out] forces The forces
  98  * \param[in] pbc Object used for computing distances accounting for PBC's
  99  * \return Computed kernel energies
 100  */
 101 template <class Buffer, class TwoCenterType, class BasicVector, class Pbc,
 102           std::enable_if_t<Contains<TwoCenterType, SupportedTwoCenterTypes>{}>* = nullptr>
 103 inline NBLIB_ALWAYS_INLINE
 104 auto dispatchInteraction(const InteractionIndex<TwoCenterType>& index,
 105                          const std::vector<TwoCenterType>& bondInstances,
 106                          gmx::ArrayRef<const BasicVector> x,
 107                          Buffer* forces,
 108                          const Pbc& pbc)
 109 {
 110     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 111
 112     int i = std::get<0>(index);
 113     int j = std::get<1>(index);
 114     const gmx::RVec& x1 = x[i];
 115     const gmx::RVec& x2 = x[j];
 116     const TwoCenterType& bond = bondInstances[std::get<2>(index)];
 117
 118     gmx::RVec dx;
 119     // calculate x1 - x2 modulo pbc
 120     pbc.dxAiuc(x1, x2, dx);
 121
 122     energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j]);
 123     return energy;
 124 }
 125
 126
 127 template<class ThreeCenterType, class BasicVector>
 128 inline NBLIB_ALWAYS_INLINE
 129 std::enable_if_t<HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
 130 addTwoCenterAggregate(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
 131                       BasicVector* fi, BasicVector* fj, BasicVector* fk)
 132 {
 133     if (parameters.manifest == ThreeCenterType::Cargo::ij)
 134     {
 135         // i-j bond
 136         return computeTwoCenter(parameters.twoCenter(), rij, fi, fj);
 137     }
 138     if (parameters.manifest == ThreeCenterType::Cargo::jk)
 139     {
 140         // j-k bond
 141         return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj);
 142     }
 143
 144     // aggregate is empty
 145     return 0.0;
 146 };
 147
 148 template<class ThreeCenterType, class BasicVector>
 149 inline NBLIB_ALWAYS_INLINE
 150 std::enable_if_t<!HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
 151 addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters,
 152                       [[maybe_unused]] const BasicVector& rij,
 153                       [[maybe_unused]] const BasicVector& rkj,
 154                       [[maybe_unused]] BasicVector* fi,
 155                       [[maybe_unused]] BasicVector* fj,
 156                       [[maybe_unused]] BasicVector* fk)
 157 {
 158     return 0.0;
 159 };
 160
 161 template<class ThreeCenterType, class BasicVector>
 162 inline NBLIB_ALWAYS_INLINE
 163 auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
 164                         BasicVector* fi, BasicVector* fj, BasicVector* fk)
 165 {
 166     using ValueType = BasicVectorValueType_t<BasicVector>;
 167     // calculate 3-center common quantities: angle between x1-x2 and x2-x3
 168     // Todo: after sufficient evaluation, switch over to atan2 based algorithm
 169     ValueType costh = cos_angle(rij, rkj); /* 25             */
 170     ValueType theta = std::acos(costh);    /* 10             */
 171
 172     // call type-specific angle kernel, e.g. harmonic, linear, quartic,  etc.
 173     auto [force, energy] = threeCenterKernel(theta, parameters);
 174
 175     spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk);
 176
 177     return energy;
 178 }
 179
 180 /*! \brief Calculate three-center interactions
 181  *
 182  * \tparam Force buffer type
 183  * \tparam Three centre interaction parameters
 184  * \tparam Cartesian vector type
 185  * \tparam PBC type
 186  * \param[in] index
 187  * \param[in] Bond parameters
 188  * \param[in] x coordinate array
 189  * \param[in/out] Force buffer
 190  * \param[in] PBC
 191  * \return Computed kernel energies
 192  */
 193 template <class Buffer, class ThreeCenterType, class BasicVector, class Pbc,
 194           std::enable_if_t<Contains<ThreeCenterType, SupportedThreeCenterTypes>{}>* = nullptr>
 195 inline NBLIB_ALWAYS_INLINE
 196 auto dispatchInteraction(const InteractionIndex<ThreeCenterType>& index,
 197                          const std::vector<ThreeCenterType>& parameters,
 198                          gmx::ArrayRef<const BasicVector> x,
 199                          Buffer* forces,
 200                          const Pbc& pbc)
 201 {
 202     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 203
 204     //! fetch input data: position vectors x1-x3 and interaction parameters
 205     int i = std::get<0>(index);
 206     int j = std::get<1>(index);
 207     int k = std::get<2>(index);
 208     const gmx::RVec& xi = x[i];
 209     const gmx::RVec& xj = x[j];
 210     const gmx::RVec& xk = x[k];
 211     const ThreeCenterType& threeCenterParameters = parameters[std::get<3>(index)];
 212
 213     gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
 214
 215     gmx::RVec rij, rkj;
 216     pbc.dxAiuc(xi, xj, rij); /*  3              */
 217     pbc.dxAiuc(xk, xj, rkj); /*  3              */
 218
 219     energy.carrier()            = computeThreeCenter(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
 220     energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
 221
 222     (*forces)[i] += fi;
 223     (*forces)[j] += fj;
 224     (*forces)[k] += fk;
 225
 226     return energy;
 227 }
 228
 229 template<class FourCenterType, class BasicVector>
 230 inline NBLIB_ALWAYS_INLINE
 231 std::enable_if_t<HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
 232 addThreeCenterAggregate(const FourCenterType& parameters,
 233                         const BasicVector& rij, const BasicVector& rkj, const BasicVector& rkl,
 234                         BasicVector* fi, BasicVector* fj, BasicVector* fk, BasicVector* fl)
 235 {
 236     using ValueType = BasicVectorValueType_t<BasicVector>;
 237     if (parameters.manifest == FourCenterType::Cargo::j)
 238     {
 239         return computeThreeCenter(parameters.threeCenter(), rij, rkj, fi, fj, fk);
 240     }
 241     if (parameters.manifest == FourCenterType::Cargo::k)
 242     {
 243         return computeThreeCenter(parameters.threeCenter(), ValueType(-1.0)*rkj, ValueType(-1.0)*rkl, fj, fk, fl);
 244         //return computeThreeCenter(parameters.threeCenter(), rkj, rkl, fj, fk, fl);
 245     }
 246
 247     // aggregate is empty
 248     return 0.0;
 249 };
 250
 251 template<class FourCenterType, class BasicVector>
 252 inline NBLIB_ALWAYS_INLINE
 253 std::enable_if_t<!HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
 254 addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
 255                         [[maybe_unused]] const BasicVector& rij,
 256                         [[maybe_unused]] const BasicVector& rkj,
 257                         [[maybe_unused]] const BasicVector& rkl,
 258                         [[maybe_unused]] BasicVector* fi,
 259                         [[maybe_unused]] BasicVector* fj,
 260                         [[maybe_unused]] BasicVector* fk,
 261                         [[maybe_unused]] BasicVector* fl)
 262 {
 263     return 0.0;
 264 };
 265
 266 /*! \brief Calculate four-center interactions
 267  *
 268  * \tparam Force buffer type
 269  * \tparam FourCenterType The bond type to compute; used for type deduction
 270  * \tparam Cartesian vector type
 271  * \tparam PBC type
 272  * \param[in] index The atom and parameter indices used for computing the interaction
 273  * \param[in] parameters The full type-specific interaction list
 274  * \param[in] x The coordinates
 275  * \param[in/out] forces The forces
 276  * \param[in] pbc Object used for computing distances accounting for PBC's
 277  * \return Computed kernel energies
 278  */
 279 template <class Buffer, class FourCenterType, class BasicVector, class Pbc,
 280           std::enable_if_t<Contains<FourCenterType, SupportedFourCenterTypes>{}>* = nullptr>
 281 inline NBLIB_ALWAYS_INLINE
 282 auto dispatchInteraction(const InteractionIndex<FourCenterType>& index,
 283                          const std::vector<FourCenterType>& parameters,
 284                          gmx::ArrayRef<const BasicVector> x,
 285                          Buffer* forces,
 286                          const Pbc& pbc)
 287 {
 288     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 289
 290     int i = std::get<0>(index);
 291     int j = std::get<1>(index);
 292     int k = std::get<2>(index);
 293     int l = std::get<3>(index);
 294
 295     const gmx::RVec& xi = x[i];
 296     const gmx::RVec& xj = x[j];
 297     const gmx::RVec& xk = x[k];
 298     const gmx::RVec& xl = x[l];
 299
 300     gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
 301
 302     gmx::RVec dxIJ, dxKJ, dxKL;
 303     pbc.dxAiuc(xi, xj, dxIJ);
 304     pbc.dxAiuc(xk, xj, dxKJ);
 305     pbc.dxAiuc(xk, xl, dxKL);
 306
 307     const FourCenterType& fourCenterTypeParams = parameters[std::get<4>(index)];
 308
 309     rvec m, n;
 310     real phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
 311
 312     auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
 313
 314     energy.carrier()              = kernelEnergy;
 315     energy.threeCenterAggregate() = addThreeCenterAggregate(fourCenterTypeParams, dxIJ, dxKJ, dxKL, &fi, &fj, &fk, &fl);
 316
 317     spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl);
 318
 319     (*forces)[i] += fi;
 320     (*forces)[j] += fj;
 321     (*forces)[k] += fk;
 322     (*forces)[l] += fl;
 323
 324     return energy;
 325 }
 326
 327 /*! \brief Calculate five-center interactions
 328  *
 329  * \tparam Force buffer type
 330  * \tparam FiveCenterType The bond type to compute; used for type deduction
 331  * \tparam Cartesian vector type
 332  * \tparam PBC type
 333  * \param[in] index The atom and parameter indices used for computing the interaction
 334  * \param[in] parameters The full type-specific interaction list
 335  * \param[in] x The coordinates
 336  * \param[in/out] forces The forces
 337  * \param[in] pbc Object used for computing distances accounting for PBC's
 338  * \return Computed kernel energies
 339  */
 340 template <class Buffer, class FiveCenterType, class BasicVector, class Pbc,
 341           std::enable_if_t<Contains<FiveCenterType, SupportedFiveCenterTypes>{}>* = nullptr>
 342 inline NBLIB_ALWAYS_INLINE
 343 auto dispatchInteraction(const InteractionIndex<FiveCenterType>& index,
 344                          const std::vector<FiveCenterType>& parameters,
 345                          gmx::ArrayRef<const BasicVector> x,
 346                          Buffer* forces,
 347                          const Pbc& pbc)
 348 {
 349     KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
 350
 351     int i = std::get<0>(index);
 352     int j = std::get<1>(index);
 353     int k = std::get<2>(index);
 354     int l = std::get<3>(index);
 355     int m = std::get<4>(index);
 356
 357     const gmx::RVec& xi = x[i];
 358     const gmx::RVec& xj = x[j];
 359     const gmx::RVec& xk = x[k];
 360     const gmx::RVec& xl = x[l];
 361     const gmx::RVec& xm = x[m];
 362
 363     gmx::RVec dxIJ, dxJK, dxKL, dxLM;
 364     pbc.dxAiuc(xi, xj, dxIJ);
 365     pbc.dxAiuc(xj, xk, dxJK);
 366     pbc.dxAiuc(xk, xl, dxKL);
 367     pbc.dxAiuc(xl, xm, dxLM);
 368
 369     const FiveCenterType& fiveCenterTypeParams = parameters[std::get<5>(index)];
 370
 371     // this dispatch function is not in use yet, because CMap is not yet implemented
 372     // we don't want to add [[maybe_unused]] in the signature
 373     // and we also don't want compiler warnings, so we cast to void
 374     (void)fiveCenterTypeParams;
 375     (void)forces;
 376
 377     return energy;
 378 }
 379
 380
 381 /*! \brief implement a loop over bonds for a given BondType and Kernel
 382  *  corresponds to e.g. the "bonds" function at Gromacs:bonded.cpp@450
 383  *
 384  * \param[in] indices interaction atom pair indices + bond parameter index
 385  * \param[in] interactionParameters bond/interaction parameters
 386  * \param[in] x coordinate input
 387  * \param[in/out] forces The forces
 388  * \param[in] pbc Object used for computing distances accounting for PBC's
 389  * \return Computed kernel energies
 390  */
 391 template <class Index, class InteractionType, class Buffer, class Pbc>
 392 auto computeForces(const std::vector<Index>& indices,
 393                    const std::vector<InteractionType>& interactionParameters,
 394                    gmx::ArrayRef<const Vec3> x,
 395                    Buffer* forces,
 396                    const Pbc& pbc)
 397 {
 398     KernelEnergy<BasicVectorValueType_t<Vec3>> energy;
 399
 400     for (const auto& index : indices)
 401     {
 402         energy += dispatchInteraction(index, interactionParameters, x, forces, pbc);
 403     }
 404
 405     return energy;
 406 }
 407
 408 /*! \brief implement a loop over bond types and accumulate their force contributions
 409  *
 410  * \param[in] interactions interaction pairs and bond parameters
 411  * \param[in] x coordinate input
 412  * \param[in/out] forces output force buffer
 413  * \param[in] pbc Object used for computing distances accounting for PBC's
 414  * \return Computed kernel energies
 415  */
 416 template<class Buffer, class Pbc>
 417 auto reduceListedForces(const ListedInteractionData& interactions,
 418                         gmx::ArrayRef<const Vec3> x,
 419                         Buffer* forces,
 420                         const Pbc& pbc)
 421 {
 422     using ValueType = BasicVectorValueType_t<Vec3>;
 423     std::array<ValueType, std::tuple_size<ListedInteractionData>::value> energies{0};
 424     energies.fill(0);
 425
 426     // calculate one bond type
 427     auto computeForceType = [forces, &x, &energies, &pbc](const auto& interactionElement) {
 428         using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
 429
 430         KernelEnergy<ValueType> energy = computeForces(interactionElement.indices, interactionElement.parameters, x, forces, pbc);
 431
 432         energies[CarrierIndex<InteractionType>{}] += energy.carrier();
 433         energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();
 434         energies[ThreeCenterAggregateIndex<InteractionType>{}] += energy.threeCenterAggregate();
 435         // energies[FepIndex{}] += energy.freeEnergyDerivative();
 436     };
 437
 438     // calculate all bond types, returns a tuple with the energies for each type
 439     for_each_tuple(computeForceType, interactions);
 440
 441     return energies;
 442 }
 443
 444 } // namespace nblib
 445
 446 #undef NBLIB_ALWAYS_INLINE
 447
 448 #endif // NBLIB_LISTEDFORCES_DATAFLOW_HPP