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