api/nblib/listed_forces/calculator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,2021, 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 a bonded force calculator
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  * \author Artem Zhmurov <zhmurov@gmail.com>
  44  */
  45 #include <algorithm>
  46
  47 #include "nblib/box.h"
  48 #include "nblib/exception.h"
  49 #include "nblib/pbc.hpp"
  50 #include "nblib/listed_forces/calculator.h"
  51 #include "nblib/listed_forces/dataflow.hpp"
  52 #include "nblib/listed_forces/helpers.hpp"
  53
  54 namespace nblib
  55 {
  56
  57 ListedForceCalculator::~ListedForceCalculator() = default;
  58
  59 ListedForceCalculator::ListedForceCalculator(const ListedInteractionData& interactions,
  60                                              size_t                       bufferSize,
  61                                              int                          nthr,
  62                                              const Box&                   box) :
  63     numThreads(nthr),
  64     threadedForceBuffers_(numThreads),
  65     threadedShiftForceBuffers_(numThreads),
  66     pbcHolder_(std::make_unique<PbcHolder>(PbcType::Xyz, box))
  67 {
  68     // split up the work
  69     threadedInteractions_ = splitListedWork(interactions, bufferSize, numThreads);
  70
  71     // set up the reduction buffers
  72     int threadRange = (bufferSize + numThreads - 1) / numThreads;
  73 #pragma omp parallel for num_threads(numThreads) schedule(static)
  74     for (int i = 0; i < numThreads; ++i)
  75     {
  76         int rangeStart = i * threadRange;
  77         int rangeEnd   = rangeStart + threadRange;
  78         // the last range is a bit bigger due to integer division
  79         if (i == numThreads - 1)
  80         {
  81             rangeEnd = bufferSize;
  82         }
  83
  84         threadedForceBuffers_[i]      = ForceBufferProxy<Vec3>(rangeStart, rangeEnd);
  85         threadedShiftForceBuffers_[i] = std::vector<Vec3>(gmx::c_numShiftVectors);
  86     }
  87 }
  88
  89 template<class ShiftForce>
  90 void ListedForceCalculator::computeForcesAndEnergies(gmx::ArrayRef<const Vec3> x,
  91                                                      gmx::ArrayRef<Vec3>       forces,
  92                                                      [[maybe_unused]] gmx::ArrayRef<ShiftForce> shiftForces,
  93                                                      bool usePbc)
  94 {
  95     if (x.size() != forces.size())
  96     {
  97         throw InputException("Coordinates array and force buffer size mismatch");
  98     }
  99
 100     energyBuffer_.fill(0);
 101     std::vector<std::array<real, std::tuple_size<ListedInteractionData>::value>> energiesPerThread(numThreads);
 102
 103     constexpr bool haveShiftForces = !std::is_same_v<ShiftForce, std::nullptr_t>;
 104     if constexpr (haveShiftForces)
 105     {
 106         if (shiftForces.size() != gmx::c_numShiftVectors)
 107         {
 108             throw InputException("Shift vectors array size mismatch");
 109         }
 110     }
 111
 112 #pragma omp parallel for num_threads(numThreads) schedule(static)
 113     for (int thread = 0; thread < numThreads; ++thread)
 114     {
 115         std::conditional_t<haveShiftForces, gmx::ArrayRef<Vec3>, gmx::ArrayRef<std::nullptr_t>> shiftForceBuffer;
 116         if constexpr (haveShiftForces)
 117         {
 118             shiftForceBuffer = gmx::ArrayRef<Vec3>(threadedShiftForceBuffers_[thread]);
 119             std::fill(shiftForceBuffer.begin(), shiftForceBuffer.end(), Vec3{ 0, 0, 0 });
 120         }
 121
 122         ForceBufferProxy<Vec3>* threadBuffer = &threadedForceBuffers_[thread];
 123
 124         // forces in range of this thread are directly written into the output buffer
 125         threadBuffer->setMasterBuffer(forces);
 126
 127         // zero out the outliers in the thread buffer
 128         threadBuffer->clearOutliers();
 129
 130         if (usePbc)
 131         {
 132             energiesPerThread[thread] = reduceListedForces(
 133                     threadedInteractions_[thread], x, threadBuffer, shiftForceBuffer, *pbcHolder_);
 134         }
 135         else
 136         {
 137             energiesPerThread[thread] = reduceListedForces(
 138                     threadedInteractions_[thread], x, threadBuffer, shiftForceBuffer, NoPbc{});
 139         }
 140     }
 141
 142     // reduce shift forces
 143     // This is a potential candidate for OMP parallelization, but attention should be paid to the
 144     // relative costs of thread synchronization overhead vs reduction cost in contexts where the
 145     // number of threads could be large vs where number of threads could be small
 146     if constexpr (haveShiftForces)
 147     {
 148         for (int i = 0; i < gmx::c_numShiftVectors; ++i)
 149         {
 150             Vec3 threadSum{ 0, 0, 0 };
 151             for (int thread = 0; thread < numThreads; ++thread)
 152             {
 153                 threadSum += (threadedShiftForceBuffers_[thread])[i];
 154             }
 155             shiftForces[i] += threadSum;
 156         }
 157     }
 158
 159     // reduce energies
 160     for (int thread = 0; thread < numThreads; ++thread)
 161     {
 162         for (int type = 0; type < int(energyBuffer_.size()); ++type)
 163         {
 164             energyBuffer_[type] += energiesPerThread[thread][type];
 165         }
 166     }
 167     // reduce forces
 168 #pragma omp parallel for num_threads(numThreads) schedule(static)
 169     for (int thread = 0; thread < numThreads; ++thread)
 170     {
 171         auto& thisBuffer = threadedForceBuffers_[thread];
 172         // access outliers from other threads
 173         for (int otherThread = 0; otherThread < numThreads; ++otherThread)
 174         {
 175             auto& otherBuffer = threadedForceBuffers_[otherThread];
 176             for (const auto& outlier : otherBuffer)
 177             {
 178                 int index = outlier.first;
 179                 // check whether this outlier falls within the range of <thread>
 180                 if (thisBuffer.inRange(index))
 181                 {
 182                     auto force = outlier.second;
 183                     forces[index] += force;
 184                 }
 185             }
 186         }
 187     }
 188 }
 189
 190 void ListedForceCalculator::compute(gmx::ArrayRef<const Vec3> coordinates,
 191                                     gmx::ArrayRef<Vec3>       forces,
 192                                     gmx::ArrayRef<Vec3>       shiftForces,
 193                                     gmx::ArrayRef<real>       energies,
 194                                     bool                      usePbc)
 195 {
 196     computeForcesAndEnergies(coordinates, forces, shiftForces, usePbc);
 197     if (!energies.empty())
 198     {
 199         std::copy(energyBuffer_.begin(), energyBuffer_.end(), energies.begin());
 200     }
 201 }
 202
 203 void ListedForceCalculator::compute(gmx::ArrayRef<const Vec3> coordinates,
 204                                     gmx::ArrayRef<Vec3>       forces,
 205                                     gmx::ArrayRef<real>       energies,
 206                                     bool                      usePbc)
 207 {
 208     computeForcesAndEnergies(coordinates, forces, gmx::ArrayRef<std::nullptr_t>{}, usePbc);
 209     if (!energies.empty())
 210     {
 211         std::copy(energyBuffer_.begin(), energyBuffer_.end(), energies.begin());
 212     }
 213 }
 214
 215 } // namespace nblib