libhoomd/cuda/TwoStepNVTMTKGPU.cu

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 // Maintainer: jglaser
  51
  52 #include "TwoStepNVTMTKGPU.cuh"
  53
  54 #ifdef WIN32
  55 #include <cassert>
  56 #else
  57 #include <assert.h>
  58 #endif
  59
  60 /*! \file TwoStepNVTGPU.cu
  61     \brief Defines GPU kernel code for NVT integration on the GPU. Used by TwoStepNVTGPU.
  62 */
  63
  64 //! Takes the first 1/2 step forward in the NVT integration step
  65 /*! \param d_pos array of particle positions
  66     \param d_vel array of particle velocities
  67     \param d_accel array of particle accelerations
  68     \param d_image array of particle images
  69     \param d_group_members Device array listing the indicies of the mebers of the group to integrate
  70     \param group_size Number of members in the group
  71     \param box Box dimensions for periodic boundary condition handling
  72     \param exp_fac Velocity rescaling factor from thermostat
  73     \param deltaT Amount of real time to step forward in one time step
  74
  75     Take the first half step forward in the NVT integration.
  76
  77     See gpu_nve_step_one_kernel() for some performance notes on how to handle the group data reads efficiently.
  78 */
  79 extern "C" __global__
  80 void gpu_nvt_mtk_step_one_kernel(Scalar4 *d_pos,
  81                              Scalar4 *d_vel,
  82                              const Scalar3 *d_accel,
  83                              int3 *d_image,
  84                              unsigned int *d_group_members,
  85                              unsigned int group_size,
  86                              BoxDim box,
  87                              Scalar exp_fac,
  88                              Scalar deltaT)
  89     {
  90     // determine which particle this thread works on
  91     int group_idx = blockIdx.x * blockDim.x + threadIdx.x;
  92
  93     if (group_idx < group_size)
  94         {
  95         unsigned int idx = d_group_members[group_idx];
  96
  97         // update positions to the next timestep and update velocities to the next half step
  98         Scalar4 postype = d_pos[idx];
  99         Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z);
 100
 101         Scalar4 velmass = d_vel[idx];
 102         Scalar3 vel = make_scalar3(velmass.x, velmass.y, velmass.z);
 103         Scalar3 accel = d_accel[idx];
 104
 105         // perform update computation
 106         vel = vel*exp_fac + Scalar(1.0/2.0) * accel * deltaT;
 107         pos += vel * deltaT;
 108
 109         // read in the image flags
 110         int3 image = d_image[idx];
 111
 112         // time to fix the periodic boundary conditions
 113         box.wrap(pos, image);
 114
 115         // write out the results
 116         d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, postype.w);
 117         d_vel[idx] = make_scalar4(vel.x, vel.y, vel.z, velmass.w);
 118         d_image[idx] = image;
 119         }
 120     }
 121
 122 /*! \param d_pos array of particle positions
 123     \param d_vel array of particle velocities
 124     \param d_accel array of particle accelerations
 125     \param d_image array of particle images
 126     \param d_group_members Device array listing the indicies of the mebers of the group to integrate
 127     \param group_size Number of members in the group
 128     \param box Box dimensions for periodic boundary condition handling
 129     \param block_size Size of the block to run
 130     \param exp_fac Thermostat rescaling factor
 131     \param deltaT Amount of real time to step forward in one time step
 132 */
 133 cudaError_t gpu_nvt_mtk_step_one(Scalar4 *d_pos,
 134                              Scalar4 *d_vel,
 135                              const Scalar3 *d_accel,
 136                              int3 *d_image,
 137                              unsigned int *d_group_members,
 138                              unsigned int group_size,
 139                              const BoxDim& box,
 140                              unsigned int block_size,
 141                              Scalar exp_fac,
 142                              Scalar deltaT)
 143     {
 144     static unsigned int max_block_size = UINT_MAX;
 145     if (max_block_size == UINT_MAX)
 146         {
 147         cudaFuncAttributes attr;
 148         cudaFuncGetAttributes(&attr, (const void *)gpu_nvt_mtk_step_one_kernel);
 149         max_block_size = attr.maxThreadsPerBlock;
 150         }
 151
 152     unsigned int run_block_size = min(block_size, max_block_size);
 153
 154     // setup the grid to run the kernel
 155     dim3 grid( (group_size/run_block_size) + 1, 1, 1);
 156     dim3 threads(run_block_size, 1, 1);
 157
 158     // run the kernel
 159     gpu_nvt_mtk_step_one_kernel<<< grid, threads >>>(d_pos,
 160                          d_vel,
 161                          d_accel,
 162                          d_image,
 163                          d_group_members,
 164                          group_size,
 165                          box,
 166                          exp_fac,
 167                          deltaT);
 168     return cudaSuccess;
 169     }
 170
 171 //! Takes the second 1/2 step forward in the NVT integration step
 172 /*! \param d_vel array of particle velocities
 173     \param d_accel array of particle accelerations
 174     \param d_group_members Device array listing the indicies of the mebers of the group to integrate
 175     \param group_size Number of members in the group
 176     \param d_net_force Net force on each particle
 177     \param deltaT Amount of real time to step forward in one time step
 178 */
 179 extern "C" __global__
 180 void gpu_nvt_mtk_step_two_kernel(Scalar4 *d_vel,
 181                              Scalar3 *d_accel,
 182                              unsigned int *d_group_members,
 183                              unsigned int group_size,
 184                              Scalar4 *d_net_force,
 185                              Scalar deltaT)
 186     {
 187     // determine which particle this thread works on
 188     int group_idx = blockIdx.x * blockDim.x + threadIdx.x;
 189
 190     if (group_idx < group_size)
 191         {
 192         unsigned int idx = d_group_members[group_idx];
 193
 194         // read in the net force and calculate the acceleration
 195         Scalar4 net_force = d_net_force[idx];
 196         Scalar3 accel = make_scalar3(net_force.x,net_force.y,net_force.z);
 197
 198         Scalar4 vel = d_vel[idx];
 199         Scalar3 v = make_scalar3(vel.x,vel.y,vel.z);
 200
 201         Scalar mass = vel.w;
 202         accel = accel/mass;
 203         v += Scalar(1.0/2.0) * deltaT * accel;
 204
 205         // write out data
 206         d_vel[idx] = make_scalar4(v.x,v.y,v.z,vel.w);
 207
 208         // since we calculate the acceleration, we need to write it for the next step
 209         d_accel[idx] = accel;
 210         }
 211     }
 212
 213 /*! \param d_vel array of particle velocities
 214     \param d_accel array of particle accelerations
 215     \param d_group_members Device array listing the indicies of the mebers of the group to integrate
 216     \param group_size Number of members in the group
 217     \param d_net_force Net force on each particle
 218     \param block_size Size of the block to execute on the device
 219     \param deltaT Amount of real time to step forward in one time step
 220 */
 221 cudaError_t gpu_nvt_mtk_step_two(Scalar4 *d_vel,
 222                              Scalar3 *d_accel,
 223                              unsigned int *d_group_members,
 224                              unsigned int group_size,
 225                              Scalar4 *d_net_force,
 226                              unsigned int block_size,
 227                              Scalar deltaT)
 228     {
 229     static unsigned int max_block_size = UINT_MAX;
 230     if (max_block_size == UINT_MAX)
 231         {
 232         cudaFuncAttributes attr;
 233         cudaFuncGetAttributes(&attr, (const void *)gpu_nvt_mtk_step_two_kernel);
 234         max_block_size = attr.maxThreadsPerBlock;
 235         }
 236
 237     unsigned int run_block_size = min(block_size, max_block_size);
 238
 239     // setup the grid to run the kernel
 240     dim3 grid( (group_size/run_block_size) + 1, 1, 1);
 241     dim3 threads(run_block_size, 1, 1);
 242
 243     // run the kernel
 244     gpu_nvt_mtk_step_two_kernel<<< grid, threads >>>(d_vel, d_accel, d_group_members, group_size, d_net_force, deltaT);
 245
 246     return cudaSuccess;
 247     }
 248
 249 // vim:syntax=cpp