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/virials.h"
  61
  62 namespace nblib
  63 {
  64
  65 class GmxNBForceCalculatorCpu::CpuImpl final
  66 {
  67 public:
  68     CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
  69             gmx::ArrayRef<real>    nonBondedParams,
  70             gmx::ArrayRef<real>    charges,
  71             gmx::ArrayRef<int64_t> particleInteractionFlags,
  72             gmx::ArrayRef<int>     exclusionRanges,
  73             gmx::ArrayRef<int>     exclusionElements,
  74             const NBKernelOptions& options);
  75
  76     //! calculates a new pair list based on new coordinates (for every NS step)
  77     void updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box);
  78
  79     //! Compute forces and return
  80     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  81                  const Box&                     box,
  82                  gmx::ArrayRef<gmx::RVec>       forceOutput);
  83
  84     //! Compute forces and virial tensor
  85     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  86                  const Box&                     box,
  87                  gmx::ArrayRef<gmx::RVec>       forceOutput,
  88                  gmx::ArrayRef<real>            virialOutput);
  89
  90     //! Compute forces, virial tensor and potential energies
  91     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
  92                  const Box&                     box,
  93                  gmx::ArrayRef<gmx::RVec>       forceOutput,
  94                  gmx::ArrayRef<real>            virialOutput,
  95                  gmx::ArrayRef<real>            energyOutput);
  96
  97 private:
  98     //! \brief client-side provided system description data
  99     SystemDescription system_;
 100
 101     //! \brief Gmx backend objects, employed for calculating the forces
 102     GmxBackendData backend_;
 103 };
 104
 105 GmxNBForceCalculatorCpu::CpuImpl::CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
 106                                           gmx::ArrayRef<real>    nonBondedParams,
 107                                           gmx::ArrayRef<real>    charges,
 108                                           gmx::ArrayRef<int64_t> particleInteractionFlags,
 109                                           gmx::ArrayRef<int>     exclusionRanges,
 110                                           gmx::ArrayRef<int>     exclusionElements,
 111                                           const NBKernelOptions& options) :
 112     system_(SystemDescription(particleTypeIdOfAllParticles, nonBondedParams, charges, particleInteractionFlags)),
 113     backend_(GmxBackendData(options, findNumEnergyGroups(particleInteractionFlags), exclusionRanges, exclusionElements))
 114 {
 115     // Set up non-bonded verlet in the backend
 116     backend_.nbv_ = createNbnxmCPU(system_.numParticleTypes_,
 117                                    options,
 118                                    findNumEnergyGroups(particleInteractionFlags),
 119                                    system_.nonBondedParams_);
 120 }
 121
 122 void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates,
 123                                                       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     // Put particles on a grid based on bounds specified by the box
 146     nbnxn_put_on_grid(backend_.nbv_.get(),
 147                       legacyBox,
 148                       0,
 149                       lowerCorner,
 150                       upperCorner,
 151                       nullptr,
 152                       { 0, int(coordinates.size()) },
 153                       particleDensity,
 154                       system_.particleInfo_,
 155                       coordinates,
 156                       0,
 157                       nullptr);
 158
 159     backend_.nbv_->constructPairlist(
 160             gmx::InteractionLocality::Local, backend_.exclusions_, 0, &backend_.nrnb_);
 161
 162     // Set Particle Types and Charges and VdW params
 163     backend_.nbv_->setAtomProperties(
 164             system_.particleTypeIdOfAllParticles_, system_.charges_, system_.particleInfo_);
 165     backend_.updatePairlistCalled = true;
 166 }
 167
 168 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 169                                                const Box&                     box,
 170                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 171                                                gmx::ArrayRef<real>            virialOutput,
 172                                                gmx::ArrayRef<real>            energyOutput)
 173 {
 174     if (coordinateInput.size() != forceOutput.size())
 175     {
 176         throw InputException("coordinate array and force buffer size mismatch");
 177     }
 178
 179     if (!backend_.updatePairlistCalled)
 180     {
 181         throw InputException("compute called without updating pairlist at least once");
 182     }
 183
 184     // update the box if changed
 185     if (!(system_.box_ == box))
 186     {
 187         system_.box_ = box;
 188         updateForcerec(&backend_.forcerec_, box.legacyMatrix());
 189     }
 190
 191     bool computeVirial               = !virialOutput.empty();
 192     bool computeEnergies             = !energyOutput.empty();
 193     backend_.stepWork_.computeVirial = computeVirial;
 194     backend_.stepWork_.computeEnergy = computeEnergies;
 195
 196     // update the coordinates in the backend
 197     backend_.nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput);
 198
 199     backend_.nbv_->dispatchNonbondedKernel(
 200             gmx::InteractionLocality::Local,
 201             backend_.interactionConst_,
 202             backend_.stepWork_,
 203             enbvClearFYes,
 204             backend_.forcerec_.shift_vec,
 205             backend_.enerd_.grpp.energyGroupPairTerms[backend_.forcerec_.haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR
 206                                                                                         : NonBondedEnergyTerms::LJSR],
 207             backend_.enerd_.grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
 208             &backend_.nrnb_);
 209
 210     backend_.nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput);
 211
 212     if (computeVirial)
 213     {
 214         // calculate shift forces and turn into an array ref
 215         std::vector<Vec3> shiftForcesVector(gmx::c_numShiftVectors, Vec3(0.0, 0.0, 0.0));
 216         nbnxn_atomdata_add_nbat_fshift_to_fshift(*backend_.nbv_->nbat, shiftForcesVector);
 217         auto shiftForcesRef = constArrayRefFromArray(shiftForcesVector.data(), shiftForcesVector.size());
 218
 219         std::vector<Vec3> shiftVectorsArray(gmx::c_numShiftVectors);
 220
 221         // copy shift vectors from ForceRec
 222         std::copy(backend_.forcerec_.shift_vec.begin(),
 223                   backend_.forcerec_.shift_vec.end(),
 224                   shiftVectorsArray.begin());
 225
 226         computeVirialTensor(
 227                 coordinateInput, forceOutput, shiftVectorsArray, shiftForcesRef, box, virialOutput);
 228     }
 229
 230     // extract term energies (per interaction type)
 231     if (computeEnergies)
 232     {
 233         int nGroupPairs = backend_.enerd_.grpp.nener;
 234         if (int(energyOutput.size()) != int(NonBondedEnergyTerms::Count) * nGroupPairs)
 235         {
 236             throw InputException("Array size for energy output is wrong\n");
 237         }
 238
 239         for (int eg = 0; eg < int(NonBondedEnergyTerms::Count); ++eg)
 240         {
 241             std::copy(begin(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 242                       end(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
 243                       energyOutput.begin() + eg * nGroupPairs);
 244         }
 245     }
 246 }
 247
 248 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 249                                                const Box&                     box,
 250                                                gmx::ArrayRef<gmx::RVec>       forceOutput)
 251 {
 252     // compute forces and fill in force buffer
 253     compute(coordinateInput, box, forceOutput, gmx::ArrayRef<real>{}, gmx::ArrayRef<real>{});
 254 }
 255
 256 void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 257                                                const Box&                     box,
 258                                                gmx::ArrayRef<gmx::RVec>       forceOutput,
 259                                                gmx::ArrayRef<real>            virialOutput)
 260 {
 261     // compute forces and fill in force buffer
 262     compute(coordinateInput, box, forceOutput, virialOutput, gmx::ArrayRef<real>{});
 263 }
 264
 265 GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef<int>  particleTypeIdOfAllParticles,
 266                                                  gmx::ArrayRef<real> nonBondedParams,
 267                                                  gmx::ArrayRef<real> charges,
 268                                                  gmx::ArrayRef<int64_t> particleInteractionFlags,
 269                                                  gmx::ArrayRef<int>     exclusionRanges,
 270                                                  gmx::ArrayRef<int>     exclusionElements,
 271                                                  const NBKernelOptions& options)
 272 {
 273     if (options.useGpu)
 274     {
 275         throw InputException("Use GmxNBForceCalculatorGpu for GPU support");
 276     }
 277
 278     impl_ = std::make_unique<CpuImpl>(particleTypeIdOfAllParticles,
 279                                       nonBondedParams,
 280                                       charges,
 281                                       particleInteractionFlags,
 282                                       exclusionRanges,
 283                                       exclusionElements,
 284                                       options);
 285 }
 286
 287 GmxNBForceCalculatorCpu::~GmxNBForceCalculatorCpu() = default;
 288
 289 //! calculates a new pair list based on new coordinates (for every NS step)
 290 void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box)
 291 {
 292     impl_->updatePairlist(coordinates, box);
 293 }
 294
 295 //! Compute forces and return
 296 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 297                                       const Box&                     box,
 298                                       gmx::ArrayRef<gmx::RVec>       forceOutput)
 299 {
 300     impl_->compute(coordinateInput, box, forceOutput);
 301 }
 302
 303 //! Compute forces and virial tensor
 304 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 305                                       const Box&                     box,
 306                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 307                                       gmx::ArrayRef<real>            virialOutput)
 308 {
 309     impl_->compute(coordinateInput, box, forceOutput, virialOutput);
 310 }
 311
 312 //! Compute forces, virial tensor and potential energies
 313 void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
 314                                       const Box&                     box,
 315                                       gmx::ArrayRef<gmx::RVec>       forceOutput,
 316                                       gmx::ArrayRef<real>            virialOutput,
 317                                       gmx::ArrayRef<real>            energyOutput)
 318 {
 319     impl_->compute(coordinateInput, box, forceOutput, virialOutput, energyOutput);
 320 }
 321
 322 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
 323                                                                     const NBKernelOptions& options)
 324 {
 325     std::vector<real> nonBondedParameters = createNonBondedParameters(
 326             topology.getParticleTypes(), topology.getNonBondedInteractionMap());
 327
 328     std::vector<int64_t> particleInteractionFlags = createParticleInfoAllVdw(topology.numParticles());
 329
 330     return std::make_unique<GmxNBForceCalculatorCpu>(topology.getParticleTypeIdOfAllParticles(),
 331                                                      nonBondedParameters,
 332                                                      topology.getCharges(),
 333                                                      particleInteractionFlags,
 334                                                      topology.exclusionLists().ListRanges,
 335                                                      topology.exclusionLists().ListElements,
 336                                                      options);
 337 }
 338
 339 } // namespace nblib