src/gromacs/domdec/redistribute.cpp

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