?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Partial conversion of gmx_domdec_t to C++
[gromacs.git] / src / gromacs / domdec / distribute.cpp
blobfd374eef215c5d249178de2d5a5293bcfe20db25
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2018, 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 distribution functions.
38 *
39 * \author Berk Hess <hess@kth.se>
40 * \ingroup module_domdec
41 */
42
43 #include "gmxpre.h"
44
45 #include "distribute.h"
46
47 #include "config.h"
48
49 #include <vector>
50
51 #include "gromacs/domdec/domdec_network.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/commrec.h"
54 #include "gromacs/mdtypes/df_history.h"
55 #include "gromacs/mdtypes/state.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
58
59 #include "atomdistribution.h"
60 #include "cellsizes.h"
61 #include "domdec_internal.h"
62 #include "utility.h"
63
64 static void distributeVecSendrecv(gmx_domdec_t *dd,
65 const t_block *cgs,
66 gmx::ArrayRef<const gmx::RVec> globalVec,
67 gmx::ArrayRef<gmx::RVec> localVec)
68 {
69 if (DDMASTER(dd))
70 {
71 std::vector<gmx::RVec> buffer;
72
73 for (int rank = 0; rank < dd->nnodes; rank++)
74 {
75 if (rank != dd->rank)
76 {
77 const auto &domainGroups = dd->ma->domainGroups[rank];
78
79 buffer.resize(domainGroups.numAtoms);
80
81 int localAtom = 0;
82 for (const int &cg : domainGroups.atomGroups)
83 {
84 for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
85 {
86 buffer[localAtom++] = globalVec[globalAtom];
87 }
88 }
89 GMX_RELEASE_ASSERT(localAtom == domainGroups.numAtoms, "The index count and number of indices should match");
90
91 #if GMX_MPI
92 MPI_Send(buffer.data(), domainGroups.numAtoms*sizeof(gmx::RVec), MPI_BYTE,
93 rank, rank, dd->mpi_comm_all);
94 #endif
95 }
96 }
97
98 const auto &domainGroups = dd->ma->domainGroups[dd->masterrank];
99 int localAtom = 0;
100 for (const int &cg : domainGroups.atomGroups)
101 {
102 for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
103 {
104 localVec[localAtom++] = globalVec[globalAtom];
105 }
106 }
107 }
108 else
109 {
110 #if GMX_MPI
111 MPI_Recv(localVec.data(), dd->nat_home*sizeof(gmx::RVec), MPI_BYTE, dd->masterrank,
112 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
113 #endif
114 }
115 }
116
117 static void distributeVecScatterv(gmx_domdec_t *dd,
118 const t_block *cgs,
119 gmx::ArrayRef<const gmx::RVec> globalVec,
120 gmx::ArrayRef<gmx::RVec> localVec)
121 {
122 int *sendCounts = nullptr;
123 int *displacements = nullptr;
124
125 if (DDMASTER(dd))
126 {
127 AtomDistribution &ma = *dd->ma.get();
128
129 get_commbuffer_counts(&ma, &sendCounts, &displacements);
130
131 gmx::ArrayRef<gmx::RVec> buffer = ma.rvecBuffer;
132 int localAtom = 0;
133 for (int rank = 0; rank < dd->nnodes; rank++)
134 {
135 const auto &domainGroups = ma.domainGroups[rank];
136 for (const int &cg : domainGroups.atomGroups)
137 {
138 for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
139 {
140 buffer[localAtom++] = globalVec[globalAtom];
141 }
142 }
143 }
144 }
145
146 dd_scatterv(dd, sendCounts, displacements,
147 DDMASTER(dd) ? dd->ma->rvecBuffer.data() : nullptr,
148 dd->nat_home*sizeof(gmx::RVec), localVec.data());
149 }
150
151 static void distributeVec(gmx_domdec_t *dd,
152 const t_block *cgs,
153 gmx::ArrayRef<const gmx::RVec> globalVec,
154 gmx::ArrayRef<gmx::RVec> localVec)
155 {
156 if (dd->nnodes <= c_maxNumRanksUseSendRecvForScatterAndGather)
157 {
158 distributeVecSendrecv(dd, cgs, globalVec, localVec);
159 }
160 else
161 {
162 distributeVecScatterv(dd, cgs, globalVec, localVec);
163 }
164 }
165
166 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
167 {
168 if (dfhist == nullptr)
169 {
170 return;
171 }
172
173 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
174 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
175 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
176
177 if (dfhist->nlambda > 0)
178 {
179 int nlam = dfhist->nlambda;
180 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
181 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
182 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
183 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
184 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
185 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
186
187 for (int i = 0; i < nlam; i++)
188 {
189 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
190 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
191 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
192 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
193 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
194 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
195 }
196 }
197 }
198
199 static void dd_distribute_state(gmx_domdec_t *dd, const t_block *cgs,
200 const t_state *state, t_state *state_local,
201 PaddedRVecVector *f)
202 {
203 int nh = state_local->nhchainlength;
204
205 if (DDMASTER(dd))
206 {
207 GMX_RELEASE_ASSERT(state->nhchainlength == nh, "The global and local Nose-Hoover chain lengths should match");
208
209 for (int i = 0; i < efptNR; i++)
210 {
211 state_local->lambda[i] = state->lambda[i];
212 }
213 state_local->fep_state = state->fep_state;
214 state_local->veta = state->veta;
215 state_local->vol0 = state->vol0;
216 copy_mat(state->box, state_local->box);
217 copy_mat(state->box_rel, state_local->box_rel);
218 copy_mat(state->boxv, state_local->boxv);
219 copy_mat(state->svir_prev, state_local->svir_prev);
220 copy_mat(state->fvir_prev, state_local->fvir_prev);
221 if (state->dfhist != nullptr)
222 {
223 copy_df_history(state_local->dfhist, state->dfhist);
224 }
225 for (int i = 0; i < state_local->ngtc; i++)
226 {
227 for (int j = 0; j < nh; j++)
228 {
229 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
230 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
231 }
232 state_local->therm_integral[i] = state->therm_integral[i];
233 }
234 for (int i = 0; i < state_local->nnhpres; i++)
235 {
236 for (int j = 0; j < nh; j++)
237 {
238 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
239 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
240 }
241 }
242 state_local->baros_integral = state->baros_integral;
243 }
244 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda.data());
245 dd_bcast(dd, sizeof(int), &state_local->fep_state);
246 dd_bcast(dd, sizeof(real), &state_local->veta);
247 dd_bcast(dd, sizeof(real), &state_local->vol0);
248 dd_bcast(dd, sizeof(state_local->box), state_local->box);
249 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
250 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
251 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
252 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
253 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi.data());
254 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi.data());
255 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral.data());
256 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi.data());
257 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi.data());
258
259 /* communicate df_history -- required for restarting from checkpoint */
260 dd_distribute_dfhist(dd, state_local->dfhist);
261
262 dd_resize_state(state_local, f, dd->nat_home);
263
264 if (state_local->flags & (1 << estX))
265 {
266 distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->x) : gmx::EmptyArrayRef(), state_local->x);
267 }
268 if (state_local->flags & (1 << estV))
269 {
270 distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->v) : gmx::EmptyArrayRef(), state_local->v);
271 }
272 if (state_local->flags & (1 << estCGP))
273 {
274 distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->cg_p) : gmx::EmptyArrayRef(), state_local->cg_p);
275 }
276 }
277
278 static std::vector < std::vector < int>>
279 getAtomGroupDistribution(FILE *fplog,
280 const matrix box, const gmx_ddbox_t &ddbox,
281 const t_block *cgs, rvec pos[],
282 gmx_domdec_t *dd)
283 {
284 AtomDistribution &ma = *dd->ma.get();
285
286 /* Clear the count */
287 for (int rank = 0; rank < dd->nnodes; rank++)
288 {
289 ma.domainGroups[rank].numAtoms = 0;
290 }
291
292 matrix tcm;
293 make_tric_corr_matrix(dd->npbcdim, box, tcm);
294
295 ivec npulse;
296 const auto cellBoundaries =
297 set_dd_cell_sizes_slb(dd, &ddbox, setcellsizeslbMASTER, npulse);
298
299 const int *cgindex = cgs->index;
300
301 std::vector < std::vector < int>> indices(dd->nnodes);
302
303 /* Compute the center of geometry for all charge groups */
304 for (int icg = 0; icg < cgs->nr; icg++)
305 {
306 /* Set the reference location cg_cm for assigning the group */
307 rvec cg_cm;
308 int k0 = cgindex[icg];
309 int k1 = cgindex[icg + 1];
310 int nrcg = k1 - k0;
311 if (nrcg == 1)
312 {
313 copy_rvec(pos[k0], cg_cm);
314 }
315 else
316 {
317 real inv_ncg = 1/static_cast<real>(nrcg);
318
319 clear_rvec(cg_cm);
320 for (int k = k0; k < k1; k++)
321 {
322 rvec_inc(cg_cm, pos[k]);
323 }
324 for (int d = 0; d < DIM; d++)
325 {
326 cg_cm[d] *= inv_ncg;
327 }
328 }
329 /* Put the charge group in the box and determine the cell index ind */
330 ivec ind;
331 for (int d = DIM - 1; d >= 0; d--)
332 {
333 real pos_d = cg_cm[d];
334 if (d < dd->npbcdim)
335 {
336 bool bScrew = (dd->bScrewPBC && d == XX);
337 if (ddbox.tric_dir[d] && dd->nc[d] > 1)
338 {
339 /* Use triclinic coordinates for this dimension */
340 for (int j = d + 1; j < DIM; j++)
341 {
342 pos_d += cg_cm[j]*tcm[j][d];
343 }
344 }
345 while (pos_d >= box[d][d])
346 {
347 pos_d -= box[d][d];
348 rvec_dec(cg_cm, box[d]);
349 if (bScrew)
350 {
351 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
352 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
353 }
354 for (int k = k0; k < k1; k++)
355 {
356 rvec_dec(pos[k], box[d]);
357 if (bScrew)
358 {
359 pos[k][YY] = box[YY][YY] - pos[k][YY];
360 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
361 }
362 }
363 }
364 while (pos_d < 0)
365 {
366 pos_d += box[d][d];
367 rvec_inc(cg_cm, box[d]);
368 if (bScrew)
369 {
370 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
371 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
372 }
373 for (int k = k0; k < k1; k++)
374 {
375 rvec_inc(pos[k], box[d]);
376 if (bScrew)
377 {
378 pos[k][YY] = box[YY][YY] - pos[k][YY];
379 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
380 }
381 }
382 }
383 }
384 /* This could be done more efficiently */
385 ind[d] = 0;
386 while (ind[d]+1 < dd->nc[d] && pos_d >= cellBoundaries[d][ind[d] + 1])
387 {
388 ind[d]++;
389 }
390 }
391 int domainIndex = dd_index(dd->nc, ind);
392 indices[domainIndex].push_back(icg);
393 ma.domainGroups[domainIndex].numAtoms += cgindex[icg + 1] - cgindex[icg];
394 }
395
396 if (fplog)
397 {
398 // Use double for the sums to avoid natoms^2 overflowing
399 // (65537^2 > 2^32)
400 int nat_sum, nat_min, nat_max;
401 double nat2_sum;
402
403 nat_sum = 0;
404 nat2_sum = 0;
405 nat_min = ma.domainGroups[0].numAtoms;
406 nat_max = ma.domainGroups[0].numAtoms;
407 for (int rank = 0; rank < dd->nnodes; rank++)
408 {
409 int numAtoms = ma.domainGroups[rank].numAtoms;
410 nat_sum += numAtoms;
411 // cast to double to avoid integer overflows when squaring
412 nat2_sum += gmx::square(static_cast<double>(numAtoms));
413 nat_min = std::min(nat_min, numAtoms);
414 nat_max = std::max(nat_max, numAtoms);
415 }
416 nat_sum /= dd->nnodes;
417 nat2_sum /= dd->nnodes;
418
419 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
420 dd->nnodes,
421 nat_sum,
422 static_cast<int>(std::sqrt(nat2_sum - gmx::square(static_cast<double>(nat_sum)) + 0.5)),
423 nat_min, nat_max);
424 }
425
426 return indices;
427 }
428
429 static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
430 const t_block *cgs,
431 const matrix box, const gmx_ddbox_t *ddbox,
432 rvec pos[])
433 {
434 AtomDistribution *ma = dd->ma.get();
435 int *ibuf, buf2[2] = { 0, 0 };
436 gmx_bool bMaster = DDMASTER(dd);
437
438 std::vector < std::vector < int>> groupIndices;
439
440 if (bMaster)
441 {
442 if (dd->bScrewPBC)
443 {
444 check_screw_box(box);
445 }
446
447 groupIndices = getAtomGroupDistribution(fplog, box, *ddbox, cgs, pos, dd);
448
449 for (int rank = 0; rank < dd->nnodes; rank++)
450 {
451 ma->intBuffer[rank*2] = groupIndices[rank].size();
452 ma->intBuffer[rank*2 + 1] = ma->domainGroups[rank].numAtoms;
453 }
454 ibuf = ma->intBuffer.data();
455 }
456 else
457 {
458 ibuf = nullptr;
459 }
460 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
461
462 dd->ncg_home = buf2[0];
463 dd->nat_home = buf2[1];
464 dd->globalAtomGroupIndices.resize(dd->ncg_home);
465 dd->globalAtomIndices.resize(dd->nat_home);
466
467 if (bMaster)
468 {
469 ma->atomGroups.clear();
470
471 int groupOffset = 0;
472 for (int rank = 0; rank < dd->nnodes; rank++)
473 {
474 ma->intBuffer[rank] = groupIndices[rank].size()*sizeof(int);
475 ma->intBuffer[dd->nnodes + rank] = groupOffset*sizeof(int);
476
477 ma->atomGroups.insert(ma->atomGroups.end(),
478 groupIndices[rank].begin(), groupIndices[rank].end());
479
480 ma->domainGroups[rank].atomGroups = gmx::constArrayRefFromArray(ma->atomGroups.data() + groupOffset, groupIndices[rank].size());
481
482 groupOffset += groupIndices[rank].size();
483 }
484 }
485
486 dd_scatterv(dd,
487 bMaster ? ma->intBuffer.data() : nullptr,
488 bMaster ? ma->intBuffer.data() + dd->nnodes : nullptr,
489 bMaster ? ma->atomGroups.data() : nullptr,
490 dd->ncg_home*sizeof(int), dd->globalAtomGroupIndices.data());
491
492 /* Determine the home charge group sizes */
493 dd->atomGroups_.index.resize(dd->ncg_home + 1);
494 dd->atomGroups_.index[0] = 0;
495 for (int i = 0; i < dd->ncg_home; i++)
496 {
497 int cg_gl = dd->globalAtomGroupIndices[i];
498 dd->atomGroups_.index[i+1] =
499 dd->atomGroups_.index[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
500 }
501
502 if (debug)
503 {
504 fprintf(debug, "Home charge groups:\n");
505 for (int i = 0; i < dd->ncg_home; i++)
506 {
507 fprintf(debug, " %d", dd->globalAtomGroupIndices[i]);
508 if (i % 10 == 9)
509 {
510 fprintf(debug, "\n");
511 }
512 }
513 fprintf(debug, "\n");
514 }
515 }
516
517 void distributeState(FILE *fplog,
518 gmx_domdec_t *dd,
519 t_state *state_global,
520 const t_block &cgs_gl,
521 const gmx_ddbox_t &ddbox,
522 t_state *state_local,
523 PaddedRVecVector *f)
524 {
525 rvec *xGlobal = (DDMASTER(dd) ? as_rvec_array(state_global->x.data()) : nullptr);
526
527 distributeAtomGroups(fplog, dd, &cgs_gl,
528 DDMASTER(dd) ? state_global->box : nullptr,
529 &ddbox, xGlobal);
530
531 dd_distribute_state(dd, &cgs_gl,
532 state_global, state_local, f);
533 }
The ultimate molecular dynamics simulation package