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<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<gmx::RVec> coordinates, const Box& box)
 124 {
 125     if (coordinates.size() != system_.numParticles_)
 126     {
 127         throw InputException(
 128                 "Coordinate array containing different number of entries than particles in the "
 129                 "system");
 130     }
 131
 132     const auto* legacyBox = box.legacyMatrix();
 133     system_.box_          = box;
 134     updateForcerec(&backend_.forcerec_, box.legacyMatrix());
 135     if (TRICLINIC(legacyBox))
 136     {
 137         throw InputException("Only rectangular unit-cells are supported here");
 138     }
 139
 140     const rvec lowerCorner = { 0, 0, 0 };
 141     const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] };
 142
 143     const real particleDensity = static_cast<real>(coordinates.size()) / det(legacyBox);
 144
 145     // If particles are too far outside the box, the grid setup can fail
 146     put_atoms_in_box_omp(PbcType::Xyz, box.legacyMatrix(), coordinates, backend_.numThreads_);
 147
 148     // Put particles on a grid based on bounds specified by the box
 149     nbnxn_put_on_grid(backend_.nbv_.get(),
 150                       legacyBox,
 151                       0,
 152                       lowerCorner,
 153                       upperCorner,
 154                       nullptr,
 155                       { 0, int(coordinates.size()) },
 156                       particleDensity,
 157                       system_.particleInfo_,
 158                       coordinates,
 159                       0,
 160                       nullptr);
 161
 162     backend_.nbv_->constructPairlist(
 163             gmx::InteractionLocality::Local, backend_.exclusions_, 0, &backend_.nrnb_);
 164
 165     // Set Particle Types and Charges and VdW params
 166     backend_.nbv_->setAtomProperties(
 167             system_.particleTypeIdOfAllParticles_, system_.charges_, system_.particleInfo_);
 168     backend_.updatePairlistCalled = true;
 169 }
 170
 171 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 172                                                const Box&                     box,
 173                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 174                                                gmx::ArrayRef<real>            virialOutput,
 175                                                gmx::ArrayRef<real>            energyOutput)
 176 {
 177     if (coordinateInput.size() != forceOutput.size())
 178     {
 179         throw InputException("coordinate array and force buffer size mismatch");
 180     }
 181
 182     if (!backend_.updatePairlistCalled)
 183     {
 184         throw InputException("compute called without updating pairlist at least once");
 185     }
 186
 187     // update the box if changed
 188     if (!(system_.box_ == box))
 189     {
 190         system_.box_ = box;
 191         updateForcerec(&backend_.forcerec_, box.legacyMatrix());
 192     }
 193
 194     bool computeVirial               = !virialOutput.empty();
 195     bool computeEnergies             = !energyOutput.empty();
 196     backend_.stepWork_.computeVirial = computeVirial;
 197     backend_.stepWork_.computeEnergy = computeEnergies;
 198
 199     // update the coordinates in the backend
 200     backend_.nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput);
 201
 202     backend_.nbv_->dispatchNonbondedKernel(
 203             gmx::InteractionLocality::Local,
 204             backend_.interactionConst_,
 205             backend_.stepWork_,
 206             enbvClearFYes,
 207             backend_.forcerec_.shift_vec,
 208             backend_.enerd_.grpp.energyGroupPairTerms[backend_.forcerec_.haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR
 209                                                                                         : NonBondedEnergyTerms::LJSR],
 210             backend_.enerd_.grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
 211             &backend_.nrnb_);
 212
 213     backend_.nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput);
 214
 215     if (computeVirial)
 216     {
 217         // calculate shift forces and turn into an array ref
 218         std::vector<Vec3> shiftForcesVector(gmx::c_numShiftVectors, Vec3(0.0, 0.0, 0.0));
 219         nbnxn_atomdata_add_nbat_fshift_to_fshift(*backend_.nbv_->nbat, shiftForcesVector);
 220         auto shiftForcesRef = constArrayRefFromArray(shiftForcesVector.data(), shiftForcesVector.size());
 221
 222         std::vector<Vec3> shiftVectorsArray(gmx::c_numShiftVectors);
 223
 224         // copy shift vectors from ForceRec
 225         std::copy(backend_.forcerec_.shift_vec.begin(),
 226                   backend_.forcerec_.shift_vec.end(),
 227                   shiftVectorsArray.begin());
 228
 229         computeVirialTensor(
 230                 coordinateInput, forceOutput, shiftVectorsArray, shiftForcesRef, box, virialOutput);
 231     }
 232
 233     // extract term energies (per interaction type)
 234     if (computeEnergies)
 235     {
 236         int nGroupPairs = backend_.enerd_.grpp.nener;
 237         if (int(energyOutput.size()) != int(NonBondedEnergyTerms::Count) * nGroupPairs)
 238         {
 239             throw InputException("Array size for energy output is wrong\n");
 240         }
 241
 242         for (int eg = 0; eg < int(NonBondedEnergyTerms::Count); ++eg)
 243         {
 244             std::copy(begin(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 245                       end(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 246                       energyOutput.begin() + eg * nGroupPairs);
 247         }
 248     }
 249 }
 250
 251 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 252                                                const Box&                     box,
 253                                                gmx::ArrayRef<gmx::RVec>       forceOutput)
 254 {
 255     // compute forces and fill in force buffer
 256     compute(coordinateInput, box, forceOutput, gmx::ArrayRef<real>{}, gmx::ArrayRef<real>{});
 257 }
 258
 259 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 260                                                const Box&                     box,
 261                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 262                                                gmx::ArrayRef<real>            virialOutput)
 263 {
 264     // compute forces and fill in force buffer
 265     compute(coordinateInput, box, forceOutput, virialOutput, gmx::ArrayRef<real>{});
 266 }
 267
 268 GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef<int>  particleTypeIdOfAllParticles,
 269                                                  gmx::ArrayRef<real> nonBondedParams,
 270                                                  gmx::ArrayRef<real> charges,
 271                                                  gmx::ArrayRef<int64_t> particleInteractionFlags,
 272                                                  gmx::ArrayRef<int>     exclusionRanges,
 273                                                  gmx::ArrayRef<int>     exclusionElements,
 274                                                  const NBKernelOptions& options)
 275 {
 276     if (options.useGpu)
 277     {
 278         throw InputException("Use GmxNBForceCalculatorGpu for GPU support");
 279     }
 280
 281     impl_ = std::make_unique<CpuImpl>(particleTypeIdOfAllParticles,
 282                                       nonBondedParams,
 283                                       charges,
 284                                       particleInteractionFlags,
 285                                       exclusionRanges,
 286                                       exclusionElements,
 287                                       options);
 288 }
 289
 290 GmxNBForceCalculatorCpu::~GmxNBForceCalculatorCpu() = default;
 291
 292 //! calculates a new pair list based on new coordinates (for every NS step)
 293 void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef<gmx::RVec> coordinates, const Box& box)
 294 {
 295     impl_->updatePairlist(coordinates, box);
 296 }
 297
 298 //! Compute forces and return
 299 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 300                                       const Box&                     box,
 301                                       gmx::ArrayRef<gmx::RVec>       forceOutput)
 302 {
 303     impl_->compute(coordinateInput, box, forceOutput);
 304 }
 305
 306 //! Compute forces and virial tensor
 307 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 308                                       const Box&                     box,
 309                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 310                                       gmx::ArrayRef<real>            virialOutput)
 311 {
 312     impl_->compute(coordinateInput, box, forceOutput, virialOutput);
 313 }
 314
 315 //! Compute forces, virial tensor and potential energies
 316 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 317                                       const Box&                     box,
 318                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 319                                       gmx::ArrayRef<real>            virialOutput,
 320                                       gmx::ArrayRef<real>            energyOutput)
 321 {
 322     impl_->compute(coordinateInput, box, forceOutput, virialOutput, energyOutput);
 323 }
 324
 325 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
 326                                                                     const NBKernelOptions& options)
 327 {
 328     std::vector<real> nonBondedParameters = createNonBondedParameters(
 329             topology.getParticleTypes(), topology.getNonBondedInteractionMap());
 330
 331     std::vector<int64_t> particleInteractionFlags = createParticleInfoAllVdw(topology.numParticles());
 332
 333     return std::make_unique<GmxNBForceCalculatorCpu>(topology.getParticleTypeIdOfAllParticles(),
 334                                                      nonBondedParameters,
 335                                                      topology.getCharges(),
 336                                                      particleInteractionFlags,
 337                                                      topology.exclusionLists().ListRanges,
 338                                                      topology.exclusionLists().ListElements,
 339                                                      options);
 340 }
 341
 342 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(TprReader& tprReader,
 343                                                                     const NBKernelOptions& options)
 344 {
 345     return std::make_unique<GmxNBForceCalculatorCpu>(tprReader.particleTypeIdOfAllParticles_,
 346                                                      tprReader.nonbondedParameters_,
 347                                                      tprReader.charges_,
 348                                                      tprReader.particleInteractionFlags_,
 349                                                      tprReader.exclusionListRanges_,
 350                                                      tprReader.exclusionListElements_,
 351                                                      options);
 352 }
 353
 354 } // namespace nblib