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