2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * \brief Implements atom redistribution functions.
39 * \author Berk Hess <hess@kth.se>
40 * \ingroup module_domdec
45 #include "redistribute.h"
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"
62 #include "domdec_internal.h"
65 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
66 static constexpr int DD_CGIBS = 2;
68 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
70 /* The lower 16 bits are reserved for the charge group size */
71 static constexpr int DD_FLAG_NRCG = 65535;
73 /* Returns which bit tells whether to move a group forward along dimension d */
74 static inline int DD_FLAG_FW(int d)
76 return 1 << (16 + d*2);
79 /* Returns which bit tells whether to move a group backward along dimension d */
80 static inline int DD_FLAG_BW(int d)
82 return 1 << (16 + d*2 + 1);
86 copyMovedAtomsToBufferPerAtom(gmx::ArrayRef<const int> move,
87 const gmx::RangePartitioning &atomGroups,
89 rvec *src, gmx_domdec_comm_t *comm)
91 int pos_vec[DIM*2] = { 0 };
93 for (int g = 0; g < move.ssize(); g++)
95 const auto atomGroup = atomGroups.block(g);
96 /* Skip moved atoms */
100 /* Copy to the communication buffer */
101 int numAtoms = atomGroup.size();
102 pos_vec[m] += 1 + vec*numAtoms;
103 for (int i : atomGroup)
105 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
107 pos_vec[m] += (nvec - vec - 1)*numAtoms;
113 copyMovedChargeGroupCogs(gmx::ArrayRef<const int> move,
114 const gmx::RangePartitioning &atomGroups,
115 int nvec, const rvec *src,
116 gmx_domdec_comm_t *comm)
118 int pos_vec[DIM*2] = { 0 };
120 for (int g = 0; g < move.ssize(); g++)
122 const auto atomGroup = atomGroups.block(g);
123 /* Skip moved atoms */
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;
135 copyMovedUpdateGroupCogs(gmx::ArrayRef<const int> move,
136 int nvec, gmx::ArrayRef<const gmx::RVec> coordinates,
137 gmx_domdec_comm_t *comm)
139 int pos_vec[DIM*2] = { 0 };
141 for (int g = 0; g < move.ssize(); g++)
143 /* Skip moved atoms */
144 const int m = move[g];
147 /* Copy to the communication buffer */
148 const gmx::RVec &cog = (comm->useUpdateGroups ?
149 comm->updateGroupsCog->cogForAtom(g) :
151 copy_rvec(cog, comm->cgcm_state[m][pos_vec[m]]);
152 pos_vec[m] += 1 + nvec;
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,
165 for (int g = 0; g < move.ssize(); g++)
169 /* Clear the global indices */
170 for (int a : atomGroups.block(g))
172 ga2la->erase(globalAtomIndices[a]);
176 bLocalCG[globalAtomGroupIndices[g]] = FALSE;
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.
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)
193 const gmx_domdec_comm_t *comm = dd->comm;
196 fprintf(fplog, "\nStep %" PRId64 ":\n", step);
200 mesg += "The charge group starting at atom";
202 else if (comm->useUpdateGroups)
204 mesg += "The update group starting at atom";
210 mesg += gmx::formatString(" %d moved more than the distance allowed by the domain decomposition",
211 ddglatnr(dd, dd->atomGrouping().block(cg).begin()));
214 mesg += gmx::formatString(" (%f)", limitd);
216 mesg += gmx::formatString(" in direction %c\n", dim2char(dim));
217 fprintf(fplog, "%s", mesg.c_str());
219 fprintf(fplog, "distance out of cell %f\n",
220 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
223 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
224 cm_old[XX], cm_old[YY], cm_old[ZZ]);
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",
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",
233 comm->cell_x0[dim], comm->cell_x1[dim]);
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)
245 print_cg_move(fplog, dd, step, cg, dim, dir,
246 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
248 print_cg_move(stderr, dd, step, cg, dim, dir,
249 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
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");
255 static void rotate_state_atom(t_state *state, int a)
257 if (state->flags & (1 << estX))
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];
264 if (state->flags & (1 << estV))
266 auto v = makeArrayRef(state->v);
267 v[a][YY] = -v[a][YY];
268 v[a][ZZ] = -v[a][ZZ];
270 if (state->flags & (1 << estCGP))
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];
278 /* Returns a pointer to a buffer of size \p numAtomsNew with contents 0 from numAtomsOld upwards
280 * Note: numAtomsOld should either be 0 or match the current buffer size.
282 static int *getMovedBuffer(gmx_domdec_comm_t *comm,
286 std::vector<int> &movedBuffer = comm->movedBuffer;
288 GMX_RELEASE_ASSERT(numAtomsOld == 0 || movedBuffer.size() == numAtomsOld, "numAtomsOld should either be 0 or match the current size");
290 /* Contents up to numAtomsOld should be preserved, clear from numAtomsOld */
291 if (numAtomsOld == 0)
295 movedBuffer.resize(numAtomsNew);
297 return movedBuffer.data();
300 /* Bounds to determine whether an atom group moved to far between DD steps */
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 */
308 /* Compute flag that tells where to move an atom group
310 * The input dev tells whether the group moves fw,remains,bw (-1,0,1) along
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
318 static int computeMoveFlag(const gmx_domdec_t &dd,
322 int firstMoveDimValue = -1;
323 for (int d = 0; d < dd.ndim; d++)
325 const int dim = dd.dim[d];
328 flag |= DD_FLAG_FW(d);
329 if (firstMoveDimValue == -1)
331 firstMoveDimValue = d*2;
334 else if (dev[dim] == -1)
336 flag |= DD_FLAG_BW(d);
337 if (firstMoveDimValue == -1)
341 firstMoveDimValue = d*2 + 1;
345 firstMoveDimValue = d*2;
351 return firstMoveDimValue + flag;
354 /* Determine to which domains atomGroups in the range \p cg_start, \p cg_end should go.
356 * Returns in the move array where the groups should go.
357 * Also updates \p cg_cm for jumps over periodic boundaries.
359 * \TODO Rename cg to atomGroup.
361 static void calc_cg_move(FILE *fplog, int64_t step,
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,
370 gmx::ArrayRef<int> move)
372 const int npbcdim = dd->npbcdim;
373 auto x = makeArrayRef(state->x);
375 for (int g = cg_start; g < cg_end; g++)
377 const auto atomGroup = atomGroups.block(g);
378 const int numAtoms = atomGroup.size();
379 // TODO: Rename this center of geometry variable to cogNew
383 copy_rvec(x[atomGroup.begin()], cm_new);
387 real invNumAtoms = 1/static_cast<real>(numAtoms);
389 for (int k : atomGroup)
391 rvec_inc(cm_new, x[k]);
393 for (int d = 0; d < DIM; d++)
395 cm_new[d] = invNumAtoms*cm_new[d];
400 /* Do pbc and check DD cell boundary crossings */
401 for (int d = DIM - 1; d >= 0; d--)
405 bool bScrew = (dd->bScrewPBC && d == XX);
406 /* Determine the location of this cg in lattice coordinates */
407 real pos_d = cm_new[d];
410 for (int d2 = d + 1; d2 < DIM; d2++)
412 pos_d += cm_new[d2]*tcm[d2][d];
415 /* Put the charge group in the triclinic unit-cell */
416 if (pos_d >= cell_x1[d])
418 if (pos_d >= moveLimits.upper[d])
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);
425 if (dd->ci[d] == dd->nc[d] - 1)
427 rvec_dec(cm_new, state->box[d]);
430 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
431 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
433 for (int k : atomGroup)
435 rvec_dec(x[k], state->box[d]);
438 rotate_state_atom(state, k);
443 else if (pos_d < cell_x0[d])
445 if (pos_d < moveLimits.lower[d])
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);
454 rvec_inc(cm_new, state->box[d]);
457 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
458 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
460 for (int k : atomGroup)
462 rvec_inc(x[k], state->box[d]);
465 rotate_state_atom(state, k);
471 else if (d < npbcdim)
473 /* Put the charge group in the rectangular unit-cell */
474 while (cm_new[d] >= state->box[d][d])
476 rvec_dec(cm_new, state->box[d]);
477 for (int k : atomGroup)
479 rvec_dec(x[k], state->box[d]);
482 while (cm_new[d] < 0)
484 rvec_inc(cm_new, state->box[d]);
485 for (int k : atomGroup)
487 rvec_inc(x[k], state->box[d]);
493 copy_rvec(cm_new, cg_cm[g]);
495 /* Temporarily store the flag in move */
496 move[g] = computeMoveFlag(*dd, dev);
502 /* Constructor that purposely does not initialize anything */
511 /* Determine to which domains update groups in the range \p groupBegin, \p groupEnd should go.
513 * Returns in the move array where the groups should go.
514 * Also updates the COGs and coordinates for jumps over periodic boundaries.
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)
525 GMX_RELEASE_ASSERT(!dd->bScrewPBC, "Screw PBC is not supported here");
527 const int npbcdim = dd->npbcdim;
529 gmx::UpdateGroupsCog *updateGroupsCog = dd->comm->updateGroupsCog.get();
531 for (int g = groupBegin; g < groupEnd; g++)
534 gmx::RVec &cog = updateGroupsCog->cog(g);
535 gmx::RVec cogOld = cog;
538 /* Do pbc and check DD cell boundary crossings */
539 for (int d = DIM - 1; d >= 0; d--)
543 /* Determine the location of this COG in lattice coordinates */
547 for (int d2 = d + 1; d2 < DIM; d2++)
549 pos_d += cog[d2]*tcm[d2][d];
552 /* Put the COG in the triclinic unit-cell */
553 if (pos_d >= cell_x1[d])
555 if (pos_d >= moveLimits.upper[d])
557 cg_move_error(fplog, dd, step, g, d, 1,
558 true, moveLimits.distance[d],
562 if (dd->ci[d] == dd->nc[d] - 1)
564 rvec_dec(cog, state->box[d]);
567 else if (pos_d < cell_x0[d])
569 if (pos_d < moveLimits.lower[d])
571 cg_move_error(fplog, dd, step, g, d, -1,
572 true, moveLimits.distance[d],
578 rvec_inc(cog, state->box[d]);
582 else if (d < npbcdim)
584 /* Put the COG in the rectangular unit-cell */
585 while (cog[d] >= state->box[d][d])
587 rvec_dec(cog, state->box[d]);
591 rvec_inc(cog, state->box[d]);
596 /* Store the PBC and move flag, so we can later apply them to the atoms */
597 PbcAndFlag &pbcAndFlag = pbcAndFlags[g];
599 rvec_sub(cog, cogOld, pbcAndFlag.pbcShift);
600 pbcAndFlag.moveFlag = computeMoveFlag(*dd, dev);
605 applyPbcAndSetMoveFlags(const gmx::UpdateGroupsCog &updateGroupsCog,
606 gmx::ArrayRef<const PbcAndFlag> pbcAndFlags,
609 gmx::ArrayRef<gmx::RVec> atomCoords,
610 gmx::ArrayRef<int> move)
612 for (int a = atomBegin; a < atomEnd; a++)
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;
621 void dd_redistribute_cg(FILE *fplog, int64_t step,
622 gmx_domdec_t *dd, ivec tric_dir,
624 PaddedVector<gmx::RVec> *f,
629 gmx_domdec_comm_t *comm = dd->comm;
633 check_screw_box(state->box);
636 rvec *cg_cm = nullptr;
637 if (fr->cutoff_scheme == ecutsGROUP)
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;
646 DDBufferAccess<int> moveBuffer(comm->intBuffer, dd->ncg_home);
647 gmx::ArrayRef<int> move = moveBuffer.buffer;
649 const int npbcdim = dd->npbcdim;
651 rvec cell_x0, cell_x1;
652 MoveLimits moveLimits;
653 for (int d = 0; (d < DIM); d++)
655 moveLimits.distance[d] = dd->comm->cellsize_min[d];
656 if (d >= npbcdim && dd->ci[d] == 0)
658 cell_x0[d] = -GMX_FLOAT_MAX;
662 cell_x0[d] = comm->cell_x0[d];
664 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
666 cell_x1[d] = GMX_FLOAT_MAX;
670 cell_x1[d] = comm->cell_x1[d];
674 moveLimits.lower[d] = comm->old_cell_x0[d] - moveLimits.distance[d];
675 moveLimits.upper[d] = comm->old_cell_x1[d] + moveLimits.distance[d];
679 /* We check after communication if a charge group moved
680 * more than one cell. Set the pre-comm check limit to float_max.
682 moveLimits.lower[d] = -GMX_FLOAT_MAX;
683 moveLimits.upper[d] = GMX_FLOAT_MAX;
688 make_tric_corr_matrix(npbcdim, state->box, tcm);
690 const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
692 const int nthread = gmx_omp_nthreads_get(emntDomdec);
694 /* Compute the center of geometry for all home charge groups
695 * and put them in the box and determine where they should go.
697 std::vector<PbcAndFlag> pbcAndFlags(comm->useUpdateGroups ? comm->updateGroupsCog->numCogs() : 0);
699 #pragma omp parallel num_threads(nthread)
703 const int thread = gmx_omp_get_thread_num();
705 if (comm->useUpdateGroups)
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,
714 /* We need a barrier as atoms below can be in a COG of a different thread */
716 const int numHomeAtoms = comm->atomRanges.numHomeAtoms();
717 applyPbcAndSetMoveFlags(updateGroupsCog, pbcAndFlags,
718 ( thread *numHomeAtoms)/nthread,
719 ((thread+1)*numHomeAtoms)/nthread,
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,
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()),
735 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
738 int ncg[DIM*2] = { 0 };
739 int nat[DIM*2] = { 0 };
740 for (int cg = 0; cg < dd->ncg_home; cg++)
744 const int flag = move[cg] & ~DD_FLAG_NRCG;
745 const int mc = move[cg] & DD_FLAG_NRCG;
748 std::vector<int> &cggl_flag = comm->cggl_flag[mc];
750 /* TODO: See if we can use push_back instead */
751 if ((ncg[mc] + 1)*DD_CGIBS > gmx::index(cggl_flag.size()))
753 cggl_flag.resize((ncg[mc] + 1)*DD_CGIBS);
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.
760 const int nrcg = atomGrouping.block(cg).size();
761 cggl_flag[ncg[mc]*DD_CGIBS+1] = nrcg | flag;
767 inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
768 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
771 for (int i = 0; i < dd->ndim*2; i++)
773 *ncg_moved += ncg[i];
786 /* Make sure the communication buffers are large enough */
787 for (int mc = 0; mc < dd->ndim*2; mc++)
789 size_t nvr = ncg[mc] + nat[mc]*nvec;
790 if (nvr > comm->cgcm_state[mc].size())
792 comm->cgcm_state[mc].resize(nvr);
796 switch (fr->cutoff_scheme)
799 /* Recalculating cg_cm might be cheaper than communicating,
800 * but that could give rise to rounding issues.
802 copyMovedChargeGroupCogs(move, dd->atomGrouping(),
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.
811 copyMovedUpdateGroupCogs(move, nvec, state->x, comm);
814 gmx_incons("unimplemented");
818 copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
820 state->x.rvec_array(),
824 copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
826 state->v.rvec_array(),
831 copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
833 state->cg_p.rvec_array(),
837 int *moved = getMovedBuffer(comm, 0, dd->ncg_home);
839 clear_and_mark_ind(move,
840 dd->globalAtomGroupIndices, dd->atomGrouping(), dd->globalAtomIndices,
841 dd->ga2la, comm->bLocalCG,
844 /* Now we can remove the excess global atom-group indices from the list */
845 dd->globalAtomGroupIndices.resize(dd->ncg_home);
846 dd->atomGrouping_.reduceNumBlocks(dd->ncg_home);
848 /* We reuse the intBuffer without reacquiring since we are in the same scope */
849 DDBufferAccess<int> &flagBuffer = moveBuffer;
851 const cginfo_mb_t *cginfo_mb = fr->cginfo_mb;
853 /* Temporarily store atoms passed to our rank at the end of the range */
854 int home_pos_cg = dd->ncg_home;
855 int home_pos_at = dd->atomGrouping().subRange(0, dd->ncg_home).end();
856 for (int d = 0; d < dd->ndim; d++)
858 DDBufferAccess<gmx::RVec> rvecBuffer(comm->rvecBuffer, 0);
860 const int dim = dd->dim[d];
863 for (int dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
865 const int cdd = d*2 + dir;
866 /* Communicate the cg and atom counts */
867 int sbuf[2] = { ncg[cdd], nat[cdd] };
870 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
871 d, dir, sbuf[0], sbuf[1]);
874 ddSendrecv(dd, d, dir, sbuf, 2, rbuf, 2);
876 flagBuffer.resize((ncg_recv + rbuf[0])*DD_CGIBS);
878 /* Communicate the charge group indices, sizes and flags */
879 ddSendrecv(dd, d, dir,
880 comm->cggl_flag[cdd].data(), sbuf[0]*DD_CGIBS,
881 flagBuffer.buffer.data() + ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
883 const int nvs = ncg[cdd] + nat[cdd]*nvec;
884 const int i = rbuf[0] + rbuf[1] *nvec;
885 rvecBuffer.resize(nvr + i);
887 /* Communicate cgcm and state */
888 ddSendrecv(dd, d, dir,
889 as_rvec_array(comm->cgcm_state[cdd].data()), nvs,
890 as_rvec_array(rvecBuffer.buffer.data()) + nvr, i);
895 dd_check_alloc_ncg(fr, state, f, home_pos_cg + ncg_recv);
896 if (fr->cutoff_scheme == ecutsGROUP)
898 /* Here we resize to more than necessary and shrink later */
899 dd_resize_state(state, f, home_pos_at + ncg_recv*MAX_CGCGSIZE);
902 /* Process the received charge or update groups */
904 for (int cg = 0; cg < ncg_recv; cg++)
906 /* Extract the move flags and COG for the charge or update group */
907 int flag = flagBuffer.buffer[cg*DD_CGIBS + 1];
908 const gmx::RVec &cog = rvecBuffer.buffer[buf_pos];
910 if (dim >= npbcdim && dd->nc[dim] > 2)
912 /* No pbc in this dim and more than one domain boundary.
913 * We do a separate check if a charge group didn't move too far.
915 if (((flag & DD_FLAG_FW(d)) &&
916 cog[dim] > cell_x1[dim]) ||
917 ((flag & DD_FLAG_BW(d)) &&
918 cog[dim] < cell_x0[dim]))
920 rvec pos = { cog[0], cog[1], cog[2] };
921 cg_move_error(fplog, dd, step, cg, dim,
922 (flag & DD_FLAG_FW(d)) ? 1 : 0,
923 fr->cutoff_scheme == ecutsGROUP, 0,
931 /* Check which direction this cg should go */
932 for (int d2 = d + 1; (d2 < dd->ndim && mc == -1); d2++)
934 if (isDlbOn(dd->comm))
936 /* The cell boundaries for dimension d2 are not equal
937 * for each cell row of the lower dimension(s),
938 * therefore we might need to redetermine where
941 const int dim2 = dd->dim[d2];
942 /* The DD grid is not staggered at the box boundaries,
943 * so we do not need to handle boundary crossings.
944 * This also means we do not have to handle PBC here.
946 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
947 (flag & DD_FLAG_FW(d2))) ||
948 (dd->ci[dim2] == 0 &&
949 (flag & DD_FLAG_BW(d2)))))
951 /* Clear the two flags for this dimension */
952 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
953 /* Determine the location of this cg
954 * in lattice coordinates
956 real pos_d = cog[dim2];
959 for (int d3 = dim2+1; d3 < DIM; d3++)
961 pos_d += cog[d3]*tcm[d3][dim2];
965 GMX_ASSERT(dim2 >= 0 && dim2 < DIM, "Keep the static analyzer happy");
967 /* Check if we need to move this group
968 * to an adjacent cell because of the
971 if (pos_d >= cell_x1[dim2] &&
972 dd->ci[dim2] != dd->nc[dim2]-1)
974 flag |= DD_FLAG_FW(d2);
976 else if (pos_d < cell_x0[dim2] &&
979 flag |= DD_FLAG_BW(d2);
982 flagBuffer.buffer[cg*DD_CGIBS + 1] = flag;
985 /* Set to which neighboring cell this cg should go */
986 if (flag & DD_FLAG_FW(d2))
990 else if (flag & DD_FLAG_BW(d2))
992 if (dd->nc[dd->dim[d2]] > 2)
1004 const int nrcg = flag & DD_FLAG_NRCG;
1007 /* Set the global charge group index and size */
1008 const int globalAtomGroupIndex = flagBuffer.buffer[cg*DD_CGIBS];
1009 dd->globalAtomGroupIndices.push_back(globalAtomGroupIndex);
1010 dd->atomGrouping_.appendBlock(nrcg);
1011 /* Copy the state from the buffer */
1012 if (fr->cutoff_scheme == ecutsGROUP)
1015 copy_rvec(cog, cg_cm[home_pos_cg]);
1019 /* Set the cginfo */
1020 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
1021 globalAtomGroupIndex);
1024 comm->bLocalCG[globalAtomGroupIndex] = TRUE;
1027 auto x = makeArrayRef(state->x);
1028 auto v = makeArrayRef(state->v);
1029 auto cg_p = makeArrayRef(state->cg_p);
1030 rvec *rvecPtr = as_rvec_array(rvecBuffer.buffer.data());
1031 for (int i = 0; i < nrcg; i++)
1033 copy_rvec(rvecPtr[buf_pos++],
1038 for (int i = 0; i < nrcg; i++)
1040 copy_rvec(rvecPtr[buf_pos++],
1046 for (int i = 0; i < nrcg; i++)
1048 copy_rvec(rvecPtr[buf_pos++],
1049 cg_p[home_pos_at+i]);
1053 home_pos_at += nrcg;
1057 /* Reallocate the buffers if necessary */
1058 if ((ncg[mc] + 1)*DD_CGIBS > gmx::index(comm->cggl_flag[mc].size()))
1060 comm->cggl_flag[mc].resize((ncg[mc] + 1)*DD_CGIBS);
1062 size_t nvr = ncg[mc] + nat[mc]*nvec;
1063 if (nvr + 1 + nrcg*nvec > comm->cgcm_state[mc].size())
1065 comm->cgcm_state[mc].resize(nvr + 1 + nrcg*nvec);
1067 /* Copy from the receive to the send buffers */
1068 memcpy(comm->cggl_flag[mc].data() + ncg[mc]*DD_CGIBS,
1069 flagBuffer.buffer.data() + cg*DD_CGIBS,
1070 DD_CGIBS*sizeof(int));
1071 memcpy(comm->cgcm_state[mc][nvr],
1072 rvecBuffer.buffer.data() + buf_pos,
1073 (1 + nrcg*nvec)*sizeof(rvec));
1074 buf_pos += 1 + nrcg*nvec;
1081 /* Note that the indices are now only partially up to date
1082 * and ncg_home and nat_home are not the real count, since there are
1083 * "holes" in the arrays for the charge groups that moved to neighbors.
1085 if (fr->cutoff_scheme == ecutsVERLET)
1087 /* We need to clear the moved flags for the received atoms,
1088 * because the moved buffer will be passed to the nbnxn gridding call.
1090 int *moved = getMovedBuffer(comm, dd->ncg_home, home_pos_cg);
1092 for (int i = dd->ncg_home; i < home_pos_cg; i++)
1098 dd->ncg_home = home_pos_cg;
1099 comm->atomRanges.setEnd(DDAtomRanges::Type::Home, home_pos_at);
1101 if (fr->cutoff_scheme == ecutsGROUP)
1103 /* We overallocated before, we need to set the right size here */
1104 dd_resize_state(state, f, comm->atomRanges.numHomeAtoms());
1110 "Finished repartitioning: cgs moved out %d, new home %d\n",
1111 *ncg_moved, dd->ncg_home-*ncg_moved);