src/gromacs/nbnxm/grid.h

   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
  36 /*! \internal \file
  37  *
  38  * \brief
  39  * Declares the Grid class.
  40  *
  41  * This class provides functionality for setting up and accessing atoms
  42  * on a grid for one domain decomposition zone. This grid is used for
  43  * generating cluster pair lists for computing non-bonded pair interactions.
  44  * The grid consists of a regular array of columns along dimensions x and y.
  45  * Along z the number of cells and their boundaries vary between the columns.
  46  * Each cell can hold one or more clusters of atoms, depending on the grid
  47  * geometry, which is set by the pair-list type.
  48  *
  49  * \author Berk Hess <hess@kth.se>
  50  * \ingroup module_nbnxm
  51  */
  52
  53 #ifndef GMX_NBNXM_GRID_H
  54 #define GMX_NBNXM_GRID_H
  55
  56 #include <memory>
  57 #include <vector>
  58
  59 #include "gromacs/gpu_utils/hostallocator.h"
  60 #include "gromacs/math/vectypes.h"
  61 #include "gromacs/utility/alignedallocator.h"
  62 #include "gromacs/utility/arrayref.h"
  63 #include "gromacs/utility/range.h"
  64
  65 struct gmx_domdec_zones_t;
  66 struct nbnxn_atomdata_t;
  67 struct nbnxn_search;
  68 enum class PairlistType;
  69
  70 namespace gmx
  71 {
  72 class UpdateGroupsCog;
  73 } // namespace gmx
  74
  75 namespace Nbnxm
  76 {
  77
  78 struct GridSetData;
  79 struct GridWork;
  80
  81 /*! \internal
  82  * \brief Bounding box for a nbnxm atom cluster
  83  *
  84  * \note Should be aligned in memory to enable 4-wide SIMD operations.
  85  */
  86 struct BoundingBox
  87 {
  88     /*! \internal
  89      * \brief Corner for the bounding box, padded with one element to enable 4-wide SIMD operations
  90      */
  91     struct Corner
  92     {
  93         //! Returns a corner with the minimum coordinates along each dimension
  94         static Corner min(const Corner &c1,
  95                           const Corner &c2)
  96         {
  97             Corner cMin;
  98
  99             cMin.x       = std::min(c1.x, c2.x);
 100             cMin.y       = std::min(c1.y, c2.y);
 101             cMin.z       = std::min(c1.z, c2.z);
 102             /* This value of the padding is irrelevant, as long as it
 103              * is initialized. We use min to allow auto-vectorization.
 104              */
 105             cMin.padding = std::min(c1.padding, c2.padding);
 106
 107             return cMin;
 108         }
 109
 110         //! Returns a corner with the maximum coordinates along each dimension
 111         static Corner max(const Corner &c1,
 112                           const Corner &c2)
 113         {
 114             Corner cMax;
 115
 116             cMax.x       = std::max(c1.x, c2.x);
 117             cMax.y       = std::max(c1.y, c2.y);
 118             cMax.z       = std::max(c1.z, c2.z);
 119             cMax.padding = std::max(c1.padding, c2.padding);
 120
 121             return cMax;
 122         }
 123
 124         //! Returns a pointer for SIMD loading of a Corner object
 125         const float *ptr() const
 126         {
 127             return &x;
 128         }
 129
 130         //! Returns a pointer for SIMD storing of a Corner object
 131         float *ptr()
 132         {
 133             return &x;
 134         }
 135
 136         float x;       //!< x coordinate
 137         float y;       //!< y coordinate
 138         float z;       //!< z coordinate
 139         float padding; //!< padding, unused, but should be set to avoid operations on unitialized data
 140     };
 141
 142     Corner lower; //!< lower, along x and y and z, corner
 143     Corner upper; //!< upper, along x and y and z, corner
 144 };
 145
 146 /*! \internal
 147  * \brief Bounding box for one dimension of a grid cell
 148  */
 149 struct BoundingBox1D
 150 {
 151     float lower; //!< lower bound
 152     float upper; //!< upper bound
 153 };
 154
 155 } // namespace Nbnxm
 156
 157 namespace Nbnxm
 158 {
 159
 160 /*! \internal
 161  * \brief A pair-search grid object for one domain decomposition zone
 162  *
 163  * This is a rectangular 3D grid covering a potentially non-rectangular
 164  * volume which is either the whole unit cell or the local zone or part
 165  * of a non-local zone when using domain decomposition. The grid cells
 166  * are even spaced along x/y and irregular along z. Each cell is sub-divided
 167  * into atom clusters. With a CPU geometry, each cell contains 1 or 2 clusters.
 168  * With a GPU geometry, each cell contains up to 8 clusters. The geometry is
 169  * set by the pairlist type which is the only argument of the constructor.
 170  *
 171  * When multiple grids are used, i.e. with domain decomposition, we want
 172  * to avoid the overhead of multiple coordinate arrays or extra indexing.
 173  * Therefore each grid stores a cell offset, so a contiguous cell index
 174  * can be used to index atom arrays. All methods returning atom indices
 175  * return indices which index into a full atom array.
 176  *
 177  * Note that when atom groups, instead of individual atoms, are assigned
 178  * to grid cells, individual atoms can be geometrically outside the cell
 179  * and grid that they have been assigned to (as determined by the center
 180  * or geometry of the atom group they belong to).
 181  */
 182 class Grid
 183 {
 184     public:
 185         /*! \internal
 186          * \brief The cluster and cell geometry of a grid
 187          */
 188         struct Geometry
 189         {
 190             //! Constructs the cluster/cell geometry given the type of pairlist
 191             Geometry(PairlistType pairlistType);
 192
 193             bool isSimple;             //!< Is this grid simple (CPU) or hierarchical (GPU)
 194             int  numAtomsICluster;     //!< Number of atoms per cluster
 195             int  numAtomsJCluster;     //!< Number of atoms for list j-clusters
 196             int  numAtomsPerCell;      //!< Number of atoms per cell
 197             int  numAtomsICluster2Log; //!< 2log of na_c
 198         };
 199
 200         // The physical dimensions of a grid
 201         struct Dimensions
 202         {
 203             //! The lower corner of the (local) grid
 204             rvec lowerCorner;
 205             //! The upper corner of the (local) grid
 206             rvec upperCorner;
 207             //! The physical grid size: upperCorner - lowerCorner
 208             rvec gridSize;
 209             //! An estimate for the atom number density of the region targeted by the grid
 210             real atomDensity;
 211             //! The maximum distance an atom can be outside of a cell and outside of the grid
 212             real maxAtomGroupRadius;
 213             //! Size of cell along dimension x and y
 214             real cellSize[DIM - 1];
 215             //! 1/size of a cell along dimensions x and y
 216             real invCellSize[DIM - 1];
 217             //! The number of grid cells along dimensions x and y
 218             int  numCells[DIM - 1];
 219         };
 220
 221         //! Constructs a grid given the type of pairlist
 222         Grid(PairlistType  pairlistType,
 223              const bool   &haveFep);
 224
 225         //! Returns the geometry of the grid cells
 226         const Geometry &geometry() const
 227         {
 228             return geometry_;
 229         }
 230
 231         //! Returns the dimensions of the grid
 232         const Dimensions &dimensions() const
 233         {
 234             return dimensions_;
 235         }
 236
 237         //! Returns the total number of grid columns
 238         int numColumns() const
 239         {
 240             return dimensions_.numCells[XX]*dimensions_.numCells[YY];
 241         }
 242
 243         //! Returns the total number of grid cells
 244         int numCells() const
 245         {
 246             return numCellsTotal_;
 247         }
 248
 249         //! Returns the cell offset of (the first cell of) this grid in the list of cells combined over all grids
 250         int cellOffset() const
 251         {
 252             return cellOffset_;
 253         }
 254
 255         //! Returns the maximum number of grid cells in a column
 256         int numCellsColumnMax() const
 257         {
 258             return numCellsColumnMax_;
 259         }
 260
 261         //! Returns the start of the source atom range mapped to this grid
 262         int srcAtomBegin() const
 263         {
 264             return srcAtomBegin_;
 265         }
 266
 267         //! Returns the end of the source atom range mapped to this grid
 268         int srcAtomEnd() const
 269         {
 270             return srcAtomEnd_;
 271         }
 272
 273         //! Returns the first cell index in the grid, starting at 0 in this grid
 274         int firstCellInColumn(int columnIndex) const
 275         {
 276             return cxy_ind_[columnIndex];
 277         }
 278
 279         //! Returns the number of cells in the column
 280         int numCellsInColumn(int columnIndex) const
 281         {
 282             return cxy_ind_[columnIndex + 1] - cxy_ind_[columnIndex];
 283         }
 284
 285         //! Returns the index of the first atom in the column
 286         int firstAtomInColumn(int columnIndex) const
 287         {
 288             return (cellOffset_ + cxy_ind_[columnIndex])*geometry_.numAtomsPerCell;
 289         }
 290
 291         //! Returns the number of real atoms in the column
 292         int numAtomsInColumn(int columnIndex) const
 293         {
 294             return cxy_na_[columnIndex];
 295         }
 296
 297         /*! \brief Returns a view of the number of non-filler, atoms for each grid column
 298          *
 299          * \todo Needs a useful name. */
 300         gmx::ArrayRef<const int> cxy_na() const
 301         {
 302             return cxy_na_;
 303         }
 304         /*! \brief Returns a view of the grid-local cell index for each grid column
 305          *
 306          * \todo Needs a useful name. */
 307         gmx::ArrayRef<const int> cxy_ind() const
 308         {
 309             return cxy_ind_;
 310         }
 311
 312         //! Returns the number of real atoms in the column
 313         int numAtomsPerCell() const
 314         {
 315             return geometry_.numAtomsPerCell;
 316         }
 317
 318         //! Returns the number of atoms in the column including padding
 319         int paddedNumAtomsInColumn(int columnIndex) const
 320         {
 321             return numCellsInColumn(columnIndex)*geometry_.numAtomsPerCell;
 322         }
 323
 324         //! Returns the end of the atom index range on the grid, including padding
 325         int atomIndexEnd() const
 326         {
 327             return (cellOffset_ + numCellsTotal_)*geometry_.numAtomsPerCell;
 328         }
 329
 330         //! Returns whether any atom in the cluster is perturbed
 331         bool clusterIsPerturbed(int clusterIndex) const
 332         {
 333             return fep_[clusterIndex] != 0U;
 334         }
 335
 336         //! Returns whether the given atom in the cluster is perturbed
 337         bool atomIsPerturbed(int clusterIndex,
 338                              int atomIndexInCluster) const
 339         {
 340             return (fep_[clusterIndex] & (1 << atomIndexInCluster)) != 0U;
 341         }
 342
 343         //! Returns the free-energy perturbation bits for the cluster
 344         unsigned int fepBits(int clusterIndex) const
 345         {
 346             return fep_[clusterIndex];
 347         }
 348
 349         //! Returns the i-bounding boxes for all clusters on the grid
 350         gmx::ArrayRef<const BoundingBox> iBoundingBoxes() const
 351         {
 352             return bb_;
 353         }
 354
 355         //! Returns the j-bounding boxes for all clusters on the grid
 356         gmx::ArrayRef<const BoundingBox> jBoundingBoxes() const
 357         {
 358             return bbj_;
 359         }
 360
 361         //! Returns the packed bounding boxes for all clusters on the grid, empty with a CPU list
 362         gmx::ArrayRef<const float> packedBoundingBoxes() const
 363         {
 364             return pbb_;
 365         }
 366
 367         //! Returns the bounding boxes along z for all cells on the grid
 368         gmx::ArrayRef<const BoundingBox1D> zBoundingBoxes() const
 369         {
 370             return bbcz_;
 371         }
 372
 373         //! Returns the flags for all clusters on the grid
 374         gmx::ArrayRef<const int> clusterFlags() const
 375         {
 376             return flags_;
 377         }
 378
 379         //! Returns the number of clusters for all cells on the grid, empty with a CPU geometry
 380         gmx::ArrayRef<const int> numClustersPerCell() const
 381         {
 382             return numClusters_;
 383         }
 384
 385         //! Returns the cluster index for an atom
 386         int atomToCluster(int atomIndex) const
 387         {
 388             return (atomIndex >> geometry_.numAtomsICluster2Log);
 389         }
 390
 391         //! Returns the total number of clusters on the grid
 392         int numClusters() const
 393         {
 394             if (geometry_.isSimple)
 395             {
 396                 return numCellsTotal_;
 397             }
 398             else
 399             {
 400                 return numClustersTotal_;
 401             }
 402         }
 403
 404         //! Sets the grid dimensions
 405         void setDimensions(int                   ddZone,
 406                            int                   numAtoms,
 407                            gmx::RVec             lowerCorner,
 408                            gmx::RVec             upperCorner,
 409                            real                  atomDensity,
 410                            real                  maxAtomGroupRadius,
 411                            bool                  haveFep,
 412                            gmx::PinningPolicy    pinningPolicy);
 413
 414         //! Sets the cell indices using indices in \p gridSetData and \p gridWork
 415         void setCellIndices(int                             ddZone,
 416                             int                             cellOffset,
 417                             GridSetData                    *gridSetData,
 418                             gmx::ArrayRef<GridWork>         gridWork,
 419                             gmx::Range<int>                 atomRange,
 420                             const int                      *atinfo,
 421                             gmx::ArrayRef<const gmx::RVec>  x,
 422                             int                             numAtomsMoved,
 423                             nbnxn_atomdata_t               *nbat);
 424
 425         //! Determine in which grid columns atoms should go, store cells and atom counts in \p cell and \p cxy_na
 426         static void calcColumnIndices(const Grid::Dimensions         &gridDims,
 427                                       const gmx::UpdateGroupsCog     *updateGroupsCog,
 428                                       gmx::Range<int>                 atomRange,
 429                                       gmx::ArrayRef<const gmx::RVec>  x,
 430                                       int                             dd_zone,
 431                                       const int                      *move,
 432                                       int                             thread,
 433                                       int                             nthread,
 434                                       gmx::ArrayRef<int>              cell,
 435                                       gmx::ArrayRef<int>              cxy_na);
 436
 437     private:
 438         /*! \brief Fill a pair search cell with atoms
 439          *
 440          * Potentially sorts atoms and sets the interaction flags.
 441          */
 442         void fillCell(GridSetData                    *gridSetData,
 443                       nbnxn_atomdata_t               *nbat,
 444                       int                             atomStart,
 445                       int                             atomEnd,
 446                       const int                      *atinfo,
 447                       gmx::ArrayRef<const gmx::RVec>  x,
 448                       BoundingBox gmx_unused         *bb_work_aligned);
 449
 450         //! Spatially sort the atoms within the given column range, for CPU geometry
 451         void sortColumnsCpuGeometry(GridSetData                   *gridSetData,
 452                                     int                            dd_zone,
 453                                     const int                     *atinfo,
 454                                     gmx::ArrayRef<const gmx::RVec> x,
 455                                     nbnxn_atomdata_t              *nbat,
 456                                     gmx::Range<int>                columnRange,
 457                                     gmx::ArrayRef<int>             sort_work);
 458
 459         //! Spatially sort the atoms within the given column range, for GPU geometry
 460         void sortColumnsGpuGeometry(GridSetData                   *gridSetData,
 461                                     int                            dd_zone,
 462                                     const int                     *atinfo,
 463                                     gmx::ArrayRef<const gmx::RVec> x,
 464                                     nbnxn_atomdata_t              *nbat,
 465                                     gmx::Range<int>                columnRange,
 466                                     gmx::ArrayRef<int>             sort_work);
 467
 468         /* Data members */
 469         //! The geometry of the grid clusters and cells
 470         Geometry   geometry_;
 471         //! The physical dimensions of the grid
 472         Dimensions dimensions_;
 473
 474         //! The total number of cells in this grid
 475         int        numCellsTotal_;
 476         //! Index in nbs->cell corresponding to cell 0
 477         int        cellOffset_;
 478         //! The maximum number of cells in a column
 479         int        numCellsColumnMax_;
 480
 481         //! The start of the source atom range mapped to this grid
 482         int        srcAtomBegin_;
 483         //! The end of the source atom range mapped to this grid
 484         int        srcAtomEnd_;
 485
 486         /* Grid data */
 487         /*! \brief The number of, non-filler, atoms for each grid column.
 488          *
 489          * \todo Needs a useful name. */
 490         gmx::HostVector<int>    cxy_na_;
 491         /*! \brief The grid-local cell index for each grid column
 492          *
 493          * \todo Needs a useful name. */
 494         gmx::HostVector<int>    cxy_ind_;
 495
 496         //! The number of cluster for each cell
 497         std::vector<int> numClusters_;
 498
 499         /* Bounding boxes */
 500         //! Bounding boxes in z for the cells
 501         std::vector<BoundingBox1D>                          bbcz_;
 502         //! 3D bounding boxes for the sub cells
 503         std::vector < BoundingBox, gmx::AlignedAllocator < BoundingBox>> bb_;
 504         //! 3D j-bounding boxes for the case where the i- and j-cluster sizes are different
 505         std::vector < BoundingBox, gmx::AlignedAllocator < BoundingBox>> bbjStorage_;
 506         //! 3D j-bounding boxes
 507         gmx::ArrayRef<BoundingBox>                           bbj_;
 508         //! 3D bounding boxes in packed xxxx format per cell
 509         std::vector < float, gmx::AlignedAllocator < float>> pbb_;
 510
 511         //! Tells whether we have perturbed interactions, authorative source is in GridSet (never modified)
 512         const bool               &haveFep_;
 513
 514         /* Bit-flag information */
 515         //! Flags for properties of clusters in each cell
 516         std::vector<int>          flags_;
 517         //! Signal bits for atoms in each cell that tell whether an atom is perturbed
 518         std::vector<unsigned int> fep_;
 519
 520         /* Statistics */
 521         //! Total number of clusters, used for printing
 522         int numClustersTotal_;
 523 };
 524
 525 } // namespace Nbnxm
 526
 527 #endif