2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * Declares the Grid class.
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.
49 * \author Berk Hess <hess@kth.se>
50 * \ingroup module_nbnxm
53 #ifndef GMX_NBNXM_GRID_H
54 #define GMX_NBNXM_GRID_H
59 #include "gromacs/math/vectypes.h"
60 #include "gromacs/simd/simd.h"
61 #include "gromacs/utility/alignedallocator.h"
62 #include "gromacs/utility/arrayref.h"
65 struct gmx_domdec_zones_t;
66 struct nbnxn_atomdata_t;
68 enum class PairlistType;
72 class UpdateGroupsCog;
79 * \brief Bounding box for a nbnxm atom cluster
81 * \note Should be aligned in memory to enable 4-wide SIMD operations.
86 * \brief Corner for the bounding box, padded with one element to enable 4-wide SIMD operations
90 //! Returns a corner with the minimum coordinates along each dimension
91 static Corner min(const Corner &c1,
96 cMin.x = std::min(c1.x, c2.x);
97 cMin.y = std::min(c1.y, c2.y);
98 cMin.z = std::min(c1.z, c2.z);
99 /* This value of the padding is irrelevant, as long as it
100 * is initialized. We use min to allow auto-vectorization.
102 cMin.padding = std::min(c1.padding, c2.padding);
107 //! Returns a corner with the maximum coordinates along each dimension
108 static Corner max(const Corner &c1,
113 cMax.x = std::max(c1.x, c2.x);
114 cMax.y = std::max(c1.y, c2.y);
115 cMax.z = std::max(c1.z, c2.z);
116 cMax.padding = std::max(c1.padding, c2.padding);
121 //! Returns a pointer for SIMD loading of a Corner object
122 const float *ptr() const
127 //! Returns a pointer for SIMD storing of a Corner object
133 float x; //!< x coordinate
134 float y; //!< y coordinate
135 float z; //!< z coordinate
136 float padding; //!< padding, unused, but should be set to avoid operations on unitialized data
139 Corner lower; //!< lower, along x and y and z, corner
140 Corner upper; //!< upper, along x and y and z, corner
146 // TODO: Convert macros to constexpr int
148 /* Pair search box lower and upper bound in z only. */
151 /* Bounding box calculations are (currently) always in single precision, so
152 * we only need to check for single precision support here.
153 * This uses less (cache-)memory and SIMD is faster, at least on x86.
155 #if GMX_SIMD4_HAVE_FLOAT
156 # define NBNXN_SEARCH_BB_SIMD4 1
157 /* Memory alignment in bytes as required by SIMD aligned loads/stores */
158 # define NBNXN_SEARCH_BB_MEM_ALIGN (GMX_SIMD4_WIDTH*sizeof(float))
160 # define NBNXN_SEARCH_BB_SIMD4 0
161 /* No alignment required, but set it so we can call the same routines */
162 # define NBNXN_SEARCH_BB_MEM_ALIGN 32
166 #if NBNXN_SEARCH_BB_SIMD4
167 /* Always use 4-wide SIMD for bounding box calculations */
170 /* Single precision BBs + coordinates, we can also load coordinates with SIMD */
171 # define NBNXN_SEARCH_SIMD4_FLOAT_X_BB 1
173 # define NBNXN_SEARCH_SIMD4_FLOAT_X_BB 0
176 /* The packed bounding box coordinate stride is always set to 4.
177 * With AVX we could use 8, but that turns out not to be faster.
179 # define STRIDE_PBB GMX_SIMD4_WIDTH
180 # define STRIDE_PBB_2LOG 2
182 /* Store bounding boxes corners as quadruplets: xxxxyyyyzzzz */
183 # define NBNXN_BBXXXX 1
184 /* Size of a quadruplet of bounding boxes, each 2 corners, stored packed */
185 # define NNBSBB_XXXX (2*DIM*STRIDE_PBB)
187 #else /* NBNXN_SEARCH_BB_SIMD4 */
189 # define NBNXN_SEARCH_SIMD4_FLOAT_X_BB 0
190 # define NBNXN_BBXXXX 0
192 #endif /* NBNXN_SEARCH_BB_SIMD4 */
198 * \brief A pair-search grid object for one domain decomposition zone
200 * This is a rectangular 3D grid covering a potentially non-rectangular
201 * volume which is either the whole unit cell or the local zone or part
202 * of a non-local zone when using domain decomposition. The grid cells
203 * are even spaced along x/y and irregular along z. Each cell is sub-divided
204 * into atom clusters. With a CPU geometry, each cell contains 1 or 2 clusters.
205 * With a GPU geometry, each cell contains up to 8 clusters. The geometry is
206 * set by the pairlist type which is the only argument of the constructor.
208 * When multiple grids are used, i.e. with domain decomposition, we want
209 * to avoid the overhead of multiple coordinate arrays or extra indexing.
210 * Therefore each grid stores a cell offset, so a contiguous cell index
211 * can be used to index atom arrays. All methods returning atom indices
212 * return indices which index into a full atom array.
214 * Note that when atom groups, instead of individual atoms, are assigned
215 * to grid cells, individual atoms can be geometrically outside the cell
216 * and grid that they have been assigned to (as determined by the center
217 * or geometry of the atom group they belong to).
223 * \brief The cluster and cell geometry of a grid
227 //! Constructs the cluster/cell geometry given the type of pairlist
228 Geometry(PairlistType pairlistType);
230 bool isSimple; //!< Is this grid simple (CPU) or hierarchical (GPU)
231 int numAtomsICluster; //!< Number of atoms per cluster
232 int numAtomsJCluster; //!< Number of atoms for list j-clusters
233 int numAtomsPerCell; //!< Number of atoms per cell
234 int numAtomsICluster2Log; //!< 2log of na_c
237 // The physical dimensions of a grid
240 //! The lower corner of the (local) grid
242 //! The upper corner of the (local) grid
244 //! The physical grid size: upperCorner - lowerCorner
246 //! An estimate for the atom number density of the region targeted by the grid
248 //! The maximum distance an atom can be outside of a cell and outside of the grid
249 real maxAtomGroupRadius;
250 //! Size of cell along dimension x and y
251 real cellSize[DIM - 1];
252 //! 1/size of a cell along dimensions x and y
253 real invCellSize[DIM - 1];
254 //! The number of grid cells along dimensions x and y
255 int numCells[DIM - 1];
258 //! Constructs a grid given the type of pairlist
259 Grid(PairlistType pairlistType);
261 //! Returns the geometry of the grid cells
262 const Geometry &geometry() const
267 //! Returns the dimensions of the grid
268 const Dimensions &dimensions() const
273 //! Returns the total number of grid columns
274 int numColumns() const
276 return dimensions_.numCells[XX]*dimensions_.numCells[YY];
279 //! Returns the total number of grid cells
282 return numCellsTotal_;
285 //! Returns the cell offset of (the first cell of) this grid in the list of cells combined over all grids
286 int cellOffset() const
291 //! Returns the first cell index in the grid, starting at 0 in this grid
292 int firstCellInColumn(int columnIndex) const
294 return cxy_ind_[columnIndex];
297 //! Returns the number of cells in the column
298 int numCellsInColumn(int columnIndex) const
300 return cxy_ind_[columnIndex + 1] - cxy_ind_[columnIndex];
303 //! Returns the index of the first atom in the column
304 int firstAtomInColumn(int columnIndex) const
306 return (cellOffset_ + cxy_ind_[columnIndex])*geometry_.numAtomsPerCell;
309 //! Returns the number of real atoms in the column
310 int numAtomsInColumn(int columnIndex) const
312 return cxy_na_[columnIndex];
315 //! Returns the number of atoms in the column including padding
316 int paddedNumAtomsInColumn(int columnIndex) const
318 return numCellsInColumn(columnIndex)*geometry_.numAtomsPerCell;
321 //! Returns the end of the atom index range on the grid, including padding
322 int atomIndexEnd() const
324 return (cellOffset_ + numCellsTotal_)*geometry_.numAtomsPerCell;
327 //! Returns whether any atom in the cluster is perturbed
328 bool clusterIsPerturbed(int clusterIndex) const
330 return fep_[clusterIndex] != 0u;
333 //! Returns whether the given atom in the cluster is perturbed
334 bool atomIsPerturbed(int clusterIndex,
335 int atomIndexInCluster) const
337 return (fep_[clusterIndex] & (1 << atomIndexInCluster)) != 0u;
340 //! Returns the free-energy perturbation bits for the cluster
341 unsigned int fepBits(int clusterIndex) const
343 return fep_[clusterIndex];
346 //! Returns the i-bounding boxes for all clusters on the grid
347 gmx::ArrayRef<const BoundingBox> iBoundingBoxes() const
352 //! Returns the j-bounding boxes for all clusters on the grid
353 gmx::ArrayRef<const BoundingBox> jBoundingBoxes() const
358 //! Returns the packed bounding boxes for all clusters on the grid, empty with a CPU list
359 gmx::ArrayRef<const float> packedBoundingBoxes() const
364 //! Returns the bounding boxes along z for all cells on the grid
365 gmx::ArrayRef<const float> zBoundingBoxes() const
370 //! Returns the flags for all clusters on the grid
371 gmx::ArrayRef<const int> clusterFlags() const
376 //! Returns the number of clusters for all cells on the grid, empty with a CPU geometry
377 gmx::ArrayRef<const int> numClustersPerCell() const
382 //! Returns the cluster index for an atom
383 int atomToCluster(int atomIndex) const
385 return (atomIndex >> geometry_.numAtomsICluster2Log);
388 //! Returns the total number of clusters on the grid
389 int numClusters() const
391 if (geometry_.isSimple)
393 return numCellsTotal_;
397 return numClustersTotal_;
401 //! Sets the grid dimensions
402 void setDimensions(const nbnxn_search *nbs,
405 const rvec lowerCorner,
406 const rvec upperCorner,
408 real maxAtomGroupRadius);
410 //! Determine in which grid cells the atoms should go
411 void calcCellIndices(nbnxn_search *nbs,
414 const gmx::UpdateGroupsCog *updateGroupsCog,
418 gmx::ArrayRef<const gmx::RVec> x,
421 nbnxn_atomdata_t *nbat);
424 /*! \brief Fill a pair search cell with atoms
426 * Potentially sorts atoms and sets the interaction flags.
428 void fillCell(nbnxn_search *nbs,
429 nbnxn_atomdata_t *nbat,
433 gmx::ArrayRef<const gmx::RVec> x,
434 BoundingBox gmx_unused *bb_work_aligned);
436 //! Spatially sort the atoms within one grid column
437 void sortColumnsCpuGeometry(nbnxn_search *nbs,
439 int atomStart, int atomEnd,
441 gmx::ArrayRef<const gmx::RVec> x,
442 nbnxn_atomdata_t *nbat,
443 int cxy_start, int cxy_end,
444 gmx::ArrayRef<int> sort_work);
446 //! Spatially sort the atoms within one grid column
447 void sortColumnsGpuGeometry(nbnxn_search *nbs,
449 int atomStart, int atomEnd,
451 gmx::ArrayRef<const gmx::RVec> x,
452 nbnxn_atomdata_t *nbat,
453 int cxy_start, int cxy_end,
454 gmx::ArrayRef<int> sort_work);
457 //! The geometry of the grid clusters and cells
459 //! The physical dimensions of the grid
460 Dimensions dimensions_;
462 //! The total number of cells in this grid
464 //! Index in nbs->cell corresponding to cell 0 */
468 //! The number of, non-filler, atoms for each grid column
469 std::vector<int> cxy_na_;
470 //! The grid-local cell index for each grid column
471 std::vector<int> cxy_ind_;
473 //! The number of cluster for each cell
474 std::vector<int> numClusters_;
477 //! Bounding boxes in z for the cells
478 std::vector<float> bbcz_;
479 //! 3D bounding boxes for the sub cells
480 std::vector < BoundingBox, gmx::AlignedAllocator < BoundingBox>> bb_;
481 //! 3D j-bounding boxes for the case where the i- and j-cluster sizes are different
482 std::vector < BoundingBox, gmx::AlignedAllocator < BoundingBox>> bbjStorage_;
483 //! 3D j-bounding boxes
484 gmx::ArrayRef<BoundingBox> bbj_;
485 //! 3D bounding boxes in packed xxxx format per cell
486 std::vector < float, gmx::AlignedAllocator < float>> pbb_;
488 /* Bit-flag information */
489 //! Flags for properties of clusters in each cell
490 std::vector<int> flags_;
491 //! Signal bits for atoms in each cell that tell whether an atom is perturbed
492 std::vector<unsigned int> fep_;
495 //! Total number of clusters, used for printing
496 int numClustersTotal_;