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 "gromacs/utility/arrayref.h"
  53
  54 #include "nblib/pbc.hpp"
  55 #include "definitions.h"
  56 #include "nblib/util/util.hpp"
  57
  58 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
  59
  60 namespace nblib
  61 {
  62
  63 namespace detail
  64 {
  65 template<class T>
  66 inline void gmxRVecZeroWorkaround([[maybe_unused]] T& value)
  67 {
  68 }
  69
  70 template<>
  71 inline void gmxRVecZeroWorkaround<gmx::RVec>(gmx::RVec& value)
  72 {
  73     for (int i = 0; i < dimSize; ++i)
  74     {
  75         value[i] = 0;
  76     }
  77 }
  78 } // namespace detail
  79
  80 /*! \internal \brief proxy object to access forces in an underlying buffer
  81  *
  82  * Depending on the index, either the underlying master buffer, or local
  83  * storage for outliers is accessed. This object does not own the master buffer.
  84  *
  85  */
  86 template<class T>
  87 class ForceBufferProxy
  88 {
  89     using HashMap = std::unordered_map<int, T>;
  90
  91 public:
  92     ForceBufferProxy() : rangeStart_(0), rangeEnd_(0) { }
  93
  94     ForceBufferProxy(int rangeStart, int rangeEnd) : rangeStart_(rangeStart), rangeEnd_(rangeEnd)
  95     {
  96     }
  97
  98     void clearOutliers() { outliers.clear(); }
  99
 100     inline NBLIB_ALWAYS_INLINE T& operator[](int i)
 101     {
 102         if (i >= rangeStart_ && i < rangeEnd_)
 103         {
 104             return masterForceBuffer[i];
 105         }
 106         else
 107         {
 108             if (outliers.count(i) == 0)
 109             {
 110                 T zero = T();
 111                 // if T = gmx::RVec, need to explicitly initialize it to zeros
 112                 detail::gmxRVecZeroWorkaround(zero);
 113                 outliers[i] = zero;
 114             }
 115             return outliers[i];
 116         }
 117     }
 118
 119     typename HashMap::const_iterator begin() { return outliers.begin(); }
 120     typename HashMap::const_iterator end() { return outliers.end(); }
 121
 122     [[nodiscard]] bool inRange(int index) const { return (index >= rangeStart_ && index < rangeEnd_); }
 123
 124     void setMasterBuffer(gmx::ArrayRef<T> buffer) { masterForceBuffer = buffer; }
 125
 126 private:
 127     gmx::ArrayRef<T> masterForceBuffer;
 128     int rangeStart_;
 129     int rangeEnd_;
 130
 131     HashMap outliers;
 132 };
 133
 134 namespace detail
 135 {
 136
 137 static int computeChunkIndex(int index, int totalRange, int nSplits)
 138 {
 139     if (totalRange < nSplits)
 140     {
 141         // if there's more threads than particles
 142         return index;
 143     }
 144
 145     int splitLength = totalRange / nSplits;
 146     return std::min(index / splitLength, nSplits - 1);
 147 }
 148
 149 } // namespace detail
 150
 151
 152 /*! \internal \brief splits an interaction tuple into nSplits interaction tuples
 153  *
 154  * \param interactions
 155  * \param totalRange the number of particle sequence coordinates
 156  * \param nSplits number to divide the total work by
 157  * \return
 158  */
 159 inline
 160 std::vector<ListedInteractionData> splitListedWork(const ListedInteractionData& interactions,
 161                                                    int                          totalRange,
 162                                                    int                          nSplits)
 163 {
 164     std::vector<ListedInteractionData> workDivision(nSplits);
 165
 166     auto splitOneElement = [totalRange, nSplits, &workDivision](const auto& inputElement) {
 167         // the index of inputElement in the ListedInteractionsTuple
 168         constexpr int elementIndex =
 169                 FindIndex<std::decay_t<decltype(inputElement)>, ListedInteractionData>{};
 170
 171         // for now, copy all parameters to each split
 172         // Todo: extract only the parameters needed for this split
 173         for (auto& workDivisionSplit : workDivision)
 174         {
 175             std::get<elementIndex>(workDivisionSplit).parameters = inputElement.parameters;
 176         }
 177
 178         // loop over all interactions in inputElement
 179         for (const auto& interactionIndex : inputElement.indices)
 180         {
 181             // each interaction has multiple coordinate indices
 182             // we must pick one of them to assign this interaction to one of the output index ranges
 183             // Todo: count indices outside the current split range in order to minimize the buffer size
 184             int representativeIndex =
 185                     *std::min_element(begin(interactionIndex), end(interactionIndex) - 1);
 186             int splitIndex = detail::computeChunkIndex(representativeIndex, totalRange, nSplits);
 187
 188             std::get<elementIndex>(workDivision[splitIndex]).indices.push_back(interactionIndex);
 189         }
 190     };
 191
 192     // split each interaction type in the input interaction tuple
 193     for_each_tuple(splitOneElement, interactions);
 194
 195     return workDivision;
 196 }
 197
 198 } // namespace nblib
 199
 200 #undef NBLIB_ALWAYS_INLINE
 201
 202 #endif // NBLIB_LISTEDFORCSES_HELPERS_HPP