api/nblib/listed_forces/helpers.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, 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  * Helper data structures and utility functions for the nblib force calculator.
  38  * Intended for internal use.
  39  *
  40  * \author Victor Holanda <victor.holanda@cscs.ch>
  41  * \author Joe Jordan <ejjordan@kth.se>
  42  * \author Prashanth Kanduri <kanduri@cscs.ch>
  43  * \author Sebastian Keller <keller@cscs.ch>
  44  * \author Artem Zhmurov <zhmurov@gmail.com>
  45  */
  46
  47 #ifndef NBLIB_LISTEDFORCSES_HELPERS_HPP
  48 #define NBLIB_LISTEDFORCSES_HELPERS_HPP
  49
  50 #include <unordered_map>
  51
  52 #include "nblib/pbc.hpp"
  53 #include "definitions.h"
  54 #include "nblib/util/internal.h"
  55
  56 namespace nblib
  57 {
  58
  59 namespace detail
  60 {
  61 template<class T>
  62 inline void gmxRVecZeroWorkaround([[maybe_unused]] T& value)
  63 {
  64 }
  65
  66 template<>
  67 inline void gmxRVecZeroWorkaround<gmx::RVec>(gmx::RVec& value)
  68 {
  69     for (int i = 0; i < dimSize; ++i)
  70     {
  71         value[i] = 0;
  72     }
  73 }
  74 } // namespace detail
  75
  76 /*! \internal \brief object to store forces for multithreaded listed forces computation
  77  *
  78  */
  79 template<class T>
  80 class ForceBuffer
  81 {
  82     using HashMap = std::unordered_map<int, T>;
  83
  84 public:
  85     ForceBuffer() : rangeStart(0), rangeEnd(0) { }
  86
  87     ForceBuffer(T* mbuf, int rs, int re) :
  88         masterForceBuffer(mbuf),
  89         rangeStart(rs),
  90         rangeEnd(re)
  91     {
  92     }
  93
  94     void clear() { outliers.clear(); }
  95
  96     inline NBLIB_ALWAYS_INLINE T& operator[](int i)
  97     {
  98         if (i >= rangeStart && i < rangeEnd)
  99         {
 100             return masterForceBuffer[i];
 101         }
 102         else
 103         {
 104             if (outliers.count(i) == 0)
 105             {
 106                 T zero = T();
 107                 // if T = gmx::RVec, need to explicitly initialize it to zeros
 108                 detail::gmxRVecZeroWorkaround(zero);
 109                 outliers[i] = zero;
 110             }
 111             return outliers[i];
 112         }
 113     }
 114
 115     typename HashMap::const_iterator begin() { return outliers.begin(); }
 116     typename HashMap::const_iterator end() { return outliers.end(); }
 117
 118     [[nodiscard]] bool inRange(int index) const { return (index >= rangeStart && index < rangeEnd); }
 119
 120 private:
 121     T*  masterForceBuffer;
 122     int rangeStart;
 123     int rangeEnd;
 124
 125     HashMap outliers;
 126 };
 127
 128 namespace detail
 129 {
 130
 131 static int computeChunkIndex(int index, int totalRange, int nSplits)
 132 {
 133     if (totalRange < nSplits)
 134     {
 135         // if there's more threads than particles
 136         return index;
 137     }
 138
 139     int splitLength = totalRange / nSplits;
 140     return std::min(index / splitLength, nSplits - 1);
 141 }
 142
 143 } // namespace detail
 144
 145
 146 /*! \internal \brief splits an interaction tuple into nSplits interaction tuples
 147  *
 148  * \param interactions
 149  * \param totalRange the number of particle sequence coordinates
 150  * \param nSplits number to divide the total work by
 151  * \return
 152  */
 153 inline
 154 std::vector<ListedInteractionData> splitListedWork(const ListedInteractionData& interactions,
 155                                                    int                          totalRange,
 156                                                    int                          nSplits)
 157 {
 158     std::vector<ListedInteractionData> workDivision(nSplits);
 159
 160     auto splitOneElement = [totalRange, nSplits, &workDivision](const auto& inputElement) {
 161         // the index of inputElement in the ListedInteractionsTuple
 162         constexpr int elementIndex =
 163                 FindIndex<std::decay_t<decltype(inputElement)>, ListedInteractionData>{};
 164
 165         // for now, copy all parameters to each split
 166         // Todo: extract only the parameters needed for this split
 167         for (auto& workDivisionSplit : workDivision)
 168         {
 169             std::get<elementIndex>(workDivisionSplit).parameters = inputElement.parameters;
 170         }
 171
 172         // loop over all interactions in inputElement
 173         for (const auto& interactionIndex : inputElement.indices)
 174         {
 175             // each interaction has multiple coordinate indices
 176             // we must pick one of them to assign this interaction to one of the output index ranges
 177             // Todo: count indices outside the current split range in order to minimize the buffer size
 178             int representativeIndex =
 179                     *std::min_element(begin(interactionIndex), end(interactionIndex) - 1);
 180             int splitIndex = detail::computeChunkIndex(representativeIndex, totalRange, nSplits);
 181
 182             std::get<elementIndex>(workDivision[splitIndex]).indices.push_back(interactionIndex);
 183         }
 184     };
 185
 186     // split each interaction type in the input interaction tuple
 187     for_each_tuple(splitOneElement, interactions);
 188
 189     return workDivision;
 190 }
 191
 192 } // namespace nblib
 193
 194 #endif // NBLIB_LISTEDFORCSES_HELPERS_HPP