api/nblib/gmxcalculatorcpu.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 Implements a force calculator based on GROMACS data structures.
  37  *
  38  * \author Victor Holanda <victor.holanda@cscs.ch>
  39  * \author Joe Jordan <ejjordan@kth.se>
  40  * \author Prashanth Kanduri <kanduri@cscs.ch>
  41  * \author Sebastian Keller <keller@cscs.ch>
  42  */
  43 #include "gromacs/ewald/ewald_utils.h"
  44 #include "gromacs/mdtypes/enerdata.h"
  45 #include "gromacs/mdtypes/interaction_const.h"
  46 #include "gromacs/nbnxm/atomdata.h"
  47 #include "gromacs/nbnxm/nbnxm.h"
  48 #include "gromacs/nbnxm/pairlistset.h"
  49 #include "gromacs/nbnxm/pairlistsets.h"
  50 #include "gromacs/nbnxm/pairsearch.h"
  51 #include "gromacs/utility/listoflists.h"
  52 #include "gromacs/utility/range.h"
  53 #include "nblib/exception.h"
  54 #include "nblib/gmxbackenddata.h"
  55 #include "nblib/gmxcalculatorcpu.h"
  56 #include "nblib/nbnxmsetuphelpers.h"
  57 #include "nblib/pbc.hpp"
  58 #include "nblib/systemdescription.h"
  59 #include "nblib/topology.h"
  60 #include "nblib/tpr.h"
  61 #include "nblib/virials.h"
  62
  63 namespace nblib
  64 {
  65
  66 class GmxNBForceCalculatorCpu::CpuImpl final
  67 {
  68 public:
  69     CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
  70             gmx::ArrayRef<real>    nonBondedParams,
  71             gmx::ArrayRef<real>    charges,
  72             gmx::ArrayRef<int64_t> particleInteractionFlags,
  73             gmx::ArrayRef<int>     exclusionRanges,
  74             gmx::ArrayRef<int>     exclusionElements,
  75             const NBKernelOptions& options);
  76
  77     //! calculates a new pair list based on new coordinates (for every NS step)
  78     void updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box);
  79
  80     //! Compute forces and return
  81     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  82                  const Box&                     box,
  83                  gmx::ArrayRef<gmx::RVec>       forceOutput);
  84
  85     //! Compute forces and virial tensor
  86     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  87                  const Box&                     box,
  88                  gmx::ArrayRef<gmx::RVec>       forceOutput,
  89                  gmx::ArrayRef<real>            virialOutput);
  90
  91     //! Compute forces, virial tensor and potential energies
  92     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  93                  const Box&                     box,
  94                  gmx::ArrayRef<gmx::RVec>       forceOutput,
  95                  gmx::ArrayRef<real>            virialOutput,
  96                  gmx::ArrayRef<real>            energyOutput);
  97
  98 private:
  99     //! \brief client-side provided system description data
 100     SystemDescription system_;
 101
 102     //! \brief Gmx backend objects, employed for calculating the forces
 103     GmxBackendData backend_;
 104 };
 105
 106 GmxNBForceCalculatorCpu::CpuImpl::CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
 107                                           gmx::ArrayRef<real>    nonBondedParams,
 108                                           gmx::ArrayRef<real>    charges,
 109                                           gmx::ArrayRef<int64_t> particleInteractionFlags,
 110                                           gmx::ArrayRef<int>     exclusionRanges,
 111                                           gmx::ArrayRef<int>     exclusionElements,
 112                                           const NBKernelOptions& options) :
 113     system_(SystemDescription(particleTypeIdOfAllParticles, nonBondedParams, charges, particleInteractionFlags)),
 114     backend_(GmxBackendData(options, findNumEnergyGroups(particleInteractionFlags), exclusionRanges, exclusionElements))
 115 {
 116     // Set up non-bonded verlet in the backend
 117     backend_.nbv_ = createNbnxmCPU(system_.numParticleTypes_,
 118                                    options,
 119                                    findNumEnergyGroups(particleInteractionFlags),
 120                                    system_.nonBondedParams_);
 121 }
 122
 123 void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates,
 124                                                       const Box&                     box)
 125 {
 126     if (coordinates.size() != system_.numParticles_)
 127     {
 128         throw InputException(
 129                 "Coordinate array containing different number of entries than particles in the "
 130                 "system");
 131     }
 132
 133     const auto* legacyBox = box.legacyMatrix();
 134     system_.box_          = box;
 135     updateForcerec(&backend_.forcerec_, box.legacyMatrix());
 136     if (TRICLINIC(legacyBox))
 137     {
 138         throw InputException("Only rectangular unit-cells are supported here");
 139     }
 140
 141     const rvec lowerCorner = { 0, 0, 0 };
 142     const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] };
 143
 144     const real particleDensity = static_cast<real>(coordinates.size()) / det(legacyBox);
 145
 146     // Put particles on a grid based on bounds specified by the box
 147     nbnxn_put_on_grid(backend_.nbv_.get(),
 148                       legacyBox,
 149                       0,
 150                       lowerCorner,
 151                       upperCorner,
 152                       nullptr,
 153                       { 0, int(coordinates.size()) },
 154                       particleDensity,
 155                       system_.particleInfo_,
 156                       coordinates,
 157                       0,
 158                       nullptr);
 159
 160     backend_.nbv_->constructPairlist(
 161             gmx::InteractionLocality::Local, backend_.exclusions_, 0, &backend_.nrnb_);
 162
 163     // Set Particle Types and Charges and VdW params
 164     backend_.nbv_->setAtomProperties(
 165             system_.particleTypeIdOfAllParticles_, system_.charges_, system_.particleInfo_);
 166     backend_.updatePairlistCalled = true;
 167 }
 168
 169 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 170                                                const Box&                     box,
 171                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 172                                                gmx::ArrayRef<real>            virialOutput,
 173                                                gmx::ArrayRef<real>            energyOutput)
 174 {
 175     if (coordinateInput.size() != forceOutput.size())
 176     {
 177         throw InputException("coordinate array and force buffer size mismatch");
 178     }
 179
 180     if (!backend_.updatePairlistCalled)
 181     {
 182         throw InputException("compute called without updating pairlist at least once");
 183     }
 184
 185     // update the box if changed
 186     if (!(system_.box_ == box))
 187     {
 188         system_.box_ = box;
 189         updateForcerec(&backend_.forcerec_, box.legacyMatrix());
 190     }
 191
 192     bool computeVirial               = !virialOutput.empty();
 193     bool computeEnergies             = !energyOutput.empty();
 194     backend_.stepWork_.computeVirial = computeVirial;
 195     backend_.stepWork_.computeEnergy = computeEnergies;
 196
 197     // update the coordinates in the backend
 198     backend_.nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput);
 199
 200     backend_.nbv_->dispatchNonbondedKernel(
 201             gmx::InteractionLocality::Local,
 202             backend_.interactionConst_,
 203             backend_.stepWork_,
 204             enbvClearFYes,
 205             backend_.forcerec_.shift_vec,
 206             backend_.enerd_.grpp.energyGroupPairTerms[backend_.forcerec_.haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR
 207                                                                                         : NonBondedEnergyTerms::LJSR],
 208             backend_.enerd_.grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
 209             &backend_.nrnb_);
 210
 211     backend_.nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput);
 212
 213     if (computeVirial)
 214     {
 215         // calculate shift forces and turn into an array ref
 216         std::vector<Vec3> shiftForcesVector(gmx::c_numShiftVectors, Vec3(0.0, 0.0, 0.0));
 217         nbnxn_atomdata_add_nbat_fshift_to_fshift(*backend_.nbv_->nbat, shiftForcesVector);
 218         auto shiftForcesRef = constArrayRefFromArray(shiftForcesVector.data(), shiftForcesVector.size());
 219
 220         std::vector<Vec3> shiftVectorsArray(gmx::c_numShiftVectors);
 221
 222         // copy shift vectors from ForceRec
 223         std::copy(backend_.forcerec_.shift_vec.begin(),
 224                   backend_.forcerec_.shift_vec.end(),
 225                   shiftVectorsArray.begin());
 226
 227         computeVirialTensor(
 228                 coordinateInput, forceOutput, shiftVectorsArray, shiftForcesRef, box, virialOutput);
 229     }
 230
 231     // extract term energies (per interaction type)
 232     if (computeEnergies)
 233     {
 234         int nGroupPairs = backend_.enerd_.grpp.nener;
 235         if (int(energyOutput.size()) != int(NonBondedEnergyTerms::Count) * nGroupPairs)
 236         {
 237             throw InputException("Array size for energy output is wrong\n");
 238         }
 239
 240         for (int eg = 0; eg < int(NonBondedEnergyTerms::Count); ++eg)
 241         {
 242             std::copy(begin(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 243                       end(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 244                       energyOutput.begin() + eg * nGroupPairs);
 245         }
 246     }
 247 }
 248
 249 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 250                                                const Box&                     box,
 251                                                gmx::ArrayRef<gmx::RVec>       forceOutput)
 252 {
 253     // compute forces and fill in force buffer
 254     compute(coordinateInput, box, forceOutput, gmx::ArrayRef<real>{}, gmx::ArrayRef<real>{});
 255 }
 256
 257 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 258                                                const Box&                     box,
 259                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 260                                                gmx::ArrayRef<real>            virialOutput)
 261 {
 262     // compute forces and fill in force buffer
 263     compute(coordinateInput, box, forceOutput, virialOutput, gmx::ArrayRef<real>{});
 264 }
 265
 266 GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef<int>  particleTypeIdOfAllParticles,
 267                                                  gmx::ArrayRef<real> nonBondedParams,
 268                                                  gmx::ArrayRef<real> charges,
 269                                                  gmx::ArrayRef<int64_t> particleInteractionFlags,
 270                                                  gmx::ArrayRef<int>     exclusionRanges,
 271                                                  gmx::ArrayRef<int>     exclusionElements,
 272                                                  const NBKernelOptions& options)
 273 {
 274     if (options.useGpu)
 275     {
 276         throw InputException("Use GmxNBForceCalculatorGpu for GPU support");
 277     }
 278
 279     impl_ = std::make_unique<CpuImpl>(particleTypeIdOfAllParticles,
 280                                       nonBondedParams,
 281                                       charges,
 282                                       particleInteractionFlags,
 283                                       exclusionRanges,
 284                                       exclusionElements,
 285                                       options);
 286 }
 287
 288 GmxNBForceCalculatorCpu::~GmxNBForceCalculatorCpu() = default;
 289
 290 //! calculates a new pair list based on new coordinates (for every NS step)
 291 void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box)
 292 {
 293     impl_->updatePairlist(coordinates, box);
 294 }
 295
 296 //! Compute forces and return
 297 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 298                                       const Box&                     box,
 299                                       gmx::ArrayRef<gmx::RVec>       forceOutput)
 300 {
 301     impl_->compute(coordinateInput, box, forceOutput);
 302 }
 303
 304 //! Compute forces and virial tensor
 305 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 306                                       const Box&                     box,
 307                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 308                                       gmx::ArrayRef<real>            virialOutput)
 309 {
 310     impl_->compute(coordinateInput, box, forceOutput, virialOutput);
 311 }
 312
 313 //! Compute forces, virial tensor and potential energies
 314 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 315                                       const Box&                     box,
 316                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 317                                       gmx::ArrayRef<real>            virialOutput,
 318                                       gmx::ArrayRef<real>            energyOutput)
 319 {
 320     impl_->compute(coordinateInput, box, forceOutput, virialOutput, energyOutput);
 321 }
 322
 323 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
 324                                                                     const NBKernelOptions& options)
 325 {
 326     std::vector<real> nonBondedParameters = createNonBondedParameters(
 327             topology.getParticleTypes(), topology.getNonBondedInteractionMap());
 328
 329     std::vector<int64_t> particleInteractionFlags = createParticleInfoAllVdw(topology.numParticles());
 330
 331     return std::make_unique<GmxNBForceCalculatorCpu>(topology.getParticleTypeIdOfAllParticles(),
 332                                                      nonBondedParameters,
 333                                                      topology.getCharges(),
 334                                                      particleInteractionFlags,
 335                                                      topology.exclusionLists().ListRanges,
 336                                                      topology.exclusionLists().ListElements,
 337                                                      options);
 338 }
 339
 340 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(TprReader& tprReader,
 341                                                                     const NBKernelOptions& options)
 342 {
 343     return std::make_unique<GmxNBForceCalculatorCpu>(tprReader.particleTypeIdOfAllParticles_,
 344                                                      tprReader.nonbondedParameters_,
 345                                                      tprReader.charges_,
 346                                                      tprReader.particleInteractionFlags_,
 347                                                      tprReader.exclusionListRanges_,
 348                                                      tprReader.exclusionListElements_,
 349                                                      options);
 350 }
 351
 352 } // namespace nblib