From 893eb8e85ad94831eed37cbb8aa10930c014b0fe Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Mon, 21 May 2018 22:42:49 +0200 Subject: [PATCH] Make nbnxn_grid_t counts and sizes arrays This enables templating on the dimension. Also put the cell index storing and atom counting together in a separate function, since these need to be consistent. Change-Id: I9931ff258a90a6ea5af788b7ec3e56d32859e47b --- src/gromacs/mdlib/nbnxn_atomdata.cpp | 24 +++--- src/gromacs/mdlib/nbnxn_grid.cpp | 156 ++++++++++++++++++----------------- src/gromacs/mdlib/nbnxn_internal.h | 39 ++++----- src/gromacs/mdlib/nbnxn_search.cpp | 84 ++++++++++--------- 4 files changed, 155 insertions(+), 148 deletions(-) diff --git a/src/gromacs/mdlib/nbnxn_atomdata.cpp b/src/gromacs/mdlib/nbnxn_atomdata.cpp index fa5a923cb2..6cb43f2744 100644 --- a/src/gromacs/mdlib/nbnxn_atomdata.cpp +++ b/src/gromacs/mdlib/nbnxn_atomdata.cpp @@ -798,7 +798,7 @@ static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t *nbat, for (const nbnxn_grid_t &grid : nbs->grid) { /* Loop over all columns and copy and fill */ - for (int i = 0; i < grid.ncx*grid.ncy; i++) + for (int i = 0; i < grid.numCells[XX]*grid.numCells[YY]; i++) { int ncz = grid.cxy_ind[i+1] - grid.cxy_ind[i]; int ash = (grid.cell0 + grid.cxy_ind[i])*grid.na_sc; @@ -818,7 +818,7 @@ static void nbnxn_atomdata_set_ljcombparams(nbnxn_atomdata_t *nbat, for (const nbnxn_grid_t &grid : nbs->grid) { /* Loop over all columns and copy and fill */ - for (int i = 0; i < grid.ncx*grid.ncy; i++) + for (int i = 0; i < grid.numCells[XX]*grid.numCells[YY]; i++) { int ncz = grid.cxy_ind[i+1] - grid.cxy_ind[i]; int ash = (grid.cell0 + grid.cxy_ind[i])*grid.na_sc; @@ -857,7 +857,7 @@ static void nbnxn_atomdata_set_charges(nbnxn_atomdata_t *nbat, for (const nbnxn_grid_t &grid : nbs->grid) { /* Loop over all columns and copy and fill */ - for (int cxy = 0; cxy < grid.ncx*grid.ncy; cxy++) + for (int cxy = 0; cxy < grid.numCells[XX]*grid.numCells[YY]; cxy++) { int ash = (grid.cell0 + grid.cxy_ind[cxy])*grid.na_sc; int na = grid.cxy_na[cxy]; @@ -1000,7 +1000,7 @@ static void nbnxn_atomdata_set_energygroups(nbnxn_atomdata_t *nbat, for (const nbnxn_grid_t &grid : nbs->grid) { /* Loop over all columns and copy and fill */ - for (int i = 0; i < grid.ncx*grid.ncy; i++) + for (int i = 0; i < grid.numCells[XX]*grid.numCells[YY]; i++) { int ncz = grid.cxy_ind[i+1] - grid.cxy_ind[i]; int ash = (grid.cell0 + grid.cxy_ind[i])*grid.na_sc; @@ -1091,25 +1091,23 @@ void nbnxn_atomdata_copy_x_to_nbat_x(const nbnxn_search_t nbs, { for (int g = g0; g < g1; g++) { - const nbnxn_grid_t *grid; - int cxy0, cxy1; + const nbnxn_grid_t &grid = nbs->grid[g]; + const int numCellsXY = grid.numCells[XX]*grid.numCells[YY]; - grid = &nbs->grid[g]; - - cxy0 = (grid->ncx*grid->ncy* th +nth-1)/nth; - cxy1 = (grid->ncx*grid->ncy*(th+1)+nth-1)/nth; + const int cxy0 = (numCellsXY* th + nth - 1)/nth; + const int cxy1 = (numCellsXY*(th + 1) + nth - 1)/nth; for (int cxy = cxy0; cxy < cxy1; cxy++) { int na, ash, na_fill; - na = grid->cxy_na[cxy]; - ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc; + na = grid.cxy_na[cxy]; + ash = (grid.cell0 + grid.cxy_ind[cxy])*grid.na_sc; if (g == 0 && FillLocal) { na_fill = - (grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy])*grid->na_sc; + (grid.cxy_ind[cxy+1] - grid.cxy_ind[cxy])*grid.na_sc; } else { diff --git a/src/gromacs/mdlib/nbnxn_grid.cpp b/src/gromacs/mdlib/nbnxn_grid.cpp index 821ffb5710..13d1c75e16 100644 --- a/src/gromacs/mdlib/nbnxn_grid.cpp +++ b/src/gromacs/mdlib/nbnxn_grid.cpp @@ -119,19 +119,20 @@ static void set_grid_size_xy(const nbnxn_search_t nbs, * in the nbsist when the fixed cell dimensions (x,y) are * larger than the variable one (z) than the other way around. */ - grid->ncx = std::max(1, static_cast(size[XX]/tlen_x)); - grid->ncy = std::max(1, static_cast(size[YY]/tlen_y)); + grid->numCells[XX] = std::max(1, static_cast(size[XX]/tlen_x)); + grid->numCells[YY] = std::max(1, static_cast(size[YY]/tlen_y)); } else { - grid->ncx = 1; - grid->ncy = 1; + grid->numCells[XX] = 1; + grid->numCells[YY] = 1; } - grid->sx = size[XX]/grid->ncx; - grid->sy = size[YY]/grid->ncy; - grid->inv_sx = 1/grid->sx; - grid->inv_sy = 1/grid->sy; + for (int d = 0; d < DIM - 1; d++) + { + grid->cellSize[d] = size[d]/grid->numCells[d]; + grid->invCellSize[d] = 1/grid->cellSize[d]; + } if (ddZone > 0) { @@ -141,28 +142,29 @@ static void set_grid_size_xy(const nbnxn_search_t nbs, * they end up on the grid, but for performance it's better * if they don't end up in cells that can be within cut-off range. */ - grid->ncx++; - grid->ncy++; + grid->numCells[XX]++; + grid->numCells[YY]++; } /* We need one additional cell entry for particles moved by DD */ - grid->cxy_na.resize(grid->ncx*grid->ncy + 1); - grid->cxy_ind.resize(grid->ncx*grid->ncy + 2); + int numCellsXY = grid->numCells[XX]*grid->numCells[YY]; + grid->cxy_na.resize(numCellsXY + 1); + grid->cxy_ind.resize(numCellsXY + 2); for (nbnxn_search_work_t &work : nbs->work) { - work.cxy_na.resize(grid->ncx*grid->ncy + 1); + work.cxy_na.resize(numCellsXY + 1); } /* Worst case scenario of 1 atom in each last cell */ int maxNumCells; if (grid->na_cj <= grid->na_c) { - maxNumCells = numAtoms/grid->na_sc + grid->ncx*grid->ncy; + maxNumCells = numAtoms/grid->na_sc + numCellsXY; } else { - maxNumCells = numAtoms/grid->na_sc + grid->ncx*grid->ncy*grid->na_cj/grid->na_c; + maxNumCells = numAtoms/grid->na_sc + numCellsXY*grid->na_cj/grid->na_c; } grid->nsubc.resize(maxNumCells); @@ -594,7 +596,7 @@ static void combine_bounding_box_pairs(nbnxn_grid_t *grid, // TODO: During SIMDv2 transition only some archs use namespace (remove when done) using namespace gmx; - for (int i = 0; i < grid->ncx*grid->ncy; i++) + for (int i = 0; i < grid->numCells[XX]*grid->numCells[YY]; i++) { /* Starting bb in a column is expected to be 2-aligned */ int sc2 = grid->cxy_ind[i]>>1; @@ -651,16 +653,16 @@ static void print_bbsizes_simple(FILE *fp, } dsvmul(1.0/grid->nc, ba, ba); + real avgCellSizeZ = (grid->atom_density > 0 ? grid->na_c/(grid->atom_density*grid->cellSize[XX]*grid->cellSize[YY]) : 0.0); + fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n", - grid->sx, - grid->sy, - grid->atom_density > 0 ? - grid->na_c/(grid->atom_density*grid->sx*grid->sy) : 0.0, + grid->cellSize[XX], + grid->cellSize[YY], + avgCellSizeZ, ba[XX], ba[YY], ba[ZZ], - ba[XX]/grid->sx, - ba[YY]/grid->sy, - grid->atom_density > 0 ? - ba[ZZ]/(grid->na_c/(grid->atom_density*grid->sx*grid->sy)) : 0.0); + ba[XX]*grid->invCellSize[XX], + ba[YY]*grid->invCellSize[YY], + grid->atom_density > 0 ? ba[ZZ]/avgCellSizeZ : 0.0); } /* Prints the average bb size, used for debug output */ @@ -702,16 +704,16 @@ static void print_bbsizes_supersub(FILE *fp, } dsvmul(1.0/ns, ba, ba); + real avgClusterSizeZ = (grid->atom_density > 0 ? grid->na_sc/(grid->atom_density*grid->cellSize[XX]*grid->cellSize[YY]*c_gpuNumClusterPerCellZ) : 0.0); + fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n", - grid->sx/c_gpuNumClusterPerCellX, - grid->sy/c_gpuNumClusterPerCellY, - grid->atom_density > 0 ? - grid->na_sc/(grid->atom_density*grid->sx*grid->sy*c_gpuNumClusterPerCellZ) : 0.0, + grid->cellSize[XX]/c_gpuNumClusterPerCellX, + grid->cellSize[YY]/c_gpuNumClusterPerCellY, + avgClusterSizeZ, ba[XX], ba[YY], ba[ZZ], - ba[XX]*c_gpuNumClusterPerCellX/grid->sx, - ba[YY]*c_gpuNumClusterPerCellY/grid->sy, - grid->atom_density > 0 ? - ba[ZZ]/(grid->na_sc/(grid->atom_density*grid->sx*grid->sy*c_gpuNumClusterPerCellZ)) : 0.0); + ba[XX]*c_gpuNumClusterPerCellX*grid->invCellSize[XX], + ba[YY]*c_gpuNumClusterPerCellY*grid->invCellSize[YY], + grid->atom_density > 0 ? ba[ZZ]/avgClusterSizeZ : 0.0); } /* Set non-bonded interaction flags for the current cluster. @@ -1013,8 +1015,8 @@ static void sort_columns_supersub(const nbnxn_search_t nbs, */ for (int cxy = cxy_start; cxy < cxy_end; cxy++) { - int gridX = cxy/grid->ncy; - int gridY = cxy - gridX*grid->ncy; + int gridX = cxy/grid->numCells[YY]; + int gridY = cxy - gridX*grid->numCells[YY]; int numAtomsInColumn = grid->cxy_na[cxy]; int numCellsInColumn = grid->cxy_ind[cxy + 1] - grid->cxy_ind[cxy]; @@ -1056,8 +1058,8 @@ static void sort_columns_supersub(const nbnxn_search_t nbs, /* Sort the atoms along y */ sort_atoms(YY, (sub_z & 1), dd_zone, nbs->a.data() + ash_z, na_z, x, - grid->c0[YY] + gridY*grid->sy, - grid->inv_sy, subdiv_z, + grid->c0[YY] + gridY*grid->cellSize[YY], + grid->invCellSize[YY], subdiv_z, sort_work); } @@ -1071,8 +1073,8 @@ static void sort_columns_supersub(const nbnxn_search_t nbs, /* Sort the atoms along x */ sort_atoms(XX, ((cz*c_gpuNumClusterPerCellY + sub_y) & 1), dd_zone, nbs->a.data() + ash_y, na_y, x, - grid->c0[XX] + gridX*grid->sx, - grid->inv_sx, subdiv_y, + grid->c0[XX] + gridX*grid->cellSize[XX], + grid->invCellSize[XX], subdiv_y, sort_work); } @@ -1096,6 +1098,16 @@ static void sort_columns_supersub(const nbnxn_search_t nbs, } } +/* Sets the cell index in the cell array for atom \p atomIndex and increments the atom count for the grid column */ +static void setCellAndAtomCount(gmx::ArrayRef cell, + int cellIndex, + gmx::ArrayRef cxy_na, + int atomIndex) +{ + cell[atomIndex] = cellIndex; + cxy_na[cellIndex] += 1; +} + /* Determine in which grid column atoms should go */ static void calc_column_indices(nbnxn_grid_t *grid, int atomStart, int atomEnd, @@ -1106,7 +1118,7 @@ static void calc_column_indices(nbnxn_grid_t *grid, gmx::ArrayRef cxy_na) { /* We add one extra cell for particles which moved during DD */ - for (int i = 0; i < grid->ncx*grid->ncy+1; i++) + for (int i = 0; i < grid->numCells[XX]*grid->numCells[YY] + 1; i++) { cxy_na[i] = 0; } @@ -1125,38 +1137,36 @@ static void calc_column_indices(nbnxn_grid_t *grid, * The int cast takes care of the lower bound, * we will explicitly take care of the upper bound. */ - int cx = static_cast((x[i][XX] - grid->c0[XX])*grid->inv_sx); - int cy = static_cast((x[i][YY] - grid->c0[YY])*grid->inv_sy); + int cx = static_cast((x[i][XX] - grid->c0[XX])*grid->invCellSize[XX]); + int cy = static_cast((x[i][YY] - grid->c0[YY])*grid->invCellSize[YY]); #ifndef NDEBUG - if (cx < 0 || cx > grid->ncx || - cy < 0 || cy > grid->ncy) + if (cx < 0 || cx > grid->numCells[XX] || + cy < 0 || cy > grid->numCells[YY]) { gmx_fatal(FARGS, "grid cell cx %d cy %d out of range (max %d %d)\n" "atom %f %f %f, grid->c0 %f %f", - cx, cy, grid->ncx, grid->ncy, + cx, cy, grid->numCells[XX], grid->numCells[YY], x[i][XX], x[i][YY], x[i][ZZ], grid->c0[XX], grid->c0[YY]); } #endif /* Take care of potential rouding issues */ - cx = std::min(cx, grid->ncx - 1); - cy = std::min(cy, grid->ncy - 1); + cx = std::min(cx, grid->numCells[XX] - 1); + cy = std::min(cy, grid->numCells[YY] - 1); /* For the moment cell will contain only the, grid local, * x and y indices, not z. */ - cell[i] = cx*grid->ncy + cy; + setCellAndAtomCount(cell, cx*grid->numCells[YY] + cy, cxy_na, i); } else { /* Put this moved particle after the end of the grid, * so we can process it later without using conditionals. */ - cell[i] = grid->ncx*grid->ncy; + setCellAndAtomCount(cell, grid->numCells[XX]*grid->numCells[YY], cxy_na, i); } - - cxy_na[cell[i]]++; } } else @@ -1164,8 +1174,8 @@ static void calc_column_indices(nbnxn_grid_t *grid, /* Non-home zone */ for (int i = taskAtomStart; i < taskAtomEnd; i++) { - int cx = static_cast((x[i][XX] - grid->c0[XX])*grid->inv_sx); - int cy = static_cast((x[i][YY] - grid->c0[YY])*grid->inv_sy); + int cx = static_cast((x[i][XX] - grid->c0[XX])*grid->invCellSize[XX]); + int cy = static_cast((x[i][YY] - grid->c0[YY])*grid->invCellSize[YY]); /* For non-home zones there could be particles outside * the non-bonded cut-off range, which have been communicated @@ -1174,16 +1184,14 @@ static void calc_column_indices(nbnxn_grid_t *grid, * we put them in an extra row at the border. */ cx = std::max(cx, 0); - cx = std::min(cx, grid->ncx - 1); + cx = std::min(cx, grid->numCells[XX] - 1); cy = std::max(cy, 0); - cy = std::min(cy, grid->ncy - 1); + cy = std::min(cy, grid->numCells[YY] - 1); /* For the moment cell will contain only the, grid local, * x and y indices, not z. */ - cell[i] = cx*grid->ncy + cy; - - cxy_na[cell[i]]++; + setCellAndAtomCount(cell, cx*grid->numCells[YY] + cy, cxy_na, i); } } } @@ -1244,10 +1252,10 @@ static void calc_cell_indices(const nbnxn_search_t nbs, int ncz_max = 0; int ncz = 0; grid->cxy_ind[0] = 0; - for (int i = 0; i < grid->ncx*grid->ncy + 1; i++) + for (int i = 0; i < grid->numCells[XX]*grid->numCells[YY] + 1; i++) { /* We set ncz_max at the beginning of the loop iso at the end - * to skip i=grid->ncx*grid->ncy which are moved particles + * to skip i=grid->ncx*grid->numCells[YY] which are moved particles * that do not need to be ordered on the grid. */ if (ncz > ncz_max) @@ -1269,7 +1277,7 @@ static void calc_cell_indices(const nbnxn_search_t nbs, /* Clear cxy_na, so we can reuse the array below */ grid->cxy_na[i] = 0; } - grid->nc = grid->cxy_ind[grid->ncx*grid->ncy] - grid->cxy_ind[0]; + grid->nc = grid->cxy_ind[grid->numCells[XX]*grid->numCells[YY]] - grid->cxy_ind[0]; resizeForNumberOfCells(*grid, numAtomsMoved, nbs, nbat); @@ -1277,14 +1285,14 @@ static void calc_cell_indices(const nbnxn_search_t nbs, { fprintf(debug, "ns na_sc %d na_c %d super-cells: %d x %d y %d z %.1f maxz %d\n", grid->na_sc, grid->na_c, grid->nc, - grid->ncx, grid->ncy, grid->nc/((double)(grid->ncx*grid->ncy)), + grid->numCells[XX], grid->numCells[YY], grid->nc/((double)(grid->numCells[XX]*grid->numCells[YY])), ncz_max); if (gmx_debug_at) { int i = 0; - for (int cy = 0; cy < grid->ncy; cy++) + for (int cy = 0; cy < grid->numCells[YY]; cy++) { - for (int cx = 0; cx < grid->ncx; cx++) + for (int cx = 0; cx < grid->numCells[XX]; cx++) { fprintf(debug, " %2d", grid->cxy_ind[i+1]-grid->cxy_ind[i]); i++; @@ -1318,7 +1326,7 @@ static void calc_cell_indices(const nbnxn_search_t nbs, { /* Set the cell indices for the moved particles */ int n0 = grid->nc*grid->na_sc; - int n1 = grid->nc*grid->na_sc+grid->cxy_na[grid->ncx*grid->ncy]; + int n1 = grid->nc*grid->na_sc+grid->cxy_na[grid->numCells[XX]*grid->numCells[YY]]; if (ddZone == 0) { for (int i = n0; i < n1; i++) @@ -1334,18 +1342,18 @@ static void calc_cell_indices(const nbnxn_search_t nbs, { try { + int columnStart = ((thread + 0)*grid->numCells[XX]*grid->numCells[YY])/nthread; + int columnEnd = ((thread + 1)*grid->numCells[XX]*grid->numCells[YY])/nthread; if (grid->bSimple) { sort_columns_simple(nbs, ddZone, grid, atomStart, atomEnd, atinfo, x, nbat, - ((thread+0)*grid->ncx*grid->ncy)/nthread, - ((thread+1)*grid->ncx*grid->ncy)/nthread, + columnStart, columnEnd, nbs->work[thread].sortBuffer); } else { sort_columns_supersub(nbs, ddZone, grid, atomStart, atomEnd, atinfo, x, nbat, - ((thread+0)*grid->ncx*grid->ncy)/nthread, - ((thread+1)*grid->ncx*grid->ncy)/nthread, + columnStart, columnEnd, nbs->work[thread].sortBuffer); } } @@ -1512,8 +1520,8 @@ void nbnxn_put_on_grid_nonlocal(nbnxn_search_t nbs, void nbnxn_get_ncells(nbnxn_search_t nbs, int *ncx, int *ncy) { - *ncx = nbs->grid[0].ncx; - *ncy = nbs->grid[0].ncy; + *ncx = nbs->grid[0].numCells[XX]; + *ncy = nbs->grid[0].numCells[YY]; } void nbnxn_get_atomorder(const nbnxn_search_t nbs, const int **a, int *n) @@ -1525,7 +1533,7 @@ void nbnxn_get_atomorder(const nbnxn_search_t nbs, const int **a, int *n) /* Return the atom order for the home cell (index 0) */ *a = nbs->a.data(); - *n = grid->cxy_ind[grid->ncx*grid->ncy]*grid->na_sc; + *n = grid->cxy_ind[grid->numCells[XX]*grid->numCells[YY]]*grid->na_sc; } void nbnxn_set_atomorder(nbnxn_search_t nbs) @@ -1534,11 +1542,11 @@ void nbnxn_set_atomorder(nbnxn_search_t nbs) nbnxn_grid_t *grid = &nbs->grid[0]; int ao = 0; - for (int cx = 0; cx < grid->ncx; cx++) + for (int cx = 0; cx < grid->numCells[XX]; cx++) { - for (int cy = 0; cy < grid->ncy; cy++) + for (int cy = 0; cy < grid->numCells[YY]; cy++) { - int cxy = cx*grid->ncy + cy; + int cxy = cx*grid->numCells[YY] + cy; int j = grid->cxy_ind[cxy]*grid->na_sc; for (int cz = 0; cz < grid->cxy_na[cxy]; cz++) { diff --git a/src/gromacs/mdlib/nbnxn_internal.h b/src/gromacs/mdlib/nbnxn_internal.h index fe4a5a3448..f63d96e04c 100644 --- a/src/gromacs/mdlib/nbnxn_internal.h +++ b/src/gromacs/mdlib/nbnxn_internal.h @@ -153,27 +153,24 @@ typedef struct { /* A pair-search grid struct for one domain decomposition zone */ struct nbnxn_grid_t { - rvec c0; /* The lower corner of the (local) grid */ - rvec c1; /* The upper corner of the (local) grid */ - rvec size; /* c1 - c0 */ - real atom_density; /* The atom number density for the local grid */ - - gmx_bool bSimple; /* Is this grid simple or super/sub */ - int na_c; /* Number of atoms per cluster */ - int na_cj; /* Number of atoms for list j-clusters */ - int na_sc; /* Number of atoms per super-cluster */ - int na_c_2log; /* 2log of na_c */ - - int ncx; /* Number of (super-)cells along x */ - int ncy; /* Number of (super-)cells along y */ - int nc; /* Total number of (super-)cells */ - - real sx; /* x-size of a (super-)cell */ - real sy; /* y-size of a (super-)cell */ - real inv_sx; /* 1/sx */ - real inv_sy; /* 1/sy */ - - int cell0; /* Index in nbs->cell corresponding to cell 0 */ + rvec c0; /* The lower corner of the (local) grid */ + rvec c1; /* The upper corner of the (local) grid */ + rvec size; /* c1 - c0 */ + real atom_density; /* The atom number density for the local grid */ + + gmx_bool bSimple; /* Is this grid simple or super/sub */ + int na_c; /* Number of atoms per cluster */ + int na_cj; /* Number of atoms for list j-clusters */ + int na_sc; /* Number of atoms per super-cluster */ + int na_c_2log; /* 2log of na_c */ + + int numCells[DIM - 1]; /* Number of cells along x/y */ + int nc; /* Total number of cells */ + + real cellSize[DIM - 1]; /* size of a cell */ + real invCellSize[DIM - 1]; /* 1/cellSize */ + + int cell0; /* Index in nbs->cell corresponding to cell 0 */ /* Grid data */ std::vector cxy_na; /* The number of atoms for each column in x,y */ diff --git a/src/gromacs/mdlib/nbnxn_search.cpp b/src/gromacs/mdlib/nbnxn_search.cpp index 237ee4ebd7..8c4cefad68 100644 --- a/src/gromacs/mdlib/nbnxn_search.cpp +++ b/src/gromacs/mdlib/nbnxn_search.cpp @@ -361,19 +361,23 @@ static void init_buffer_flags(nbnxn_buffer_flags_t *flags, /* Determines the cell range along one dimension that * the bounding box b0 - b1 sees. */ +template static void get_cell_range(real b0, real b1, - int nc, real c0, real s, real invs, + const nbnxn_grid_t &gridj, real d2, real r2, int *cf, int *cl) { - *cf = std::max(static_cast((b0 - c0)*invs), 0); + real distanceInCells = (b0 - gridj.c0[dim])*gridj.invCellSize[dim]; + *cf = std::max(static_cast(distanceInCells), 0); - while (*cf > 0 && d2 + gmx::square((b0 - c0) - (*cf-1+1)*s) < r2) + while (*cf > 0 && + d2 + gmx::square((b0 - gridj.c0[dim]) - (*cf - 1 + 1)*gridj.cellSize[dim]) < r2) { (*cf)--; } - *cl = std::min(static_cast((b1 - c0)*invs), nc-1); - while (*cl < nc-1 && d2 + gmx::square((*cl+1)*s - (b1 - c0)) < r2) + *cl = std::min(static_cast((b1 - gridj.c0[dim])*gridj.invCellSize[dim]), gridj.numCells[dim] - 1); + while (*cl < gridj.numCells[dim] - 1 && + d2 + gmx::square((*cl + 1)*gridj.cellSize[dim] - (b1 - gridj.c0[dim])) < r2) { (*cl)++; } @@ -2570,12 +2574,12 @@ static real minimum_subgrid_size_xy(const nbnxn_grid_t *grid) { if (grid->bSimple) { - return std::min(grid->sx, grid->sy); + return std::min(grid->cellSize[XX], grid->cellSize[YY]); } else { - return std::min(grid->sx/c_gpuNumClusterPerCellX, - grid->sy/c_gpuNumClusterPerCellY); + return std::min(grid->cellSize[XX]/c_gpuNumClusterPerCellX, + grid->cellSize[YY]/c_gpuNumClusterPerCellY); } } @@ -2707,8 +2711,8 @@ static void get_nsubpair_target(const nbnxn_search_t nbs, return; } - ls[XX] = (grid->c1[XX] - grid->c0[XX])/(grid->ncx*c_gpuNumClusterPerCellX); - ls[YY] = (grid->c1[YY] - grid->c0[YY])/(grid->ncy*c_gpuNumClusterPerCellY); + ls[XX] = (grid->c1[XX] - grid->c0[XX])/(grid->numCells[XX]*c_gpuNumClusterPerCellX); + ls[YY] = (grid->c1[YY] - grid->c0[YY])/(grid->numCells[YY]*c_gpuNumClusterPerCellY); ls[ZZ] = grid->na_c/(grid->atom_density*ls[XX]*ls[YY]); /* The average length of the diagonal of a sub cell */ @@ -3075,10 +3079,10 @@ static gmx_bool next_ci(const nbnxn_grid_t *grid, return FALSE; } - while (*ci >= grid->cxy_ind[*ci_x*grid->ncy + *ci_y + 1]) + while (*ci >= grid->cxy_ind[*ci_x*grid->numCells[YY] + *ci_y + 1]) { *ci_y += 1; - if (*ci_y == grid->ncy) + if (*ci_y == grid->numCells[YY]) { *ci_x += 1; *ci_y = 0; @@ -3112,8 +3116,8 @@ static float boundingbox_only_distance2(const nbnxn_grid_t *gridi, real bbx, bby; real rbb2; - bbx = 0.5*(gridi->sx + gridj->sx); - bby = 0.5*(gridi->sy + gridj->sy); + bbx = 0.5*(gridi->cellSize[XX] + gridj->cellSize[XX]); + bby = 0.5*(gridi->cellSize[YY] + gridj->cellSize[YY]); if (!simple) { bbx /= c_gpuNumClusterPerCellX; @@ -3149,7 +3153,7 @@ static int get_ci_block_size(const nbnxn_grid_t *gridi, * zone boundaries with 3D domain decomposition. At the same time * the blocks will not become too small. */ - ci_block = (gridi->nc*ci_block_enum)/(ci_block_denom*gridi->ncx*nth); + ci_block = (gridi->nc*ci_block_enum)/(ci_block_denom*gridi->numCells[XX]*nth); /* Ensure the blocks are not too small: avoids cache invalidation */ if (ci_block*gridi->na_sc < ci_block_min_atoms) @@ -3329,7 +3333,7 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, if (debug) { fprintf(debug, "nbl nc_i %d col.av. %.1f ci_block %d\n", - gridi->nc, gridi->nc/(double)(gridi->ncx*gridi->ncy), ci_block); + gridi->nc, gridi->nc/(double)(gridi->numCells[XX]*gridi->numCells[YY]), ci_block); } numDistanceChecks = 0; @@ -3359,7 +3363,7 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, } else { - bx1 = gridi->c0[XX] + (ci_x+1)*gridi->sx; + bx1 = gridi->c0[XX] + (ci_x+1)*gridi->cellSize[XX]; } if (bx1 < gridj->c0[XX]) { @@ -3372,7 +3376,7 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, } } - ci_xy = ci_x*gridi->ncy + ci_y; + ci_xy = ci_x*gridi->numCells[YY] + ci_y; /* Loop over shift vectors in three dimensions */ for (int tz = -shp[ZZ]; tz <= shp[ZZ]; tz++) @@ -3420,14 +3424,14 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, } else { - by0 = gridi->c0[YY] + (ci_y )*gridi->sy + shy; - by1 = gridi->c0[YY] + (ci_y+1)*gridi->sy + shy; + by0 = gridi->c0[YY] + (ci_y )*gridi->cellSize[YY] + shy; + by1 = gridi->c0[YY] + (ci_y+1)*gridi->cellSize[YY] + shy; } - get_cell_range(by0, by1, - gridj->ncy, gridj->c0[YY], gridj->sy, gridj->inv_sy, - d2z_cx, rlist2, - &cyf, &cyl); + get_cell_range(by0, by1, + *gridj, + d2z_cx, rlist2, + &cyf, &cyl); if (cyf > cyl) { @@ -3462,14 +3466,14 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, } else { - bx0 = gridi->c0[XX] + (ci_x )*gridi->sx + shx; - bx1 = gridi->c0[XX] + (ci_x+1)*gridi->sx + shx; + bx0 = gridi->c0[XX] + (ci_x )*gridi->cellSize[XX] + shx; + bx1 = gridi->c0[XX] + (ci_x+1)*gridi->cellSize[XX] + shx; } - get_cell_range(bx0, bx1, - gridj->ncx, gridj->c0[XX], gridj->sx, gridj->inv_sx, - d2z_cy, rlist2, - &cxf, &cxl); + get_cell_range(bx0, bx1, + *gridj, + d2z_cy, rlist2, + &cxf, &cxl); if (cxf > cxl) { @@ -3518,13 +3522,13 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, for (int cx = cxf; cx <= cxl; cx++) { d2zx = d2z; - if (gridj->c0[XX] + cx*gridj->sx > bx1) + if (gridj->c0[XX] + cx*gridj->cellSize[XX] > bx1) { - d2zx += gmx::square(gridj->c0[XX] + cx*gridj->sx - bx1); + d2zx += gmx::square(gridj->c0[XX] + cx*gridj->cellSize[XX] - bx1); } - else if (gridj->c0[XX] + (cx+1)*gridj->sx < bx0) + else if (gridj->c0[XX] + (cx+1)*gridj->cellSize[XX] < bx0) { - d2zx += gmx::square(gridj->c0[XX] + (cx+1)*gridj->sx - bx0); + d2zx += gmx::square(gridj->c0[XX] + (cx+1)*gridj->cellSize[XX] - bx0); } if (gridi == gridj && @@ -3544,17 +3548,17 @@ static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs, for (int cy = cyf_x; cy <= cyl; cy++) { - const int columnStart = gridj->cxy_ind[cx*gridj->ncy + cy]; - const int columnEnd = gridj->cxy_ind[cx*gridj->ncy + cy + 1]; + const int columnStart = gridj->cxy_ind[cx*gridj->numCells[YY] + cy]; + const int columnEnd = gridj->cxy_ind[cx*gridj->numCells[YY] + cy + 1]; d2zxy = d2zx; - if (gridj->c0[YY] + cy*gridj->sy > by1) + if (gridj->c0[YY] + cy*gridj->cellSize[YY] > by1) { - d2zxy += gmx::square(gridj->c0[YY] + cy*gridj->sy - by1); + d2zxy += gmx::square(gridj->c0[YY] + cy*gridj->cellSize[YY] - by1); } - else if (gridj->c0[YY] + (cy+1)*gridj->sy < by0) + else if (gridj->c0[YY] + (cy+1)*gridj->cellSize[YY] < by0) { - d2zxy += gmx::square(gridj->c0[YY] + (cy+1)*gridj->sy - by0); + d2zxy += gmx::square(gridj->c0[YY] + (cy+1)*gridj->cellSize[YY] - by0); } if (columnStart < columnEnd && d2zxy < rlist2) { -- 2.11.4.GIT