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