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.
5 * Copyright (c) 2019, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * \brief Implements atom redistribution functions.
40 * \author Berk Hess <hess@kth.se>
41 * \ingroup module_domdec
46 #include "redistribute.h"
50 #include "gromacs/domdec/domdec_network.h"
51 #include "gromacs/domdec/ga2la.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/gmx_omp_nthreads.h"
55 #include "gromacs/mdlib/nsgrid.h"
56 #include "gromacs/mdtypes/forcerec.h"
57 #include "gromacs/mdtypes/nblist.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/utility/cstringutil.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/gmxomp.h"
63 #include "domdec_internal.h"
66 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
67 static constexpr int DD_CGIBS
= 2;
69 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
71 /* The lower 16 bits are reserved for the charge group size */
72 static constexpr int DD_FLAG_NRCG
= 65535;
74 /* Returns which bit tells whether to move a group forward along dimension d */
75 static inline int DD_FLAG_FW(int d
)
77 return 1 << (16 + d
* 2);
80 /* Returns which bit tells whether to move a group backward along dimension d */
81 static inline int DD_FLAG_BW(int d
)
83 return 1 << (16 + d
* 2 + 1);
86 static void copyMovedAtomsToBufferPerAtom(gmx::ArrayRef
<const int> move
,
90 gmx_domdec_comm_t
* comm
)
92 int pos_vec
[DIM
* 2] = { 0 };
94 for (gmx::index i
= 0; i
< move
.ssize(); i
++)
96 /* Skip moved atoms */
97 const int m
= move
[i
];
100 /* Copy to the communication buffer */
101 pos_vec
[m
] += 1 + vec
;
102 copy_rvec(src
[i
], comm
->cgcm_state
[m
][pos_vec
[m
]++]);
103 pos_vec
[m
] += nvec
- vec
- 1;
108 static void copyMovedUpdateGroupCogs(gmx::ArrayRef
<const int> move
,
110 gmx::ArrayRef
<const gmx::RVec
> coordinates
,
111 gmx_domdec_comm_t
* comm
)
113 int pos_vec
[DIM
* 2] = { 0 };
115 for (gmx::index g
= 0; g
< move
.ssize(); g
++)
117 /* Skip moved atoms */
118 const int m
= move
[g
];
121 /* Copy to the communication buffer */
122 const gmx::RVec
& cog
=
123 (comm
->systemInfo
.useUpdateGroups
? comm
->updateGroupsCog
->cogForAtom(g
)
125 copy_rvec(cog
, comm
->cgcm_state
[m
][pos_vec
[m
]]);
126 pos_vec
[m
] += 1 + nvec
;
131 static void clear_and_mark_ind(gmx::ArrayRef
<const int> move
,
132 gmx::ArrayRef
<const int> globalAtomIndices
,
136 for (gmx::index a
= 0; a
< move
.ssize(); a
++)
140 /* Clear the global indices */
141 ga2la
->erase(globalAtomIndices
[a
]);
142 /* Signal that this atom has moved using the ns cell index.
143 * Here we set it to -1. fill_grid will change it
144 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
151 static void print_cg_move(FILE* fplog
,
152 const gmx_domdec_t
* dd
,
157 gmx_bool bHaveCgcmOld
,
163 const gmx_domdec_comm_t
* comm
= dd
->comm
;
166 fprintf(fplog
, "\nStep %" PRId64
":\n", step
);
168 if (comm
->systemInfo
.useUpdateGroups
)
170 mesg
+= "The update group starting at atom";
176 mesg
+= gmx::formatString(
177 " %d moved more than the distance allowed by the domain decomposition", ddglatnr(dd
, cg
));
180 mesg
+= gmx::formatString(" (%f)", limitd
);
182 mesg
+= gmx::formatString(" in direction %c\n", dim2char(dim
));
183 fprintf(fplog
, "%s", mesg
.c_str());
185 fprintf(fplog
, "distance out of cell %f\n",
186 dir
== 1 ? pos_d
- comm
->cell_x1
[dim
] : pos_d
- comm
->cell_x0
[dim
]);
189 fprintf(fplog
, "Old coordinates: %8.3f %8.3f %8.3f\n", cm_old
[XX
], cm_old
[YY
], cm_old
[ZZ
]);
191 fprintf(fplog
, "New coordinates: %8.3f %8.3f %8.3f\n", cm_new
[XX
], cm_new
[YY
], cm_new
[ZZ
]);
192 fprintf(fplog
, "Old cell boundaries in direction %c: %8.3f %8.3f\n", dim2char(dim
),
193 comm
->old_cell_x0
[dim
], comm
->old_cell_x1
[dim
]);
194 fprintf(fplog
, "New cell boundaries in direction %c: %8.3f %8.3f\n", dim2char(dim
),
195 comm
->cell_x0
[dim
], comm
->cell_x1
[dim
]);
198 [[noreturn
]] static void cg_move_error(FILE* fplog
,
199 const gmx_domdec_t
* dd
,
204 gmx_bool bHaveCgcmOld
,
212 print_cg_move(fplog
, dd
, step
, cg
, dim
, dir
, bHaveCgcmOld
, limitd
, cm_old
, cm_new
, pos_d
);
214 print_cg_move(stderr
, dd
, step
, cg
, dim
, dir
, bHaveCgcmOld
, limitd
, cm_old
, cm_new
, pos_d
);
216 "One or more atoms moved too far between two domain decomposition steps.\n"
217 "This usually means that your system is not well equilibrated");
220 static void rotate_state_atom(t_state
* state
, int a
)
222 if (state
->flags
& (1 << estX
))
224 auto x
= makeArrayRef(state
->x
);
225 /* Rotate the complete state; for a rectangular box only */
226 x
[a
][YY
] = state
->box
[YY
][YY
] - x
[a
][YY
];
227 x
[a
][ZZ
] = state
->box
[ZZ
][ZZ
] - x
[a
][ZZ
];
229 if (state
->flags
& (1 << estV
))
231 auto v
= makeArrayRef(state
->v
);
232 v
[a
][YY
] = -v
[a
][YY
];
233 v
[a
][ZZ
] = -v
[a
][ZZ
];
235 if (state
->flags
& (1 << estCGP
))
237 auto cg_p
= makeArrayRef(state
->cg_p
);
238 cg_p
[a
][YY
] = -cg_p
[a
][YY
];
239 cg_p
[a
][ZZ
] = -cg_p
[a
][ZZ
];
243 /* Returns a pointer to a buffer of size \p numAtomsNew with contents 0 from numAtomsOld upwards
245 * Note: numAtomsOld should either be 0 or match the current buffer size.
247 static int* getMovedBuffer(gmx_domdec_comm_t
* comm
, size_t numAtomsOld
, size_t numAtomsNew
)
249 std::vector
<int>& movedBuffer
= comm
->movedBuffer
;
251 GMX_RELEASE_ASSERT(numAtomsOld
== 0 || movedBuffer
.size() == numAtomsOld
,
252 "numAtomsOld should either be 0 or match the current size");
254 /* Contents up to numAtomsOld should be preserved, clear from numAtomsOld */
255 if (numAtomsOld
== 0)
259 movedBuffer
.resize(numAtomsNew
);
261 return movedBuffer
.data();
264 /* Bounds to determine whether an atom group moved to far between DD steps */
267 rvec distance
; /* Maximum distance over which a group can move */
268 rvec lower
; /* Lower bound for group location */
269 rvec upper
; /* Upper bound for group location */
272 /* Compute flag that tells where to move an atom group
274 * The input dev tells whether the group moves fw,remains,bw (-1,0,1) along
276 * The return value is -1 when the group does not move.
277 * The return has move flags set when the group does move and the lower 4 bits
278 * are (mis)used to tell which is the first dimension (bit 1,2,3) the group
279 * needs to be moved along and in which direction (bit 0 not set for fw
282 static int computeMoveFlag(const gmx_domdec_t
& dd
, const ivec
& dev
)
285 int firstMoveDimValue
= -1;
286 for (int d
= 0; d
< dd
.ndim
; d
++)
288 const int dim
= dd
.dim
[d
];
291 flag
|= DD_FLAG_FW(d
);
292 if (firstMoveDimValue
== -1)
294 firstMoveDimValue
= d
* 2;
297 else if (dev
[dim
] == -1)
299 flag
|= DD_FLAG_BW(d
);
300 if (firstMoveDimValue
== -1)
302 if (dd
.numCells
[dim
] > 2)
304 firstMoveDimValue
= d
* 2 + 1;
308 firstMoveDimValue
= d
* 2;
314 return firstMoveDimValue
+ flag
;
317 /* Determine to which atoms in the range \p cg_start, \p cg_end should go.
319 * Returns in the move array where the atoms should go.
321 static void calc_cg_move(FILE* fplog
,
329 const MoveLimits
& moveLimits
,
332 gmx::ArrayRef
<int> move
)
334 const int npbcdim
= dd
->unitCellInfo
.npbcdim
;
335 auto x
= makeArrayRef(state
->x
);
337 for (int a
= cg_start
; a
< cg_end
; a
++)
339 // TODO: Rename this center of geometry variable to cogNew
341 copy_rvec(x
[a
], cm_new
);
344 /* Do pbc and check DD cell boundary crossings */
345 for (int d
= DIM
- 1; d
>= 0; d
--)
347 if (dd
->numCells
[d
] > 1)
349 bool bScrew
= (dd
->unitCellInfo
.haveScrewPBC
&& d
== XX
);
350 /* Determine the location of this cg in lattice coordinates */
351 real pos_d
= cm_new
[d
];
354 for (int d2
= d
+ 1; d2
< DIM
; d2
++)
356 pos_d
+= cm_new
[d2
] * tcm
[d2
][d
];
359 /* Put the charge group in the triclinic unit-cell */
360 if (pos_d
>= cell_x1
[d
])
362 if (pos_d
>= moveLimits
.upper
[d
])
364 cg_move_error(fplog
, dd
, step
, a
, d
, 1, false, moveLimits
.distance
[d
],
365 cm_new
, cm_new
, pos_d
);
368 if (dd
->ci
[d
] == dd
->numCells
[d
] - 1)
370 rvec_dec(cm_new
, state
->box
[d
]);
373 cm_new
[YY
] = state
->box
[YY
][YY
] - cm_new
[YY
];
374 cm_new
[ZZ
] = state
->box
[ZZ
][ZZ
] - cm_new
[ZZ
];
376 rvec_dec(x
[a
], state
->box
[d
]);
379 rotate_state_atom(state
, a
);
383 else if (pos_d
< cell_x0
[d
])
385 if (pos_d
< moveLimits
.lower
[d
])
387 cg_move_error(fplog
, dd
, step
, a
, d
, -1, false, moveLimits
.distance
[d
],
388 cm_new
, cm_new
, pos_d
);
393 rvec_inc(cm_new
, state
->box
[d
]);
396 cm_new
[YY
] = state
->box
[YY
][YY
] - cm_new
[YY
];
397 cm_new
[ZZ
] = state
->box
[ZZ
][ZZ
] - cm_new
[ZZ
];
399 rvec_inc(x
[a
], state
->box
[d
]);
402 rotate_state_atom(state
, a
);
407 else if (d
< npbcdim
)
409 /* Put the charge group in the rectangular unit-cell */
410 while (cm_new
[d
] >= state
->box
[d
][d
])
412 rvec_dec(cm_new
, state
->box
[d
]);
413 rvec_dec(x
[a
], state
->box
[d
]);
415 while (cm_new
[d
] < 0)
417 rvec_inc(cm_new
, state
->box
[d
]);
418 rvec_inc(x
[a
], state
->box
[d
]);
423 /* Temporarily store the flag in move */
424 move
[a
] = computeMoveFlag(*dd
, dev
);
430 /* Constructor that purposely does not initialize anything */
437 /* Determine to which domains update groups in the range \p groupBegin, \p groupEnd should go.
439 * Returns in the move array where the groups should go.
440 * Also updates the COGs and coordinates for jumps over periodic boundaries.
442 static void calcGroupMove(FILE* fplog
,
444 const gmx_domdec_t
* dd
,
445 const t_state
* state
,
450 const MoveLimits
& moveLimits
,
453 gmx::ArrayRef
<PbcAndFlag
> pbcAndFlags
)
455 GMX_RELEASE_ASSERT(!dd
->unitCellInfo
.haveScrewPBC
, "Screw PBC is not supported here");
457 const int npbcdim
= dd
->unitCellInfo
.npbcdim
;
459 gmx::UpdateGroupsCog
* updateGroupsCog
= dd
->comm
->updateGroupsCog
.get();
461 for (int g
= groupBegin
; g
< groupEnd
; g
++)
464 gmx::RVec
& cog
= updateGroupsCog
->cog(g
);
465 gmx::RVec cogOld
= cog
;
468 /* Do pbc and check DD cell boundary crossings */
469 for (int d
= DIM
- 1; d
>= 0; d
--)
471 if (dd
->numCells
[d
] > 1)
473 /* Determine the location of this COG in lattice coordinates */
477 for (int d2
= d
+ 1; d2
< DIM
; d2
++)
479 pos_d
+= cog
[d2
] * tcm
[d2
][d
];
482 /* Put the COG in the triclinic unit-cell */
483 if (pos_d
>= cell_x1
[d
])
485 if (pos_d
>= moveLimits
.upper
[d
])
487 cg_move_error(fplog
, dd
, step
, g
, d
, 1, true, moveLimits
.distance
[d
],
491 if (dd
->ci
[d
] == dd
->numCells
[d
] - 1)
493 rvec_dec(cog
, state
->box
[d
]);
496 else if (pos_d
< cell_x0
[d
])
498 if (pos_d
< moveLimits
.lower
[d
])
500 cg_move_error(fplog
, dd
, step
, g
, d
, -1, true, moveLimits
.distance
[d
],
506 rvec_inc(cog
, state
->box
[d
]);
510 else if (d
< npbcdim
)
512 /* Put the COG in the rectangular unit-cell */
513 while (cog
[d
] >= state
->box
[d
][d
])
515 rvec_dec(cog
, state
->box
[d
]);
519 rvec_inc(cog
, state
->box
[d
]);
524 /* Store the PBC and move flag, so we can later apply them to the atoms */
525 PbcAndFlag
& pbcAndFlag
= pbcAndFlags
[g
];
527 rvec_sub(cog
, cogOld
, pbcAndFlag
.pbcShift
);
528 pbcAndFlag
.moveFlag
= computeMoveFlag(*dd
, dev
);
532 static void applyPbcAndSetMoveFlags(const gmx::UpdateGroupsCog
& updateGroupsCog
,
533 gmx::ArrayRef
<const PbcAndFlag
> pbcAndFlags
,
536 gmx::ArrayRef
<gmx::RVec
> atomCoords
,
537 gmx::ArrayRef
<int> move
)
539 for (int a
= atomBegin
; a
< atomEnd
; a
++)
541 const PbcAndFlag
& pbcAndFlag
= pbcAndFlags
[updateGroupsCog
.cogIndex(a
)];
542 rvec_inc(atomCoords
[a
], pbcAndFlag
.pbcShift
);
543 /* Temporarily store the flag in move */
544 move
[a
] = pbcAndFlag
.moveFlag
;
548 void dd_redistribute_cg(FILE* fplog
,
553 PaddedHostVector
<gmx::RVec
>* f
,
558 gmx_domdec_comm_t
* comm
= dd
->comm
;
560 if (dd
->unitCellInfo
.haveScrewPBC
)
562 check_screw_box(state
->box
);
565 // Positions are always present, so there's nothing to flag
566 bool bV
= (state
->flags
& (1 << estV
)) != 0;
567 bool bCGP
= (state
->flags
& (1 << estCGP
)) != 0;
569 DDBufferAccess
<int> moveBuffer(comm
->intBuffer
, dd
->ncg_home
);
570 gmx::ArrayRef
<int> move
= moveBuffer
.buffer
;
572 const int npbcdim
= dd
->unitCellInfo
.npbcdim
;
574 rvec cell_x0
, cell_x1
;
575 MoveLimits moveLimits
;
576 for (int d
= 0; (d
< DIM
); d
++)
578 moveLimits
.distance
[d
] = dd
->comm
->cellsize_min
[d
];
579 if (d
>= npbcdim
&& dd
->ci
[d
] == 0)
581 cell_x0
[d
] = -GMX_FLOAT_MAX
;
585 cell_x0
[d
] = comm
->cell_x0
[d
];
587 if (d
>= npbcdim
&& dd
->ci
[d
] == dd
->numCells
[d
] - 1)
589 cell_x1
[d
] = GMX_FLOAT_MAX
;
593 cell_x1
[d
] = comm
->cell_x1
[d
];
597 moveLimits
.lower
[d
] = comm
->old_cell_x0
[d
] - moveLimits
.distance
[d
];
598 moveLimits
.upper
[d
] = comm
->old_cell_x1
[d
] + moveLimits
.distance
[d
];
602 /* We check after communication if a charge group moved
603 * more than one cell. Set the pre-comm check limit to float_max.
605 moveLimits
.lower
[d
] = -GMX_FLOAT_MAX
;
606 moveLimits
.upper
[d
] = GMX_FLOAT_MAX
;
611 make_tric_corr_matrix(npbcdim
, state
->box
, tcm
);
613 const int nthread
= gmx_omp_nthreads_get(emntDomdec
);
615 /* Compute the center of geometry for all home charge groups
616 * and put them in the box and determine where they should go.
618 std::vector
<PbcAndFlag
> pbcAndFlags(
619 comm
->systemInfo
.useUpdateGroups
? comm
->updateGroupsCog
->numCogs() : 0);
621 #pragma omp parallel num_threads(nthread)
625 const int thread
= gmx_omp_get_thread_num();
627 if (comm
->systemInfo
.useUpdateGroups
)
629 const auto& updateGroupsCog
= *comm
->updateGroupsCog
;
630 const int numGroups
= updateGroupsCog
.numCogs();
631 calcGroupMove(fplog
, step
, dd
, state
, tric_dir
, tcm
, cell_x0
, cell_x1
, moveLimits
,
632 (thread
* numGroups
) / nthread
, ((thread
+ 1) * numGroups
) / nthread
,
634 /* We need a barrier as atoms below can be in a COG of a different thread */
636 const int numHomeAtoms
= comm
->atomRanges
.numHomeAtoms();
637 applyPbcAndSetMoveFlags(updateGroupsCog
, pbcAndFlags
, (thread
* numHomeAtoms
) / nthread
,
638 ((thread
+ 1) * numHomeAtoms
) / nthread
, state
->x
, move
);
642 /* Here we handle single atoms or charge groups */
643 calc_cg_move(fplog
, step
, dd
, state
, tric_dir
, tcm
, cell_x0
, cell_x1
, moveLimits
,
644 (thread
* dd
->ncg_home
) / nthread
,
645 ((thread
+ 1) * dd
->ncg_home
) / nthread
, move
);
648 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
651 int ncg
[DIM
* 2] = { 0 };
652 int nat
[DIM
* 2] = { 0 };
653 for (int cg
= 0; cg
< dd
->ncg_home
; cg
++)
657 const int flag
= move
[cg
] & ~DD_FLAG_NRCG
;
658 const int mc
= move
[cg
] & DD_FLAG_NRCG
;
661 std::vector
<int>& cggl_flag
= comm
->cggl_flag
[mc
];
663 /* TODO: See if we can use push_back instead */
664 if ((ncg
[mc
] + 1) * DD_CGIBS
> gmx::index(cggl_flag
.size()))
666 cggl_flag
.resize((ncg
[mc
] + 1) * DD_CGIBS
);
668 cggl_flag
[ncg
[mc
] * DD_CGIBS
] = dd
->globalAtomGroupIndices
[cg
];
669 /* We store the cg size in the lower 16 bits
670 * and the place where the charge group should go
671 * in the next 6 bits. This saves some communication volume.
673 * TODO: Remove the size, as it is now always 1.
675 const int numAtomsInGroup
= 1;
676 cggl_flag
[ncg
[mc
] * DD_CGIBS
+ 1] = numAtomsInGroup
| flag
;
678 nat
[mc
] += numAtomsInGroup
;
682 inc_nrnb(nrnb
, eNR_CGCM
, comm
->atomRanges
.numHomeAtoms());
683 inc_nrnb(nrnb
, eNR_RESETX
, dd
->ncg_home
);
686 for (int i
= 0; i
< dd
->ndim
* 2; i
++)
688 *ncg_moved
+= ncg
[i
];
701 /* Make sure the communication buffers are large enough */
702 for (int mc
= 0; mc
< dd
->ndim
* 2; mc
++)
704 size_t nvr
= ncg
[mc
] + nat
[mc
] * nvec
;
705 if (nvr
> comm
->cgcm_state
[mc
].size())
707 comm
->cgcm_state
[mc
].resize(nvr
);
711 /* With update groups we send over their COGs.
712 * Without update groups we send the moved atom coordinates
713 * over twice. This is so the code further down can be used
714 * without many conditionals both with and without update groups.
716 copyMovedUpdateGroupCogs(move
, nvec
, state
->x
, comm
);
719 copyMovedAtomsToBufferPerAtom(move
, nvec
, vectorIndex
++, state
->x
.rvec_array(), comm
);
722 copyMovedAtomsToBufferPerAtom(move
, nvec
, vectorIndex
++, state
->v
.rvec_array(), comm
);
726 copyMovedAtomsToBufferPerAtom(move
, nvec
, vectorIndex
++, state
->cg_p
.rvec_array(), comm
);
729 int* moved
= getMovedBuffer(comm
, 0, dd
->ncg_home
);
731 clear_and_mark_ind(move
, dd
->globalAtomIndices
, dd
->ga2la
, moved
);
733 /* Now we can remove the excess global atom-group indices from the list */
734 dd
->globalAtomGroupIndices
.resize(dd
->ncg_home
);
736 /* We reuse the intBuffer without reacquiring since we are in the same scope */
737 DDBufferAccess
<int>& flagBuffer
= moveBuffer
;
739 gmx::ArrayRef
<const cginfo_mb_t
> cginfo_mb
= fr
->cginfo_mb
;
741 /* Temporarily store atoms passed to our rank at the end of the range */
742 int home_pos_cg
= dd
->ncg_home
;
743 int home_pos_at
= dd
->ncg_home
;
744 for (int d
= 0; d
< dd
->ndim
; d
++)
746 DDBufferAccess
<gmx::RVec
> rvecBuffer(comm
->rvecBuffer
, 0);
748 const int dim
= dd
->dim
[d
];
751 for (int dir
= 0; dir
< (dd
->numCells
[dim
] == 2 ? 1 : 2); dir
++)
753 const int cdd
= d
* 2 + dir
;
754 /* Communicate the cg and atom counts */
755 int sbuf
[2] = { ncg
[cdd
], nat
[cdd
] };
758 fprintf(debug
, "Sending ddim %d dir %d: ncg %d nat %d\n", d
, dir
, sbuf
[0], sbuf
[1]);
761 ddSendrecv(dd
, d
, dir
, sbuf
, 2, rbuf
, 2);
763 flagBuffer
.resize((ncg_recv
+ rbuf
[0]) * DD_CGIBS
);
765 /* Communicate the charge group indices, sizes and flags */
766 ddSendrecv(dd
, d
, dir
, comm
->cggl_flag
[cdd
].data(), sbuf
[0] * DD_CGIBS
,
767 flagBuffer
.buffer
.data() + ncg_recv
* DD_CGIBS
, rbuf
[0] * DD_CGIBS
);
769 const int nvs
= ncg
[cdd
] + nat
[cdd
] * nvec
;
770 const int i
= rbuf
[0] + rbuf
[1] * nvec
;
771 rvecBuffer
.resize(nvr
+ i
);
773 /* Communicate cgcm and state */
774 ddSendrecv(dd
, d
, dir
, as_rvec_array(comm
->cgcm_state
[cdd
].data()), nvs
,
775 as_rvec_array(rvecBuffer
.buffer
.data()) + nvr
, i
);
780 dd_check_alloc_ncg(fr
, state
, f
, home_pos_cg
+ ncg_recv
);
782 /* Process the received charge or update groups */
784 for (int cg
= 0; cg
< ncg_recv
; cg
++)
786 /* Extract the move flags and COG for the charge or update group */
787 int flag
= flagBuffer
.buffer
[cg
* DD_CGIBS
+ 1];
788 const gmx::RVec
& cog
= rvecBuffer
.buffer
[buf_pos
];
790 if (dim
>= npbcdim
&& dd
->numCells
[dim
] > 2)
792 /* No pbc in this dim and more than one domain boundary.
793 * We do a separate check if a charge group didn't move too far.
795 if (((flag
& DD_FLAG_FW(d
)) && cog
[dim
] > cell_x1
[dim
])
796 || ((flag
& DD_FLAG_BW(d
)) && cog
[dim
] < cell_x0
[dim
]))
798 rvec pos
= { cog
[0], cog
[1], cog
[2] };
799 cg_move_error(fplog
, dd
, step
, cg
, dim
, (flag
& DD_FLAG_FW(d
)) ? 1 : 0, false,
800 0, pos
, pos
, pos
[dim
]);
805 if (d
< dd
->ndim
- 1)
807 /* Check which direction this cg should go */
808 for (int d2
= d
+ 1; (d2
< dd
->ndim
&& mc
== -1); d2
++)
810 if (isDlbOn(dd
->comm
))
812 /* The cell boundaries for dimension d2 are not equal
813 * for each cell row of the lower dimension(s),
814 * therefore we might need to redetermine where
817 const int dim2
= dd
->dim
[d2
];
818 /* The DD grid is not staggered at the box boundaries,
819 * so we do not need to handle boundary crossings.
820 * This also means we do not have to handle PBC here.
822 if (!((dd
->ci
[dim2
] == dd
->numCells
[dim2
] - 1 && (flag
& DD_FLAG_FW(d2
)))
823 || (dd
->ci
[dim2
] == 0 && (flag
& DD_FLAG_BW(d2
)))))
825 /* Clear the two flags for this dimension */
826 flag
&= ~(DD_FLAG_FW(d2
) | DD_FLAG_BW(d2
));
827 /* Determine the location of this cg
828 * in lattice coordinates
830 real pos_d
= cog
[dim2
];
833 for (int d3
= dim2
+ 1; d3
< DIM
; d3
++)
835 pos_d
+= cog
[d3
] * tcm
[d3
][dim2
];
839 GMX_ASSERT(dim2
>= 0 && dim2
< DIM
, "Keep the static analyzer happy");
841 /* Check if we need to move this group
842 * to an adjacent cell because of the
845 if (pos_d
>= cell_x1
[dim2
] && dd
->ci
[dim2
] != dd
->numCells
[dim2
] - 1)
847 flag
|= DD_FLAG_FW(d2
);
849 else if (pos_d
< cell_x0
[dim2
] && dd
->ci
[dim2
] != 0)
851 flag
|= DD_FLAG_BW(d2
);
854 flagBuffer
.buffer
[cg
* DD_CGIBS
+ 1] = flag
;
857 /* Set to which neighboring cell this cg should go */
858 if (flag
& DD_FLAG_FW(d2
))
862 else if (flag
& DD_FLAG_BW(d2
))
864 if (dd
->numCells
[dd
->dim
[d2
]] > 2)
876 GMX_ASSERT((flag
& DD_FLAG_NRCG
) == 1,
877 "Charge groups are gone, so all groups should have size 1");
878 constexpr int nrcg
= 1;
881 /* Set the global charge group index and size */
882 const int globalAtomGroupIndex
= flagBuffer
.buffer
[cg
* DD_CGIBS
];
883 dd
->globalAtomGroupIndices
.push_back(globalAtomGroupIndex
);
884 /* Skip the COG entry in the buffer */
888 fr
->cginfo
[home_pos_cg
] = ddcginfo(cginfo_mb
, globalAtomGroupIndex
);
890 auto x
= makeArrayRef(state
->x
);
891 auto v
= makeArrayRef(state
->v
);
892 auto cg_p
= makeArrayRef(state
->cg_p
);
893 rvec
* rvecPtr
= as_rvec_array(rvecBuffer
.buffer
.data());
894 for (int i
= 0; i
< nrcg
; i
++)
896 copy_rvec(rvecPtr
[buf_pos
++], x
[home_pos_at
+ i
]);
900 for (int i
= 0; i
< nrcg
; i
++)
902 copy_rvec(rvecPtr
[buf_pos
++], v
[home_pos_at
+ i
]);
907 for (int i
= 0; i
< nrcg
; i
++)
909 copy_rvec(rvecPtr
[buf_pos
++], cg_p
[home_pos_at
+ i
]);
917 /* Reallocate the buffers if necessary */
918 if ((ncg
[mc
] + 1) * DD_CGIBS
> gmx::index(comm
->cggl_flag
[mc
].size()))
920 comm
->cggl_flag
[mc
].resize((ncg
[mc
] + 1) * DD_CGIBS
);
922 size_t nvr
= ncg
[mc
] + nat
[mc
] * nvec
;
923 if (nvr
+ 1 + nrcg
* nvec
> comm
->cgcm_state
[mc
].size())
925 comm
->cgcm_state
[mc
].resize(nvr
+ 1 + nrcg
* nvec
);
927 /* Copy from the receive to the send buffers */
928 memcpy(comm
->cggl_flag
[mc
].data() + ncg
[mc
] * DD_CGIBS
,
929 flagBuffer
.buffer
.data() + cg
* DD_CGIBS
, DD_CGIBS
* sizeof(int));
930 memcpy(comm
->cgcm_state
[mc
][nvr
], rvecBuffer
.buffer
.data() + buf_pos
,
931 (1 + nrcg
* nvec
) * sizeof(rvec
));
932 buf_pos
+= 1 + nrcg
* nvec
;
939 /* Note that the indices are now only partially up to date
940 * and ncg_home and nat_home are not the real count, since there are
941 * "holes" in the arrays for the charge groups that moved to neighbors.
944 /* We need to clear the moved flags for the received atoms,
945 * because the moved buffer will be passed to the nbnxm gridding call.
947 moved
= getMovedBuffer(comm
, dd
->ncg_home
, home_pos_cg
);
949 for (int i
= dd
->ncg_home
; i
< home_pos_cg
; i
++)
954 dd
->ncg_home
= home_pos_cg
;
955 comm
->atomRanges
.setEnd(DDAtomRanges::Type::Home
, home_pos_at
);
959 fprintf(debug
, "Finished repartitioning: cgs moved out %d, new home %d\n", *ncg_moved
,
960 dd
->ncg_home
- *ncg_moved
);