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/lincs_cuda.cuh"
  68 #include "gromacs/mdlib/update_constrain_cuda.h"
  69 #include "gromacs/pbcutil/pbc.h"
  70 #include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
  71
  72 #include "leapfrog_cuda_impl.h"
  73 #include "settle_cuda_impl.h"
  74
  75 namespace gmx
  76 {
  77
  78 void UpdateConstrainCuda::Impl::integrate(const real  dt,
  79                                           const bool  updateVelocities,
  80                                           const bool  computeVirial,
  81                                           tensor      virial)
  82 {
  83     // Clearing virial matrix
  84     // TODO There is no point in having saparate virial matrix for constraints
  85     clear_mat(virial);
  86
  87     integrator_->integrate(d_x_, d_xp_, d_v_, d_f_, dt);
  88     lincsCuda_->apply(d_x_, d_xp_,
  89                       updateVelocities, d_v_, 1.0/dt,
  90                       computeVirial, virial);
  91     settleCuda_->apply(d_x_, d_xp_,
  92                        updateVelocities, d_v_, 1.0/dt,
  93                        computeVirial, virial);
  94
  95     // scaledVirial -> virial (methods above returns scaled values)
  96     float scaleFactor = 0.5f/(dt*dt);
  97     for (int i = 0; i < DIM; i++)
  98     {
  99         for (int j = 0; j < DIM; j++)
 100         {
 101             virial[i][j] = scaleFactor*virial[i][j];
 102         }
 103     }
 104
 105     return;
 106 }
 107
 108 UpdateConstrainCuda::Impl::Impl(int                numAtoms,
 109                                 const t_inputrec  &ir,
 110                                 const gmx_mtop_t  &mtop)
 111     : numAtoms_(numAtoms)
 112 {
 113     allocateDeviceBuffer(&d_x_,              numAtoms, nullptr);
 114     allocateDeviceBuffer(&d_xp_,             numAtoms, nullptr);
 115     allocateDeviceBuffer(&d_v_,              numAtoms, nullptr);
 116     allocateDeviceBuffer(&d_f_,              numAtoms, nullptr);
 117     allocateDeviceBuffer(&d_inverseMasses_,  numAtoms, nullptr);
 118
 119     // TODO When the code will be integrated into the schedule, it will be assigned non-default stream.
 120     stream_ = nullptr;
 121
 122     GMX_RELEASE_ASSERT(numAtoms == mtop.natoms, "State and topology number of atoms should be the same.");
 123     integrator_ = std::make_unique<LeapFrogCuda::Impl>();
 124     lincsCuda_  = std::make_unique<LincsCuda>(ir.nLincsIter, ir.nProjOrder);
 125     settleCuda_ = std::make_unique<SettleCuda::Impl>(mtop);
 126
 127 }
 128
 129 UpdateConstrainCuda::Impl::~Impl()
 130 {
 131 }
 132
 133 void UpdateConstrainCuda::Impl::set(const t_idef      &idef,
 134                                     const t_mdatoms   &md)
 135 {
 136     // Integrator should also update something, but it does not even have a method yet
 137     integrator_->set(md);
 138     lincsCuda_->set(idef, md);
 139     settleCuda_->set(idef, md);
 140 }
 141
 142 void UpdateConstrainCuda::Impl::setPbc(const t_pbc *pbc)
 143 {
 144     setPbcAiuc(pbc->ndim_ePBC, pbc->box, &pbcAiuc_);
 145     integrator_->setPbc(pbc);
 146     lincsCuda_->setPbc(pbc);
 147     settleCuda_->setPbc(pbc);
 148 }
 149
 150 void UpdateConstrainCuda::Impl::copyCoordinatesToGpu(const rvec *h_x)
 151 {
 152     copyToDeviceBuffer(&d_x_, (float3*)h_x, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 153 }
 154
 155 void UpdateConstrainCuda::Impl::copyVelocitiesToGpu(const rvec *h_v)
 156 {
 157     copyToDeviceBuffer(&d_v_, (float3*)h_v, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 158 }
 159
 160 void UpdateConstrainCuda::Impl::copyForcesToGpu(const rvec *h_f)
 161 {
 162     copyToDeviceBuffer(&d_f_, (float3*)h_f, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 163 }
 164
 165 void UpdateConstrainCuda::Impl::copyCoordinatesFromGpu(rvec *h_xp)
 166 {
 167     copyFromDeviceBuffer((float3*)h_xp, &d_xp_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 168 }
 169
 170 void UpdateConstrainCuda::Impl::copyVelocitiesFromGpu(rvec *h_v)
 171 {
 172     copyFromDeviceBuffer((float3*)h_v, &d_v_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 173 }
 174
 175 void UpdateConstrainCuda::Impl::copyForcesFromGpu(rvec *h_f)
 176 {
 177     copyFromDeviceBuffer((float3*)h_f, &d_f_, 0, numAtoms_, stream_, GpuApiCallBehavior::Sync, nullptr);
 178 }
 179
 180 void UpdateConstrainCuda::Impl::setXVFPointers(rvec *d_x, rvec *d_xp, rvec *d_v, rvec *d_f)
 181 {
 182     d_x_  = (float3*)d_x;
 183     d_xp_ = (float3*)d_xp;
 184     d_v_  = (float3*)d_v;
 185     d_f_  = (float3*)d_f;
 186 }
 187
 188
 189 UpdateConstrainCuda::UpdateConstrainCuda(int                numAtoms,
 190                                          const t_inputrec  &ir,
 191                                          const gmx_mtop_t  &mtop)
 192     : impl_(new Impl(numAtoms, ir, mtop))
 193 {
 194 }
 195
 196 UpdateConstrainCuda::~UpdateConstrainCuda() = default;
 197
 198 void UpdateConstrainCuda::integrate(const real  dt,
 199                                     const bool  updateVelocities,
 200                                     const bool  computeVirial,
 201                                     tensor      virialScaled)
 202 {
 203     impl_->integrate(dt, updateVelocities, computeVirial, virialScaled);
 204 }
 205
 206 void UpdateConstrainCuda::set(const t_idef               &idef,
 207                               const t_mdatoms gmx_unused &md)
 208 {
 209     impl_->set(idef, md);
 210 }
 211
 212 void UpdateConstrainCuda::setPbc(const t_pbc *pbc)
 213 {
 214     impl_->setPbc(pbc);
 215 }
 216
 217 void UpdateConstrainCuda::copyCoordinatesToGpu(const rvec *h_x)
 218 {
 219     impl_->copyCoordinatesToGpu(h_x);
 220 }
 221
 222 void UpdateConstrainCuda::copyVelocitiesToGpu(const rvec *h_v)
 223 {
 224     impl_->copyVelocitiesToGpu(h_v);
 225 }
 226
 227 void UpdateConstrainCuda::copyForcesToGpu(const rvec *h_f)
 228 {
 229     impl_->copyForcesToGpu(h_f);
 230 }
 231
 232 void UpdateConstrainCuda::copyCoordinatesFromGpu(rvec *h_xp)
 233 {
 234     impl_->copyCoordinatesFromGpu(h_xp);
 235 }
 236
 237 void UpdateConstrainCuda::copyVelocitiesFromGpu(rvec *h_v)
 238 {
 239     impl_->copyVelocitiesFromGpu(h_v);
 240 }
 241
 242 void UpdateConstrainCuda::copyForcesFromGpu(rvec *h_f)
 243 {
 244     impl_->copyForcesFromGpu(h_f);
 245 }
 246
 247 void UpdateConstrainCuda::setXVFPointers(rvec *d_x, rvec *d_xp, rvec *d_v, rvec *d_f)
 248 {
 249     impl_->setXVFPointers(d_x, d_xp, d_v, d_f);
 250 }
 251
 252 } //namespace gmx