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.
   5  * Copyright (c) 2019, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /* \internal \file
  37  *
  38  * \brief Implements atom redistribution functions.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \ingroup module_domdec
  42  */
  43
  44 #include "gmxpre.h"
  45
  46 #include "redistribute.h"
  47
  48 #include <cstring>
  49
  50 #include "gromacs/domdec/domdec_network.h"
  51 #include "gromacs/domdec/ga2la.h"
  52 #include "gromacs/gmxlib/nrnb.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  55 #include "gromacs/mdlib/nsgrid.h"
  56 #include "gromacs/mdtypes/forcerec.h"
  57 #include "gromacs/mdtypes/nblist.h"
  58 #include "gromacs/mdtypes/state.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/gmxomp.h"
  62
  63 #include "domdec_internal.h"
  64 #include "utility.h"
  65
  66 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
  67 static constexpr int DD_CGIBS = 2;
  68
  69 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
  70
  71 /* The lower 16 bits are reserved for the charge group size */
  72 static constexpr int DD_FLAG_NRCG = 65535;
  73
  74 /* Returns which bit tells whether to move a group forward along dimension d */
  75 static inline int DD_FLAG_FW(int d)
  76 {
  77     return 1 << (16 + d * 2);
  78 }
  79
  80 /* Returns which bit tells whether to move a group backward along dimension d */
  81 static inline int DD_FLAG_BW(int d)
  82 {
  83     return 1 << (16 + d * 2 + 1);
  84 }
  85
  86 static void copyMovedAtomsToBufferPerAtom(gmx::ArrayRef<const int> move,
  87                                           int                      nvec,
  88                                           int                      vec,
  89                                           rvec*                    src,
  90                                           gmx_domdec_comm_t*       comm)
  91 {
  92     int pos_vec[DIM * 2] = { 0 };
  93
  94     for (gmx::index i = 0; i < move.ssize(); i++)
  95     {
  96         /* Skip moved atoms */
  97         const int m = move[i];
  98         if (m >= 0)
  99         {
 100             /* Copy to the communication buffer */
 101             pos_vec[m] += 1 + vec;
 102             copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
 103             pos_vec[m] += nvec - vec - 1;
 104         }
 105     }
 106 }
 107
 108 static void copyMovedUpdateGroupCogs(gmx::ArrayRef<const int>       move,
 109                                      int                            nvec,
 110                                      gmx::ArrayRef<const gmx::RVec> coordinates,
 111                                      gmx_domdec_comm_t*             comm)
 112 {
 113     int pos_vec[DIM * 2] = { 0 };
 114
 115     for (gmx::index g = 0; g < move.ssize(); g++)
 116     {
 117         /* Skip moved atoms */
 118         const int m = move[g];
 119         if (m >= 0)
 120         {
 121             /* Copy to the communication buffer */
 122             const gmx::RVec& cog =
 123                     (comm->systemInfo.useUpdateGroups ? comm->updateGroupsCog->cogForAtom(g)
 124                                                       : coordinates[g]);
 125             copy_rvec(cog, comm->cgcm_state[m][pos_vec[m]]);
 126             pos_vec[m] += 1 + nvec;
 127         }
 128     }
 129 }
 130
 131 static void clear_and_mark_ind(gmx::ArrayRef<const int> move,
 132                                gmx::ArrayRef<const int> globalAtomIndices,
 133                                gmx_ga2la_t*             ga2la,
 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             /* Signal that this atom has moved using the ns cell index.
 143              * Here we set it to -1. fill_grid will change it
 144              * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
 145              */
 146             cell_index[a] = -1;
 147         }
 148     }
 149 }
 150
 151 static void print_cg_move(FILE*               fplog,
 152                           const gmx_domdec_t* dd,
 153                           int64_t             step,
 154                           int                 cg,
 155                           int                 dim,
 156                           int                 dir,
 157                           gmx_bool            bHaveCgcmOld,
 158                           real                limitd,
 159                           rvec                cm_old,
 160                           rvec                cm_new,
 161                           real                pos_d)
 162 {
 163     const gmx_domdec_comm_t* comm = dd->comm;
 164     std::string              mesg;
 165
 166     fprintf(fplog, "\nStep %" PRId64 ":\n", step);
 167
 168     if (comm->systemInfo.useUpdateGroups)
 169     {
 170         mesg += "The update group starting at atom";
 171     }
 172     else
 173     {
 174         mesg += "Atom";
 175     }
 176     mesg += gmx::formatString(
 177             " %d moved more than the distance allowed by the domain decomposition", ddglatnr(dd, cg));
 178     if (limitd > 0)
 179     {
 180         mesg += gmx::formatString(" (%f)", limitd);
 181     }
 182     mesg += gmx::formatString(" in direction %c\n", dim2char(dim));
 183     fprintf(fplog, "%s", mesg.c_str());
 184
 185     fprintf(fplog, "distance out of cell %f\n",
 186             dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
 187     if (bHaveCgcmOld)
 188     {
 189         fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n", cm_old[XX], cm_old[YY], cm_old[ZZ]);
 190     }
 191     fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n", cm_new[XX], cm_new[YY], cm_new[ZZ]);
 192     fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n", dim2char(dim),
 193             comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
 194     fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n", dim2char(dim),
 195             comm->cell_x0[dim], comm->cell_x1[dim]);
 196 }
 197
 198 [[noreturn]] static void cg_move_error(FILE*               fplog,
 199                                        const gmx_domdec_t* dd,
 200                                        int64_t             step,
 201                                        int                 cg,
 202                                        int                 dim,
 203                                        int                 dir,
 204                                        gmx_bool            bHaveCgcmOld,
 205                                        real                limitd,
 206                                        rvec                cm_old,
 207                                        rvec                cm_new,
 208                                        real                pos_d)
 209 {
 210     if (fplog)
 211     {
 212         print_cg_move(fplog, dd, step, cg, dim, dir, bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
 213     }
 214     print_cg_move(stderr, dd, step, cg, dim, dir, bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
 215     gmx_fatal(FARGS,
 216               "One or more atoms moved too far between two domain decomposition steps.\n"
 217               "This usually means that your system is not well equilibrated");
 218 }
 219
 220 static void rotate_state_atom(t_state* state, int a)
 221 {
 222     if (state->flags & (1 << estX))
 223     {
 224         auto x = makeArrayRef(state->x);
 225         /* Rotate the complete state; for a rectangular box only */
 226         x[a][YY] = state->box[YY][YY] - x[a][YY];
 227         x[a][ZZ] = state->box[ZZ][ZZ] - x[a][ZZ];
 228     }
 229     if (state->flags & (1 << estV))
 230     {
 231         auto v   = makeArrayRef(state->v);
 232         v[a][YY] = -v[a][YY];
 233         v[a][ZZ] = -v[a][ZZ];
 234     }
 235     if (state->flags & (1 << estCGP))
 236     {
 237         auto cg_p   = makeArrayRef(state->cg_p);
 238         cg_p[a][YY] = -cg_p[a][YY];
 239         cg_p[a][ZZ] = -cg_p[a][ZZ];
 240     }
 241 }
 242
 243 /* Returns a pointer to a buffer of size \p numAtomsNew with contents 0 from numAtomsOld upwards
 244  *
 245  * Note: numAtomsOld should either be 0 or match the current buffer size.
 246  */
 247 static int* getMovedBuffer(gmx_domdec_comm_t* comm, size_t numAtomsOld, size_t numAtomsNew)
 248 {
 249     std::vector<int>& movedBuffer = comm->movedBuffer;
 250
 251     GMX_RELEASE_ASSERT(numAtomsOld == 0 || movedBuffer.size() == numAtomsOld,
 252                        "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, const ivec& dev)
 283 {
 284     int flag              = 0;
 285     int firstMoveDimValue = -1;
 286     for (int d = 0; d < dd.ndim; d++)
 287     {
 288         const int dim = dd.dim[d];
 289         if (dev[dim] == 1)
 290         {
 291             flag |= DD_FLAG_FW(d);
 292             if (firstMoveDimValue == -1)
 293             {
 294                 firstMoveDimValue = d * 2;
 295             }
 296         }
 297         else if (dev[dim] == -1)
 298         {
 299             flag |= DD_FLAG_BW(d);
 300             if (firstMoveDimValue == -1)
 301             {
 302                 if (dd.nc[dim] > 2)
 303                 {
 304                     firstMoveDimValue = d * 2 + 1;
 305                 }
 306                 else
 307                 {
 308                     firstMoveDimValue = d * 2;
 309                 }
 310             }
 311         }
 312     }
 313
 314     return firstMoveDimValue + flag;
 315 }
 316
 317 /* Determine to which atoms in the range \p cg_start, \p cg_end should go.
 318  *
 319  * Returns in the move array where the atoms should go.
 320  */
 321 static void calc_cg_move(FILE*              fplog,
 322                          int64_t            step,
 323                          gmx_domdec_t*      dd,
 324                          t_state*           state,
 325                          const ivec         tric_dir,
 326                          matrix             tcm,
 327                          const rvec         cell_x0,
 328                          const rvec         cell_x1,
 329                          const MoveLimits&  moveLimits,
 330                          int                cg_start,
 331                          int                cg_end,
 332                          gmx::ArrayRef<int> move)
 333 {
 334     const int npbcdim = dd->unitCellInfo.npbcdim;
 335     auto      x       = makeArrayRef(state->x);
 336
 337     for (int a = cg_start; a < cg_end; a++)
 338     {
 339         // TODO: Rename this center of geometry variable to cogNew
 340         rvec cm_new;
 341         copy_rvec(x[a], cm_new);
 342
 343         ivec dev = { 0 };
 344         /* Do pbc and check DD cell boundary crossings */
 345         for (int d = DIM - 1; d >= 0; d--)
 346         {
 347             if (dd->nc[d] > 1)
 348             {
 349                 bool bScrew = (dd->unitCellInfo.haveScrewPBC && d == XX);
 350                 /* Determine the location of this cg in lattice coordinates */
 351                 real pos_d = cm_new[d];
 352                 if (tric_dir[d])
 353                 {
 354                     for (int d2 = d + 1; d2 < DIM; d2++)
 355                     {
 356                         pos_d += cm_new[d2] * tcm[d2][d];
 357                     }
 358                 }
 359                 /* Put the charge group in the triclinic unit-cell */
 360                 if (pos_d >= cell_x1[d])
 361                 {
 362                     if (pos_d >= moveLimits.upper[d])
 363                     {
 364                         cg_move_error(fplog, dd, step, a, d, 1, false, moveLimits.distance[d],
 365                                       cm_new, cm_new, pos_d);
 366                     }
 367                     dev[d] = 1;
 368                     if (dd->ci[d] == dd->nc[d] - 1)
 369                     {
 370                         rvec_dec(cm_new, state->box[d]);
 371                         if (bScrew)
 372                         {
 373                             cm_new[YY] = state->box[YY][YY] - cm_new[YY];
 374                             cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
 375                         }
 376                         rvec_dec(x[a], state->box[d]);
 377                         if (bScrew)
 378                         {
 379                             rotate_state_atom(state, a);
 380                         }
 381                     }
 382                 }
 383                 else if (pos_d < cell_x0[d])
 384                 {
 385                     if (pos_d < moveLimits.lower[d])
 386                     {
 387                         cg_move_error(fplog, dd, step, a, d, -1, false, moveLimits.distance[d],
 388                                       cm_new, cm_new, pos_d);
 389                     }
 390                     dev[d] = -1;
 391                     if (dd->ci[d] == 0)
 392                     {
 393                         rvec_inc(cm_new, state->box[d]);
 394                         if (bScrew)
 395                         {
 396                             cm_new[YY] = state->box[YY][YY] - cm_new[YY];
 397                             cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
 398                         }
 399                         rvec_inc(x[a], state->box[d]);
 400                         if (bScrew)
 401                         {
 402                             rotate_state_atom(state, a);
 403                         }
 404                     }
 405                 }
 406             }
 407             else if (d < npbcdim)
 408             {
 409                 /* Put the charge group in the rectangular unit-cell */
 410                 while (cm_new[d] >= state->box[d][d])
 411                 {
 412                     rvec_dec(cm_new, state->box[d]);
 413                     rvec_dec(x[a], state->box[d]);
 414                 }
 415                 while (cm_new[d] < 0)
 416                 {
 417                     rvec_inc(cm_new, state->box[d]);
 418                     rvec_inc(x[a], state->box[d]);
 419                 }
 420             }
 421         }
 422
 423         /* Temporarily store the flag in move */
 424         move[a] = computeMoveFlag(*dd, dev);
 425     }
 426 }
 427
 428 struct PbcAndFlag
 429 {
 430     /* Constructor that purposely does not initialize anything */
 431     PbcAndFlag() {}
 432
 433     gmx::RVec pbcShift;
 434     int       moveFlag;
 435 };
 436
 437 /* Determine to which domains update groups in the range \p groupBegin, \p groupEnd should go.
 438  *
 439  * Returns in the move array where the groups should go.
 440  * Also updates the COGs and coordinates for jumps over periodic boundaries.
 441  */
 442 static void calcGroupMove(FILE*                     fplog,
 443                           int64_t                   step,
 444                           const gmx_domdec_t*       dd,
 445                           const t_state*            state,
 446                           const ivec                tric_dir,
 447                           matrix                    tcm,
 448                           const rvec                cell_x0,
 449                           const rvec                cell_x1,
 450                           const MoveLimits&         moveLimits,
 451                           int                       groupBegin,
 452                           int                       groupEnd,
 453                           gmx::ArrayRef<PbcAndFlag> pbcAndFlags)
 454 {
 455     GMX_RELEASE_ASSERT(!dd->unitCellInfo.haveScrewPBC, "Screw PBC is not supported here");
 456
 457     const int npbcdim = dd->unitCellInfo.npbcdim;
 458
 459     gmx::UpdateGroupsCog* updateGroupsCog = dd->comm->updateGroupsCog.get();
 460
 461     for (int g = groupBegin; g < groupEnd; g++)
 462     {
 463
 464         gmx::RVec& cog    = updateGroupsCog->cog(g);
 465         gmx::RVec  cogOld = cog;
 466
 467         ivec dev = { 0 };
 468         /* Do pbc and check DD cell boundary crossings */
 469         for (int d = DIM - 1; d >= 0; d--)
 470         {
 471             if (dd->nc[d] > 1)
 472             {
 473                 /* Determine the location of this COG in lattice coordinates */
 474                 real pos_d = cog[d];
 475                 if (tric_dir[d])
 476                 {
 477                     for (int d2 = d + 1; d2 < DIM; d2++)
 478                     {
 479                         pos_d += cog[d2] * tcm[d2][d];
 480                     }
 481                 }
 482                 /* Put the COG in the triclinic unit-cell */
 483                 if (pos_d >= cell_x1[d])
 484                 {
 485                     if (pos_d >= moveLimits.upper[d])
 486                     {
 487                         cg_move_error(fplog, dd, step, g, d, 1, true, moveLimits.distance[d],
 488                                       cogOld, cog, pos_d);
 489                     }
 490                     dev[d] = 1;
 491                     if (dd->ci[d] == dd->nc[d] - 1)
 492                     {
 493                         rvec_dec(cog, state->box[d]);
 494                     }
 495                 }
 496                 else if (pos_d < cell_x0[d])
 497                 {
 498                     if (pos_d < moveLimits.lower[d])
 499                     {
 500                         cg_move_error(fplog, dd, step, g, d, -1, true, moveLimits.distance[d],
 501                                       cogOld, cog, pos_d);
 502                     }
 503                     dev[d] = -1;
 504                     if (dd->ci[d] == 0)
 505                     {
 506                         rvec_inc(cog, state->box[d]);
 507                     }
 508                 }
 509             }
 510             else if (d < npbcdim)
 511             {
 512                 /* Put the COG in the rectangular unit-cell */
 513                 while (cog[d] >= state->box[d][d])
 514                 {
 515                     rvec_dec(cog, state->box[d]);
 516                 }
 517                 while (cog[d] < 0)
 518                 {
 519                     rvec_inc(cog, state->box[d]);
 520                 }
 521             }
 522         }
 523
 524         /* Store the PBC and move flag, so we can later apply them to the atoms */
 525         PbcAndFlag& pbcAndFlag = pbcAndFlags[g];
 526
 527         rvec_sub(cog, cogOld, pbcAndFlag.pbcShift);
 528         pbcAndFlag.moveFlag = computeMoveFlag(*dd, dev);
 529     }
 530 }
 531
 532 static void 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,
 549                         int64_t                      step,
 550                         gmx_domdec_t*                dd,
 551                         ivec                         tric_dir,
 552                         t_state*                     state,
 553                         PaddedHostVector<gmx::RVec>* f,
 554                         t_forcerec*                  fr,
 555                         t_nrnb*                      nrnb,
 556                         int*                         ncg_moved)
 557 {
 558     gmx_domdec_comm_t* comm = dd->comm;
 559
 560     if (dd->unitCellInfo.haveScrewPBC)
 561     {
 562         check_screw_box(state->box);
 563     }
 564
 565     // Positions are always present, so there's nothing to flag
 566     bool bV   = (state->flags & (1 << estV)) != 0;
 567     bool bCGP = (state->flags & (1 << estCGP)) != 0;
 568
 569     DDBufferAccess<int> moveBuffer(comm->intBuffer, dd->ncg_home);
 570     gmx::ArrayRef<int>  move = moveBuffer.buffer;
 571
 572     const int npbcdim = dd->unitCellInfo.npbcdim;
 573
 574     rvec       cell_x0, cell_x1;
 575     MoveLimits moveLimits;
 576     for (int d = 0; (d < DIM); d++)
 577     {
 578         moveLimits.distance[d] = dd->comm->cellsize_min[d];
 579         if (d >= npbcdim && dd->ci[d] == 0)
 580         {
 581             cell_x0[d] = -GMX_FLOAT_MAX;
 582         }
 583         else
 584         {
 585             cell_x0[d] = comm->cell_x0[d];
 586         }
 587         if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
 588         {
 589             cell_x1[d] = GMX_FLOAT_MAX;
 590         }
 591         else
 592         {
 593             cell_x1[d] = comm->cell_x1[d];
 594         }
 595         if (d < npbcdim)
 596         {
 597             moveLimits.lower[d] = comm->old_cell_x0[d] - moveLimits.distance[d];
 598             moveLimits.upper[d] = comm->old_cell_x1[d] + moveLimits.distance[d];
 599         }
 600         else
 601         {
 602             /* We check after communication if a charge group moved
 603              * more than one cell. Set the pre-comm check limit to float_max.
 604              */
 605             moveLimits.lower[d] = -GMX_FLOAT_MAX;
 606             moveLimits.upper[d] = GMX_FLOAT_MAX;
 607         }
 608     }
 609
 610     matrix tcm;
 611     make_tric_corr_matrix(npbcdim, state->box, tcm);
 612
 613     const int nthread = gmx_omp_nthreads_get(emntDomdec);
 614
 615     /* Compute the center of geometry for all home charge groups
 616      * and put them in the box and determine where they should go.
 617      */
 618     std::vector<PbcAndFlag> pbcAndFlags(
 619             comm->systemInfo.useUpdateGroups ? comm->updateGroupsCog->numCogs() : 0);
 620
 621 #pragma omp parallel num_threads(nthread)
 622     {
 623         try
 624         {
 625             const int thread = gmx_omp_get_thread_num();
 626
 627             if (comm->systemInfo.useUpdateGroups)
 628             {
 629                 const auto& updateGroupsCog = *comm->updateGroupsCog;
 630                 const int   numGroups       = updateGroupsCog.numCogs();
 631                 calcGroupMove(fplog, step, dd, state, tric_dir, tcm, cell_x0, cell_x1, moveLimits,
 632                               (thread * numGroups) / nthread, ((thread + 1) * numGroups) / nthread,
 633                               pbcAndFlags);
 634                 /* We need a barrier as atoms below can be in a COG of a different thread */
 635 #pragma omp barrier
 636                 const int numHomeAtoms = comm->atomRanges.numHomeAtoms();
 637                 applyPbcAndSetMoveFlags(updateGroupsCog, pbcAndFlags, (thread * numHomeAtoms) / nthread,
 638                                         ((thread + 1) * numHomeAtoms) / nthread, state->x, move);
 639             }
 640             else
 641             {
 642                 /* Here we handle single atoms or charge groups */
 643                 calc_cg_move(fplog, step, dd, state, tric_dir, tcm, cell_x0, cell_x1, moveLimits,
 644                              (thread * dd->ncg_home) / nthread,
 645                              ((thread + 1) * dd->ncg_home) / nthread, move);
 646             }
 647         }
 648         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 649     }
 650
 651     int ncg[DIM * 2] = { 0 };
 652     int nat[DIM * 2] = { 0 };
 653     for (int cg = 0; cg < dd->ncg_home; cg++)
 654     {
 655         if (move[cg] >= 0)
 656         {
 657             const int flag = move[cg] & ~DD_FLAG_NRCG;
 658             const int mc   = move[cg] & DD_FLAG_NRCG;
 659             move[cg]       = mc;
 660
 661             std::vector<int>& cggl_flag = comm->cggl_flag[mc];
 662
 663             /* TODO: See if we can use push_back instead */
 664             if ((ncg[mc] + 1) * DD_CGIBS > gmx::index(cggl_flag.size()))
 665             {
 666                 cggl_flag.resize((ncg[mc] + 1) * DD_CGIBS);
 667             }
 668             cggl_flag[ncg[mc] * DD_CGIBS] = dd->globalAtomGroupIndices[cg];
 669             /* We store the cg size in the lower 16 bits
 670              * and the place where the charge group should go
 671              * in the next 6 bits. This saves some communication volume.
 672              *
 673              * TODO: Remove the size, as it is now always 1.
 674              */
 675             const int numAtomsInGroup         = 1;
 676             cggl_flag[ncg[mc] * DD_CGIBS + 1] = numAtomsInGroup | flag;
 677             ncg[mc] += 1;
 678             nat[mc] += numAtomsInGroup;
 679         }
 680     }
 681
 682     inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
 683     inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
 684
 685     *ncg_moved = 0;
 686     for (int i = 0; i < dd->ndim * 2; i++)
 687     {
 688         *ncg_moved += ncg[i];
 689     }
 690
 691     int nvec = 1;
 692     if (bV)
 693     {
 694         nvec++;
 695     }
 696     if (bCGP)
 697     {
 698         nvec++;
 699     }
 700
 701     /* Make sure the communication buffers are large enough */
 702     for (int mc = 0; mc < dd->ndim * 2; mc++)
 703     {
 704         size_t nvr = ncg[mc] + nat[mc] * nvec;
 705         if (nvr > comm->cgcm_state[mc].size())
 706         {
 707             comm->cgcm_state[mc].resize(nvr);
 708         }
 709     }
 710
 711     /* With update groups we send over their COGs.
 712      * Without update groups we send the moved atom coordinates
 713      * over twice. This is so the code further down can be used
 714      * without many conditionals both with and without update groups.
 715      */
 716     copyMovedUpdateGroupCogs(move, nvec, state->x, comm);
 717
 718     int vectorIndex = 0;
 719     copyMovedAtomsToBufferPerAtom(move, nvec, vectorIndex++, state->x.rvec_array(), comm);
 720     if (bV)
 721     {
 722         copyMovedAtomsToBufferPerAtom(move, nvec, vectorIndex++, state->v.rvec_array(), comm);
 723     }
 724     if (bCGP)
 725     {
 726         copyMovedAtomsToBufferPerAtom(move, nvec, vectorIndex++, state->cg_p.rvec_array(), comm);
 727     }
 728
 729     int* moved = getMovedBuffer(comm, 0, dd->ncg_home);
 730
 731     clear_and_mark_ind(move, dd->globalAtomIndices, dd->ga2la, moved);
 732
 733     /* Now we can remove the excess global atom-group indices from the list */
 734     dd->globalAtomGroupIndices.resize(dd->ncg_home);
 735
 736     /* We reuse the intBuffer without reacquiring since we are in the same scope */
 737     DDBufferAccess<int>& flagBuffer = moveBuffer;
 738
 739     gmx::ArrayRef<const cginfo_mb_t> cginfo_mb = fr->cginfo_mb;
 740
 741     /* Temporarily store atoms passed to our rank at the end of the range */
 742     int home_pos_cg = dd->ncg_home;
 743     int home_pos_at = dd->ncg_home;
 744     for (int d = 0; d < dd->ndim; d++)
 745     {
 746         DDBufferAccess<gmx::RVec> rvecBuffer(comm->rvecBuffer, 0);
 747
 748         const int dim      = dd->dim[d];
 749         int       ncg_recv = 0;
 750         int       nvr      = 0;
 751         for (int dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
 752         {
 753             const int cdd = d * 2 + dir;
 754             /* Communicate the cg and atom counts */
 755             int sbuf[2] = { ncg[cdd], nat[cdd] };
 756             if (debug)
 757             {
 758                 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n", d, dir, sbuf[0], sbuf[1]);
 759             }
 760             int rbuf[2];
 761             ddSendrecv(dd, d, dir, sbuf, 2, rbuf, 2);
 762
 763             flagBuffer.resize((ncg_recv + rbuf[0]) * DD_CGIBS);
 764
 765             /* Communicate the charge group indices, sizes and flags */
 766             ddSendrecv(dd, d, dir, comm->cggl_flag[cdd].data(), sbuf[0] * DD_CGIBS,
 767                        flagBuffer.buffer.data() + ncg_recv * DD_CGIBS, rbuf[0] * DD_CGIBS);
 768
 769             const int nvs = ncg[cdd] + nat[cdd] * nvec;
 770             const int i   = rbuf[0] + rbuf[1] * nvec;
 771             rvecBuffer.resize(nvr + i);
 772
 773             /* Communicate cgcm and state */
 774             ddSendrecv(dd, d, dir, as_rvec_array(comm->cgcm_state[cdd].data()), nvs,
 775                        as_rvec_array(rvecBuffer.buffer.data()) + nvr, i);
 776             ncg_recv += rbuf[0];
 777             nvr += i;
 778         }
 779
 780         dd_check_alloc_ncg(fr, state, f, home_pos_cg + ncg_recv);
 781
 782         /* Process the received charge or update groups */
 783         int buf_pos = 0;
 784         for (int cg = 0; cg < ncg_recv; cg++)
 785         {
 786             /* Extract the move flags and COG for the charge or update group */
 787             int              flag = flagBuffer.buffer[cg * DD_CGIBS + 1];
 788             const gmx::RVec& cog  = rvecBuffer.buffer[buf_pos];
 789
 790             if (dim >= npbcdim && dd->nc[dim] > 2)
 791             {
 792                 /* No pbc in this dim and more than one domain boundary.
 793                  * We do a separate check if a charge group didn't move too far.
 794                  */
 795                 if (((flag & DD_FLAG_FW(d)) && cog[dim] > cell_x1[dim])
 796                     || ((flag & DD_FLAG_BW(d)) && cog[dim] < cell_x0[dim]))
 797                 {
 798                     rvec pos = { cog[0], cog[1], cog[2] };
 799                     cg_move_error(fplog, dd, step, cg, dim, (flag & DD_FLAG_FW(d)) ? 1 : 0, false,
 800                                   0, pos, pos, pos[dim]);
 801                 }
 802             }
 803
 804             int mc = -1;
 805             if (d < dd->ndim - 1)
 806             {
 807                 /* Check which direction this cg should go */
 808                 for (int d2 = d + 1; (d2 < dd->ndim && mc == -1); d2++)
 809                 {
 810                     if (isDlbOn(dd->comm))
 811                     {
 812                         /* The cell boundaries for dimension d2 are not equal
 813                          * for each cell row of the lower dimension(s),
 814                          * therefore we might need to redetermine where
 815                          * this cg should go.
 816                          */
 817                         const int dim2 = dd->dim[d2];
 818                         /* The DD grid is not staggered at the box boundaries,
 819                          * so we do not need to handle boundary crossings.
 820                          * This also means we do not have to handle PBC here.
 821                          */
 822                         if (!((dd->ci[dim2] == dd->nc[dim2] - 1 && (flag & DD_FLAG_FW(d2)))
 823                               || (dd->ci[dim2] == 0 && (flag & DD_FLAG_BW(d2)))))
 824                         {
 825                             /* Clear the two flags for this dimension */
 826                             flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
 827                             /* Determine the location of this cg
 828                              * in lattice coordinates
 829                              */
 830                             real pos_d = cog[dim2];
 831                             if (tric_dir[dim2])
 832                             {
 833                                 for (int d3 = dim2 + 1; d3 < DIM; d3++)
 834                                 {
 835                                     pos_d += cog[d3] * tcm[d3][dim2];
 836                                 }
 837                             }
 838
 839                             GMX_ASSERT(dim2 >= 0 && dim2 < DIM, "Keep the static analyzer happy");
 840
 841                             /* Check if we need to move this group
 842                              * to an adjacent cell because of the
 843                              * staggering.
 844                              */
 845                             if (pos_d >= cell_x1[dim2] && dd->ci[dim2] != dd->nc[dim2] - 1)
 846                             {
 847                                 flag |= DD_FLAG_FW(d2);
 848                             }
 849                             else if (pos_d < cell_x0[dim2] && dd->ci[dim2] != 0)
 850                             {
 851                                 flag |= DD_FLAG_BW(d2);
 852                             }
 853
 854                             flagBuffer.buffer[cg * DD_CGIBS + 1] = flag;
 855                         }
 856                     }
 857                     /* Set to which neighboring cell this cg should go */
 858                     if (flag & DD_FLAG_FW(d2))
 859                     {
 860                         mc = d2 * 2;
 861                     }
 862                     else if (flag & DD_FLAG_BW(d2))
 863                     {
 864                         if (dd->nc[dd->dim[d2]] > 2)
 865                         {
 866                             mc = d2 * 2 + 1;
 867                         }
 868                         else
 869                         {
 870                             mc = d2 * 2;
 871                         }
 872                     }
 873                 }
 874             }
 875
 876             GMX_ASSERT((flag & DD_FLAG_NRCG) == 1,
 877                        "Charge groups are gone, so all groups should have size 1");
 878             constexpr int nrcg = 1;
 879             if (mc == -1)
 880             {
 881                 /* Set the global charge group index and size */
 882                 const int globalAtomGroupIndex = flagBuffer.buffer[cg * DD_CGIBS];
 883                 dd->globalAtomGroupIndices.push_back(globalAtomGroupIndex);
 884                 /* Skip the COG entry in the buffer */
 885                 buf_pos++;
 886
 887                 /* Set the cginfo */
 888                 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb, globalAtomGroupIndex);
 889
 890                 auto  x       = makeArrayRef(state->x);
 891                 auto  v       = makeArrayRef(state->v);
 892                 auto  cg_p    = makeArrayRef(state->cg_p);
 893                 rvec* rvecPtr = as_rvec_array(rvecBuffer.buffer.data());
 894                 for (int i = 0; i < nrcg; i++)
 895                 {
 896                     copy_rvec(rvecPtr[buf_pos++], x[home_pos_at + i]);
 897                 }
 898                 if (bV)
 899                 {
 900                     for (int i = 0; i < nrcg; i++)
 901                     {
 902                         copy_rvec(rvecPtr[buf_pos++], v[home_pos_at + i]);
 903                     }
 904                 }
 905                 if (bCGP)
 906                 {
 907                     for (int i = 0; i < nrcg; i++)
 908                     {
 909                         copy_rvec(rvecPtr[buf_pos++], cg_p[home_pos_at + i]);
 910                     }
 911                 }
 912                 home_pos_cg += 1;
 913                 home_pos_at += nrcg;
 914             }
 915             else
 916             {
 917                 /* Reallocate the buffers if necessary  */
 918                 if ((ncg[mc] + 1) * DD_CGIBS > gmx::index(comm->cggl_flag[mc].size()))
 919                 {
 920                     comm->cggl_flag[mc].resize((ncg[mc] + 1) * DD_CGIBS);
 921                 }
 922                 size_t nvr = ncg[mc] + nat[mc] * nvec;
 923                 if (nvr + 1 + nrcg * nvec > comm->cgcm_state[mc].size())
 924                 {
 925                     comm->cgcm_state[mc].resize(nvr + 1 + nrcg * nvec);
 926                 }
 927                 /* Copy from the receive to the send buffers */
 928                 memcpy(comm->cggl_flag[mc].data() + ncg[mc] * DD_CGIBS,
 929                        flagBuffer.buffer.data() + cg * DD_CGIBS, DD_CGIBS * sizeof(int));
 930                 memcpy(comm->cgcm_state[mc][nvr], rvecBuffer.buffer.data() + buf_pos,
 931                        (1 + nrcg * nvec) * sizeof(rvec));
 932                 buf_pos += 1 + nrcg * nvec;
 933                 ncg[mc] += 1;
 934                 nat[mc] += nrcg;
 935             }
 936         }
 937     }
 938
 939     /* Note that the indices are now only partially up to date
 940      * and ncg_home and nat_home are not the real count, since there are
 941      * "holes" in the arrays for the charge groups that moved to neighbors.
 942      */
 943
 944     /* We need to clear the moved flags for the received atoms,
 945      * because the moved buffer will be passed to the nbnxm gridding call.
 946      */
 947     moved = getMovedBuffer(comm, dd->ncg_home, home_pos_cg);
 948
 949     for (int i = dd->ncg_home; i < home_pos_cg; i++)
 950     {
 951         moved[i] = 0;
 952     }
 953
 954     dd->ncg_home = home_pos_cg;
 955     comm->atomRanges.setEnd(DDAtomRanges::Type::Home, home_pos_at);
 956
 957     if (debug)
 958     {
 959         fprintf(debug, "Finished repartitioning: cgs moved out %d, new home %d\n", *ncg_moved,
 960                 dd->ncg_home - *ncg_moved);
 961     }
 962 }