src/gromacs/domdec/redistribute.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, 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 atom redistribution functions.
  38  *
  39  * \author Berk Hess <hess@kth.se>
  40  * \ingroup module_domdec
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "redistribute.h"
  46
  47 #include <cstring>
  48
  49 #include "gromacs/domdec/domdec_network.h"
  50 #include "gromacs/domdec/ga2la.h"
  51 #include "gromacs/gmxlib/nrnb.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  54 #include "gromacs/mdlib/nsgrid.h"
  55 #include "gromacs/mdtypes/forcerec.h"
  56 #include "gromacs/mdtypes/nblist.h"
  57 #include "gromacs/mdtypes/state.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/gmxomp.h"
  61
  62 #include "domdec_internal.h"
  63 #include "utility.h"
  64
  65 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
  66 static constexpr int DD_CGIBS = 2;
  67
  68 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
  69
  70 /* The lower 16 bits are reserved for the charge group size */
  71 static constexpr int DD_FLAG_NRCG = 65535;
  72
  73 /* Returns which bit tells whether to move a group forward along dimension d */
  74 static inline int DD_FLAG_FW(int d)
  75 {
  76     return 1 << (16 + d*2);
  77 }
  78
  79 /* Returns which bit tells whether to move a group backward along dimension d */
  80 static inline int DD_FLAG_BW(int d)
  81 {
  82     return 1 << (16 + d*2 + 1);
  83 }
  84
  85 static void
  86 copyMovedAtomsToBufferPerAtom(gmx::ArrayRef<const int> move,
  87                               const gmx::RangePartitioning &atomGroups,
  88                               int nvec, int vec,
  89                               rvec *src, gmx_domdec_comm_t *comm)
  90 {
  91     int pos_vec[DIM*2] = { 0 };
  92
  93     for (int g = 0; g < move.size(); g++)
  94     {
  95         const auto atomGroup = atomGroups.block(g);
  96         /* Skip moved atoms */
  97         int        m         = move[g];
  98         if (m >= 0)
  99         {
 100             /* Copy to the communication buffer */
 101             int numAtoms  = atomGroup.size();
 102             pos_vec[m]   += 1 + vec*numAtoms;
 103             for (int i : atomGroup)
 104             {
 105                 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
 106             }
 107             pos_vec[m] += (nvec - vec - 1)*numAtoms;
 108         }
 109     }
 110 }
 111
 112 static void
 113 copyMovedAtomsToBufferPerChargeGroup(gmx::ArrayRef<const int> move,
 114                                      const gmx::RangePartitioning &atomGroups,
 115                                      int nvec, rvec *src, gmx_domdec_comm_t *comm)
 116 {
 117     int pos_vec[DIM*2] = { 0 };
 118
 119     for (int g = 0; g < move.size(); g++)
 120     {
 121         const auto atomGroup = atomGroups.block(g);
 122         /* Skip moved atoms */
 123         int        m         = move[g];
 124         if (m >= 0)
 125         {
 126             /* Copy to the communication buffer */
 127             copy_rvec(src[g], comm->cgcm_state[m][pos_vec[m]]);
 128             pos_vec[m] += 1 + atomGroup.size()*nvec;
 129         }
 130     }
 131 }
 132
 133 static void clear_and_mark_ind(gmx::ArrayRef<const int>      move,
 134                                gmx::ArrayRef<const int>      globalAtomGroupIndices,
 135                                const gmx::RangePartitioning &atomGroups,
 136                                gmx::ArrayRef<const int>      globalAtomIndices,
 137                                gmx_ga2la_t                  *ga2la,
 138                                char                         *bLocalCG,
 139                                int                          *cell_index)
 140 {
 141     for (int g = 0; g < move.size(); g++)
 142     {
 143         if (move[g] >= 0)
 144         {
 145             /* Clear the global indices */
 146             for (int a : atomGroups.block(g))
 147             {
 148                 ga2la->erase(globalAtomIndices[a]);
 149             }
 150             if (bLocalCG)
 151             {
 152                 bLocalCG[globalAtomGroupIndices[g]] = FALSE;
 153             }
 154             /* Signal that this group has moved using the ns cell index.
 155              * Here we set it to -1. fill_grid will change it
 156              * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
 157              */
 158             cell_index[g] = -1;
 159         }
 160     }
 161 }
 162
 163 static void print_cg_move(FILE *fplog,
 164                           gmx_domdec_t *dd,
 165                           int64_t step, int cg, int dim, int dir,
 166                           gmx_bool bHaveCgcmOld, real limitd,
 167                           rvec cm_old, rvec cm_new, real pos_d)
 168 {
 169     gmx_domdec_comm_t *comm;
 170     char               buf[22];
 171
 172     comm = dd->comm;
 173
 174     fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
 175     if (limitd > 0)
 176     {
 177         fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
 178                 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
 179                 ddglatnr(dd, dd->atomGrouping().block(cg).begin()), limitd, dim2char(dim));
 180     }
 181     else
 182     {
 183         /* We don't have a limiting distance available: don't print it */
 184         fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
 185                 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
 186                 ddglatnr(dd, dd->atomGrouping().block(cg).begin()), dim2char(dim));
 187     }
 188     fprintf(fplog, "distance out of cell %f\n",
 189             dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
 190     if (bHaveCgcmOld)
 191     {
 192         fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
 193                 cm_old[XX], cm_old[YY], cm_old[ZZ]);
 194     }
 195     fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
 196             cm_new[XX], cm_new[YY], cm_new[ZZ]);
 197     fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
 198             dim2char(dim),
 199             comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
 200     fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
 201             dim2char(dim),
 202             comm->cell_x0[dim], comm->cell_x1[dim]);
 203 }
 204
 205 [[ noreturn ]]
 206 static void cg_move_error(FILE *fplog,
 207                           gmx_domdec_t *dd,
 208                           int64_t step, int cg, int dim, int dir,
 209                           gmx_bool bHaveCgcmOld, real limitd,
 210                           rvec cm_old, rvec cm_new, real pos_d)
 211 {
 212     if (fplog)
 213     {
 214         print_cg_move(fplog, dd, step, cg, dim, dir,
 215                       bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
 216     }
 217     print_cg_move(stderr, dd, step, cg, dim, dir,
 218                   bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
 219     gmx_fatal(FARGS,
 220               "%s moved too far between two domain decomposition steps\n"
 221               "This usually means that your system is not well equilibrated",
 222               dd->comm->bCGs ? "An atom group" : "An atom");
 223 }
 224
 225 static void rotate_state_atom(t_state *state, int a)
 226 {
 227     if (state->flags & (1 << estX))
 228     {
 229         auto x = makeArrayRef(state->x);
 230         /* Rotate the complete state; for a rectangular box only */
 231         x[a][YY] = state->box[YY][YY] - x[a][YY];
 232         x[a][ZZ] = state->box[ZZ][ZZ] - x[a][ZZ];
 233     }
 234     if (state->flags & (1 << estV))
 235     {
 236         auto v = makeArrayRef(state->v);
 237         v[a][YY] = -v[a][YY];
 238         v[a][ZZ] = -v[a][ZZ];
 239     }
 240     if (state->flags & (1 << estCGP))
 241     {
 242         auto cg_p = makeArrayRef(state->cg_p);
 243         cg_p[a][YY] = -cg_p[a][YY];
 244         cg_p[a][ZZ] = -cg_p[a][ZZ];
 245     }
 246 }
 247
 248 /* Returns a pointer to a buffer of size \p numAtomsNew with contents 0 from numAtomsOld upwards
 249  *
 250  * Note: numAtomsOld should either be 0 or match the current buffer size.
 251  */
 252 static int *getMovedBuffer(gmx_domdec_comm_t *comm,
 253                            size_t             numAtomsOld,
 254                            size_t             numAtomsNew)
 255 {
 256     std::vector<int> &movedBuffer = comm->movedBuffer;
 257
 258     GMX_RELEASE_ASSERT(numAtomsOld == 0 || movedBuffer.size() == numAtomsOld, "numAtomsOld should either be 0 or match the current size");
 259
 260     /* Contents up to numAtomsOld should be preserved, clear from numAtomsOld */
 261     if (numAtomsOld == 0)
 262     {
 263         movedBuffer.clear();
 264     }
 265     movedBuffer.resize(numAtomsNew);
 266
 267     return movedBuffer.data();
 268 }
 269
 270 /* Bounds to determine whether an atom group moved to far between DD steps */
 271 struct MoveLimits
 272 {
 273     rvec distance; /* Maximum distance over which a group can move */
 274     rvec lower;    /* Lower bound for group location */
 275     rvec upper;    /* Upper bound for group location */
 276 };
 277
 278 /* Compute flag that tells where to move an atom group
 279  *
 280  * The input dev tells whether the group moves fw,remains,bw (-1,0,1) along
 281  * dimensions.
 282  * The return value is -1 when the group does not move.
 283  * The return has move flags set when the group does move and the lower 4 bits
 284  * are (mis)used to tell which is the first dimension (bit 1,2,3) the group
 285  * needs to be moved along and in which direction (bit 0 not set for fw
 286  * and set for bw).
 287  */
 288 static int computeMoveFlag(const gmx_domdec_t &dd,
 289                            const ivec         &dev)
 290 {
 291     int flag              = 0;
 292     int firstMoveDimValue = -1;
 293     for (int d = 0; d < dd.ndim; d++)
 294     {
 295         const int dim = dd.dim[d];
 296         if (dev[dim] == 1)
 297         {
 298             flag |= DD_FLAG_FW(d);
 299             if (firstMoveDimValue == -1)
 300             {
 301                 firstMoveDimValue = d*2;
 302             }
 303         }
 304         else if (dev[dim] == -1)
 305         {
 306             flag |= DD_FLAG_BW(d);
 307             if (firstMoveDimValue == -1)
 308             {
 309                 if (dd.nc[dim] > 2)
 310                 {
 311                     firstMoveDimValue = d*2 + 1;
 312                 }
 313                 else
 314                 {
 315                     firstMoveDimValue = d*2;
 316                 }
 317             }
 318         }
 319     }
 320
 321     return firstMoveDimValue + flag;
 322 }
 323
 324 /* Determine to which domains atomGroups in the range \p cg_start, \p cg_end should go.
 325  *
 326  * Returns in the move array where the groups should go.
 327  * Also updates \p cg_cm for jumps over periodic boundaries.
 328  *
 329  * \TODO Rename cg to atomGroup.
 330  */
 331 static void calc_cg_move(FILE *fplog, int64_t step,
 332                          gmx_domdec_t *dd,
 333                          t_state *state,
 334                          const ivec tric_dir, matrix tcm,
 335                          const rvec cell_x0, const rvec cell_x1,
 336                          const MoveLimits &moveLimits,
 337                          const gmx::RangePartitioning &atomGroups,
 338                          int cg_start, int cg_end,
 339                          rvec *cg_cm,
 340                          gmx::ArrayRef<int> move)
 341 {
 342     const int npbcdim = dd->npbcdim;
 343     auto      x       = makeArrayRef(state->x);
 344
 345     for (int g = cg_start; g < cg_end; g++)
 346     {
 347         const auto atomGroup = atomGroups.block(g);
 348         const int  numAtoms  = atomGroup.size();
 349         // TODO: Rename this center of geometry variable to cogNew
 350         rvec       cm_new;
 351         if (numAtoms == 1)
 352         {
 353             copy_rvec(x[atomGroup.begin()], cm_new);
 354         }
 355         else
 356         {
 357             real invNumAtoms = 1/static_cast<real>(numAtoms);
 358             clear_rvec(cm_new);
 359             for (int k : atomGroup)
 360             {
 361                 rvec_inc(cm_new, x[k]);
 362             }
 363             for (int d = 0; d < DIM; d++)
 364             {
 365                 cm_new[d] = invNumAtoms*cm_new[d];
 366             }
 367         }
 368
 369         ivec dev = { 0 };
 370         /* Do pbc and check DD cell boundary crossings */
 371         for (int d = DIM - 1; d >= 0; d--)
 372         {
 373             if (dd->nc[d] > 1)
 374             {
 375                 bool bScrew = (dd->bScrewPBC && d == XX);
 376                 /* Determine the location of this cg in lattice coordinates */
 377                 real pos_d = cm_new[d];
 378                 if (tric_dir[d])
 379                 {
 380                     for (int d2 = d + 1; d2 < DIM; d2++)
 381                     {
 382                         pos_d += cm_new[d2]*tcm[d2][d];
 383                     }
 384                 }
 385                 /* Put the charge group in the triclinic unit-cell */
 386                 if (pos_d >= cell_x1[d])
 387                 {
 388                     if (pos_d >= moveLimits.upper[d])
 389                     {
 390                         cg_move_error(fplog, dd, step, g, d, 1,
 391                                       cg_cm != state->x.rvec_array(), moveLimits.distance[d],
 392                                       cg_cm[g], cm_new, pos_d);
 393                     }
 394                     dev[d] = 1;
 395                     if (dd->ci[d] == dd->nc[d] - 1)
 396                     {
 397                         rvec_dec(cm_new, state->box[d]);
 398                         if (bScrew)
 399                         {
 400                             cm_new[YY] = state->box[YY][YY] - cm_new[YY];
 401                             cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
 402                         }
 403                         for (int k : atomGroup)
 404                         {
 405                             rvec_dec(x[k], state->box[d]);
 406                             if (bScrew)
 407                             {
 408                                 rotate_state_atom(state, k);
 409                             }
 410                         }
 411                     }
 412                 }
 413                 else if (pos_d < cell_x0[d])
 414                 {
 415                     if (pos_d < moveLimits.lower[d])
 416                     {
 417                         cg_move_error(fplog, dd, step, g, d, -1,
 418                                       cg_cm != state->x.rvec_array(), moveLimits.distance[d],
 419                                       cg_cm[g], cm_new, pos_d);
 420                     }
 421                     dev[d] = -1;
 422                     if (dd->ci[d] == 0)
 423                     {
 424                         rvec_inc(cm_new, state->box[d]);
 425                         if (bScrew)
 426                         {
 427                             cm_new[YY] = state->box[YY][YY] - cm_new[YY];
 428                             cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
 429                         }
 430                         for (int k : atomGroup)
 431                         {
 432                             rvec_inc(x[k], state->box[d]);
 433                             if (bScrew)
 434                             {
 435                                 rotate_state_atom(state, k);
 436                             }
 437                         }
 438                     }
 439                 }
 440             }
 441             else if (d < npbcdim)
 442             {
 443                 /* Put the charge group in the rectangular unit-cell */
 444                 while (cm_new[d] >= state->box[d][d])
 445                 {
 446                     rvec_dec(cm_new, state->box[d]);
 447                     for (int k : atomGroup)
 448                     {
 449                         rvec_dec(x[k], state->box[d]);
 450                     }
 451                 }
 452                 while (cm_new[d] < 0)
 453                 {
 454                     rvec_inc(cm_new, state->box[d]);
 455                     for (int k : atomGroup)
 456                     {
 457                         rvec_inc(x[k], state->box[d]);
 458                     }
 459                 }
 460             }
 461         }
 462
 463         copy_rvec(cm_new, cg_cm[g]);
 464
 465         /* Temporarily store the flag in move */
 466         move[g] = computeMoveFlag(*dd, dev);
 467     }
 468 }
 469
 470 struct PbcAndFlag
 471 {
 472     /* Constructor that purposely does not initialize anything */
 473     PbcAndFlag()
 474     {
 475     }
 476
 477     gmx::RVec pbcShift;
 478     int       moveFlag;
 479 };
 480
 481 /* Determine to which domains update groups in the range \p groupBegin, \p groupEnd should go.
 482  *
 483  * Returns in the move array where the groups should go.
 484  * Also updates the COGs and coordinates for jumps over periodic boundaries.
 485  */
 486 static void calcGroupMove(FILE *fplog, int64_t step,
 487                           gmx_domdec_t *dd,
 488                           t_state *state,
 489                           const ivec tric_dir, matrix tcm,
 490                           const rvec cell_x0, const rvec cell_x1,
 491                           const MoveLimits &moveLimits,
 492                           int groupBegin, int groupEnd,
 493                           gmx::ArrayRef<PbcAndFlag> pbcAndFlags)
 494 {
 495     GMX_RELEASE_ASSERT(!dd->bScrewPBC, "Screw PBC is not supported here");
 496
 497     const int             npbcdim         = dd->npbcdim;
 498
 499     gmx::UpdateGroupsCog *updateGroupsCog = dd->comm->updateGroupsCog.get();
 500
 501     for (int g = groupBegin; g < groupEnd; g++)
 502     {
 503
 504         gmx::RVec &cog    = updateGroupsCog->cog(g);
 505         gmx::RVec  cogOld = cog;
 506
 507         ivec       dev = { 0 };
 508         /* Do pbc and check DD cell boundary crossings */
 509         for (int d = DIM - 1; d >= 0; d--)
 510         {
 511             if (dd->nc[d] > 1)
 512             {
 513                 /* Determine the location of this COG in lattice coordinates */
 514                 real pos_d = cog[d];
 515                 if (tric_dir[d])
 516                 {
 517                     for (int d2 = d + 1; d2 < DIM; d2++)
 518                     {
 519                         pos_d += cog[d2]*tcm[d2][d];
 520                     }
 521                 }
 522                 /* Put the COG in the triclinic unit-cell */
 523                 if (pos_d >= cell_x1[d])
 524                 {
 525                     if (pos_d >= moveLimits.upper[d])
 526                     {
 527                         cg_move_error(fplog, dd, step, g, d, 1,
 528                                       true, moveLimits.distance[d],
 529                                       cogOld, cog, pos_d);
 530                     }
 531                     dev[d] = 1;
 532                     if (dd->ci[d] == dd->nc[d] - 1)
 533                     {
 534                         rvec_dec(cog, state->box[d]);
 535                     }
 536                 }
 537                 else if (pos_d < cell_x0[d])
 538                 {
 539                     if (pos_d < moveLimits.lower[d])
 540                     {
 541                         cg_move_error(fplog, dd, step, g, d, -1,
 542                                       true, moveLimits.distance[d],
 543                                       cogOld, cog, pos_d);
 544                     }
 545                     dev[d] = -1;
 546                     if (dd->ci[d] == 0)
 547                     {
 548                         rvec_inc(cog, state->box[d]);
 549                     }
 550                 }
 551             }
 552             else if (d < npbcdim)
 553             {
 554                 /* Put the COG in the rectangular unit-cell */
 555                 while (cog[d] >= state->box[d][d])
 556                 {
 557                     rvec_dec(cog, state->box[d]);
 558                 }
 559                 while (cog[d] < 0)
 560                 {
 561                     rvec_inc(cog, state->box[d]);
 562                 }
 563             }
 564         }
 565
 566         /* Store the PBC and move flag, so we can later apply them to the atoms */
 567         PbcAndFlag &pbcAndFlag = pbcAndFlags[g];
 568
 569         rvec_sub(cog, cogOld, pbcAndFlag.pbcShift);
 570         pbcAndFlag.moveFlag = computeMoveFlag(*dd, dev);
 571     }
 572 }
 573
 574 static void
 575 applyPbcAndSetMoveFlags(const gmx::UpdateGroupsCog      &updateGroupsCog,
 576                         gmx::ArrayRef<const PbcAndFlag>  pbcAndFlags,
 577                         int                              atomBegin,
 578                         int                              atomEnd,
 579                         gmx::ArrayRef<gmx::RVec>         atomCoords,
 580                         gmx::ArrayRef<int>               move)
 581 {
 582     for (int a = atomBegin; a < atomEnd; a++)
 583     {
 584         const PbcAndFlag &pbcAndFlag = pbcAndFlags[updateGroupsCog.cogIndex(a)];
 585         rvec_inc(atomCoords[a], pbcAndFlag.pbcShift);
 586         /* Temporarily store the flag in move */
 587         move[a] = pbcAndFlag.moveFlag;
 588     }
 589 }
 590
 591 void dd_redistribute_cg(FILE *fplog, int64_t step,
 592                         gmx_domdec_t *dd, ivec tric_dir,
 593                         t_state *state,
 594                         PaddedVector<gmx::RVec> *f,
 595                         t_forcerec *fr,
 596                         t_nrnb *nrnb,
 597                         int *ncg_moved)
 598 {
 599     gmx_domdec_comm_t *comm = dd->comm;
 600
 601     if (dd->bScrewPBC)
 602     {
 603         check_screw_box(state->box);
 604     }
 605
 606     rvec *cg_cm = nullptr;
 607     if (fr->cutoff_scheme == ecutsGROUP)
 608     {
 609         cg_cm = fr->cg_cm;
 610     }
 611
 612     // Positions are always present, so there's nothing to flag
 613     bool                bV   = (state->flags & (1<<estV)) != 0;
 614     bool                bCGP = (state->flags & (1<<estCGP)) != 0;
 615
 616     DDBufferAccess<int> moveBuffer(comm->intBuffer, dd->ncg_home);
 617     gmx::ArrayRef<int>  move = moveBuffer.buffer;
 618
 619     const int           npbcdim = dd->npbcdim;
 620
 621     rvec                cell_x0, cell_x1;
 622     MoveLimits          moveLimits;
 623     for (int d = 0; (d < DIM); d++)
 624     {
 625         moveLimits.distance[d] = dd->comm->cellsize_min[d];
 626         if (d >= npbcdim && dd->ci[d] == 0)
 627         {
 628             cell_x0[d] = -GMX_FLOAT_MAX;
 629         }
 630         else
 631         {
 632             cell_x0[d] = comm->cell_x0[d];
 633         }
 634         if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
 635         {
 636             cell_x1[d] = GMX_FLOAT_MAX;
 637         }
 638         else
 639         {
 640             cell_x1[d] = comm->cell_x1[d];
 641         }
 642         if (d < npbcdim)
 643         {
 644             moveLimits.lower[d] = comm->old_cell_x0[d] - moveLimits.distance[d];
 645             moveLimits.upper[d] = comm->old_cell_x1[d] + moveLimits.distance[d];
 646         }
 647         else
 648         {
 649             /* We check after communication if a charge group moved
 650              * more than one cell. Set the pre-comm check limit to float_max.
 651              */
 652             moveLimits.lower[d] = -GMX_FLOAT_MAX;
 653             moveLimits.upper[d] =  GMX_FLOAT_MAX;
 654         }
 655     }
 656
 657     matrix tcm;
 658     make_tric_corr_matrix(npbcdim, state->box, tcm);
 659
 660     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
 661
 662     const int                     nthread = gmx_omp_nthreads_get(emntDomdec);
 663
 664     /* Compute the center of geometry for all home charge groups
 665      * and put them in the box and determine where they should go.
 666      */
 667     std::vector<PbcAndFlag>  pbcAndFlags(comm->useUpdateGroups ? comm->updateGroupsCog->numCogs() : 0);
 668
 669 #pragma omp parallel num_threads(nthread)
 670     {
 671         try
 672         {
 673             const int thread = gmx_omp_get_thread_num();
 674
 675             if (comm->useUpdateGroups)
 676             {
 677                 const auto &updateGroupsCog = *comm->updateGroupsCog;
 678                 const int   numGroups       = updateGroupsCog.numCogs();
 679                 calcGroupMove(fplog, step, dd, state, tric_dir, tcm,
 680                               cell_x0, cell_x1, moveLimits,
 681                               ( thread   *numGroups)/nthread,
 682                               ((thread+1)*numGroups)/nthread,
 683                               pbcAndFlags);
 684                 /* We need a barrier as atoms below can be in a COG of a different thread */
 685 #pragma omp barrier
 686                 const int numHomeAtoms = comm->atomRanges.numHomeAtoms();
 687                 applyPbcAndSetMoveFlags(updateGroupsCog, pbcAndFlags,
 688                                         ( thread   *numHomeAtoms)/nthread,
 689                                         ((thread+1)*numHomeAtoms)/nthread,
 690                                         state->x,
 691                                         move);
 692             }
 693             else
 694             {
 695                 /* Here we handle single atoms or charge groups */
 696                 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
 697                              cell_x0, cell_x1, moveLimits,
 698                              atomGrouping,
 699                              ( thread   *dd->ncg_home)/nthread,
 700                              ((thread+1)*dd->ncg_home)/nthread,
 701                              fr->cutoff_scheme == ecutsGROUP ? cg_cm : as_rvec_array(state->x.data()),
 702                              move);
 703             }
 704         }
 705         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 706     }
 707
 708     int ncg[DIM*2] = { 0 };
 709     int nat[DIM*2] = { 0 };
 710     for (int cg = 0; cg < dd->ncg_home; cg++)
 711     {
 712         if (move[cg] >= 0)
 713         {
 714             const int flag = move[cg] & ~DD_FLAG_NRCG;
 715             const int mc   = move[cg] & DD_FLAG_NRCG;
 716             move[cg]       = mc;
 717
 718             std::vector<int> &cggl_flag = comm->cggl_flag[mc];
 719
 720             /* TODO: See if we can use push_back instead */
 721             if ((ncg[mc] + 1)*DD_CGIBS > gmx::index(cggl_flag.size()))
 722             {
 723                 cggl_flag.resize((ncg[mc] + 1)*DD_CGIBS);
 724             }
 725             cggl_flag[ncg[mc]*DD_CGIBS  ] = dd->globalAtomGroupIndices[cg];
 726             /* We store the cg size in the lower 16 bits
 727              * and the place where the charge group should go
 728              * in the next 6 bits. This saves some communication volume.
 729              */
 730             const int nrcg = atomGrouping.block(cg).size();
 731             cggl_flag[ncg[mc]*DD_CGIBS+1] = nrcg | flag;
 732             ncg[mc] += 1;
 733             nat[mc] += nrcg;
 734         }
 735     }
 736
 737     inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
 738     inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
 739
 740     *ncg_moved = 0;
 741     for (int i = 0; i < dd->ndim*2; i++)
 742     {
 743         *ncg_moved += ncg[i];
 744     }
 745
 746     int nvec = 1;
 747     if (bV)
 748     {
 749         nvec++;
 750     }
 751     if (bCGP)
 752     {
 753         nvec++;
 754     }
 755
 756     /* Make sure the communication buffers are large enough */
 757     for (int mc = 0; mc < dd->ndim*2; mc++)
 758     {
 759         size_t nvr = ncg[mc] + nat[mc]*nvec;
 760         if (nvr > comm->cgcm_state[mc].size())
 761         {
 762             comm->cgcm_state[mc].resize(nvr);
 763         }
 764     }
 765
 766     switch (fr->cutoff_scheme)
 767     {
 768         case ecutsGROUP:
 769             /* Recalculating cg_cm might be cheaper than communicating,
 770              * but that could give rise to rounding issues.
 771              */
 772             copyMovedAtomsToBufferPerChargeGroup(move, dd->atomGrouping(),
 773                                                  nvec, cg_cm, comm);
 774             break;
 775         case ecutsVERLET:
 776             /* Without charge groups we send the moved atom coordinates
 777              * over twice. This is so the code below can be used without
 778              * many conditionals for both for with and without charge groups.
 779              */
 780             copyMovedAtomsToBufferPerChargeGroup(move, dd->atomGrouping(),
 781                                                  nvec, state->x.rvec_array(), comm);
 782             break;
 783         default:
 784             gmx_incons("unimplemented");
 785     }
 786
 787     int vectorIndex = 0;
 788     copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
 789                                   nvec, vectorIndex++,
 790                                   state->x.rvec_array(),
 791                                   comm);
 792     if (bV)
 793     {
 794         copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
 795                                       nvec, vectorIndex++,
 796                                       state->v.rvec_array(),
 797                                       comm);
 798     }
 799     if (bCGP)
 800     {
 801         copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
 802                                       nvec, vectorIndex++,
 803                                       state->cg_p.rvec_array(),
 804                                       comm);
 805     }
 806
 807     int *moved;
 808     if (fr->cutoff_scheme == ecutsVERLET)
 809     {
 810         moved = getMovedBuffer(comm, 0, dd->ncg_home);
 811     }
 812     else
 813     {
 814         moved = fr->ns->grid->cell_index;
 815     }
 816
 817     clear_and_mark_ind(move,
 818                        dd->globalAtomGroupIndices, dd->atomGrouping(), dd->globalAtomIndices,
 819                        dd->ga2la, comm->bLocalCG,
 820                        moved);
 821
 822     /* Now we can remove the excess global atom-group indices from the list */
 823     dd->globalAtomGroupIndices.resize(dd->ncg_home);
 824     dd->atomGrouping_.reduceNumBlocks(dd->ncg_home);
 825
 826     /* We reuse the intBuffer without reacquiring since we are in the same scope */
 827     DDBufferAccess<int> &flagBuffer = moveBuffer;
 828
 829     const cginfo_mb_t   *cginfo_mb  = fr->cginfo_mb;
 830
 831     /* Temporarily store atoms passed to our rank at the end of the range */
 832     int home_pos_cg = dd->ncg_home;
 833     int home_pos_at = dd->atomGrouping().subRange(0, dd->ncg_home).end();
 834     for (int d = 0; d < dd->ndim; d++)
 835     {
 836         DDBufferAccess<gmx::RVec> rvecBuffer(comm->rvecBuffer, 0);
 837
 838         const int                 dim = dd->dim[d];
 839         int ncg_recv                  = 0;
 840         int                       nvr = 0;
 841         for (int dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
 842         {
 843             const int cdd = d*2 + dir;
 844             /* Communicate the cg and atom counts */
 845             int       sbuf[2] = { ncg[cdd], nat[cdd] };
 846             if (debug)
 847             {
 848                 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
 849                         d, dir, sbuf[0], sbuf[1]);
 850             }
 851             int rbuf[2];
 852             ddSendrecv(dd, d, dir, sbuf, 2, rbuf, 2);
 853
 854             flagBuffer.resize((ncg_recv + rbuf[0])*DD_CGIBS);
 855
 856             /* Communicate the charge group indices, sizes and flags */
 857             ddSendrecv(dd, d, dir,
 858                        comm->cggl_flag[cdd].data(), sbuf[0]*DD_CGIBS,
 859                        flagBuffer.buffer.data() + ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
 860
 861             const int nvs = ncg[cdd] + nat[cdd]*nvec;
 862             const int i   = rbuf[0]  + rbuf[1] *nvec;
 863             rvecBuffer.resize(nvr + i);
 864
 865             /* Communicate cgcm and state */
 866             ddSendrecv(dd, d, dir,
 867                        as_rvec_array(comm->cgcm_state[cdd].data()), nvs,
 868                        as_rvec_array(rvecBuffer.buffer.data()) + nvr, i);
 869             ncg_recv += rbuf[0];
 870             nvr      += i;
 871         }
 872
 873         dd_check_alloc_ncg(fr, state, f, home_pos_cg + ncg_recv);
 874         if (fr->cutoff_scheme == ecutsGROUP)
 875         {
 876             /* Here we resize to more than necessary and shrink later */
 877             dd_resize_state(state, f, home_pos_at + ncg_recv*MAX_CGCGSIZE);
 878         }
 879
 880         /* Process the received charge groups */
 881         int buf_pos = 0;
 882         for (int cg = 0; cg < ncg_recv; cg++)
 883         {
 884             int   flag = flagBuffer.buffer[cg*DD_CGIBS + 1];
 885             rvec &pos  = as_rvec_array(rvecBuffer.buffer.data())[buf_pos];
 886
 887             if (dim >= npbcdim && dd->nc[dim] > 2)
 888             {
 889                 rvec &pos = as_rvec_array(rvecBuffer.buffer.data())[buf_pos];
 890
 891                 /* No pbc in this dim and more than one domain boundary.
 892                  * We do a separate check if a charge group didn't move too far.
 893                  */
 894                 if (((flag & DD_FLAG_FW(d)) &&
 895                      pos[dim] > cell_x1[dim]) ||
 896                     ((flag & DD_FLAG_BW(d)) &&
 897                      pos[dim] < cell_x0[dim]))
 898                 {
 899                     cg_move_error(fplog, dd, step, cg, dim,
 900                                   (flag & DD_FLAG_FW(d)) ? 1 : 0,
 901                                   fr->cutoff_scheme == ecutsGROUP, 0,
 902                                   pos, pos, pos[dim]);
 903                 }
 904             }
 905
 906             int mc = -1;
 907             if (d < dd->ndim-1)
 908             {
 909                 /* Check which direction this cg should go */
 910                 for (int d2 = d + 1; (d2 < dd->ndim && mc == -1); d2++)
 911                 {
 912                     if (isDlbOn(dd->comm))
 913                     {
 914                         /* The cell boundaries for dimension d2 are not equal
 915                          * for each cell row of the lower dimension(s),
 916                          * therefore we might need to redetermine where
 917                          * this cg should go.
 918                          */
 919                         const int dim2 = dd->dim[d2];
 920                         /* If this cg crosses the box boundary in dimension d2
 921                          * we can use the communicated flag, so we do not
 922                          * have to worry about pbc.
 923                          */
 924                         if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
 925                                (flag & DD_FLAG_FW(d2))) ||
 926                               (dd->ci[dim2] == 0 &&
 927                                (flag & DD_FLAG_BW(d2)))))
 928                         {
 929                             /* Clear the two flags for this dimension */
 930                             flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
 931                             /* Determine the location of this cg
 932                              * in lattice coordinates
 933                              */
 934                             real pos_d = rvecBuffer.buffer[buf_pos][dim2];
 935                             if (tric_dir[dim2])
 936                             {
 937                                 for (int d3 = dim2+1; d3 < DIM; d3++)
 938                                 {
 939                                     pos_d += pos[d3]*tcm[d3][dim2];
 940                                 }
 941                             }
 942
 943                             GMX_ASSERT(dim2 >= 0 && dim2 < DIM, "Keep the static analyzer happy");
 944
 945                             /* Check of we are not at the box edge.
 946                              * pbc is only handled in the first step above,
 947                              * but this check could move over pbc while
 948                              * the first step did not due to different rounding.
 949                              */
 950                             if (pos_d >= cell_x1[dim2] &&
 951                                 dd->ci[dim2] != dd->nc[dim2]-1)
 952                             {
 953                                 flag |= DD_FLAG_FW(d2);
 954                             }
 955                             else if (pos_d < cell_x0[dim2] &&
 956                                      dd->ci[dim2] != 0)
 957                             {
 958                                 flag |= DD_FLAG_BW(d2);
 959                             }
 960                             flagBuffer.buffer[cg*DD_CGIBS + 1] = flag;
 961                         }
 962                     }
 963                     /* Set to which neighboring cell this cg should go */
 964                     if (flag & DD_FLAG_FW(d2))
 965                     {
 966                         mc = d2*2;
 967                     }
 968                     else if (flag & DD_FLAG_BW(d2))
 969                     {
 970                         if (dd->nc[dd->dim[d2]] > 2)
 971                         {
 972                             mc = d2*2+1;
 973                         }
 974                         else
 975                         {
 976                             mc = d2*2;
 977                         }
 978                     }
 979                 }
 980             }
 981
 982             const int nrcg = flag & DD_FLAG_NRCG;
 983             if (mc == -1)
 984             {
 985                 /* Set the global charge group index and size */
 986                 const int globalAtomGroupIndex = flagBuffer.buffer[cg*DD_CGIBS];
 987                 dd->globalAtomGroupIndices.push_back(globalAtomGroupIndex);
 988                 dd->atomGrouping_.appendBlock(nrcg);
 989                 /* Copy the state from the buffer */
 990                 if (fr->cutoff_scheme == ecutsGROUP)
 991                 {
 992                     cg_cm = fr->cg_cm;
 993                     copy_rvec(pos, cg_cm[home_pos_cg]);
 994                 }
 995                 buf_pos++;
 996
 997                 /* Set the cginfo */
 998                 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
 999                                                    globalAtomGroupIndex);
1000                 if (comm->bLocalCG)
1001                 {
1002                     comm->bLocalCG[globalAtomGroupIndex] = TRUE;
1003                 }
1004
1005                 auto  x       = makeArrayRef(state->x);
1006                 auto  v       = makeArrayRef(state->v);
1007                 auto  cg_p    = makeArrayRef(state->cg_p);
1008                 rvec *rvecPtr = as_rvec_array(rvecBuffer.buffer.data());
1009                 for (int i = 0; i < nrcg; i++)
1010                 {
1011                     copy_rvec(rvecPtr[buf_pos++],
1012                               x[home_pos_at+i]);
1013                 }
1014                 if (bV)
1015                 {
1016                     for (int i = 0; i < nrcg; i++)
1017                     {
1018                         copy_rvec(rvecPtr[buf_pos++],
1019                                   v[home_pos_at+i]);
1020                     }
1021                 }
1022                 if (bCGP)
1023                 {
1024                     for (int i = 0; i < nrcg; i++)
1025                     {
1026                         copy_rvec(rvecPtr[buf_pos++],
1027                                   cg_p[home_pos_at+i]);
1028                     }
1029                 }
1030                 home_pos_cg += 1;
1031                 home_pos_at += nrcg;
1032             }
1033             else
1034             {
1035                 /* Reallocate the buffers if necessary  */
1036                 if ((ncg[mc] + 1)*DD_CGIBS > gmx::index(comm->cggl_flag[mc].size()))
1037                 {
1038                     comm->cggl_flag[mc].resize((ncg[mc] + 1)*DD_CGIBS);
1039                 }
1040                 size_t nvr = ncg[mc] + nat[mc]*nvec;
1041                 if (nvr + 1 + nrcg*nvec > comm->cgcm_state[mc].size())
1042                 {
1043                     comm->cgcm_state[mc].resize(nvr + 1 + nrcg*nvec);
1044                 }
1045                 /* Copy from the receive to the send buffers */
1046                 memcpy(comm->cggl_flag[mc].data() + ncg[mc]*DD_CGIBS,
1047                        flagBuffer.buffer.data() + cg*DD_CGIBS,
1048                        DD_CGIBS*sizeof(int));
1049                 memcpy(comm->cgcm_state[mc][nvr],
1050                        rvecBuffer.buffer.data() + buf_pos,
1051                        (1 + nrcg*nvec)*sizeof(rvec));
1052                 buf_pos += 1 + nrcg*nvec;
1053                 ncg[mc] += 1;
1054                 nat[mc] += nrcg;
1055             }
1056         }
1057     }
1058
1059     /* Note that the indices are now only partially up to date
1060      * and ncg_home and nat_home are not the real count, since there are
1061      * "holes" in the arrays for the charge groups that moved to neighbors.
1062      */
1063     if (fr->cutoff_scheme == ecutsVERLET)
1064     {
1065         /* We need to clear the moved flags for the received atoms,
1066          * because the moved buffer will be passed to the nbnxn gridding call.
1067          */
1068         int *moved = getMovedBuffer(comm, dd->ncg_home, home_pos_cg);
1069
1070         for (int i =  dd->ncg_home; i < home_pos_cg; i++)
1071         {
1072             moved[i] = 0;
1073         }
1074     }
1075
1076     dd->ncg_home = home_pos_cg;
1077     comm->atomRanges.setEnd(DDAtomRanges::Type::Home, home_pos_at);
1078
1079     if (fr->cutoff_scheme == ecutsGROUP)
1080     {
1081         /* We overallocated before, we need to set the right size here */
1082         dd_resize_state(state, f, comm->atomRanges.numHomeAtoms());
1083     }
1084
1085     if (debug)
1086     {
1087         fprintf(debug,
1088                 "Finished repartitioning: cgs moved out %d, new home %d\n",
1089                 *ncg_moved, dd->ncg_home-*ncg_moved);
1090
1091     }
1092 }