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