2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014, 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.
38 #include "gromacs/legacyheaders/domdec.h"
51 #include "gromacs/fileio/gmxfio.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/imd/imd.h"
54 #include "gromacs/legacyheaders/bonded-threading.h"
55 #include "gromacs/legacyheaders/chargegroup.h"
56 #include "gromacs/legacyheaders/constr.h"
57 #include "gromacs/legacyheaders/domdec_network.h"
58 #include "gromacs/legacyheaders/force.h"
59 #include "gromacs/legacyheaders/gmx_ga2la.h"
60 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
61 #include "gromacs/legacyheaders/gpu_utils.h"
62 #include "gromacs/legacyheaders/macros.h"
63 #include "gromacs/legacyheaders/mdatoms.h"
64 #include "gromacs/legacyheaders/mdrun.h"
65 #include "gromacs/legacyheaders/names.h"
66 #include "gromacs/legacyheaders/network.h"
67 #include "gromacs/legacyheaders/nrnb.h"
68 #include "gromacs/legacyheaders/nsgrid.h"
69 #include "gromacs/legacyheaders/pme.h"
70 #include "gromacs/legacyheaders/shellfc.h"
71 #include "gromacs/legacyheaders/typedefs.h"
72 #include "gromacs/listed-forces/bonded.h"
73 #include "gromacs/math/vec.h"
74 #include "gromacs/mdlib/nb_verlet.h"
75 #include "gromacs/mdlib/nbnxn_search.h"
76 #include "gromacs/pbcutil/ishift.h"
77 #include "gromacs/pbcutil/pbc.h"
78 #include "gromacs/pulling/pull.h"
79 #include "gromacs/pulling/pull_rotation.h"
80 #include "gromacs/swap/swapcoords.h"
81 #include "gromacs/timing/wallcycle.h"
82 #include "gromacs/topology/mtop_util.h"
83 #include "gromacs/utility/basenetwork.h"
84 #include "gromacs/utility/fatalerror.h"
85 #include "gromacs/utility/futil.h"
86 #include "gromacs/utility/gmxmpi.h"
87 #include "gromacs/utility/qsort_threadsafe.h"
88 #include "gromacs/utility/smalloc.h"
90 #define DDRANK(dd, rank) (rank)
91 #define DDMASTERRANK(dd) (dd->masterrank)
93 typedef struct gmx_domdec_master
95 /* The cell boundaries */
97 /* The global charge group division */
98 int *ncg; /* Number of home charge groups for each node */
99 int *index; /* Index of nnodes+1 into cg */
100 int *cg; /* Global charge group index */
101 int *nat; /* Number of home atoms for each node. */
102 int *ibuf; /* Buffer for communication */
103 rvec *vbuf; /* Buffer for state scattering and gathering */
104 } gmx_domdec_master_t;
108 /* The numbers of charge groups to send and receive for each cell
109 * that requires communication, the last entry contains the total
110 * number of atoms that needs to be communicated.
112 int nsend[DD_MAXIZONE+2];
113 int nrecv[DD_MAXIZONE+2];
114 /* The charge groups to send */
117 /* The atom range for non-in-place communication */
118 int cell2at0[DD_MAXIZONE];
119 int cell2at1[DD_MAXIZONE];
124 int np; /* Number of grid pulses in this dimension */
125 int np_dlb; /* For dlb, for use with edlbAUTO */
126 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
128 gmx_bool bInPlace; /* Can we communicate in place? */
129 } gmx_domdec_comm_dim_t;
133 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
134 real *cell_f; /* State var.: cell boundaries, box relative */
135 real *old_cell_f; /* Temp. var.: old cell size */
136 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
137 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
138 real *bound_min; /* Temp. var.: lower limit for cell boundary */
139 real *bound_max; /* Temp. var.: upper limit for cell boundary */
140 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
141 real *buf_ncd; /* Temp. var. */
144 #define DD_NLOAD_MAX 9
146 /* Here floats are accurate enough, since these variables
147 * only influence the load balancing, not the actual MD results.
174 gmx_cgsort_t *sort_new;
186 /* This enum determines the order of the coordinates.
187 * ddnatHOME and ddnatZONE should be first and second,
188 * the others can be ordered as wanted.
191 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
195 edlbAUTO, edlbNO, edlbYES, edlbNR
197 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
201 int dim; /* The dimension */
202 gmx_bool dim_match; /* Tells if DD and PME dims match */
203 int nslab; /* The number of PME slabs in this dimension */
204 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
205 int *pp_min; /* The minimum pp node location, size nslab */
206 int *pp_max; /* The maximum pp node location,size nslab */
207 int maxshift; /* The maximum shift for coordinate redistribution in PME */
212 real min0; /* The minimum bottom of this zone */
213 real max1; /* The maximum top of this zone */
214 real min1; /* The minimum top of this zone */
215 real mch0; /* The maximum bottom communicaton height for this zone */
216 real mch1; /* The maximum top communicaton height for this zone */
217 real p1_0; /* The bottom value of the first cell in this zone */
218 real p1_1; /* The top value of the first cell in this zone */
223 gmx_domdec_ind_t ind;
230 } dd_comm_setup_work_t;
232 typedef struct gmx_domdec_comm
234 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
235 * unless stated otherwise.
238 /* The number of decomposition dimensions for PME, 0: no PME */
240 /* The number of nodes doing PME (PP/PME or only PME) */
244 /* The communication setup including the PME only nodes */
245 gmx_bool bCartesianPP_PME;
248 int *pmenodes; /* size npmenodes */
249 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
250 * but with bCartesianPP_PME */
251 gmx_ddpme_t ddpme[2];
253 /* The DD particle-particle nodes only */
254 gmx_bool bCartesianPP;
255 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
257 /* The global charge groups */
260 /* Should we sort the cgs */
262 gmx_domdec_sort_t *sort;
264 /* Are there charge groups? */
267 /* Are there bonded and multi-body interactions between charge groups? */
268 gmx_bool bInterCGBondeds;
269 gmx_bool bInterCGMultiBody;
271 /* Data for the optional bonded interaction atom communication range */
278 /* Is eDLB=edlbAUTO locked such that we currently can't turn it on? */
279 gmx_bool bDLB_locked;
280 /* Are we actually using DLB? */
281 gmx_bool bDynLoadBal;
283 /* Cell sizes for static load balancing, first index cartesian */
286 /* The width of the communicated boundaries */
289 /* The minimum cell size (including triclinic correction) */
291 /* For dlb, for use with edlbAUTO */
292 rvec cellsize_min_dlb;
293 /* The lower limit for the DD cell size with DLB */
295 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
296 gmx_bool bVacDLBNoLimit;
298 /* With PME load balancing we set limits on DLB */
299 gmx_bool bPMELoadBalDLBLimits;
300 /* DLB needs to take into account that we want to allow this maximum
301 * cut-off (for PME load balancing), this could limit cell boundaries.
303 real PMELoadBal_max_cutoff;
305 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
307 /* box0 and box_size are required with dim's without pbc and -gcom */
311 /* The cell boundaries */
315 /* The old location of the cell boundaries, to check cg displacements */
319 /* The communication setup and charge group boundaries for the zones */
320 gmx_domdec_zones_t zones;
322 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
323 * cell boundaries of neighboring cells for dynamic load balancing.
325 gmx_ddzone_t zone_d1[2];
326 gmx_ddzone_t zone_d2[2][2];
328 /* The coordinate/force communication setup and indices */
329 gmx_domdec_comm_dim_t cd[DIM];
330 /* The maximum number of cells to communicate with in one dimension */
333 /* Which cg distribution is stored on the master node */
334 int master_cg_ddp_count;
336 /* The number of cg's received from the direct neighbors */
337 int zone_ncg1[DD_MAXZONE];
339 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
342 /* Array for signalling if atoms have moved to another domain */
346 /* Communication buffer for general use */
350 /* Communication buffer for general use */
353 /* Temporary storage for thread parallel communication setup */
355 dd_comm_setup_work_t *dth;
357 /* Communication buffers only used with multiple grid pulses */
362 /* Communication buffers for local redistribution */
364 int cggl_flag_nalloc[DIM*2];
366 int cgcm_state_nalloc[DIM*2];
368 /* Cell sizes for dynamic load balancing */
369 gmx_domdec_root_t **root;
373 real cell_f_max0[DIM];
374 real cell_f_min1[DIM];
376 /* Stuff for load communication */
377 gmx_bool bRecordLoad;
378 gmx_domdec_load_t *load;
379 int nrank_gpu_shared;
381 MPI_Comm *mpi_comm_load;
382 MPI_Comm mpi_comm_gpu_shared;
385 /* Maximum DLB scaling per load balancing step in percent */
389 float cycl[ddCyclNr];
390 int cycl_n[ddCyclNr];
391 float cycl_max[ddCyclNr];
392 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
396 /* How many times have did we have load measurements */
398 /* How many times have we collected the load measurements */
402 double sum_nat[ddnatNR-ddnatZONE];
412 /* The last partition step */
413 gmx_int64_t partition_step;
421 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
424 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
425 #define DD_FLAG_NRCG 65535
426 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
427 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
429 /* Zone permutation required to obtain consecutive charge groups
430 * for neighbor searching.
432 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
434 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
435 * components see only j zones with that component 0.
438 /* The DD zone order */
439 static const ivec dd_zo[DD_MAXZONE] =
440 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
445 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
450 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
455 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
457 /* Factors used to avoid problems due to rounding issues */
458 #define DD_CELL_MARGIN 1.0001
459 #define DD_CELL_MARGIN2 1.00005
460 /* Factor to account for pressure scaling during nstlist steps */
461 #define DD_PRES_SCALE_MARGIN 1.02
463 /* Turn on DLB when the load imbalance causes this amount of total loss.
464 * There is a bit of overhead with DLB and it's difficult to achieve
465 * a load imbalance of less than 2% with DLB.
467 #define DD_PERF_LOSS_DLB_ON 0.02
469 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
470 #define DD_PERF_LOSS_WARN 0.05
472 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
474 /* Use separate MPI send and receive commands
475 * when nnodes <= GMX_DD_NNODES_SENDRECV.
476 * This saves memory (and some copying for small nnodes).
477 * For high parallelization scatter and gather calls are used.
479 #define GMX_DD_NNODES_SENDRECV 4
483 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
485 static void index2xyz(ivec nc,int ind,ivec xyz)
487 xyz[XX] = ind % nc[XX];
488 xyz[YY] = (ind / nc[XX]) % nc[YY];
489 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
493 /* This order is required to minimize the coordinate communication in PME
494 * which uses decomposition in the x direction.
496 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
498 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
500 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
501 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
502 xyz[ZZ] = ind % nc[ZZ];
505 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
510 ddindex = dd_index(dd->nc, c);
511 if (dd->comm->bCartesianPP_PME)
513 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
515 else if (dd->comm->bCartesianPP)
518 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
529 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
531 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
534 int ddglatnr(gmx_domdec_t *dd, int i)
544 if (i >= dd->comm->nat[ddnatNR-1])
546 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
548 atnr = dd->gatindex[i] + 1;
554 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
556 return &dd->comm->cgs_gl;
559 static void vec_rvec_init(vec_rvec_t *v)
565 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
569 v->nalloc = over_alloc_dd(n);
570 srenew(v->v, v->nalloc);
574 void dd_store_state(gmx_domdec_t *dd, t_state *state)
578 if (state->ddp_count != dd->ddp_count)
580 gmx_incons("The state does not the domain decomposition state");
583 state->ncg_gl = dd->ncg_home;
584 if (state->ncg_gl > state->cg_gl_nalloc)
586 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
587 srenew(state->cg_gl, state->cg_gl_nalloc);
589 for (i = 0; i < state->ncg_gl; i++)
591 state->cg_gl[i] = dd->index_gl[i];
594 state->ddp_count_cg_gl = dd->ddp_count;
597 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
599 return &dd->comm->zones;
602 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
603 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
605 gmx_domdec_zones_t *zones;
608 zones = &dd->comm->zones;
611 while (icg >= zones->izone[izone].cg1)
620 else if (izone < zones->nizone)
622 *jcg0 = zones->izone[izone].jcg0;
626 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
627 icg, izone, zones->nizone);
630 *jcg1 = zones->izone[izone].jcg1;
632 for (d = 0; d < dd->ndim; d++)
635 shift0[dim] = zones->izone[izone].shift0[dim];
636 shift1[dim] = zones->izone[izone].shift1[dim];
637 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
639 /* A conservative approach, this can be optimized */
646 int dd_natoms_vsite(gmx_domdec_t *dd)
648 return dd->comm->nat[ddnatVSITE];
651 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
653 *at_start = dd->comm->nat[ddnatCON-1];
654 *at_end = dd->comm->nat[ddnatCON];
657 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
659 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
660 int *index, *cgindex;
661 gmx_domdec_comm_t *comm;
662 gmx_domdec_comm_dim_t *cd;
663 gmx_domdec_ind_t *ind;
664 rvec shift = {0, 0, 0}, *buf, *rbuf;
665 gmx_bool bPBC, bScrew;
669 cgindex = dd->cgindex;
674 nat_tot = dd->nat_home;
675 for (d = 0; d < dd->ndim; d++)
677 bPBC = (dd->ci[dd->dim[d]] == 0);
678 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
681 copy_rvec(box[dd->dim[d]], shift);
684 for (p = 0; p < cd->np; p++)
691 for (i = 0; i < ind->nsend[nzone]; i++)
693 at0 = cgindex[index[i]];
694 at1 = cgindex[index[i]+1];
695 for (j = at0; j < at1; j++)
697 copy_rvec(x[j], buf[n]);
704 for (i = 0; i < ind->nsend[nzone]; i++)
706 at0 = cgindex[index[i]];
707 at1 = cgindex[index[i]+1];
708 for (j = at0; j < at1; j++)
710 /* We need to shift the coordinates */
711 rvec_add(x[j], shift, buf[n]);
718 for (i = 0; i < ind->nsend[nzone]; i++)
720 at0 = cgindex[index[i]];
721 at1 = cgindex[index[i]+1];
722 for (j = at0; j < at1; j++)
725 buf[n][XX] = x[j][XX] + shift[XX];
727 * This operation requires a special shift force
728 * treatment, which is performed in calc_vir.
730 buf[n][YY] = box[YY][YY] - x[j][YY];
731 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
743 rbuf = comm->vbuf2.v;
745 /* Send and receive the coordinates */
746 dd_sendrecv_rvec(dd, d, dddirBackward,
747 buf, ind->nsend[nzone+1],
748 rbuf, ind->nrecv[nzone+1]);
752 for (zone = 0; zone < nzone; zone++)
754 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
756 copy_rvec(rbuf[j], x[i]);
761 nat_tot += ind->nrecv[nzone+1];
767 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
769 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
770 int *index, *cgindex;
771 gmx_domdec_comm_t *comm;
772 gmx_domdec_comm_dim_t *cd;
773 gmx_domdec_ind_t *ind;
777 gmx_bool bShiftForcesNeedPbc, bScrew;
781 cgindex = dd->cgindex;
785 nzone = comm->zones.n/2;
786 nat_tot = dd->nat_tot;
787 for (d = dd->ndim-1; d >= 0; d--)
789 /* Only forces in domains near the PBC boundaries need to
790 consider PBC in the treatment of fshift */
791 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
792 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
793 if (fshift == NULL && !bScrew)
795 bShiftForcesNeedPbc = FALSE;
797 /* Determine which shift vector we need */
803 for (p = cd->np-1; p >= 0; p--)
806 nat_tot -= ind->nrecv[nzone+1];
813 sbuf = comm->vbuf2.v;
815 for (zone = 0; zone < nzone; zone++)
817 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
819 copy_rvec(f[i], sbuf[j]);
824 /* Communicate the forces */
825 dd_sendrecv_rvec(dd, d, dddirForward,
826 sbuf, ind->nrecv[nzone+1],
827 buf, ind->nsend[nzone+1]);
829 /* Add the received forces */
831 if (!bShiftForcesNeedPbc)
833 for (i = 0; i < ind->nsend[nzone]; i++)
835 at0 = cgindex[index[i]];
836 at1 = cgindex[index[i]+1];
837 for (j = at0; j < at1; j++)
839 rvec_inc(f[j], buf[n]);
846 /* fshift should always be defined if this function is
847 * called when bShiftForcesNeedPbc is true */
848 assert(NULL != fshift);
849 for (i = 0; i < ind->nsend[nzone]; i++)
851 at0 = cgindex[index[i]];
852 at1 = cgindex[index[i]+1];
853 for (j = at0; j < at1; j++)
855 rvec_inc(f[j], buf[n]);
856 /* Add this force to the shift force */
857 rvec_inc(fshift[is], buf[n]);
864 for (i = 0; i < ind->nsend[nzone]; i++)
866 at0 = cgindex[index[i]];
867 at1 = cgindex[index[i]+1];
868 for (j = at0; j < at1; j++)
870 /* Rotate the force */
871 f[j][XX] += buf[n][XX];
872 f[j][YY] -= buf[n][YY];
873 f[j][ZZ] -= buf[n][ZZ];
876 /* Add this force to the shift force */
877 rvec_inc(fshift[is], buf[n]);
888 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
890 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
891 int *index, *cgindex;
892 gmx_domdec_comm_t *comm;
893 gmx_domdec_comm_dim_t *cd;
894 gmx_domdec_ind_t *ind;
899 cgindex = dd->cgindex;
901 buf = &comm->vbuf.v[0][0];
904 nat_tot = dd->nat_home;
905 for (d = 0; d < dd->ndim; d++)
908 for (p = 0; p < cd->np; p++)
913 for (i = 0; i < ind->nsend[nzone]; i++)
915 at0 = cgindex[index[i]];
916 at1 = cgindex[index[i]+1];
917 for (j = at0; j < at1; j++)
930 rbuf = &comm->vbuf2.v[0][0];
932 /* Send and receive the coordinates */
933 dd_sendrecv_real(dd, d, dddirBackward,
934 buf, ind->nsend[nzone+1],
935 rbuf, ind->nrecv[nzone+1]);
939 for (zone = 0; zone < nzone; zone++)
941 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
948 nat_tot += ind->nrecv[nzone+1];
954 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
956 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
957 int *index, *cgindex;
958 gmx_domdec_comm_t *comm;
959 gmx_domdec_comm_dim_t *cd;
960 gmx_domdec_ind_t *ind;
965 cgindex = dd->cgindex;
967 buf = &comm->vbuf.v[0][0];
969 nzone = comm->zones.n/2;
970 nat_tot = dd->nat_tot;
971 for (d = dd->ndim-1; d >= 0; d--)
974 for (p = cd->np-1; p >= 0; p--)
977 nat_tot -= ind->nrecv[nzone+1];
984 sbuf = &comm->vbuf2.v[0][0];
986 for (zone = 0; zone < nzone; zone++)
988 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
995 /* Communicate the forces */
996 dd_sendrecv_real(dd, d, dddirForward,
997 sbuf, ind->nrecv[nzone+1],
998 buf, ind->nsend[nzone+1]);
1000 /* Add the received forces */
1002 for (i = 0; i < ind->nsend[nzone]; i++)
1004 at0 = cgindex[index[i]];
1005 at1 = cgindex[index[i]+1];
1006 for (j = at0; j < at1; j++)
1017 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1019 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",
1021 zone->min0, zone->max1,
1022 zone->mch0, zone->mch0,
1023 zone->p1_0, zone->p1_1);
1027 #define DDZONECOMM_MAXZONE 5
1028 #define DDZONECOMM_BUFSIZE 3
1030 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1031 int ddimind, int direction,
1032 gmx_ddzone_t *buf_s, int n_s,
1033 gmx_ddzone_t *buf_r, int n_r)
1035 #define ZBS DDZONECOMM_BUFSIZE
1036 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1037 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1040 for (i = 0; i < n_s; i++)
1042 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1043 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1044 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1045 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1046 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1047 vbuf_s[i*ZBS+1][2] = 0;
1048 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1049 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1050 vbuf_s[i*ZBS+2][2] = 0;
1053 dd_sendrecv_rvec(dd, ddimind, direction,
1057 for (i = 0; i < n_r; i++)
1059 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1060 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1061 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1062 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1063 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1064 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1065 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1071 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1072 rvec cell_ns_x0, rvec cell_ns_x1)
1074 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1076 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1077 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1078 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1079 rvec extr_s[2], extr_r[2];
1081 real dist_d, c = 0, det;
1082 gmx_domdec_comm_t *comm;
1083 gmx_bool bPBC, bUse;
1087 for (d = 1; d < dd->ndim; d++)
1090 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1091 zp->min0 = cell_ns_x0[dim];
1092 zp->max1 = cell_ns_x1[dim];
1093 zp->min1 = cell_ns_x1[dim];
1094 zp->mch0 = cell_ns_x0[dim];
1095 zp->mch1 = cell_ns_x1[dim];
1096 zp->p1_0 = cell_ns_x0[dim];
1097 zp->p1_1 = cell_ns_x1[dim];
1100 for (d = dd->ndim-2; d >= 0; d--)
1103 bPBC = (dim < ddbox->npbcdim);
1105 /* Use an rvec to store two reals */
1106 extr_s[d][0] = comm->cell_f0[d+1];
1107 extr_s[d][1] = comm->cell_f1[d+1];
1108 extr_s[d][2] = comm->cell_f1[d+1];
1111 /* Store the extremes in the backward sending buffer,
1112 * so the get updated separately from the forward communication.
1114 for (d1 = d; d1 < dd->ndim-1; d1++)
1116 /* We invert the order to be able to use the same loop for buf_e */
1117 buf_s[pos].min0 = extr_s[d1][1];
1118 buf_s[pos].max1 = extr_s[d1][0];
1119 buf_s[pos].min1 = extr_s[d1][2];
1120 buf_s[pos].mch0 = 0;
1121 buf_s[pos].mch1 = 0;
1122 /* Store the cell corner of the dimension we communicate along */
1123 buf_s[pos].p1_0 = comm->cell_x0[dim];
1124 buf_s[pos].p1_1 = 0;
1128 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1131 if (dd->ndim == 3 && d == 0)
1133 buf_s[pos] = comm->zone_d2[0][1];
1135 buf_s[pos] = comm->zone_d1[0];
1139 /* We only need to communicate the extremes
1140 * in the forward direction
1142 npulse = comm->cd[d].np;
1145 /* Take the minimum to avoid double communication */
1146 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1150 /* Without PBC we should really not communicate over
1151 * the boundaries, but implementing that complicates
1152 * the communication setup and therefore we simply
1153 * do all communication, but ignore some data.
1155 npulse_min = npulse;
1157 for (p = 0; p < npulse_min; p++)
1159 /* Communicate the extremes forward */
1160 bUse = (bPBC || dd->ci[dim] > 0);
1162 dd_sendrecv_rvec(dd, d, dddirForward,
1163 extr_s+d, dd->ndim-d-1,
1164 extr_r+d, dd->ndim-d-1);
1168 for (d1 = d; d1 < dd->ndim-1; d1++)
1170 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1171 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1172 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1178 for (p = 0; p < npulse; p++)
1180 /* Communicate all the zone information backward */
1181 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1183 dd_sendrecv_ddzone(dd, d, dddirBackward,
1190 for (d1 = d+1; d1 < dd->ndim; d1++)
1192 /* Determine the decrease of maximum required
1193 * communication height along d1 due to the distance along d,
1194 * this avoids a lot of useless atom communication.
1196 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1198 if (ddbox->tric_dir[dim])
1200 /* c is the off-diagonal coupling between the cell planes
1201 * along directions d and d1.
1203 c = ddbox->v[dim][dd->dim[d1]][dim];
1209 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1212 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1216 /* A negative value signals out of range */
1222 /* Accumulate the extremes over all pulses */
1223 for (i = 0; i < buf_size; i++)
1227 buf_e[i] = buf_r[i];
1233 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1234 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1235 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1238 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1246 if (bUse && dh[d1] >= 0)
1248 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1249 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1252 /* Copy the received buffer to the send buffer,
1253 * to pass the data through with the next pulse.
1255 buf_s[i] = buf_r[i];
1257 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1258 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1260 /* Store the extremes */
1263 for (d1 = d; d1 < dd->ndim-1; d1++)
1265 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1266 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1267 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1271 if (d == 1 || (d == 0 && dd->ndim == 3))
1273 for (i = d; i < 2; i++)
1275 comm->zone_d2[1-d][i] = buf_e[pos];
1281 comm->zone_d1[1] = buf_e[pos];
1291 for (i = 0; i < 2; i++)
1295 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1297 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1298 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1304 for (i = 0; i < 2; i++)
1306 for (j = 0; j < 2; j++)
1310 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1312 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1313 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1317 for (d = 1; d < dd->ndim; d++)
1319 comm->cell_f_max0[d] = extr_s[d-1][0];
1320 comm->cell_f_min1[d] = extr_s[d-1][1];
1323 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1324 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1329 static void dd_collect_cg(gmx_domdec_t *dd,
1330 t_state *state_local)
1332 gmx_domdec_master_t *ma = NULL;
1333 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1335 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1337 /* The master has the correct distribution */
1341 if (state_local->ddp_count == dd->ddp_count)
1343 /* The local state and DD are in sync, use the DD indices */
1344 ncg_home = dd->ncg_home;
1346 nat_home = dd->nat_home;
1348 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1350 /* The DD is out of sync with the local state, but we have stored
1351 * the cg indices with the local state, so we can use those.
1355 cgs_gl = &dd->comm->cgs_gl;
1357 ncg_home = state_local->ncg_gl;
1358 cg = state_local->cg_gl;
1360 for (i = 0; i < ncg_home; i++)
1362 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1367 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1381 /* Collect the charge group and atom counts on the master */
1382 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1387 for (i = 0; i < dd->nnodes; i++)
1389 ma->ncg[i] = ma->ibuf[2*i];
1390 ma->nat[i] = ma->ibuf[2*i+1];
1391 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1394 /* Make byte counts and indices */
1395 for (i = 0; i < dd->nnodes; i++)
1397 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1398 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1402 fprintf(debug, "Initial charge group distribution: ");
1403 for (i = 0; i < dd->nnodes; i++)
1405 fprintf(debug, " %d", ma->ncg[i]);
1407 fprintf(debug, "\n");
1411 /* Collect the charge group indices on the master */
1413 ncg_home*sizeof(int), cg,
1414 DDMASTER(dd) ? ma->ibuf : NULL,
1415 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1416 DDMASTER(dd) ? ma->cg : NULL);
1418 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1421 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1424 gmx_domdec_master_t *ma;
1425 int n, i, c, a, nalloc = 0;
1434 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1435 dd->rank, dd->mpi_comm_all);
1440 /* Copy the master coordinates to the global array */
1441 cgs_gl = &dd->comm->cgs_gl;
1443 n = DDMASTERRANK(dd);
1445 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1447 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1449 copy_rvec(lv[a++], v[c]);
1453 for (n = 0; n < dd->nnodes; n++)
1457 if (ma->nat[n] > nalloc)
1459 nalloc = over_alloc_dd(ma->nat[n]);
1460 srenew(buf, nalloc);
1463 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1464 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1467 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1469 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1471 copy_rvec(buf[a++], v[c]);
1480 static void get_commbuffer_counts(gmx_domdec_t *dd,
1481 int **counts, int **disps)
1483 gmx_domdec_master_t *ma;
1488 /* Make the rvec count and displacment arrays */
1490 *disps = ma->ibuf + dd->nnodes;
1491 for (n = 0; n < dd->nnodes; n++)
1493 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1494 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1498 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1501 gmx_domdec_master_t *ma;
1502 int *rcounts = NULL, *disps = NULL;
1511 get_commbuffer_counts(dd, &rcounts, &disps);
1516 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1520 cgs_gl = &dd->comm->cgs_gl;
1523 for (n = 0; n < dd->nnodes; n++)
1525 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1527 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1529 copy_rvec(buf[a++], v[c]);
1536 void dd_collect_vec(gmx_domdec_t *dd,
1537 t_state *state_local, rvec *lv, rvec *v)
1539 dd_collect_cg(dd, state_local);
1541 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1543 dd_collect_vec_sendrecv(dd, lv, v);
1547 dd_collect_vec_gatherv(dd, lv, v);
1552 void dd_collect_state(gmx_domdec_t *dd,
1553 t_state *state_local, t_state *state)
1557 nh = state->nhchainlength;
1561 for (i = 0; i < efptNR; i++)
1563 state->lambda[i] = state_local->lambda[i];
1565 state->fep_state = state_local->fep_state;
1566 state->veta = state_local->veta;
1567 state->vol0 = state_local->vol0;
1568 copy_mat(state_local->box, state->box);
1569 copy_mat(state_local->boxv, state->boxv);
1570 copy_mat(state_local->svir_prev, state->svir_prev);
1571 copy_mat(state_local->fvir_prev, state->fvir_prev);
1572 copy_mat(state_local->pres_prev, state->pres_prev);
1574 for (i = 0; i < state_local->ngtc; i++)
1576 for (j = 0; j < nh; j++)
1578 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1579 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1581 state->therm_integral[i] = state_local->therm_integral[i];
1583 for (i = 0; i < state_local->nnhpres; i++)
1585 for (j = 0; j < nh; j++)
1587 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1588 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1592 for (est = 0; est < estNR; est++)
1594 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1599 dd_collect_vec(dd, state_local, state_local->x, state->x);
1602 dd_collect_vec(dd, state_local, state_local->v, state->v);
1605 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1608 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1610 case estDISRE_INITF:
1611 case estDISRE_RM3TAV:
1612 case estORIRE_INITF:
1616 gmx_incons("Unknown state entry encountered in dd_collect_state");
1622 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1628 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1631 state->nalloc = over_alloc_dd(nalloc);
1633 for (est = 0; est < estNR; est++)
1635 if (EST_DISTR(est) && (state->flags & (1<<est)))
1640 srenew(state->x, state->nalloc);
1643 srenew(state->v, state->nalloc);
1646 srenew(state->sd_X, state->nalloc);
1649 srenew(state->cg_p, state->nalloc);
1651 case estDISRE_INITF:
1652 case estDISRE_RM3TAV:
1653 case estORIRE_INITF:
1655 /* No reallocation required */
1658 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1665 srenew(*f, state->nalloc);
1669 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1672 if (nalloc > fr->cg_nalloc)
1676 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1678 fr->cg_nalloc = over_alloc_dd(nalloc);
1679 srenew(fr->cginfo, fr->cg_nalloc);
1680 if (fr->cutoff_scheme == ecutsGROUP)
1682 srenew(fr->cg_cm, fr->cg_nalloc);
1685 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1687 /* We don't use charge groups, we use x in state to set up
1688 * the atom communication.
1690 dd_realloc_state(state, f, nalloc);
1694 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1697 gmx_domdec_master_t *ma;
1698 int n, i, c, a, nalloc = 0;
1705 for (n = 0; n < dd->nnodes; n++)
1709 if (ma->nat[n] > nalloc)
1711 nalloc = over_alloc_dd(ma->nat[n]);
1712 srenew(buf, nalloc);
1714 /* Use lv as a temporary buffer */
1716 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1718 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1720 copy_rvec(v[c], buf[a++]);
1723 if (a != ma->nat[n])
1725 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1730 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1731 DDRANK(dd, n), n, dd->mpi_comm_all);
1736 n = DDMASTERRANK(dd);
1738 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1740 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1742 copy_rvec(v[c], lv[a++]);
1749 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1750 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1755 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1758 gmx_domdec_master_t *ma;
1759 int *scounts = NULL, *disps = NULL;
1767 get_commbuffer_counts(dd, &scounts, &disps);
1771 for (n = 0; n < dd->nnodes; n++)
1773 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1775 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1777 copy_rvec(v[c], buf[a++]);
1783 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1786 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1788 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1790 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1794 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1798 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1801 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1802 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1803 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1805 if (dfhist->nlambda > 0)
1807 int nlam = dfhist->nlambda;
1808 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1809 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1810 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1811 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1812 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1813 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1815 for (i = 0; i < nlam; i++)
1817 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1818 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1819 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1820 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1821 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1822 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1827 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1828 t_state *state, t_state *state_local,
1833 nh = state->nhchainlength;
1837 for (i = 0; i < efptNR; i++)
1839 state_local->lambda[i] = state->lambda[i];
1841 state_local->fep_state = state->fep_state;
1842 state_local->veta = state->veta;
1843 state_local->vol0 = state->vol0;
1844 copy_mat(state->box, state_local->box);
1845 copy_mat(state->box_rel, state_local->box_rel);
1846 copy_mat(state->boxv, state_local->boxv);
1847 copy_mat(state->svir_prev, state_local->svir_prev);
1848 copy_mat(state->fvir_prev, state_local->fvir_prev);
1849 copy_df_history(&state_local->dfhist, &state->dfhist);
1850 for (i = 0; i < state_local->ngtc; i++)
1852 for (j = 0; j < nh; j++)
1854 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1855 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1857 state_local->therm_integral[i] = state->therm_integral[i];
1859 for (i = 0; i < state_local->nnhpres; i++)
1861 for (j = 0; j < nh; j++)
1863 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1864 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1868 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1869 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1870 dd_bcast(dd, sizeof(real), &state_local->veta);
1871 dd_bcast(dd, sizeof(real), &state_local->vol0);
1872 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1873 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1874 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1875 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1876 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1877 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1878 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1879 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1880 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1881 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1883 /* communicate df_history -- required for restarting from checkpoint */
1884 dd_distribute_dfhist(dd, &state_local->dfhist);
1886 if (dd->nat_home > state_local->nalloc)
1888 dd_realloc_state(state_local, f, dd->nat_home);
1890 for (i = 0; i < estNR; i++)
1892 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1897 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1900 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1903 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1906 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1908 case estDISRE_INITF:
1909 case estDISRE_RM3TAV:
1910 case estORIRE_INITF:
1912 /* Not implemented yet */
1915 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1921 static char dim2char(int dim)
1927 case XX: c = 'X'; break;
1928 case YY: c = 'Y'; break;
1929 case ZZ: c = 'Z'; break;
1930 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1936 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1937 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1939 rvec grid_s[2], *grid_r = NULL, cx, r;
1940 char fname[STRLEN], buf[22];
1942 int a, i, d, z, y, x;
1946 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1947 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1951 snew(grid_r, 2*dd->nnodes);
1954 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1958 for (d = 0; d < DIM; d++)
1960 for (i = 0; i < DIM; i++)
1968 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1970 tric[d][i] = box[i][d]/box[i][i];
1979 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1980 out = gmx_fio_fopen(fname, "w");
1981 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1983 for (i = 0; i < dd->nnodes; i++)
1985 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1986 for (d = 0; d < DIM; d++)
1988 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1990 for (z = 0; z < 2; z++)
1992 for (y = 0; y < 2; y++)
1994 for (x = 0; x < 2; x++)
1996 cx[XX] = grid_r[i*2+x][XX];
1997 cx[YY] = grid_r[i*2+y][YY];
1998 cx[ZZ] = grid_r[i*2+z][ZZ];
2000 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2001 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2005 for (d = 0; d < DIM; d++)
2007 for (x = 0; x < 4; x++)
2011 case 0: y = 1 + i*8 + 2*x; break;
2012 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2013 case 2: y = 1 + i*8 + x; break;
2015 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2019 gmx_fio_fclose(out);
2024 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2025 gmx_mtop_t *mtop, t_commrec *cr,
2026 int natoms, rvec x[], matrix box)
2028 char fname[STRLEN], buf[22];
2030 int i, ii, resnr, c;
2031 char *atomname, *resname;
2038 natoms = dd->comm->nat[ddnatVSITE];
2041 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2043 out = gmx_fio_fopen(fname, "w");
2045 fprintf(out, "TITLE %s\n", title);
2046 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2047 for (i = 0; i < natoms; i++)
2049 ii = dd->gatindex[i];
2050 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2051 if (i < dd->comm->nat[ddnatZONE])
2054 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2060 else if (i < dd->comm->nat[ddnatVSITE])
2062 b = dd->comm->zones.n;
2066 b = dd->comm->zones.n + 1;
2068 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2069 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2071 fprintf(out, "TER\n");
2073 gmx_fio_fclose(out);
2076 real dd_cutoff_mbody(gmx_domdec_t *dd)
2078 gmx_domdec_comm_t *comm;
2085 if (comm->bInterCGBondeds)
2087 if (comm->cutoff_mbody > 0)
2089 r = comm->cutoff_mbody;
2093 /* cutoff_mbody=0 means we do not have DLB */
2094 r = comm->cellsize_min[dd->dim[0]];
2095 for (di = 1; di < dd->ndim; di++)
2097 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2099 if (comm->bBondComm)
2101 r = std::max(r, comm->cutoff_mbody);
2105 r = std::min(r, comm->cutoff);
2113 real dd_cutoff_twobody(gmx_domdec_t *dd)
2117 r_mb = dd_cutoff_mbody(dd);
2119 return std::max(dd->comm->cutoff, r_mb);
2123 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2127 nc = dd->nc[dd->comm->cartpmedim];
2128 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2129 copy_ivec(coord, coord_pme);
2130 coord_pme[dd->comm->cartpmedim] =
2131 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2134 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2136 /* Here we assign a PME node to communicate with this DD node
2137 * by assuming that the major index of both is x.
2138 * We add cr->npmenodes/2 to obtain an even distribution.
2140 return (ddindex*npme + npme/2)/ndd;
2143 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2145 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2148 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2150 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2153 static int *dd_pmenodes(t_commrec *cr)
2158 snew(pmenodes, cr->npmenodes);
2160 for (i = 0; i < cr->dd->nnodes; i++)
2162 p0 = cr_ddindex2pmeindex(cr, i);
2163 p1 = cr_ddindex2pmeindex(cr, i+1);
2164 if (i+1 == cr->dd->nnodes || p1 > p0)
2168 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2170 pmenodes[n] = i + 1 + n;
2178 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2186 if (dd->comm->bCartesian) {
2187 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2188 dd_coords2pmecoords(dd,coords,coords_pme);
2189 copy_ivec(dd->ntot,nc);
2190 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2191 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2193 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2195 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2201 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2206 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2208 gmx_domdec_comm_t *comm;
2210 int ddindex, nodeid = -1;
2212 comm = cr->dd->comm;
2217 if (comm->bCartesianPP_PME)
2220 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2225 ddindex = dd_index(cr->dd->nc, coords);
2226 if (comm->bCartesianPP)
2228 nodeid = comm->ddindex2simnodeid[ddindex];
2234 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2246 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2249 gmx_domdec_comm_t *comm;
2256 /* This assumes a uniform x domain decomposition grid cell size */
2257 if (comm->bCartesianPP_PME)
2260 ivec coord, coord_pme;
2261 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2262 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2264 /* This is a PP node */
2265 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2266 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2270 else if (comm->bCartesianPP)
2272 if (sim_nodeid < dd->nnodes)
2274 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2279 /* This assumes DD cells with identical x coordinates
2280 * are numbered sequentially.
2282 if (dd->comm->pmenodes == NULL)
2284 if (sim_nodeid < dd->nnodes)
2286 /* The DD index equals the nodeid */
2287 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2293 while (sim_nodeid > dd->comm->pmenodes[i])
2297 if (sim_nodeid < dd->comm->pmenodes[i])
2299 pmenode = dd->comm->pmenodes[i];
2307 void get_pme_nnodes(const gmx_domdec_t *dd,
2308 int *npmenodes_x, int *npmenodes_y)
2312 *npmenodes_x = dd->comm->npmenodes_x;
2313 *npmenodes_y = dd->comm->npmenodes_y;
2322 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2324 gmx_bool bPMEOnlyNode;
2326 if (DOMAINDECOMP(cr))
2328 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2332 bPMEOnlyNode = FALSE;
2335 return bPMEOnlyNode;
2338 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2339 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2343 ivec coord, coord_pme;
2347 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2350 for (x = 0; x < dd->nc[XX]; x++)
2352 for (y = 0; y < dd->nc[YY]; y++)
2354 for (z = 0; z < dd->nc[ZZ]; z++)
2356 if (dd->comm->bCartesianPP_PME)
2361 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2362 if (dd->ci[XX] == coord_pme[XX] &&
2363 dd->ci[YY] == coord_pme[YY] &&
2364 dd->ci[ZZ] == coord_pme[ZZ])
2366 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2371 /* The slab corresponds to the nodeid in the PME group */
2372 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2374 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2381 /* The last PP-only node is the peer node */
2382 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2386 fprintf(debug, "Receive coordinates from PP ranks:");
2387 for (x = 0; x < *nmy_ddnodes; x++)
2389 fprintf(debug, " %d", (*my_ddnodes)[x]);
2391 fprintf(debug, "\n");
2395 static gmx_bool receive_vir_ener(t_commrec *cr)
2397 gmx_domdec_comm_t *comm;
2402 if (cr->npmenodes < cr->dd->nnodes)
2404 comm = cr->dd->comm;
2405 if (comm->bCartesianPP_PME)
2407 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2410 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2411 coords[comm->cartpmedim]++;
2412 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2415 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2416 if (dd_simnode2pmenode(cr, rank) == pmenode)
2418 /* This is not the last PP node for pmenode */
2426 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2427 if (cr->sim_nodeid+1 < cr->nnodes &&
2428 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2430 /* This is not the last PP node for pmenode */
2439 static void set_zones_ncg_home(gmx_domdec_t *dd)
2441 gmx_domdec_zones_t *zones;
2444 zones = &dd->comm->zones;
2446 zones->cg_range[0] = 0;
2447 for (i = 1; i < zones->n+1; i++)
2449 zones->cg_range[i] = dd->ncg_home;
2451 /* zone_ncg1[0] should always be equal to ncg_home */
2452 dd->comm->zone_ncg1[0] = dd->ncg_home;
2455 static void rebuild_cgindex(gmx_domdec_t *dd,
2456 const int *gcgs_index, t_state *state)
2458 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2461 dd_cg_gl = dd->index_gl;
2462 cgindex = dd->cgindex;
2465 for (i = 0; i < state->ncg_gl; i++)
2469 dd_cg_gl[i] = cg_gl;
2470 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2474 dd->ncg_home = state->ncg_gl;
2477 set_zones_ncg_home(dd);
2480 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2482 while (cg >= cginfo_mb->cg_end)
2487 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2490 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2491 t_forcerec *fr, char *bLocalCG)
2493 cginfo_mb_t *cginfo_mb;
2499 cginfo_mb = fr->cginfo_mb;
2500 cginfo = fr->cginfo;
2502 for (cg = cg0; cg < cg1; cg++)
2504 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2508 if (bLocalCG != NULL)
2510 for (cg = cg0; cg < cg1; cg++)
2512 bLocalCG[index_gl[cg]] = TRUE;
2517 static void make_dd_indices(gmx_domdec_t *dd,
2518 const int *gcgs_index, int cg_start)
2520 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2521 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2524 if (dd->nat_tot > dd->gatindex_nalloc)
2526 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2527 srenew(dd->gatindex, dd->gatindex_nalloc);
2530 nzone = dd->comm->zones.n;
2531 zone2cg = dd->comm->zones.cg_range;
2532 zone_ncg1 = dd->comm->zone_ncg1;
2533 index_gl = dd->index_gl;
2534 gatindex = dd->gatindex;
2535 bCGs = dd->comm->bCGs;
2537 if (zone2cg[1] != dd->ncg_home)
2539 gmx_incons("dd->ncg_zone is not up to date");
2542 /* Make the local to global and global to local atom index */
2543 a = dd->cgindex[cg_start];
2544 for (zone = 0; zone < nzone; zone++)
2552 cg0 = zone2cg[zone];
2554 cg1 = zone2cg[zone+1];
2555 cg1_p1 = cg0 + zone_ncg1[zone];
2557 for (cg = cg0; cg < cg1; cg++)
2562 /* Signal that this cg is from more than one pulse away */
2565 cg_gl = index_gl[cg];
2568 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2571 ga2la_set(dd->ga2la, a_gl, a, zone1);
2577 gatindex[a] = cg_gl;
2578 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2585 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2591 if (bLocalCG == NULL)
2595 for (i = 0; i < dd->ncg_tot; i++)
2597 if (!bLocalCG[dd->index_gl[i]])
2600 "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);
2605 for (i = 0; i < ncg_sys; i++)
2612 if (ngl != dd->ncg_tot)
2614 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);
2621 static void check_index_consistency(gmx_domdec_t *dd,
2622 int natoms_sys, int ncg_sys,
2625 int nerr, ngl, i, a, cell;
2630 if (dd->comm->DD_debug > 1)
2632 snew(have, natoms_sys);
2633 for (a = 0; a < dd->nat_tot; a++)
2635 if (have[dd->gatindex[a]] > 0)
2637 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);
2641 have[dd->gatindex[a]] = a + 1;
2647 snew(have, dd->nat_tot);
2650 for (i = 0; i < natoms_sys; i++)
2652 if (ga2la_get(dd->ga2la, i, &a, &cell))
2654 if (a >= dd->nat_tot)
2656 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);
2662 if (dd->gatindex[a] != i)
2664 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);
2671 if (ngl != dd->nat_tot)
2674 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2675 dd->rank, where, ngl, dd->nat_tot);
2677 for (a = 0; a < dd->nat_tot; a++)
2682 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2683 dd->rank, where, a+1, dd->gatindex[a]+1);
2688 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2692 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2693 dd->rank, where, nerr);
2697 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2704 /* Clear the whole list without searching */
2705 ga2la_clear(dd->ga2la);
2709 for (i = a_start; i < dd->nat_tot; i++)
2711 ga2la_del(dd->ga2la, dd->gatindex[i]);
2715 bLocalCG = dd->comm->bLocalCG;
2718 for (i = cg_start; i < dd->ncg_tot; i++)
2720 bLocalCG[dd->index_gl[i]] = FALSE;
2724 dd_clear_local_vsite_indices(dd);
2726 if (dd->constraints)
2728 dd_clear_local_constraint_indices(dd);
2732 /* This function should be used for moving the domain boudaries during DLB,
2733 * for obtaining the minimum cell size. It checks the initially set limit
2734 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2735 * and, possibly, a longer cut-off limit set for PME load balancing.
2737 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2741 cellsize_min = comm->cellsize_min[dim];
2743 if (!comm->bVacDLBNoLimit)
2745 /* The cut-off might have changed, e.g. by PME load balacning,
2746 * from the value used to set comm->cellsize_min, so check it.
2748 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2750 if (comm->bPMELoadBalDLBLimits)
2752 /* Check for the cut-off limit set by the PME load balancing */
2753 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2757 return cellsize_min;
2760 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2763 real grid_jump_limit;
2765 /* The distance between the boundaries of cells at distance
2766 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2767 * and by the fact that cells should not be shifted by more than
2768 * half their size, such that cg's only shift by one cell
2769 * at redecomposition.
2771 grid_jump_limit = comm->cellsize_limit;
2772 if (!comm->bVacDLBNoLimit)
2774 if (comm->bPMELoadBalDLBLimits)
2776 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2778 grid_jump_limit = std::max(grid_jump_limit,
2779 cutoff/comm->cd[dim_ind].np);
2782 return grid_jump_limit;
2785 static gmx_bool check_grid_jump(gmx_int64_t step,
2791 gmx_domdec_comm_t *comm;
2800 for (d = 1; d < dd->ndim; d++)
2803 limit = grid_jump_limit(comm, cutoff, d);
2804 bfac = ddbox->box_size[dim];
2805 if (ddbox->tric_dir[dim])
2807 bfac *= ddbox->skew_fac[dim];
2809 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2810 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2818 /* This error should never be triggered under normal
2819 * circumstances, but you never know ...
2821 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.",
2822 gmx_step_str(step, buf),
2823 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2831 static int dd_load_count(gmx_domdec_comm_t *comm)
2833 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2836 static float dd_force_load(gmx_domdec_comm_t *comm)
2843 if (comm->eFlop > 1)
2845 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2850 load = comm->cycl[ddCyclF];
2851 if (comm->cycl_n[ddCyclF] > 1)
2853 /* Subtract the maximum of the last n cycle counts
2854 * to get rid of possible high counts due to other sources,
2855 * for instance system activity, that would otherwise
2856 * affect the dynamic load balancing.
2858 load -= comm->cycl_max[ddCyclF];
2862 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2864 float gpu_wait, gpu_wait_sum;
2866 gpu_wait = comm->cycl[ddCyclWaitGPU];
2867 if (comm->cycl_n[ddCyclF] > 1)
2869 /* We should remove the WaitGPU time of the same MD step
2870 * as the one with the maximum F time, since the F time
2871 * and the wait time are not independent.
2872 * Furthermore, the step for the max F time should be chosen
2873 * the same on all ranks that share the same GPU.
2874 * But to keep the code simple, we remove the average instead.
2875 * The main reason for artificially long times at some steps
2876 * is spurious CPU activity or MPI time, so we don't expect
2877 * that changes in the GPU wait time matter a lot here.
2879 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2881 /* Sum the wait times over the ranks that share the same GPU */
2882 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2883 comm->mpi_comm_gpu_shared);
2884 /* Replace the wait time by the average over the ranks */
2885 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2893 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2895 gmx_domdec_comm_t *comm;
2900 snew(*dim_f, dd->nc[dim]+1);
2902 for (i = 1; i < dd->nc[dim]; i++)
2904 if (comm->slb_frac[dim])
2906 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2910 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2913 (*dim_f)[dd->nc[dim]] = 1;
2916 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2918 int pmeindex, slab, nso, i;
2921 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2927 ddpme->dim = dimind;
2929 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2931 ddpme->nslab = (ddpme->dim == 0 ?
2932 dd->comm->npmenodes_x :
2933 dd->comm->npmenodes_y);
2935 if (ddpme->nslab <= 1)
2940 nso = dd->comm->npmenodes/ddpme->nslab;
2941 /* Determine for each PME slab the PP location range for dimension dim */
2942 snew(ddpme->pp_min, ddpme->nslab);
2943 snew(ddpme->pp_max, ddpme->nslab);
2944 for (slab = 0; slab < ddpme->nslab; slab++)
2946 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2947 ddpme->pp_max[slab] = 0;
2949 for (i = 0; i < dd->nnodes; i++)
2951 ddindex2xyz(dd->nc, i, xyz);
2952 /* For y only use our y/z slab.
2953 * This assumes that the PME x grid size matches the DD grid size.
2955 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2957 pmeindex = ddindex2pmeindex(dd, i);
2960 slab = pmeindex/nso;
2964 slab = pmeindex % ddpme->nslab;
2966 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2967 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2971 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2974 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2976 if (dd->comm->ddpme[0].dim == XX)
2978 return dd->comm->ddpme[0].maxshift;
2986 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2988 if (dd->comm->ddpme[0].dim == YY)
2990 return dd->comm->ddpme[0].maxshift;
2992 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2994 return dd->comm->ddpme[1].maxshift;
3002 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3003 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3005 gmx_domdec_comm_t *comm;
3008 real range, pme_boundary;
3012 nc = dd->nc[ddpme->dim];
3015 if (!ddpme->dim_match)
3017 /* PP decomposition is not along dim: the worst situation */
3020 else if (ns <= 3 || (bUniform && ns == nc))
3022 /* The optimal situation */
3027 /* We need to check for all pme nodes which nodes they
3028 * could possibly need to communicate with.
3030 xmin = ddpme->pp_min;
3031 xmax = ddpme->pp_max;
3032 /* Allow for atoms to be maximally 2/3 times the cut-off
3033 * out of their DD cell. This is a reasonable balance between
3034 * between performance and support for most charge-group/cut-off
3037 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3038 /* Avoid extra communication when we are exactly at a boundary */
3042 for (s = 0; s < ns; s++)
3044 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3045 pme_boundary = (real)s/ns;
3048 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3050 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3054 pme_boundary = (real)(s+1)/ns;
3057 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3059 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3066 ddpme->maxshift = sh;
3070 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3071 ddpme->dim, ddpme->maxshift);
3075 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3079 for (d = 0; d < dd->ndim; d++)
3082 if (dim < ddbox->nboundeddim &&
3083 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3084 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3086 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",
3087 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3088 dd->nc[dim], dd->comm->cellsize_limit);
3094 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3097 /* Set the domain boundaries. Use for static (or no) load balancing,
3098 * and also for the starting state for dynamic load balancing.
3099 * setmode determine if and where the boundaries are stored, use enum above.
3100 * Returns the number communication pulses in npulse.
3102 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3103 int setmode, ivec npulse)
3105 gmx_domdec_comm_t *comm;
3108 real *cell_x, cell_dx, cellsize;
3112 for (d = 0; d < DIM; d++)
3114 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3116 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3119 cell_dx = ddbox->box_size[d]/dd->nc[d];
3122 case setcellsizeslbMASTER:
3123 for (j = 0; j < dd->nc[d]+1; j++)
3125 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3128 case setcellsizeslbLOCAL:
3129 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3130 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3135 cellsize = cell_dx*ddbox->skew_fac[d];
3136 while (cellsize*npulse[d] < comm->cutoff)
3140 cellsize_min[d] = cellsize;
3144 /* Statically load balanced grid */
3145 /* Also when we are not doing a master distribution we determine
3146 * all cell borders in a loop to obtain identical values
3147 * to the master distribution case and to determine npulse.
3149 if (setmode == setcellsizeslbMASTER)
3151 cell_x = dd->ma->cell_x[d];
3155 snew(cell_x, dd->nc[d]+1);
3157 cell_x[0] = ddbox->box0[d];
3158 for (j = 0; j < dd->nc[d]; j++)
3160 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3161 cell_x[j+1] = cell_x[j] + cell_dx;
3162 cellsize = cell_dx*ddbox->skew_fac[d];
3163 while (cellsize*npulse[d] < comm->cutoff &&
3164 npulse[d] < dd->nc[d]-1)
3168 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3170 if (setmode == setcellsizeslbLOCAL)
3172 comm->cell_x0[d] = cell_x[dd->ci[d]];
3173 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3175 if (setmode != setcellsizeslbMASTER)
3180 /* The following limitation is to avoid that a cell would receive
3181 * some of its own home charge groups back over the periodic boundary.
3182 * Double charge groups cause trouble with the global indices.
3184 if (d < ddbox->npbcdim &&
3185 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3187 char error_string[STRLEN];
3189 sprintf(error_string,
3190 "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",
3191 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3193 dd->nc[d], dd->nc[d],
3194 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3196 if (setmode == setcellsizeslbLOCAL)
3198 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3202 gmx_fatal(FARGS, error_string);
3207 if (!comm->bDynLoadBal)
3209 copy_rvec(cellsize_min, comm->cellsize_min);
3212 for (d = 0; d < comm->npmedecompdim; d++)
3214 set_pme_maxshift(dd, &comm->ddpme[d],
3215 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3216 comm->ddpme[d].slb_dim_f);
3221 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3222 int d, int dim, gmx_domdec_root_t *root,
3224 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3226 gmx_domdec_comm_t *comm;
3227 int ncd, i, j, nmin, nmin_old;
3228 gmx_bool bLimLo, bLimHi;
3230 real fac, halfway, cellsize_limit_f_i, region_size;
3231 gmx_bool bPBC, bLastHi = FALSE;
3232 int nrange[] = {range[0], range[1]};
3234 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3240 bPBC = (dim < ddbox->npbcdim);
3242 cell_size = root->buf_ncd;
3246 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3249 /* First we need to check if the scaling does not make cells
3250 * smaller than the smallest allowed size.
3251 * We need to do this iteratively, since if a cell is too small,
3252 * it needs to be enlarged, which makes all the other cells smaller,
3253 * which could in turn make another cell smaller than allowed.
3255 for (i = range[0]; i < range[1]; i++)
3257 root->bCellMin[i] = FALSE;
3263 /* We need the total for normalization */
3265 for (i = range[0]; i < range[1]; i++)
3267 if (root->bCellMin[i] == FALSE)
3269 fac += cell_size[i];
3272 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3273 /* Determine the cell boundaries */
3274 for (i = range[0]; i < range[1]; i++)
3276 if (root->bCellMin[i] == FALSE)
3278 cell_size[i] *= fac;
3279 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3281 cellsize_limit_f_i = 0;
3285 cellsize_limit_f_i = cellsize_limit_f;
3287 if (cell_size[i] < cellsize_limit_f_i)
3289 root->bCellMin[i] = TRUE;
3290 cell_size[i] = cellsize_limit_f_i;
3294 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3297 while (nmin > nmin_old);
3300 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3301 /* For this check we should not use DD_CELL_MARGIN,
3302 * but a slightly smaller factor,
3303 * since rounding could get use below the limit.
3305 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3308 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",
3309 gmx_step_str(step, buf),
3310 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3311 ncd, comm->cellsize_min[dim]);
3314 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3318 /* Check if the boundary did not displace more than halfway
3319 * each of the cells it bounds, as this could cause problems,
3320 * especially when the differences between cell sizes are large.
3321 * If changes are applied, they will not make cells smaller
3322 * than the cut-off, as we check all the boundaries which
3323 * might be affected by a change and if the old state was ok,
3324 * the cells will at most be shrunk back to their old size.
3326 for (i = range[0]+1; i < range[1]; i++)
3328 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3329 if (root->cell_f[i] < halfway)
3331 root->cell_f[i] = halfway;
3332 /* Check if the change also causes shifts of the next boundaries */
3333 for (j = i+1; j < range[1]; j++)
3335 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3337 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3341 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3342 if (root->cell_f[i] > halfway)
3344 root->cell_f[i] = halfway;
3345 /* Check if the change also causes shifts of the next boundaries */
3346 for (j = i-1; j >= range[0]+1; j--)
3348 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3350 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3357 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3358 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3359 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3360 * for a and b nrange is used */
3363 /* Take care of the staggering of the cell boundaries */
3366 for (i = range[0]; i < range[1]; i++)
3368 root->cell_f_max0[i] = root->cell_f[i];
3369 root->cell_f_min1[i] = root->cell_f[i+1];
3374 for (i = range[0]+1; i < range[1]; i++)
3376 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3377 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3378 if (bLimLo && bLimHi)
3380 /* Both limits violated, try the best we can */
3381 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3382 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3383 nrange[0] = range[0];
3385 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3388 nrange[1] = range[1];
3389 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3395 /* root->cell_f[i] = root->bound_min[i]; */
3396 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3399 else if (bLimHi && !bLastHi)
3402 if (nrange[1] < range[1]) /* found a LimLo before */
3404 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3405 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3406 nrange[0] = nrange[1];
3408 root->cell_f[i] = root->bound_max[i];
3410 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3412 nrange[1] = range[1];
3415 if (nrange[1] < range[1]) /* found last a LimLo */
3417 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3418 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3419 nrange[0] = nrange[1];
3420 nrange[1] = range[1];
3421 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3423 else if (nrange[0] > range[0]) /* found at least one LimHi */
3425 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3432 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3433 int d, int dim, gmx_domdec_root_t *root,
3434 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3435 gmx_bool bUniform, gmx_int64_t step)
3437 gmx_domdec_comm_t *comm;
3438 int ncd, d1, i, pos;
3440 real load_aver, load_i, imbalance, change, change_max, sc;
3441 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3445 int range[] = { 0, 0 };
3449 /* Convert the maximum change from the input percentage to a fraction */
3450 change_limit = comm->dlb_scale_lim*0.01;
3454 bPBC = (dim < ddbox->npbcdim);
3456 cell_size = root->buf_ncd;
3458 /* Store the original boundaries */
3459 for (i = 0; i < ncd+1; i++)
3461 root->old_cell_f[i] = root->cell_f[i];
3465 for (i = 0; i < ncd; i++)
3467 cell_size[i] = 1.0/ncd;
3470 else if (dd_load_count(comm) > 0)
3472 load_aver = comm->load[d].sum_m/ncd;
3474 for (i = 0; i < ncd; i++)
3476 /* Determine the relative imbalance of cell i */
3477 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3478 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3479 /* Determine the change of the cell size using underrelaxation */
3480 change = -relax*imbalance;
3481 change_max = std::max(change_max, std::max(change, -change));
3483 /* Limit the amount of scaling.
3484 * We need to use the same rescaling for all cells in one row,
3485 * otherwise the load balancing might not converge.
3488 if (change_max > change_limit)
3490 sc *= change_limit/change_max;
3492 for (i = 0; i < ncd; i++)
3494 /* Determine the relative imbalance of cell i */
3495 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3496 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3497 /* Determine the change of the cell size using underrelaxation */
3498 change = -sc*imbalance;
3499 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3503 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3504 cellsize_limit_f *= DD_CELL_MARGIN;
3505 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3506 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3507 if (ddbox->tric_dir[dim])
3509 cellsize_limit_f /= ddbox->skew_fac[dim];
3510 dist_min_f /= ddbox->skew_fac[dim];
3512 if (bDynamicBox && d > 0)
3514 dist_min_f *= DD_PRES_SCALE_MARGIN;
3516 if (d > 0 && !bUniform)
3518 /* Make sure that the grid is not shifted too much */
3519 for (i = 1; i < ncd; i++)
3521 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3523 gmx_incons("Inconsistent DD boundary staggering limits!");
3525 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3526 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3529 root->bound_min[i] += 0.5*space;
3531 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3532 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3535 root->bound_max[i] += 0.5*space;
3540 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3542 root->cell_f_max0[i-1] + dist_min_f,
3543 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3544 root->cell_f_min1[i] - dist_min_f);
3549 root->cell_f[0] = 0;
3550 root->cell_f[ncd] = 1;
3551 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3554 /* After the checks above, the cells should obey the cut-off
3555 * restrictions, but it does not hurt to check.
3557 for (i = 0; i < ncd; i++)
3561 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3562 dim, i, root->cell_f[i], root->cell_f[i+1]);
3565 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3566 root->cell_f[i+1] - root->cell_f[i] <
3567 cellsize_limit_f/DD_CELL_MARGIN)
3571 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3572 gmx_step_str(step, buf), dim2char(dim), i,
3573 (root->cell_f[i+1] - root->cell_f[i])
3574 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3579 /* Store the cell boundaries of the lower dimensions at the end */
3580 for (d1 = 0; d1 < d; d1++)
3582 root->cell_f[pos++] = comm->cell_f0[d1];
3583 root->cell_f[pos++] = comm->cell_f1[d1];
3586 if (d < comm->npmedecompdim)
3588 /* The master determines the maximum shift for
3589 * the coordinate communication between separate PME nodes.
3591 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3593 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3596 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3600 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3601 gmx_ddbox_t *ddbox, int dimind)
3603 gmx_domdec_comm_t *comm;
3608 /* Set the cell dimensions */
3609 dim = dd->dim[dimind];
3610 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3611 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3612 if (dim >= ddbox->nboundeddim)
3614 comm->cell_x0[dim] += ddbox->box0[dim];
3615 comm->cell_x1[dim] += ddbox->box0[dim];
3619 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3620 int d, int dim, real *cell_f_row,
3623 gmx_domdec_comm_t *comm;
3629 /* Each node would only need to know two fractions,
3630 * but it is probably cheaper to broadcast the whole array.
3632 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3633 0, comm->mpi_comm_load[d]);
3635 /* Copy the fractions for this dimension from the buffer */
3636 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3637 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3638 /* The whole array was communicated, so set the buffer position */
3639 pos = dd->nc[dim] + 1;
3640 for (d1 = 0; d1 <= d; d1++)
3644 /* Copy the cell fractions of the lower dimensions */
3645 comm->cell_f0[d1] = cell_f_row[pos++];
3646 comm->cell_f1[d1] = cell_f_row[pos++];
3648 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3650 /* Convert the communicated shift from float to int */
3651 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3654 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3658 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3659 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3660 gmx_bool bUniform, gmx_int64_t step)
3662 gmx_domdec_comm_t *comm;
3664 gmx_bool bRowMember, bRowRoot;
3669 for (d = 0; d < dd->ndim; d++)
3674 for (d1 = d; d1 < dd->ndim; d1++)
3676 if (dd->ci[dd->dim[d1]] > 0)
3689 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3690 ddbox, bDynamicBox, bUniform, step);
3691 cell_f_row = comm->root[d]->cell_f;
3695 cell_f_row = comm->cell_f_row;
3697 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3702 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3706 /* This function assumes the box is static and should therefore
3707 * not be called when the box has changed since the last
3708 * call to dd_partition_system.
3710 for (d = 0; d < dd->ndim; d++)
3712 relative_to_absolute_cell_bounds(dd, ddbox, d);
3718 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3719 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3720 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3721 gmx_wallcycle_t wcycle)
3723 gmx_domdec_comm_t *comm;
3730 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3731 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3732 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3734 else if (bDynamicBox)
3736 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3739 /* Set the dimensions for which no DD is used */
3740 for (dim = 0; dim < DIM; dim++)
3742 if (dd->nc[dim] == 1)
3744 comm->cell_x0[dim] = 0;
3745 comm->cell_x1[dim] = ddbox->box_size[dim];
3746 if (dim >= ddbox->nboundeddim)
3748 comm->cell_x0[dim] += ddbox->box0[dim];
3749 comm->cell_x1[dim] += ddbox->box0[dim];
3755 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3758 gmx_domdec_comm_dim_t *cd;
3760 for (d = 0; d < dd->ndim; d++)
3762 cd = &dd->comm->cd[d];
3763 np = npulse[dd->dim[d]];
3764 if (np > cd->np_nalloc)
3768 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3769 dim2char(dd->dim[d]), np);
3771 if (DDMASTER(dd) && cd->np_nalloc > 0)
3773 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3775 srenew(cd->ind, np);
3776 for (i = cd->np_nalloc; i < np; i++)
3778 cd->ind[i].index = NULL;
3779 cd->ind[i].nalloc = 0;
3788 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3789 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3790 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3791 gmx_wallcycle_t wcycle)
3793 gmx_domdec_comm_t *comm;
3799 /* Copy the old cell boundaries for the cg displacement check */
3800 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3801 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3803 if (comm->bDynLoadBal)
3807 check_box_size(dd, ddbox);
3809 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3813 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3814 realloc_comm_ind(dd, npulse);
3819 for (d = 0; d < DIM; d++)
3821 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3822 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3827 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3829 rvec cell_ns_x0, rvec cell_ns_x1,
3832 gmx_domdec_comm_t *comm;
3837 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3839 dim = dd->dim[dim_ind];
3841 /* Without PBC we don't have restrictions on the outer cells */
3842 if (!(dim >= ddbox->npbcdim &&
3843 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3844 comm->bDynLoadBal &&
3845 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3846 comm->cellsize_min[dim])
3849 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",
3850 gmx_step_str(step, buf), dim2char(dim),
3851 comm->cell_x1[dim] - comm->cell_x0[dim],
3852 ddbox->skew_fac[dim],
3853 dd->comm->cellsize_min[dim],
3854 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3858 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3860 /* Communicate the boundaries and update cell_ns_x0/1 */
3861 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3862 if (dd->bGridJump && dd->ndim > 1)
3864 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3869 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3873 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3881 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3882 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3891 static void check_screw_box(matrix box)
3893 /* Mathematical limitation */
3894 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3896 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3899 /* Limitation due to the asymmetry of the eighth shell method */
3900 if (box[ZZ][YY] != 0)
3902 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3906 static void distribute_cg(FILE *fplog, gmx_int64_t step,
3907 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3910 gmx_domdec_master_t *ma;
3911 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3912 int i, icg, j, k, k0, k1, d;
3916 real nrcg, inv_ncg, pos_d;
3922 if (tmp_ind == NULL)
3924 snew(tmp_nalloc, dd->nnodes);
3925 snew(tmp_ind, dd->nnodes);
3926 for (i = 0; i < dd->nnodes; i++)
3928 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3929 snew(tmp_ind[i], tmp_nalloc[i]);
3933 /* Clear the count */
3934 for (i = 0; i < dd->nnodes; i++)
3940 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3942 cgindex = cgs->index;
3944 /* Compute the center of geometry for all charge groups */
3945 for (icg = 0; icg < cgs->nr; icg++)
3948 k1 = cgindex[icg+1];
3952 copy_rvec(pos[k0], cg_cm);
3959 for (k = k0; (k < k1); k++)
3961 rvec_inc(cg_cm, pos[k]);
3963 for (d = 0; (d < DIM); d++)
3965 cg_cm[d] *= inv_ncg;
3968 /* Put the charge group in the box and determine the cell index */
3969 for (d = DIM-1; d >= 0; d--)
3972 if (d < dd->npbcdim)
3974 bScrew = (dd->bScrewPBC && d == XX);
3975 if (tric_dir[d] && dd->nc[d] > 1)
3977 /* Use triclinic coordintates for this dimension */
3978 for (j = d+1; j < DIM; j++)
3980 pos_d += cg_cm[j]*tcm[j][d];
3983 while (pos_d >= box[d][d])
3986 rvec_dec(cg_cm, box[d]);
3989 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3990 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3992 for (k = k0; (k < k1); k++)
3994 rvec_dec(pos[k], box[d]);
3997 pos[k][YY] = box[YY][YY] - pos[k][YY];
3998 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4005 rvec_inc(cg_cm, box[d]);
4008 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4009 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4011 for (k = k0; (k < k1); k++)
4013 rvec_inc(pos[k], box[d]);
4016 pos[k][YY] = box[YY][YY] - pos[k][YY];
4017 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4022 /* This could be done more efficiently */
4024 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4029 i = dd_index(dd->nc, ind);
4030 if (ma->ncg[i] == tmp_nalloc[i])
4032 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4033 srenew(tmp_ind[i], tmp_nalloc[i]);
4035 tmp_ind[i][ma->ncg[i]] = icg;
4037 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4041 for (i = 0; i < dd->nnodes; i++)
4044 for (k = 0; k < ma->ncg[i]; k++)
4046 ma->cg[k1++] = tmp_ind[i][k];
4049 ma->index[dd->nnodes] = k1;
4051 for (i = 0; i < dd->nnodes; i++)
4061 fprintf(fplog, "Charge group distribution at step %s:",
4062 gmx_step_str(step, buf));
4063 for (i = 0; i < dd->nnodes; i++)
4065 fprintf(fplog, " %d", ma->ncg[i]);
4067 fprintf(fplog, "\n");
4071 static void get_cg_distribution(FILE *fplog, gmx_int64_t step, gmx_domdec_t *dd,
4072 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4075 gmx_domdec_master_t *ma = NULL;
4078 int *ibuf, buf2[2] = { 0, 0 };
4079 gmx_bool bMaster = DDMASTER(dd);
4087 check_screw_box(box);
4090 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4092 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4093 for (i = 0; i < dd->nnodes; i++)
4095 ma->ibuf[2*i] = ma->ncg[i];
4096 ma->ibuf[2*i+1] = ma->nat[i];
4104 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4106 dd->ncg_home = buf2[0];
4107 dd->nat_home = buf2[1];
4108 dd->ncg_tot = dd->ncg_home;
4109 dd->nat_tot = dd->nat_home;
4110 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4112 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4113 srenew(dd->index_gl, dd->cg_nalloc);
4114 srenew(dd->cgindex, dd->cg_nalloc+1);
4118 for (i = 0; i < dd->nnodes; i++)
4120 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4121 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4126 bMaster ? ma->ibuf : NULL,
4127 bMaster ? ma->ibuf+dd->nnodes : NULL,
4128 bMaster ? ma->cg : NULL,
4129 dd->ncg_home*sizeof(int), dd->index_gl);
4131 /* Determine the home charge group sizes */
4133 for (i = 0; i < dd->ncg_home; i++)
4135 cg_gl = dd->index_gl[i];
4137 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4142 fprintf(debug, "Home charge groups:\n");
4143 for (i = 0; i < dd->ncg_home; i++)
4145 fprintf(debug, " %d", dd->index_gl[i]);
4148 fprintf(debug, "\n");
4151 fprintf(debug, "\n");
4155 static int compact_and_copy_vec_at(int ncg, int *move,
4158 rvec *src, gmx_domdec_comm_t *comm,
4161 int m, icg, i, i0, i1, nrcg;
4167 for (m = 0; m < DIM*2; m++)
4173 for (icg = 0; icg < ncg; icg++)
4175 i1 = cgindex[icg+1];
4181 /* Compact the home array in place */
4182 for (i = i0; i < i1; i++)
4184 copy_rvec(src[i], src[home_pos++]);
4190 /* Copy to the communication buffer */
4192 pos_vec[m] += 1 + vec*nrcg;
4193 for (i = i0; i < i1; i++)
4195 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4197 pos_vec[m] += (nvec - vec - 1)*nrcg;
4201 home_pos += i1 - i0;
4209 static int compact_and_copy_vec_cg(int ncg, int *move,
4211 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4214 int m, icg, i0, i1, nrcg;
4220 for (m = 0; m < DIM*2; m++)
4226 for (icg = 0; icg < ncg; icg++)
4228 i1 = cgindex[icg+1];
4234 /* Compact the home array in place */
4235 copy_rvec(src[icg], src[home_pos++]);
4241 /* Copy to the communication buffer */
4242 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4243 pos_vec[m] += 1 + nrcg*nvec;
4255 static int compact_ind(int ncg, int *move,
4256 int *index_gl, int *cgindex,
4258 gmx_ga2la_t ga2la, char *bLocalCG,
4261 int cg, nat, a0, a1, a, a_gl;
4266 for (cg = 0; cg < ncg; cg++)
4272 /* Compact the home arrays in place.
4273 * Anything that can be done here avoids access to global arrays.
4275 cgindex[home_pos] = nat;
4276 for (a = a0; a < a1; a++)
4279 gatindex[nat] = a_gl;
4280 /* The cell number stays 0, so we don't need to set it */
4281 ga2la_change_la(ga2la, a_gl, nat);
4284 index_gl[home_pos] = index_gl[cg];
4285 cginfo[home_pos] = cginfo[cg];
4286 /* The charge group remains local, so bLocalCG does not change */
4291 /* Clear the global indices */
4292 for (a = a0; a < a1; a++)
4294 ga2la_del(ga2la, gatindex[a]);
4298 bLocalCG[index_gl[cg]] = FALSE;
4302 cgindex[home_pos] = nat;
4307 static void clear_and_mark_ind(int ncg, int *move,
4308 int *index_gl, int *cgindex, int *gatindex,
4309 gmx_ga2la_t ga2la, char *bLocalCG,
4314 for (cg = 0; cg < ncg; cg++)
4320 /* Clear the global indices */
4321 for (a = a0; a < a1; a++)
4323 ga2la_del(ga2la, gatindex[a]);
4327 bLocalCG[index_gl[cg]] = FALSE;
4329 /* Signal that this cg has moved using the ns cell index.
4330 * Here we set it to -1. fill_grid will change it
4331 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4333 cell_index[cg] = -1;
4338 static void print_cg_move(FILE *fplog,
4340 gmx_int64_t step, int cg, int dim, int dir,
4341 gmx_bool bHaveCgcmOld, real limitd,
4342 rvec cm_old, rvec cm_new, real pos_d)
4344 gmx_domdec_comm_t *comm;
4349 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4352 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4353 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4354 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4358 /* We don't have a limiting distance available: don't print it */
4359 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4360 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4361 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4363 fprintf(fplog, "distance out of cell %f\n",
4364 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4367 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4368 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4370 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4371 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4372 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4374 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4375 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4377 comm->cell_x0[dim], comm->cell_x1[dim]);
4380 static void cg_move_error(FILE *fplog,
4382 gmx_int64_t step, int cg, int dim, int dir,
4383 gmx_bool bHaveCgcmOld, real limitd,
4384 rvec cm_old, rvec cm_new, real pos_d)
4388 print_cg_move(fplog, dd, step, cg, dim, dir,
4389 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4391 print_cg_move(stderr, dd, step, cg, dim, dir,
4392 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4394 "%s moved too far between two domain decomposition steps\n"
4395 "This usually means that your system is not well equilibrated",
4396 dd->comm->bCGs ? "A charge group" : "An atom");
4399 static void rotate_state_atom(t_state *state, int a)
4403 for (est = 0; est < estNR; est++)
4405 if (EST_DISTR(est) && (state->flags & (1<<est)))
4410 /* Rotate the complete state; for a rectangular box only */
4411 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4412 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4415 state->v[a][YY] = -state->v[a][YY];
4416 state->v[a][ZZ] = -state->v[a][ZZ];
4419 state->sd_X[a][YY] = -state->sd_X[a][YY];
4420 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4423 state->cg_p[a][YY] = -state->cg_p[a][YY];
4424 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4426 case estDISRE_INITF:
4427 case estDISRE_RM3TAV:
4428 case estORIRE_INITF:
4430 /* These are distances, so not affected by rotation */
4433 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4439 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4441 if (natoms > comm->moved_nalloc)
4443 /* Contents should be preserved here */
4444 comm->moved_nalloc = over_alloc_dd(natoms);
4445 srenew(comm->moved, comm->moved_nalloc);
4451 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4454 ivec tric_dir, matrix tcm,
4455 rvec cell_x0, rvec cell_x1,
4456 rvec limitd, rvec limit0, rvec limit1,
4458 int cg_start, int cg_end,
4463 int cg, k, k0, k1, d, dim, d2;
4468 real inv_ncg, pos_d;
4471 npbcdim = dd->npbcdim;
4473 for (cg = cg_start; cg < cg_end; cg++)
4480 copy_rvec(state->x[k0], cm_new);
4487 for (k = k0; (k < k1); k++)
4489 rvec_inc(cm_new, state->x[k]);
4491 for (d = 0; (d < DIM); d++)
4493 cm_new[d] = inv_ncg*cm_new[d];
4498 /* Do pbc and check DD cell boundary crossings */
4499 for (d = DIM-1; d >= 0; d--)
4503 bScrew = (dd->bScrewPBC && d == XX);
4504 /* Determine the location of this cg in lattice coordinates */
4508 for (d2 = d+1; d2 < DIM; d2++)
4510 pos_d += cm_new[d2]*tcm[d2][d];
4513 /* Put the charge group in the triclinic unit-cell */
4514 if (pos_d >= cell_x1[d])
4516 if (pos_d >= limit1[d])
4518 cg_move_error(fplog, dd, step, cg, d, 1,
4519 cg_cm != state->x, limitd[d],
4520 cg_cm[cg], cm_new, pos_d);
4523 if (dd->ci[d] == dd->nc[d] - 1)
4525 rvec_dec(cm_new, state->box[d]);
4528 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4529 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4531 for (k = k0; (k < k1); k++)
4533 rvec_dec(state->x[k], state->box[d]);
4536 rotate_state_atom(state, k);
4541 else if (pos_d < cell_x0[d])
4543 if (pos_d < limit0[d])
4545 cg_move_error(fplog, dd, step, cg, d, -1,
4546 cg_cm != state->x, limitd[d],
4547 cg_cm[cg], cm_new, pos_d);
4552 rvec_inc(cm_new, state->box[d]);
4555 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4556 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4558 for (k = k0; (k < k1); k++)
4560 rvec_inc(state->x[k], state->box[d]);
4563 rotate_state_atom(state, k);
4569 else if (d < npbcdim)
4571 /* Put the charge group in the rectangular unit-cell */
4572 while (cm_new[d] >= state->box[d][d])
4574 rvec_dec(cm_new, state->box[d]);
4575 for (k = k0; (k < k1); k++)
4577 rvec_dec(state->x[k], state->box[d]);
4580 while (cm_new[d] < 0)
4582 rvec_inc(cm_new, state->box[d]);
4583 for (k = k0; (k < k1); k++)
4585 rvec_inc(state->x[k], state->box[d]);
4591 copy_rvec(cm_new, cg_cm[cg]);
4593 /* Determine where this cg should go */
4596 for (d = 0; d < dd->ndim; d++)
4601 flag |= DD_FLAG_FW(d);
4607 else if (dev[dim] == -1)
4609 flag |= DD_FLAG_BW(d);
4612 if (dd->nc[dim] > 2)
4623 /* Temporarily store the flag in move */
4624 move[cg] = mc + flag;
4628 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4629 gmx_domdec_t *dd, ivec tric_dir,
4630 t_state *state, rvec **f,
4639 int ncg[DIM*2], nat[DIM*2];
4640 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4641 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4642 int sbuf[2], rbuf[2];
4643 int home_pos_cg, home_pos_at, buf_pos;
4645 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4648 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4650 cginfo_mb_t *cginfo_mb;
4651 gmx_domdec_comm_t *comm;
4653 int nthread, thread;
4657 check_screw_box(state->box);
4661 if (fr->cutoff_scheme == ecutsGROUP)
4666 for (i = 0; i < estNR; i++)
4672 case estX: /* Always present */ break;
4673 case estV: bV = (state->flags & (1<<i)); break;
4674 case estSDX: bSDX = (state->flags & (1<<i)); break;
4675 case estCGP: bCGP = (state->flags & (1<<i)); break;
4678 case estDISRE_INITF:
4679 case estDISRE_RM3TAV:
4680 case estORIRE_INITF:
4682 /* No processing required */
4685 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4690 if (dd->ncg_tot > comm->nalloc_int)
4692 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4693 srenew(comm->buf_int, comm->nalloc_int);
4695 move = comm->buf_int;
4697 /* Clear the count */
4698 for (c = 0; c < dd->ndim*2; c++)
4704 npbcdim = dd->npbcdim;
4706 for (d = 0; (d < DIM); d++)
4708 limitd[d] = dd->comm->cellsize_min[d];
4709 if (d >= npbcdim && dd->ci[d] == 0)
4711 cell_x0[d] = -GMX_FLOAT_MAX;
4715 cell_x0[d] = comm->cell_x0[d];
4717 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4719 cell_x1[d] = GMX_FLOAT_MAX;
4723 cell_x1[d] = comm->cell_x1[d];
4727 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4728 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4732 /* We check after communication if a charge group moved
4733 * more than one cell. Set the pre-comm check limit to float_max.
4735 limit0[d] = -GMX_FLOAT_MAX;
4736 limit1[d] = GMX_FLOAT_MAX;
4740 make_tric_corr_matrix(npbcdim, state->box, tcm);
4742 cgindex = dd->cgindex;
4744 nthread = gmx_omp_nthreads_get(emntDomdec);
4746 /* Compute the center of geometry for all home charge groups
4747 * and put them in the box and determine where they should go.
4749 #pragma omp parallel for num_threads(nthread) schedule(static)
4750 for (thread = 0; thread < nthread; thread++)
4752 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4753 cell_x0, cell_x1, limitd, limit0, limit1,
4755 ( thread *dd->ncg_home)/nthread,
4756 ((thread+1)*dd->ncg_home)/nthread,
4757 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4761 for (cg = 0; cg < dd->ncg_home; cg++)
4766 flag = mc & ~DD_FLAG_NRCG;
4767 mc = mc & DD_FLAG_NRCG;
4770 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4772 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4773 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4775 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4776 /* We store the cg size in the lower 16 bits
4777 * and the place where the charge group should go
4778 * in the next 6 bits. This saves some communication volume.
4780 nrcg = cgindex[cg+1] - cgindex[cg];
4781 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4787 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4788 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4791 for (i = 0; i < dd->ndim*2; i++)
4793 *ncg_moved += ncg[i];
4810 /* Make sure the communication buffers are large enough */
4811 for (mc = 0; mc < dd->ndim*2; mc++)
4813 nvr = ncg[mc] + nat[mc]*nvec;
4814 if (nvr > comm->cgcm_state_nalloc[mc])
4816 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4817 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4821 switch (fr->cutoff_scheme)
4824 /* Recalculating cg_cm might be cheaper than communicating,
4825 * but that could give rise to rounding issues.
4828 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4829 nvec, cg_cm, comm, bCompact);
4832 /* Without charge groups we send the moved atom coordinates
4833 * over twice. This is so the code below can be used without
4834 * many conditionals for both for with and without charge groups.
4837 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4838 nvec, state->x, comm, FALSE);
4841 home_pos_cg -= *ncg_moved;
4845 gmx_incons("unimplemented");
4851 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4852 nvec, vec++, state->x, comm, bCompact);
4855 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4856 nvec, vec++, state->v, comm, bCompact);
4860 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4861 nvec, vec++, state->sd_X, comm, bCompact);
4865 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4866 nvec, vec++, state->cg_p, comm, bCompact);
4871 compact_ind(dd->ncg_home, move,
4872 dd->index_gl, dd->cgindex, dd->gatindex,
4873 dd->ga2la, comm->bLocalCG,
4878 if (fr->cutoff_scheme == ecutsVERLET)
4880 moved = get_moved(comm, dd->ncg_home);
4882 for (k = 0; k < dd->ncg_home; k++)
4889 moved = fr->ns.grid->cell_index;
4892 clear_and_mark_ind(dd->ncg_home, move,
4893 dd->index_gl, dd->cgindex, dd->gatindex,
4894 dd->ga2la, comm->bLocalCG,
4898 cginfo_mb = fr->cginfo_mb;
4900 *ncg_stay_home = home_pos_cg;
4901 for (d = 0; d < dd->ndim; d++)
4906 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4909 /* Communicate the cg and atom counts */
4914 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4915 d, dir, sbuf[0], sbuf[1]);
4917 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4919 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4921 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4922 srenew(comm->buf_int, comm->nalloc_int);
4925 /* Communicate the charge group indices, sizes and flags */
4926 dd_sendrecv_int(dd, d, dir,
4927 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4928 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4930 nvs = ncg[cdd] + nat[cdd]*nvec;
4931 i = rbuf[0] + rbuf[1] *nvec;
4932 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4934 /* Communicate cgcm and state */
4935 dd_sendrecv_rvec(dd, d, dir,
4936 comm->cgcm_state[cdd], nvs,
4937 comm->vbuf.v+nvr, i);
4938 ncg_recv += rbuf[0];
4942 /* Process the received charge groups */
4944 for (cg = 0; cg < ncg_recv; cg++)
4946 flag = comm->buf_int[cg*DD_CGIBS+1];
4948 if (dim >= npbcdim && dd->nc[dim] > 2)
4950 /* No pbc in this dim and more than one domain boundary.
4951 * We do a separate check if a charge group didn't move too far.
4953 if (((flag & DD_FLAG_FW(d)) &&
4954 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4955 ((flag & DD_FLAG_BW(d)) &&
4956 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4958 cg_move_error(fplog, dd, step, cg, dim,
4959 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4960 fr->cutoff_scheme == ecutsGROUP, 0,
4961 comm->vbuf.v[buf_pos],
4962 comm->vbuf.v[buf_pos],
4963 comm->vbuf.v[buf_pos][dim]);
4970 /* Check which direction this cg should go */
4971 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4975 /* The cell boundaries for dimension d2 are not equal
4976 * for each cell row of the lower dimension(s),
4977 * therefore we might need to redetermine where
4978 * this cg should go.
4981 /* If this cg crosses the box boundary in dimension d2
4982 * we can use the communicated flag, so we do not
4983 * have to worry about pbc.
4985 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4986 (flag & DD_FLAG_FW(d2))) ||
4987 (dd->ci[dim2] == 0 &&
4988 (flag & DD_FLAG_BW(d2)))))
4990 /* Clear the two flags for this dimension */
4991 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4992 /* Determine the location of this cg
4993 * in lattice coordinates
4995 pos_d = comm->vbuf.v[buf_pos][dim2];
4998 for (d3 = dim2+1; d3 < DIM; d3++)
5001 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5004 /* Check of we are not at the box edge.
5005 * pbc is only handled in the first step above,
5006 * but this check could move over pbc while
5007 * the first step did not due to different rounding.
5009 if (pos_d >= cell_x1[dim2] &&
5010 dd->ci[dim2] != dd->nc[dim2]-1)
5012 flag |= DD_FLAG_FW(d2);
5014 else if (pos_d < cell_x0[dim2] &&
5017 flag |= DD_FLAG_BW(d2);
5019 comm->buf_int[cg*DD_CGIBS+1] = flag;
5022 /* Set to which neighboring cell this cg should go */
5023 if (flag & DD_FLAG_FW(d2))
5027 else if (flag & DD_FLAG_BW(d2))
5029 if (dd->nc[dd->dim[d2]] > 2)
5041 nrcg = flag & DD_FLAG_NRCG;
5044 if (home_pos_cg+1 > dd->cg_nalloc)
5046 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5047 srenew(dd->index_gl, dd->cg_nalloc);
5048 srenew(dd->cgindex, dd->cg_nalloc+1);
5050 /* Set the global charge group index and size */
5051 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5052 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5053 /* Copy the state from the buffer */
5054 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5055 if (fr->cutoff_scheme == ecutsGROUP)
5058 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5062 /* Set the cginfo */
5063 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5064 dd->index_gl[home_pos_cg]);
5067 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5070 if (home_pos_at+nrcg > state->nalloc)
5072 dd_realloc_state(state, f, home_pos_at+nrcg);
5074 for (i = 0; i < nrcg; i++)
5076 copy_rvec(comm->vbuf.v[buf_pos++],
5077 state->x[home_pos_at+i]);
5081 for (i = 0; i < nrcg; i++)
5083 copy_rvec(comm->vbuf.v[buf_pos++],
5084 state->v[home_pos_at+i]);
5089 for (i = 0; i < nrcg; i++)
5091 copy_rvec(comm->vbuf.v[buf_pos++],
5092 state->sd_X[home_pos_at+i]);
5097 for (i = 0; i < nrcg; i++)
5099 copy_rvec(comm->vbuf.v[buf_pos++],
5100 state->cg_p[home_pos_at+i]);
5104 home_pos_at += nrcg;
5108 /* Reallocate the buffers if necessary */
5109 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5111 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5112 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5114 nvr = ncg[mc] + nat[mc]*nvec;
5115 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5117 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5118 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5120 /* Copy from the receive to the send buffers */
5121 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5122 comm->buf_int + cg*DD_CGIBS,
5123 DD_CGIBS*sizeof(int));
5124 memcpy(comm->cgcm_state[mc][nvr],
5125 comm->vbuf.v[buf_pos],
5126 (1+nrcg*nvec)*sizeof(rvec));
5127 buf_pos += 1 + nrcg*nvec;
5134 /* With sorting (!bCompact) the indices are now only partially up to date
5135 * and ncg_home and nat_home are not the real count, since there are
5136 * "holes" in the arrays for the charge groups that moved to neighbors.
5138 if (fr->cutoff_scheme == ecutsVERLET)
5140 moved = get_moved(comm, home_pos_cg);
5142 for (i = dd->ncg_home; i < home_pos_cg; i++)
5147 dd->ncg_home = home_pos_cg;
5148 dd->nat_home = home_pos_at;
5153 "Finished repartitioning: cgs moved out %d, new home %d\n",
5154 *ncg_moved, dd->ncg_home-*ncg_moved);
5159 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5161 dd->comm->cycl[ddCycl] += cycles;
5162 dd->comm->cycl_n[ddCycl]++;
5163 if (cycles > dd->comm->cycl_max[ddCycl])
5165 dd->comm->cycl_max[ddCycl] = cycles;
5169 static double force_flop_count(t_nrnb *nrnb)
5176 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5178 /* To get closer to the real timings, we half the count
5179 * for the normal loops and again half it for water loops.
5182 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5184 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5188 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5191 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5194 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5196 sum += nrnb->n[i]*cost_nrnb(i);
5199 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5201 sum += nrnb->n[i]*cost_nrnb(i);
5207 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5209 if (dd->comm->eFlop)
5211 dd->comm->flop -= force_flop_count(nrnb);
5214 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5216 if (dd->comm->eFlop)
5218 dd->comm->flop += force_flop_count(nrnb);
5223 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5227 for (i = 0; i < ddCyclNr; i++)
5229 dd->comm->cycl[i] = 0;
5230 dd->comm->cycl_n[i] = 0;
5231 dd->comm->cycl_max[i] = 0;
5234 dd->comm->flop_n = 0;
5237 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5239 gmx_domdec_comm_t *comm;
5240 gmx_domdec_load_t *load;
5241 gmx_domdec_root_t *root = NULL;
5243 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5248 fprintf(debug, "get_load_distribution start\n");
5251 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5255 bSepPME = (dd->pme_nodeid >= 0);
5257 for (d = dd->ndim-1; d >= 0; d--)
5260 /* Check if we participate in the communication in this dimension */
5261 if (d == dd->ndim-1 ||
5262 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5264 load = &comm->load[d];
5267 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5270 if (d == dd->ndim-1)
5272 sbuf[pos++] = dd_force_load(comm);
5273 sbuf[pos++] = sbuf[0];
5276 sbuf[pos++] = sbuf[0];
5277 sbuf[pos++] = cell_frac;
5280 sbuf[pos++] = comm->cell_f_max0[d];
5281 sbuf[pos++] = comm->cell_f_min1[d];
5286 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5287 sbuf[pos++] = comm->cycl[ddCyclPME];
5292 sbuf[pos++] = comm->load[d+1].sum;
5293 sbuf[pos++] = comm->load[d+1].max;
5296 sbuf[pos++] = comm->load[d+1].sum_m;
5297 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5298 sbuf[pos++] = comm->load[d+1].flags;
5301 sbuf[pos++] = comm->cell_f_max0[d];
5302 sbuf[pos++] = comm->cell_f_min1[d];
5307 sbuf[pos++] = comm->load[d+1].mdf;
5308 sbuf[pos++] = comm->load[d+1].pme;
5312 /* Communicate a row in DD direction d.
5313 * The communicators are setup such that the root always has rank 0.
5316 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5317 load->load, load->nload*sizeof(float), MPI_BYTE,
5318 0, comm->mpi_comm_load[d]);
5320 if (dd->ci[dim] == dd->master_ci[dim])
5322 /* We are the root, process this row */
5323 if (comm->bDynLoadBal)
5325 root = comm->root[d];
5335 for (i = 0; i < dd->nc[dim]; i++)
5337 load->sum += load->load[pos++];
5338 load->max = std::max(load->max, load->load[pos]);
5344 /* This direction could not be load balanced properly,
5345 * therefore we need to use the maximum iso the average load.
5347 load->sum_m = std::max(load->sum_m, load->load[pos]);
5351 load->sum_m += load->load[pos];
5354 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5358 load->flags = (int)(load->load[pos++] + 0.5);
5362 root->cell_f_max0[i] = load->load[pos++];
5363 root->cell_f_min1[i] = load->load[pos++];
5368 load->mdf = std::max(load->mdf, load->load[pos]);
5370 load->pme = std::max(load->pme, load->load[pos]);
5374 if (comm->bDynLoadBal && root->bLimited)
5376 load->sum_m *= dd->nc[dim];
5377 load->flags |= (1<<d);
5385 comm->nload += dd_load_count(comm);
5386 comm->load_step += comm->cycl[ddCyclStep];
5387 comm->load_sum += comm->load[0].sum;
5388 comm->load_max += comm->load[0].max;
5389 if (comm->bDynLoadBal)
5391 for (d = 0; d < dd->ndim; d++)
5393 if (comm->load[0].flags & (1<<d))
5395 comm->load_lim[d]++;
5401 comm->load_mdf += comm->load[0].mdf;
5402 comm->load_pme += comm->load[0].pme;
5406 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5410 fprintf(debug, "get_load_distribution finished\n");
5414 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5416 /* Return the relative performance loss on the total run time
5417 * due to the force calculation load imbalance.
5419 if (dd->comm->nload > 0)
5422 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5423 (dd->comm->load_step*dd->nnodes);
5431 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5434 int npp, npme, nnodes, d, limp;
5435 float imbal, pme_f_ratio, lossf, lossp = 0;
5437 gmx_domdec_comm_t *comm;
5440 if (DDMASTER(dd) && comm->nload > 0)
5443 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5444 nnodes = npp + npme;
5445 imbal = comm->load_max*npp/comm->load_sum - 1;
5446 lossf = dd_force_imb_perf_loss(dd);
5447 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5448 fprintf(fplog, "%s", buf);
5449 fprintf(stderr, "\n");
5450 fprintf(stderr, "%s", buf);
5451 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5452 fprintf(fplog, "%s", buf);
5453 fprintf(stderr, "%s", buf);
5455 if (comm->bDynLoadBal)
5457 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5458 for (d = 0; d < dd->ndim; d++)
5460 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5461 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5467 sprintf(buf+strlen(buf), "\n");
5468 fprintf(fplog, "%s", buf);
5469 fprintf(stderr, "%s", buf);
5473 pme_f_ratio = comm->load_pme/comm->load_mdf;
5474 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5477 lossp *= (float)npme/(float)nnodes;
5481 lossp *= (float)npp/(float)nnodes;
5483 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5484 fprintf(fplog, "%s", buf);
5485 fprintf(stderr, "%s", buf);
5486 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5487 fprintf(fplog, "%s", buf);
5488 fprintf(stderr, "%s", buf);
5490 fprintf(fplog, "\n");
5491 fprintf(stderr, "\n");
5493 if (lossf >= DD_PERF_LOSS_WARN)
5496 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5497 " in the domain decomposition.\n", lossf*100);
5498 if (!comm->bDynLoadBal)
5500 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5504 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5506 fprintf(fplog, "%s\n", buf);
5507 fprintf(stderr, "%s\n", buf);
5509 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5512 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5513 " had %s work to do than the PP ranks.\n"
5514 " You might want to %s the number of PME ranks\n"
5515 " or %s the cut-off and the grid spacing.\n",
5517 (lossp < 0) ? "less" : "more",
5518 (lossp < 0) ? "decrease" : "increase",
5519 (lossp < 0) ? "decrease" : "increase");
5520 fprintf(fplog, "%s\n", buf);
5521 fprintf(stderr, "%s\n", buf);
5526 static float dd_vol_min(gmx_domdec_t *dd)
5528 return dd->comm->load[0].cvol_min*dd->nnodes;
5531 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5533 return dd->comm->load[0].flags;
5536 static float dd_f_imbal(gmx_domdec_t *dd)
5538 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5541 float dd_pme_f_ratio(gmx_domdec_t *dd)
5543 if (dd->comm->cycl_n[ddCyclPME] > 0)
5545 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5553 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5558 flags = dd_load_flags(dd);
5562 "DD load balancing is limited by minimum cell size in dimension");
5563 for (d = 0; d < dd->ndim; d++)
5567 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5570 fprintf(fplog, "\n");
5572 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5573 if (dd->comm->bDynLoadBal)
5575 fprintf(fplog, " vol min/aver %5.3f%c",
5576 dd_vol_min(dd), flags ? '!' : ' ');
5578 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5579 if (dd->comm->cycl_n[ddCyclPME])
5581 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5583 fprintf(fplog, "\n\n");
5586 static void dd_print_load_verbose(gmx_domdec_t *dd)
5588 if (dd->comm->bDynLoadBal)
5590 fprintf(stderr, "vol %4.2f%c ",
5591 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5593 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5594 if (dd->comm->cycl_n[ddCyclPME])
5596 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5601 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5606 gmx_domdec_root_t *root;
5607 gmx_bool bPartOfGroup = FALSE;
5609 dim = dd->dim[dim_ind];
5610 copy_ivec(loc, loc_c);
5611 for (i = 0; i < dd->nc[dim]; i++)
5614 rank = dd_index(dd->nc, loc_c);
5615 if (rank == dd->rank)
5617 /* This process is part of the group */
5618 bPartOfGroup = TRUE;
5621 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5625 dd->comm->mpi_comm_load[dim_ind] = c_row;
5626 if (dd->comm->eDLB != edlbNO)
5628 if (dd->ci[dim] == dd->master_ci[dim])
5630 /* This is the root process of this row */
5631 snew(dd->comm->root[dim_ind], 1);
5632 root = dd->comm->root[dim_ind];
5633 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5634 snew(root->old_cell_f, dd->nc[dim]+1);
5635 snew(root->bCellMin, dd->nc[dim]);
5638 snew(root->cell_f_max0, dd->nc[dim]);
5639 snew(root->cell_f_min1, dd->nc[dim]);
5640 snew(root->bound_min, dd->nc[dim]);
5641 snew(root->bound_max, dd->nc[dim]);
5643 snew(root->buf_ncd, dd->nc[dim]);
5647 /* This is not a root process, we only need to receive cell_f */
5648 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5651 if (dd->ci[dim] == dd->master_ci[dim])
5653 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5659 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5660 const gmx_hw_info_t gmx_unused *hwinfo,
5661 const gmx_hw_opt_t gmx_unused *hw_opt)
5664 int physicalnode_id_hash;
5667 MPI_Comm mpi_comm_pp_physicalnode;
5669 if (!(cr->duty & DUTY_PP) ||
5670 hw_opt->gpu_opt.ncuda_dev_use == 0)
5672 /* Only PP nodes (currently) use GPUs.
5673 * If we don't have GPUs, there are no resources to share.
5678 physicalnode_id_hash = gmx_physicalnode_id_hash();
5680 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5686 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5687 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5688 dd->rank, physicalnode_id_hash, gpu_id);
5690 /* Split the PP communicator over the physical nodes */
5691 /* TODO: See if we should store this (before), as it's also used for
5692 * for the nodecomm summution.
5694 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5695 &mpi_comm_pp_physicalnode);
5696 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5697 &dd->comm->mpi_comm_gpu_shared);
5698 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5699 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5703 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5706 /* Note that some ranks could share a GPU, while others don't */
5708 if (dd->comm->nrank_gpu_shared == 1)
5710 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5715 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5718 int dim0, dim1, i, j;
5723 fprintf(debug, "Making load communicators\n");
5726 snew(dd->comm->load, dd->ndim);
5727 snew(dd->comm->mpi_comm_load, dd->ndim);
5730 make_load_communicator(dd, 0, loc);
5734 for (i = 0; i < dd->nc[dim0]; i++)
5737 make_load_communicator(dd, 1, loc);
5743 for (i = 0; i < dd->nc[dim0]; i++)
5747 for (j = 0; j < dd->nc[dim1]; j++)
5750 make_load_communicator(dd, 2, loc);
5757 fprintf(debug, "Finished making load communicators\n");
5762 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5764 int d, dim, i, j, m;
5767 ivec dd_zp[DD_MAXIZONE];
5768 gmx_domdec_zones_t *zones;
5769 gmx_domdec_ns_ranges_t *izone;
5771 for (d = 0; d < dd->ndim; d++)
5774 copy_ivec(dd->ci, tmp);
5775 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5776 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5777 copy_ivec(dd->ci, tmp);
5778 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5779 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5782 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5785 dd->neighbor[d][1]);
5791 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5793 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5794 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5801 for (i = 0; i < nzonep; i++)
5803 copy_ivec(dd_zp3[i], dd_zp[i]);
5809 for (i = 0; i < nzonep; i++)
5811 copy_ivec(dd_zp2[i], dd_zp[i]);
5817 for (i = 0; i < nzonep; i++)
5819 copy_ivec(dd_zp1[i], dd_zp[i]);
5823 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5828 zones = &dd->comm->zones;
5830 for (i = 0; i < nzone; i++)
5833 clear_ivec(zones->shift[i]);
5834 for (d = 0; d < dd->ndim; d++)
5836 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5841 for (i = 0; i < nzone; i++)
5843 for (d = 0; d < DIM; d++)
5845 s[d] = dd->ci[d] - zones->shift[i][d];
5850 else if (s[d] >= dd->nc[d])
5856 zones->nizone = nzonep;
5857 for (i = 0; i < zones->nizone; i++)
5859 if (dd_zp[i][0] != i)
5861 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5863 izone = &zones->izone[i];
5864 izone->j0 = dd_zp[i][1];
5865 izone->j1 = dd_zp[i][2];
5866 for (dim = 0; dim < DIM; dim++)
5868 if (dd->nc[dim] == 1)
5870 /* All shifts should be allowed */
5871 izone->shift0[dim] = -1;
5872 izone->shift1[dim] = 1;
5877 izone->shift0[d] = 0;
5878 izone->shift1[d] = 0;
5879 for(j=izone->j0; j<izone->j1; j++) {
5880 if (dd->shift[j][d] > dd->shift[i][d])
5881 izone->shift0[d] = -1;
5882 if (dd->shift[j][d] < dd->shift[i][d])
5883 izone->shift1[d] = 1;
5889 /* Assume the shift are not more than 1 cell */
5890 izone->shift0[dim] = 1;
5891 izone->shift1[dim] = -1;
5892 for (j = izone->j0; j < izone->j1; j++)
5894 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5895 if (shift_diff < izone->shift0[dim])
5897 izone->shift0[dim] = shift_diff;
5899 if (shift_diff > izone->shift1[dim])
5901 izone->shift1[dim] = shift_diff;
5908 if (dd->comm->eDLB != edlbNO)
5910 snew(dd->comm->root, dd->ndim);
5913 if (dd->comm->bRecordLoad)
5915 make_load_communicators(dd);
5919 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5925 gmx_domdec_comm_t *comm;
5932 if (comm->bCartesianPP)
5934 /* Set up cartesian communication for the particle-particle part */
5937 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5938 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5941 for (int i = 0; i < DIM; i++)
5945 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5947 /* We overwrite the old communicator with the new cartesian one */
5948 cr->mpi_comm_mygroup = comm_cart;
5951 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5952 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5954 if (comm->bCartesianPP_PME)
5956 /* Since we want to use the original cartesian setup for sim,
5957 * and not the one after split, we need to make an index.
5959 snew(comm->ddindex2ddnodeid, dd->nnodes);
5960 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5961 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5962 /* Get the rank of the DD master,
5963 * above we made sure that the master node is a PP node.
5973 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5975 else if (comm->bCartesianPP)
5977 if (cr->npmenodes == 0)
5979 /* The PP communicator is also
5980 * the communicator for this simulation
5982 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5984 cr->nodeid = dd->rank;
5986 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5988 /* We need to make an index to go from the coordinates
5989 * to the nodeid of this simulation.
5991 snew(comm->ddindex2simnodeid, dd->nnodes);
5992 snew(buf, dd->nnodes);
5993 if (cr->duty & DUTY_PP)
5995 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5997 /* Communicate the ddindex to simulation nodeid index */
5998 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5999 cr->mpi_comm_mysim);
6002 /* Determine the master coordinates and rank.
6003 * The DD master should be the same node as the master of this sim.
6005 for (int i = 0; i < dd->nnodes; i++)
6007 if (comm->ddindex2simnodeid[i] == 0)
6009 ddindex2xyz(dd->nc, i, dd->master_ci);
6010 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6015 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6020 /* No Cartesian communicators */
6021 /* We use the rank in dd->comm->all as DD index */
6022 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6023 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6025 clear_ivec(dd->master_ci);
6032 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6033 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6038 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6039 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6043 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6047 gmx_domdec_comm_t *comm;
6052 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6055 snew(comm->ddindex2simnodeid, dd->nnodes);
6056 snew(buf, dd->nnodes);
6057 if (cr->duty & DUTY_PP)
6059 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6061 /* Communicate the ddindex to simulation nodeid index */
6062 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6063 cr->mpi_comm_mysim);
6069 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6070 int ncg, int natoms)
6072 gmx_domdec_master_t *ma;
6077 snew(ma->ncg, dd->nnodes);
6078 snew(ma->index, dd->nnodes+1);
6080 snew(ma->nat, dd->nnodes);
6081 snew(ma->ibuf, dd->nnodes*2);
6082 snew(ma->cell_x, DIM);
6083 for (i = 0; i < DIM; i++)
6085 snew(ma->cell_x[i], dd->nc[i]+1);
6088 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6094 snew(ma->vbuf, natoms);
6100 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6101 int gmx_unused reorder)
6104 gmx_domdec_comm_t *comm;
6114 if (comm->bCartesianPP)
6116 for (i = 1; i < DIM; i++)
6118 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6120 if (bDiv[YY] || bDiv[ZZ])
6122 comm->bCartesianPP_PME = TRUE;
6123 /* If we have 2D PME decomposition, which is always in x+y,
6124 * we stack the PME only nodes in z.
6125 * Otherwise we choose the direction that provides the thinnest slab
6126 * of PME only nodes as this will have the least effect
6127 * on the PP communication.
6128 * But for the PME communication the opposite might be better.
6130 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6132 dd->nc[YY] > dd->nc[ZZ]))
6134 comm->cartpmedim = ZZ;
6138 comm->cartpmedim = YY;
6140 comm->ntot[comm->cartpmedim]
6141 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6145 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]);
6147 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6152 if (comm->bCartesianPP_PME)
6159 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]);
6162 for (i = 0; i < DIM; i++)
6166 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6168 MPI_Comm_rank(comm_cart, &rank);
6169 if (MASTERNODE(cr) && rank != 0)
6171 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6174 /* With this assigment we loose the link to the original communicator
6175 * which will usually be MPI_COMM_WORLD, unless have multisim.
6177 cr->mpi_comm_mysim = comm_cart;
6178 cr->sim_nodeid = rank;
6180 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6184 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6185 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6188 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6192 if (cr->npmenodes == 0 ||
6193 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6195 cr->duty = DUTY_PME;
6198 /* Split the sim communicator into PP and PME only nodes */
6199 MPI_Comm_split(cr->mpi_comm_mysim,
6201 dd_index(comm->ntot, dd->ci),
6202 &cr->mpi_comm_mygroup);
6206 switch (dd_node_order)
6211 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6214 case ddnoINTERLEAVE:
6215 /* Interleave the PP-only and PME-only nodes,
6216 * as on clusters with dual-core machines this will double
6217 * the communication bandwidth of the PME processes
6218 * and thus speed up the PP <-> PME and inter PME communication.
6222 fprintf(fplog, "Interleaving PP and PME ranks\n");
6224 comm->pmenodes = dd_pmenodes(cr);
6229 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6232 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6234 cr->duty = DUTY_PME;
6241 /* Split the sim communicator into PP and PME only nodes */
6242 MPI_Comm_split(cr->mpi_comm_mysim,
6245 &cr->mpi_comm_mygroup);
6246 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6252 fprintf(fplog, "This rank does only %s work.\n\n",
6253 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6257 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6260 gmx_domdec_comm_t *comm;
6266 copy_ivec(dd->nc, comm->ntot);
6268 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6269 comm->bCartesianPP_PME = FALSE;
6271 /* Reorder the nodes by default. This might change the MPI ranks.
6272 * Real reordering is only supported on very few architectures,
6273 * Blue Gene is one of them.
6275 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6277 if (cr->npmenodes > 0)
6279 /* Split the communicator into a PP and PME part */
6280 split_communicator(fplog, cr, dd_node_order, CartReorder);
6281 if (comm->bCartesianPP_PME)
6283 /* We (possibly) reordered the nodes in split_communicator,
6284 * so it is no longer required in make_pp_communicator.
6286 CartReorder = FALSE;
6291 /* All nodes do PP and PME */
6293 /* We do not require separate communicators */
6294 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6298 if (cr->duty & DUTY_PP)
6300 /* Copy or make a new PP communicator */
6301 make_pp_communicator(fplog, cr, CartReorder);
6305 receive_ddindex2simnodeid(cr);
6308 if (!(cr->duty & DUTY_PME))
6310 /* Set up the commnuication to our PME node */
6311 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6312 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6315 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6316 dd->pme_nodeid, dd->pme_receive_vir_ener);
6321 dd->pme_nodeid = -1;
6326 dd->ma = init_gmx_domdec_master_t(dd,
6328 comm->cgs_gl.index[comm->cgs_gl.nr]);
6332 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6334 real *slb_frac, tot;
6339 if (nc > 1 && size_string != NULL)
6343 fprintf(fplog, "Using static load balancing for the %s direction\n",
6348 for (i = 0; i < nc; i++)
6351 sscanf(size_string, "%20lf%n", &dbl, &n);
6354 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6363 fprintf(fplog, "Relative cell sizes:");
6365 for (i = 0; i < nc; i++)
6370 fprintf(fplog, " %5.3f", slb_frac[i]);
6375 fprintf(fplog, "\n");
6382 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6385 gmx_mtop_ilistloop_t iloop;
6389 iloop = gmx_mtop_ilistloop_init(mtop);
6390 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6392 for (ftype = 0; ftype < F_NRE; ftype++)
6394 if ((interaction_function[ftype].flags & IF_BOND) &&
6397 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6405 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6411 val = getenv(env_var);
6414 if (sscanf(val, "%20d", &nst) <= 0)
6420 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6428 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6432 fprintf(stderr, "\n%s\n", warn_string);
6436 fprintf(fplog, "\n%s\n", warn_string);
6440 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6441 t_inputrec *ir, FILE *fplog)
6443 if (ir->ePBC == epbcSCREW &&
6444 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6446 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6449 if (ir->ns_type == ensSIMPLE)
6451 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6454 if (ir->nstlist == 0)
6456 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6459 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6461 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6465 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6470 r = ddbox->box_size[XX];
6471 for (di = 0; di < dd->ndim; di++)
6474 /* Check using the initial average cell size */
6475 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6481 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6482 const char *dlb_opt, gmx_bool bRecordLoad,
6483 unsigned long Flags, t_inputrec *ir)
6490 case 'a': eDLB = edlbAUTO; break;
6491 case 'n': eDLB = edlbNO; break;
6492 case 'y': eDLB = edlbYES; break;
6493 default: gmx_incons("Unknown dlb_opt");
6496 if (Flags & MD_RERUN)
6501 if (!EI_DYNAMICS(ir->eI))
6503 if (eDLB == edlbYES)
6505 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6506 dd_warning(cr, fplog, buf);
6514 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6519 if (Flags & MD_REPRODUCIBLE)
6526 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6530 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6533 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6541 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6546 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6548 /* Decomposition order z,y,x */
6551 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6553 for (dim = DIM-1; dim >= 0; dim--)
6555 if (dd->nc[dim] > 1)
6557 dd->dim[dd->ndim++] = dim;
6563 /* Decomposition order x,y,z */
6564 for (dim = 0; dim < DIM; dim++)
6566 if (dd->nc[dim] > 1)
6568 dd->dim[dd->ndim++] = dim;
6574 static gmx_domdec_comm_t *init_dd_comm()
6576 gmx_domdec_comm_t *comm;
6580 snew(comm->cggl_flag, DIM*2);
6581 snew(comm->cgcm_state, DIM*2);
6582 for (i = 0; i < DIM*2; i++)
6584 comm->cggl_flag_nalloc[i] = 0;
6585 comm->cgcm_state_nalloc[i] = 0;
6588 comm->nalloc_int = 0;
6589 comm->buf_int = NULL;
6591 vec_rvec_init(&comm->vbuf);
6593 comm->n_load_have = 0;
6594 comm->n_load_collect = 0;
6596 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6598 comm->sum_nat[i] = 0;
6602 comm->load_step = 0;
6605 clear_ivec(comm->load_lim);
6612 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6613 unsigned long Flags,
6615 real comm_distance_min, real rconstr,
6616 const char *dlb_opt, real dlb_scale,
6617 const char *sizex, const char *sizey, const char *sizez,
6618 gmx_mtop_t *mtop, t_inputrec *ir,
6619 matrix box, rvec *x,
6621 int *npme_x, int *npme_y)
6624 gmx_domdec_comm_t *comm;
6626 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6629 const real tenPercentMargin = 1.1;
6634 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6639 dd->comm = init_dd_comm();
6641 snew(comm->cggl_flag, DIM*2);
6642 snew(comm->cgcm_state, DIM*2);
6644 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6645 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6647 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6648 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6649 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6650 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6651 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6652 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6653 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6654 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6656 dd->pme_recv_f_alloc = 0;
6657 dd->pme_recv_f_buf = NULL;
6659 if (dd->bSendRecv2 && fplog)
6661 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");
6667 fprintf(fplog, "Will load balance based on FLOP count\n");
6669 if (comm->eFlop > 1)
6671 srand(1+cr->nodeid);
6673 comm->bRecordLoad = TRUE;
6677 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6681 /* Initialize to GPU share count to 0, might change later */
6682 comm->nrank_gpu_shared = 0;
6684 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6685 comm->bDLB_locked = FALSE;
6687 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6690 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6692 dd->bGridJump = comm->bDynLoadBal;
6693 comm->bPMELoadBalDLBLimits = FALSE;
6695 if (comm->nstSortCG)
6699 if (comm->nstSortCG == 1)
6701 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6705 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6709 snew(comm->sort, 1);
6715 fprintf(fplog, "Will not sort the charge groups\n");
6719 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6721 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6722 if (comm->bInterCGBondeds)
6724 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6728 comm->bInterCGMultiBody = FALSE;
6731 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6732 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6734 if (ir->rlistlong == 0)
6736 /* Set the cut-off to some very large value,
6737 * so we don't need if statements everywhere in the code.
6738 * We use sqrt, since the cut-off is squared in some places.
6740 comm->cutoff = GMX_CUTOFF_INF;
6744 comm->cutoff = ir->rlistlong;
6746 comm->cutoff_mbody = 0;
6748 comm->cellsize_limit = 0;
6749 comm->bBondComm = FALSE;
6751 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6752 * within nstlist steps. Since boundaries are allowed to displace by half
6753 * a cell size, DD cells should be at least the size of the list buffer.
6755 comm->cellsize_limit = std::max(comm->cellsize_limit,
6756 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6758 if (comm->bInterCGBondeds)
6760 if (comm_distance_min > 0)
6762 comm->cutoff_mbody = comm_distance_min;
6763 if (Flags & MD_DDBONDCOMM)
6765 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6769 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6771 r_bonded_limit = comm->cutoff_mbody;
6773 else if (ir->bPeriodicMols)
6775 /* Can not easily determine the required cut-off */
6776 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");
6777 comm->cutoff_mbody = comm->cutoff/2;
6778 r_bonded_limit = comm->cutoff_mbody;
6784 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6785 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6787 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6788 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6790 /* We use an initial margin of 10% for the minimum cell size,
6791 * except when we are just below the non-bonded cut-off.
6793 if (Flags & MD_DDBONDCOMM)
6795 if (std::max(r_2b, r_mb) > comm->cutoff)
6797 r_bonded = std::max(r_2b, r_mb);
6798 r_bonded_limit = tenPercentMargin*r_bonded;
6799 comm->bBondComm = TRUE;
6804 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6806 /* We determine cutoff_mbody later */
6810 /* No special bonded communication,
6811 * simply increase the DD cut-off.
6813 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6814 comm->cutoff_mbody = r_bonded_limit;
6815 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6818 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6822 "Minimum cell size due to bonded interactions: %.3f nm\n",
6823 comm->cellsize_limit);
6827 if (dd->bInterCGcons && rconstr <= 0)
6829 /* There is a cell size limit due to the constraints (P-LINCS) */
6830 rconstr = constr_r_max(fplog, mtop, ir);
6834 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6836 if (rconstr > comm->cellsize_limit)
6838 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6842 else if (rconstr > 0 && fplog)
6844 /* Here we do not check for dd->bInterCGcons,
6845 * because one can also set a cell size limit for virtual sites only
6846 * and at this point we don't know yet if there are intercg v-sites.
6849 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6852 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6854 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6858 copy_ivec(nc, dd->nc);
6859 set_dd_dim(fplog, dd);
6860 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6862 if (cr->npmenodes == -1)
6866 acs = average_cellsize_min(dd, ddbox);
6867 if (acs < comm->cellsize_limit)
6871 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6873 gmx_fatal_collective(FARGS, cr, NULL,
6874 "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",
6875 acs, comm->cellsize_limit);
6880 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6882 /* We need to choose the optimal DD grid and possibly PME nodes */
6883 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6884 comm->eDLB != edlbNO, dlb_scale,
6885 comm->cellsize_limit, comm->cutoff,
6886 comm->bInterCGBondeds);
6888 if (dd->nc[XX] == 0)
6890 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6891 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6892 !bC ? "-rdd" : "-rcon",
6893 comm->eDLB != edlbNO ? " or -dds" : "",
6894 bC ? " or your LINCS settings" : "");
6896 gmx_fatal_collective(FARGS, cr, NULL,
6897 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6899 "Look in the log file for details on the domain decomposition",
6900 cr->nnodes-cr->npmenodes, limit, buf);
6902 set_dd_dim(fplog, dd);
6908 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6909 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6912 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6913 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6915 gmx_fatal_collective(FARGS, cr, NULL,
6916 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6917 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6919 if (cr->npmenodes > dd->nnodes)
6921 gmx_fatal_collective(FARGS, cr, NULL,
6922 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6924 if (cr->npmenodes > 0)
6926 comm->npmenodes = cr->npmenodes;
6930 comm->npmenodes = dd->nnodes;
6933 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6935 /* The following choices should match those
6936 * in comm_cost_est in domdec_setup.c.
6937 * Note that here the checks have to take into account
6938 * that the decomposition might occur in a different order than xyz
6939 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6940 * in which case they will not match those in comm_cost_est,
6941 * but since that is mainly for testing purposes that's fine.
6943 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6944 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6945 getenv("GMX_PMEONEDD") == NULL)
6947 comm->npmedecompdim = 2;
6948 comm->npmenodes_x = dd->nc[XX];
6949 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6953 /* In case nc is 1 in both x and y we could still choose to
6954 * decompose pme in y instead of x, but we use x for simplicity.
6956 comm->npmedecompdim = 1;
6957 if (dd->dim[0] == YY)
6959 comm->npmenodes_x = 1;
6960 comm->npmenodes_y = comm->npmenodes;
6964 comm->npmenodes_x = comm->npmenodes;
6965 comm->npmenodes_y = 1;
6970 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6971 comm->npmenodes_x, comm->npmenodes_y, 1);
6976 comm->npmedecompdim = 0;
6977 comm->npmenodes_x = 0;
6978 comm->npmenodes_y = 0;
6981 /* Technically we don't need both of these,
6982 * but it simplifies code not having to recalculate it.
6984 *npme_x = comm->npmenodes_x;
6985 *npme_y = comm->npmenodes_y;
6987 snew(comm->slb_frac, DIM);
6988 if (comm->eDLB == edlbNO)
6990 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6991 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6992 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6995 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6997 if (comm->bBondComm || comm->eDLB != edlbNO)
6999 /* Set the bonded communication distance to halfway
7000 * the minimum and the maximum,
7001 * since the extra communication cost is nearly zero.
7003 acs = average_cellsize_min(dd, ddbox);
7004 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7005 if (comm->eDLB != edlbNO)
7007 /* Check if this does not limit the scaling */
7008 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7010 if (!comm->bBondComm)
7012 /* Without bBondComm do not go beyond the n.b. cut-off */
7013 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7014 if (comm->cellsize_limit >= comm->cutoff)
7016 /* We don't loose a lot of efficieny
7017 * when increasing it to the n.b. cut-off.
7018 * It can even be slightly faster, because we need
7019 * less checks for the communication setup.
7021 comm->cutoff_mbody = comm->cutoff;
7024 /* Check if we did not end up below our original limit */
7025 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7027 if (comm->cutoff_mbody > comm->cellsize_limit)
7029 comm->cellsize_limit = comm->cutoff_mbody;
7032 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7037 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7038 "cellsize limit %f\n",
7039 comm->bBondComm, comm->cellsize_limit);
7044 check_dd_restrictions(cr, dd, ir, fplog);
7047 comm->partition_step = INT_MIN;
7050 clear_dd_cycle_counts(dd);
7055 static void set_dlb_limits(gmx_domdec_t *dd)
7060 for (d = 0; d < dd->ndim; d++)
7062 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7063 dd->comm->cellsize_min[dd->dim[d]] =
7064 dd->comm->cellsize_min_dlb[dd->dim[d]];
7069 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7072 gmx_domdec_comm_t *comm;
7082 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);
7085 cellsize_min = comm->cellsize_min[dd->dim[0]];
7086 for (d = 1; d < dd->ndim; d++)
7088 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7091 if (cellsize_min < comm->cellsize_limit*1.05)
7093 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");
7095 /* Change DLB from "auto" to "no". */
7096 comm->eDLB = edlbNO;
7101 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7102 comm->bDynLoadBal = TRUE;
7103 dd->bGridJump = TRUE;
7107 /* We can set the required cell size info here,
7108 * so we do not need to communicate this.
7109 * The grid is completely uniform.
7111 for (d = 0; d < dd->ndim; d++)
7115 comm->load[d].sum_m = comm->load[d].sum;
7117 nc = dd->nc[dd->dim[d]];
7118 for (i = 0; i < nc; i++)
7120 comm->root[d]->cell_f[i] = i/(real)nc;
7123 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7124 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7127 comm->root[d]->cell_f[nc] = 1.0;
7132 static char *init_bLocalCG(gmx_mtop_t *mtop)
7137 ncg = ncg_mtop(mtop);
7138 snew(bLocalCG, ncg);
7139 for (cg = 0; cg < ncg; cg++)
7141 bLocalCG[cg] = FALSE;
7147 void dd_init_bondeds(FILE *fplog,
7148 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7150 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7152 gmx_domdec_comm_t *comm;
7154 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7158 if (comm->bBondComm)
7160 /* Communicate atoms beyond the cut-off for bonded interactions */
7163 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7165 comm->bLocalCG = init_bLocalCG(mtop);
7169 /* Only communicate atoms based on cut-off */
7170 comm->cglink = NULL;
7171 comm->bLocalCG = NULL;
7175 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7177 gmx_bool bDynLoadBal, real dlb_scale,
7180 gmx_domdec_comm_t *comm;
7195 fprintf(fplog, "The maximum number of communication pulses is:");
7196 for (d = 0; d < dd->ndim; d++)
7198 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7200 fprintf(fplog, "\n");
7201 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7202 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7203 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7204 for (d = 0; d < DIM; d++)
7208 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7215 comm->cellsize_min_dlb[d]/
7216 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7218 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7221 fprintf(fplog, "\n");
7225 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7226 fprintf(fplog, "The initial number of communication pulses is:");
7227 for (d = 0; d < dd->ndim; d++)
7229 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7231 fprintf(fplog, "\n");
7232 fprintf(fplog, "The initial domain decomposition cell size is:");
7233 for (d = 0; d < DIM; d++)
7237 fprintf(fplog, " %c %.2f nm",
7238 dim2char(d), dd->comm->cellsize_min[d]);
7241 fprintf(fplog, "\n\n");
7244 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7246 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7247 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7248 "non-bonded interactions", "", comm->cutoff);
7252 limit = dd->comm->cellsize_limit;
7256 if (dynamic_dd_box(ddbox, ir))
7258 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7260 limit = dd->comm->cellsize_min[XX];
7261 for (d = 1; d < DIM; d++)
7263 limit = std::min(limit, dd->comm->cellsize_min[d]);
7267 if (comm->bInterCGBondeds)
7269 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7270 "two-body bonded interactions", "(-rdd)",
7271 std::max(comm->cutoff, comm->cutoff_mbody));
7272 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7273 "multi-body bonded interactions", "(-rdd)",
7274 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7278 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7279 "virtual site constructions", "(-rcon)", limit);
7281 if (dd->constraint_comm)
7283 sprintf(buf, "atoms separated by up to %d constraints",
7285 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7286 buf, "(-rcon)", limit);
7288 fprintf(fplog, "\n");
7294 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7296 const t_inputrec *ir,
7297 const gmx_ddbox_t *ddbox)
7299 gmx_domdec_comm_t *comm;
7300 int d, dim, npulse, npulse_d_max, npulse_d;
7305 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7307 /* Determine the maximum number of comm. pulses in one dimension */
7309 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7311 /* Determine the maximum required number of grid pulses */
7312 if (comm->cellsize_limit >= comm->cutoff)
7314 /* Only a single pulse is required */
7317 else if (!bNoCutOff && comm->cellsize_limit > 0)
7319 /* We round down slightly here to avoid overhead due to the latency
7320 * of extra communication calls when the cut-off
7321 * would be only slightly longer than the cell size.
7322 * Later cellsize_limit is redetermined,
7323 * so we can not miss interactions due to this rounding.
7325 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7329 /* There is no cell size limit */
7330 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7333 if (!bNoCutOff && npulse > 1)
7335 /* See if we can do with less pulses, based on dlb_scale */
7337 for (d = 0; d < dd->ndim; d++)
7340 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7341 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7342 npulse_d_max = std::max(npulse_d_max, npulse_d);
7344 npulse = std::min(npulse, npulse_d_max);
7347 /* This env var can override npulse */
7348 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7355 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7356 for (d = 0; d < dd->ndim; d++)
7358 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7359 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7360 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7361 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7362 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7364 comm->bVacDLBNoLimit = FALSE;
7368 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7369 if (!comm->bVacDLBNoLimit)
7371 comm->cellsize_limit = std::max(comm->cellsize_limit,
7372 comm->cutoff/comm->maxpulse);
7374 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7375 /* Set the minimum cell size for each DD dimension */
7376 for (d = 0; d < dd->ndim; d++)
7378 if (comm->bVacDLBNoLimit ||
7379 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7381 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7385 comm->cellsize_min_dlb[dd->dim[d]] =
7386 comm->cutoff/comm->cd[d].np_dlb;
7389 if (comm->cutoff_mbody <= 0)
7391 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7393 if (comm->bDynLoadBal)
7399 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7401 /* If each molecule is a single charge group
7402 * or we use domain decomposition for each periodic dimension,
7403 * we do not need to take pbc into account for the bonded interactions.
7405 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7408 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7411 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7412 t_inputrec *ir, gmx_ddbox_t *ddbox)
7414 gmx_domdec_comm_t *comm;
7420 /* Initialize the thread data.
7421 * This can not be done in init_domain_decomposition,
7422 * as the numbers of threads is determined later.
7424 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7427 snew(comm->dth, comm->nth);
7430 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7432 init_ddpme(dd, &comm->ddpme[0], 0);
7433 if (comm->npmedecompdim >= 2)
7435 init_ddpme(dd, &comm->ddpme[1], 1);
7440 comm->npmenodes = 0;
7441 if (dd->pme_nodeid >= 0)
7443 gmx_fatal_collective(FARGS, NULL, dd,
7444 "Can not have separate PME ranks without PME electrostatics");
7450 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7452 if (comm->eDLB != edlbNO)
7454 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7457 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7458 if (comm->eDLB == edlbAUTO)
7462 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7464 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7467 if (ir->ePBC == epbcNONE)
7469 vol_frac = 1 - 1/(double)dd->nnodes;
7474 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7478 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7480 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7482 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7485 static gmx_bool test_dd_cutoff(t_commrec *cr,
7486 t_state *state, t_inputrec *ir,
7497 set_ddbox(dd, FALSE, cr, ir, state->box,
7498 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7502 for (d = 0; d < dd->ndim; d++)
7506 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7507 if (dynamic_dd_box(&ddbox, ir))
7509 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7512 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7514 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7515 dd->comm->cd[d].np_dlb > 0)
7517 if (np > dd->comm->cd[d].np_dlb)
7522 /* If a current local cell size is smaller than the requested
7523 * cut-off, we could still fix it, but this gets very complicated.
7524 * Without fixing here, we might actually need more checks.
7526 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7533 if (dd->comm->eDLB != edlbNO)
7535 /* If DLB is not active yet, we don't need to check the grid jumps.
7536 * Actually we shouldn't, because then the grid jump data is not set.
7538 if (dd->comm->bDynLoadBal &&
7539 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7544 gmx_sumi(1, &LocallyLimited, cr);
7546 if (LocallyLimited > 0)
7555 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7558 gmx_bool bCutoffAllowed;
7560 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7564 cr->dd->comm->cutoff = cutoff_req;
7567 return bCutoffAllowed;
7570 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7572 gmx_domdec_comm_t *comm;
7574 comm = cr->dd->comm;
7576 /* Turn on the DLB limiting (might have been on already) */
7577 comm->bPMELoadBalDLBLimits = TRUE;
7579 /* Change the cut-off limit */
7580 comm->PMELoadBal_max_cutoff = comm->cutoff;
7583 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7585 return dd->comm->bDLB_locked;
7588 void dd_dlb_set_lock(gmx_domdec_t *dd, gmx_bool bValue)
7590 /* We can only lock the DLB when it is set to auto, otherwise don't lock */
7591 if (dd->comm->eDLB == edlbAUTO)
7593 dd->comm->bDLB_locked = bValue;
7597 static void merge_cg_buffers(int ncell,
7598 gmx_domdec_comm_dim_t *cd, int pulse,
7600 int *index_gl, int *recv_i,
7601 rvec *cg_cm, rvec *recv_vr,
7603 cginfo_mb_t *cginfo_mb, int *cginfo)
7605 gmx_domdec_ind_t *ind, *ind_p;
7606 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7607 int shift, shift_at;
7609 ind = &cd->ind[pulse];
7611 /* First correct the already stored data */
7612 shift = ind->nrecv[ncell];
7613 for (cell = ncell-1; cell >= 0; cell--)
7615 shift -= ind->nrecv[cell];
7618 /* Move the cg's present from previous grid pulses */
7619 cg0 = ncg_cell[ncell+cell];
7620 cg1 = ncg_cell[ncell+cell+1];
7621 cgindex[cg1+shift] = cgindex[cg1];
7622 for (cg = cg1-1; cg >= cg0; cg--)
7624 index_gl[cg+shift] = index_gl[cg];
7625 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7626 cgindex[cg+shift] = cgindex[cg];
7627 cginfo[cg+shift] = cginfo[cg];
7629 /* Correct the already stored send indices for the shift */
7630 for (p = 1; p <= pulse; p++)
7632 ind_p = &cd->ind[p];
7634 for (c = 0; c < cell; c++)
7636 cg0 += ind_p->nsend[c];
7638 cg1 = cg0 + ind_p->nsend[cell];
7639 for (cg = cg0; cg < cg1; cg++)
7641 ind_p->index[cg] += shift;
7647 /* Merge in the communicated buffers */
7651 for (cell = 0; cell < ncell; cell++)
7653 cg1 = ncg_cell[ncell+cell+1] + shift;
7656 /* Correct the old cg indices */
7657 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7659 cgindex[cg+1] += shift_at;
7662 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7664 /* Copy this charge group from the buffer */
7665 index_gl[cg1] = recv_i[cg0];
7666 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7667 /* Add it to the cgindex */
7668 cg_gl = index_gl[cg1];
7669 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7670 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7671 cgindex[cg1+1] = cgindex[cg1] + nat;
7676 shift += ind->nrecv[cell];
7677 ncg_cell[ncell+cell+1] = cg1;
7681 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7682 int nzone, int cg0, const int *cgindex)
7686 /* Store the atom block boundaries for easy copying of communication buffers
7689 for (zone = 0; zone < nzone; zone++)
7691 for (p = 0; p < cd->np; p++)
7693 cd->ind[p].cell2at0[zone] = cgindex[cg];
7694 cg += cd->ind[p].nrecv[zone];
7695 cd->ind[p].cell2at1[zone] = cgindex[cg];
7700 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7706 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7708 if (!bLocalCG[link->a[i]])
7717 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7719 real c[DIM][4]; /* the corners for the non-bonded communication */
7720 real cr0; /* corner for rounding */
7721 real cr1[4]; /* corners for rounding */
7722 real bc[DIM]; /* corners for bounded communication */
7723 real bcr1; /* corner for rounding for bonded communication */
7726 /* Determine the corners of the domain(s) we are communicating with */
7728 set_dd_corners(const gmx_domdec_t *dd,
7729 int dim0, int dim1, int dim2,
7733 const gmx_domdec_comm_t *comm;
7734 const gmx_domdec_zones_t *zones;
7739 zones = &comm->zones;
7741 /* Keep the compiler happy */
7745 /* The first dimension is equal for all cells */
7746 c->c[0][0] = comm->cell_x0[dim0];
7749 c->bc[0] = c->c[0][0];
7754 /* This cell row is only seen from the first row */
7755 c->c[1][0] = comm->cell_x0[dim1];
7756 /* All rows can see this row */
7757 c->c[1][1] = comm->cell_x0[dim1];
7760 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7763 /* For the multi-body distance we need the maximum */
7764 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7767 /* Set the upper-right corner for rounding */
7768 c->cr0 = comm->cell_x1[dim0];
7773 for (j = 0; j < 4; j++)
7775 c->c[2][j] = comm->cell_x0[dim2];
7779 /* Use the maximum of the i-cells that see a j-cell */
7780 for (i = 0; i < zones->nizone; i++)
7782 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7787 std::max(c->c[2][j-4],
7788 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7794 /* For the multi-body distance we need the maximum */
7795 c->bc[2] = comm->cell_x0[dim2];
7796 for (i = 0; i < 2; i++)
7798 for (j = 0; j < 2; j++)
7800 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7806 /* Set the upper-right corner for rounding */
7807 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7808 * Only cell (0,0,0) can see cell 7 (1,1,1)
7810 c->cr1[0] = comm->cell_x1[dim1];
7811 c->cr1[3] = comm->cell_x1[dim1];
7814 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7817 /* For the multi-body distance we need the maximum */
7818 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7825 /* Determine which cg's we need to send in this pulse from this zone */
7827 get_zone_pulse_cgs(gmx_domdec_t *dd,
7828 int zonei, int zone,
7830 const int *index_gl,
7832 int dim, int dim_ind,
7833 int dim0, int dim1, int dim2,
7834 real r_comm2, real r_bcomm2,
7838 real skew_fac2_d, real skew_fac_01,
7839 rvec *v_d, rvec *v_0, rvec *v_1,
7840 const dd_corners_t *c,
7842 gmx_bool bDistBonded,
7848 gmx_domdec_ind_t *ind,
7849 int **ibuf, int *ibuf_nalloc,
7855 gmx_domdec_comm_t *comm;
7857 gmx_bool bDistMB_pulse;
7859 real r2, rb2, r, tric_sh;
7862 int nsend_z, nsend, nat;
7866 bScrew = (dd->bScrewPBC && dim == XX);
7868 bDistMB_pulse = (bDistMB && bDistBonded);
7874 for (cg = cg0; cg < cg1; cg++)
7878 if (tric_dist[dim_ind] == 0)
7880 /* Rectangular direction, easy */
7881 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7888 r = cg_cm[cg][dim] - c->bc[dim_ind];
7894 /* Rounding gives at most a 16% reduction
7895 * in communicated atoms
7897 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7899 r = cg_cm[cg][dim0] - c->cr0;
7900 /* This is the first dimension, so always r >= 0 */
7907 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7909 r = cg_cm[cg][dim1] - c->cr1[zone];
7916 r = cg_cm[cg][dim1] - c->bcr1;
7926 /* Triclinic direction, more complicated */
7929 /* Rounding, conservative as the skew_fac multiplication
7930 * will slightly underestimate the distance.
7932 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7934 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7935 for (i = dim0+1; i < DIM; i++)
7937 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7939 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7942 rb[dim0] = rn[dim0];
7945 /* Take care that the cell planes along dim0 might not
7946 * be orthogonal to those along dim1 and dim2.
7948 for (i = 1; i <= dim_ind; i++)
7951 if (normal[dim0][dimd] > 0)
7953 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7956 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7961 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7963 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7965 for (i = dim1+1; i < DIM; i++)
7967 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7969 rn[dim1] += tric_sh;
7972 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7973 /* Take care of coupling of the distances
7974 * to the planes along dim0 and dim1 through dim2.
7976 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7977 /* Take care that the cell planes along dim1
7978 * might not be orthogonal to that along dim2.
7980 if (normal[dim1][dim2] > 0)
7982 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7988 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7991 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7992 /* Take care of coupling of the distances
7993 * to the planes along dim0 and dim1 through dim2.
7995 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7996 /* Take care that the cell planes along dim1
7997 * might not be orthogonal to that along dim2.
7999 if (normal[dim1][dim2] > 0)
8001 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8006 /* The distance along the communication direction */
8007 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8009 for (i = dim+1; i < DIM; i++)
8011 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8016 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8017 /* Take care of coupling of the distances
8018 * to the planes along dim0 and dim1 through dim2.
8020 if (dim_ind == 1 && zonei == 1)
8022 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8028 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8031 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8032 /* Take care of coupling of the distances
8033 * to the planes along dim0 and dim1 through dim2.
8035 if (dim_ind == 1 && zonei == 1)
8037 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8045 ((bDistMB && rb2 < r_bcomm2) ||
8046 (bDist2B && r2 < r_bcomm2)) &&
8048 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8049 missing_link(comm->cglink, index_gl[cg],
8052 /* Make an index to the local charge groups */
8053 if (nsend+1 > ind->nalloc)
8055 ind->nalloc = over_alloc_large(nsend+1);
8056 srenew(ind->index, ind->nalloc);
8058 if (nsend+1 > *ibuf_nalloc)
8060 *ibuf_nalloc = over_alloc_large(nsend+1);
8061 srenew(*ibuf, *ibuf_nalloc);
8063 ind->index[nsend] = cg;
8064 (*ibuf)[nsend] = index_gl[cg];
8066 vec_rvec_check_alloc(vbuf, nsend+1);
8068 if (dd->ci[dim] == 0)
8070 /* Correct cg_cm for pbc */
8071 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8074 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8075 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8080 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8083 nat += cgindex[cg+1] - cgindex[cg];
8089 *nsend_z_ptr = nsend_z;
8092 static void setup_dd_communication(gmx_domdec_t *dd,
8093 matrix box, gmx_ddbox_t *ddbox,
8094 t_forcerec *fr, t_state *state, rvec **f)
8096 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8097 int nzone, nzone_send, zone, zonei, cg0, cg1;
8098 int c, i, cg, cg_gl, nrcg;
8099 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8100 gmx_domdec_comm_t *comm;
8101 gmx_domdec_zones_t *zones;
8102 gmx_domdec_comm_dim_t *cd;
8103 gmx_domdec_ind_t *ind;
8104 cginfo_mb_t *cginfo_mb;
8105 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8106 real r_comm2, r_bcomm2;
8107 dd_corners_t corners;
8109 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8110 real skew_fac2_d, skew_fac_01;
8117 fprintf(debug, "Setting up DD communication\n");
8122 switch (fr->cutoff_scheme)
8131 gmx_incons("unimplemented");
8135 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8137 /* Check if we need to use triclinic distances */
8138 tric_dist[dim_ind] = 0;
8139 for (i = 0; i <= dim_ind; i++)
8141 if (ddbox->tric_dir[dd->dim[i]])
8143 tric_dist[dim_ind] = 1;
8148 bBondComm = comm->bBondComm;
8150 /* Do we need to determine extra distances for multi-body bondeds? */
8151 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8153 /* Do we need to determine extra distances for only two-body bondeds? */
8154 bDist2B = (bBondComm && !bDistMB);
8156 r_comm2 = sqr(comm->cutoff);
8157 r_bcomm2 = sqr(comm->cutoff_mbody);
8161 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8164 zones = &comm->zones;
8167 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8168 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8170 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8172 /* Triclinic stuff */
8173 normal = ddbox->normal;
8177 v_0 = ddbox->v[dim0];
8178 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8180 /* Determine the coupling coefficient for the distances
8181 * to the cell planes along dim0 and dim1 through dim2.
8182 * This is required for correct rounding.
8185 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8188 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8194 v_1 = ddbox->v[dim1];
8197 zone_cg_range = zones->cg_range;
8198 index_gl = dd->index_gl;
8199 cgindex = dd->cgindex;
8200 cginfo_mb = fr->cginfo_mb;
8202 zone_cg_range[0] = 0;
8203 zone_cg_range[1] = dd->ncg_home;
8204 comm->zone_ncg1[0] = dd->ncg_home;
8205 pos_cg = dd->ncg_home;
8207 nat_tot = dd->nat_home;
8209 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8211 dim = dd->dim[dim_ind];
8212 cd = &comm->cd[dim_ind];
8214 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8216 /* No pbc in this dimension, the first node should not comm. */
8224 v_d = ddbox->v[dim];
8225 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8227 cd->bInPlace = TRUE;
8228 for (p = 0; p < cd->np; p++)
8230 /* Only atoms communicated in the first pulse are used
8231 * for multi-body bonded interactions or for bBondComm.
8233 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8238 for (zone = 0; zone < nzone_send; zone++)
8240 if (tric_dist[dim_ind] && dim_ind > 0)
8242 /* Determine slightly more optimized skew_fac's
8244 * This reduces the number of communicated atoms
8245 * by about 10% for 3D DD of rhombic dodecahedra.
8247 for (dimd = 0; dimd < dim; dimd++)
8249 sf2_round[dimd] = 1;
8250 if (ddbox->tric_dir[dimd])
8252 for (i = dd->dim[dimd]+1; i < DIM; i++)
8254 /* If we are shifted in dimension i
8255 * and the cell plane is tilted forward
8256 * in dimension i, skip this coupling.
8258 if (!(zones->shift[nzone+zone][i] &&
8259 ddbox->v[dimd][i][dimd] >= 0))
8262 sqr(ddbox->v[dimd][i][dimd]);
8265 sf2_round[dimd] = 1/sf2_round[dimd];
8270 zonei = zone_perm[dim_ind][zone];
8273 /* Here we permutate the zones to obtain a convenient order
8274 * for neighbor searching
8276 cg0 = zone_cg_range[zonei];
8277 cg1 = zone_cg_range[zonei+1];
8281 /* Look only at the cg's received in the previous grid pulse
8283 cg1 = zone_cg_range[nzone+zone+1];
8284 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8287 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8288 for (th = 0; th < comm->nth; th++)
8290 gmx_domdec_ind_t *ind_p;
8291 int **ibuf_p, *ibuf_nalloc_p;
8293 int *nsend_p, *nat_p;
8299 /* Thread 0 writes in the comm buffers */
8301 ibuf_p = &comm->buf_int;
8302 ibuf_nalloc_p = &comm->nalloc_int;
8303 vbuf_p = &comm->vbuf;
8306 nsend_zone_p = &ind->nsend[zone];
8310 /* Other threads write into temp buffers */
8311 ind_p = &comm->dth[th].ind;
8312 ibuf_p = &comm->dth[th].ibuf;
8313 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8314 vbuf_p = &comm->dth[th].vbuf;
8315 nsend_p = &comm->dth[th].nsend;
8316 nat_p = &comm->dth[th].nat;
8317 nsend_zone_p = &comm->dth[th].nsend_zone;
8319 comm->dth[th].nsend = 0;
8320 comm->dth[th].nat = 0;
8321 comm->dth[th].nsend_zone = 0;
8331 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8332 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8335 /* Get the cg's for this pulse in this zone */
8336 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8338 dim, dim_ind, dim0, dim1, dim2,
8341 normal, skew_fac2_d, skew_fac_01,
8342 v_d, v_0, v_1, &corners, sf2_round,
8343 bDistBonded, bBondComm,
8347 ibuf_p, ibuf_nalloc_p,
8353 /* Append data of threads>=1 to the communication buffers */
8354 for (th = 1; th < comm->nth; th++)
8356 dd_comm_setup_work_t *dth;
8359 dth = &comm->dth[th];
8361 ns1 = nsend + dth->nsend_zone;
8362 if (ns1 > ind->nalloc)
8364 ind->nalloc = over_alloc_dd(ns1);
8365 srenew(ind->index, ind->nalloc);
8367 if (ns1 > comm->nalloc_int)
8369 comm->nalloc_int = over_alloc_dd(ns1);
8370 srenew(comm->buf_int, comm->nalloc_int);
8372 if (ns1 > comm->vbuf.nalloc)
8374 comm->vbuf.nalloc = over_alloc_dd(ns1);
8375 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8378 for (i = 0; i < dth->nsend_zone; i++)
8380 ind->index[nsend] = dth->ind.index[i];
8381 comm->buf_int[nsend] = dth->ibuf[i];
8382 copy_rvec(dth->vbuf.v[i],
8383 comm->vbuf.v[nsend]);
8387 ind->nsend[zone] += dth->nsend_zone;
8390 /* Clear the counts in case we do not have pbc */
8391 for (zone = nzone_send; zone < nzone; zone++)
8393 ind->nsend[zone] = 0;
8395 ind->nsend[nzone] = nsend;
8396 ind->nsend[nzone+1] = nat;
8397 /* Communicate the number of cg's and atoms to receive */
8398 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8399 ind->nsend, nzone+2,
8400 ind->nrecv, nzone+2);
8402 /* The rvec buffer is also required for atom buffers of size nsend
8403 * in dd_move_x and dd_move_f.
8405 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8409 /* We can receive in place if only the last zone is not empty */
8410 for (zone = 0; zone < nzone-1; zone++)
8412 if (ind->nrecv[zone] > 0)
8414 cd->bInPlace = FALSE;
8419 /* The int buffer is only required here for the cg indices */
8420 if (ind->nrecv[nzone] > comm->nalloc_int2)
8422 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8423 srenew(comm->buf_int2, comm->nalloc_int2);
8425 /* The rvec buffer is also required for atom buffers
8426 * of size nrecv in dd_move_x and dd_move_f.
8428 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8429 vec_rvec_check_alloc(&comm->vbuf2, i);
8433 /* Make space for the global cg indices */
8434 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8435 || dd->cg_nalloc == 0)
8437 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8438 srenew(index_gl, dd->cg_nalloc);
8439 srenew(cgindex, dd->cg_nalloc+1);
8441 /* Communicate the global cg indices */
8444 recv_i = index_gl + pos_cg;
8448 recv_i = comm->buf_int2;
8450 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8451 comm->buf_int, nsend,
8452 recv_i, ind->nrecv[nzone]);
8454 /* Make space for cg_cm */
8455 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8456 if (fr->cutoff_scheme == ecutsGROUP)
8464 /* Communicate cg_cm */
8467 recv_vr = cg_cm + pos_cg;
8471 recv_vr = comm->vbuf2.v;
8473 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8474 comm->vbuf.v, nsend,
8475 recv_vr, ind->nrecv[nzone]);
8477 /* Make the charge group index */
8480 zone = (p == 0 ? 0 : nzone - 1);
8481 while (zone < nzone)
8483 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8485 cg_gl = index_gl[pos_cg];
8486 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8487 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8488 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8491 /* Update the charge group presence,
8492 * so we can use it in the next pass of the loop.
8494 comm->bLocalCG[cg_gl] = TRUE;
8500 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8503 zone_cg_range[nzone+zone] = pos_cg;
8508 /* This part of the code is never executed with bBondComm. */
8509 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8510 index_gl, recv_i, cg_cm, recv_vr,
8511 cgindex, fr->cginfo_mb, fr->cginfo);
8512 pos_cg += ind->nrecv[nzone];
8514 nat_tot += ind->nrecv[nzone+1];
8518 /* Store the atom block for easy copying of communication buffers */
8519 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8523 dd->index_gl = index_gl;
8524 dd->cgindex = cgindex;
8526 dd->ncg_tot = zone_cg_range[zones->n];
8527 dd->nat_tot = nat_tot;
8528 comm->nat[ddnatHOME] = dd->nat_home;
8529 for (i = ddnatZONE; i < ddnatNR; i++)
8531 comm->nat[i] = dd->nat_tot;
8536 /* We don't need to update cginfo, since that was alrady done above.
8537 * So we pass NULL for the forcerec.
8539 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8540 NULL, comm->bLocalCG);
8545 fprintf(debug, "Finished setting up DD communication, zones:");
8546 for (c = 0; c < zones->n; c++)
8548 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8550 fprintf(debug, "\n");
8554 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8558 for (c = 0; c < zones->nizone; c++)
8560 zones->izone[c].cg1 = zones->cg_range[c+1];
8561 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8562 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8566 static void set_zones_size(gmx_domdec_t *dd,
8567 matrix box, const gmx_ddbox_t *ddbox,
8568 int zone_start, int zone_end)
8570 gmx_domdec_comm_t *comm;
8571 gmx_domdec_zones_t *zones;
8580 zones = &comm->zones;
8582 /* Do we need to determine extra distances for multi-body bondeds? */
8583 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8585 for (z = zone_start; z < zone_end; z++)
8587 /* Copy cell limits to zone limits.
8588 * Valid for non-DD dims and non-shifted dims.
8590 copy_rvec(comm->cell_x0, zones->size[z].x0);
8591 copy_rvec(comm->cell_x1, zones->size[z].x1);
8594 for (d = 0; d < dd->ndim; d++)
8598 for (z = 0; z < zones->n; z++)
8600 /* With a staggered grid we have different sizes
8601 * for non-shifted dimensions.
8603 if (dd->bGridJump && zones->shift[z][dim] == 0)
8607 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8608 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8612 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8613 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8619 rcmbs = comm->cutoff_mbody;
8620 if (ddbox->tric_dir[dim])
8622 rcs /= ddbox->skew_fac[dim];
8623 rcmbs /= ddbox->skew_fac[dim];
8626 /* Set the lower limit for the shifted zone dimensions */
8627 for (z = zone_start; z < zone_end; z++)
8629 if (zones->shift[z][dim] > 0)
8632 if (!dd->bGridJump || d == 0)
8634 zones->size[z].x0[dim] = comm->cell_x1[dim];
8635 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8639 /* Here we take the lower limit of the zone from
8640 * the lowest domain of the zone below.
8644 zones->size[z].x0[dim] =
8645 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8651 zones->size[z].x0[dim] =
8652 zones->size[zone_perm[2][z-4]].x0[dim];
8656 zones->size[z].x0[dim] =
8657 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8660 /* A temporary limit, is updated below */
8661 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8665 for (zi = 0; zi < zones->nizone; zi++)
8667 if (zones->shift[zi][dim] == 0)
8669 /* This takes the whole zone into account.
8670 * With multiple pulses this will lead
8671 * to a larger zone then strictly necessary.
8673 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8674 zones->size[zi].x1[dim]+rcmbs);
8682 /* Loop over the i-zones to set the upper limit of each
8685 for (zi = 0; zi < zones->nizone; zi++)
8687 if (zones->shift[zi][dim] == 0)
8689 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8691 if (zones->shift[z][dim] > 0)
8693 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8694 zones->size[zi].x1[dim]+rcs);
8701 for (z = zone_start; z < zone_end; z++)
8703 /* Initialization only required to keep the compiler happy */
8704 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8707 /* To determine the bounding box for a zone we need to find
8708 * the extreme corners of 4, 2 or 1 corners.
8710 nc = 1 << (ddbox->npbcdim - 1);
8712 for (c = 0; c < nc; c++)
8714 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8718 corner[YY] = zones->size[z].x0[YY];
8722 corner[YY] = zones->size[z].x1[YY];
8726 corner[ZZ] = zones->size[z].x0[ZZ];
8730 corner[ZZ] = zones->size[z].x1[ZZ];
8732 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8734 /* With 1D domain decomposition the cg's are not in
8735 * the triclinic box, but triclinic x-y and rectangular y-z.
8736 * Shift y back, so it will later end up at 0.
8738 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8740 /* Apply the triclinic couplings */
8741 assert(ddbox->npbcdim <= DIM);
8742 for (i = YY; i < ddbox->npbcdim; i++)
8744 for (j = XX; j < i; j++)
8746 corner[j] += corner[i]*box[i][j]/box[i][i];
8751 copy_rvec(corner, corner_min);
8752 copy_rvec(corner, corner_max);
8756 for (i = 0; i < DIM; i++)
8758 corner_min[i] = std::min(corner_min[i], corner[i]);
8759 corner_max[i] = std::max(corner_max[i], corner[i]);
8763 /* Copy the extreme cornes without offset along x */
8764 for (i = 0; i < DIM; i++)
8766 zones->size[z].bb_x0[i] = corner_min[i];
8767 zones->size[z].bb_x1[i] = corner_max[i];
8769 /* Add the offset along x */
8770 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8771 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8774 if (zone_start == 0)
8777 for (dim = 0; dim < DIM; dim++)
8779 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8781 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8786 for (z = zone_start; z < zone_end; z++)
8788 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8790 zones->size[z].x0[XX], zones->size[z].x1[XX],
8791 zones->size[z].x0[YY], zones->size[z].x1[YY],
8792 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8793 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8795 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8796 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8797 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8802 static int comp_cgsort(const void *a, const void *b)
8806 gmx_cgsort_t *cga, *cgb;
8807 cga = (gmx_cgsort_t *)a;
8808 cgb = (gmx_cgsort_t *)b;
8810 comp = cga->nsc - cgb->nsc;
8813 comp = cga->ind_gl - cgb->ind_gl;
8819 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8824 /* Order the data */
8825 for (i = 0; i < n; i++)
8827 buf[i] = a[sort[i].ind];
8830 /* Copy back to the original array */
8831 for (i = 0; i < n; i++)
8837 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8842 /* Order the data */
8843 for (i = 0; i < n; i++)
8845 copy_rvec(v[sort[i].ind], buf[i]);
8848 /* Copy back to the original array */
8849 for (i = 0; i < n; i++)
8851 copy_rvec(buf[i], v[i]);
8855 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8858 int a, atot, cg, cg0, cg1, i;
8860 if (cgindex == NULL)
8862 /* Avoid the useless loop of the atoms within a cg */
8863 order_vec_cg(ncg, sort, v, buf);
8868 /* Order the data */
8870 for (cg = 0; cg < ncg; cg++)
8872 cg0 = cgindex[sort[cg].ind];
8873 cg1 = cgindex[sort[cg].ind+1];
8874 for (i = cg0; i < cg1; i++)
8876 copy_rvec(v[i], buf[a]);
8882 /* Copy back to the original array */
8883 for (a = 0; a < atot; a++)
8885 copy_rvec(buf[a], v[a]);
8889 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8890 int nsort_new, gmx_cgsort_t *sort_new,
8891 gmx_cgsort_t *sort1)
8895 /* The new indices are not very ordered, so we qsort them */
8896 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8898 /* sort2 is already ordered, so now we can merge the two arrays */
8902 while (i2 < nsort2 || i_new < nsort_new)
8906 sort1[i1++] = sort_new[i_new++];
8908 else if (i_new == nsort_new)
8910 sort1[i1++] = sort2[i2++];
8912 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8913 (sort2[i2].nsc == sort_new[i_new].nsc &&
8914 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8916 sort1[i1++] = sort2[i2++];
8920 sort1[i1++] = sort_new[i_new++];
8925 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8927 gmx_domdec_sort_t *sort;
8928 gmx_cgsort_t *cgsort, *sort_i;
8929 int ncg_new, nsort2, nsort_new, i, *a, moved;
8931 sort = dd->comm->sort;
8933 a = fr->ns.grid->cell_index;
8935 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8937 if (ncg_home_old >= 0)
8939 /* The charge groups that remained in the same ns grid cell
8940 * are completely ordered. So we can sort efficiently by sorting
8941 * the charge groups that did move into the stationary list.
8946 for (i = 0; i < dd->ncg_home; i++)
8948 /* Check if this cg did not move to another node */
8951 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8953 /* This cg is new on this node or moved ns grid cell */
8954 if (nsort_new >= sort->sort_new_nalloc)
8956 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8957 srenew(sort->sort_new, sort->sort_new_nalloc);
8959 sort_i = &(sort->sort_new[nsort_new++]);
8963 /* This cg did not move */
8964 sort_i = &(sort->sort2[nsort2++]);
8966 /* Sort on the ns grid cell indices
8967 * and the global topology index.
8968 * index_gl is irrelevant with cell ns,
8969 * but we set it here anyhow to avoid a conditional.
8972 sort_i->ind_gl = dd->index_gl[i];
8979 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8982 /* Sort efficiently */
8983 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8988 cgsort = sort->sort;
8990 for (i = 0; i < dd->ncg_home; i++)
8992 /* Sort on the ns grid cell indices
8993 * and the global topology index
8995 cgsort[i].nsc = a[i];
8996 cgsort[i].ind_gl = dd->index_gl[i];
8998 if (cgsort[i].nsc < moved)
9005 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9007 /* Determine the order of the charge groups using qsort */
9008 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9014 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9017 int ncg_new, i, *a, na;
9019 sort = dd->comm->sort->sort;
9021 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9024 for (i = 0; i < na; i++)
9028 sort[ncg_new].ind = a[i];
9036 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9039 gmx_domdec_sort_t *sort;
9040 gmx_cgsort_t *cgsort;
9042 int ncg_new, i, *ibuf, cgsize;
9045 sort = dd->comm->sort;
9047 if (dd->ncg_home > sort->sort_nalloc)
9049 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9050 srenew(sort->sort, sort->sort_nalloc);
9051 srenew(sort->sort2, sort->sort_nalloc);
9053 cgsort = sort->sort;
9055 switch (fr->cutoff_scheme)
9058 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9061 ncg_new = dd_sort_order_nbnxn(dd, fr);
9064 gmx_incons("unimplemented");
9068 /* We alloc with the old size, since cgindex is still old */
9069 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9070 vbuf = dd->comm->vbuf.v;
9074 cgindex = dd->cgindex;
9081 /* Remove the charge groups which are no longer at home here */
9082 dd->ncg_home = ncg_new;
9085 fprintf(debug, "Set the new home charge group count to %d\n",
9089 /* Reorder the state */
9090 for (i = 0; i < estNR; i++)
9092 if (EST_DISTR(i) && (state->flags & (1<<i)))
9097 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9100 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9103 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9106 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9110 case estDISRE_INITF:
9111 case estDISRE_RM3TAV:
9112 case estORIRE_INITF:
9114 /* No ordering required */
9117 gmx_incons("Unknown state entry encountered in dd_sort_state");
9122 if (fr->cutoff_scheme == ecutsGROUP)
9125 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9128 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9130 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9131 srenew(sort->ibuf, sort->ibuf_nalloc);
9134 /* Reorder the global cg index */
9135 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9136 /* Reorder the cginfo */
9137 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9138 /* Rebuild the local cg index */
9142 for (i = 0; i < dd->ncg_home; i++)
9144 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9145 ibuf[i+1] = ibuf[i] + cgsize;
9147 for (i = 0; i < dd->ncg_home+1; i++)
9149 dd->cgindex[i] = ibuf[i];
9154 for (i = 0; i < dd->ncg_home+1; i++)
9159 /* Set the home atom number */
9160 dd->nat_home = dd->cgindex[dd->ncg_home];
9162 if (fr->cutoff_scheme == ecutsVERLET)
9164 /* The atoms are now exactly in grid order, update the grid order */
9165 nbnxn_set_atomorder(fr->nbv->nbs);
9169 /* Copy the sorted ns cell indices back to the ns grid struct */
9170 for (i = 0; i < dd->ncg_home; i++)
9172 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9174 fr->ns.grid->nr = dd->ncg_home;
9178 static void add_dd_statistics(gmx_domdec_t *dd)
9180 gmx_domdec_comm_t *comm;
9185 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9187 comm->sum_nat[ddnat-ddnatZONE] +=
9188 comm->nat[ddnat] - comm->nat[ddnat-1];
9193 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9195 gmx_domdec_comm_t *comm;
9200 /* Reset all the statistics and counters for total run counting */
9201 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9203 comm->sum_nat[ddnat-ddnatZONE] = 0;
9207 comm->load_step = 0;
9210 clear_ivec(comm->load_lim);
9215 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9217 gmx_domdec_comm_t *comm;
9221 comm = cr->dd->comm;
9223 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9230 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");
9232 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9234 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9239 " av. #atoms communicated per step for force: %d x %.1f\n",
9243 if (cr->dd->vsite_comm)
9246 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9247 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9252 if (cr->dd->constraint_comm)
9255 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9256 1 + ir->nLincsIter, av);
9260 gmx_incons(" Unknown type for DD statistics");
9263 fprintf(fplog, "\n");
9265 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9267 print_dd_load_av(fplog, cr->dd);
9271 void dd_partition_system(FILE *fplog,
9274 gmx_bool bMasterState,
9276 t_state *state_global,
9277 gmx_mtop_t *top_global,
9279 t_state *state_local,
9282 gmx_localtop_t *top_local,
9285 gmx_shellfc_t shellfc,
9286 gmx_constr_t constr,
9288 gmx_wallcycle_t wcycle,
9292 gmx_domdec_comm_t *comm;
9293 gmx_ddbox_t ddbox = {0};
9295 gmx_int64_t step_pcoupl;
9296 rvec cell_ns_x0, cell_ns_x1;
9297 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9298 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9299 gmx_bool bRedist, bSortCG, bResortAll;
9300 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9307 bBoxChanged = (bMasterState || DEFORM(*ir));
9308 if (ir->epc != epcNO)
9310 /* With nstpcouple > 1 pressure coupling happens.
9311 * one step after calculating the pressure.
9312 * Box scaling happens at the end of the MD step,
9313 * after the DD partitioning.
9314 * We therefore have to do DLB in the first partitioning
9315 * after an MD step where P-coupling occured.
9316 * We need to determine the last step in which p-coupling occurred.
9317 * MRS -- need to validate this for vv?
9322 step_pcoupl = step - 1;
9326 step_pcoupl = ((step - 1)/n)*n + 1;
9328 if (step_pcoupl >= comm->partition_step)
9334 bNStGlobalComm = (step % nstglobalcomm == 0);
9336 if (!comm->bDynLoadBal)
9342 /* Should we do dynamic load balacing this step?
9343 * Since it requires (possibly expensive) global communication,
9344 * we might want to do DLB less frequently.
9346 if (bBoxChanged || ir->epc != epcNO)
9348 bDoDLB = bBoxChanged;
9352 bDoDLB = bNStGlobalComm;
9356 /* Check if we have recorded loads on the nodes */
9357 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9359 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal && !dd_dlb_is_locked(dd))
9361 /* Check if we should use DLB at the second partitioning
9362 * and every 100 partitionings,
9363 * so the extra communication cost is negligible.
9365 const int nddp_chk_dlb = 100;
9366 bCheckDLB = (comm->n_load_collect == 0 ||
9367 comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
9374 /* Print load every nstlog, first and last step to the log file */
9375 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9376 comm->n_load_collect == 0 ||
9378 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9380 /* Avoid extra communication due to verbose screen output
9381 * when nstglobalcomm is set.
9383 if (bDoDLB || bLogLoad || bCheckDLB ||
9384 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9386 get_load_distribution(dd, wcycle);
9391 dd_print_load(fplog, dd, step-1);
9395 dd_print_load_verbose(dd);
9398 comm->n_load_collect++;
9402 /* Since the timings are node dependent, the master decides */
9405 /* Here we check if the max PME rank load is more than 0.98
9406 * the max PP force load. If so, PP DLB will not help,
9407 * since we are (almost) limited by PME. Furthermore,
9408 * DLB will cause a significant extra x/f redistribution
9409 * cost on the PME ranks, which will then surely result
9410 * in lower total performance.
9411 * This check might be fragile, since one measurement
9412 * below 0.98 (although only done once every 100 DD part.)
9413 * could turn on DLB for the rest of the run.
9415 if (cr->npmenodes > 0 &&
9416 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9423 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9427 fprintf(debug, "step %s, imb loss %f\n",
9428 gmx_step_str(step, sbuf),
9429 dd_force_imb_perf_loss(dd));
9432 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9435 turn_on_dlb(fplog, cr, step);
9440 comm->n_load_have++;
9443 cgs_gl = &comm->cgs_gl;
9448 /* Clear the old state */
9449 clear_dd_indices(dd, 0, 0);
9452 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9453 TRUE, cgs_gl, state_global->x, &ddbox);
9455 get_cg_distribution(fplog, step, dd, cgs_gl,
9456 state_global->box, &ddbox, state_global->x);
9458 dd_distribute_state(dd, cgs_gl,
9459 state_global, state_local, f);
9461 dd_make_local_cgs(dd, &top_local->cgs);
9463 /* Ensure that we have space for the new distribution */
9464 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9466 if (fr->cutoff_scheme == ecutsGROUP)
9468 calc_cgcm(fplog, 0, dd->ncg_home,
9469 &top_local->cgs, state_local->x, fr->cg_cm);
9472 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9474 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9476 else if (state_local->ddp_count != dd->ddp_count)
9478 if (state_local->ddp_count > dd->ddp_count)
9480 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9483 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9485 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);
9488 /* Clear the old state */
9489 clear_dd_indices(dd, 0, 0);
9491 /* Build the new indices */
9492 rebuild_cgindex(dd, cgs_gl->index, state_local);
9493 make_dd_indices(dd, cgs_gl->index, 0);
9494 ncgindex_set = dd->ncg_home;
9496 if (fr->cutoff_scheme == ecutsGROUP)
9498 /* Redetermine the cg COMs */
9499 calc_cgcm(fplog, 0, dd->ncg_home,
9500 &top_local->cgs, state_local->x, fr->cg_cm);
9503 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9505 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9507 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9508 TRUE, &top_local->cgs, state_local->x, &ddbox);
9510 bRedist = comm->bDynLoadBal;
9514 /* We have the full state, only redistribute the cgs */
9516 /* Clear the non-home indices */
9517 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9520 /* Avoid global communication for dim's without pbc and -gcom */
9521 if (!bNStGlobalComm)
9523 copy_rvec(comm->box0, ddbox.box0 );
9524 copy_rvec(comm->box_size, ddbox.box_size);
9526 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9527 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9532 /* For dim's without pbc and -gcom */
9533 copy_rvec(ddbox.box0, comm->box0 );
9534 copy_rvec(ddbox.box_size, comm->box_size);
9536 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9539 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9541 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9544 /* Check if we should sort the charge groups */
9545 if (comm->nstSortCG > 0)
9547 bSortCG = (bMasterState ||
9548 (bRedist && (step % comm->nstSortCG == 0)));
9555 ncg_home_old = dd->ncg_home;
9560 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9562 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9564 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9566 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9569 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9571 &comm->cell_x0, &comm->cell_x1,
9572 dd->ncg_home, fr->cg_cm,
9573 cell_ns_x0, cell_ns_x1, &grid_density);
9577 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9580 switch (fr->cutoff_scheme)
9583 copy_ivec(fr->ns.grid->n, ncells_old);
9584 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9585 state_local->box, cell_ns_x0, cell_ns_x1,
9586 fr->rlistlong, grid_density);
9589 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9592 gmx_incons("unimplemented");
9594 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9595 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9599 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9601 /* Sort the state on charge group position.
9602 * This enables exact restarts from this step.
9603 * It also improves performance by about 15% with larger numbers
9604 * of atoms per node.
9607 /* Fill the ns grid with the home cell,
9608 * so we can sort with the indices.
9610 set_zones_ncg_home(dd);
9612 switch (fr->cutoff_scheme)
9615 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9617 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9619 comm->zones.size[0].bb_x0,
9620 comm->zones.size[0].bb_x1,
9622 comm->zones.dens_zone0,
9625 ncg_moved, bRedist ? comm->moved : NULL,
9626 fr->nbv->grp[eintLocal].kernel_type,
9627 fr->nbv->grp[eintLocal].nbat);
9629 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9632 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9633 0, dd->ncg_home, fr->cg_cm);
9635 copy_ivec(fr->ns.grid->n, ncells_new);
9638 gmx_incons("unimplemented");
9641 bResortAll = bMasterState;
9643 /* Check if we can user the old order and ns grid cell indices
9644 * of the charge groups to sort the charge groups efficiently.
9646 if (ncells_new[XX] != ncells_old[XX] ||
9647 ncells_new[YY] != ncells_old[YY] ||
9648 ncells_new[ZZ] != ncells_old[ZZ])
9655 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9656 gmx_step_str(step, sbuf), dd->ncg_home);
9658 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9659 bResortAll ? -1 : ncg_home_old);
9660 /* Rebuild all the indices */
9661 ga2la_clear(dd->ga2la);
9664 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9667 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9669 /* Setup up the communication and communicate the coordinates */
9670 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9672 /* Set the indices */
9673 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9675 /* Set the charge group boundaries for neighbor searching */
9676 set_cg_boundaries(&comm->zones);
9678 if (fr->cutoff_scheme == ecutsVERLET)
9680 set_zones_size(dd, state_local->box, &ddbox,
9681 bSortCG ? 1 : 0, comm->zones.n);
9684 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9687 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9688 -1,state_local->x,state_local->box);
9691 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9693 /* Extract a local topology from the global topology */
9694 for (i = 0; i < dd->ndim; i++)
9696 np[dd->dim[i]] = comm->cd[i].np;
9698 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9699 comm->cellsize_min, np,
9701 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9702 vsite, top_global, top_local);
9704 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9706 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9708 /* Set up the special atom communication */
9709 n = comm->nat[ddnatZONE];
9710 for (i = ddnatZONE+1; i < ddnatNR; i++)
9715 if (vsite && vsite->n_intercg_vsite)
9717 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9721 if (dd->bInterCGcons || dd->bInterCGsettles)
9723 /* Only for inter-cg constraints we need special code */
9724 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9725 constr, ir->nProjOrder,
9726 top_local->idef.il);
9730 gmx_incons("Unknown special atom type setup");
9735 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9737 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9739 /* Make space for the extra coordinates for virtual site
9740 * or constraint communication.
9742 state_local->natoms = comm->nat[ddnatNR-1];
9743 if (state_local->natoms > state_local->nalloc)
9745 dd_realloc_state(state_local, f, state_local->natoms);
9748 if (fr->bF_NoVirSum)
9750 if (vsite && vsite->n_intercg_vsite)
9752 nat_f_novirsum = comm->nat[ddnatVSITE];
9756 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9758 nat_f_novirsum = dd->nat_tot;
9762 nat_f_novirsum = dd->nat_home;
9771 /* Set the number of atoms required for the force calculation.
9772 * Forces need to be constrained when using a twin-range setup
9773 * or with energy minimization. For simple simulations we could
9774 * avoid some allocation, zeroing and copying, but this is
9775 * probably not worth the complications ande checking.
9777 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9778 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9780 /* We make the all mdatoms up to nat_tot_con.
9781 * We could save some work by only setting invmass
9782 * between nat_tot and nat_tot_con.
9784 /* This call also sets the new number of home particles to dd->nat_home */
9785 atoms2md(top_global, ir,
9786 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9788 /* Now we have the charges we can sort the FE interactions */
9789 dd_sort_local_top(dd, mdatoms, top_local);
9793 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9794 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9795 mdatoms, FALSE, vsite);
9800 /* Make the local shell stuff, currently no communication is done */
9801 make_local_shells(cr, mdatoms, shellfc);
9804 if (ir->implicit_solvent)
9806 make_local_gb(cr, fr->born, ir->gb_algorithm);
9809 setup_bonded_threading(fr, &top_local->idef);
9811 if (!(cr->duty & DUTY_PME))
9813 /* Send the charges and/or c6/sigmas to our PME only node */
9814 gmx_pme_send_parameters(cr,
9816 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9817 mdatoms->chargeA, mdatoms->chargeB,
9818 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9819 mdatoms->sigmaA, mdatoms->sigmaB,
9820 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9825 set_constraints(constr, top_local, ir, mdatoms, cr);
9828 if (ir->ePull != epullNO)
9830 /* Update the local pull groups */
9831 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9836 /* Update the local rotation groups */
9837 dd_make_local_rotation_groups(dd, ir->rot);
9840 if (ir->eSwapCoords != eswapNO)
9842 /* Update the local groups needed for ion swapping */
9843 dd_make_local_swap_groups(dd, ir->swap);
9846 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9847 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9849 add_dd_statistics(dd);
9851 /* Make sure we only count the cycles for this DD partitioning */
9852 clear_dd_cycle_counts(dd);
9854 /* Because the order of the atoms might have changed since
9855 * the last vsite construction, we need to communicate the constructing
9856 * atom coordinates again (for spreading the forces this MD step).
9858 dd_move_x_vsites(dd, state_local->box, state_local->x);
9860 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9862 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9864 dd_move_x(dd, state_local->box, state_local->x);
9865 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9866 -1, state_local->x, state_local->box);
9869 /* Store the partitioning step */
9870 comm->partition_step = step;
9872 /* Increase the DD partitioning counter */
9874 /* The state currently matches this DD partitioning count, store it */
9875 state_local->ddp_count = dd->ddp_count;
9878 /* The DD master node knows the complete cg distribution,
9879 * store the count so we can possibly skip the cg info communication.
9881 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9884 if (comm->DD_debug > 0)
9886 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9887 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9888 "after partitioning");