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, 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.
51 #include "gromacs/domdec/domdec_network.h"
52 #include "gromacs/domdec/ga2la.h"
53 #include "gromacs/ewald/pme.h"
54 #include "gromacs/fileio/gmxfio.h"
55 #include "gromacs/fileio/pdbio.h"
56 #include "gromacs/gmxlib/chargegroup.h"
57 #include "gromacs/gmxlib/network.h"
58 #include "gromacs/gmxlib/nrnb.h"
59 #include "gromacs/gpu_utils/gpu_utils.h"
60 #include "gromacs/hardware/hw_info.h"
61 #include "gromacs/imd/imd.h"
62 #include "gromacs/listed-forces/manage-threading.h"
63 #include "gromacs/math/functions.h"
64 #include "gromacs/math/vec.h"
65 #include "gromacs/math/vectypes.h"
66 #include "gromacs/mdlib/constr.h"
67 #include "gromacs/mdlib/force.h"
68 #include "gromacs/mdlib/forcerec.h"
69 #include "gromacs/mdlib/genborn.h"
70 #include "gromacs/mdlib/gmx_omp_nthreads.h"
71 #include "gromacs/mdlib/mdatoms.h"
72 #include "gromacs/mdlib/mdrun.h"
73 #include "gromacs/mdlib/mdsetup.h"
74 #include "gromacs/mdlib/nb_verlet.h"
75 #include "gromacs/mdlib/nbnxn_grid.h"
76 #include "gromacs/mdlib/nsgrid.h"
77 #include "gromacs/mdlib/vsite.h"
78 #include "gromacs/mdtypes/commrec.h"
79 #include "gromacs/mdtypes/df_history.h"
80 #include "gromacs/mdtypes/forcerec.h"
81 #include "gromacs/mdtypes/inputrec.h"
82 #include "gromacs/mdtypes/md_enums.h"
83 #include "gromacs/mdtypes/mdatom.h"
84 #include "gromacs/mdtypes/nblist.h"
85 #include "gromacs/mdtypes/state.h"
86 #include "gromacs/pbcutil/ishift.h"
87 #include "gromacs/pbcutil/pbc.h"
88 #include "gromacs/pulling/pull.h"
89 #include "gromacs/pulling/pull_rotation.h"
90 #include "gromacs/swap/swapcoords.h"
91 #include "gromacs/timing/wallcycle.h"
92 #include "gromacs/topology/block.h"
93 #include "gromacs/topology/idef.h"
94 #include "gromacs/topology/ifunc.h"
95 #include "gromacs/topology/mtop_util.h"
96 #include "gromacs/topology/topology.h"
97 #include "gromacs/utility/basedefinitions.h"
98 #include "gromacs/utility/basenetwork.h"
99 #include "gromacs/utility/cstringutil.h"
100 #include "gromacs/utility/exceptions.h"
101 #include "gromacs/utility/fatalerror.h"
102 #include "gromacs/utility/gmxmpi.h"
103 #include "gromacs/utility/qsort_threadsafe.h"
104 #include "gromacs/utility/real.h"
105 #include "gromacs/utility/smalloc.h"
107 #include "domdec_constraints.h"
108 #include "domdec_internal.h"
109 #include "domdec_vsite.h"
111 #define DDRANK(dd, rank) (rank)
112 #define DDMASTERRANK(dd) (dd->masterrank)
114 struct gmx_domdec_master_t
116 /* The cell boundaries */
118 /* The global charge group division */
119 int *ncg; /* Number of home charge groups for each node */
120 int *index; /* Index of nnodes+1 into cg */
121 int *cg; /* Global charge group index */
122 int *nat; /* Number of home atoms for each node. */
123 int *ibuf; /* Buffer for communication */
124 rvec *vbuf; /* Buffer for state scattering and gathering */
127 #define DD_NLOAD_MAX 9
129 const char *edlbs_names[edlbsNR] = { "off", "auto", "locked", "on" };
131 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
134 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
135 #define DD_FLAG_NRCG 65535
136 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
137 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
139 /* The DD zone order */
140 static const ivec dd_zo[DD_MAXZONE] =
141 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
143 /* The non-bonded zone-pair setup for domain decomposition
144 * The first number is the i-zone, the second number the first j-zone seen by
145 * this i-zone, the third number the last+1 j-zone seen by this i-zone.
146 * As is, this is for 3D decomposition, where there are 4 i-zones.
147 * With 2D decomposition use only the first 2 i-zones and a last+1 j-zone of 4.
148 * With 1D decomposition use only the first i-zone and a last+1 j-zone of 2.
151 ddNonbondedZonePairRanges[DD_MAXIZONE][3] = {{0, 0, 8},
156 /* Factors used to avoid problems due to rounding issues */
157 #define DD_CELL_MARGIN 1.0001
158 #define DD_CELL_MARGIN2 1.00005
159 /* Factor to account for pressure scaling during nstlist steps */
160 #define DD_PRES_SCALE_MARGIN 1.02
162 /* Turn on DLB when the load imbalance causes this amount of total loss.
163 * There is a bit of overhead with DLB and it's difficult to achieve
164 * a load imbalance of less than 2% with DLB.
166 #define DD_PERF_LOSS_DLB_ON 0.02
168 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
169 #define DD_PERF_LOSS_WARN 0.05
171 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
173 /* Use separate MPI send and receive commands
174 * when nnodes <= GMX_DD_NNODES_SENDRECV.
175 * This saves memory (and some copying for small nnodes).
176 * For high parallelization scatter and gather calls are used.
178 #define GMX_DD_NNODES_SENDRECV 4
182 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
184 static void index2xyz(ivec nc,int ind,ivec xyz)
186 xyz[XX] = ind % nc[XX];
187 xyz[YY] = (ind / nc[XX]) % nc[YY];
188 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
192 /* This order is required to minimize the coordinate communication in PME
193 * which uses decomposition in the x direction.
195 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
197 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
199 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
200 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
201 xyz[ZZ] = ind % nc[ZZ];
204 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
209 ddindex = dd_index(dd->nc, c);
210 if (dd->comm->bCartesianPP_PME)
212 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
214 else if (dd->comm->bCartesianPP)
217 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
228 static gmx_bool dynamic_dd_box(const gmx_ddbox_t *ddbox, const t_inputrec *ir)
230 return (ddbox->nboundeddim < DIM || inputrecDynamicBox(ir));
233 int ddglatnr(const gmx_domdec_t *dd, int i)
243 if (i >= dd->comm->nat[ddnatNR-1])
245 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
247 atnr = dd->gatindex[i] + 1;
253 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
255 return &dd->comm->cgs_gl;
258 static bool dlbIsOn(const gmx_domdec_comm_t *comm)
260 return (comm->dlbState == edlbsOn);
263 static void vec_rvec_init(vec_rvec_t *v)
269 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
273 v->nalloc = over_alloc_dd(n);
274 srenew(v->v, v->nalloc);
278 void dd_store_state(gmx_domdec_t *dd, t_state *state)
282 if (state->ddp_count != dd->ddp_count)
284 gmx_incons("The state does not the domain decomposition state");
287 state->cg_gl.resize(dd->ncg_home);
288 for (i = 0; i < dd->ncg_home; i++)
290 state->cg_gl[i] = dd->index_gl[i];
293 state->ddp_count_cg_gl = dd->ddp_count;
296 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
298 return &dd->comm->zones;
301 void dd_get_ns_ranges(const gmx_domdec_t *dd, int icg,
302 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
304 gmx_domdec_zones_t *zones;
307 zones = &dd->comm->zones;
310 while (icg >= zones->izone[izone].cg1)
319 else if (izone < zones->nizone)
321 *jcg0 = zones->izone[izone].jcg0;
325 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
326 icg, izone, zones->nizone);
329 *jcg1 = zones->izone[izone].jcg1;
331 for (d = 0; d < dd->ndim; d++)
334 shift0[dim] = zones->izone[izone].shift0[dim];
335 shift1[dim] = zones->izone[izone].shift1[dim];
336 if (dd->comm->tric_dir[dim] || (dlbIsOn(dd->comm) && d > 0))
338 /* A conservative approach, this can be optimized */
345 int dd_natoms_mdatoms(const gmx_domdec_t *dd)
347 /* We currently set mdatoms entries for all atoms:
348 * local + non-local + communicated for vsite + constraints
351 return dd->comm->nat[ddnatNR - 1];
354 int dd_natoms_vsite(const gmx_domdec_t *dd)
356 return dd->comm->nat[ddnatVSITE];
359 void dd_get_constraint_range(const gmx_domdec_t *dd, int *at_start, int *at_end)
361 *at_start = dd->comm->nat[ddnatCON-1];
362 *at_end = dd->comm->nat[ddnatCON];
365 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
367 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
368 int *index, *cgindex;
369 gmx_domdec_comm_t *comm;
370 gmx_domdec_comm_dim_t *cd;
371 gmx_domdec_ind_t *ind;
372 rvec shift = {0, 0, 0}, *buf, *rbuf;
373 gmx_bool bPBC, bScrew;
377 cgindex = dd->cgindex;
382 nat_tot = dd->nat_home;
383 for (d = 0; d < dd->ndim; d++)
385 bPBC = (dd->ci[dd->dim[d]] == 0);
386 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
389 copy_rvec(box[dd->dim[d]], shift);
392 for (p = 0; p < cd->np; p++)
399 for (i = 0; i < ind->nsend[nzone]; i++)
401 at0 = cgindex[index[i]];
402 at1 = cgindex[index[i]+1];
403 for (j = at0; j < at1; j++)
405 copy_rvec(x[j], buf[n]);
412 for (i = 0; i < ind->nsend[nzone]; i++)
414 at0 = cgindex[index[i]];
415 at1 = cgindex[index[i]+1];
416 for (j = at0; j < at1; j++)
418 /* We need to shift the coordinates */
419 rvec_add(x[j], shift, buf[n]);
426 for (i = 0; i < ind->nsend[nzone]; i++)
428 at0 = cgindex[index[i]];
429 at1 = cgindex[index[i]+1];
430 for (j = at0; j < at1; j++)
433 buf[n][XX] = x[j][XX] + shift[XX];
435 * This operation requires a special shift force
436 * treatment, which is performed in calc_vir.
438 buf[n][YY] = box[YY][YY] - x[j][YY];
439 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
451 rbuf = comm->vbuf2.v;
453 /* Send and receive the coordinates */
454 dd_sendrecv_rvec(dd, d, dddirBackward,
455 buf, ind->nsend[nzone+1],
456 rbuf, ind->nrecv[nzone+1]);
460 for (zone = 0; zone < nzone; zone++)
462 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
464 copy_rvec(rbuf[j], x[i]);
469 nat_tot += ind->nrecv[nzone+1];
475 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
477 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
478 int *index, *cgindex;
479 gmx_domdec_comm_t *comm;
480 gmx_domdec_comm_dim_t *cd;
481 gmx_domdec_ind_t *ind;
485 gmx_bool bShiftForcesNeedPbc, bScrew;
489 cgindex = dd->cgindex;
493 nzone = comm->zones.n/2;
494 nat_tot = dd->nat_tot;
495 for (d = dd->ndim-1; d >= 0; d--)
497 /* Only forces in domains near the PBC boundaries need to
498 consider PBC in the treatment of fshift */
499 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
500 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
501 if (fshift == NULL && !bScrew)
503 bShiftForcesNeedPbc = FALSE;
505 /* Determine which shift vector we need */
511 for (p = cd->np-1; p >= 0; p--)
514 nat_tot -= ind->nrecv[nzone+1];
521 sbuf = comm->vbuf2.v;
523 for (zone = 0; zone < nzone; zone++)
525 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
527 copy_rvec(f[i], sbuf[j]);
532 /* Communicate the forces */
533 dd_sendrecv_rvec(dd, d, dddirForward,
534 sbuf, ind->nrecv[nzone+1],
535 buf, ind->nsend[nzone+1]);
537 /* Add the received forces */
539 if (!bShiftForcesNeedPbc)
541 for (i = 0; i < ind->nsend[nzone]; i++)
543 at0 = cgindex[index[i]];
544 at1 = cgindex[index[i]+1];
545 for (j = at0; j < at1; j++)
547 rvec_inc(f[j], buf[n]);
554 /* fshift should always be defined if this function is
555 * called when bShiftForcesNeedPbc is true */
556 assert(NULL != fshift);
557 for (i = 0; i < ind->nsend[nzone]; i++)
559 at0 = cgindex[index[i]];
560 at1 = cgindex[index[i]+1];
561 for (j = at0; j < at1; j++)
563 rvec_inc(f[j], buf[n]);
564 /* Add this force to the shift force */
565 rvec_inc(fshift[is], buf[n]);
572 for (i = 0; i < ind->nsend[nzone]; i++)
574 at0 = cgindex[index[i]];
575 at1 = cgindex[index[i]+1];
576 for (j = at0; j < at1; j++)
578 /* Rotate the force */
579 f[j][XX] += buf[n][XX];
580 f[j][YY] -= buf[n][YY];
581 f[j][ZZ] -= buf[n][ZZ];
584 /* Add this force to the shift force */
585 rvec_inc(fshift[is], buf[n]);
596 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
598 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
599 int *index, *cgindex;
600 gmx_domdec_comm_t *comm;
601 gmx_domdec_comm_dim_t *cd;
602 gmx_domdec_ind_t *ind;
607 cgindex = dd->cgindex;
609 buf = &comm->vbuf.v[0][0];
612 nat_tot = dd->nat_home;
613 for (d = 0; d < dd->ndim; d++)
616 for (p = 0; p < cd->np; p++)
621 for (i = 0; i < ind->nsend[nzone]; i++)
623 at0 = cgindex[index[i]];
624 at1 = cgindex[index[i]+1];
625 for (j = at0; j < at1; j++)
638 rbuf = &comm->vbuf2.v[0][0];
640 /* Send and receive the coordinates */
641 dd_sendrecv_real(dd, d, dddirBackward,
642 buf, ind->nsend[nzone+1],
643 rbuf, ind->nrecv[nzone+1]);
647 for (zone = 0; zone < nzone; zone++)
649 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
656 nat_tot += ind->nrecv[nzone+1];
662 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
664 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
665 int *index, *cgindex;
666 gmx_domdec_comm_t *comm;
667 gmx_domdec_comm_dim_t *cd;
668 gmx_domdec_ind_t *ind;
673 cgindex = dd->cgindex;
675 buf = &comm->vbuf.v[0][0];
677 nzone = comm->zones.n/2;
678 nat_tot = dd->nat_tot;
679 for (d = dd->ndim-1; d >= 0; d--)
682 for (p = cd->np-1; p >= 0; p--)
685 nat_tot -= ind->nrecv[nzone+1];
692 sbuf = &comm->vbuf2.v[0][0];
694 for (zone = 0; zone < nzone; zone++)
696 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
703 /* Communicate the forces */
704 dd_sendrecv_real(dd, d, dddirForward,
705 sbuf, ind->nrecv[nzone+1],
706 buf, ind->nsend[nzone+1]);
708 /* Add the received forces */
710 for (i = 0; i < ind->nsend[nzone]; i++)
712 at0 = cgindex[index[i]];
713 at1 = cgindex[index[i]+1];
714 for (j = at0; j < at1; j++)
725 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
727 fprintf(fp, "zone d0 %d d1 %d d2 %d min0 %6.3f max1 %6.3f mch0 %6.3f mch1 %6.3f p1_0 %6.3f p1_1 %6.3f\n",
729 zone->min0, zone->max1,
730 zone->mch0, zone->mch0,
731 zone->p1_0, zone->p1_1);
735 #define DDZONECOMM_MAXZONE 5
736 #define DDZONECOMM_BUFSIZE 3
738 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
739 int ddimind, int direction,
740 gmx_ddzone_t *buf_s, int n_s,
741 gmx_ddzone_t *buf_r, int n_r)
743 #define ZBS DDZONECOMM_BUFSIZE
744 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
745 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
748 for (i = 0; i < n_s; i++)
750 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
751 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
752 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
753 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
754 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
755 vbuf_s[i*ZBS+1][2] = 0;
756 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
757 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
758 vbuf_s[i*ZBS+2][2] = 0;
761 dd_sendrecv_rvec(dd, ddimind, direction,
765 for (i = 0; i < n_r; i++)
767 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
768 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
769 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
770 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
771 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
772 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
773 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
779 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
780 rvec cell_ns_x0, rvec cell_ns_x1)
782 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
784 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
785 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
786 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
787 rvec extr_s[2], extr_r[2];
789 real dist_d, c = 0, det;
790 gmx_domdec_comm_t *comm;
795 for (d = 1; d < dd->ndim; d++)
798 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
799 zp->min0 = cell_ns_x0[dim];
800 zp->max1 = cell_ns_x1[dim];
801 zp->min1 = cell_ns_x1[dim];
802 zp->mch0 = cell_ns_x0[dim];
803 zp->mch1 = cell_ns_x1[dim];
804 zp->p1_0 = cell_ns_x0[dim];
805 zp->p1_1 = cell_ns_x1[dim];
808 for (d = dd->ndim-2; d >= 0; d--)
811 bPBC = (dim < ddbox->npbcdim);
813 /* Use an rvec to store two reals */
814 extr_s[d][0] = comm->cell_f0[d+1];
815 extr_s[d][1] = comm->cell_f1[d+1];
816 extr_s[d][2] = comm->cell_f1[d+1];
819 /* Store the extremes in the backward sending buffer,
820 * so the get updated separately from the forward communication.
822 for (d1 = d; d1 < dd->ndim-1; d1++)
824 /* We invert the order to be able to use the same loop for buf_e */
825 buf_s[pos].min0 = extr_s[d1][1];
826 buf_s[pos].max1 = extr_s[d1][0];
827 buf_s[pos].min1 = extr_s[d1][2];
830 /* Store the cell corner of the dimension we communicate along */
831 buf_s[pos].p1_0 = comm->cell_x0[dim];
836 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
839 if (dd->ndim == 3 && d == 0)
841 buf_s[pos] = comm->zone_d2[0][1];
843 buf_s[pos] = comm->zone_d1[0];
847 /* We only need to communicate the extremes
848 * in the forward direction
850 npulse = comm->cd[d].np;
853 /* Take the minimum to avoid double communication */
854 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
858 /* Without PBC we should really not communicate over
859 * the boundaries, but implementing that complicates
860 * the communication setup and therefore we simply
861 * do all communication, but ignore some data.
865 for (p = 0; p < npulse_min; p++)
867 /* Communicate the extremes forward */
868 bUse = (bPBC || dd->ci[dim] > 0);
870 dd_sendrecv_rvec(dd, d, dddirForward,
871 extr_s+d, dd->ndim-d-1,
872 extr_r+d, dd->ndim-d-1);
876 for (d1 = d; d1 < dd->ndim-1; d1++)
878 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
879 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
880 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
886 for (p = 0; p < npulse; p++)
888 /* Communicate all the zone information backward */
889 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
891 dd_sendrecv_ddzone(dd, d, dddirBackward,
898 for (d1 = d+1; d1 < dd->ndim; d1++)
900 /* Determine the decrease of maximum required
901 * communication height along d1 due to the distance along d,
902 * this avoids a lot of useless atom communication.
904 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
906 if (ddbox->tric_dir[dim])
908 /* c is the off-diagonal coupling between the cell planes
909 * along directions d and d1.
911 c = ddbox->v[dim][dd->dim[d1]][dim];
917 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
920 dh[d1] = comm->cutoff - (c*dist_d + std::sqrt(det))/(1 + c*c);
924 /* A negative value signals out of range */
930 /* Accumulate the extremes over all pulses */
931 for (i = 0; i < buf_size; i++)
941 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
942 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
943 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
946 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
954 if (bUse && dh[d1] >= 0)
956 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
957 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
960 /* Copy the received buffer to the send buffer,
961 * to pass the data through with the next pulse.
965 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
966 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
968 /* Store the extremes */
971 for (d1 = d; d1 < dd->ndim-1; d1++)
973 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
974 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
975 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
979 if (d == 1 || (d == 0 && dd->ndim == 3))
981 for (i = d; i < 2; i++)
983 comm->zone_d2[1-d][i] = buf_e[pos];
989 comm->zone_d1[1] = buf_e[pos];
999 for (i = 0; i < 2; i++)
1003 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1005 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1006 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1012 for (i = 0; i < 2; i++)
1014 for (j = 0; j < 2; j++)
1018 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1020 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1021 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1025 for (d = 1; d < dd->ndim; d++)
1027 comm->cell_f_max0[d] = extr_s[d-1][0];
1028 comm->cell_f_min1[d] = extr_s[d-1][1];
1031 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1032 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1037 static void dd_collect_cg(gmx_domdec_t *dd,
1038 t_state *state_local)
1040 gmx_domdec_master_t *ma = NULL;
1041 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1043 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1045 /* The master has the correct distribution */
1049 if (state_local->ddp_count == dd->ddp_count)
1051 /* The local state and DD are in sync, use the DD indices */
1052 ncg_home = dd->ncg_home;
1054 nat_home = dd->nat_home;
1056 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1058 /* The DD is out of sync with the local state, but we have stored
1059 * the cg indices with the local state, so we can use those.
1063 cgs_gl = &dd->comm->cgs_gl;
1065 ncg_home = state_local->cg_gl.size();
1066 cg = state_local->cg_gl.data();
1068 for (i = 0; i < ncg_home; i++)
1070 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1075 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1089 /* Collect the charge group and atom counts on the master */
1090 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1095 for (i = 0; i < dd->nnodes; i++)
1097 ma->ncg[i] = ma->ibuf[2*i];
1098 ma->nat[i] = ma->ibuf[2*i+1];
1099 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1102 /* Make byte counts and indices */
1103 for (i = 0; i < dd->nnodes; i++)
1105 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1106 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1110 fprintf(debug, "Initial charge group distribution: ");
1111 for (i = 0; i < dd->nnodes; i++)
1113 fprintf(debug, " %d", ma->ncg[i]);
1115 fprintf(debug, "\n");
1119 /* Collect the charge group indices on the master */
1121 ncg_home*sizeof(int), cg,
1122 DDMASTER(dd) ? ma->ibuf : NULL,
1123 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1124 DDMASTER(dd) ? ma->cg : NULL);
1126 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1129 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1130 const rvec *lv, rvec *v)
1132 gmx_domdec_master_t *ma;
1133 int n, i, c, a, nalloc = 0;
1142 MPI_Send(const_cast<void *>(static_cast<const void *>(lv)), dd->nat_home*sizeof(rvec), MPI_BYTE,
1143 DDMASTERRANK(dd), dd->rank, dd->mpi_comm_all);
1148 /* Copy the master coordinates to the global array */
1149 cgs_gl = &dd->comm->cgs_gl;
1151 n = DDMASTERRANK(dd);
1153 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1155 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1157 copy_rvec(lv[a++], v[c]);
1161 for (n = 0; n < dd->nnodes; n++)
1165 if (ma->nat[n] > nalloc)
1167 nalloc = over_alloc_dd(ma->nat[n]);
1168 srenew(buf, nalloc);
1171 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1172 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1175 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1177 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1179 copy_rvec(buf[a++], v[c]);
1188 static void get_commbuffer_counts(gmx_domdec_t *dd,
1189 int **counts, int **disps)
1191 gmx_domdec_master_t *ma;
1196 /* Make the rvec count and displacment arrays */
1198 *disps = ma->ibuf + dd->nnodes;
1199 for (n = 0; n < dd->nnodes; n++)
1201 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1202 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1206 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1207 const rvec *lv, rvec *v)
1209 gmx_domdec_master_t *ma;
1210 int *rcounts = NULL, *disps = NULL;
1219 get_commbuffer_counts(dd, &rcounts, &disps);
1224 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1228 cgs_gl = &dd->comm->cgs_gl;
1231 for (n = 0; n < dd->nnodes; n++)
1233 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1235 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1237 copy_rvec(buf[a++], v[c]);
1244 void dd_collect_vec(gmx_domdec_t *dd,
1245 t_state *state_local,
1246 const PaddedRVecVector *localVector,
1249 dd_collect_cg(dd, state_local);
1251 const rvec *lv = as_rvec_array(localVector->data());
1253 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1255 dd_collect_vec_sendrecv(dd, lv, v);
1259 dd_collect_vec_gatherv(dd, lv, v);
1263 void dd_collect_vec(gmx_domdec_t *dd,
1264 t_state *state_local,
1265 const PaddedRVecVector *localVector,
1266 PaddedRVecVector *vector)
1268 dd_collect_vec(dd, state_local, localVector, as_rvec_array(vector->data()));
1272 void dd_collect_state(gmx_domdec_t *dd,
1273 t_state *state_local, t_state *state)
1277 nh = state->nhchainlength;
1281 for (i = 0; i < efptNR; i++)
1283 state->lambda[i] = state_local->lambda[i];
1285 state->fep_state = state_local->fep_state;
1286 state->veta = state_local->veta;
1287 state->vol0 = state_local->vol0;
1288 copy_mat(state_local->box, state->box);
1289 copy_mat(state_local->boxv, state->boxv);
1290 copy_mat(state_local->svir_prev, state->svir_prev);
1291 copy_mat(state_local->fvir_prev, state->fvir_prev);
1292 copy_mat(state_local->pres_prev, state->pres_prev);
1294 for (i = 0; i < state_local->ngtc; i++)
1296 for (j = 0; j < nh; j++)
1298 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1299 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1301 state->therm_integral[i] = state_local->therm_integral[i];
1303 for (i = 0; i < state_local->nnhpres; i++)
1305 for (j = 0; j < nh; j++)
1307 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1308 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1312 for (est = 0; est < estNR; est++)
1314 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1319 dd_collect_vec(dd, state_local, &state_local->x, &state->x);
1322 dd_collect_vec(dd, state_local, &state_local->v, &state->v);
1324 case est_SDX_NOTSUPPORTED:
1327 dd_collect_vec(dd, state_local, &state_local->cg_p, &state->cg_p);
1329 case estDISRE_INITF:
1330 case estDISRE_RM3TAV:
1331 case estORIRE_INITF:
1335 gmx_incons("Unknown state entry encountered in dd_collect_state");
1341 static void dd_resize_state(t_state *state, PaddedRVecVector *f, int natoms)
1347 fprintf(debug, "Resizing state: currently %d, required %d\n", state->natoms, natoms);
1350 for (est = 0; est < estNR; est++)
1352 if (EST_DISTR(est) && (state->flags & (1<<est)))
1354 /* We need to allocate one element extra, since we might use
1355 * (unaligned) 4-wide SIMD loads to access rvec entries.
1360 state->x.resize(natoms + 1);
1363 state->v.resize(natoms + 1);
1365 case est_SDX_NOTSUPPORTED:
1368 state->cg_p.resize(natoms + 1);
1370 case estDISRE_INITF:
1371 case estDISRE_RM3TAV:
1372 case estORIRE_INITF:
1374 /* No reallocation required */
1377 gmx_incons("Unknown state entry encountered in dd_resize_state");
1384 (*f).resize(natoms + 1);
1388 static void dd_check_alloc_ncg(t_forcerec *fr,
1390 PaddedRVecVector *f,
1391 int numChargeGroups)
1393 if (numChargeGroups > fr->cg_nalloc)
1397 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, numChargeGroups, over_alloc_dd(numChargeGroups));
1399 fr->cg_nalloc = over_alloc_dd(numChargeGroups);
1400 srenew(fr->cginfo, fr->cg_nalloc);
1401 if (fr->cutoff_scheme == ecutsGROUP)
1403 srenew(fr->cg_cm, fr->cg_nalloc);
1406 if (fr->cutoff_scheme == ecutsVERLET)
1408 /* We don't use charge groups, we use x in state to set up
1409 * the atom communication.
1411 dd_resize_state(state, f, numChargeGroups);
1415 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1418 gmx_domdec_master_t *ma;
1419 int n, i, c, a, nalloc = 0;
1426 for (n = 0; n < dd->nnodes; n++)
1430 if (ma->nat[n] > nalloc)
1432 nalloc = over_alloc_dd(ma->nat[n]);
1433 srenew(buf, nalloc);
1435 /* Use lv as a temporary buffer */
1437 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1439 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1441 copy_rvec(v[c], buf[a++]);
1444 if (a != ma->nat[n])
1446 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1451 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1452 DDRANK(dd, n), n, dd->mpi_comm_all);
1457 n = DDMASTERRANK(dd);
1459 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1461 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1463 copy_rvec(v[c], lv[a++]);
1470 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1471 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1476 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1479 gmx_domdec_master_t *ma;
1480 int *scounts = NULL, *disps = NULL;
1488 get_commbuffer_counts(dd, &scounts, &disps);
1492 for (n = 0; n < dd->nnodes; n++)
1494 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1496 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1498 copy_rvec(v[c], buf[a++]);
1504 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1507 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1509 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1511 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1515 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1519 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1526 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1527 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1528 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1530 if (dfhist->nlambda > 0)
1532 int nlam = dfhist->nlambda;
1533 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1534 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1535 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1536 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1537 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1538 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1540 for (int i = 0; i < nlam; i++)
1542 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1543 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1544 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1545 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1546 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1547 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1552 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1553 t_state *state, t_state *state_local,
1554 PaddedRVecVector *f)
1558 nh = state->nhchainlength;
1562 for (i = 0; i < efptNR; i++)
1564 state_local->lambda[i] = state->lambda[i];
1566 state_local->fep_state = state->fep_state;
1567 state_local->veta = state->veta;
1568 state_local->vol0 = state->vol0;
1569 copy_mat(state->box, state_local->box);
1570 copy_mat(state->box_rel, state_local->box_rel);
1571 copy_mat(state->boxv, state_local->boxv);
1572 copy_mat(state->svir_prev, state_local->svir_prev);
1573 copy_mat(state->fvir_prev, state_local->fvir_prev);
1574 if (state->dfhist != NULL)
1576 copy_df_history(state_local->dfhist, state->dfhist);
1578 for (i = 0; i < state_local->ngtc; i++)
1580 for (j = 0; j < nh; j++)
1582 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1583 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1585 state_local->therm_integral[i] = state->therm_integral[i];
1587 for (i = 0; i < state_local->nnhpres; i++)
1589 for (j = 0; j < nh; j++)
1591 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1592 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1596 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda.data());
1597 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1598 dd_bcast(dd, sizeof(real), &state_local->veta);
1599 dd_bcast(dd, sizeof(real), &state_local->vol0);
1600 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1601 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1602 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1603 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1604 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1605 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi.data());
1606 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi.data());
1607 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral.data());
1608 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi.data());
1609 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi.data());
1611 /* communicate df_history -- required for restarting from checkpoint */
1612 dd_distribute_dfhist(dd, state_local->dfhist);
1614 dd_resize_state(state_local, f, dd->nat_home);
1616 for (i = 0; i < estNR; i++)
1618 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1623 dd_distribute_vec(dd, cgs, as_rvec_array(state->x.data()), as_rvec_array(state_local->x.data()));
1626 dd_distribute_vec(dd, cgs, as_rvec_array(state->v.data()), as_rvec_array(state_local->v.data()));
1628 case est_SDX_NOTSUPPORTED:
1631 dd_distribute_vec(dd, cgs, as_rvec_array(state->cg_p.data()), as_rvec_array(state_local->cg_p.data()));
1633 case estDISRE_INITF:
1634 case estDISRE_RM3TAV:
1635 case estORIRE_INITF:
1637 /* Not implemented yet */
1640 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1646 static char dim2char(int dim)
1652 case XX: c = 'X'; break;
1653 case YY: c = 'Y'; break;
1654 case ZZ: c = 'Z'; break;
1655 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1661 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1662 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1664 rvec grid_s[2], *grid_r = NULL, cx, r;
1665 char fname[STRLEN], buf[22];
1667 int a, i, d, z, y, x;
1671 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1672 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1676 snew(grid_r, 2*dd->nnodes);
1679 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1683 for (d = 0; d < DIM; d++)
1685 for (i = 0; i < DIM; i++)
1693 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1695 tric[d][i] = box[i][d]/box[i][i];
1704 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1705 out = gmx_fio_fopen(fname, "w");
1706 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1708 for (i = 0; i < dd->nnodes; i++)
1710 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1711 for (d = 0; d < DIM; d++)
1713 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1715 for (z = 0; z < 2; z++)
1717 for (y = 0; y < 2; y++)
1719 for (x = 0; x < 2; x++)
1721 cx[XX] = grid_r[i*2+x][XX];
1722 cx[YY] = grid_r[i*2+y][YY];
1723 cx[ZZ] = grid_r[i*2+z][ZZ];
1725 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
1726 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
1730 for (d = 0; d < DIM; d++)
1732 for (x = 0; x < 4; x++)
1736 case 0: y = 1 + i*8 + 2*x; break;
1737 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
1738 case 2: y = 1 + i*8 + x; break;
1740 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
1744 gmx_fio_fclose(out);
1749 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
1750 const gmx_mtop_t *mtop, t_commrec *cr,
1751 int natoms, rvec x[], matrix box)
1753 char fname[STRLEN], buf[22];
1755 int i, ii, resnr, c;
1756 char *atomname, *resname;
1763 natoms = dd->comm->nat[ddnatVSITE];
1766 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
1768 out = gmx_fio_fopen(fname, "w");
1770 fprintf(out, "TITLE %s\n", title);
1771 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1772 for (i = 0; i < natoms; i++)
1774 ii = dd->gatindex[i];
1775 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
1776 if (i < dd->comm->nat[ddnatZONE])
1779 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
1785 else if (i < dd->comm->nat[ddnatVSITE])
1787 b = dd->comm->zones.n;
1791 b = dd->comm->zones.n + 1;
1793 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
1794 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
1796 fprintf(out, "TER\n");
1798 gmx_fio_fclose(out);
1801 real dd_cutoff_multibody(const gmx_domdec_t *dd)
1803 gmx_domdec_comm_t *comm;
1810 if (comm->bInterCGBondeds)
1812 if (comm->cutoff_mbody > 0)
1814 r = comm->cutoff_mbody;
1818 /* cutoff_mbody=0 means we do not have DLB */
1819 r = comm->cellsize_min[dd->dim[0]];
1820 for (di = 1; di < dd->ndim; di++)
1822 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
1824 if (comm->bBondComm)
1826 r = std::max(r, comm->cutoff_mbody);
1830 r = std::min(r, comm->cutoff);
1838 real dd_cutoff_twobody(const gmx_domdec_t *dd)
1842 r_mb = dd_cutoff_multibody(dd);
1844 return std::max(dd->comm->cutoff, r_mb);
1848 static void dd_cart_coord2pmecoord(const gmx_domdec_t *dd, const ivec coord,
1853 nc = dd->nc[dd->comm->cartpmedim];
1854 ntot = dd->comm->ntot[dd->comm->cartpmedim];
1855 copy_ivec(coord, coord_pme);
1856 coord_pme[dd->comm->cartpmedim] =
1857 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
1860 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
1865 npme = dd->comm->npmenodes;
1867 /* Here we assign a PME node to communicate with this DD node
1868 * by assuming that the major index of both is x.
1869 * We add cr->npmenodes/2 to obtain an even distribution.
1871 return (ddindex*npme + npme/2)/npp;
1874 static int *dd_interleaved_pme_ranks(const gmx_domdec_t *dd)
1879 snew(pme_rank, dd->comm->npmenodes);
1881 for (i = 0; i < dd->nnodes; i++)
1883 p0 = ddindex2pmeindex(dd, i);
1884 p1 = ddindex2pmeindex(dd, i+1);
1885 if (i+1 == dd->nnodes || p1 > p0)
1889 fprintf(debug, "pme_rank[%d] = %d\n", n, i+1+n);
1891 pme_rank[n] = i + 1 + n;
1899 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
1907 if (dd->comm->bCartesian) {
1908 gmx_ddindex2xyz(dd->nc,ddindex,coords);
1909 dd_coords2pmecoords(dd,coords,coords_pme);
1910 copy_ivec(dd->ntot,nc);
1911 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
1912 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
1914 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
1916 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
1922 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
1927 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
1929 gmx_domdec_comm_t *comm;
1931 int ddindex, nodeid = -1;
1933 comm = cr->dd->comm;
1938 if (comm->bCartesianPP_PME)
1941 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
1946 ddindex = dd_index(cr->dd->nc, coords);
1947 if (comm->bCartesianPP)
1949 nodeid = comm->ddindex2simnodeid[ddindex];
1955 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
1967 static int dd_simnode2pmenode(const gmx_domdec_t *dd,
1968 const t_commrec gmx_unused *cr,
1973 const gmx_domdec_comm_t *comm = dd->comm;
1975 /* This assumes a uniform x domain decomposition grid cell size */
1976 if (comm->bCartesianPP_PME)
1979 ivec coord, coord_pme;
1980 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
1981 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
1983 /* This is a PP node */
1984 dd_cart_coord2pmecoord(dd, coord, coord_pme);
1985 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
1989 else if (comm->bCartesianPP)
1991 if (sim_nodeid < dd->nnodes)
1993 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
1998 /* This assumes DD cells with identical x coordinates
1999 * are numbered sequentially.
2001 if (dd->comm->pmenodes == NULL)
2003 if (sim_nodeid < dd->nnodes)
2005 /* The DD index equals the nodeid */
2006 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2012 while (sim_nodeid > dd->comm->pmenodes[i])
2016 if (sim_nodeid < dd->comm->pmenodes[i])
2018 pmenode = dd->comm->pmenodes[i];
2026 void get_pme_nnodes(const gmx_domdec_t *dd,
2027 int *npmenodes_x, int *npmenodes_y)
2031 *npmenodes_x = dd->comm->npmenodes_x;
2032 *npmenodes_y = dd->comm->npmenodes_y;
2041 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2042 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2046 ivec coord, coord_pme;
2050 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2053 for (x = 0; x < dd->nc[XX]; x++)
2055 for (y = 0; y < dd->nc[YY]; y++)
2057 for (z = 0; z < dd->nc[ZZ]; z++)
2059 if (dd->comm->bCartesianPP_PME)
2064 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2065 if (dd->ci[XX] == coord_pme[XX] &&
2066 dd->ci[YY] == coord_pme[YY] &&
2067 dd->ci[ZZ] == coord_pme[ZZ])
2069 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2074 /* The slab corresponds to the nodeid in the PME group */
2075 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2077 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2084 /* The last PP-only node is the peer node */
2085 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2089 fprintf(debug, "Receive coordinates from PP ranks:");
2090 for (x = 0; x < *nmy_ddnodes; x++)
2092 fprintf(debug, " %d", (*my_ddnodes)[x]);
2094 fprintf(debug, "\n");
2098 static gmx_bool receive_vir_ener(const gmx_domdec_t *dd, const t_commrec *cr)
2100 gmx_bool bReceive = TRUE;
2102 if (cr->npmenodes < dd->nnodes)
2104 gmx_domdec_comm_t *comm = dd->comm;
2105 if (comm->bCartesianPP_PME)
2108 int pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
2110 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2111 coords[comm->cartpmedim]++;
2112 if (coords[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2115 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2116 if (dd_simnode2pmenode(dd, cr, rank) == pmenode)
2118 /* This is not the last PP node for pmenode */
2123 GMX_RELEASE_ASSERT(false, "Without MPI we should not have Cartesian PP-PME with #PMEnodes < #DDnodes");
2128 int pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
2129 if (cr->sim_nodeid+1 < cr->nnodes &&
2130 dd_simnode2pmenode(dd, cr, cr->sim_nodeid+1) == pmenode)
2132 /* This is not the last PP node for pmenode */
2141 static void set_zones_ncg_home(gmx_domdec_t *dd)
2143 gmx_domdec_zones_t *zones;
2146 zones = &dd->comm->zones;
2148 zones->cg_range[0] = 0;
2149 for (i = 1; i < zones->n+1; i++)
2151 zones->cg_range[i] = dd->ncg_home;
2153 /* zone_ncg1[0] should always be equal to ncg_home */
2154 dd->comm->zone_ncg1[0] = dd->ncg_home;
2157 static void rebuild_cgindex(gmx_domdec_t *dd,
2158 const int *gcgs_index, const t_state *state)
2160 int * gmx_restrict dd_cg_gl = dd->index_gl;
2161 int * gmx_restrict cgindex = dd->cgindex;
2164 /* Copy back the global charge group indices from state
2165 * and rebuild the local charge group to atom index.
2168 for (unsigned int i = 0; i < state->cg_gl.size(); i++)
2171 int cg_gl = state->cg_gl[i];
2172 dd_cg_gl[i] = cg_gl;
2173 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2175 cgindex[state->cg_gl.size()] = nat;
2177 dd->ncg_home = state->cg_gl.size();
2180 set_zones_ncg_home(dd);
2183 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2185 while (cg >= cginfo_mb->cg_end)
2190 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2193 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2194 t_forcerec *fr, char *bLocalCG)
2196 cginfo_mb_t *cginfo_mb;
2202 cginfo_mb = fr->cginfo_mb;
2203 cginfo = fr->cginfo;
2205 for (cg = cg0; cg < cg1; cg++)
2207 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2211 if (bLocalCG != NULL)
2213 for (cg = cg0; cg < cg1; cg++)
2215 bLocalCG[index_gl[cg]] = TRUE;
2220 static void make_dd_indices(gmx_domdec_t *dd,
2221 const int *gcgs_index, int cg_start)
2223 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2224 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2227 if (dd->nat_tot > dd->gatindex_nalloc)
2229 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2230 srenew(dd->gatindex, dd->gatindex_nalloc);
2233 nzone = dd->comm->zones.n;
2234 zone2cg = dd->comm->zones.cg_range;
2235 zone_ncg1 = dd->comm->zone_ncg1;
2236 index_gl = dd->index_gl;
2237 gatindex = dd->gatindex;
2238 bCGs = dd->comm->bCGs;
2240 if (zone2cg[1] != dd->ncg_home)
2242 gmx_incons("dd->ncg_zone is not up to date");
2245 /* Make the local to global and global to local atom index */
2246 a = dd->cgindex[cg_start];
2247 for (zone = 0; zone < nzone; zone++)
2255 cg0 = zone2cg[zone];
2257 cg1 = zone2cg[zone+1];
2258 cg1_p1 = cg0 + zone_ncg1[zone];
2260 for (cg = cg0; cg < cg1; cg++)
2265 /* Signal that this cg is from more than one pulse away */
2268 cg_gl = index_gl[cg];
2271 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2274 ga2la_set(dd->ga2la, a_gl, a, zone1);
2280 gatindex[a] = cg_gl;
2281 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2288 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2294 if (bLocalCG == NULL)
2298 for (i = 0; i < dd->ncg_tot; i++)
2300 if (!bLocalCG[dd->index_gl[i]])
2303 "DD rank %d, %s: cg %d, global cg %d is not marked in bLocalCG (ncg_home %d)\n", dd->rank, where, i+1, dd->index_gl[i]+1, dd->ncg_home);
2308 for (i = 0; i < ncg_sys; i++)
2315 if (ngl != dd->ncg_tot)
2317 fprintf(stderr, "DD rank %d, %s: In bLocalCG %d cgs are marked as local, whereas there are %d\n", dd->rank, where, ngl, dd->ncg_tot);
2324 static void check_index_consistency(gmx_domdec_t *dd,
2325 int natoms_sys, int ncg_sys,
2328 int nerr, ngl, i, a, cell;
2333 if (dd->comm->DD_debug > 1)
2335 snew(have, natoms_sys);
2336 for (a = 0; a < dd->nat_tot; a++)
2338 if (have[dd->gatindex[a]] > 0)
2340 fprintf(stderr, "DD rank %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
2344 have[dd->gatindex[a]] = a + 1;
2350 snew(have, dd->nat_tot);
2353 for (i = 0; i < natoms_sys; i++)
2355 if (ga2la_get(dd->ga2la, i, &a, &cell))
2357 if (a >= dd->nat_tot)
2359 fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which is larger than nat_tot (%d)\n", dd->rank, i+1, a+1, dd->nat_tot);
2365 if (dd->gatindex[a] != i)
2367 fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which has global atom index %d\n", dd->rank, i+1, a+1, dd->gatindex[a]+1);
2374 if (ngl != dd->nat_tot)
2377 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2378 dd->rank, where, ngl, dd->nat_tot);
2380 for (a = 0; a < dd->nat_tot; a++)
2385 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2386 dd->rank, where, a+1, dd->gatindex[a]+1);
2391 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2395 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2396 dd->rank, where, nerr);
2400 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2407 /* Clear the whole list without searching */
2408 ga2la_clear(dd->ga2la);
2412 for (i = a_start; i < dd->nat_tot; i++)
2414 ga2la_del(dd->ga2la, dd->gatindex[i]);
2418 bLocalCG = dd->comm->bLocalCG;
2421 for (i = cg_start; i < dd->ncg_tot; i++)
2423 bLocalCG[dd->index_gl[i]] = FALSE;
2427 dd_clear_local_vsite_indices(dd);
2429 if (dd->constraints)
2431 dd_clear_local_constraint_indices(dd);
2435 /* This function should be used for moving the domain boudaries during DLB,
2436 * for obtaining the minimum cell size. It checks the initially set limit
2437 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2438 * and, possibly, a longer cut-off limit set for PME load balancing.
2440 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2444 cellsize_min = comm->cellsize_min[dim];
2446 if (!comm->bVacDLBNoLimit)
2448 /* The cut-off might have changed, e.g. by PME load balacning,
2449 * from the value used to set comm->cellsize_min, so check it.
2451 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2453 if (comm->bPMELoadBalDLBLimits)
2455 /* Check for the cut-off limit set by the PME load balancing */
2456 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2460 return cellsize_min;
2463 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2466 real grid_jump_limit;
2468 /* The distance between the boundaries of cells at distance
2469 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2470 * and by the fact that cells should not be shifted by more than
2471 * half their size, such that cg's only shift by one cell
2472 * at redecomposition.
2474 grid_jump_limit = comm->cellsize_limit;
2475 if (!comm->bVacDLBNoLimit)
2477 if (comm->bPMELoadBalDLBLimits)
2479 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2481 grid_jump_limit = std::max(grid_jump_limit,
2482 cutoff/comm->cd[dim_ind].np);
2485 return grid_jump_limit;
2488 static gmx_bool check_grid_jump(gmx_int64_t step,
2494 gmx_domdec_comm_t *comm;
2503 for (d = 1; d < dd->ndim; d++)
2506 limit = grid_jump_limit(comm, cutoff, d);
2507 bfac = ddbox->box_size[dim];
2508 if (ddbox->tric_dir[dim])
2510 bfac *= ddbox->skew_fac[dim];
2512 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2513 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2521 /* This error should never be triggered under normal
2522 * circumstances, but you never know ...
2524 gmx_fatal(FARGS, "Step %s: The domain decomposition grid has shifted too much in the %c-direction around cell %d %d %d. This should not have happened. Running with fewer ranks might avoid this issue.",
2525 gmx_step_str(step, buf),
2526 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2534 static int dd_load_count(gmx_domdec_comm_t *comm)
2536 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2539 static float dd_force_load(gmx_domdec_comm_t *comm)
2546 if (comm->eFlop > 1)
2548 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2553 load = comm->cycl[ddCyclF];
2554 if (comm->cycl_n[ddCyclF] > 1)
2556 /* Subtract the maximum of the last n cycle counts
2557 * to get rid of possible high counts due to other sources,
2558 * for instance system activity, that would otherwise
2559 * affect the dynamic load balancing.
2561 load -= comm->cycl_max[ddCyclF];
2565 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2567 float gpu_wait, gpu_wait_sum;
2569 gpu_wait = comm->cycl[ddCyclWaitGPU];
2570 if (comm->cycl_n[ddCyclF] > 1)
2572 /* We should remove the WaitGPU time of the same MD step
2573 * as the one with the maximum F time, since the F time
2574 * and the wait time are not independent.
2575 * Furthermore, the step for the max F time should be chosen
2576 * the same on all ranks that share the same GPU.
2577 * But to keep the code simple, we remove the average instead.
2578 * The main reason for artificially long times at some steps
2579 * is spurious CPU activity or MPI time, so we don't expect
2580 * that changes in the GPU wait time matter a lot here.
2582 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2584 /* Sum the wait times over the ranks that share the same GPU */
2585 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2586 comm->mpi_comm_gpu_shared);
2587 /* Replace the wait time by the average over the ranks */
2588 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2596 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2598 gmx_domdec_comm_t *comm;
2603 snew(*dim_f, dd->nc[dim]+1);
2605 for (i = 1; i < dd->nc[dim]; i++)
2607 if (comm->slb_frac[dim])
2609 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2613 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2616 (*dim_f)[dd->nc[dim]] = 1;
2619 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2621 int pmeindex, slab, nso, i;
2624 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2630 ddpme->dim = dimind;
2632 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2634 ddpme->nslab = (ddpme->dim == 0 ?
2635 dd->comm->npmenodes_x :
2636 dd->comm->npmenodes_y);
2638 if (ddpme->nslab <= 1)
2643 nso = dd->comm->npmenodes/ddpme->nslab;
2644 /* Determine for each PME slab the PP location range for dimension dim */
2645 snew(ddpme->pp_min, ddpme->nslab);
2646 snew(ddpme->pp_max, ddpme->nslab);
2647 for (slab = 0; slab < ddpme->nslab; slab++)
2649 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2650 ddpme->pp_max[slab] = 0;
2652 for (i = 0; i < dd->nnodes; i++)
2654 ddindex2xyz(dd->nc, i, xyz);
2655 /* For y only use our y/z slab.
2656 * This assumes that the PME x grid size matches the DD grid size.
2658 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2660 pmeindex = ddindex2pmeindex(dd, i);
2663 slab = pmeindex/nso;
2667 slab = pmeindex % ddpme->nslab;
2669 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2670 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2674 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2677 int dd_pme_maxshift_x(const gmx_domdec_t *dd)
2679 if (dd->comm->ddpme[0].dim == XX)
2681 return dd->comm->ddpme[0].maxshift;
2689 int dd_pme_maxshift_y(const gmx_domdec_t *dd)
2691 if (dd->comm->ddpme[0].dim == YY)
2693 return dd->comm->ddpme[0].maxshift;
2695 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2697 return dd->comm->ddpme[1].maxshift;
2705 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2706 gmx_bool bUniform, const gmx_ddbox_t *ddbox,
2709 gmx_domdec_comm_t *comm;
2712 real range, pme_boundary;
2716 nc = dd->nc[ddpme->dim];
2719 if (!ddpme->dim_match)
2721 /* PP decomposition is not along dim: the worst situation */
2724 else if (ns <= 3 || (bUniform && ns == nc))
2726 /* The optimal situation */
2731 /* We need to check for all pme nodes which nodes they
2732 * could possibly need to communicate with.
2734 xmin = ddpme->pp_min;
2735 xmax = ddpme->pp_max;
2736 /* Allow for atoms to be maximally 2/3 times the cut-off
2737 * out of their DD cell. This is a reasonable balance between
2738 * between performance and support for most charge-group/cut-off
2741 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2742 /* Avoid extra communication when we are exactly at a boundary */
2746 for (s = 0; s < ns; s++)
2748 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2749 pme_boundary = (real)s/ns;
2752 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
2754 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
2758 pme_boundary = (real)(s+1)/ns;
2761 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
2763 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
2770 ddpme->maxshift = sh;
2774 fprintf(debug, "PME slab communication range for dim %d is %d\n",
2775 ddpme->dim, ddpme->maxshift);
2779 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
2783 for (d = 0; d < dd->ndim; d++)
2786 if (dim < ddbox->nboundeddim &&
2787 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
2788 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
2790 gmx_fatal(FARGS, "The %c-size of the box (%f) times the triclinic skew factor (%f) is smaller than the number of DD cells (%d) times the smallest allowed cell size (%f)\n",
2791 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
2792 dd->nc[dim], dd->comm->cellsize_limit);
2798 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
2801 /* Set the domain boundaries. Use for static (or no) load balancing,
2802 * and also for the starting state for dynamic load balancing.
2803 * setmode determine if and where the boundaries are stored, use enum above.
2804 * Returns the number communication pulses in npulse.
2806 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, const gmx_ddbox_t *ddbox,
2807 int setmode, ivec npulse)
2809 gmx_domdec_comm_t *comm;
2812 real *cell_x, cell_dx, cellsize;
2816 for (d = 0; d < DIM; d++)
2818 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
2820 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
2823 cell_dx = ddbox->box_size[d]/dd->nc[d];
2826 case setcellsizeslbMASTER:
2827 for (j = 0; j < dd->nc[d]+1; j++)
2829 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
2832 case setcellsizeslbLOCAL:
2833 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
2834 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
2839 cellsize = cell_dx*ddbox->skew_fac[d];
2840 while (cellsize*npulse[d] < comm->cutoff)
2844 cellsize_min[d] = cellsize;
2848 /* Statically load balanced grid */
2849 /* Also when we are not doing a master distribution we determine
2850 * all cell borders in a loop to obtain identical values
2851 * to the master distribution case and to determine npulse.
2853 if (setmode == setcellsizeslbMASTER)
2855 cell_x = dd->ma->cell_x[d];
2859 snew(cell_x, dd->nc[d]+1);
2861 cell_x[0] = ddbox->box0[d];
2862 for (j = 0; j < dd->nc[d]; j++)
2864 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
2865 cell_x[j+1] = cell_x[j] + cell_dx;
2866 cellsize = cell_dx*ddbox->skew_fac[d];
2867 while (cellsize*npulse[d] < comm->cutoff &&
2868 npulse[d] < dd->nc[d]-1)
2872 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
2874 if (setmode == setcellsizeslbLOCAL)
2876 comm->cell_x0[d] = cell_x[dd->ci[d]];
2877 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
2879 if (setmode != setcellsizeslbMASTER)
2884 /* The following limitation is to avoid that a cell would receive
2885 * some of its own home charge groups back over the periodic boundary.
2886 * Double charge groups cause trouble with the global indices.
2888 if (d < ddbox->npbcdim &&
2889 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
2891 char error_string[STRLEN];
2893 sprintf(error_string,
2894 "The box size in direction %c (%f) times the triclinic skew factor (%f) is too small for a cut-off of %f with %d domain decomposition cells, use 1 or more than %d %s or increase the box size in this direction",
2895 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
2897 dd->nc[d], dd->nc[d],
2898 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
2900 if (setmode == setcellsizeslbLOCAL)
2902 gmx_fatal_collective(FARGS, dd->mpi_comm_all, DDMASTER(dd),
2907 gmx_fatal(FARGS, error_string);
2914 copy_rvec(cellsize_min, comm->cellsize_min);
2917 for (d = 0; d < comm->npmedecompdim; d++)
2919 set_pme_maxshift(dd, &comm->ddpme[d],
2920 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
2921 comm->ddpme[d].slb_dim_f);
2926 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
2927 int d, int dim, domdec_root_t *root,
2928 const gmx_ddbox_t *ddbox,
2929 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
2931 gmx_domdec_comm_t *comm;
2932 int ncd, i, j, nmin, nmin_old;
2933 gmx_bool bLimLo, bLimHi;
2935 real fac, halfway, cellsize_limit_f_i, region_size;
2936 gmx_bool bPBC, bLastHi = FALSE;
2937 int nrange[] = {range[0], range[1]};
2939 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
2945 bPBC = (dim < ddbox->npbcdim);
2947 cell_size = root->buf_ncd;
2951 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
2954 /* First we need to check if the scaling does not make cells
2955 * smaller than the smallest allowed size.
2956 * We need to do this iteratively, since if a cell is too small,
2957 * it needs to be enlarged, which makes all the other cells smaller,
2958 * which could in turn make another cell smaller than allowed.
2960 for (i = range[0]; i < range[1]; i++)
2962 root->bCellMin[i] = FALSE;
2968 /* We need the total for normalization */
2970 for (i = range[0]; i < range[1]; i++)
2972 if (root->bCellMin[i] == FALSE)
2974 fac += cell_size[i];
2977 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
2978 /* Determine the cell boundaries */
2979 for (i = range[0]; i < range[1]; i++)
2981 if (root->bCellMin[i] == FALSE)
2983 cell_size[i] *= fac;
2984 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
2986 cellsize_limit_f_i = 0;
2990 cellsize_limit_f_i = cellsize_limit_f;
2992 if (cell_size[i] < cellsize_limit_f_i)
2994 root->bCellMin[i] = TRUE;
2995 cell_size[i] = cellsize_limit_f_i;
2999 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3002 while (nmin > nmin_old);
3005 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3006 /* For this check we should not use DD_CELL_MARGIN,
3007 * but a slightly smaller factor,
3008 * since rounding could get use below the limit.
3010 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3013 gmx_fatal(FARGS, "Step %s: the dynamic load balancing could not balance dimension %c: box size %f, triclinic skew factor %f, #cells %d, minimum cell size %f\n",
3014 gmx_step_str(step, buf),
3015 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3016 ncd, comm->cellsize_min[dim]);
3019 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3023 /* Check if the boundary did not displace more than halfway
3024 * each of the cells it bounds, as this could cause problems,
3025 * especially when the differences between cell sizes are large.
3026 * If changes are applied, they will not make cells smaller
3027 * than the cut-off, as we check all the boundaries which
3028 * might be affected by a change and if the old state was ok,
3029 * the cells will at most be shrunk back to their old size.
3031 for (i = range[0]+1; i < range[1]; i++)
3033 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3034 if (root->cell_f[i] < halfway)
3036 root->cell_f[i] = halfway;
3037 /* Check if the change also causes shifts of the next boundaries */
3038 for (j = i+1; j < range[1]; j++)
3040 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3042 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3046 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3047 if (root->cell_f[i] > halfway)
3049 root->cell_f[i] = halfway;
3050 /* Check if the change also causes shifts of the next boundaries */
3051 for (j = i-1; j >= range[0]+1; j--)
3053 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3055 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3062 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3063 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3064 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3065 * for a and b nrange is used */
3068 /* Take care of the staggering of the cell boundaries */
3071 for (i = range[0]; i < range[1]; i++)
3073 root->cell_f_max0[i] = root->cell_f[i];
3074 root->cell_f_min1[i] = root->cell_f[i+1];
3079 for (i = range[0]+1; i < range[1]; i++)
3081 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3082 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3083 if (bLimLo && bLimHi)
3085 /* Both limits violated, try the best we can */
3086 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3087 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3088 nrange[0] = range[0];
3090 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3093 nrange[1] = range[1];
3094 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3100 /* root->cell_f[i] = root->bound_min[i]; */
3101 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3104 else if (bLimHi && !bLastHi)
3107 if (nrange[1] < range[1]) /* found a LimLo before */
3109 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3110 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3111 nrange[0] = nrange[1];
3113 root->cell_f[i] = root->bound_max[i];
3115 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3117 nrange[1] = range[1];
3120 if (nrange[1] < range[1]) /* found last a LimLo */
3122 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3123 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3124 nrange[0] = nrange[1];
3125 nrange[1] = range[1];
3126 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3128 else if (nrange[0] > range[0]) /* found at least one LimHi */
3130 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3137 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3138 int d, int dim, domdec_root_t *root,
3139 const gmx_ddbox_t *ddbox,
3140 gmx_bool bDynamicBox,
3141 gmx_bool bUniform, gmx_int64_t step)
3143 gmx_domdec_comm_t *comm;
3144 int ncd, d1, i, pos;
3146 real load_aver, load_i, imbalance, change, change_max, sc;
3147 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3151 int range[] = { 0, 0 };
3155 /* Convert the maximum change from the input percentage to a fraction */
3156 change_limit = comm->dlb_scale_lim*0.01;
3160 bPBC = (dim < ddbox->npbcdim);
3162 cell_size = root->buf_ncd;
3164 /* Store the original boundaries */
3165 for (i = 0; i < ncd+1; i++)
3167 root->old_cell_f[i] = root->cell_f[i];
3171 for (i = 0; i < ncd; i++)
3173 cell_size[i] = 1.0/ncd;
3176 else if (dd_load_count(comm) > 0)
3178 load_aver = comm->load[d].sum_m/ncd;
3180 for (i = 0; i < ncd; i++)
3182 /* Determine the relative imbalance of cell i */
3183 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3184 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3185 /* Determine the change of the cell size using underrelaxation */
3186 change = -relax*imbalance;
3187 change_max = std::max(change_max, std::max(change, -change));
3189 /* Limit the amount of scaling.
3190 * We need to use the same rescaling for all cells in one row,
3191 * otherwise the load balancing might not converge.
3194 if (change_max > change_limit)
3196 sc *= change_limit/change_max;
3198 for (i = 0; i < ncd; i++)
3200 /* Determine the relative imbalance of cell i */
3201 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3202 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3203 /* Determine the change of the cell size using underrelaxation */
3204 change = -sc*imbalance;
3205 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3209 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3210 cellsize_limit_f *= DD_CELL_MARGIN;
3211 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3212 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3213 if (ddbox->tric_dir[dim])
3215 cellsize_limit_f /= ddbox->skew_fac[dim];
3216 dist_min_f /= ddbox->skew_fac[dim];
3218 if (bDynamicBox && d > 0)
3220 dist_min_f *= DD_PRES_SCALE_MARGIN;
3222 if (d > 0 && !bUniform)
3224 /* Make sure that the grid is not shifted too much */
3225 for (i = 1; i < ncd; i++)
3227 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3229 gmx_incons("Inconsistent DD boundary staggering limits!");
3231 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3232 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3235 root->bound_min[i] += 0.5*space;
3237 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3238 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3241 root->bound_max[i] += 0.5*space;
3246 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3248 root->cell_f_max0[i-1] + dist_min_f,
3249 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3250 root->cell_f_min1[i] - dist_min_f);
3255 root->cell_f[0] = 0;
3256 root->cell_f[ncd] = 1;
3257 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3260 /* After the checks above, the cells should obey the cut-off
3261 * restrictions, but it does not hurt to check.
3263 for (i = 0; i < ncd; i++)
3267 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3268 dim, i, root->cell_f[i], root->cell_f[i+1]);
3271 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3272 root->cell_f[i+1] - root->cell_f[i] <
3273 cellsize_limit_f/DD_CELL_MARGIN)
3277 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3278 gmx_step_str(step, buf), dim2char(dim), i,
3279 (root->cell_f[i+1] - root->cell_f[i])
3280 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3285 /* Store the cell boundaries of the lower dimensions at the end */
3286 for (d1 = 0; d1 < d; d1++)
3288 root->cell_f[pos++] = comm->cell_f0[d1];
3289 root->cell_f[pos++] = comm->cell_f1[d1];
3292 if (d < comm->npmedecompdim)
3294 /* The master determines the maximum shift for
3295 * the coordinate communication between separate PME nodes.
3297 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3299 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3302 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3306 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3307 const gmx_ddbox_t *ddbox,
3310 gmx_domdec_comm_t *comm;
3315 /* Set the cell dimensions */
3316 dim = dd->dim[dimind];
3317 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3318 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3319 if (dim >= ddbox->nboundeddim)
3321 comm->cell_x0[dim] += ddbox->box0[dim];
3322 comm->cell_x1[dim] += ddbox->box0[dim];
3326 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3327 int d, int dim, real *cell_f_row,
3328 const gmx_ddbox_t *ddbox)
3330 gmx_domdec_comm_t *comm;
3336 /* Each node would only need to know two fractions,
3337 * but it is probably cheaper to broadcast the whole array.
3339 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3340 0, comm->mpi_comm_load[d]);
3342 /* Copy the fractions for this dimension from the buffer */
3343 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3344 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3345 /* The whole array was communicated, so set the buffer position */
3346 pos = dd->nc[dim] + 1;
3347 for (d1 = 0; d1 <= d; d1++)
3351 /* Copy the cell fractions of the lower dimensions */
3352 comm->cell_f0[d1] = cell_f_row[pos++];
3353 comm->cell_f1[d1] = cell_f_row[pos++];
3355 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3357 /* Convert the communicated shift from float to int */
3358 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3361 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3365 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3366 const gmx_ddbox_t *ddbox,
3367 gmx_bool bDynamicBox,
3368 gmx_bool bUniform, gmx_int64_t step)
3370 gmx_domdec_comm_t *comm;
3372 gmx_bool bRowMember, bRowRoot;
3377 for (d = 0; d < dd->ndim; d++)
3382 for (d1 = d; d1 < dd->ndim; d1++)
3384 if (dd->ci[dd->dim[d1]] > 0)
3397 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3398 ddbox, bDynamicBox, bUniform, step);
3399 cell_f_row = comm->root[d]->cell_f;
3403 cell_f_row = comm->cell_f_row;
3405 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3410 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd,
3411 const gmx_ddbox_t *ddbox)
3415 /* This function assumes the box is static and should therefore
3416 * not be called when the box has changed since the last
3417 * call to dd_partition_system.
3419 for (d = 0; d < dd->ndim; d++)
3421 relative_to_absolute_cell_bounds(dd, ddbox, d);
3427 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3428 const gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3429 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3430 gmx_wallcycle_t wcycle)
3432 gmx_domdec_comm_t *comm;
3439 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3440 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3441 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3443 else if (bDynamicBox)
3445 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3448 /* Set the dimensions for which no DD is used */
3449 for (dim = 0; dim < DIM; dim++)
3451 if (dd->nc[dim] == 1)
3453 comm->cell_x0[dim] = 0;
3454 comm->cell_x1[dim] = ddbox->box_size[dim];
3455 if (dim >= ddbox->nboundeddim)
3457 comm->cell_x0[dim] += ddbox->box0[dim];
3458 comm->cell_x1[dim] += ddbox->box0[dim];
3464 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3467 gmx_domdec_comm_dim_t *cd;
3469 for (d = 0; d < dd->ndim; d++)
3471 cd = &dd->comm->cd[d];
3472 np = npulse[dd->dim[d]];
3473 if (np > cd->np_nalloc)
3477 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3478 dim2char(dd->dim[d]), np);
3480 if (DDMASTER(dd) && cd->np_nalloc > 0)
3482 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3484 srenew(cd->ind, np);
3485 for (i = cd->np_nalloc; i < np; i++)
3487 cd->ind[i].index = NULL;
3488 cd->ind[i].nalloc = 0;
3497 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3498 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3499 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3500 gmx_wallcycle_t wcycle)
3502 gmx_domdec_comm_t *comm;
3508 /* Copy the old cell boundaries for the cg displacement check */
3509 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3510 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3516 check_box_size(dd, ddbox);
3518 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3522 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3523 realloc_comm_ind(dd, npulse);
3528 for (d = 0; d < DIM; d++)
3530 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3531 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3536 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3538 rvec cell_ns_x0, rvec cell_ns_x1,
3541 gmx_domdec_comm_t *comm;
3546 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3548 dim = dd->dim[dim_ind];
3550 /* Without PBC we don't have restrictions on the outer cells */
3551 if (!(dim >= ddbox->npbcdim &&
3552 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3554 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3555 comm->cellsize_min[dim])
3558 gmx_fatal(FARGS, "Step %s: The %c-size (%f) times the triclinic skew factor (%f) is smaller than the smallest allowed cell size (%f) for domain decomposition grid cell %d %d %d",
3559 gmx_step_str(step, buf), dim2char(dim),
3560 comm->cell_x1[dim] - comm->cell_x0[dim],
3561 ddbox->skew_fac[dim],
3562 dd->comm->cellsize_min[dim],
3563 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3567 if ((dlbIsOn(dd->comm) && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3569 /* Communicate the boundaries and update cell_ns_x0/1 */
3570 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3571 if (dlbIsOn(dd->comm) && dd->ndim > 1)
3573 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3578 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3582 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3590 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3591 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3600 static void check_screw_box(matrix box)
3602 /* Mathematical limitation */
3603 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3605 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3608 /* Limitation due to the asymmetry of the eighth shell method */
3609 if (box[ZZ][YY] != 0)
3611 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3615 static void distribute_cg(FILE *fplog,
3616 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3619 gmx_domdec_master_t *ma;
3620 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3621 int i, icg, j, k, k0, k1, d;
3625 real nrcg, inv_ncg, pos_d;
3631 snew(tmp_nalloc, dd->nnodes);
3632 snew(tmp_ind, dd->nnodes);
3633 for (i = 0; i < dd->nnodes; i++)
3635 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3636 snew(tmp_ind[i], tmp_nalloc[i]);
3639 /* Clear the count */
3640 for (i = 0; i < dd->nnodes; i++)
3646 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3648 cgindex = cgs->index;
3650 /* Compute the center of geometry for all charge groups */
3651 for (icg = 0; icg < cgs->nr; icg++)
3654 k1 = cgindex[icg+1];
3658 copy_rvec(pos[k0], cg_cm);
3665 for (k = k0; (k < k1); k++)
3667 rvec_inc(cg_cm, pos[k]);
3669 for (d = 0; (d < DIM); d++)
3671 cg_cm[d] *= inv_ncg;
3674 /* Put the charge group in the box and determine the cell index */
3675 for (d = DIM-1; d >= 0; d--)
3678 if (d < dd->npbcdim)
3680 bScrew = (dd->bScrewPBC && d == XX);
3681 if (tric_dir[d] && dd->nc[d] > 1)
3683 /* Use triclinic coordintates for this dimension */
3684 for (j = d+1; j < DIM; j++)
3686 pos_d += cg_cm[j]*tcm[j][d];
3689 while (pos_d >= box[d][d])
3692 rvec_dec(cg_cm, box[d]);
3695 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3696 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3698 for (k = k0; (k < k1); k++)
3700 rvec_dec(pos[k], box[d]);
3703 pos[k][YY] = box[YY][YY] - pos[k][YY];
3704 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3711 rvec_inc(cg_cm, box[d]);
3714 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3715 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3717 for (k = k0; (k < k1); k++)
3719 rvec_inc(pos[k], box[d]);
3722 pos[k][YY] = box[YY][YY] - pos[k][YY];
3723 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3728 /* This could be done more efficiently */
3730 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3735 i = dd_index(dd->nc, ind);
3736 if (ma->ncg[i] == tmp_nalloc[i])
3738 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3739 srenew(tmp_ind[i], tmp_nalloc[i]);
3741 tmp_ind[i][ma->ncg[i]] = icg;
3743 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3747 for (i = 0; i < dd->nnodes; i++)
3750 for (k = 0; k < ma->ncg[i]; k++)
3752 ma->cg[k1++] = tmp_ind[i][k];
3755 ma->index[dd->nnodes] = k1;
3757 for (i = 0; i < dd->nnodes; i++)
3766 // Use double for the sums to avoid natoms^2 overflowing
3768 int nat_sum, nat_min, nat_max;
3773 nat_min = ma->nat[0];
3774 nat_max = ma->nat[0];
3775 for (i = 0; i < dd->nnodes; i++)
3777 nat_sum += ma->nat[i];
3778 // cast to double to avoid integer overflows when squaring
3779 nat2_sum += gmx::square(static_cast<double>(ma->nat[i]));
3780 nat_min = std::min(nat_min, ma->nat[i]);
3781 nat_max = std::max(nat_max, ma->nat[i]);
3783 nat_sum /= dd->nnodes;
3784 nat2_sum /= dd->nnodes;
3786 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
3789 static_cast<int>(std::sqrt(nat2_sum - gmx::square(static_cast<double>(nat_sum)) + 0.5)),
3794 static void get_cg_distribution(FILE *fplog, gmx_domdec_t *dd,
3795 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
3798 gmx_domdec_master_t *ma = NULL;
3801 int *ibuf, buf2[2] = { 0, 0 };
3802 gmx_bool bMaster = DDMASTER(dd);
3810 check_screw_box(box);
3813 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
3815 distribute_cg(fplog, box, ddbox->tric_dir, cgs, pos, dd);
3816 for (i = 0; i < dd->nnodes; i++)
3818 ma->ibuf[2*i] = ma->ncg[i];
3819 ma->ibuf[2*i+1] = ma->nat[i];
3827 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
3829 dd->ncg_home = buf2[0];
3830 dd->nat_home = buf2[1];
3831 dd->ncg_tot = dd->ncg_home;
3832 dd->nat_tot = dd->nat_home;
3833 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
3835 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
3836 srenew(dd->index_gl, dd->cg_nalloc);
3837 srenew(dd->cgindex, dd->cg_nalloc+1);
3841 for (i = 0; i < dd->nnodes; i++)
3843 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
3844 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
3849 bMaster ? ma->ibuf : NULL,
3850 bMaster ? ma->ibuf+dd->nnodes : NULL,
3851 bMaster ? ma->cg : NULL,
3852 dd->ncg_home*sizeof(int), dd->index_gl);
3854 /* Determine the home charge group sizes */
3856 for (i = 0; i < dd->ncg_home; i++)
3858 cg_gl = dd->index_gl[i];
3860 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
3865 fprintf(debug, "Home charge groups:\n");
3866 for (i = 0; i < dd->ncg_home; i++)
3868 fprintf(debug, " %d", dd->index_gl[i]);
3871 fprintf(debug, "\n");
3874 fprintf(debug, "\n");
3878 static int compact_and_copy_vec_at(int ncg, int *move,
3881 rvec *src, gmx_domdec_comm_t *comm,
3884 int m, icg, i, i0, i1, nrcg;
3890 for (m = 0; m < DIM*2; m++)
3896 for (icg = 0; icg < ncg; icg++)
3898 i1 = cgindex[icg+1];
3904 /* Compact the home array in place */
3905 for (i = i0; i < i1; i++)
3907 copy_rvec(src[i], src[home_pos++]);
3913 /* Copy to the communication buffer */
3915 pos_vec[m] += 1 + vec*nrcg;
3916 for (i = i0; i < i1; i++)
3918 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
3920 pos_vec[m] += (nvec - vec - 1)*nrcg;
3924 home_pos += i1 - i0;
3932 static int compact_and_copy_vec_cg(int ncg, int *move,
3934 int nvec, rvec *src, gmx_domdec_comm_t *comm,
3937 int m, icg, i0, i1, nrcg;
3943 for (m = 0; m < DIM*2; m++)
3949 for (icg = 0; icg < ncg; icg++)
3951 i1 = cgindex[icg+1];
3957 /* Compact the home array in place */
3958 copy_rvec(src[icg], src[home_pos++]);
3964 /* Copy to the communication buffer */
3965 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
3966 pos_vec[m] += 1 + nrcg*nvec;
3978 static int compact_ind(int ncg, int *move,
3979 int *index_gl, int *cgindex,
3981 gmx_ga2la_t *ga2la, char *bLocalCG,
3984 int cg, nat, a0, a1, a, a_gl;
3989 for (cg = 0; cg < ncg; cg++)
3995 /* Compact the home arrays in place.
3996 * Anything that can be done here avoids access to global arrays.
3998 cgindex[home_pos] = nat;
3999 for (a = a0; a < a1; a++)
4002 gatindex[nat] = a_gl;
4003 /* The cell number stays 0, so we don't need to set it */
4004 ga2la_change_la(ga2la, a_gl, nat);
4007 index_gl[home_pos] = index_gl[cg];
4008 cginfo[home_pos] = cginfo[cg];
4009 /* The charge group remains local, so bLocalCG does not change */
4014 /* Clear the global indices */
4015 for (a = a0; a < a1; a++)
4017 ga2la_del(ga2la, gatindex[a]);
4021 bLocalCG[index_gl[cg]] = FALSE;
4025 cgindex[home_pos] = nat;
4030 static void clear_and_mark_ind(int ncg, int *move,
4031 int *index_gl, int *cgindex, int *gatindex,
4032 gmx_ga2la_t *ga2la, char *bLocalCG,
4037 for (cg = 0; cg < ncg; cg++)
4043 /* Clear the global indices */
4044 for (a = a0; a < a1; a++)
4046 ga2la_del(ga2la, gatindex[a]);
4050 bLocalCG[index_gl[cg]] = FALSE;
4052 /* Signal that this cg has moved using the ns cell index.
4053 * Here we set it to -1. fill_grid will change it
4054 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4056 cell_index[cg] = -1;
4061 static void print_cg_move(FILE *fplog,
4063 gmx_int64_t step, int cg, int dim, int dir,
4064 gmx_bool bHaveCgcmOld, real limitd,
4065 rvec cm_old, rvec cm_new, real pos_d)
4067 gmx_domdec_comm_t *comm;
4072 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4075 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4076 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4077 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4081 /* We don't have a limiting distance available: don't print it */
4082 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4083 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4084 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4086 fprintf(fplog, "distance out of cell %f\n",
4087 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4090 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4091 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4093 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4094 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4095 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4097 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4098 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4100 comm->cell_x0[dim], comm->cell_x1[dim]);
4103 static void cg_move_error(FILE *fplog,
4105 gmx_int64_t step, int cg, int dim, int dir,
4106 gmx_bool bHaveCgcmOld, real limitd,
4107 rvec cm_old, rvec cm_new, real pos_d)
4111 print_cg_move(fplog, dd, step, cg, dim, dir,
4112 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4114 print_cg_move(stderr, dd, step, cg, dim, dir,
4115 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4117 "%s moved too far between two domain decomposition steps\n"
4118 "This usually means that your system is not well equilibrated",
4119 dd->comm->bCGs ? "A charge group" : "An atom");
4122 static void rotate_state_atom(t_state *state, int a)
4126 for (est = 0; est < estNR; est++)
4128 if (EST_DISTR(est) && (state->flags & (1<<est)))
4133 /* Rotate the complete state; for a rectangular box only */
4134 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4135 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4138 state->v[a][YY] = -state->v[a][YY];
4139 state->v[a][ZZ] = -state->v[a][ZZ];
4141 case est_SDX_NOTSUPPORTED:
4144 state->cg_p[a][YY] = -state->cg_p[a][YY];
4145 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4147 case estDISRE_INITF:
4148 case estDISRE_RM3TAV:
4149 case estORIRE_INITF:
4151 /* These are distances, so not affected by rotation */
4154 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4160 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4162 if (natoms > comm->moved_nalloc)
4164 /* Contents should be preserved here */
4165 comm->moved_nalloc = over_alloc_dd(natoms);
4166 srenew(comm->moved, comm->moved_nalloc);
4172 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4175 ivec tric_dir, matrix tcm,
4176 rvec cell_x0, rvec cell_x1,
4177 rvec limitd, rvec limit0, rvec limit1,
4179 int cg_start, int cg_end,
4184 int cg, k, k0, k1, d, dim, d2;
4189 real inv_ncg, pos_d;
4192 npbcdim = dd->npbcdim;
4194 for (cg = cg_start; cg < cg_end; cg++)
4201 copy_rvec(state->x[k0], cm_new);
4208 for (k = k0; (k < k1); k++)
4210 rvec_inc(cm_new, state->x[k]);
4212 for (d = 0; (d < DIM); d++)
4214 cm_new[d] = inv_ncg*cm_new[d];
4219 /* Do pbc and check DD cell boundary crossings */
4220 for (d = DIM-1; d >= 0; d--)
4224 bScrew = (dd->bScrewPBC && d == XX);
4225 /* Determine the location of this cg in lattice coordinates */
4229 for (d2 = d+1; d2 < DIM; d2++)
4231 pos_d += cm_new[d2]*tcm[d2][d];
4234 /* Put the charge group in the triclinic unit-cell */
4235 if (pos_d >= cell_x1[d])
4237 if (pos_d >= limit1[d])
4239 cg_move_error(fplog, dd, step, cg, d, 1,
4240 cg_cm != as_rvec_array(state->x.data()), limitd[d],
4241 cg_cm[cg], cm_new, pos_d);
4244 if (dd->ci[d] == dd->nc[d] - 1)
4246 rvec_dec(cm_new, state->box[d]);
4249 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4250 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4252 for (k = k0; (k < k1); k++)
4254 rvec_dec(state->x[k], state->box[d]);
4257 rotate_state_atom(state, k);
4262 else if (pos_d < cell_x0[d])
4264 if (pos_d < limit0[d])
4266 cg_move_error(fplog, dd, step, cg, d, -1,
4267 cg_cm != as_rvec_array(state->x.data()), limitd[d],
4268 cg_cm[cg], cm_new, pos_d);
4273 rvec_inc(cm_new, state->box[d]);
4276 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4277 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4279 for (k = k0; (k < k1); k++)
4281 rvec_inc(state->x[k], state->box[d]);
4284 rotate_state_atom(state, k);
4290 else if (d < npbcdim)
4292 /* Put the charge group in the rectangular unit-cell */
4293 while (cm_new[d] >= state->box[d][d])
4295 rvec_dec(cm_new, state->box[d]);
4296 for (k = k0; (k < k1); k++)
4298 rvec_dec(state->x[k], state->box[d]);
4301 while (cm_new[d] < 0)
4303 rvec_inc(cm_new, state->box[d]);
4304 for (k = k0; (k < k1); k++)
4306 rvec_inc(state->x[k], state->box[d]);
4312 copy_rvec(cm_new, cg_cm[cg]);
4314 /* Determine where this cg should go */
4317 for (d = 0; d < dd->ndim; d++)
4322 flag |= DD_FLAG_FW(d);
4328 else if (dev[dim] == -1)
4330 flag |= DD_FLAG_BW(d);
4333 if (dd->nc[dim] > 2)
4344 /* Temporarily store the flag in move */
4345 move[cg] = mc + flag;
4349 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4350 gmx_domdec_t *dd, ivec tric_dir,
4351 t_state *state, PaddedRVecVector *f,
4360 int ncg[DIM*2], nat[DIM*2];
4361 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4362 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4363 int sbuf[2], rbuf[2];
4364 int home_pos_cg, home_pos_at, buf_pos;
4366 gmx_bool bV = FALSE, bCGP = FALSE;
4369 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4371 cginfo_mb_t *cginfo_mb;
4372 gmx_domdec_comm_t *comm;
4374 int nthread, thread;
4378 check_screw_box(state->box);
4382 if (fr->cutoff_scheme == ecutsGROUP)
4387 for (i = 0; i < estNR; i++)
4393 case estX: /* Always present */ break;
4394 case estV: bV = (state->flags & (1<<i)); break;
4395 case est_SDX_NOTSUPPORTED: break;
4396 case estCGP: bCGP = (state->flags & (1<<i)); break;
4399 case estDISRE_INITF:
4400 case estDISRE_RM3TAV:
4401 case estORIRE_INITF:
4403 /* No processing required */
4406 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4411 if (dd->ncg_tot > comm->nalloc_int)
4413 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4414 srenew(comm->buf_int, comm->nalloc_int);
4416 move = comm->buf_int;
4418 /* Clear the count */
4419 for (c = 0; c < dd->ndim*2; c++)
4425 npbcdim = dd->npbcdim;
4427 for (d = 0; (d < DIM); d++)
4429 limitd[d] = dd->comm->cellsize_min[d];
4430 if (d >= npbcdim && dd->ci[d] == 0)
4432 cell_x0[d] = -GMX_FLOAT_MAX;
4436 cell_x0[d] = comm->cell_x0[d];
4438 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4440 cell_x1[d] = GMX_FLOAT_MAX;
4444 cell_x1[d] = comm->cell_x1[d];
4448 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4449 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4453 /* We check after communication if a charge group moved
4454 * more than one cell. Set the pre-comm check limit to float_max.
4456 limit0[d] = -GMX_FLOAT_MAX;
4457 limit1[d] = GMX_FLOAT_MAX;
4461 make_tric_corr_matrix(npbcdim, state->box, tcm);
4463 cgindex = dd->cgindex;
4465 nthread = gmx_omp_nthreads_get(emntDomdec);
4467 /* Compute the center of geometry for all home charge groups
4468 * and put them in the box and determine where they should go.
4470 #pragma omp parallel for num_threads(nthread) schedule(static)
4471 for (thread = 0; thread < nthread; thread++)
4475 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4476 cell_x0, cell_x1, limitd, limit0, limit1,
4478 ( thread *dd->ncg_home)/nthread,
4479 ((thread+1)*dd->ncg_home)/nthread,
4480 fr->cutoff_scheme == ecutsGROUP ? cg_cm : as_rvec_array(state->x.data()),
4483 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
4486 for (cg = 0; cg < dd->ncg_home; cg++)
4491 flag = mc & ~DD_FLAG_NRCG;
4492 mc = mc & DD_FLAG_NRCG;
4495 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4497 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4498 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4500 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4501 /* We store the cg size in the lower 16 bits
4502 * and the place where the charge group should go
4503 * in the next 6 bits. This saves some communication volume.
4505 nrcg = cgindex[cg+1] - cgindex[cg];
4506 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4512 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4513 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4516 for (i = 0; i < dd->ndim*2; i++)
4518 *ncg_moved += ncg[i];
4531 /* Make sure the communication buffers are large enough */
4532 for (mc = 0; mc < dd->ndim*2; mc++)
4534 nvr = ncg[mc] + nat[mc]*nvec;
4535 if (nvr > comm->cgcm_state_nalloc[mc])
4537 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4538 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4542 switch (fr->cutoff_scheme)
4545 /* Recalculating cg_cm might be cheaper than communicating,
4546 * but that could give rise to rounding issues.
4549 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4550 nvec, cg_cm, comm, bCompact);
4553 /* Without charge groups we send the moved atom coordinates
4554 * over twice. This is so the code below can be used without
4555 * many conditionals for both for with and without charge groups.
4558 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4559 nvec, as_rvec_array(state->x.data()), comm, FALSE);
4562 home_pos_cg -= *ncg_moved;
4566 gmx_incons("unimplemented");
4572 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4573 nvec, vec++, as_rvec_array(state->x.data()),
4577 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4578 nvec, vec++, as_rvec_array(state->v.data()),
4583 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4584 nvec, vec++, as_rvec_array(state->cg_p.data()),
4590 compact_ind(dd->ncg_home, move,
4591 dd->index_gl, dd->cgindex, dd->gatindex,
4592 dd->ga2la, comm->bLocalCG,
4597 if (fr->cutoff_scheme == ecutsVERLET)
4599 moved = get_moved(comm, dd->ncg_home);
4601 for (k = 0; k < dd->ncg_home; k++)
4608 moved = fr->ns->grid->cell_index;
4611 clear_and_mark_ind(dd->ncg_home, move,
4612 dd->index_gl, dd->cgindex, dd->gatindex,
4613 dd->ga2la, comm->bLocalCG,
4617 cginfo_mb = fr->cginfo_mb;
4619 *ncg_stay_home = home_pos_cg;
4620 for (d = 0; d < dd->ndim; d++)
4625 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4628 /* Communicate the cg and atom counts */
4633 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4634 d, dir, sbuf[0], sbuf[1]);
4636 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4638 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4640 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4641 srenew(comm->buf_int, comm->nalloc_int);
4644 /* Communicate the charge group indices, sizes and flags */
4645 dd_sendrecv_int(dd, d, dir,
4646 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4647 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4649 nvs = ncg[cdd] + nat[cdd]*nvec;
4650 i = rbuf[0] + rbuf[1] *nvec;
4651 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4653 /* Communicate cgcm and state */
4654 dd_sendrecv_rvec(dd, d, dir,
4655 comm->cgcm_state[cdd], nvs,
4656 comm->vbuf.v+nvr, i);
4657 ncg_recv += rbuf[0];
4661 dd_check_alloc_ncg(fr, state, f, home_pos_cg + ncg_recv);
4662 if (fr->cutoff_scheme == ecutsGROUP)
4664 /* Here we resize to more than necessary and shrink later */
4665 dd_resize_state(state, f, home_pos_at + ncg_recv*MAX_CGCGSIZE);
4668 /* Process the received charge groups */
4670 for (cg = 0; cg < ncg_recv; cg++)
4672 flag = comm->buf_int[cg*DD_CGIBS+1];
4674 if (dim >= npbcdim && dd->nc[dim] > 2)
4676 /* No pbc in this dim and more than one domain boundary.
4677 * We do a separate check if a charge group didn't move too far.
4679 if (((flag & DD_FLAG_FW(d)) &&
4680 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4681 ((flag & DD_FLAG_BW(d)) &&
4682 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4684 cg_move_error(fplog, dd, step, cg, dim,
4685 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4686 fr->cutoff_scheme == ecutsGROUP, 0,
4687 comm->vbuf.v[buf_pos],
4688 comm->vbuf.v[buf_pos],
4689 comm->vbuf.v[buf_pos][dim]);
4696 /* Check which direction this cg should go */
4697 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4699 if (dlbIsOn(dd->comm))
4701 /* The cell boundaries for dimension d2 are not equal
4702 * for each cell row of the lower dimension(s),
4703 * therefore we might need to redetermine where
4704 * this cg should go.
4707 /* If this cg crosses the box boundary in dimension d2
4708 * we can use the communicated flag, so we do not
4709 * have to worry about pbc.
4711 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4712 (flag & DD_FLAG_FW(d2))) ||
4713 (dd->ci[dim2] == 0 &&
4714 (flag & DD_FLAG_BW(d2)))))
4716 /* Clear the two flags for this dimension */
4717 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4718 /* Determine the location of this cg
4719 * in lattice coordinates
4721 pos_d = comm->vbuf.v[buf_pos][dim2];
4724 for (d3 = dim2+1; d3 < DIM; d3++)
4727 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4730 /* Check of we are not at the box edge.
4731 * pbc is only handled in the first step above,
4732 * but this check could move over pbc while
4733 * the first step did not due to different rounding.
4735 if (pos_d >= cell_x1[dim2] &&
4736 dd->ci[dim2] != dd->nc[dim2]-1)
4738 flag |= DD_FLAG_FW(d2);
4740 else if (pos_d < cell_x0[dim2] &&
4743 flag |= DD_FLAG_BW(d2);
4745 comm->buf_int[cg*DD_CGIBS+1] = flag;
4748 /* Set to which neighboring cell this cg should go */
4749 if (flag & DD_FLAG_FW(d2))
4753 else if (flag & DD_FLAG_BW(d2))
4755 if (dd->nc[dd->dim[d2]] > 2)
4767 nrcg = flag & DD_FLAG_NRCG;
4770 if (home_pos_cg+1 > dd->cg_nalloc)
4772 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4773 srenew(dd->index_gl, dd->cg_nalloc);
4774 srenew(dd->cgindex, dd->cg_nalloc+1);
4776 /* Set the global charge group index and size */
4777 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4778 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4779 /* Copy the state from the buffer */
4780 if (fr->cutoff_scheme == ecutsGROUP)
4783 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4787 /* Set the cginfo */
4788 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4789 dd->index_gl[home_pos_cg]);
4792 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4795 for (i = 0; i < nrcg; i++)
4797 copy_rvec(comm->vbuf.v[buf_pos++],
4798 state->x[home_pos_at+i]);
4802 for (i = 0; i < nrcg; i++)
4804 copy_rvec(comm->vbuf.v[buf_pos++],
4805 state->v[home_pos_at+i]);
4810 for (i = 0; i < nrcg; i++)
4812 copy_rvec(comm->vbuf.v[buf_pos++],
4813 state->cg_p[home_pos_at+i]);
4817 home_pos_at += nrcg;
4821 /* Reallocate the buffers if necessary */
4822 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4824 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4825 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4827 nvr = ncg[mc] + nat[mc]*nvec;
4828 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
4830 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
4831 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4833 /* Copy from the receive to the send buffers */
4834 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
4835 comm->buf_int + cg*DD_CGIBS,
4836 DD_CGIBS*sizeof(int));
4837 memcpy(comm->cgcm_state[mc][nvr],
4838 comm->vbuf.v[buf_pos],
4839 (1+nrcg*nvec)*sizeof(rvec));
4840 buf_pos += 1 + nrcg*nvec;
4847 /* With sorting (!bCompact) the indices are now only partially up to date
4848 * and ncg_home and nat_home are not the real count, since there are
4849 * "holes" in the arrays for the charge groups that moved to neighbors.
4851 if (fr->cutoff_scheme == ecutsVERLET)
4853 moved = get_moved(comm, home_pos_cg);
4855 for (i = dd->ncg_home; i < home_pos_cg; i++)
4860 dd->ncg_home = home_pos_cg;
4861 dd->nat_home = home_pos_at;
4863 if (fr->cutoff_scheme == ecutsGROUP && !bCompact)
4865 /* We overallocated before, we need to set the right size here */
4866 dd_resize_state(state, f, dd->nat_home);
4872 "Finished repartitioning: cgs moved out %d, new home %d\n",
4873 *ncg_moved, dd->ncg_home-*ncg_moved);
4878 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
4880 /* Note that the cycles value can be incorrect, either 0 or some
4881 * extremely large value, when our thread migrated to another core
4882 * with an unsynchronized cycle counter. If this happens less often
4883 * that once per nstlist steps, this will not cause issues, since
4884 * we later subtract the maximum value from the sum over nstlist steps.
4885 * A zero count will slightly lower the total, but that's a small effect.
4886 * Note that the main purpose of the subtraction of the maximum value
4887 * is to avoid throwing off the load balancing when stalls occur due
4888 * e.g. system activity or network congestion.
4890 dd->comm->cycl[ddCycl] += cycles;
4891 dd->comm->cycl_n[ddCycl]++;
4892 if (cycles > dd->comm->cycl_max[ddCycl])
4894 dd->comm->cycl_max[ddCycl] = cycles;
4898 static double force_flop_count(t_nrnb *nrnb)
4905 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
4907 /* To get closer to the real timings, we half the count
4908 * for the normal loops and again half it for water loops.
4911 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
4913 sum += nrnb->n[i]*0.25*cost_nrnb(i);
4917 sum += nrnb->n[i]*0.50*cost_nrnb(i);
4920 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
4923 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
4925 sum += nrnb->n[i]*cost_nrnb(i);
4928 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
4930 sum += nrnb->n[i]*cost_nrnb(i);
4936 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
4938 if (dd->comm->eFlop)
4940 dd->comm->flop -= force_flop_count(nrnb);
4943 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
4945 if (dd->comm->eFlop)
4947 dd->comm->flop += force_flop_count(nrnb);
4952 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
4956 for (i = 0; i < ddCyclNr; i++)
4958 dd->comm->cycl[i] = 0;
4959 dd->comm->cycl_n[i] = 0;
4960 dd->comm->cycl_max[i] = 0;
4963 dd->comm->flop_n = 0;
4966 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
4968 gmx_domdec_comm_t *comm;
4969 domdec_load_t *load;
4970 domdec_root_t *root = NULL;
4972 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
4977 fprintf(debug, "get_load_distribution start\n");
4980 wallcycle_start(wcycle, ewcDDCOMMLOAD);
4984 bSepPME = (dd->pme_nodeid >= 0);
4986 if (dd->ndim == 0 && bSepPME)
4988 /* Without decomposition, but with PME nodes, we need the load */
4989 comm->load[0].mdf = comm->cycl[ddCyclPPduringPME];
4990 comm->load[0].pme = comm->cycl[ddCyclPME];
4993 for (d = dd->ndim-1; d >= 0; d--)
4996 /* Check if we participate in the communication in this dimension */
4997 if (d == dd->ndim-1 ||
4998 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5000 load = &comm->load[d];
5001 if (dlbIsOn(dd->comm))
5003 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5006 if (d == dd->ndim-1)
5008 sbuf[pos++] = dd_force_load(comm);
5009 sbuf[pos++] = sbuf[0];
5010 if (dlbIsOn(dd->comm))
5012 sbuf[pos++] = sbuf[0];
5013 sbuf[pos++] = cell_frac;
5016 sbuf[pos++] = comm->cell_f_max0[d];
5017 sbuf[pos++] = comm->cell_f_min1[d];
5022 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5023 sbuf[pos++] = comm->cycl[ddCyclPME];
5028 sbuf[pos++] = comm->load[d+1].sum;
5029 sbuf[pos++] = comm->load[d+1].max;
5030 if (dlbIsOn(dd->comm))
5032 sbuf[pos++] = comm->load[d+1].sum_m;
5033 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5034 sbuf[pos++] = comm->load[d+1].flags;
5037 sbuf[pos++] = comm->cell_f_max0[d];
5038 sbuf[pos++] = comm->cell_f_min1[d];
5043 sbuf[pos++] = comm->load[d+1].mdf;
5044 sbuf[pos++] = comm->load[d+1].pme;
5048 /* Communicate a row in DD direction d.
5049 * The communicators are setup such that the root always has rank 0.
5052 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5053 load->load, load->nload*sizeof(float), MPI_BYTE,
5054 0, comm->mpi_comm_load[d]);
5056 if (dd->ci[dim] == dd->master_ci[dim])
5058 /* We are the root, process this row */
5061 root = comm->root[d];
5071 for (i = 0; i < dd->nc[dim]; i++)
5073 load->sum += load->load[pos++];
5074 load->max = std::max(load->max, load->load[pos]);
5076 if (dlbIsOn(dd->comm))
5080 /* This direction could not be load balanced properly,
5081 * therefore we need to use the maximum iso the average load.
5083 load->sum_m = std::max(load->sum_m, load->load[pos]);
5087 load->sum_m += load->load[pos];
5090 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5094 load->flags = (int)(load->load[pos++] + 0.5);
5098 root->cell_f_max0[i] = load->load[pos++];
5099 root->cell_f_min1[i] = load->load[pos++];
5104 load->mdf = std::max(load->mdf, load->load[pos]);
5106 load->pme = std::max(load->pme, load->load[pos]);
5110 if (dlbIsOn(comm) && root->bLimited)
5112 load->sum_m *= dd->nc[dim];
5113 load->flags |= (1<<d);
5121 comm->nload += dd_load_count(comm);
5122 comm->load_step += comm->cycl[ddCyclStep];
5123 comm->load_sum += comm->load[0].sum;
5124 comm->load_max += comm->load[0].max;
5127 for (d = 0; d < dd->ndim; d++)
5129 if (comm->load[0].flags & (1<<d))
5131 comm->load_lim[d]++;
5137 comm->load_mdf += comm->load[0].mdf;
5138 comm->load_pme += comm->load[0].pme;
5142 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5146 fprintf(debug, "get_load_distribution finished\n");
5150 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5152 /* Return the relative performance loss on the total run time
5153 * due to the force calculation load imbalance.
5155 if (dd->comm->nload > 0 && dd->comm->load_step > 0)
5158 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5159 (dd->comm->load_step*dd->nnodes);
5167 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5170 int npp, npme, nnodes, d, limp;
5171 float imbal, pme_f_ratio, lossf = 0, lossp = 0;
5173 gmx_domdec_comm_t *comm;
5176 if (DDMASTER(dd) && comm->nload > 0)
5179 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5180 nnodes = npp + npme;
5181 if (dd->nnodes > 1 && comm->load_sum > 0)
5183 imbal = comm->load_max*npp/comm->load_sum - 1;
5184 lossf = dd_force_imb_perf_loss(dd);
5185 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5186 fprintf(fplog, "%s", buf);
5187 fprintf(stderr, "\n");
5188 fprintf(stderr, "%s", buf);
5189 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5190 fprintf(fplog, "%s", buf);
5191 fprintf(stderr, "%s", buf);
5196 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5197 for (d = 0; d < dd->ndim; d++)
5199 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5200 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5206 sprintf(buf+strlen(buf), "\n");
5207 fprintf(fplog, "%s", buf);
5208 fprintf(stderr, "%s", buf);
5210 if (npme > 0 && comm->load_mdf > 0 && comm->load_step > 0)
5212 pme_f_ratio = comm->load_pme/comm->load_mdf;
5213 lossp = (comm->load_pme - comm->load_mdf)/comm->load_step;
5216 lossp *= (float)npme/(float)nnodes;
5220 lossp *= (float)npp/(float)nnodes;
5222 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5223 fprintf(fplog, "%s", buf);
5224 fprintf(stderr, "%s", buf);
5225 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5226 fprintf(fplog, "%s", buf);
5227 fprintf(stderr, "%s", buf);
5229 fprintf(fplog, "\n");
5230 fprintf(stderr, "\n");
5232 if (lossf >= DD_PERF_LOSS_WARN)
5235 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5236 " in the domain decomposition.\n", lossf*100);
5239 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5243 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5245 fprintf(fplog, "%s\n", buf);
5246 fprintf(stderr, "%s\n", buf);
5248 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5251 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5252 " had %s work to do than the PP ranks.\n"
5253 " You might want to %s the number of PME ranks\n"
5254 " or %s the cut-off and the grid spacing.\n",
5256 (lossp < 0) ? "less" : "more",
5257 (lossp < 0) ? "decrease" : "increase",
5258 (lossp < 0) ? "decrease" : "increase");
5259 fprintf(fplog, "%s\n", buf);
5260 fprintf(stderr, "%s\n", buf);
5265 static float dd_vol_min(gmx_domdec_t *dd)
5267 return dd->comm->load[0].cvol_min*dd->nnodes;
5270 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5272 return dd->comm->load[0].flags;
5275 static float dd_f_imbal(gmx_domdec_t *dd)
5277 if (dd->comm->load[0].sum > 0)
5279 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1.0f;
5283 /* Something is wrong in the cycle counting, report no load imbalance */
5288 float dd_pme_f_ratio(gmx_domdec_t *dd)
5290 /* Should only be called on the DD master rank */
5291 assert(DDMASTER(dd));
5293 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5295 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5303 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5308 flags = dd_load_flags(dd);
5312 "DD load balancing is limited by minimum cell size in dimension");
5313 for (d = 0; d < dd->ndim; d++)
5317 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5320 fprintf(fplog, "\n");
5322 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5323 if (dlbIsOn(dd->comm))
5325 fprintf(fplog, " vol min/aver %5.3f%c",
5326 dd_vol_min(dd), flags ? '!' : ' ');
5330 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5332 if (dd->comm->cycl_n[ddCyclPME])
5334 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5336 fprintf(fplog, "\n\n");
5339 static void dd_print_load_verbose(gmx_domdec_t *dd)
5341 if (dlbIsOn(dd->comm))
5343 fprintf(stderr, "vol %4.2f%c ",
5344 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5348 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5350 if (dd->comm->cycl_n[ddCyclPME])
5352 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5357 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5362 domdec_root_t *root;
5363 gmx_bool bPartOfGroup = FALSE;
5365 dim = dd->dim[dim_ind];
5366 copy_ivec(loc, loc_c);
5367 for (i = 0; i < dd->nc[dim]; i++)
5370 rank = dd_index(dd->nc, loc_c);
5371 if (rank == dd->rank)
5373 /* This process is part of the group */
5374 bPartOfGroup = TRUE;
5377 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5381 dd->comm->mpi_comm_load[dim_ind] = c_row;
5382 if (dd->comm->dlbState != edlbsOffForever)
5384 if (dd->ci[dim] == dd->master_ci[dim])
5386 /* This is the root process of this row */
5387 snew(dd->comm->root[dim_ind], 1);
5388 root = dd->comm->root[dim_ind];
5389 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5390 snew(root->old_cell_f, dd->nc[dim]+1);
5391 snew(root->bCellMin, dd->nc[dim]);
5394 snew(root->cell_f_max0, dd->nc[dim]);
5395 snew(root->cell_f_min1, dd->nc[dim]);
5396 snew(root->bound_min, dd->nc[dim]);
5397 snew(root->bound_max, dd->nc[dim]);
5399 snew(root->buf_ncd, dd->nc[dim]);
5403 /* This is not a root process, we only need to receive cell_f */
5404 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5407 if (dd->ci[dim] == dd->master_ci[dim])
5409 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5415 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5416 const gmx_hw_info_t gmx_unused *hwinfo,
5417 const gmx_hw_opt_t gmx_unused *hw_opt)
5420 int physicalnode_id_hash;
5423 MPI_Comm mpi_comm_pp_physicalnode;
5425 if (!(cr->duty & DUTY_PP) || hw_opt->gpu_opt.n_dev_use == 0)
5427 /* Only PP nodes (currently) use GPUs.
5428 * If we don't have GPUs, there are no resources to share.
5433 physicalnode_id_hash = gmx_physicalnode_id_hash();
5435 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5441 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5442 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5443 dd->rank, physicalnode_id_hash, gpu_id);
5445 /* Split the PP communicator over the physical nodes */
5446 /* TODO: See if we should store this (before), as it's also used for
5447 * for the nodecomm summution.
5449 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5450 &mpi_comm_pp_physicalnode);
5451 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5452 &dd->comm->mpi_comm_gpu_shared);
5453 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5454 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5458 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5461 /* Note that some ranks could share a GPU, while others don't */
5463 if (dd->comm->nrank_gpu_shared == 1)
5465 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5470 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5473 int dim0, dim1, i, j;
5478 fprintf(debug, "Making load communicators\n");
5481 snew(dd->comm->load, std::max(dd->ndim, 1));
5482 snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
5490 make_load_communicator(dd, 0, loc);
5494 for (i = 0; i < dd->nc[dim0]; i++)
5497 make_load_communicator(dd, 1, loc);
5503 for (i = 0; i < dd->nc[dim0]; i++)
5507 for (j = 0; j < dd->nc[dim1]; j++)
5510 make_load_communicator(dd, 2, loc);
5517 fprintf(debug, "Finished making load communicators\n");
5522 /*! \brief Sets up the relation between neighboring domains and zones */
5523 static void setup_neighbor_relations(gmx_domdec_t *dd)
5525 int d, dim, i, j, m;
5527 gmx_domdec_zones_t *zones;
5528 gmx_domdec_ns_ranges_t *izone;
5530 for (d = 0; d < dd->ndim; d++)
5533 copy_ivec(dd->ci, tmp);
5534 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5535 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5536 copy_ivec(dd->ci, tmp);
5537 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5538 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5541 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5544 dd->neighbor[d][1]);
5548 int nzone = (1 << dd->ndim);
5549 int nizone = (1 << std::max(dd->ndim - 1, 0));
5550 assert(nizone >= 1 && nizone <= DD_MAXIZONE);
5552 zones = &dd->comm->zones;
5554 for (i = 0; i < nzone; i++)
5557 clear_ivec(zones->shift[i]);
5558 for (d = 0; d < dd->ndim; d++)
5560 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5565 for (i = 0; i < nzone; i++)
5567 for (d = 0; d < DIM; d++)
5569 s[d] = dd->ci[d] - zones->shift[i][d];
5574 else if (s[d] >= dd->nc[d])
5580 zones->nizone = nizone;
5581 for (i = 0; i < zones->nizone; i++)
5583 assert(ddNonbondedZonePairRanges[i][0] == i);
5585 izone = &zones->izone[i];
5586 /* dd_zp3 is for 3D decomposition, for fewer dimensions use only
5587 * j-zones up to nzone.
5589 izone->j0 = std::min(ddNonbondedZonePairRanges[i][1], nzone);
5590 izone->j1 = std::min(ddNonbondedZonePairRanges[i][2], nzone);
5591 for (dim = 0; dim < DIM; dim++)
5593 if (dd->nc[dim] == 1)
5595 /* All shifts should be allowed */
5596 izone->shift0[dim] = -1;
5597 izone->shift1[dim] = 1;
5601 /* Determine the min/max j-zone shift wrt the i-zone */
5602 izone->shift0[dim] = 1;
5603 izone->shift1[dim] = -1;
5604 for (j = izone->j0; j < izone->j1; j++)
5606 int shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5607 if (shift_diff < izone->shift0[dim])
5609 izone->shift0[dim] = shift_diff;
5611 if (shift_diff > izone->shift1[dim])
5613 izone->shift1[dim] = shift_diff;
5620 if (dd->comm->dlbState != edlbsOffForever)
5622 snew(dd->comm->root, dd->ndim);
5625 if (dd->comm->bRecordLoad)
5627 make_load_communicators(dd);
5631 static void make_pp_communicator(FILE *fplog,
5633 t_commrec gmx_unused *cr,
5634 int gmx_unused reorder)
5637 gmx_domdec_comm_t *comm;
5644 if (comm->bCartesianPP)
5646 /* Set up cartesian communication for the particle-particle part */
5649 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5650 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5653 for (int i = 0; i < DIM; i++)
5657 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5659 /* We overwrite the old communicator with the new cartesian one */
5660 cr->mpi_comm_mygroup = comm_cart;
5663 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5664 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5666 if (comm->bCartesianPP_PME)
5668 /* Since we want to use the original cartesian setup for sim,
5669 * and not the one after split, we need to make an index.
5671 snew(comm->ddindex2ddnodeid, dd->nnodes);
5672 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5673 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5674 /* Get the rank of the DD master,
5675 * above we made sure that the master node is a PP node.
5685 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5687 else if (comm->bCartesianPP)
5689 if (cr->npmenodes == 0)
5691 /* The PP communicator is also
5692 * the communicator for this simulation
5694 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5696 cr->nodeid = dd->rank;
5698 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5700 /* We need to make an index to go from the coordinates
5701 * to the nodeid of this simulation.
5703 snew(comm->ddindex2simnodeid, dd->nnodes);
5704 snew(buf, dd->nnodes);
5705 if (cr->duty & DUTY_PP)
5707 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5709 /* Communicate the ddindex to simulation nodeid index */
5710 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5711 cr->mpi_comm_mysim);
5714 /* Determine the master coordinates and rank.
5715 * The DD master should be the same node as the master of this sim.
5717 for (int i = 0; i < dd->nnodes; i++)
5719 if (comm->ddindex2simnodeid[i] == 0)
5721 ddindex2xyz(dd->nc, i, dd->master_ci);
5722 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5727 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5732 /* No Cartesian communicators */
5733 /* We use the rank in dd->comm->all as DD index */
5734 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5735 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5737 clear_ivec(dd->master_ci);
5744 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
5745 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5750 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
5751 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5755 static void receive_ddindex2simnodeid(gmx_domdec_t *dd,
5759 gmx_domdec_comm_t *comm = dd->comm;
5761 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5764 snew(comm->ddindex2simnodeid, dd->nnodes);
5765 snew(buf, dd->nnodes);
5766 if (cr->duty & DUTY_PP)
5768 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5770 /* Communicate the ddindex to simulation nodeid index */
5771 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5772 cr->mpi_comm_mysim);
5776 GMX_UNUSED_VALUE(dd);
5777 GMX_UNUSED_VALUE(cr);
5781 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5782 int ncg, int natoms)
5784 gmx_domdec_master_t *ma;
5789 snew(ma->ncg, dd->nnodes);
5790 snew(ma->index, dd->nnodes+1);
5792 snew(ma->nat, dd->nnodes);
5793 snew(ma->ibuf, dd->nnodes*2);
5794 snew(ma->cell_x, DIM);
5795 for (i = 0; i < DIM; i++)
5797 snew(ma->cell_x[i], dd->nc[i]+1);
5800 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5806 snew(ma->vbuf, natoms);
5812 static void split_communicator(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd,
5813 int gmx_unused dd_rank_order,
5814 int gmx_unused reorder)
5816 gmx_domdec_comm_t *comm;
5825 if (comm->bCartesianPP)
5827 for (i = 1; i < DIM; i++)
5829 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5831 if (bDiv[YY] || bDiv[ZZ])
5833 comm->bCartesianPP_PME = TRUE;
5834 /* If we have 2D PME decomposition, which is always in x+y,
5835 * we stack the PME only nodes in z.
5836 * Otherwise we choose the direction that provides the thinnest slab
5837 * of PME only nodes as this will have the least effect
5838 * on the PP communication.
5839 * But for the PME communication the opposite might be better.
5841 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
5843 dd->nc[YY] > dd->nc[ZZ]))
5845 comm->cartpmedim = ZZ;
5849 comm->cartpmedim = YY;
5851 comm->ntot[comm->cartpmedim]
5852 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
5856 fprintf(fplog, "Number of PME-only ranks (%d) is not a multiple of nx*ny (%d*%d) or nx*nz (%d*%d)\n", cr->npmenodes, dd->nc[XX], dd->nc[YY], dd->nc[XX], dd->nc[ZZ]);
5858 "Will not use a Cartesian communicator for PP <-> PME\n\n");
5863 if (comm->bCartesianPP_PME)
5870 fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
5873 for (i = 0; i < DIM; i++)
5877 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
5879 MPI_Comm_rank(comm_cart, &rank);
5880 if (MASTER(cr) && rank != 0)
5882 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
5885 /* With this assigment we loose the link to the original communicator
5886 * which will usually be MPI_COMM_WORLD, unless have multisim.
5888 cr->mpi_comm_mysim = comm_cart;
5889 cr->sim_nodeid = rank;
5891 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
5895 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
5896 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5899 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
5903 if (cr->npmenodes == 0 ||
5904 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
5906 cr->duty = DUTY_PME;
5909 /* Split the sim communicator into PP and PME only nodes */
5910 MPI_Comm_split(cr->mpi_comm_mysim,
5912 dd_index(comm->ntot, dd->ci),
5913 &cr->mpi_comm_mygroup);
5917 switch (dd_rank_order)
5919 case ddrankorderPP_PME:
5922 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
5925 case ddrankorderINTERLEAVE:
5926 /* Interleave the PP-only and PME-only ranks */
5929 fprintf(fplog, "Interleaving PP and PME ranks\n");
5931 comm->pmenodes = dd_interleaved_pme_ranks(dd);
5933 case ddrankorderCARTESIAN:
5936 gmx_fatal(FARGS, "Unknown dd_rank_order=%d", dd_rank_order);
5939 if (dd_simnode2pmenode(dd, cr, cr->sim_nodeid) == -1)
5941 cr->duty = DUTY_PME;
5948 /* Split the sim communicator into PP and PME only nodes */
5949 MPI_Comm_split(cr->mpi_comm_mysim,
5952 &cr->mpi_comm_mygroup);
5953 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
5959 fprintf(fplog, "This rank does only %s work.\n\n",
5960 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
5964 /*! \brief Generates the MPI communicators for domain decomposition */
5965 static void make_dd_communicators(FILE *fplog, t_commrec *cr,
5966 gmx_domdec_t *dd, int dd_rank_order)
5968 gmx_domdec_comm_t *comm;
5973 copy_ivec(dd->nc, comm->ntot);
5975 comm->bCartesianPP = (dd_rank_order == ddrankorderCARTESIAN);
5976 comm->bCartesianPP_PME = FALSE;
5978 /* Reorder the nodes by default. This might change the MPI ranks.
5979 * Real reordering is only supported on very few architectures,
5980 * Blue Gene is one of them.
5982 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
5984 if (cr->npmenodes > 0)
5986 /* Split the communicator into a PP and PME part */
5987 split_communicator(fplog, cr, dd, dd_rank_order, CartReorder);
5988 if (comm->bCartesianPP_PME)
5990 /* We (possibly) reordered the nodes in split_communicator,
5991 * so it is no longer required in make_pp_communicator.
5993 CartReorder = FALSE;
5998 /* All nodes do PP and PME */
6000 /* We do not require separate communicators */
6001 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6005 if (cr->duty & DUTY_PP)
6007 /* Copy or make a new PP communicator */
6008 make_pp_communicator(fplog, dd, cr, CartReorder);
6012 receive_ddindex2simnodeid(dd, cr);
6015 if (!(cr->duty & DUTY_PME))
6017 /* Set up the commnuication to our PME node */
6018 dd->pme_nodeid = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
6019 dd->pme_receive_vir_ener = receive_vir_ener(dd, cr);
6022 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6023 dd->pme_nodeid, dd->pme_receive_vir_ener);
6028 dd->pme_nodeid = -1;
6033 dd->ma = init_gmx_domdec_master_t(dd,
6035 comm->cgs_gl.index[comm->cgs_gl.nr]);
6039 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6041 real *slb_frac, tot;
6046 if (nc > 1 && size_string != NULL)
6050 fprintf(fplog, "Using static load balancing for the %s direction\n",
6055 for (i = 0; i < nc; i++)
6058 sscanf(size_string, "%20lf%n", &dbl, &n);
6061 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6070 fprintf(fplog, "Relative cell sizes:");
6072 for (i = 0; i < nc; i++)
6077 fprintf(fplog, " %5.3f", slb_frac[i]);
6082 fprintf(fplog, "\n");
6089 static int multi_body_bondeds_count(const gmx_mtop_t *mtop)
6092 gmx_mtop_ilistloop_t iloop;
6096 iloop = gmx_mtop_ilistloop_init(mtop);
6097 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6099 for (ftype = 0; ftype < F_NRE; ftype++)
6101 if ((interaction_function[ftype].flags & IF_BOND) &&
6104 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6112 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6118 val = getenv(env_var);
6121 if (sscanf(val, "%20d", &nst) <= 0)
6127 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6135 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6139 fprintf(stderr, "\n%s\n", warn_string);
6143 fprintf(fplog, "\n%s\n", warn_string);
6147 static void check_dd_restrictions(t_commrec *cr, const gmx_domdec_t *dd,
6148 const t_inputrec *ir, FILE *fplog)
6150 if (ir->ePBC == epbcSCREW &&
6151 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6153 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6156 if (ir->ns_type == ensSIMPLE)
6158 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6161 if (ir->nstlist == 0)
6163 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6166 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6168 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6172 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6177 r = ddbox->box_size[XX];
6178 for (di = 0; di < dd->ndim; di++)
6181 /* Check using the initial average cell size */
6182 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6188 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6189 const char *dlb_opt, gmx_bool bRecordLoad,
6190 unsigned long Flags, const t_inputrec *ir)
6197 case 'a': dlbState = edlbsOffCanTurnOn; break;
6198 case 'n': dlbState = edlbsOffForever; break;
6199 case 'y': dlbState = edlbsOn; break;
6200 default: gmx_incons("Unknown dlb_opt");
6203 if (Flags & MD_RERUN)
6205 return edlbsOffForever;
6208 if (!EI_DYNAMICS(ir->eI))
6210 if (dlbState == edlbsOn)
6212 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6213 dd_warning(cr, fplog, buf);
6216 return edlbsOffForever;
6221 dd_warning(cr, fplog, "NOTE: Cycle counters unsupported or not enabled in kernel. Cannot use dynamic load balancing.\n");
6222 return edlbsOffForever;
6225 if (Flags & MD_REPRODUCIBLE)
6229 case edlbsOffForever:
6231 case edlbsOffCanTurnOn:
6232 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6233 dlbState = edlbsOffForever;
6236 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6239 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", dlbState);
6247 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6252 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6254 /* Decomposition order z,y,x */
6257 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6259 for (dim = DIM-1; dim >= 0; dim--)
6261 if (dd->nc[dim] > 1)
6263 dd->dim[dd->ndim++] = dim;
6269 /* Decomposition order x,y,z */
6270 for (dim = 0; dim < DIM; dim++)
6272 if (dd->nc[dim] > 1)
6274 dd->dim[dd->ndim++] = dim;
6280 static gmx_domdec_comm_t *init_dd_comm()
6282 gmx_domdec_comm_t *comm;
6286 snew(comm->cggl_flag, DIM*2);
6287 snew(comm->cgcm_state, DIM*2);
6288 for (i = 0; i < DIM*2; i++)
6290 comm->cggl_flag_nalloc[i] = 0;
6291 comm->cgcm_state_nalloc[i] = 0;
6294 comm->nalloc_int = 0;
6295 comm->buf_int = NULL;
6297 vec_rvec_init(&comm->vbuf);
6299 comm->n_load_have = 0;
6300 comm->n_load_collect = 0;
6302 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6304 comm->sum_nat[i] = 0;
6308 comm->load_step = 0;
6311 clear_ivec(comm->load_lim);
6318 /*! \brief Set the cell size and interaction limits, as well as the DD grid */
6319 static void set_dd_limits_and_grid(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd,
6320 unsigned long Flags,
6321 ivec nc, int nPmeRanks,
6322 real comm_distance_min, real rconstr,
6323 const char *dlb_opt, real dlb_scale,
6324 const char *sizex, const char *sizey, const char *sizez,
6325 const gmx_mtop_t *mtop,
6326 const t_inputrec *ir,
6327 matrix box, const rvec *x,
6329 int *npme_x, int *npme_y)
6332 real r_bonded_limit = -1;
6333 const real tenPercentMargin = 1.1;
6334 gmx_domdec_comm_t *comm = dd->comm;
6336 snew(comm->cggl_flag, DIM*2);
6337 snew(comm->cgcm_state, DIM*2);
6339 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6340 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6342 dd->pme_recv_f_alloc = 0;
6343 dd->pme_recv_f_buf = NULL;
6345 /* Initialize to GPU share count to 0, might change later */
6346 comm->nrank_gpu_shared = 0;
6348 comm->dlbState = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6349 comm->bCheckWhetherToTurnDlbOn = TRUE;
6353 fprintf(fplog, "Dynamic load balancing: %s\n",
6354 edlbs_names[comm->dlbState]);
6356 comm->bPMELoadBalDLBLimits = FALSE;
6358 /* Allocate the charge group/atom sorting struct */
6359 snew(comm->sort, 1);
6361 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6363 comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
6364 mtop->bIntermolecularInteractions);
6365 if (comm->bInterCGBondeds)
6367 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6371 comm->bInterCGMultiBody = FALSE;
6374 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6375 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6379 /* Set the cut-off to some very large value,
6380 * so we don't need if statements everywhere in the code.
6381 * We use sqrt, since the cut-off is squared in some places.
6383 comm->cutoff = GMX_CUTOFF_INF;
6387 comm->cutoff = ir->rlist;
6389 comm->cutoff_mbody = 0;
6391 comm->cellsize_limit = 0;
6392 comm->bBondComm = FALSE;
6394 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6395 * within nstlist steps. Since boundaries are allowed to displace by half
6396 * a cell size, DD cells should be at least the size of the list buffer.
6398 comm->cellsize_limit = std::max(comm->cellsize_limit,
6399 ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
6401 if (comm->bInterCGBondeds)
6403 if (comm_distance_min > 0)
6405 comm->cutoff_mbody = comm_distance_min;
6406 if (Flags & MD_DDBONDCOMM)
6408 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6412 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6414 r_bonded_limit = comm->cutoff_mbody;
6416 else if (ir->bPeriodicMols)
6418 /* Can not easily determine the required cut-off */
6419 dd_warning(cr, fplog, "NOTE: Periodic molecules are present in this system. Because of this, the domain decomposition algorithm cannot easily determine the minimum cell size that it requires for treating bonded interactions. Instead, domain decomposition will assume that half the non-bonded cut-off will be a suitable lower bound.\n");
6420 comm->cutoff_mbody = comm->cutoff/2;
6421 r_bonded_limit = comm->cutoff_mbody;
6429 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6430 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6432 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6433 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6435 /* We use an initial margin of 10% for the minimum cell size,
6436 * except when we are just below the non-bonded cut-off.
6438 if (Flags & MD_DDBONDCOMM)
6440 if (std::max(r_2b, r_mb) > comm->cutoff)
6442 r_bonded = std::max(r_2b, r_mb);
6443 r_bonded_limit = tenPercentMargin*r_bonded;
6444 comm->bBondComm = TRUE;
6449 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6451 /* We determine cutoff_mbody later */
6455 /* No special bonded communication,
6456 * simply increase the DD cut-off.
6458 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6459 comm->cutoff_mbody = r_bonded_limit;
6460 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6466 "Minimum cell size due to bonded interactions: %.3f nm\n",
6469 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6472 if (dd->bInterCGcons && rconstr <= 0)
6474 /* There is a cell size limit due to the constraints (P-LINCS) */
6475 rconstr = constr_r_max(fplog, mtop, ir);
6479 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6481 if (rconstr > comm->cellsize_limit)
6483 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6487 else if (rconstr > 0 && fplog)
6489 /* Here we do not check for dd->bInterCGcons,
6490 * because one can also set a cell size limit for virtual sites only
6491 * and at this point we don't know yet if there are intercg v-sites.
6494 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6497 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6499 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6503 copy_ivec(nc, dd->nc);
6504 set_dd_dim(fplog, dd);
6505 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6509 cr->npmenodes = nPmeRanks;
6513 /* When the DD grid is set explicitly and -npme is set to auto,
6514 * don't use PME ranks. We check later if the DD grid is
6515 * compatible with the total number of ranks.
6520 real acs = average_cellsize_min(dd, ddbox);
6521 if (acs < comm->cellsize_limit)
6525 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6527 gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
6528 "The initial cell size (%f) is smaller than the cell size limit (%f), change options -dd, -rdd or -rcon, see the log file for details",
6529 acs, comm->cellsize_limit);
6534 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6536 /* We need to choose the optimal DD grid and possibly PME nodes */
6538 dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6540 comm->dlbState != edlbsOffForever, dlb_scale,
6541 comm->cellsize_limit, comm->cutoff,
6542 comm->bInterCGBondeds);
6544 if (dd->nc[XX] == 0)
6547 gmx_bool bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6548 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6549 !bC ? "-rdd" : "-rcon",
6550 comm->dlbState != edlbsOffForever ? " or -dds" : "",
6551 bC ? " or your LINCS settings" : "");
6553 gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
6554 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6556 "Look in the log file for details on the domain decomposition",
6557 cr->nnodes-cr->npmenodes, limit, buf);
6559 set_dd_dim(fplog, dd);
6565 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6566 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6569 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6570 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6572 gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
6573 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6574 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6576 if (cr->npmenodes > dd->nnodes)
6578 gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
6579 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6581 if (cr->npmenodes > 0)
6583 comm->npmenodes = cr->npmenodes;
6587 comm->npmenodes = dd->nnodes;
6590 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6592 /* The following choices should match those
6593 * in comm_cost_est in domdec_setup.c.
6594 * Note that here the checks have to take into account
6595 * that the decomposition might occur in a different order than xyz
6596 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6597 * in which case they will not match those in comm_cost_est,
6598 * but since that is mainly for testing purposes that's fine.
6600 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6601 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6602 getenv("GMX_PMEONEDD") == NULL)
6604 comm->npmedecompdim = 2;
6605 comm->npmenodes_x = dd->nc[XX];
6606 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6610 /* In case nc is 1 in both x and y we could still choose to
6611 * decompose pme in y instead of x, but we use x for simplicity.
6613 comm->npmedecompdim = 1;
6614 if (dd->dim[0] == YY)
6616 comm->npmenodes_x = 1;
6617 comm->npmenodes_y = comm->npmenodes;
6621 comm->npmenodes_x = comm->npmenodes;
6622 comm->npmenodes_y = 1;
6627 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6628 comm->npmenodes_x, comm->npmenodes_y, 1);
6633 comm->npmedecompdim = 0;
6634 comm->npmenodes_x = 0;
6635 comm->npmenodes_y = 0;
6638 /* Technically we don't need both of these,
6639 * but it simplifies code not having to recalculate it.
6641 *npme_x = comm->npmenodes_x;
6642 *npme_y = comm->npmenodes_y;
6644 snew(comm->slb_frac, DIM);
6645 if (comm->dlbState == edlbsOffForever)
6647 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6648 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6649 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6652 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6654 if (comm->bBondComm || comm->dlbState != edlbsOffForever)
6656 /* Set the bonded communication distance to halfway
6657 * the minimum and the maximum,
6658 * since the extra communication cost is nearly zero.
6660 real acs = average_cellsize_min(dd, ddbox);
6661 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6662 if (comm->dlbState != edlbsOffForever)
6664 /* Check if this does not limit the scaling */
6665 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
6667 if (!comm->bBondComm)
6669 /* Without bBondComm do not go beyond the n.b. cut-off */
6670 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
6671 if (comm->cellsize_limit >= comm->cutoff)
6673 /* We don't loose a lot of efficieny
6674 * when increasing it to the n.b. cut-off.
6675 * It can even be slightly faster, because we need
6676 * less checks for the communication setup.
6678 comm->cutoff_mbody = comm->cutoff;
6681 /* Check if we did not end up below our original limit */
6682 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
6684 if (comm->cutoff_mbody > comm->cellsize_limit)
6686 comm->cellsize_limit = comm->cutoff_mbody;
6689 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6694 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6695 "cellsize limit %f\n",
6696 comm->bBondComm, comm->cellsize_limit);
6701 check_dd_restrictions(cr, dd, ir, fplog);
6705 static void set_dlb_limits(gmx_domdec_t *dd)
6710 for (d = 0; d < dd->ndim; d++)
6712 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6713 dd->comm->cellsize_min[dd->dim[d]] =
6714 dd->comm->cellsize_min_dlb[dd->dim[d]];
6719 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
6722 gmx_domdec_comm_t *comm;
6732 fprintf(fplog, "At step %s the performance loss due to force load imbalance is %.1f %%\n", gmx_step_str(step, buf), dd_force_imb_perf_loss(dd)*100);
6735 cellsize_min = comm->cellsize_min[dd->dim[0]];
6736 for (d = 1; d < dd->ndim; d++)
6738 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6741 if (cellsize_min < comm->cellsize_limit*1.05)
6743 dd_warning(cr, fplog, "NOTE: the minimum cell size is smaller than 1.05 times the cell size limit, will not turn on dynamic load balancing\n");
6745 /* Change DLB from "auto" to "no". */
6746 comm->dlbState = edlbsOffForever;
6751 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6752 comm->dlbState = edlbsOn;
6756 /* We can set the required cell size info here,
6757 * so we do not need to communicate this.
6758 * The grid is completely uniform.
6760 for (d = 0; d < dd->ndim; d++)
6764 comm->load[d].sum_m = comm->load[d].sum;
6766 nc = dd->nc[dd->dim[d]];
6767 for (i = 0; i < nc; i++)
6769 comm->root[d]->cell_f[i] = i/(real)nc;
6772 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6773 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6776 comm->root[d]->cell_f[nc] = 1.0;
6781 static char *init_bLocalCG(const gmx_mtop_t *mtop)
6786 ncg = ncg_mtop(mtop);
6787 snew(bLocalCG, ncg);
6788 for (cg = 0; cg < ncg; cg++)
6790 bLocalCG[cg] = FALSE;
6796 void dd_init_bondeds(FILE *fplog,
6798 const gmx_mtop_t *mtop,
6799 const gmx_vsite_t *vsite,
6800 const t_inputrec *ir,
6801 gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
6803 gmx_domdec_comm_t *comm;
6805 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
6809 if (comm->bBondComm)
6811 /* Communicate atoms beyond the cut-off for bonded interactions */
6814 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
6816 comm->bLocalCG = init_bLocalCG(mtop);
6820 /* Only communicate atoms based on cut-off */
6821 comm->cglink = NULL;
6822 comm->bLocalCG = NULL;
6826 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
6827 const gmx_mtop_t *mtop, const t_inputrec *ir,
6828 gmx_bool bDynLoadBal, real dlb_scale,
6829 const gmx_ddbox_t *ddbox)
6831 gmx_domdec_comm_t *comm;
6846 fprintf(fplog, "The maximum number of communication pulses is:");
6847 for (d = 0; d < dd->ndim; d++)
6849 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
6851 fprintf(fplog, "\n");
6852 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
6853 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
6854 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
6855 for (d = 0; d < DIM; d++)
6859 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
6866 comm->cellsize_min_dlb[d]/
6867 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6869 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
6872 fprintf(fplog, "\n");
6876 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
6877 fprintf(fplog, "The initial number of communication pulses is:");
6878 for (d = 0; d < dd->ndim; d++)
6880 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
6882 fprintf(fplog, "\n");
6883 fprintf(fplog, "The initial domain decomposition cell size is:");
6884 for (d = 0; d < DIM; d++)
6888 fprintf(fplog, " %c %.2f nm",
6889 dim2char(d), dd->comm->cellsize_min[d]);
6892 fprintf(fplog, "\n\n");
6895 gmx_bool bInterCGVsites = count_intercg_vsites(mtop);
6897 if (comm->bInterCGBondeds ||
6899 dd->bInterCGcons || dd->bInterCGsettles)
6901 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
6902 fprintf(fplog, "%40s %-7s %6.3f nm\n",
6903 "non-bonded interactions", "", comm->cutoff);
6907 limit = dd->comm->cellsize_limit;
6911 if (dynamic_dd_box(ddbox, ir))
6913 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
6915 limit = dd->comm->cellsize_min[XX];
6916 for (d = 1; d < DIM; d++)
6918 limit = std::min(limit, dd->comm->cellsize_min[d]);
6922 if (comm->bInterCGBondeds)
6924 fprintf(fplog, "%40s %-7s %6.3f nm\n",
6925 "two-body bonded interactions", "(-rdd)",
6926 std::max(comm->cutoff, comm->cutoff_mbody));
6927 fprintf(fplog, "%40s %-7s %6.3f nm\n",
6928 "multi-body bonded interactions", "(-rdd)",
6929 (comm->bBondComm || dlbIsOn(dd->comm)) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
6933 fprintf(fplog, "%40s %-7s %6.3f nm\n",
6934 "virtual site constructions", "(-rcon)", limit);
6936 if (dd->bInterCGcons || dd->bInterCGsettles)
6938 sprintf(buf, "atoms separated by up to %d constraints",
6940 fprintf(fplog, "%40s %-7s %6.3f nm\n",
6941 buf, "(-rcon)", limit);
6943 fprintf(fplog, "\n");
6949 static void set_cell_limits_dlb(gmx_domdec_t *dd,
6951 const t_inputrec *ir,
6952 const gmx_ddbox_t *ddbox)
6954 gmx_domdec_comm_t *comm;
6955 int d, dim, npulse, npulse_d_max, npulse_d;
6960 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
6962 /* Determine the maximum number of comm. pulses in one dimension */
6964 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
6966 /* Determine the maximum required number of grid pulses */
6967 if (comm->cellsize_limit >= comm->cutoff)
6969 /* Only a single pulse is required */
6972 else if (!bNoCutOff && comm->cellsize_limit > 0)
6974 /* We round down slightly here to avoid overhead due to the latency
6975 * of extra communication calls when the cut-off
6976 * would be only slightly longer than the cell size.
6977 * Later cellsize_limit is redetermined,
6978 * so we can not miss interactions due to this rounding.
6980 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
6984 /* There is no cell size limit */
6985 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
6988 if (!bNoCutOff && npulse > 1)
6990 /* See if we can do with less pulses, based on dlb_scale */
6992 for (d = 0; d < dd->ndim; d++)
6995 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
6996 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
6997 npulse_d_max = std::max(npulse_d_max, npulse_d);
6999 npulse = std::min(npulse, npulse_d_max);
7002 /* This env var can override npulse */
7003 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7010 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7011 for (d = 0; d < dd->ndim; d++)
7013 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7014 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7015 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7016 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7017 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7019 comm->bVacDLBNoLimit = FALSE;
7023 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7024 if (!comm->bVacDLBNoLimit)
7026 comm->cellsize_limit = std::max(comm->cellsize_limit,
7027 comm->cutoff/comm->maxpulse);
7029 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7030 /* Set the minimum cell size for each DD dimension */
7031 for (d = 0; d < dd->ndim; d++)
7033 if (comm->bVacDLBNoLimit ||
7034 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7036 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7040 comm->cellsize_min_dlb[dd->dim[d]] =
7041 comm->cutoff/comm->cd[d].np_dlb;
7044 if (comm->cutoff_mbody <= 0)
7046 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7054 gmx_bool dd_bonded_molpbc(const gmx_domdec_t *dd, int ePBC)
7056 /* If each molecule is a single charge group
7057 * or we use domain decomposition for each periodic dimension,
7058 * we do not need to take pbc into account for the bonded interactions.
7060 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7063 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7066 /*! \brief Sets grid size limits and PP-PME setup, prints settings to log */
7067 static void set_ddgrid_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7068 const gmx_mtop_t *mtop, const t_inputrec *ir,
7069 const gmx_ddbox_t *ddbox)
7071 gmx_domdec_comm_t *comm;
7077 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7079 init_ddpme(dd, &comm->ddpme[0], 0);
7080 if (comm->npmedecompdim >= 2)
7082 init_ddpme(dd, &comm->ddpme[1], 1);
7087 comm->npmenodes = 0;
7088 if (dd->pme_nodeid >= 0)
7090 gmx_fatal_collective(FARGS, dd->mpi_comm_all, DDMASTER(dd),
7091 "Can not have separate PME ranks without PME electrostatics");
7097 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7099 if (comm->dlbState != edlbsOffForever)
7101 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7104 print_dd_settings(fplog, dd, mtop, ir, dlbIsOn(comm), dlb_scale, ddbox);
7105 if (comm->dlbState == edlbsOffCanTurnOn)
7109 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7111 print_dd_settings(fplog, dd, mtop, ir, TRUE, dlb_scale, ddbox);
7114 if (ir->ePBC == epbcNONE)
7116 vol_frac = 1 - 1/(double)dd->nnodes;
7121 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7125 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7127 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7129 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7132 /*! \brief Set some important DD parameters that can be modified by env.vars */
7133 static void set_dd_envvar_options(FILE *fplog, gmx_domdec_t *dd, int rank_mysim)
7135 gmx_domdec_comm_t *comm = dd->comm;
7137 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
7138 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
7139 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
7140 int recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
7141 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
7142 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
7143 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
7145 if (dd->bSendRecv2 && fplog)
7147 fprintf(fplog, "Will use two sequential MPI_Sendrecv calls instead of two simultaneous non-blocking MPI_Irecv and MPI_Isend pairs for constraint and vsite communication\n");
7154 fprintf(fplog, "Will load balance based on FLOP count\n");
7156 if (comm->eFlop > 1)
7158 srand(1 + rank_mysim);
7160 comm->bRecordLoad = TRUE;
7164 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
7168 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
7169 unsigned long Flags,
7170 ivec nc, int nPmeRanks,
7172 real comm_distance_min, real rconstr,
7173 const char *dlb_opt, real dlb_scale,
7174 const char *sizex, const char *sizey, const char *sizez,
7175 const gmx_mtop_t *mtop,
7176 const t_inputrec *ir,
7177 matrix box, rvec *x,
7179 int *npme_x, int *npme_y)
7186 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
7191 dd->comm = init_dd_comm();
7193 set_dd_envvar_options(fplog, dd, cr->nodeid);
7195 set_dd_limits_and_grid(fplog, cr, dd, Flags,
7197 comm_distance_min, rconstr,
7199 sizex, sizey, sizez,
7205 make_dd_communicators(fplog, cr, dd, dd_rank_order);
7207 if (cr->duty & DUTY_PP)
7209 set_ddgrid_parameters(fplog, dd, dlb_scale, mtop, ir, ddbox);
7211 setup_neighbor_relations(dd);
7214 /* Set overallocation to avoid frequent reallocation of arrays */
7215 set_over_alloc_dd(TRUE);
7217 /* Initialize DD paritioning counters */
7218 dd->comm->partition_step = INT_MIN;
7221 /* We currently don't know the number of threads yet, we set this later */
7224 clear_dd_cycle_counts(dd);
7229 static gmx_bool test_dd_cutoff(t_commrec *cr,
7230 t_state *state, const t_inputrec *ir,
7241 set_ddbox(dd, FALSE, cr, ir, state->box,
7242 TRUE, &dd->comm->cgs_gl, as_rvec_array(state->x.data()), &ddbox);
7246 for (d = 0; d < dd->ndim; d++)
7250 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7251 if (dynamic_dd_box(&ddbox, ir))
7253 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7256 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7258 if (dd->comm->dlbState != edlbsOffForever && dim < ddbox.npbcdim &&
7259 dd->comm->cd[d].np_dlb > 0)
7261 if (np > dd->comm->cd[d].np_dlb)
7266 /* If a current local cell size is smaller than the requested
7267 * cut-off, we could still fix it, but this gets very complicated.
7268 * Without fixing here, we might actually need more checks.
7270 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7277 if (dd->comm->dlbState != edlbsOffForever)
7279 /* If DLB is not active yet, we don't need to check the grid jumps.
7280 * Actually we shouldn't, because then the grid jump data is not set.
7282 if (dlbIsOn(dd->comm) &&
7283 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7288 gmx_sumi(1, &LocallyLimited, cr);
7290 if (LocallyLimited > 0)
7299 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, const t_inputrec *ir,
7302 gmx_bool bCutoffAllowed;
7304 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7308 cr->dd->comm->cutoff = cutoff_req;
7311 return bCutoffAllowed;
7314 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
7316 gmx_domdec_comm_t *comm;
7318 comm = cr->dd->comm;
7320 /* Turn on the DLB limiting (might have been on already) */
7321 comm->bPMELoadBalDLBLimits = TRUE;
7323 /* Change the cut-off limit */
7324 comm->PMELoadBal_max_cutoff = cutoff;
7328 fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
7329 comm->PMELoadBal_max_cutoff);
7333 /* Sets whether we should later check the load imbalance data, so that
7334 * we can trigger dynamic load balancing if enough imbalance has
7337 * Used after PME load balancing unlocks DLB, so that the check
7338 * whether DLB will be useful can happen immediately.
7340 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
7342 if (dd->comm->dlbState == edlbsOffCanTurnOn)
7344 dd->comm->bCheckWhetherToTurnDlbOn = bValue;
7348 /* Returns if we should check whether there has been enough load
7349 * imbalance to trigger dynamic load balancing.
7351 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
7353 const int nddp_chk_dlb = 100;
7355 if (dd->comm->dlbState != edlbsOffCanTurnOn)
7360 /* We should check whether we should use DLB directly after
7362 if (dd->comm->bCheckWhetherToTurnDlbOn)
7364 /* This flag was set when the PME load-balancing routines
7365 unlocked DLB, and should now be cleared. */
7366 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
7369 /* We should also check whether we should use DLB every 100
7370 * partitionings (we do not do this every partioning, so that we
7371 * avoid excessive communication). */
7372 if (dd->comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1)
7380 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
7382 return (dd->comm->dlbState == edlbsOn);
7385 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7387 return (dd->comm->dlbState == edlbsOffTemporarilyLocked);
7390 void dd_dlb_lock(gmx_domdec_t *dd)
7392 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7393 if (dd->comm->dlbState == edlbsOffCanTurnOn)
7395 dd->comm->dlbState = edlbsOffTemporarilyLocked;
7399 void dd_dlb_unlock(gmx_domdec_t *dd)
7401 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7402 if (dd->comm->dlbState == edlbsOffTemporarilyLocked)
7404 dd->comm->dlbState = edlbsOffCanTurnOn;
7405 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
7409 static void merge_cg_buffers(int ncell,
7410 gmx_domdec_comm_dim_t *cd, int pulse,
7412 int *index_gl, int *recv_i,
7413 rvec *cg_cm, rvec *recv_vr,
7415 cginfo_mb_t *cginfo_mb, int *cginfo)
7417 gmx_domdec_ind_t *ind, *ind_p;
7418 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7419 int shift, shift_at;
7421 ind = &cd->ind[pulse];
7423 /* First correct the already stored data */
7424 shift = ind->nrecv[ncell];
7425 for (cell = ncell-1; cell >= 0; cell--)
7427 shift -= ind->nrecv[cell];
7430 /* Move the cg's present from previous grid pulses */
7431 cg0 = ncg_cell[ncell+cell];
7432 cg1 = ncg_cell[ncell+cell+1];
7433 cgindex[cg1+shift] = cgindex[cg1];
7434 for (cg = cg1-1; cg >= cg0; cg--)
7436 index_gl[cg+shift] = index_gl[cg];
7437 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7438 cgindex[cg+shift] = cgindex[cg];
7439 cginfo[cg+shift] = cginfo[cg];
7441 /* Correct the already stored send indices for the shift */
7442 for (p = 1; p <= pulse; p++)
7444 ind_p = &cd->ind[p];
7446 for (c = 0; c < cell; c++)
7448 cg0 += ind_p->nsend[c];
7450 cg1 = cg0 + ind_p->nsend[cell];
7451 for (cg = cg0; cg < cg1; cg++)
7453 ind_p->index[cg] += shift;
7459 /* Merge in the communicated buffers */
7463 for (cell = 0; cell < ncell; cell++)
7465 cg1 = ncg_cell[ncell+cell+1] + shift;
7468 /* Correct the old cg indices */
7469 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7471 cgindex[cg+1] += shift_at;
7474 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7476 /* Copy this charge group from the buffer */
7477 index_gl[cg1] = recv_i[cg0];
7478 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7479 /* Add it to the cgindex */
7480 cg_gl = index_gl[cg1];
7481 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7482 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7483 cgindex[cg1+1] = cgindex[cg1] + nat;
7488 shift += ind->nrecv[cell];
7489 ncg_cell[ncell+cell+1] = cg1;
7493 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7494 int nzone, int cg0, const int *cgindex)
7498 /* Store the atom block boundaries for easy copying of communication buffers
7501 for (zone = 0; zone < nzone; zone++)
7503 for (p = 0; p < cd->np; p++)
7505 cd->ind[p].cell2at0[zone] = cgindex[cg];
7506 cg += cd->ind[p].nrecv[zone];
7507 cd->ind[p].cell2at1[zone] = cgindex[cg];
7512 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7518 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7520 if (!bLocalCG[link->a[i]])
7529 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7531 real c[DIM][4]; /* the corners for the non-bonded communication */
7532 real cr0; /* corner for rounding */
7533 real cr1[4]; /* corners for rounding */
7534 real bc[DIM]; /* corners for bounded communication */
7535 real bcr1; /* corner for rounding for bonded communication */
7538 /* Determine the corners of the domain(s) we are communicating with */
7540 set_dd_corners(const gmx_domdec_t *dd,
7541 int dim0, int dim1, int dim2,
7545 const gmx_domdec_comm_t *comm;
7546 const gmx_domdec_zones_t *zones;
7551 zones = &comm->zones;
7553 /* Keep the compiler happy */
7557 /* The first dimension is equal for all cells */
7558 c->c[0][0] = comm->cell_x0[dim0];
7561 c->bc[0] = c->c[0][0];
7566 /* This cell row is only seen from the first row */
7567 c->c[1][0] = comm->cell_x0[dim1];
7568 /* All rows can see this row */
7569 c->c[1][1] = comm->cell_x0[dim1];
7570 if (dlbIsOn(dd->comm))
7572 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7575 /* For the multi-body distance we need the maximum */
7576 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7579 /* Set the upper-right corner for rounding */
7580 c->cr0 = comm->cell_x1[dim0];
7585 for (j = 0; j < 4; j++)
7587 c->c[2][j] = comm->cell_x0[dim2];
7589 if (dlbIsOn(dd->comm))
7591 /* Use the maximum of the i-cells that see a j-cell */
7592 for (i = 0; i < zones->nizone; i++)
7594 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7599 std::max(c->c[2][j-4],
7600 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7606 /* For the multi-body distance we need the maximum */
7607 c->bc[2] = comm->cell_x0[dim2];
7608 for (i = 0; i < 2; i++)
7610 for (j = 0; j < 2; j++)
7612 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7618 /* Set the upper-right corner for rounding */
7619 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7620 * Only cell (0,0,0) can see cell 7 (1,1,1)
7622 c->cr1[0] = comm->cell_x1[dim1];
7623 c->cr1[3] = comm->cell_x1[dim1];
7624 if (dlbIsOn(dd->comm))
7626 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7629 /* For the multi-body distance we need the maximum */
7630 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7637 /* Determine which cg's we need to send in this pulse from this zone */
7639 get_zone_pulse_cgs(gmx_domdec_t *dd,
7640 int zonei, int zone,
7642 const int *index_gl,
7644 int dim, int dim_ind,
7645 int dim0, int dim1, int dim2,
7646 real r_comm2, real r_bcomm2,
7650 real skew_fac2_d, real skew_fac_01,
7651 rvec *v_d, rvec *v_0, rvec *v_1,
7652 const dd_corners_t *c,
7654 gmx_bool bDistBonded,
7660 gmx_domdec_ind_t *ind,
7661 int **ibuf, int *ibuf_nalloc,
7667 gmx_domdec_comm_t *comm;
7669 gmx_bool bDistMB_pulse;
7671 real r2, rb2, r, tric_sh;
7674 int nsend_z, nsend, nat;
7678 bScrew = (dd->bScrewPBC && dim == XX);
7680 bDistMB_pulse = (bDistMB && bDistBonded);
7686 for (cg = cg0; cg < cg1; cg++)
7690 if (tric_dist[dim_ind] == 0)
7692 /* Rectangular direction, easy */
7693 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7700 r = cg_cm[cg][dim] - c->bc[dim_ind];
7706 /* Rounding gives at most a 16% reduction
7707 * in communicated atoms
7709 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7711 r = cg_cm[cg][dim0] - c->cr0;
7712 /* This is the first dimension, so always r >= 0 */
7719 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7721 r = cg_cm[cg][dim1] - c->cr1[zone];
7728 r = cg_cm[cg][dim1] - c->bcr1;
7738 /* Triclinic direction, more complicated */
7741 /* Rounding, conservative as the skew_fac multiplication
7742 * will slightly underestimate the distance.
7744 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7746 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7747 for (i = dim0+1; i < DIM; i++)
7749 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7751 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7754 rb[dim0] = rn[dim0];
7757 /* Take care that the cell planes along dim0 might not
7758 * be orthogonal to those along dim1 and dim2.
7760 for (i = 1; i <= dim_ind; i++)
7763 if (normal[dim0][dimd] > 0)
7765 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7768 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7773 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7775 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7777 for (i = dim1+1; i < DIM; i++)
7779 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7781 rn[dim1] += tric_sh;
7784 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7785 /* Take care of coupling of the distances
7786 * to the planes along dim0 and dim1 through dim2.
7788 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7789 /* Take care that the cell planes along dim1
7790 * might not be orthogonal to that along dim2.
7792 if (normal[dim1][dim2] > 0)
7794 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7800 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7803 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7804 /* Take care of coupling of the distances
7805 * to the planes along dim0 and dim1 through dim2.
7807 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7808 /* Take care that the cell planes along dim1
7809 * might not be orthogonal to that along dim2.
7811 if (normal[dim1][dim2] > 0)
7813 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7818 /* The distance along the communication direction */
7819 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7821 for (i = dim+1; i < DIM; i++)
7823 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7828 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7829 /* Take care of coupling of the distances
7830 * to the planes along dim0 and dim1 through dim2.
7832 if (dim_ind == 1 && zonei == 1)
7834 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7840 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7843 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7844 /* Take care of coupling of the distances
7845 * to the planes along dim0 and dim1 through dim2.
7847 if (dim_ind == 1 && zonei == 1)
7849 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7857 ((bDistMB && rb2 < r_bcomm2) ||
7858 (bDist2B && r2 < r_bcomm2)) &&
7860 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7861 missing_link(comm->cglink, index_gl[cg],
7864 /* Make an index to the local charge groups */
7865 if (nsend+1 > ind->nalloc)
7867 ind->nalloc = over_alloc_large(nsend+1);
7868 srenew(ind->index, ind->nalloc);
7870 if (nsend+1 > *ibuf_nalloc)
7872 *ibuf_nalloc = over_alloc_large(nsend+1);
7873 srenew(*ibuf, *ibuf_nalloc);
7875 ind->index[nsend] = cg;
7876 (*ibuf)[nsend] = index_gl[cg];
7878 vec_rvec_check_alloc(vbuf, nsend+1);
7880 if (dd->ci[dim] == 0)
7882 /* Correct cg_cm for pbc */
7883 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7886 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7887 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7892 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7895 nat += cgindex[cg+1] - cgindex[cg];
7901 *nsend_z_ptr = nsend_z;
7904 static void setup_dd_communication(gmx_domdec_t *dd,
7905 matrix box, gmx_ddbox_t *ddbox,
7907 t_state *state, PaddedRVecVector *f)
7909 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7910 int nzone, nzone_send, zone, zonei, cg0, cg1;
7911 int c, i, cg, cg_gl, nrcg;
7912 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7913 gmx_domdec_comm_t *comm;
7914 gmx_domdec_zones_t *zones;
7915 gmx_domdec_comm_dim_t *cd;
7916 gmx_domdec_ind_t *ind;
7917 cginfo_mb_t *cginfo_mb;
7918 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7919 real r_comm2, r_bcomm2;
7920 dd_corners_t corners;
7922 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7923 real skew_fac2_d, skew_fac_01;
7930 fprintf(debug, "Setting up DD communication\n");
7937 /* Initialize the thread data.
7938 * This can not be done in init_domain_decomposition,
7939 * as the numbers of threads is determined later.
7941 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7944 snew(comm->dth, comm->nth);
7948 switch (fr->cutoff_scheme)
7954 cg_cm = as_rvec_array(state->x.data());
7957 gmx_incons("unimplemented");
7961 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7963 /* Check if we need to use triclinic distances */
7964 tric_dist[dim_ind] = 0;
7965 for (i = 0; i <= dim_ind; i++)
7967 if (ddbox->tric_dir[dd->dim[i]])
7969 tric_dist[dim_ind] = 1;
7974 bBondComm = comm->bBondComm;
7976 /* Do we need to determine extra distances for multi-body bondeds? */
7977 bDistMB = (comm->bInterCGMultiBody && dlbIsOn(dd->comm) && dd->ndim > 1);
7979 /* Do we need to determine extra distances for only two-body bondeds? */
7980 bDist2B = (bBondComm && !bDistMB);
7982 r_comm2 = gmx::square(comm->cutoff);
7983 r_bcomm2 = gmx::square(comm->cutoff_mbody);
7987 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, std::sqrt(r_bcomm2));
7990 zones = &comm->zones;
7993 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
7994 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
7996 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
7998 /* Triclinic stuff */
7999 normal = ddbox->normal;
8003 v_0 = ddbox->v[dim0];
8004 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8006 /* Determine the coupling coefficient for the distances
8007 * to the cell planes along dim0 and dim1 through dim2.
8008 * This is required for correct rounding.
8011 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8014 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8020 v_1 = ddbox->v[dim1];
8023 zone_cg_range = zones->cg_range;
8024 index_gl = dd->index_gl;
8025 cgindex = dd->cgindex;
8026 cginfo_mb = fr->cginfo_mb;
8028 zone_cg_range[0] = 0;
8029 zone_cg_range[1] = dd->ncg_home;
8030 comm->zone_ncg1[0] = dd->ncg_home;
8031 pos_cg = dd->ncg_home;
8033 nat_tot = dd->nat_home;
8035 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8037 dim = dd->dim[dim_ind];
8038 cd = &comm->cd[dim_ind];
8040 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8042 /* No pbc in this dimension, the first node should not comm. */
8050 v_d = ddbox->v[dim];
8051 skew_fac2_d = gmx::square(ddbox->skew_fac[dim]);
8053 cd->bInPlace = TRUE;
8054 for (p = 0; p < cd->np; p++)
8056 /* Only atoms communicated in the first pulse are used
8057 * for multi-body bonded interactions or for bBondComm.
8059 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8064 for (zone = 0; zone < nzone_send; zone++)
8066 if (tric_dist[dim_ind] && dim_ind > 0)
8068 /* Determine slightly more optimized skew_fac's
8070 * This reduces the number of communicated atoms
8071 * by about 10% for 3D DD of rhombic dodecahedra.
8073 for (dimd = 0; dimd < dim; dimd++)
8075 sf2_round[dimd] = 1;
8076 if (ddbox->tric_dir[dimd])
8078 for (i = dd->dim[dimd]+1; i < DIM; i++)
8080 /* If we are shifted in dimension i
8081 * and the cell plane is tilted forward
8082 * in dimension i, skip this coupling.
8084 if (!(zones->shift[nzone+zone][i] &&
8085 ddbox->v[dimd][i][dimd] >= 0))
8088 gmx::square(ddbox->v[dimd][i][dimd]);
8091 sf2_round[dimd] = 1/sf2_round[dimd];
8096 zonei = zone_perm[dim_ind][zone];
8099 /* Here we permutate the zones to obtain a convenient order
8100 * for neighbor searching
8102 cg0 = zone_cg_range[zonei];
8103 cg1 = zone_cg_range[zonei+1];
8107 /* Look only at the cg's received in the previous grid pulse
8109 cg1 = zone_cg_range[nzone+zone+1];
8110 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8113 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8114 for (th = 0; th < comm->nth; th++)
8118 gmx_domdec_ind_t *ind_p;
8119 int **ibuf_p, *ibuf_nalloc_p;
8121 int *nsend_p, *nat_p;
8127 /* Thread 0 writes in the comm buffers */
8129 ibuf_p = &comm->buf_int;
8130 ibuf_nalloc_p = &comm->nalloc_int;
8131 vbuf_p = &comm->vbuf;
8134 nsend_zone_p = &ind->nsend[zone];
8138 /* Other threads write into temp buffers */
8139 ind_p = &comm->dth[th].ind;
8140 ibuf_p = &comm->dth[th].ibuf;
8141 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8142 vbuf_p = &comm->dth[th].vbuf;
8143 nsend_p = &comm->dth[th].nsend;
8144 nat_p = &comm->dth[th].nat;
8145 nsend_zone_p = &comm->dth[th].nsend_zone;
8147 comm->dth[th].nsend = 0;
8148 comm->dth[th].nat = 0;
8149 comm->dth[th].nsend_zone = 0;
8159 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8160 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8163 /* Get the cg's for this pulse in this zone */
8164 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8166 dim, dim_ind, dim0, dim1, dim2,
8169 normal, skew_fac2_d, skew_fac_01,
8170 v_d, v_0, v_1, &corners, sf2_round,
8171 bDistBonded, bBondComm,
8175 ibuf_p, ibuf_nalloc_p,
8180 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
8183 /* Append data of threads>=1 to the communication buffers */
8184 for (th = 1; th < comm->nth; th++)
8186 dd_comm_setup_work_t *dth;
8189 dth = &comm->dth[th];
8191 ns1 = nsend + dth->nsend_zone;
8192 if (ns1 > ind->nalloc)
8194 ind->nalloc = over_alloc_dd(ns1);
8195 srenew(ind->index, ind->nalloc);
8197 if (ns1 > comm->nalloc_int)
8199 comm->nalloc_int = over_alloc_dd(ns1);
8200 srenew(comm->buf_int, comm->nalloc_int);
8202 if (ns1 > comm->vbuf.nalloc)
8204 comm->vbuf.nalloc = over_alloc_dd(ns1);
8205 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8208 for (i = 0; i < dth->nsend_zone; i++)
8210 ind->index[nsend] = dth->ind.index[i];
8211 comm->buf_int[nsend] = dth->ibuf[i];
8212 copy_rvec(dth->vbuf.v[i],
8213 comm->vbuf.v[nsend]);
8217 ind->nsend[zone] += dth->nsend_zone;
8220 /* Clear the counts in case we do not have pbc */
8221 for (zone = nzone_send; zone < nzone; zone++)
8223 ind->nsend[zone] = 0;
8225 ind->nsend[nzone] = nsend;
8226 ind->nsend[nzone+1] = nat;
8227 /* Communicate the number of cg's and atoms to receive */
8228 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8229 ind->nsend, nzone+2,
8230 ind->nrecv, nzone+2);
8232 /* The rvec buffer is also required for atom buffers of size nsend
8233 * in dd_move_x and dd_move_f.
8235 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8239 /* We can receive in place if only the last zone is not empty */
8240 for (zone = 0; zone < nzone-1; zone++)
8242 if (ind->nrecv[zone] > 0)
8244 cd->bInPlace = FALSE;
8249 /* The int buffer is only required here for the cg indices */
8250 if (ind->nrecv[nzone] > comm->nalloc_int2)
8252 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8253 srenew(comm->buf_int2, comm->nalloc_int2);
8255 /* The rvec buffer is also required for atom buffers
8256 * of size nrecv in dd_move_x and dd_move_f.
8258 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8259 vec_rvec_check_alloc(&comm->vbuf2, i);
8263 /* Make space for the global cg indices */
8264 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8265 || dd->cg_nalloc == 0)
8267 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8268 srenew(index_gl, dd->cg_nalloc);
8269 srenew(cgindex, dd->cg_nalloc+1);
8271 /* Communicate the global cg indices */
8274 recv_i = index_gl + pos_cg;
8278 recv_i = comm->buf_int2;
8280 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8281 comm->buf_int, nsend,
8282 recv_i, ind->nrecv[nzone]);
8284 /* Make space for cg_cm */
8285 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8286 if (fr->cutoff_scheme == ecutsGROUP)
8292 cg_cm = as_rvec_array(state->x.data());
8294 /* Communicate cg_cm */
8297 recv_vr = cg_cm + pos_cg;
8301 recv_vr = comm->vbuf2.v;
8303 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8304 comm->vbuf.v, nsend,
8305 recv_vr, ind->nrecv[nzone]);
8307 /* Make the charge group index */
8310 zone = (p == 0 ? 0 : nzone - 1);
8311 while (zone < nzone)
8313 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8315 cg_gl = index_gl[pos_cg];
8316 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8317 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8318 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8321 /* Update the charge group presence,
8322 * so we can use it in the next pass of the loop.
8324 comm->bLocalCG[cg_gl] = TRUE;
8330 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8333 zone_cg_range[nzone+zone] = pos_cg;
8338 /* This part of the code is never executed with bBondComm. */
8339 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8340 index_gl, recv_i, cg_cm, recv_vr,
8341 cgindex, fr->cginfo_mb, fr->cginfo);
8342 pos_cg += ind->nrecv[nzone];
8344 nat_tot += ind->nrecv[nzone+1];
8348 /* Store the atom block for easy copying of communication buffers */
8349 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8353 dd->index_gl = index_gl;
8354 dd->cgindex = cgindex;
8356 dd->ncg_tot = zone_cg_range[zones->n];
8357 dd->nat_tot = nat_tot;
8358 comm->nat[ddnatHOME] = dd->nat_home;
8359 for (i = ddnatZONE; i < ddnatNR; i++)
8361 comm->nat[i] = dd->nat_tot;
8366 /* We don't need to update cginfo, since that was alrady done above.
8367 * So we pass NULL for the forcerec.
8369 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8370 NULL, comm->bLocalCG);
8375 fprintf(debug, "Finished setting up DD communication, zones:");
8376 for (c = 0; c < zones->n; c++)
8378 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8380 fprintf(debug, "\n");
8384 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8388 for (c = 0; c < zones->nizone; c++)
8390 zones->izone[c].cg1 = zones->cg_range[c+1];
8391 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8392 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8396 static void set_zones_size(gmx_domdec_t *dd,
8397 matrix box, const gmx_ddbox_t *ddbox,
8398 int zone_start, int zone_end)
8400 gmx_domdec_comm_t *comm;
8401 gmx_domdec_zones_t *zones;
8410 zones = &comm->zones;
8412 /* Do we need to determine extra distances for multi-body bondeds? */
8413 bDistMB = (comm->bInterCGMultiBody && dlbIsOn(dd->comm) && dd->ndim > 1);
8415 for (z = zone_start; z < zone_end; z++)
8417 /* Copy cell limits to zone limits.
8418 * Valid for non-DD dims and non-shifted dims.
8420 copy_rvec(comm->cell_x0, zones->size[z].x0);
8421 copy_rvec(comm->cell_x1, zones->size[z].x1);
8424 for (d = 0; d < dd->ndim; d++)
8428 for (z = 0; z < zones->n; z++)
8430 /* With a staggered grid we have different sizes
8431 * for non-shifted dimensions.
8433 if (dlbIsOn(dd->comm) && zones->shift[z][dim] == 0)
8437 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8438 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8442 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8443 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8449 rcmbs = comm->cutoff_mbody;
8450 if (ddbox->tric_dir[dim])
8452 rcs /= ddbox->skew_fac[dim];
8453 rcmbs /= ddbox->skew_fac[dim];
8456 /* Set the lower limit for the shifted zone dimensions */
8457 for (z = zone_start; z < zone_end; z++)
8459 if (zones->shift[z][dim] > 0)
8462 if (!dlbIsOn(dd->comm) || d == 0)
8464 zones->size[z].x0[dim] = comm->cell_x1[dim];
8465 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8469 /* Here we take the lower limit of the zone from
8470 * the lowest domain of the zone below.
8474 zones->size[z].x0[dim] =
8475 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8481 zones->size[z].x0[dim] =
8482 zones->size[zone_perm[2][z-4]].x0[dim];
8486 zones->size[z].x0[dim] =
8487 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8490 /* A temporary limit, is updated below */
8491 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8495 for (zi = 0; zi < zones->nizone; zi++)
8497 if (zones->shift[zi][dim] == 0)
8499 /* This takes the whole zone into account.
8500 * With multiple pulses this will lead
8501 * to a larger zone then strictly necessary.
8503 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8504 zones->size[zi].x1[dim]+rcmbs);
8512 /* Loop over the i-zones to set the upper limit of each
8515 for (zi = 0; zi < zones->nizone; zi++)
8517 if (zones->shift[zi][dim] == 0)
8519 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8521 if (zones->shift[z][dim] > 0)
8523 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8524 zones->size[zi].x1[dim]+rcs);
8531 for (z = zone_start; z < zone_end; z++)
8533 /* Initialization only required to keep the compiler happy */
8534 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8537 /* To determine the bounding box for a zone we need to find
8538 * the extreme corners of 4, 2 or 1 corners.
8540 nc = 1 << (ddbox->nboundeddim - 1);
8542 for (c = 0; c < nc; c++)
8544 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8548 corner[YY] = zones->size[z].x0[YY];
8552 corner[YY] = zones->size[z].x1[YY];
8556 corner[ZZ] = zones->size[z].x0[ZZ];
8560 corner[ZZ] = zones->size[z].x1[ZZ];
8562 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8563 box[ZZ][1 - dd->dim[0]] != 0)
8565 /* With 1D domain decomposition the cg's are not in
8566 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8567 * Shift the corner of the z-vector back to along the box
8568 * vector of dimension d, so it will later end up at 0 along d.
8569 * This can affect the location of this corner along dd->dim[0]
8570 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8572 int d = 1 - dd->dim[0];
8574 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8576 /* Apply the triclinic couplings */
8577 assert(ddbox->npbcdim <= DIM);
8578 for (i = YY; i < ddbox->npbcdim; i++)
8580 for (j = XX; j < i; j++)
8582 corner[j] += corner[i]*box[i][j]/box[i][i];
8587 copy_rvec(corner, corner_min);
8588 copy_rvec(corner, corner_max);
8592 for (i = 0; i < DIM; i++)
8594 corner_min[i] = std::min(corner_min[i], corner[i]);
8595 corner_max[i] = std::max(corner_max[i], corner[i]);
8599 /* Copy the extreme cornes without offset along x */
8600 for (i = 0; i < DIM; i++)
8602 zones->size[z].bb_x0[i] = corner_min[i];
8603 zones->size[z].bb_x1[i] = corner_max[i];
8605 /* Add the offset along x */
8606 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8607 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8610 if (zone_start == 0)
8613 for (dim = 0; dim < DIM; dim++)
8615 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8617 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8622 for (z = zone_start; z < zone_end; z++)
8624 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8626 zones->size[z].x0[XX], zones->size[z].x1[XX],
8627 zones->size[z].x0[YY], zones->size[z].x1[YY],
8628 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8629 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8631 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8632 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8633 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8638 static int comp_cgsort(const void *a, const void *b)
8642 gmx_cgsort_t *cga, *cgb;
8643 cga = (gmx_cgsort_t *)a;
8644 cgb = (gmx_cgsort_t *)b;
8646 comp = cga->nsc - cgb->nsc;
8649 comp = cga->ind_gl - cgb->ind_gl;
8655 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8660 /* Order the data */
8661 for (i = 0; i < n; i++)
8663 buf[i] = a[sort[i].ind];
8666 /* Copy back to the original array */
8667 for (i = 0; i < n; i++)
8673 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8678 /* Order the data */
8679 for (i = 0; i < n; i++)
8681 copy_rvec(v[sort[i].ind], buf[i]);
8684 /* Copy back to the original array */
8685 for (i = 0; i < n; i++)
8687 copy_rvec(buf[i], v[i]);
8691 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8694 int a, atot, cg, cg0, cg1, i;
8696 if (cgindex == NULL)
8698 /* Avoid the useless loop of the atoms within a cg */
8699 order_vec_cg(ncg, sort, v, buf);
8704 /* Order the data */
8706 for (cg = 0; cg < ncg; cg++)
8708 cg0 = cgindex[sort[cg].ind];
8709 cg1 = cgindex[sort[cg].ind+1];
8710 for (i = cg0; i < cg1; i++)
8712 copy_rvec(v[i], buf[a]);
8718 /* Copy back to the original array */
8719 for (a = 0; a < atot; a++)
8721 copy_rvec(buf[a], v[a]);
8725 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8726 int nsort_new, gmx_cgsort_t *sort_new,
8727 gmx_cgsort_t *sort1)
8731 /* The new indices are not very ordered, so we qsort them */
8732 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8734 /* sort2 is already ordered, so now we can merge the two arrays */
8738 while (i2 < nsort2 || i_new < nsort_new)
8742 sort1[i1++] = sort_new[i_new++];
8744 else if (i_new == nsort_new)
8746 sort1[i1++] = sort2[i2++];
8748 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8749 (sort2[i2].nsc == sort_new[i_new].nsc &&
8750 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8752 sort1[i1++] = sort2[i2++];
8756 sort1[i1++] = sort_new[i_new++];
8761 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8763 gmx_domdec_sort_t *sort;
8764 gmx_cgsort_t *cgsort, *sort_i;
8765 int ncg_new, nsort2, nsort_new, i, *a, moved;
8767 sort = dd->comm->sort;
8769 a = fr->ns->grid->cell_index;
8771 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns->grid->ncells;
8773 if (ncg_home_old >= 0)
8775 /* The charge groups that remained in the same ns grid cell
8776 * are completely ordered. So we can sort efficiently by sorting
8777 * the charge groups that did move into the stationary list.
8782 for (i = 0; i < dd->ncg_home; i++)
8784 /* Check if this cg did not move to another node */
8787 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8789 /* This cg is new on this node or moved ns grid cell */
8790 if (nsort_new >= sort->sort_new_nalloc)
8792 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8793 srenew(sort->sort_new, sort->sort_new_nalloc);
8795 sort_i = &(sort->sort_new[nsort_new++]);
8799 /* This cg did not move */
8800 sort_i = &(sort->sort2[nsort2++]);
8802 /* Sort on the ns grid cell indices
8803 * and the global topology index.
8804 * index_gl is irrelevant with cell ns,
8805 * but we set it here anyhow to avoid a conditional.
8808 sort_i->ind_gl = dd->index_gl[i];
8815 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8818 /* Sort efficiently */
8819 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8824 cgsort = sort->sort;
8826 for (i = 0; i < dd->ncg_home; i++)
8828 /* Sort on the ns grid cell indices
8829 * and the global topology index
8831 cgsort[i].nsc = a[i];
8832 cgsort[i].ind_gl = dd->index_gl[i];
8834 if (cgsort[i].nsc < moved)
8841 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8843 /* Determine the order of the charge groups using qsort */
8844 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8850 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8856 sort = dd->comm->sort->sort;
8858 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8861 for (i = 0; i < na; i++)
8865 sort[ncg_new].ind = a[i];
8873 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
8876 gmx_domdec_sort_t *sort;
8877 gmx_cgsort_t *cgsort;
8879 int ncg_new, i, *ibuf, cgsize;
8882 sort = dd->comm->sort;
8884 if (dd->ncg_home > sort->sort_nalloc)
8886 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8887 srenew(sort->sort, sort->sort_nalloc);
8888 srenew(sort->sort2, sort->sort_nalloc);
8890 cgsort = sort->sort;
8892 switch (fr->cutoff_scheme)
8895 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8898 ncg_new = dd_sort_order_nbnxn(dd, fr);
8901 gmx_incons("unimplemented");
8905 /* We alloc with the old size, since cgindex is still old */
8906 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8907 vbuf = dd->comm->vbuf.v;
8911 cgindex = dd->cgindex;
8918 /* Remove the charge groups which are no longer at home here */
8919 dd->ncg_home = ncg_new;
8922 fprintf(debug, "Set the new home charge group count to %d\n",
8926 /* Reorder the state */
8927 for (i = 0; i < estNR; i++)
8929 if (EST_DISTR(i) && (state->flags & (1<<i)))
8934 order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->x.data()), vbuf);
8937 order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->v.data()), vbuf);
8939 case est_SDX_NOTSUPPORTED:
8942 order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->cg_p.data()), vbuf);
8946 case estDISRE_INITF:
8947 case estDISRE_RM3TAV:
8948 case estORIRE_INITF:
8950 /* No ordering required */
8953 gmx_incons("Unknown state entry encountered in dd_sort_state");
8958 if (fr->cutoff_scheme == ecutsGROUP)
8961 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8964 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8966 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8967 srenew(sort->ibuf, sort->ibuf_nalloc);
8970 /* Reorder the global cg index */
8971 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
8972 /* Reorder the cginfo */
8973 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
8974 /* Rebuild the local cg index */
8978 for (i = 0; i < dd->ncg_home; i++)
8980 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8981 ibuf[i+1] = ibuf[i] + cgsize;
8983 for (i = 0; i < dd->ncg_home+1; i++)
8985 dd->cgindex[i] = ibuf[i];
8990 for (i = 0; i < dd->ncg_home+1; i++)
8995 /* Set the home atom number */
8996 dd->nat_home = dd->cgindex[dd->ncg_home];
8998 if (fr->cutoff_scheme == ecutsVERLET)
9000 /* The atoms are now exactly in grid order, update the grid order */
9001 nbnxn_set_atomorder(fr->nbv->nbs);
9005 /* Copy the sorted ns cell indices back to the ns grid struct */
9006 for (i = 0; i < dd->ncg_home; i++)
9008 fr->ns->grid->cell_index[i] = cgsort[i].nsc;
9010 fr->ns->grid->nr = dd->ncg_home;
9014 static void add_dd_statistics(gmx_domdec_t *dd)
9016 gmx_domdec_comm_t *comm;
9021 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9023 comm->sum_nat[ddnat-ddnatZONE] +=
9024 comm->nat[ddnat] - comm->nat[ddnat-1];
9029 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9031 gmx_domdec_comm_t *comm;
9036 /* Reset all the statistics and counters for total run counting */
9037 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9039 comm->sum_nat[ddnat-ddnatZONE] = 0;
9043 comm->load_step = 0;
9046 clear_ivec(comm->load_lim);
9051 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9053 gmx_domdec_comm_t *comm;
9057 comm = cr->dd->comm;
9059 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9066 fprintf(fplog, "\n D O M A I N D E C O M P O S I T I O N S T A T I S T I C S\n\n");
9068 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9070 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9075 " av. #atoms communicated per step for force: %d x %.1f\n",
9079 if (cr->dd->vsite_comm)
9082 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9083 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9088 if (cr->dd->constraint_comm)
9091 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9092 1 + ir->nLincsIter, av);
9096 gmx_incons(" Unknown type for DD statistics");
9099 fprintf(fplog, "\n");
9101 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9103 print_dd_load_av(fplog, cr->dd);
9107 void dd_partition_system(FILE *fplog,
9110 gmx_bool bMasterState,
9112 t_state *state_global,
9113 const gmx_mtop_t *top_global,
9114 const t_inputrec *ir,
9115 t_state *state_local,
9116 PaddedRVecVector *f,
9118 gmx_localtop_t *top_local,
9121 gmx_constr_t constr,
9123 gmx_wallcycle_t wcycle,
9127 gmx_domdec_comm_t *comm;
9128 gmx_ddbox_t ddbox = {0};
9130 gmx_int64_t step_pcoupl;
9131 rvec cell_ns_x0, cell_ns_x1;
9132 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9133 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bTurnOnDLB, bLogLoad;
9134 gmx_bool bRedist, bSortCG, bResortAll;
9135 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9139 wallcycle_start(wcycle, ewcDOMDEC);
9144 bBoxChanged = (bMasterState || inputrecDeform(ir));
9145 if (ir->epc != epcNO)
9147 /* With nstpcouple > 1 pressure coupling happens.
9148 * one step after calculating the pressure.
9149 * Box scaling happens at the end of the MD step,
9150 * after the DD partitioning.
9151 * We therefore have to do DLB in the first partitioning
9152 * after an MD step where P-coupling occurred.
9153 * We need to determine the last step in which p-coupling occurred.
9154 * MRS -- need to validate this for vv?
9159 step_pcoupl = step - 1;
9163 step_pcoupl = ((step - 1)/n)*n + 1;
9165 if (step_pcoupl >= comm->partition_step)
9171 bNStGlobalComm = (step % nstglobalcomm == 0);
9179 /* Should we do dynamic load balacing this step?
9180 * Since it requires (possibly expensive) global communication,
9181 * we might want to do DLB less frequently.
9183 if (bBoxChanged || ir->epc != epcNO)
9185 bDoDLB = bBoxChanged;
9189 bDoDLB = bNStGlobalComm;
9193 /* Check if we have recorded loads on the nodes */
9194 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9196 bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
9198 /* Print load every nstlog, first and last step to the log file */
9199 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9200 comm->n_load_collect == 0 ||
9202 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9204 /* Avoid extra communication due to verbose screen output
9205 * when nstglobalcomm is set.
9207 if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
9208 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9210 get_load_distribution(dd, wcycle);
9215 dd_print_load(fplog, dd, step-1);
9219 dd_print_load_verbose(dd);
9222 comm->n_load_collect++;
9224 if (bCheckWhetherToTurnDlbOn)
9226 /* Since the timings are node dependent, the master decides */
9229 /* Here we check if the max PME rank load is more than 0.98
9230 * the max PP force load. If so, PP DLB will not help,
9231 * since we are (almost) limited by PME. Furthermore,
9232 * DLB will cause a significant extra x/f redistribution
9233 * cost on the PME ranks, which will then surely result
9234 * in lower total performance.
9235 * This check might be fragile, since one measurement
9236 * below 0.98 (although only done once every 100 DD part.)
9237 * could turn on DLB for the rest of the run.
9239 if (cr->npmenodes > 0 &&
9240 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9247 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9251 fprintf(debug, "step %s, imb loss %f\n",
9252 gmx_step_str(step, sbuf),
9253 dd_force_imb_perf_loss(dd));
9256 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9259 turn_on_dlb(fplog, cr, step);
9264 comm->n_load_have++;
9267 cgs_gl = &comm->cgs_gl;
9272 /* Clear the old state */
9273 clear_dd_indices(dd, 0, 0);
9276 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9277 TRUE, cgs_gl, as_rvec_array(state_global->x.data()), &ddbox);
9279 get_cg_distribution(fplog, dd, cgs_gl,
9280 state_global->box, &ddbox, as_rvec_array(state_global->x.data()));
9282 dd_distribute_state(dd, cgs_gl,
9283 state_global, state_local, f);
9285 dd_make_local_cgs(dd, &top_local->cgs);
9287 /* Ensure that we have space for the new distribution */
9288 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9290 if (fr->cutoff_scheme == ecutsGROUP)
9292 calc_cgcm(fplog, 0, dd->ncg_home,
9293 &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
9296 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9298 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9300 else if (state_local->ddp_count != dd->ddp_count)
9302 if (state_local->ddp_count > dd->ddp_count)
9304 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9307 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9309 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count_cg_gl (%d) != state_local->ddp_count (%d)", state_local->ddp_count_cg_gl, state_local->ddp_count);
9312 /* Clear the old state */
9313 clear_dd_indices(dd, 0, 0);
9315 /* Build the new indices */
9316 rebuild_cgindex(dd, cgs_gl->index, state_local);
9317 make_dd_indices(dd, cgs_gl->index, 0);
9318 ncgindex_set = dd->ncg_home;
9320 if (fr->cutoff_scheme == ecutsGROUP)
9322 /* Redetermine the cg COMs */
9323 calc_cgcm(fplog, 0, dd->ncg_home,
9324 &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
9327 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9329 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9331 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9332 TRUE, &top_local->cgs, as_rvec_array(state_local->x.data()), &ddbox);
9334 bRedist = dlbIsOn(comm);
9338 /* We have the full state, only redistribute the cgs */
9340 /* Clear the non-home indices */
9341 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9344 /* Avoid global communication for dim's without pbc and -gcom */
9345 if (!bNStGlobalComm)
9347 copy_rvec(comm->box0, ddbox.box0 );
9348 copy_rvec(comm->box_size, ddbox.box_size);
9350 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9351 bNStGlobalComm, &top_local->cgs, as_rvec_array(state_local->x.data()), &ddbox);
9356 /* For dim's without pbc and -gcom */
9357 copy_rvec(ddbox.box0, comm->box0 );
9358 copy_rvec(ddbox.box_size, comm->box_size);
9360 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9363 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9365 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9368 /* Check if we should sort the charge groups */
9369 bSortCG = (bMasterState || bRedist);
9371 ncg_home_old = dd->ncg_home;
9376 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9378 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9380 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9382 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9385 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9387 &comm->cell_x0, &comm->cell_x1,
9388 dd->ncg_home, fr->cg_cm,
9389 cell_ns_x0, cell_ns_x1, &grid_density);
9393 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9396 switch (fr->cutoff_scheme)
9399 copy_ivec(fr->ns->grid->n, ncells_old);
9400 grid_first(fplog, fr->ns->grid, dd, &ddbox,
9401 state_local->box, cell_ns_x0, cell_ns_x1,
9402 fr->rlist, grid_density);
9405 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9408 gmx_incons("unimplemented");
9410 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9411 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9415 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9417 /* Sort the state on charge group position.
9418 * This enables exact restarts from this step.
9419 * It also improves performance by about 15% with larger numbers
9420 * of atoms per node.
9423 /* Fill the ns grid with the home cell,
9424 * so we can sort with the indices.
9426 set_zones_ncg_home(dd);
9428 switch (fr->cutoff_scheme)
9431 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9433 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9435 comm->zones.size[0].bb_x0,
9436 comm->zones.size[0].bb_x1,
9438 comm->zones.dens_zone0,
9440 as_rvec_array(state_local->x.data()),
9441 ncg_moved, bRedist ? comm->moved : NULL,
9442 fr->nbv->grp[eintLocal].kernel_type,
9443 fr->nbv->grp[eintLocal].nbat);
9445 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9448 fill_grid(&comm->zones, fr->ns->grid, dd->ncg_home,
9449 0, dd->ncg_home, fr->cg_cm);
9451 copy_ivec(fr->ns->grid->n, ncells_new);
9454 gmx_incons("unimplemented");
9457 bResortAll = bMasterState;
9459 /* Check if we can user the old order and ns grid cell indices
9460 * of the charge groups to sort the charge groups efficiently.
9462 if (ncells_new[XX] != ncells_old[XX] ||
9463 ncells_new[YY] != ncells_old[YY] ||
9464 ncells_new[ZZ] != ncells_old[ZZ])
9471 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9472 gmx_step_str(step, sbuf), dd->ncg_home);
9474 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9475 bResortAll ? -1 : ncg_home_old);
9477 /* After sorting and compacting we set the correct size */
9478 dd_resize_state(state_local, f, dd->nat_home);
9480 /* Rebuild all the indices */
9481 ga2la_clear(dd->ga2la);
9484 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9487 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9489 /* Setup up the communication and communicate the coordinates */
9490 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9492 /* Set the indices */
9493 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9495 /* Set the charge group boundaries for neighbor searching */
9496 set_cg_boundaries(&comm->zones);
9498 if (fr->cutoff_scheme == ecutsVERLET)
9500 set_zones_size(dd, state_local->box, &ddbox,
9501 bSortCG ? 1 : 0, comm->zones.n);
9504 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9507 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9508 -1,as_rvec_array(state_local->x.data()),state_local->box);
9511 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9513 /* Extract a local topology from the global topology */
9514 for (i = 0; i < dd->ndim; i++)
9516 np[dd->dim[i]] = comm->cd[i].np;
9518 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9519 comm->cellsize_min, np,
9521 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : as_rvec_array(state_local->x.data()),
9522 vsite, top_global, top_local);
9524 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9526 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9528 /* Set up the special atom communication */
9529 n = comm->nat[ddnatZONE];
9530 for (i = ddnatZONE+1; i < ddnatNR; i++)
9535 if (vsite && vsite->n_intercg_vsite)
9537 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9541 if (dd->bInterCGcons || dd->bInterCGsettles)
9543 /* Only for inter-cg constraints we need special code */
9544 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9545 constr, ir->nProjOrder,
9546 top_local->idef.il);
9550 gmx_incons("Unknown special atom type setup");
9555 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9557 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9559 /* Make space for the extra coordinates for virtual site
9560 * or constraint communication.
9562 state_local->natoms = comm->nat[ddnatNR-1];
9564 dd_resize_state(state_local, f, state_local->natoms);
9566 if (fr->bF_NoVirSum)
9568 if (vsite && vsite->n_intercg_vsite)
9570 nat_f_novirsum = comm->nat[ddnatVSITE];
9574 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9576 nat_f_novirsum = dd->nat_tot;
9580 nat_f_novirsum = dd->nat_home;
9589 /* Set the number of atoms required for the force calculation.
9590 * Forces need to be constrained when doing energy
9591 * minimization. For simple simulations we could avoid some
9592 * allocation, zeroing and copying, but this is probably not worth
9593 * the complications and checking.
9595 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9596 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9598 /* Update atom data for mdatoms and several algorithms */
9599 mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr,
9600 NULL, mdatoms, vsite, NULL);
9602 if (ir->implicit_solvent)
9604 make_local_gb(cr, fr->born, ir->gb_algorithm);
9607 if (!(cr->duty & DUTY_PME))
9609 /* Send the charges and/or c6/sigmas to our PME only node */
9610 gmx_pme_send_parameters(cr,
9612 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9613 mdatoms->chargeA, mdatoms->chargeB,
9614 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9615 mdatoms->sigmaA, mdatoms->sigmaB,
9616 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9621 set_constraints(constr, top_local, ir, mdatoms, cr);
9626 /* Update the local pull groups */
9627 dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
9632 /* Update the local rotation groups */
9633 dd_make_local_rotation_groups(dd, ir->rot);
9636 if (ir->eSwapCoords != eswapNO)
9638 /* Update the local groups needed for ion swapping */
9639 dd_make_local_swap_groups(dd, ir->swap);
9642 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9643 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9645 add_dd_statistics(dd);
9647 /* Make sure we only count the cycles for this DD partitioning */
9648 clear_dd_cycle_counts(dd);
9650 /* Because the order of the atoms might have changed since
9651 * the last vsite construction, we need to communicate the constructing
9652 * atom coordinates again (for spreading the forces this MD step).
9654 dd_move_x_vsites(dd, state_local->box, as_rvec_array(state_local->x.data()));
9656 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9658 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9660 dd_move_x(dd, state_local->box, as_rvec_array(state_local->x.data()));
9661 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9662 -1, as_rvec_array(state_local->x.data()), state_local->box);
9665 /* Store the partitioning step */
9666 comm->partition_step = step;
9668 /* Increase the DD partitioning counter */
9670 /* The state currently matches this DD partitioning count, store it */
9671 state_local->ddp_count = dd->ddp_count;
9674 /* The DD master node knows the complete cg distribution,
9675 * store the count so we can possibly skip the cg info communication.
9677 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9680 if (comm->DD_debug > 0)
9682 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9683 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9684 "after partitioning");
9687 wallcycle_stop(wcycle, ewcDOMDEC);