src/gromacs/mdlib/update_constrain_cuda_impl.cu

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, 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  *
  37  * \brief Implements update and constraints class using CUDA.
  38  *
  39  * The class combines Leap-Frog integrator with LINCS and SETTLE constraints.
  40  *
  41  * \todo This class should take over the management of coordinates, velocities
  42  *       forces, virial, and PBC from its members (i.e. from Leap-Frog, LINCS
  43  *       and SETTLE).
  44  * \todo The computational procedures in members should be integrated to improve
  45  *       computational performance.
  46  *
  47  * \author Artem Zhmurov <zhmurov@gmail.com>
  48  *
  49  * \ingroup module_mdlib
  50  */
  51 #include "gmxpre.h"
  52
  53 #include "update_constrain_cuda_impl.h"
  54
  55 #include <assert.h>
  56 #include <stdio.h>
  57
  58 #include <cmath>
  59
  60 #include <algorithm>
  61
  62 #include "gromacs/gpu_utils/cudautils.cuh"
  63 #include "gromacs/gpu_utils/devicebuffer.cuh"
  64 #include "gromacs/gpu_utils/gputraits.cuh"
  65 #include "gromacs/gpu_utils/vectype_ops.cuh"
  66 #include "gromacs/math/vec.h"
  67 #include "gromacs/mdlib/leapfrog_cuda.h"
  68 #include "gromacs/mdlib/lincs_cuda.h"
  69 #include "gromacs/mdlib/settle_cuda.h"
  70 #include "gromacs/mdlib/update_constrain_cuda.h"
  71 #include "gromacs/pbcutil/pbc.h"
  72 #include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
  73
  74 namespace gmx
  75 {
  76
  77 /*! \brief Integrate
  78  *
  79  * Integrates the equation of motion using Leap-Frog algorithm and applies
  80  * LINCS and SETTLE constraints.
  81  * Updates d_xp_ and d_v_ fields of this object.
  82  *
  83  * \param[in] dt                Timestep
  84  * \param[in] updateVelocities  If the velocities should be constrained.
  85  * \param[in] computeVirial     If virial should be updated.
  86  * \param[out] virial           Place to save virial tensor.
  87  */
  88 void UpdateConstrainCuda::Impl::integrate(const real  dt,
  89                                           const bool  updateVelocities,
  90                                           const bool  computeVirial,
  91                                           tensor      virial)
  92 {
  93     // Clearing virial matrix
  94     // TODO There is no point in having saparate virial matrix for constraints
  95     clear_mat(virial);
  96
  97     integrator_->integrate(dt);
  98     lincsCuda_->apply(updateVelocities, 1.0/dt, computeVirial, virial);
  99     settleCuda_->apply(updateVelocities, 1.0/dt, computeVirial, virial);
 100
 101     // scaledVirial -> virial (methods above returns scaled values)
 102     float scaleFactor = 0.5f/(dt*dt);
 103     for (int i = 0; i < DIM; i++)
 104     {
 105         for (int j = 0; j < DIM; j++)
 106         {
 107             virial[i][j] = scaleFactor*virial[i][j];
 108         }
 109     }
 110
 111     return;
 112 }
 113
 114 /*! \brief Create Update-Constrain object
 115  *
 116  * \param[in] numAtoms  Number of atoms.
 117  * \param[in] ir        Input record data: LINCS takes number of iterations and order of
 118  *                      projection from it.
 119  * \param[in] mtop      Topology of the system: SETTLE gets the masses for O and H atoms
 120  *                      and target O-H and H-H distances from this object.
 121  */
 122 UpdateConstrainCuda::Impl::Impl(int                numAtoms,
 123                                 const t_inputrec  &ir,
 124                                 const gmx_mtop_t  &mtop)
 125     : numAtoms_(numAtoms)
 126 {
 127     allocateDeviceBuffer(&d_x_,              numAtoms, nullptr);
 128     allocateDeviceBuffer(&d_xp_,             numAtoms, nullptr);
 129     allocateDeviceBuffer(&d_v_,              numAtoms, nullptr);
 130     allocateDeviceBuffer(&d_f_,              numAtoms, nullptr);
 131     allocateDeviceBuffer(&d_inverseMasses_,  numAtoms, nullptr);
 132
 133     // TODO When the code will be integrated into the schedule, it will be assigned non-default stream.
 134     stream_ = nullptr;
 135
 136     GMX_RELEASE_ASSERT(numAtoms == mtop.natoms, "State and topology number of atoms should be the same.");
 137     integrator_ = std::make_unique<LeapFrogCuda>(numAtoms);
 138     lincsCuda_  = std::make_unique<LincsCuda>(mtop.natoms, ir.nLincsIter, ir.nProjOrder);
 139     settleCuda_ = std::make_unique<SettleCuda>(mtop.natoms, mtop);
 140
 141     integrator_->setXVFPointers((rvec*)d_x_, (rvec*)d_xp_, (rvec*)d_v_, (rvec*)d_f_);
 142     lincsCuda_->setXVPointers((rvec*)d_x_, (rvec*)d_xp_, (rvec*)d_v_);
 143     settleCuda_->setXVPointers((rvec*)d_x_, (rvec*)d_xp_,  (rvec*)d_v_);
 144 }
 145
 146 UpdateConstrainCuda::Impl::~Impl()
 147 {
 148 }
 149
 150 /*! \brief
 151  * Update data-structures (e.g. after NB search step).
 152  *
 153  * \param[in] idef    System topology
 154  * \param[in] md      Atoms data.
 155  */
 156 void UpdateConstrainCuda::Impl::set(const t_idef      &idef,
 157                                     const t_mdatoms   &md)
 158 {
 159     // Integrator should also update something, but it does not even have a method yet
 160     integrator_->set(md);
 161     lincsCuda_->set(idef, md);
 162     settleCuda_->set(idef, md);
 163 }
 164
 165 /*! \brief
 166  * Update PBC data.
 167  *
 168  * Converts pbc data from t_pbc into the PbcAiuc format and stores the latter.
 169  *
 170  * \param[in] pbc The PBC data in t_pbc format.
 171  */
 172 void UpdateConstrainCuda::Impl::setPbc(const t_pbc *pbc)
 173 {
 174     setPbcAiuc(pbc->ndim_ePBC, pbc->box, &pbcAiuc_);
 175     integrator_->setPbc(pbc);
 176     lincsCuda_->setPbc(pbc);
 177     settleCuda_->setPbc(pbc);
 178 }
 179
 180 /*! \brief
 181  * Copy coordinates from CPU to GPU.
 182  *
 183  * The data are assumed to be in float3/fvec format (single precision).
 184  *
 185  * \param[in] h_x  CPU pointer where coordinates should be copied from.
 186  */
 187 void UpdateConstrainCuda::Impl::copyCoordinatesToGpu(const rvec *h_x)
 188 {
 189     copyToDeviceBuffer(&d_x_, (float3*)h_x, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 190 }
 191
 192 /*! \brief
 193  * Copy velocities from CPU to GPU.
 194  *
 195  * The data are assumed to be in float3/fvec format (single precision).
 196  *
 197  * \param[in] h_v  CPU pointer where velocities should be copied from.
 198  */
 199 void UpdateConstrainCuda::Impl::copyVelocitiesToGpu(const rvec *h_v)
 200 {
 201     copyToDeviceBuffer(&d_v_, (float3*)h_v, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 202 }
 203
 204 /*! \brief
 205  * Copy forces from CPU to GPU.
 206  *
 207  * The data are assumed to be in float3/fvec format (single precision).
 208  *
 209  * \param[in] h_f  CPU pointer where forces should be copied from.
 210  */
 211 void UpdateConstrainCuda::Impl::copyForcesToGpu(const rvec *h_f)
 212 {
 213     copyToDeviceBuffer(&d_f_, (float3*)h_f, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 214 }
 215
 216 /*! \brief
 217  * Copy coordinates from GPU to CPU.
 218  *
 219  * The data are assumed to be in float3/fvec format (single precision).
 220  *
 221  * \param[out] h_xp CPU pointer where coordinates should be copied to.
 222  */
 223 void UpdateConstrainCuda::Impl::copyCoordinatesFromGpu(rvec *h_xp)
 224 {
 225     copyFromDeviceBuffer((float3*)h_xp, &d_xp_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 226 }
 227
 228 /*! \brief
 229  * Copy velocities from GPU to CPU.
 230  *
 231  * The velocities are assumed to be in float3/fvec format (single precision).
 232  *
 233  * \param[in] h_v  Pointer to velocities data.
 234  */
 235 void UpdateConstrainCuda::Impl::copyVelocitiesFromGpu(rvec *h_v)
 236 {
 237     copyFromDeviceBuffer((float3*)h_v, &d_v_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 238 }
 239
 240 /*! \brief
 241  * Copy forces from GPU to CPU.
 242  *
 243  * The forces are assumed to be in float3/fvec format (single precision).
 244  *
 245  * \param[in] h_f  Pointer to forces data.
 246  */
 247 void UpdateConstrainCuda::Impl::copyForcesFromGpu(rvec *h_f)
 248 {
 249     copyFromDeviceBuffer((float3*)h_f, &d_f_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 250 }
 251
 252 /*! \brief
 253  * Set the internal GPU-memory x, xprime and v pointers.
 254  *
 255  * Data is not copied. The data are assumed to be in float3/fvec format
 256  * (float3 is used internally, but the data layout should be identical).
 257  *
 258  * \param[in] d_x  Pointer to the coordinates for the input (on GPU)
 259  * \param[in] d_xp Pointer to the coordinates for the output (on GPU)
 260  * \param[in] d_v  Pointer to the velocities (on GPU)
 261  * \param[in] d_f  Pointer to the forces (on GPU)
 262  */
 263 void UpdateConstrainCuda::Impl::setXVFPointers(rvec *d_x, rvec *d_xp, rvec *d_v, rvec *d_f)
 264 {
 265     d_x_  = (float3*)d_x;
 266     d_xp_ = (float3*)d_xp;
 267     d_v_  = (float3*)d_v;
 268     d_f_  = (float3*)d_f;
 269 }
 270
 271
 272 UpdateConstrainCuda::UpdateConstrainCuda(int                numAtoms,
 273                                          const t_inputrec  &ir,
 274                                          const gmx_mtop_t  &mtop)
 275     : impl_(new Impl(numAtoms, ir, mtop))
 276 {
 277 }
 278
 279 UpdateConstrainCuda::~UpdateConstrainCuda() = default;
 280
 281 void UpdateConstrainCuda::integrate(const real  dt,
 282                                     const bool  updateVelocities,
 283                                     const bool  computeVirial,
 284                                     tensor      virialScaled)
 285 {
 286     impl_->integrate(dt, updateVelocities, computeVirial, virialScaled);
 287 }
 288
 289 void UpdateConstrainCuda::set(const t_idef               &idef,
 290                               const t_mdatoms gmx_unused &md)
 291 {
 292     impl_->set(idef, md);
 293 }
 294
 295 void UpdateConstrainCuda::setPbc(const t_pbc *pbc)
 296 {
 297     impl_->setPbc(pbc);
 298 }
 299
 300 void UpdateConstrainCuda::copyCoordinatesToGpu(const rvec *h_x)
 301 {
 302     impl_->copyCoordinatesToGpu(h_x);
 303 }
 304
 305 void UpdateConstrainCuda::copyVelocitiesToGpu(const rvec *h_v)
 306 {
 307     impl_->copyVelocitiesToGpu(h_v);
 308 }
 309
 310 void UpdateConstrainCuda::copyForcesToGpu(const rvec *h_f)
 311 {
 312     impl_->copyForcesToGpu(h_f);
 313 }
 314
 315 void UpdateConstrainCuda::copyCoordinatesFromGpu(rvec *h_xp)
 316 {
 317     impl_->copyCoordinatesFromGpu(h_xp);
 318 }
 319
 320 void UpdateConstrainCuda::copyVelocitiesFromGpu(rvec *h_v)
 321 {
 322     impl_->copyVelocitiesFromGpu(h_v);
 323 }
 324
 325 void UpdateConstrainCuda::copyForcesFromGpu(rvec *h_f)
 326 {
 327     impl_->copyForcesFromGpu(h_f);
 328 }
 329
 330 void UpdateConstrainCuda::setXVFPointers(rvec *d_x, rvec *d_xp, rvec *d_v, rvec *d_f)
 331 {
 332     impl_->setXVFPointers(d_x, d_xp, d_v, d_f);
 333 }
 334
 335 } //namespace gmx