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