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