1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
30 #include "gmx_fatal.h"
31 #include "gmx_fatal_collective.h"
34 #include "domdec_network.h"
37 #include "chargegroup.h"
46 #include "pull_rotation.h"
47 #include "gmx_wallcycle.h"
51 #include "mtop_util.h"
53 #include "gmx_ga2la.h"
56 #include "nbnxn_search.h"
58 #include "gmx_omp_nthreads.h"
60 #include "gromacs/utility/gmxmpi.h"
62 #define DDRANK(dd, rank) (rank)
63 #define DDMASTERRANK(dd) (dd->masterrank)
65 typedef struct gmx_domdec_master
67 /* The cell boundaries */
69 /* The global charge group division */
70 int *ncg; /* Number of home charge groups for each node */
71 int *index; /* Index of nnodes+1 into cg */
72 int *cg; /* Global charge group index */
73 int *nat; /* Number of home atoms for each node. */
74 int *ibuf; /* Buffer for communication */
75 rvec *vbuf; /* Buffer for state scattering and gathering */
76 } gmx_domdec_master_t;
80 /* The numbers of charge groups to send and receive for each cell
81 * that requires communication, the last entry contains the total
82 * number of atoms that needs to be communicated.
84 int nsend[DD_MAXIZONE+2];
85 int nrecv[DD_MAXIZONE+2];
86 /* The charge groups to send */
89 /* The atom range for non-in-place communication */
90 int cell2at0[DD_MAXIZONE];
91 int cell2at1[DD_MAXIZONE];
96 int np; /* Number of grid pulses in this dimension */
97 int np_dlb; /* For dlb, for use with edlbAUTO */
98 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
100 gmx_bool bInPlace; /* Can we communicate in place? */
101 } gmx_domdec_comm_dim_t;
105 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
106 real *cell_f; /* State var.: cell boundaries, box relative */
107 real *old_cell_f; /* Temp. var.: old cell size */
108 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
109 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
110 real *bound_min; /* Temp. var.: lower limit for cell boundary */
111 real *bound_max; /* Temp. var.: upper limit for cell boundary */
112 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
113 real *buf_ncd; /* Temp. var. */
116 #define DD_NLOAD_MAX 9
118 /* Here floats are accurate enough, since these variables
119 * only influence the load balancing, not the actual MD results.
146 gmx_cgsort_t *sort_new;
158 /* This enum determines the order of the coordinates.
159 * ddnatHOME and ddnatZONE should be first and second,
160 * the others can be ordered as wanted.
163 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
167 edlbAUTO, edlbNO, edlbYES, edlbNR
169 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
173 int dim; /* The dimension */
174 gmx_bool dim_match; /* Tells if DD and PME dims match */
175 int nslab; /* The number of PME slabs in this dimension */
176 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
177 int *pp_min; /* The minimum pp node location, size nslab */
178 int *pp_max; /* The maximum pp node location,size nslab */
179 int maxshift; /* The maximum shift for coordinate redistribution in PME */
184 real min0; /* The minimum bottom of this zone */
185 real max1; /* The maximum top of this zone */
186 real min1; /* The minimum top of this zone */
187 real mch0; /* The maximum bottom communicaton height for this zone */
188 real mch1; /* The maximum top communicaton height for this zone */
189 real p1_0; /* The bottom value of the first cell in this zone */
190 real p1_1; /* The top value of the first cell in this zone */
195 gmx_domdec_ind_t ind;
202 } dd_comm_setup_work_t;
204 typedef struct gmx_domdec_comm
206 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
207 * unless stated otherwise.
210 /* The number of decomposition dimensions for PME, 0: no PME */
212 /* The number of nodes doing PME (PP/PME or only PME) */
216 /* The communication setup including the PME only nodes */
217 gmx_bool bCartesianPP_PME;
220 int *pmenodes; /* size npmenodes */
221 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
222 * but with bCartesianPP_PME */
223 gmx_ddpme_t ddpme[2];
225 /* The DD particle-particle nodes only */
226 gmx_bool bCartesianPP;
227 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
229 /* The global charge groups */
232 /* Should we sort the cgs */
234 gmx_domdec_sort_t *sort;
236 /* Are there charge groups? */
239 /* Are there bonded and multi-body interactions between charge groups? */
240 gmx_bool bInterCGBondeds;
241 gmx_bool bInterCGMultiBody;
243 /* Data for the optional bonded interaction atom communication range */
250 /* Are we actually using DLB? */
251 gmx_bool bDynLoadBal;
253 /* Cell sizes for static load balancing, first index cartesian */
256 /* The width of the communicated boundaries */
259 /* The minimum cell size (including triclinic correction) */
261 /* For dlb, for use with edlbAUTO */
262 rvec cellsize_min_dlb;
263 /* The lower limit for the DD cell size with DLB */
265 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
266 gmx_bool bVacDLBNoLimit;
268 /* With PME load balancing we set limits on DLB */
269 gmx_bool bPMELoadBalDLBLimits;
270 /* DLB needs to take into account that we want to allow this maximum
271 * cut-off (for PME load balancing), this could limit cell boundaries.
273 real PMELoadBal_max_cutoff;
275 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
277 /* box0 and box_size are required with dim's without pbc and -gcom */
281 /* The cell boundaries */
285 /* The old location of the cell boundaries, to check cg displacements */
289 /* The communication setup and charge group boundaries for the zones */
290 gmx_domdec_zones_t zones;
292 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
293 * cell boundaries of neighboring cells for dynamic load balancing.
295 gmx_ddzone_t zone_d1[2];
296 gmx_ddzone_t zone_d2[2][2];
298 /* The coordinate/force communication setup and indices */
299 gmx_domdec_comm_dim_t cd[DIM];
300 /* The maximum number of cells to communicate with in one dimension */
303 /* Which cg distribution is stored on the master node */
304 int master_cg_ddp_count;
306 /* The number of cg's received from the direct neighbors */
307 int zone_ncg1[DD_MAXZONE];
309 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
312 /* Array for signalling if atoms have moved to another domain */
316 /* Communication buffer for general use */
320 /* Communication buffer for general use */
323 /* Temporary storage for thread parallel communication setup */
325 dd_comm_setup_work_t *dth;
327 /* Communication buffers only used with multiple grid pulses */
332 /* Communication buffers for local redistribution */
334 int cggl_flag_nalloc[DIM*2];
336 int cgcm_state_nalloc[DIM*2];
338 /* Cell sizes for dynamic load balancing */
339 gmx_domdec_root_t **root;
343 real cell_f_max0[DIM];
344 real cell_f_min1[DIM];
346 /* Stuff for load communication */
347 gmx_bool bRecordLoad;
348 gmx_domdec_load_t *load;
350 MPI_Comm *mpi_comm_load;
353 /* Maximum DLB scaling per load balancing step in percent */
357 float cycl[ddCyclNr];
358 int cycl_n[ddCyclNr];
359 float cycl_max[ddCyclNr];
360 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
364 /* Have often have did we have load measurements */
366 /* Have often have we collected the load measurements */
370 double sum_nat[ddnatNR-ddnatZONE];
380 /* The last partition step */
381 gmx_large_int_t partition_step;
389 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
392 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
393 #define DD_FLAG_NRCG 65535
394 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
395 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
397 /* Zone permutation required to obtain consecutive charge groups
398 * for neighbor searching.
400 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
402 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
403 * components see only j zones with that component 0.
406 /* The DD zone order */
407 static const ivec dd_zo[DD_MAXZONE] =
408 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
413 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
418 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
423 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
425 /* Factors used to avoid problems due to rounding issues */
426 #define DD_CELL_MARGIN 1.0001
427 #define DD_CELL_MARGIN2 1.00005
428 /* Factor to account for pressure scaling during nstlist steps */
429 #define DD_PRES_SCALE_MARGIN 1.02
431 /* Allowed performance loss before we DLB or warn */
432 #define DD_PERF_LOSS 0.05
434 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
436 /* Use separate MPI send and receive commands
437 * when nnodes <= GMX_DD_NNODES_SENDRECV.
438 * This saves memory (and some copying for small nnodes).
439 * For high parallelization scatter and gather calls are used.
441 #define GMX_DD_NNODES_SENDRECV 4
445 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
447 static void index2xyz(ivec nc,int ind,ivec xyz)
449 xyz[XX] = ind % nc[XX];
450 xyz[YY] = (ind / nc[XX]) % nc[YY];
451 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
455 /* This order is required to minimize the coordinate communication in PME
456 * which uses decomposition in the x direction.
458 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
460 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
462 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
463 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
464 xyz[ZZ] = ind % nc[ZZ];
467 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
472 ddindex = dd_index(dd->nc, c);
473 if (dd->comm->bCartesianPP_PME)
475 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
477 else if (dd->comm->bCartesianPP)
480 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
491 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
493 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
496 int ddglatnr(gmx_domdec_t *dd, int i)
506 if (i >= dd->comm->nat[ddnatNR-1])
508 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
510 atnr = dd->gatindex[i] + 1;
516 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
518 return &dd->comm->cgs_gl;
521 static void vec_rvec_init(vec_rvec_t *v)
527 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
531 v->nalloc = over_alloc_dd(n);
532 srenew(v->v, v->nalloc);
536 void dd_store_state(gmx_domdec_t *dd, t_state *state)
540 if (state->ddp_count != dd->ddp_count)
542 gmx_incons("The state does not the domain decomposition state");
545 state->ncg_gl = dd->ncg_home;
546 if (state->ncg_gl > state->cg_gl_nalloc)
548 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
549 srenew(state->cg_gl, state->cg_gl_nalloc);
551 for (i = 0; i < state->ncg_gl; i++)
553 state->cg_gl[i] = dd->index_gl[i];
556 state->ddp_count_cg_gl = dd->ddp_count;
559 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
561 return &dd->comm->zones;
564 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
565 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
567 gmx_domdec_zones_t *zones;
570 zones = &dd->comm->zones;
573 while (icg >= zones->izone[izone].cg1)
582 else if (izone < zones->nizone)
584 *jcg0 = zones->izone[izone].jcg0;
588 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
589 icg, izone, zones->nizone);
592 *jcg1 = zones->izone[izone].jcg1;
594 for (d = 0; d < dd->ndim; d++)
597 shift0[dim] = zones->izone[izone].shift0[dim];
598 shift1[dim] = zones->izone[izone].shift1[dim];
599 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
601 /* A conservative approach, this can be optimized */
608 int dd_natoms_vsite(gmx_domdec_t *dd)
610 return dd->comm->nat[ddnatVSITE];
613 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
615 *at_start = dd->comm->nat[ddnatCON-1];
616 *at_end = dd->comm->nat[ddnatCON];
619 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
621 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
622 int *index, *cgindex;
623 gmx_domdec_comm_t *comm;
624 gmx_domdec_comm_dim_t *cd;
625 gmx_domdec_ind_t *ind;
626 rvec shift = {0, 0, 0}, *buf, *rbuf;
627 gmx_bool bPBC, bScrew;
631 cgindex = dd->cgindex;
636 nat_tot = dd->nat_home;
637 for (d = 0; d < dd->ndim; d++)
639 bPBC = (dd->ci[dd->dim[d]] == 0);
640 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
643 copy_rvec(box[dd->dim[d]], shift);
646 for (p = 0; p < cd->np; p++)
653 for (i = 0; i < ind->nsend[nzone]; i++)
655 at0 = cgindex[index[i]];
656 at1 = cgindex[index[i]+1];
657 for (j = at0; j < at1; j++)
659 copy_rvec(x[j], buf[n]);
666 for (i = 0; i < ind->nsend[nzone]; i++)
668 at0 = cgindex[index[i]];
669 at1 = cgindex[index[i]+1];
670 for (j = at0; j < at1; j++)
672 /* We need to shift the coordinates */
673 rvec_add(x[j], shift, buf[n]);
680 for (i = 0; i < ind->nsend[nzone]; i++)
682 at0 = cgindex[index[i]];
683 at1 = cgindex[index[i]+1];
684 for (j = at0; j < at1; j++)
687 buf[n][XX] = x[j][XX] + shift[XX];
689 * This operation requires a special shift force
690 * treatment, which is performed in calc_vir.
692 buf[n][YY] = box[YY][YY] - x[j][YY];
693 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
705 rbuf = comm->vbuf2.v;
707 /* Send and receive the coordinates */
708 dd_sendrecv_rvec(dd, d, dddirBackward,
709 buf, ind->nsend[nzone+1],
710 rbuf, ind->nrecv[nzone+1]);
714 for (zone = 0; zone < nzone; zone++)
716 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
718 copy_rvec(rbuf[j], x[i]);
723 nat_tot += ind->nrecv[nzone+1];
729 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
731 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
732 int *index, *cgindex;
733 gmx_domdec_comm_t *comm;
734 gmx_domdec_comm_dim_t *cd;
735 gmx_domdec_ind_t *ind;
739 gmx_bool bPBC, bScrew;
743 cgindex = dd->cgindex;
748 nzone = comm->zones.n/2;
749 nat_tot = dd->nat_tot;
750 for (d = dd->ndim-1; d >= 0; d--)
752 bPBC = (dd->ci[dd->dim[d]] == 0);
753 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
754 if (fshift == NULL && !bScrew)
758 /* Determine which shift vector we need */
764 for (p = cd->np-1; p >= 0; p--)
767 nat_tot -= ind->nrecv[nzone+1];
774 sbuf = comm->vbuf2.v;
776 for (zone = 0; zone < nzone; zone++)
778 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
780 copy_rvec(f[i], sbuf[j]);
785 /* Communicate the forces */
786 dd_sendrecv_rvec(dd, d, dddirForward,
787 sbuf, ind->nrecv[nzone+1],
788 buf, ind->nsend[nzone+1]);
790 /* Add the received forces */
794 for (i = 0; i < ind->nsend[nzone]; i++)
796 at0 = cgindex[index[i]];
797 at1 = cgindex[index[i]+1];
798 for (j = at0; j < at1; j++)
800 rvec_inc(f[j], buf[n]);
807 for (i = 0; i < ind->nsend[nzone]; i++)
809 at0 = cgindex[index[i]];
810 at1 = cgindex[index[i]+1];
811 for (j = at0; j < at1; j++)
813 rvec_inc(f[j], buf[n]);
814 /* Add this force to the shift force */
815 rvec_inc(fshift[is], buf[n]);
822 for (i = 0; i < ind->nsend[nzone]; i++)
824 at0 = cgindex[index[i]];
825 at1 = cgindex[index[i]+1];
826 for (j = at0; j < at1; j++)
828 /* Rotate the force */
829 f[j][XX] += buf[n][XX];
830 f[j][YY] -= buf[n][YY];
831 f[j][ZZ] -= buf[n][ZZ];
834 /* Add this force to the shift force */
835 rvec_inc(fshift[is], buf[n]);
846 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
848 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
849 int *index, *cgindex;
850 gmx_domdec_comm_t *comm;
851 gmx_domdec_comm_dim_t *cd;
852 gmx_domdec_ind_t *ind;
857 cgindex = dd->cgindex;
859 buf = &comm->vbuf.v[0][0];
862 nat_tot = dd->nat_home;
863 for (d = 0; d < dd->ndim; d++)
866 for (p = 0; p < cd->np; p++)
871 for (i = 0; i < ind->nsend[nzone]; i++)
873 at0 = cgindex[index[i]];
874 at1 = cgindex[index[i]+1];
875 for (j = at0; j < at1; j++)
888 rbuf = &comm->vbuf2.v[0][0];
890 /* Send and receive the coordinates */
891 dd_sendrecv_real(dd, d, dddirBackward,
892 buf, ind->nsend[nzone+1],
893 rbuf, ind->nrecv[nzone+1]);
897 for (zone = 0; zone < nzone; zone++)
899 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
906 nat_tot += ind->nrecv[nzone+1];
912 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
914 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
915 int *index, *cgindex;
916 gmx_domdec_comm_t *comm;
917 gmx_domdec_comm_dim_t *cd;
918 gmx_domdec_ind_t *ind;
923 cgindex = dd->cgindex;
925 buf = &comm->vbuf.v[0][0];
928 nzone = comm->zones.n/2;
929 nat_tot = dd->nat_tot;
930 for (d = dd->ndim-1; d >= 0; d--)
933 for (p = cd->np-1; p >= 0; p--)
936 nat_tot -= ind->nrecv[nzone+1];
943 sbuf = &comm->vbuf2.v[0][0];
945 for (zone = 0; zone < nzone; zone++)
947 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
954 /* Communicate the forces */
955 dd_sendrecv_real(dd, d, dddirForward,
956 sbuf, ind->nrecv[nzone+1],
957 buf, ind->nsend[nzone+1]);
959 /* Add the received forces */
961 for (i = 0; i < ind->nsend[nzone]; i++)
963 at0 = cgindex[index[i]];
964 at1 = cgindex[index[i]+1];
965 for (j = at0; j < at1; j++)
976 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
978 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",
980 zone->min0, zone->max1,
981 zone->mch0, zone->mch0,
982 zone->p1_0, zone->p1_1);
986 #define DDZONECOMM_MAXZONE 5
987 #define DDZONECOMM_BUFSIZE 3
989 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
990 int ddimind, int direction,
991 gmx_ddzone_t *buf_s, int n_s,
992 gmx_ddzone_t *buf_r, int n_r)
994 #define ZBS DDZONECOMM_BUFSIZE
995 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
996 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
999 for (i = 0; i < n_s; i++)
1001 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1002 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1003 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1004 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1005 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1006 vbuf_s[i*ZBS+1][2] = 0;
1007 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1008 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1009 vbuf_s[i*ZBS+2][2] = 0;
1012 dd_sendrecv_rvec(dd, ddimind, direction,
1016 for (i = 0; i < n_r; i++)
1018 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1019 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1020 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1021 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1022 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1023 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1024 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1030 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1031 rvec cell_ns_x0, rvec cell_ns_x1)
1033 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1035 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1036 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1037 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1038 rvec extr_s[2], extr_r[2];
1040 real dist_d, c = 0, det;
1041 gmx_domdec_comm_t *comm;
1042 gmx_bool bPBC, bUse;
1046 for (d = 1; d < dd->ndim; d++)
1049 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1050 zp->min0 = cell_ns_x0[dim];
1051 zp->max1 = cell_ns_x1[dim];
1052 zp->min1 = cell_ns_x1[dim];
1053 zp->mch0 = cell_ns_x0[dim];
1054 zp->mch1 = cell_ns_x1[dim];
1055 zp->p1_0 = cell_ns_x0[dim];
1056 zp->p1_1 = cell_ns_x1[dim];
1059 for (d = dd->ndim-2; d >= 0; d--)
1062 bPBC = (dim < ddbox->npbcdim);
1064 /* Use an rvec to store two reals */
1065 extr_s[d][0] = comm->cell_f0[d+1];
1066 extr_s[d][1] = comm->cell_f1[d+1];
1067 extr_s[d][2] = comm->cell_f1[d+1];
1070 /* Store the extremes in the backward sending buffer,
1071 * so the get updated separately from the forward communication.
1073 for (d1 = d; d1 < dd->ndim-1; d1++)
1075 /* We invert the order to be able to use the same loop for buf_e */
1076 buf_s[pos].min0 = extr_s[d1][1];
1077 buf_s[pos].max1 = extr_s[d1][0];
1078 buf_s[pos].min1 = extr_s[d1][2];
1079 buf_s[pos].mch0 = 0;
1080 buf_s[pos].mch1 = 0;
1081 /* Store the cell corner of the dimension we communicate along */
1082 buf_s[pos].p1_0 = comm->cell_x0[dim];
1083 buf_s[pos].p1_1 = 0;
1087 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1090 if (dd->ndim == 3 && d == 0)
1092 buf_s[pos] = comm->zone_d2[0][1];
1094 buf_s[pos] = comm->zone_d1[0];
1098 /* We only need to communicate the extremes
1099 * in the forward direction
1101 npulse = comm->cd[d].np;
1104 /* Take the minimum to avoid double communication */
1105 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1109 /* Without PBC we should really not communicate over
1110 * the boundaries, but implementing that complicates
1111 * the communication setup and therefore we simply
1112 * do all communication, but ignore some data.
1114 npulse_min = npulse;
1116 for (p = 0; p < npulse_min; p++)
1118 /* Communicate the extremes forward */
1119 bUse = (bPBC || dd->ci[dim] > 0);
1121 dd_sendrecv_rvec(dd, d, dddirForward,
1122 extr_s+d, dd->ndim-d-1,
1123 extr_r+d, dd->ndim-d-1);
1127 for (d1 = d; d1 < dd->ndim-1; d1++)
1129 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1130 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1131 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1137 for (p = 0; p < npulse; p++)
1139 /* Communicate all the zone information backward */
1140 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1142 dd_sendrecv_ddzone(dd, d, dddirBackward,
1149 for (d1 = d+1; d1 < dd->ndim; d1++)
1151 /* Determine the decrease of maximum required
1152 * communication height along d1 due to the distance along d,
1153 * this avoids a lot of useless atom communication.
1155 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1157 if (ddbox->tric_dir[dim])
1159 /* c is the off-diagonal coupling between the cell planes
1160 * along directions d and d1.
1162 c = ddbox->v[dim][dd->dim[d1]][dim];
1168 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1171 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1175 /* A negative value signals out of range */
1181 /* Accumulate the extremes over all pulses */
1182 for (i = 0; i < buf_size; i++)
1186 buf_e[i] = buf_r[i];
1192 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1193 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1194 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1197 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1205 if (bUse && dh[d1] >= 0)
1207 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1208 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1211 /* Copy the received buffer to the send buffer,
1212 * to pass the data through with the next pulse.
1214 buf_s[i] = buf_r[i];
1216 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1217 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1219 /* Store the extremes */
1222 for (d1 = d; d1 < dd->ndim-1; d1++)
1224 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1225 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1226 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1230 if (d == 1 || (d == 0 && dd->ndim == 3))
1232 for (i = d; i < 2; i++)
1234 comm->zone_d2[1-d][i] = buf_e[pos];
1240 comm->zone_d1[1] = buf_e[pos];
1250 for (i = 0; i < 2; i++)
1254 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1256 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1257 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1263 for (i = 0; i < 2; i++)
1265 for (j = 0; j < 2; j++)
1269 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1271 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1272 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1276 for (d = 1; d < dd->ndim; d++)
1278 comm->cell_f_max0[d] = extr_s[d-1][0];
1279 comm->cell_f_min1[d] = extr_s[d-1][1];
1282 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1283 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1288 static void dd_collect_cg(gmx_domdec_t *dd,
1289 t_state *state_local)
1291 gmx_domdec_master_t *ma = NULL;
1292 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1295 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1297 /* The master has the correct distribution */
1301 if (state_local->ddp_count == dd->ddp_count)
1303 ncg_home = dd->ncg_home;
1305 nat_home = dd->nat_home;
1307 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1309 cgs_gl = &dd->comm->cgs_gl;
1311 ncg_home = state_local->ncg_gl;
1312 cg = state_local->cg_gl;
1314 for (i = 0; i < ncg_home; i++)
1316 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1321 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1324 buf2[0] = dd->ncg_home;
1325 buf2[1] = dd->nat_home;
1335 /* Collect the charge group and atom counts on the master */
1336 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1341 for (i = 0; i < dd->nnodes; i++)
1343 ma->ncg[i] = ma->ibuf[2*i];
1344 ma->nat[i] = ma->ibuf[2*i+1];
1345 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1348 /* Make byte counts and indices */
1349 for (i = 0; i < dd->nnodes; i++)
1351 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1352 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1356 fprintf(debug, "Initial charge group distribution: ");
1357 for (i = 0; i < dd->nnodes; i++)
1359 fprintf(debug, " %d", ma->ncg[i]);
1361 fprintf(debug, "\n");
1365 /* Collect the charge group indices on the master */
1367 dd->ncg_home*sizeof(int), dd->index_gl,
1368 DDMASTER(dd) ? ma->ibuf : NULL,
1369 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1370 DDMASTER(dd) ? ma->cg : NULL);
1372 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1375 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1378 gmx_domdec_master_t *ma;
1379 int n, i, c, a, nalloc = 0;
1388 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1389 dd->rank, dd->mpi_comm_all);
1394 /* Copy the master coordinates to the global array */
1395 cgs_gl = &dd->comm->cgs_gl;
1397 n = DDMASTERRANK(dd);
1399 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1401 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1403 copy_rvec(lv[a++], v[c]);
1407 for (n = 0; n < dd->nnodes; n++)
1411 if (ma->nat[n] > nalloc)
1413 nalloc = over_alloc_dd(ma->nat[n]);
1414 srenew(buf, nalloc);
1417 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1418 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1421 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1423 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1425 copy_rvec(buf[a++], v[c]);
1434 static void get_commbuffer_counts(gmx_domdec_t *dd,
1435 int **counts, int **disps)
1437 gmx_domdec_master_t *ma;
1442 /* Make the rvec count and displacment arrays */
1444 *disps = ma->ibuf + dd->nnodes;
1445 for (n = 0; n < dd->nnodes; n++)
1447 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1448 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1452 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1455 gmx_domdec_master_t *ma;
1456 int *rcounts = NULL, *disps = NULL;
1465 get_commbuffer_counts(dd, &rcounts, &disps);
1470 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1474 cgs_gl = &dd->comm->cgs_gl;
1477 for (n = 0; n < dd->nnodes; n++)
1479 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1481 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1483 copy_rvec(buf[a++], v[c]);
1490 void dd_collect_vec(gmx_domdec_t *dd,
1491 t_state *state_local, rvec *lv, rvec *v)
1493 gmx_domdec_master_t *ma;
1494 int n, i, c, a, nalloc = 0;
1497 dd_collect_cg(dd, state_local);
1499 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1501 dd_collect_vec_sendrecv(dd, lv, v);
1505 dd_collect_vec_gatherv(dd, lv, v);
1510 void dd_collect_state(gmx_domdec_t *dd,
1511 t_state *state_local, t_state *state)
1515 nh = state->nhchainlength;
1519 for (i = 0; i < efptNR; i++)
1521 state->lambda[i] = state_local->lambda[i];
1523 state->fep_state = state_local->fep_state;
1524 state->veta = state_local->veta;
1525 state->vol0 = state_local->vol0;
1526 copy_mat(state_local->box, state->box);
1527 copy_mat(state_local->boxv, state->boxv);
1528 copy_mat(state_local->svir_prev, state->svir_prev);
1529 copy_mat(state_local->fvir_prev, state->fvir_prev);
1530 copy_mat(state_local->pres_prev, state->pres_prev);
1533 for (i = 0; i < state_local->ngtc; i++)
1535 for (j = 0; j < nh; j++)
1537 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1538 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1540 state->therm_integral[i] = state_local->therm_integral[i];
1542 for (i = 0; i < state_local->nnhpres; i++)
1544 for (j = 0; j < nh; j++)
1546 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1547 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1551 for (est = 0; est < estNR; est++)
1553 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1558 dd_collect_vec(dd, state_local, state_local->x, state->x);
1561 dd_collect_vec(dd, state_local, state_local->v, state->v);
1564 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1567 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1570 if (state->nrngi == 1)
1574 for (i = 0; i < state_local->nrng; i++)
1576 state->ld_rng[i] = state_local->ld_rng[i];
1582 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1583 state_local->ld_rng, state->ld_rng);
1587 if (state->nrngi == 1)
1591 state->ld_rngi[0] = state_local->ld_rngi[0];
1596 dd_gather(dd, sizeof(state->ld_rngi[0]),
1597 state_local->ld_rngi, state->ld_rngi);
1600 case estDISRE_INITF:
1601 case estDISRE_RM3TAV:
1602 case estORIRE_INITF:
1606 gmx_incons("Unknown state entry encountered in dd_collect_state");
1612 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1618 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1621 state->nalloc = over_alloc_dd(nalloc);
1623 for (est = 0; est < estNR; est++)
1625 if (EST_DISTR(est) && (state->flags & (1<<est)))
1630 srenew(state->x, state->nalloc);
1633 srenew(state->v, state->nalloc);
1636 srenew(state->sd_X, state->nalloc);
1639 srenew(state->cg_p, state->nalloc);
1643 case estDISRE_INITF:
1644 case estDISRE_RM3TAV:
1645 case estORIRE_INITF:
1647 /* No reallocation required */
1650 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1657 srenew(*f, state->nalloc);
1661 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1664 if (nalloc > fr->cg_nalloc)
1668 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1670 fr->cg_nalloc = over_alloc_dd(nalloc);
1671 srenew(fr->cginfo, fr->cg_nalloc);
1672 if (fr->cutoff_scheme == ecutsGROUP)
1674 srenew(fr->cg_cm, fr->cg_nalloc);
1677 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1679 /* We don't use charge groups, we use x in state to set up
1680 * the atom communication.
1682 dd_realloc_state(state, f, nalloc);
1686 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1689 gmx_domdec_master_t *ma;
1690 int n, i, c, a, nalloc = 0;
1697 for (n = 0; n < dd->nnodes; n++)
1701 if (ma->nat[n] > nalloc)
1703 nalloc = over_alloc_dd(ma->nat[n]);
1704 srenew(buf, nalloc);
1706 /* Use lv as a temporary buffer */
1708 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1710 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1712 copy_rvec(v[c], buf[a++]);
1715 if (a != ma->nat[n])
1717 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1722 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1723 DDRANK(dd, n), n, dd->mpi_comm_all);
1728 n = DDMASTERRANK(dd);
1730 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1732 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1734 copy_rvec(v[c], lv[a++]);
1741 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1742 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1747 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1750 gmx_domdec_master_t *ma;
1751 int *scounts = NULL, *disps = NULL;
1752 int n, i, c, a, nalloc = 0;
1759 get_commbuffer_counts(dd, &scounts, &disps);
1763 for (n = 0; n < dd->nnodes; n++)
1765 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1767 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1769 copy_rvec(v[c], buf[a++]);
1775 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1778 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1780 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1782 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1786 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1790 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1791 t_state *state, t_state *state_local,
1796 nh = state->nhchainlength;
1800 for (i = 0; i < efptNR; i++)
1802 state_local->lambda[i] = state->lambda[i];
1804 state_local->fep_state = state->fep_state;
1805 state_local->veta = state->veta;
1806 state_local->vol0 = state->vol0;
1807 copy_mat(state->box, state_local->box);
1808 copy_mat(state->box_rel, state_local->box_rel);
1809 copy_mat(state->boxv, state_local->boxv);
1810 copy_mat(state->svir_prev, state_local->svir_prev);
1811 copy_mat(state->fvir_prev, state_local->fvir_prev);
1812 for (i = 0; i < state_local->ngtc; i++)
1814 for (j = 0; j < nh; j++)
1816 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1817 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1819 state_local->therm_integral[i] = state->therm_integral[i];
1821 for (i = 0; i < state_local->nnhpres; i++)
1823 for (j = 0; j < nh; j++)
1825 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1826 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1830 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1831 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1832 dd_bcast(dd, sizeof(real), &state_local->veta);
1833 dd_bcast(dd, sizeof(real), &state_local->vol0);
1834 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1835 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1836 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1837 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1838 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1839 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1840 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1841 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1842 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1843 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1845 if (dd->nat_home > state_local->nalloc)
1847 dd_realloc_state(state_local, f, dd->nat_home);
1849 for (i = 0; i < estNR; i++)
1851 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1856 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1859 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1862 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1865 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1868 if (state->nrngi == 1)
1871 state_local->nrng*sizeof(state_local->ld_rng[0]),
1872 state->ld_rng, state_local->ld_rng);
1877 state_local->nrng*sizeof(state_local->ld_rng[0]),
1878 state->ld_rng, state_local->ld_rng);
1882 if (state->nrngi == 1)
1884 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1885 state->ld_rngi, state_local->ld_rngi);
1889 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1890 state->ld_rngi, state_local->ld_rngi);
1893 case estDISRE_INITF:
1894 case estDISRE_RM3TAV:
1895 case estORIRE_INITF:
1897 /* Not implemented yet */
1900 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1906 static char dim2char(int dim)
1912 case XX: c = 'X'; break;
1913 case YY: c = 'Y'; break;
1914 case ZZ: c = 'Z'; break;
1915 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1921 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1922 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1924 rvec grid_s[2], *grid_r = NULL, cx, r;
1925 char fname[STRLEN], format[STRLEN], buf[22];
1927 int a, i, d, z, y, x;
1931 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1932 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1936 snew(grid_r, 2*dd->nnodes);
1939 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1943 for (d = 0; d < DIM; d++)
1945 for (i = 0; i < DIM; i++)
1953 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1955 tric[d][i] = box[i][d]/box[i][i];
1964 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1965 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
1966 out = gmx_fio_fopen(fname, "w");
1967 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1969 for (i = 0; i < dd->nnodes; i++)
1971 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1972 for (d = 0; d < DIM; d++)
1974 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1976 for (z = 0; z < 2; z++)
1978 for (y = 0; y < 2; y++)
1980 for (x = 0; x < 2; x++)
1982 cx[XX] = grid_r[i*2+x][XX];
1983 cx[YY] = grid_r[i*2+y][YY];
1984 cx[ZZ] = grid_r[i*2+z][ZZ];
1986 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
1987 ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
1991 for (d = 0; d < DIM; d++)
1993 for (x = 0; x < 4; x++)
1997 case 0: y = 1 + i*8 + 2*x; break;
1998 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
1999 case 2: y = 1 + i*8 + x; break;
2001 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2005 gmx_fio_fclose(out);
2010 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2011 gmx_mtop_t *mtop, t_commrec *cr,
2012 int natoms, rvec x[], matrix box)
2014 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2016 int i, ii, resnr, c;
2017 char *atomname, *resname;
2024 natoms = dd->comm->nat[ddnatVSITE];
2027 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2029 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
2030 sprintf(format4, "%s%s\n", get_pdbformat4(), "%6.2f%6.2f");
2032 out = gmx_fio_fopen(fname, "w");
2034 fprintf(out, "TITLE %s\n", title);
2035 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2036 for (i = 0; i < natoms; i++)
2038 ii = dd->gatindex[i];
2039 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2040 if (i < dd->comm->nat[ddnatZONE])
2043 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2049 else if (i < dd->comm->nat[ddnatVSITE])
2051 b = dd->comm->zones.n;
2055 b = dd->comm->zones.n + 1;
2057 fprintf(out, strlen(atomname) < 4 ? format : format4,
2058 "ATOM", (ii+1)%100000,
2059 atomname, resname, ' ', resnr%10000, ' ',
2060 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2062 fprintf(out, "TER\n");
2064 gmx_fio_fclose(out);
2067 real dd_cutoff_mbody(gmx_domdec_t *dd)
2069 gmx_domdec_comm_t *comm;
2076 if (comm->bInterCGBondeds)
2078 if (comm->cutoff_mbody > 0)
2080 r = comm->cutoff_mbody;
2084 /* cutoff_mbody=0 means we do not have DLB */
2085 r = comm->cellsize_min[dd->dim[0]];
2086 for (di = 1; di < dd->ndim; di++)
2088 r = min(r, comm->cellsize_min[dd->dim[di]]);
2090 if (comm->bBondComm)
2092 r = max(r, comm->cutoff_mbody);
2096 r = min(r, comm->cutoff);
2104 real dd_cutoff_twobody(gmx_domdec_t *dd)
2108 r_mb = dd_cutoff_mbody(dd);
2110 return max(dd->comm->cutoff, r_mb);
2114 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2118 nc = dd->nc[dd->comm->cartpmedim];
2119 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2120 copy_ivec(coord, coord_pme);
2121 coord_pme[dd->comm->cartpmedim] =
2122 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2125 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2127 /* Here we assign a PME node to communicate with this DD node
2128 * by assuming that the major index of both is x.
2129 * We add cr->npmenodes/2 to obtain an even distribution.
2131 return (ddindex*npme + npme/2)/ndd;
2134 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2136 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2139 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2141 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2144 static int *dd_pmenodes(t_commrec *cr)
2149 snew(pmenodes, cr->npmenodes);
2151 for (i = 0; i < cr->dd->nnodes; i++)
2153 p0 = cr_ddindex2pmeindex(cr, i);
2154 p1 = cr_ddindex2pmeindex(cr, i+1);
2155 if (i+1 == cr->dd->nnodes || p1 > p0)
2159 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2161 pmenodes[n] = i + 1 + n;
2169 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2172 ivec coords, coords_pme, nc;
2177 if (dd->comm->bCartesian) {
2178 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2179 dd_coords2pmecoords(dd,coords,coords_pme);
2180 copy_ivec(dd->ntot,nc);
2181 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2182 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2184 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2186 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2192 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2197 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2199 gmx_domdec_comm_t *comm;
2201 int ddindex, nodeid = -1;
2203 comm = cr->dd->comm;
2208 if (comm->bCartesianPP_PME)
2211 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2216 ddindex = dd_index(cr->dd->nc, coords);
2217 if (comm->bCartesianPP)
2219 nodeid = comm->ddindex2simnodeid[ddindex];
2225 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2237 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2240 gmx_domdec_comm_t *comm;
2241 ivec coord, coord_pme;
2248 /* This assumes a uniform x domain decomposition grid cell size */
2249 if (comm->bCartesianPP_PME)
2252 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2253 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2255 /* This is a PP node */
2256 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2257 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2261 else if (comm->bCartesianPP)
2263 if (sim_nodeid < dd->nnodes)
2265 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2270 /* This assumes DD cells with identical x coordinates
2271 * are numbered sequentially.
2273 if (dd->comm->pmenodes == NULL)
2275 if (sim_nodeid < dd->nnodes)
2277 /* The DD index equals the nodeid */
2278 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2284 while (sim_nodeid > dd->comm->pmenodes[i])
2288 if (sim_nodeid < dd->comm->pmenodes[i])
2290 pmenode = dd->comm->pmenodes[i];
2298 void get_pme_nnodes(const gmx_domdec_t *dd,
2299 int *npmenodes_x, int *npmenodes_y)
2303 *npmenodes_x = dd->comm->npmenodes_x;
2304 *npmenodes_y = dd->comm->npmenodes_y;
2313 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2315 gmx_bool bPMEOnlyNode;
2317 if (DOMAINDECOMP(cr))
2319 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2323 bPMEOnlyNode = FALSE;
2326 return bPMEOnlyNode;
2329 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2330 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2334 ivec coord, coord_pme;
2338 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2341 for (x = 0; x < dd->nc[XX]; x++)
2343 for (y = 0; y < dd->nc[YY]; y++)
2345 for (z = 0; z < dd->nc[ZZ]; z++)
2347 if (dd->comm->bCartesianPP_PME)
2352 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2353 if (dd->ci[XX] == coord_pme[XX] &&
2354 dd->ci[YY] == coord_pme[YY] &&
2355 dd->ci[ZZ] == coord_pme[ZZ])
2357 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2362 /* The slab corresponds to the nodeid in the PME group */
2363 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2365 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2372 /* The last PP-only node is the peer node */
2373 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2377 fprintf(debug, "Receive coordinates from PP nodes:");
2378 for (x = 0; x < *nmy_ddnodes; x++)
2380 fprintf(debug, " %d", (*my_ddnodes)[x]);
2382 fprintf(debug, "\n");
2386 static gmx_bool receive_vir_ener(t_commrec *cr)
2388 gmx_domdec_comm_t *comm;
2389 int pmenode, coords[DIM], rank;
2393 if (cr->npmenodes < cr->dd->nnodes)
2395 comm = cr->dd->comm;
2396 if (comm->bCartesianPP_PME)
2398 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2400 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2401 coords[comm->cartpmedim]++;
2402 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2404 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2405 if (dd_simnode2pmenode(cr, rank) == pmenode)
2407 /* This is not the last PP node for pmenode */
2415 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2416 if (cr->sim_nodeid+1 < cr->nnodes &&
2417 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2419 /* This is not the last PP node for pmenode */
2428 static void set_zones_ncg_home(gmx_domdec_t *dd)
2430 gmx_domdec_zones_t *zones;
2433 zones = &dd->comm->zones;
2435 zones->cg_range[0] = 0;
2436 for (i = 1; i < zones->n+1; i++)
2438 zones->cg_range[i] = dd->ncg_home;
2440 /* zone_ncg1[0] should always be equal to ncg_home */
2441 dd->comm->zone_ncg1[0] = dd->ncg_home;
2444 static void rebuild_cgindex(gmx_domdec_t *dd,
2445 const int *gcgs_index, t_state *state)
2447 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2450 dd_cg_gl = dd->index_gl;
2451 cgindex = dd->cgindex;
2454 for (i = 0; i < state->ncg_gl; i++)
2458 dd_cg_gl[i] = cg_gl;
2459 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2463 dd->ncg_home = state->ncg_gl;
2466 set_zones_ncg_home(dd);
2469 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2471 while (cg >= cginfo_mb->cg_end)
2476 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2479 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2480 t_forcerec *fr, char *bLocalCG)
2482 cginfo_mb_t *cginfo_mb;
2488 cginfo_mb = fr->cginfo_mb;
2489 cginfo = fr->cginfo;
2491 for (cg = cg0; cg < cg1; cg++)
2493 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2497 if (bLocalCG != NULL)
2499 for (cg = cg0; cg < cg1; cg++)
2501 bLocalCG[index_gl[cg]] = TRUE;
2506 static void make_dd_indices(gmx_domdec_t *dd,
2507 const int *gcgs_index, int cg_start)
2509 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2510 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2515 bLocalCG = dd->comm->bLocalCG;
2517 if (dd->nat_tot > dd->gatindex_nalloc)
2519 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2520 srenew(dd->gatindex, dd->gatindex_nalloc);
2523 nzone = dd->comm->zones.n;
2524 zone2cg = dd->comm->zones.cg_range;
2525 zone_ncg1 = dd->comm->zone_ncg1;
2526 index_gl = dd->index_gl;
2527 gatindex = dd->gatindex;
2528 bCGs = dd->comm->bCGs;
2530 if (zone2cg[1] != dd->ncg_home)
2532 gmx_incons("dd->ncg_zone is not up to date");
2535 /* Make the local to global and global to local atom index */
2536 a = dd->cgindex[cg_start];
2537 for (zone = 0; zone < nzone; zone++)
2545 cg0 = zone2cg[zone];
2547 cg1 = zone2cg[zone+1];
2548 cg1_p1 = cg0 + zone_ncg1[zone];
2550 for (cg = cg0; cg < cg1; cg++)
2555 /* Signal that this cg is from more than one pulse away */
2558 cg_gl = index_gl[cg];
2561 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2564 ga2la_set(dd->ga2la, a_gl, a, zone1);
2570 gatindex[a] = cg_gl;
2571 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2578 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2581 int ncg, i, ngl, nerr;
2584 if (bLocalCG == NULL)
2588 for (i = 0; i < dd->ncg_tot; i++)
2590 if (!bLocalCG[dd->index_gl[i]])
2593 "DD node %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);
2598 for (i = 0; i < ncg_sys; i++)
2605 if (ngl != dd->ncg_tot)
2607 fprintf(stderr, "DD node %d, %s: In bLocalCG %d cgs are marked as local, whereas there are %d\n", dd->rank, where, ngl, dd->ncg_tot);
2614 static void check_index_consistency(gmx_domdec_t *dd,
2615 int natoms_sys, int ncg_sys,
2618 int nerr, ngl, i, a, cell;
2623 if (dd->comm->DD_debug > 1)
2625 snew(have, natoms_sys);
2626 for (a = 0; a < dd->nat_tot; a++)
2628 if (have[dd->gatindex[a]] > 0)
2630 fprintf(stderr, "DD node %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
2634 have[dd->gatindex[a]] = a + 1;
2640 snew(have, dd->nat_tot);
2643 for (i = 0; i < natoms_sys; i++)
2645 if (ga2la_get(dd->ga2la, i, &a, &cell))
2647 if (a >= dd->nat_tot)
2649 fprintf(stderr, "DD node %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);
2655 if (dd->gatindex[a] != i)
2657 fprintf(stderr, "DD node %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);
2664 if (ngl != dd->nat_tot)
2667 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2668 dd->rank, where, ngl, dd->nat_tot);
2670 for (a = 0; a < dd->nat_tot; a++)
2675 "DD node %d, %s: local atom %d, global %d has no global index\n",
2676 dd->rank, where, a+1, dd->gatindex[a]+1);
2681 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2685 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2686 dd->rank, where, nerr);
2690 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2697 /* Clear the whole list without searching */
2698 ga2la_clear(dd->ga2la);
2702 for (i = a_start; i < dd->nat_tot; i++)
2704 ga2la_del(dd->ga2la, dd->gatindex[i]);
2708 bLocalCG = dd->comm->bLocalCG;
2711 for (i = cg_start; i < dd->ncg_tot; i++)
2713 bLocalCG[dd->index_gl[i]] = FALSE;
2717 dd_clear_local_vsite_indices(dd);
2719 if (dd->constraints)
2721 dd_clear_local_constraint_indices(dd);
2725 /* This function should be used for moving the domain boudaries during DLB,
2726 * for obtaining the minimum cell size. It checks the initially set limit
2727 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2728 * and, possibly, a longer cut-off limit set for PME load balancing.
2730 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2734 cellsize_min = comm->cellsize_min[dim];
2736 if (!comm->bVacDLBNoLimit)
2738 /* The cut-off might have changed, e.g. by PME load balacning,
2739 * from the value used to set comm->cellsize_min, so check it.
2741 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2743 if (comm->bPMELoadBalDLBLimits)
2745 /* Check for the cut-off limit set by the PME load balancing */
2746 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2750 return cellsize_min;
2753 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2756 real grid_jump_limit;
2758 /* The distance between the boundaries of cells at distance
2759 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2760 * and by the fact that cells should not be shifted by more than
2761 * half their size, such that cg's only shift by one cell
2762 * at redecomposition.
2764 grid_jump_limit = comm->cellsize_limit;
2765 if (!comm->bVacDLBNoLimit)
2767 if (comm->bPMELoadBalDLBLimits)
2769 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2771 grid_jump_limit = max(grid_jump_limit,
2772 cutoff/comm->cd[dim_ind].np);
2775 return grid_jump_limit;
2778 static gmx_bool check_grid_jump(gmx_large_int_t step,
2784 gmx_domdec_comm_t *comm;
2793 for (d = 1; d < dd->ndim; d++)
2796 limit = grid_jump_limit(comm, cutoff, d);
2797 bfac = ddbox->box_size[dim];
2798 if (ddbox->tric_dir[dim])
2800 bfac *= ddbox->skew_fac[dim];
2802 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2803 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2811 /* This error should never be triggered under normal
2812 * circumstances, but you never know ...
2814 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 less nodes might avoid this issue.",
2815 gmx_step_str(step, buf),
2816 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2824 static int dd_load_count(gmx_domdec_comm_t *comm)
2826 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2829 static float dd_force_load(gmx_domdec_comm_t *comm)
2836 if (comm->eFlop > 1)
2838 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2843 load = comm->cycl[ddCyclF];
2844 if (comm->cycl_n[ddCyclF] > 1)
2846 /* Subtract the maximum of the last n cycle counts
2847 * to get rid of possible high counts due to other soures,
2848 * for instance system activity, that would otherwise
2849 * affect the dynamic load balancing.
2851 load -= comm->cycl_max[ddCyclF];
2858 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2860 gmx_domdec_comm_t *comm;
2865 snew(*dim_f, dd->nc[dim]+1);
2867 for (i = 1; i < dd->nc[dim]; i++)
2869 if (comm->slb_frac[dim])
2871 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2875 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2878 (*dim_f)[dd->nc[dim]] = 1;
2881 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2883 int pmeindex, slab, nso, i;
2886 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2892 ddpme->dim = dimind;
2894 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2896 ddpme->nslab = (ddpme->dim == 0 ?
2897 dd->comm->npmenodes_x :
2898 dd->comm->npmenodes_y);
2900 if (ddpme->nslab <= 1)
2905 nso = dd->comm->npmenodes/ddpme->nslab;
2906 /* Determine for each PME slab the PP location range for dimension dim */
2907 snew(ddpme->pp_min, ddpme->nslab);
2908 snew(ddpme->pp_max, ddpme->nslab);
2909 for (slab = 0; slab < ddpme->nslab; slab++)
2911 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2912 ddpme->pp_max[slab] = 0;
2914 for (i = 0; i < dd->nnodes; i++)
2916 ddindex2xyz(dd->nc, i, xyz);
2917 /* For y only use our y/z slab.
2918 * This assumes that the PME x grid size matches the DD grid size.
2920 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2922 pmeindex = ddindex2pmeindex(dd, i);
2925 slab = pmeindex/nso;
2929 slab = pmeindex % ddpme->nslab;
2931 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2932 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2936 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2939 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2941 if (dd->comm->ddpme[0].dim == XX)
2943 return dd->comm->ddpme[0].maxshift;
2951 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2953 if (dd->comm->ddpme[0].dim == YY)
2955 return dd->comm->ddpme[0].maxshift;
2957 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2959 return dd->comm->ddpme[1].maxshift;
2967 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2968 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2970 gmx_domdec_comm_t *comm;
2973 real range, pme_boundary;
2977 nc = dd->nc[ddpme->dim];
2980 if (!ddpme->dim_match)
2982 /* PP decomposition is not along dim: the worst situation */
2985 else if (ns <= 3 || (bUniform && ns == nc))
2987 /* The optimal situation */
2992 /* We need to check for all pme nodes which nodes they
2993 * could possibly need to communicate with.
2995 xmin = ddpme->pp_min;
2996 xmax = ddpme->pp_max;
2997 /* Allow for atoms to be maximally 2/3 times the cut-off
2998 * out of their DD cell. This is a reasonable balance between
2999 * between performance and support for most charge-group/cut-off
3002 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3003 /* Avoid extra communication when we are exactly at a boundary */
3007 for (s = 0; s < ns; s++)
3009 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3010 pme_boundary = (real)s/ns;
3013 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3015 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3019 pme_boundary = (real)(s+1)/ns;
3022 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3024 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3031 ddpme->maxshift = sh;
3035 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3036 ddpme->dim, ddpme->maxshift);
3040 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3044 for (d = 0; d < dd->ndim; d++)
3047 if (dim < ddbox->nboundeddim &&
3048 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3049 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3051 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",
3052 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3053 dd->nc[dim], dd->comm->cellsize_limit);
3058 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3059 gmx_bool bMaster, ivec npulse)
3061 gmx_domdec_comm_t *comm;
3064 real *cell_x, cell_dx, cellsize;
3068 for (d = 0; d < DIM; d++)
3070 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3072 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3075 cell_dx = ddbox->box_size[d]/dd->nc[d];
3078 for (j = 0; j < dd->nc[d]+1; j++)
3080 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3085 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3086 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3088 cellsize = cell_dx*ddbox->skew_fac[d];
3089 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3093 cellsize_min[d] = cellsize;
3097 /* Statically load balanced grid */
3098 /* Also when we are not doing a master distribution we determine
3099 * all cell borders in a loop to obtain identical values
3100 * to the master distribution case and to determine npulse.
3104 cell_x = dd->ma->cell_x[d];
3108 snew(cell_x, dd->nc[d]+1);
3110 cell_x[0] = ddbox->box0[d];
3111 for (j = 0; j < dd->nc[d]; j++)
3113 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3114 cell_x[j+1] = cell_x[j] + cell_dx;
3115 cellsize = cell_dx*ddbox->skew_fac[d];
3116 while (cellsize*npulse[d] < comm->cutoff &&
3117 npulse[d] < dd->nc[d]-1)
3121 cellsize_min[d] = min(cellsize_min[d], cellsize);
3125 comm->cell_x0[d] = cell_x[dd->ci[d]];
3126 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3130 /* The following limitation is to avoid that a cell would receive
3131 * some of its own home charge groups back over the periodic boundary.
3132 * Double charge groups cause trouble with the global indices.
3134 if (d < ddbox->npbcdim &&
3135 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3137 gmx_fatal_collective(FARGS, NULL, dd,
3138 "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",
3139 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3141 dd->nc[d], dd->nc[d],
3142 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3146 if (!comm->bDynLoadBal)
3148 copy_rvec(cellsize_min, comm->cellsize_min);
3151 for (d = 0; d < comm->npmedecompdim; d++)
3153 set_pme_maxshift(dd, &comm->ddpme[d],
3154 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3155 comm->ddpme[d].slb_dim_f);
3160 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3161 int d, int dim, gmx_domdec_root_t *root,
3163 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3165 gmx_domdec_comm_t *comm;
3166 int ncd, i, j, nmin, nmin_old;
3167 gmx_bool bLimLo, bLimHi;
3169 real fac, halfway, cellsize_limit_f_i, region_size;
3170 gmx_bool bPBC, bLastHi = FALSE;
3171 int nrange[] = {range[0], range[1]};
3173 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3179 bPBC = (dim < ddbox->npbcdim);
3181 cell_size = root->buf_ncd;
3185 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3188 /* First we need to check if the scaling does not make cells
3189 * smaller than the smallest allowed size.
3190 * We need to do this iteratively, since if a cell is too small,
3191 * it needs to be enlarged, which makes all the other cells smaller,
3192 * which could in turn make another cell smaller than allowed.
3194 for (i = range[0]; i < range[1]; i++)
3196 root->bCellMin[i] = FALSE;
3202 /* We need the total for normalization */
3204 for (i = range[0]; i < range[1]; i++)
3206 if (root->bCellMin[i] == FALSE)
3208 fac += cell_size[i];
3211 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3212 /* Determine the cell boundaries */
3213 for (i = range[0]; i < range[1]; i++)
3215 if (root->bCellMin[i] == FALSE)
3217 cell_size[i] *= fac;
3218 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3220 cellsize_limit_f_i = 0;
3224 cellsize_limit_f_i = cellsize_limit_f;
3226 if (cell_size[i] < cellsize_limit_f_i)
3228 root->bCellMin[i] = TRUE;
3229 cell_size[i] = cellsize_limit_f_i;
3233 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3236 while (nmin > nmin_old);
3239 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3240 /* For this check we should not use DD_CELL_MARGIN,
3241 * but a slightly smaller factor,
3242 * since rounding could get use below the limit.
3244 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3247 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",
3248 gmx_step_str(step, buf),
3249 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3250 ncd, comm->cellsize_min[dim]);
3253 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3257 /* Check if the boundary did not displace more than halfway
3258 * each of the cells it bounds, as this could cause problems,
3259 * especially when the differences between cell sizes are large.
3260 * If changes are applied, they will not make cells smaller
3261 * than the cut-off, as we check all the boundaries which
3262 * might be affected by a change and if the old state was ok,
3263 * the cells will at most be shrunk back to their old size.
3265 for (i = range[0]+1; i < range[1]; i++)
3267 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3268 if (root->cell_f[i] < halfway)
3270 root->cell_f[i] = halfway;
3271 /* Check if the change also causes shifts of the next boundaries */
3272 for (j = i+1; j < range[1]; j++)
3274 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3276 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3280 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3281 if (root->cell_f[i] > halfway)
3283 root->cell_f[i] = halfway;
3284 /* Check if the change also causes shifts of the next boundaries */
3285 for (j = i-1; j >= range[0]+1; j--)
3287 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3289 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3296 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3297 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3298 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3299 * for a and b nrange is used */
3302 /* Take care of the staggering of the cell boundaries */
3305 for (i = range[0]; i < range[1]; i++)
3307 root->cell_f_max0[i] = root->cell_f[i];
3308 root->cell_f_min1[i] = root->cell_f[i+1];
3313 for (i = range[0]+1; i < range[1]; i++)
3315 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3316 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3317 if (bLimLo && bLimHi)
3319 /* Both limits violated, try the best we can */
3320 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3321 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3322 nrange[0] = range[0];
3324 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3327 nrange[1] = range[1];
3328 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3334 /* root->cell_f[i] = root->bound_min[i]; */
3335 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3338 else if (bLimHi && !bLastHi)
3341 if (nrange[1] < range[1]) /* found a LimLo before */
3343 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3344 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3345 nrange[0] = nrange[1];
3347 root->cell_f[i] = root->bound_max[i];
3349 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3351 nrange[1] = range[1];
3354 if (nrange[1] < range[1]) /* found last a LimLo */
3356 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3357 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3358 nrange[0] = nrange[1];
3359 nrange[1] = range[1];
3360 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3362 else if (nrange[0] > range[0]) /* found at least one LimHi */
3364 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3371 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3372 int d, int dim, gmx_domdec_root_t *root,
3373 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3374 gmx_bool bUniform, gmx_large_int_t step)
3376 gmx_domdec_comm_t *comm;
3377 int ncd, d1, i, j, pos;
3379 real load_aver, load_i, imbalance, change, change_max, sc;
3380 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3384 int range[] = { 0, 0 };
3388 /* Convert the maximum change from the input percentage to a fraction */
3389 change_limit = comm->dlb_scale_lim*0.01;
3393 bPBC = (dim < ddbox->npbcdim);
3395 cell_size = root->buf_ncd;
3397 /* Store the original boundaries */
3398 for (i = 0; i < ncd+1; i++)
3400 root->old_cell_f[i] = root->cell_f[i];
3404 for (i = 0; i < ncd; i++)
3406 cell_size[i] = 1.0/ncd;
3409 else if (dd_load_count(comm))
3411 load_aver = comm->load[d].sum_m/ncd;
3413 for (i = 0; i < ncd; i++)
3415 /* Determine the relative imbalance of cell i */
3416 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3417 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3418 /* Determine the change of the cell size using underrelaxation */
3419 change = -relax*imbalance;
3420 change_max = max(change_max, max(change, -change));
3422 /* Limit the amount of scaling.
3423 * We need to use the same rescaling for all cells in one row,
3424 * otherwise the load balancing might not converge.
3427 if (change_max > change_limit)
3429 sc *= change_limit/change_max;
3431 for (i = 0; i < ncd; i++)
3433 /* Determine the relative imbalance of cell i */
3434 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3435 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3436 /* Determine the change of the cell size using underrelaxation */
3437 change = -sc*imbalance;
3438 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3442 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3443 cellsize_limit_f *= DD_CELL_MARGIN;
3444 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3445 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3446 if (ddbox->tric_dir[dim])
3448 cellsize_limit_f /= ddbox->skew_fac[dim];
3449 dist_min_f /= ddbox->skew_fac[dim];
3451 if (bDynamicBox && d > 0)
3453 dist_min_f *= DD_PRES_SCALE_MARGIN;
3455 if (d > 0 && !bUniform)
3457 /* Make sure that the grid is not shifted too much */
3458 for (i = 1; i < ncd; i++)
3460 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3462 gmx_incons("Inconsistent DD boundary staggering limits!");
3464 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3465 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3468 root->bound_min[i] += 0.5*space;
3470 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3471 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3474 root->bound_max[i] += 0.5*space;
3479 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3481 root->cell_f_max0[i-1] + dist_min_f,
3482 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3483 root->cell_f_min1[i] - dist_min_f);
3488 root->cell_f[0] = 0;
3489 root->cell_f[ncd] = 1;
3490 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3493 /* After the checks above, the cells should obey the cut-off
3494 * restrictions, but it does not hurt to check.
3496 for (i = 0; i < ncd; i++)
3500 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3501 dim, i, root->cell_f[i], root->cell_f[i+1]);
3504 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3505 root->cell_f[i+1] - root->cell_f[i] <
3506 cellsize_limit_f/DD_CELL_MARGIN)
3510 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3511 gmx_step_str(step, buf), dim2char(dim), i,
3512 (root->cell_f[i+1] - root->cell_f[i])
3513 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3518 /* Store the cell boundaries of the lower dimensions at the end */
3519 for (d1 = 0; d1 < d; d1++)
3521 root->cell_f[pos++] = comm->cell_f0[d1];
3522 root->cell_f[pos++] = comm->cell_f1[d1];
3525 if (d < comm->npmedecompdim)
3527 /* The master determines the maximum shift for
3528 * the coordinate communication between separate PME nodes.
3530 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3532 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3535 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3539 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3540 gmx_ddbox_t *ddbox, int dimind)
3542 gmx_domdec_comm_t *comm;
3547 /* Set the cell dimensions */
3548 dim = dd->dim[dimind];
3549 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3550 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3551 if (dim >= ddbox->nboundeddim)
3553 comm->cell_x0[dim] += ddbox->box0[dim];
3554 comm->cell_x1[dim] += ddbox->box0[dim];
3558 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3559 int d, int dim, real *cell_f_row,
3562 gmx_domdec_comm_t *comm;
3568 /* Each node would only need to know two fractions,
3569 * but it is probably cheaper to broadcast the whole array.
3571 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3572 0, comm->mpi_comm_load[d]);
3574 /* Copy the fractions for this dimension from the buffer */
3575 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3576 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3577 /* The whole array was communicated, so set the buffer position */
3578 pos = dd->nc[dim] + 1;
3579 for (d1 = 0; d1 <= d; d1++)
3583 /* Copy the cell fractions of the lower dimensions */
3584 comm->cell_f0[d1] = cell_f_row[pos++];
3585 comm->cell_f1[d1] = cell_f_row[pos++];
3587 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3589 /* Convert the communicated shift from float to int */
3590 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3593 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3597 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3598 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3599 gmx_bool bUniform, gmx_large_int_t step)
3601 gmx_domdec_comm_t *comm;
3603 gmx_bool bRowMember, bRowRoot;
3608 for (d = 0; d < dd->ndim; d++)
3613 for (d1 = d; d1 < dd->ndim; d1++)
3615 if (dd->ci[dd->dim[d1]] > 0)
3628 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3629 ddbox, bDynamicBox, bUniform, step);
3630 cell_f_row = comm->root[d]->cell_f;
3634 cell_f_row = comm->cell_f_row;
3636 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3641 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3645 /* This function assumes the box is static and should therefore
3646 * not be called when the box has changed since the last
3647 * call to dd_partition_system.
3649 for (d = 0; d < dd->ndim; d++)
3651 relative_to_absolute_cell_bounds(dd, ddbox, d);
3657 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3658 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3659 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3660 gmx_wallcycle_t wcycle)
3662 gmx_domdec_comm_t *comm;
3669 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3670 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3671 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3673 else if (bDynamicBox)
3675 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3678 /* Set the dimensions for which no DD is used */
3679 for (dim = 0; dim < DIM; dim++)
3681 if (dd->nc[dim] == 1)
3683 comm->cell_x0[dim] = 0;
3684 comm->cell_x1[dim] = ddbox->box_size[dim];
3685 if (dim >= ddbox->nboundeddim)
3687 comm->cell_x0[dim] += ddbox->box0[dim];
3688 comm->cell_x1[dim] += ddbox->box0[dim];
3694 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3697 gmx_domdec_comm_dim_t *cd;
3699 for (d = 0; d < dd->ndim; d++)
3701 cd = &dd->comm->cd[d];
3702 np = npulse[dd->dim[d]];
3703 if (np > cd->np_nalloc)
3707 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3708 dim2char(dd->dim[d]), np);
3710 if (DDMASTER(dd) && cd->np_nalloc > 0)
3712 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3714 srenew(cd->ind, np);
3715 for (i = cd->np_nalloc; i < np; i++)
3717 cd->ind[i].index = NULL;
3718 cd->ind[i].nalloc = 0;
3727 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3728 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3729 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3730 gmx_wallcycle_t wcycle)
3732 gmx_domdec_comm_t *comm;
3738 /* Copy the old cell boundaries for the cg displacement check */
3739 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3740 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3742 if (comm->bDynLoadBal)
3746 check_box_size(dd, ddbox);
3748 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3752 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3753 realloc_comm_ind(dd, npulse);
3758 for (d = 0; d < DIM; d++)
3760 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3761 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3766 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3768 rvec cell_ns_x0, rvec cell_ns_x1,
3769 gmx_large_int_t step)
3771 gmx_domdec_comm_t *comm;
3776 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3778 dim = dd->dim[dim_ind];
3780 /* Without PBC we don't have restrictions on the outer cells */
3781 if (!(dim >= ddbox->npbcdim &&
3782 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3783 comm->bDynLoadBal &&
3784 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3785 comm->cellsize_min[dim])
3788 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",
3789 gmx_step_str(step, buf), dim2char(dim),
3790 comm->cell_x1[dim] - comm->cell_x0[dim],
3791 ddbox->skew_fac[dim],
3792 dd->comm->cellsize_min[dim],
3793 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3797 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3799 /* Communicate the boundaries and update cell_ns_x0/1 */
3800 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3801 if (dd->bGridJump && dd->ndim > 1)
3803 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3808 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3812 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3820 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3821 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3830 static void check_screw_box(matrix box)
3832 /* Mathematical limitation */
3833 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3835 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3838 /* Limitation due to the asymmetry of the eighth shell method */
3839 if (box[ZZ][YY] != 0)
3841 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3845 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3846 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3849 gmx_domdec_master_t *ma;
3850 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3851 int i, icg, j, k, k0, k1, d, npbcdim;
3853 rvec box_size, cg_cm;
3855 real nrcg, inv_ncg, pos_d;
3857 gmx_bool bUnbounded, bScrew;
3861 if (tmp_ind == NULL)
3863 snew(tmp_nalloc, dd->nnodes);
3864 snew(tmp_ind, dd->nnodes);
3865 for (i = 0; i < dd->nnodes; i++)
3867 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3868 snew(tmp_ind[i], tmp_nalloc[i]);
3872 /* Clear the count */
3873 for (i = 0; i < dd->nnodes; i++)
3879 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3881 cgindex = cgs->index;
3883 /* Compute the center of geometry for all charge groups */
3884 for (icg = 0; icg < cgs->nr; icg++)
3887 k1 = cgindex[icg+1];
3891 copy_rvec(pos[k0], cg_cm);
3898 for (k = k0; (k < k1); k++)
3900 rvec_inc(cg_cm, pos[k]);
3902 for (d = 0; (d < DIM); d++)
3904 cg_cm[d] *= inv_ncg;
3907 /* Put the charge group in the box and determine the cell index */
3908 for (d = DIM-1; d >= 0; d--)
3911 if (d < dd->npbcdim)
3913 bScrew = (dd->bScrewPBC && d == XX);
3914 if (tric_dir[d] && dd->nc[d] > 1)
3916 /* Use triclinic coordintates for this dimension */
3917 for (j = d+1; j < DIM; j++)
3919 pos_d += cg_cm[j]*tcm[j][d];
3922 while (pos_d >= box[d][d])
3925 rvec_dec(cg_cm, box[d]);
3928 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3929 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3931 for (k = k0; (k < k1); k++)
3933 rvec_dec(pos[k], box[d]);
3936 pos[k][YY] = box[YY][YY] - pos[k][YY];
3937 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3944 rvec_inc(cg_cm, box[d]);
3947 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3948 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3950 for (k = k0; (k < k1); k++)
3952 rvec_inc(pos[k], box[d]);
3955 pos[k][YY] = box[YY][YY] - pos[k][YY];
3956 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3961 /* This could be done more efficiently */
3963 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3968 i = dd_index(dd->nc, ind);
3969 if (ma->ncg[i] == tmp_nalloc[i])
3971 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3972 srenew(tmp_ind[i], tmp_nalloc[i]);
3974 tmp_ind[i][ma->ncg[i]] = icg;
3976 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3980 for (i = 0; i < dd->nnodes; i++)
3983 for (k = 0; k < ma->ncg[i]; k++)
3985 ma->cg[k1++] = tmp_ind[i][k];
3988 ma->index[dd->nnodes] = k1;
3990 for (i = 0; i < dd->nnodes; i++)
4000 fprintf(fplog, "Charge group distribution at step %s:",
4001 gmx_step_str(step, buf));
4002 for (i = 0; i < dd->nnodes; i++)
4004 fprintf(fplog, " %d", ma->ncg[i]);
4006 fprintf(fplog, "\n");
4010 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4011 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4014 gmx_domdec_master_t *ma = NULL;
4017 int *ibuf, buf2[2] = { 0, 0 };
4018 gmx_bool bMaster = DDMASTER(dd);
4025 check_screw_box(box);
4028 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4030 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4031 for (i = 0; i < dd->nnodes; i++)
4033 ma->ibuf[2*i] = ma->ncg[i];
4034 ma->ibuf[2*i+1] = ma->nat[i];
4042 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4044 dd->ncg_home = buf2[0];
4045 dd->nat_home = buf2[1];
4046 dd->ncg_tot = dd->ncg_home;
4047 dd->nat_tot = dd->nat_home;
4048 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4050 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4051 srenew(dd->index_gl, dd->cg_nalloc);
4052 srenew(dd->cgindex, dd->cg_nalloc+1);
4056 for (i = 0; i < dd->nnodes; i++)
4058 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4059 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4064 DDMASTER(dd) ? ma->ibuf : NULL,
4065 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4066 DDMASTER(dd) ? ma->cg : NULL,
4067 dd->ncg_home*sizeof(int), dd->index_gl);
4069 /* Determine the home charge group sizes */
4071 for (i = 0; i < dd->ncg_home; i++)
4073 cg_gl = dd->index_gl[i];
4075 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4080 fprintf(debug, "Home charge groups:\n");
4081 for (i = 0; i < dd->ncg_home; i++)
4083 fprintf(debug, " %d", dd->index_gl[i]);
4086 fprintf(debug, "\n");
4089 fprintf(debug, "\n");
4093 static int compact_and_copy_vec_at(int ncg, int *move,
4096 rvec *src, gmx_domdec_comm_t *comm,
4099 int m, icg, i, i0, i1, nrcg;
4105 for (m = 0; m < DIM*2; m++)
4111 for (icg = 0; icg < ncg; icg++)
4113 i1 = cgindex[icg+1];
4119 /* Compact the home array in place */
4120 for (i = i0; i < i1; i++)
4122 copy_rvec(src[i], src[home_pos++]);
4128 /* Copy to the communication buffer */
4130 pos_vec[m] += 1 + vec*nrcg;
4131 for (i = i0; i < i1; i++)
4133 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4135 pos_vec[m] += (nvec - vec - 1)*nrcg;
4139 home_pos += i1 - i0;
4147 static int compact_and_copy_vec_cg(int ncg, int *move,
4149 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4152 int m, icg, i0, i1, nrcg;
4158 for (m = 0; m < DIM*2; m++)
4164 for (icg = 0; icg < ncg; icg++)
4166 i1 = cgindex[icg+1];
4172 /* Compact the home array in place */
4173 copy_rvec(src[icg], src[home_pos++]);
4179 /* Copy to the communication buffer */
4180 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4181 pos_vec[m] += 1 + nrcg*nvec;
4193 static int compact_ind(int ncg, int *move,
4194 int *index_gl, int *cgindex,
4196 gmx_ga2la_t ga2la, char *bLocalCG,
4199 int cg, nat, a0, a1, a, a_gl;
4204 for (cg = 0; cg < ncg; cg++)
4210 /* Compact the home arrays in place.
4211 * Anything that can be done here avoids access to global arrays.
4213 cgindex[home_pos] = nat;
4214 for (a = a0; a < a1; a++)
4217 gatindex[nat] = a_gl;
4218 /* The cell number stays 0, so we don't need to set it */
4219 ga2la_change_la(ga2la, a_gl, nat);
4222 index_gl[home_pos] = index_gl[cg];
4223 cginfo[home_pos] = cginfo[cg];
4224 /* The charge group remains local, so bLocalCG does not change */
4229 /* Clear the global indices */
4230 for (a = a0; a < a1; a++)
4232 ga2la_del(ga2la, gatindex[a]);
4236 bLocalCG[index_gl[cg]] = FALSE;
4240 cgindex[home_pos] = nat;
4245 static void clear_and_mark_ind(int ncg, int *move,
4246 int *index_gl, int *cgindex, int *gatindex,
4247 gmx_ga2la_t ga2la, char *bLocalCG,
4252 for (cg = 0; cg < ncg; cg++)
4258 /* Clear the global indices */
4259 for (a = a0; a < a1; a++)
4261 ga2la_del(ga2la, gatindex[a]);
4265 bLocalCG[index_gl[cg]] = FALSE;
4267 /* Signal that this cg has moved using the ns cell index.
4268 * Here we set it to -1. fill_grid will change it
4269 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4271 cell_index[cg] = -1;
4276 static void print_cg_move(FILE *fplog,
4278 gmx_large_int_t step, int cg, int dim, int dir,
4279 gmx_bool bHaveLimitdAndCMOld, real limitd,
4280 rvec cm_old, rvec cm_new, real pos_d)
4282 gmx_domdec_comm_t *comm;
4287 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4288 if (bHaveLimitdAndCMOld)
4290 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4291 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4295 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4296 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4298 fprintf(fplog, "distance out of cell %f\n",
4299 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4300 if (bHaveLimitdAndCMOld)
4302 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4303 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4305 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4306 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4307 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4309 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4310 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4312 comm->cell_x0[dim], comm->cell_x1[dim]);
4315 static void cg_move_error(FILE *fplog,
4317 gmx_large_int_t step, int cg, int dim, int dir,
4318 gmx_bool bHaveLimitdAndCMOld, real limitd,
4319 rvec cm_old, rvec cm_new, real pos_d)
4323 print_cg_move(fplog, dd, step, cg, dim, dir,
4324 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4326 print_cg_move(stderr, dd, step, cg, dim, dir,
4327 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4329 "A charge group moved too far between two domain decomposition steps\n"
4330 "This usually means that your system is not well equilibrated");
4333 static void rotate_state_atom(t_state *state, int a)
4337 for (est = 0; est < estNR; est++)
4339 if (EST_DISTR(est) && (state->flags & (1<<est)))
4344 /* Rotate the complete state; for a rectangular box only */
4345 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4346 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4349 state->v[a][YY] = -state->v[a][YY];
4350 state->v[a][ZZ] = -state->v[a][ZZ];
4353 state->sd_X[a][YY] = -state->sd_X[a][YY];
4354 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4357 state->cg_p[a][YY] = -state->cg_p[a][YY];
4358 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4360 case estDISRE_INITF:
4361 case estDISRE_RM3TAV:
4362 case estORIRE_INITF:
4364 /* These are distances, so not affected by rotation */
4367 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4373 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4375 if (natoms > comm->moved_nalloc)
4377 /* Contents should be preserved here */
4378 comm->moved_nalloc = over_alloc_dd(natoms);
4379 srenew(comm->moved, comm->moved_nalloc);
4385 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4388 ivec tric_dir, matrix tcm,
4389 rvec cell_x0, rvec cell_x1,
4390 rvec limitd, rvec limit0, rvec limit1,
4392 int cg_start, int cg_end,
4397 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4398 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4402 real inv_ncg, pos_d;
4405 npbcdim = dd->npbcdim;
4407 for (cg = cg_start; cg < cg_end; cg++)
4414 copy_rvec(state->x[k0], cm_new);
4421 for (k = k0; (k < k1); k++)
4423 rvec_inc(cm_new, state->x[k]);
4425 for (d = 0; (d < DIM); d++)
4427 cm_new[d] = inv_ncg*cm_new[d];
4432 /* Do pbc and check DD cell boundary crossings */
4433 for (d = DIM-1; d >= 0; d--)
4437 bScrew = (dd->bScrewPBC && d == XX);
4438 /* Determine the location of this cg in lattice coordinates */
4442 for (d2 = d+1; d2 < DIM; d2++)
4444 pos_d += cm_new[d2]*tcm[d2][d];
4447 /* Put the charge group in the triclinic unit-cell */
4448 if (pos_d >= cell_x1[d])
4450 if (pos_d >= limit1[d])
4452 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4453 cg_cm[cg], cm_new, pos_d);
4456 if (dd->ci[d] == dd->nc[d] - 1)
4458 rvec_dec(cm_new, state->box[d]);
4461 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4462 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4464 for (k = k0; (k < k1); k++)
4466 rvec_dec(state->x[k], state->box[d]);
4469 rotate_state_atom(state, k);
4474 else if (pos_d < cell_x0[d])
4476 if (pos_d < limit0[d])
4478 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4479 cg_cm[cg], cm_new, pos_d);
4484 rvec_inc(cm_new, state->box[d]);
4487 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4488 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4490 for (k = k0; (k < k1); k++)
4492 rvec_inc(state->x[k], state->box[d]);
4495 rotate_state_atom(state, k);
4501 else if (d < npbcdim)
4503 /* Put the charge group in the rectangular unit-cell */
4504 while (cm_new[d] >= state->box[d][d])
4506 rvec_dec(cm_new, state->box[d]);
4507 for (k = k0; (k < k1); k++)
4509 rvec_dec(state->x[k], state->box[d]);
4512 while (cm_new[d] < 0)
4514 rvec_inc(cm_new, state->box[d]);
4515 for (k = k0; (k < k1); k++)
4517 rvec_inc(state->x[k], state->box[d]);
4523 copy_rvec(cm_new, cg_cm[cg]);
4525 /* Determine where this cg should go */
4528 for (d = 0; d < dd->ndim; d++)
4533 flag |= DD_FLAG_FW(d);
4539 else if (dev[dim] == -1)
4541 flag |= DD_FLAG_BW(d);
4544 if (dd->nc[dim] > 2)
4555 /* Temporarily store the flag in move */
4556 move[cg] = mc + flag;
4560 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4561 gmx_domdec_t *dd, ivec tric_dir,
4562 t_state *state, rvec **f,
4571 int ncg[DIM*2], nat[DIM*2];
4572 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4573 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4574 int sbuf[2], rbuf[2];
4575 int home_pos_cg, home_pos_at, buf_pos;
4577 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4580 real inv_ncg, pos_d;
4582 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4584 cginfo_mb_t *cginfo_mb;
4585 gmx_domdec_comm_t *comm;
4587 int nthread, thread;
4591 check_screw_box(state->box);
4595 if (fr->cutoff_scheme == ecutsGROUP)
4600 for (i = 0; i < estNR; i++)
4606 case estX: /* Always present */ break;
4607 case estV: bV = (state->flags & (1<<i)); break;
4608 case estSDX: bSDX = (state->flags & (1<<i)); break;
4609 case estCGP: bCGP = (state->flags & (1<<i)); break;
4612 case estDISRE_INITF:
4613 case estDISRE_RM3TAV:
4614 case estORIRE_INITF:
4616 /* No processing required */
4619 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4624 if (dd->ncg_tot > comm->nalloc_int)
4626 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4627 srenew(comm->buf_int, comm->nalloc_int);
4629 move = comm->buf_int;
4631 /* Clear the count */
4632 for (c = 0; c < dd->ndim*2; c++)
4638 npbcdim = dd->npbcdim;
4640 for (d = 0; (d < DIM); d++)
4642 limitd[d] = dd->comm->cellsize_min[d];
4643 if (d >= npbcdim && dd->ci[d] == 0)
4645 cell_x0[d] = -GMX_FLOAT_MAX;
4649 cell_x0[d] = comm->cell_x0[d];
4651 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4653 cell_x1[d] = GMX_FLOAT_MAX;
4657 cell_x1[d] = comm->cell_x1[d];
4661 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4662 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4666 /* We check after communication if a charge group moved
4667 * more than one cell. Set the pre-comm check limit to float_max.
4669 limit0[d] = -GMX_FLOAT_MAX;
4670 limit1[d] = GMX_FLOAT_MAX;
4674 make_tric_corr_matrix(npbcdim, state->box, tcm);
4676 cgindex = dd->cgindex;
4678 nthread = gmx_omp_nthreads_get(emntDomdec);
4680 /* Compute the center of geometry for all home charge groups
4681 * and put them in the box and determine where they should go.
4683 #pragma omp parallel for num_threads(nthread) schedule(static)
4684 for (thread = 0; thread < nthread; thread++)
4686 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4687 cell_x0, cell_x1, limitd, limit0, limit1,
4689 ( thread *dd->ncg_home)/nthread,
4690 ((thread+1)*dd->ncg_home)/nthread,
4691 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4695 for (cg = 0; cg < dd->ncg_home; cg++)
4700 flag = mc & ~DD_FLAG_NRCG;
4701 mc = mc & DD_FLAG_NRCG;
4704 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4706 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4707 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4709 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4710 /* We store the cg size in the lower 16 bits
4711 * and the place where the charge group should go
4712 * in the next 6 bits. This saves some communication volume.
4714 nrcg = cgindex[cg+1] - cgindex[cg];
4715 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4721 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4722 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4725 for (i = 0; i < dd->ndim*2; i++)
4727 *ncg_moved += ncg[i];
4744 /* Make sure the communication buffers are large enough */
4745 for (mc = 0; mc < dd->ndim*2; mc++)
4747 nvr = ncg[mc] + nat[mc]*nvec;
4748 if (nvr > comm->cgcm_state_nalloc[mc])
4750 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4751 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4755 switch (fr->cutoff_scheme)
4758 /* Recalculating cg_cm might be cheaper than communicating,
4759 * but that could give rise to rounding issues.
4762 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4763 nvec, cg_cm, comm, bCompact);
4766 /* Without charge groups we send the moved atom coordinates
4767 * over twice. This is so the code below can be used without
4768 * many conditionals for both for with and without charge groups.
4771 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4772 nvec, state->x, comm, FALSE);
4775 home_pos_cg -= *ncg_moved;
4779 gmx_incons("unimplemented");
4785 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4786 nvec, vec++, state->x, comm, bCompact);
4789 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4790 nvec, vec++, state->v, comm, bCompact);
4794 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4795 nvec, vec++, state->sd_X, comm, bCompact);
4799 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4800 nvec, vec++, state->cg_p, comm, bCompact);
4805 compact_ind(dd->ncg_home, move,
4806 dd->index_gl, dd->cgindex, dd->gatindex,
4807 dd->ga2la, comm->bLocalCG,
4812 if (fr->cutoff_scheme == ecutsVERLET)
4814 moved = get_moved(comm, dd->ncg_home);
4816 for (k = 0; k < dd->ncg_home; k++)
4823 moved = fr->ns.grid->cell_index;
4826 clear_and_mark_ind(dd->ncg_home, move,
4827 dd->index_gl, dd->cgindex, dd->gatindex,
4828 dd->ga2la, comm->bLocalCG,
4832 cginfo_mb = fr->cginfo_mb;
4834 *ncg_stay_home = home_pos_cg;
4835 for (d = 0; d < dd->ndim; d++)
4841 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4844 /* Communicate the cg and atom counts */
4849 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4850 d, dir, sbuf[0], sbuf[1]);
4852 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4854 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4856 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4857 srenew(comm->buf_int, comm->nalloc_int);
4860 /* Communicate the charge group indices, sizes and flags */
4861 dd_sendrecv_int(dd, d, dir,
4862 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4863 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4865 nvs = ncg[cdd] + nat[cdd]*nvec;
4866 i = rbuf[0] + rbuf[1] *nvec;
4867 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4869 /* Communicate cgcm and state */
4870 dd_sendrecv_rvec(dd, d, dir,
4871 comm->cgcm_state[cdd], nvs,
4872 comm->vbuf.v+nvr, i);
4873 ncg_recv += rbuf[0];
4874 nat_recv += rbuf[1];
4878 /* Process the received charge groups */
4880 for (cg = 0; cg < ncg_recv; cg++)
4882 flag = comm->buf_int[cg*DD_CGIBS+1];
4884 if (dim >= npbcdim && dd->nc[dim] > 2)
4886 /* No pbc in this dim and more than one domain boundary.
4887 * We do a separate check if a charge group didn't move too far.
4889 if (((flag & DD_FLAG_FW(d)) &&
4890 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4891 ((flag & DD_FLAG_BW(d)) &&
4892 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4894 cg_move_error(fplog, dd, step, cg, dim,
4895 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4897 comm->vbuf.v[buf_pos],
4898 comm->vbuf.v[buf_pos],
4899 comm->vbuf.v[buf_pos][dim]);
4906 /* Check which direction this cg should go */
4907 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4911 /* The cell boundaries for dimension d2 are not equal
4912 * for each cell row of the lower dimension(s),
4913 * therefore we might need to redetermine where
4914 * this cg should go.
4917 /* If this cg crosses the box boundary in dimension d2
4918 * we can use the communicated flag, so we do not
4919 * have to worry about pbc.
4921 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4922 (flag & DD_FLAG_FW(d2))) ||
4923 (dd->ci[dim2] == 0 &&
4924 (flag & DD_FLAG_BW(d2)))))
4926 /* Clear the two flags for this dimension */
4927 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4928 /* Determine the location of this cg
4929 * in lattice coordinates
4931 pos_d = comm->vbuf.v[buf_pos][dim2];
4934 for (d3 = dim2+1; d3 < DIM; d3++)
4937 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4940 /* Check of we are not at the box edge.
4941 * pbc is only handled in the first step above,
4942 * but this check could move over pbc while
4943 * the first step did not due to different rounding.
4945 if (pos_d >= cell_x1[dim2] &&
4946 dd->ci[dim2] != dd->nc[dim2]-1)
4948 flag |= DD_FLAG_FW(d2);
4950 else if (pos_d < cell_x0[dim2] &&
4953 flag |= DD_FLAG_BW(d2);
4955 comm->buf_int[cg*DD_CGIBS+1] = flag;
4958 /* Set to which neighboring cell this cg should go */
4959 if (flag & DD_FLAG_FW(d2))
4963 else if (flag & DD_FLAG_BW(d2))
4965 if (dd->nc[dd->dim[d2]] > 2)
4977 nrcg = flag & DD_FLAG_NRCG;
4980 if (home_pos_cg+1 > dd->cg_nalloc)
4982 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4983 srenew(dd->index_gl, dd->cg_nalloc);
4984 srenew(dd->cgindex, dd->cg_nalloc+1);
4986 /* Set the global charge group index and size */
4987 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4988 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4989 /* Copy the state from the buffer */
4990 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4991 if (fr->cutoff_scheme == ecutsGROUP)
4994 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4998 /* Set the cginfo */
4999 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5000 dd->index_gl[home_pos_cg]);
5003 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5006 if (home_pos_at+nrcg > state->nalloc)
5008 dd_realloc_state(state, f, home_pos_at+nrcg);
5010 for (i = 0; i < nrcg; i++)
5012 copy_rvec(comm->vbuf.v[buf_pos++],
5013 state->x[home_pos_at+i]);
5017 for (i = 0; i < nrcg; i++)
5019 copy_rvec(comm->vbuf.v[buf_pos++],
5020 state->v[home_pos_at+i]);
5025 for (i = 0; i < nrcg; i++)
5027 copy_rvec(comm->vbuf.v[buf_pos++],
5028 state->sd_X[home_pos_at+i]);
5033 for (i = 0; i < nrcg; i++)
5035 copy_rvec(comm->vbuf.v[buf_pos++],
5036 state->cg_p[home_pos_at+i]);
5040 home_pos_at += nrcg;
5044 /* Reallocate the buffers if necessary */
5045 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5047 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5048 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5050 nvr = ncg[mc] + nat[mc]*nvec;
5051 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5053 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5054 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5056 /* Copy from the receive to the send buffers */
5057 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5058 comm->buf_int + cg*DD_CGIBS,
5059 DD_CGIBS*sizeof(int));
5060 memcpy(comm->cgcm_state[mc][nvr],
5061 comm->vbuf.v[buf_pos],
5062 (1+nrcg*nvec)*sizeof(rvec));
5063 buf_pos += 1 + nrcg*nvec;
5070 /* With sorting (!bCompact) the indices are now only partially up to date
5071 * and ncg_home and nat_home are not the real count, since there are
5072 * "holes" in the arrays for the charge groups that moved to neighbors.
5074 if (fr->cutoff_scheme == ecutsVERLET)
5076 moved = get_moved(comm, home_pos_cg);
5078 for (i = dd->ncg_home; i < home_pos_cg; i++)
5083 dd->ncg_home = home_pos_cg;
5084 dd->nat_home = home_pos_at;
5089 "Finished repartitioning: cgs moved out %d, new home %d\n",
5090 *ncg_moved, dd->ncg_home-*ncg_moved);
5095 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5097 dd->comm->cycl[ddCycl] += cycles;
5098 dd->comm->cycl_n[ddCycl]++;
5099 if (cycles > dd->comm->cycl_max[ddCycl])
5101 dd->comm->cycl_max[ddCycl] = cycles;
5105 static double force_flop_count(t_nrnb *nrnb)
5112 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5114 /* To get closer to the real timings, we half the count
5115 * for the normal loops and again half it for water loops.
5118 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5120 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5124 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5127 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5130 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5132 sum += nrnb->n[i]*cost_nrnb(i);
5135 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5137 sum += nrnb->n[i]*cost_nrnb(i);
5143 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5145 if (dd->comm->eFlop)
5147 dd->comm->flop -= force_flop_count(nrnb);
5150 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5152 if (dd->comm->eFlop)
5154 dd->comm->flop += force_flop_count(nrnb);
5159 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5163 for (i = 0; i < ddCyclNr; i++)
5165 dd->comm->cycl[i] = 0;
5166 dd->comm->cycl_n[i] = 0;
5167 dd->comm->cycl_max[i] = 0;
5170 dd->comm->flop_n = 0;
5173 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5175 gmx_domdec_comm_t *comm;
5176 gmx_domdec_load_t *load;
5177 gmx_domdec_root_t *root = NULL;
5178 int d, dim, cid, i, pos;
5179 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5184 fprintf(debug, "get_load_distribution start\n");
5187 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5191 bSepPME = (dd->pme_nodeid >= 0);
5193 for (d = dd->ndim-1; d >= 0; d--)
5196 /* Check if we participate in the communication in this dimension */
5197 if (d == dd->ndim-1 ||
5198 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5200 load = &comm->load[d];
5203 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5206 if (d == dd->ndim-1)
5208 sbuf[pos++] = dd_force_load(comm);
5209 sbuf[pos++] = sbuf[0];
5212 sbuf[pos++] = sbuf[0];
5213 sbuf[pos++] = cell_frac;
5216 sbuf[pos++] = comm->cell_f_max0[d];
5217 sbuf[pos++] = comm->cell_f_min1[d];
5222 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5223 sbuf[pos++] = comm->cycl[ddCyclPME];
5228 sbuf[pos++] = comm->load[d+1].sum;
5229 sbuf[pos++] = comm->load[d+1].max;
5232 sbuf[pos++] = comm->load[d+1].sum_m;
5233 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5234 sbuf[pos++] = comm->load[d+1].flags;
5237 sbuf[pos++] = comm->cell_f_max0[d];
5238 sbuf[pos++] = comm->cell_f_min1[d];
5243 sbuf[pos++] = comm->load[d+1].mdf;
5244 sbuf[pos++] = comm->load[d+1].pme;
5248 /* Communicate a row in DD direction d.
5249 * The communicators are setup such that the root always has rank 0.
5252 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5253 load->load, load->nload*sizeof(float), MPI_BYTE,
5254 0, comm->mpi_comm_load[d]);
5256 if (dd->ci[dim] == dd->master_ci[dim])
5258 /* We are the root, process this row */
5259 if (comm->bDynLoadBal)
5261 root = comm->root[d];
5271 for (i = 0; i < dd->nc[dim]; i++)
5273 load->sum += load->load[pos++];
5274 load->max = max(load->max, load->load[pos]);
5280 /* This direction could not be load balanced properly,
5281 * therefore we need to use the maximum iso the average load.
5283 load->sum_m = max(load->sum_m, load->load[pos]);
5287 load->sum_m += load->load[pos];
5290 load->cvol_min = min(load->cvol_min, load->load[pos]);
5294 load->flags = (int)(load->load[pos++] + 0.5);
5298 root->cell_f_max0[i] = load->load[pos++];
5299 root->cell_f_min1[i] = load->load[pos++];
5304 load->mdf = max(load->mdf, load->load[pos]);
5306 load->pme = max(load->pme, load->load[pos]);
5310 if (comm->bDynLoadBal && root->bLimited)
5312 load->sum_m *= dd->nc[dim];
5313 load->flags |= (1<<d);
5321 comm->nload += dd_load_count(comm);
5322 comm->load_step += comm->cycl[ddCyclStep];
5323 comm->load_sum += comm->load[0].sum;
5324 comm->load_max += comm->load[0].max;
5325 if (comm->bDynLoadBal)
5327 for (d = 0; d < dd->ndim; d++)
5329 if (comm->load[0].flags & (1<<d))
5331 comm->load_lim[d]++;
5337 comm->load_mdf += comm->load[0].mdf;
5338 comm->load_pme += comm->load[0].pme;
5342 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5346 fprintf(debug, "get_load_distribution finished\n");
5350 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5352 /* Return the relative performance loss on the total run time
5353 * due to the force calculation load imbalance.
5355 if (dd->comm->nload > 0)
5358 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5359 (dd->comm->load_step*dd->nnodes);
5367 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5370 int npp, npme, nnodes, d, limp;
5371 float imbal, pme_f_ratio, lossf, lossp = 0;
5373 gmx_domdec_comm_t *comm;
5376 if (DDMASTER(dd) && comm->nload > 0)
5379 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5380 nnodes = npp + npme;
5381 imbal = comm->load_max*npp/comm->load_sum - 1;
5382 lossf = dd_force_imb_perf_loss(dd);
5383 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5384 fprintf(fplog, "%s", buf);
5385 fprintf(stderr, "\n");
5386 fprintf(stderr, "%s", buf);
5387 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5388 fprintf(fplog, "%s", buf);
5389 fprintf(stderr, "%s", buf);
5391 if (comm->bDynLoadBal)
5393 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5394 for (d = 0; d < dd->ndim; d++)
5396 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5397 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5403 sprintf(buf+strlen(buf), "\n");
5404 fprintf(fplog, "%s", buf);
5405 fprintf(stderr, "%s", buf);
5409 pme_f_ratio = comm->load_pme/comm->load_mdf;
5410 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5413 lossp *= (float)npme/(float)nnodes;
5417 lossp *= (float)npp/(float)nnodes;
5419 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5420 fprintf(fplog, "%s", buf);
5421 fprintf(stderr, "%s", buf);
5422 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5423 fprintf(fplog, "%s", buf);
5424 fprintf(stderr, "%s", buf);
5426 fprintf(fplog, "\n");
5427 fprintf(stderr, "\n");
5429 if (lossf >= DD_PERF_LOSS)
5432 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5433 " in the domain decomposition.\n", lossf*100);
5434 if (!comm->bDynLoadBal)
5436 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5440 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5442 fprintf(fplog, "%s\n", buf);
5443 fprintf(stderr, "%s\n", buf);
5445 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5448 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5449 " had %s work to do than the PP nodes.\n"
5450 " You might want to %s the number of PME nodes\n"
5451 " or %s the cut-off and the grid spacing.\n",
5453 (lossp < 0) ? "less" : "more",
5454 (lossp < 0) ? "decrease" : "increase",
5455 (lossp < 0) ? "decrease" : "increase");
5456 fprintf(fplog, "%s\n", buf);
5457 fprintf(stderr, "%s\n", buf);
5462 static float dd_vol_min(gmx_domdec_t *dd)
5464 return dd->comm->load[0].cvol_min*dd->nnodes;
5467 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5469 return dd->comm->load[0].flags;
5472 static float dd_f_imbal(gmx_domdec_t *dd)
5474 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5477 float dd_pme_f_ratio(gmx_domdec_t *dd)
5479 if (dd->comm->cycl_n[ddCyclPME] > 0)
5481 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5489 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5494 flags = dd_load_flags(dd);
5498 "DD load balancing is limited by minimum cell size in dimension");
5499 for (d = 0; d < dd->ndim; d++)
5503 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5506 fprintf(fplog, "\n");
5508 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5509 if (dd->comm->bDynLoadBal)
5511 fprintf(fplog, " vol min/aver %5.3f%c",
5512 dd_vol_min(dd), flags ? '!' : ' ');
5514 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5515 if (dd->comm->cycl_n[ddCyclPME])
5517 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5519 fprintf(fplog, "\n\n");
5522 static void dd_print_load_verbose(gmx_domdec_t *dd)
5524 if (dd->comm->bDynLoadBal)
5526 fprintf(stderr, "vol %4.2f%c ",
5527 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5529 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5530 if (dd->comm->cycl_n[ddCyclPME])
5532 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5537 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5542 gmx_domdec_root_t *root;
5543 gmx_bool bPartOfGroup = FALSE;
5545 dim = dd->dim[dim_ind];
5546 copy_ivec(loc, loc_c);
5547 for (i = 0; i < dd->nc[dim]; i++)
5550 rank = dd_index(dd->nc, loc_c);
5551 if (rank == dd->rank)
5553 /* This process is part of the group */
5554 bPartOfGroup = TRUE;
5557 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5561 dd->comm->mpi_comm_load[dim_ind] = c_row;
5562 if (dd->comm->eDLB != edlbNO)
5564 if (dd->ci[dim] == dd->master_ci[dim])
5566 /* This is the root process of this row */
5567 snew(dd->comm->root[dim_ind], 1);
5568 root = dd->comm->root[dim_ind];
5569 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5570 snew(root->old_cell_f, dd->nc[dim]+1);
5571 snew(root->bCellMin, dd->nc[dim]);
5574 snew(root->cell_f_max0, dd->nc[dim]);
5575 snew(root->cell_f_min1, dd->nc[dim]);
5576 snew(root->bound_min, dd->nc[dim]);
5577 snew(root->bound_max, dd->nc[dim]);
5579 snew(root->buf_ncd, dd->nc[dim]);
5583 /* This is not a root process, we only need to receive cell_f */
5584 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5587 if (dd->ci[dim] == dd->master_ci[dim])
5589 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5595 static void make_load_communicators(gmx_domdec_t *dd)
5598 int dim0, dim1, i, j;
5603 fprintf(debug, "Making load communicators\n");
5606 snew(dd->comm->load, dd->ndim);
5607 snew(dd->comm->mpi_comm_load, dd->ndim);
5610 make_load_communicator(dd, 0, loc);
5614 for (i = 0; i < dd->nc[dim0]; i++)
5617 make_load_communicator(dd, 1, loc);
5623 for (i = 0; i < dd->nc[dim0]; i++)
5627 for (j = 0; j < dd->nc[dim1]; j++)
5630 make_load_communicator(dd, 2, loc);
5637 fprintf(debug, "Finished making load communicators\n");
5642 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5645 int d, dim, i, j, m;
5648 ivec dd_zp[DD_MAXIZONE];
5649 gmx_domdec_zones_t *zones;
5650 gmx_domdec_ns_ranges_t *izone;
5652 for (d = 0; d < dd->ndim; d++)
5655 copy_ivec(dd->ci, tmp);
5656 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5657 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5658 copy_ivec(dd->ci, tmp);
5659 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5660 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5663 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5666 dd->neighbor[d][1]);
5672 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5674 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5675 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5682 for (i = 0; i < nzonep; i++)
5684 copy_ivec(dd_zp3[i], dd_zp[i]);
5690 for (i = 0; i < nzonep; i++)
5692 copy_ivec(dd_zp2[i], dd_zp[i]);
5698 for (i = 0; i < nzonep; i++)
5700 copy_ivec(dd_zp1[i], dd_zp[i]);
5704 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5709 zones = &dd->comm->zones;
5711 for (i = 0; i < nzone; i++)
5714 clear_ivec(zones->shift[i]);
5715 for (d = 0; d < dd->ndim; d++)
5717 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5722 for (i = 0; i < nzone; i++)
5724 for (d = 0; d < DIM; d++)
5726 s[d] = dd->ci[d] - zones->shift[i][d];
5731 else if (s[d] >= dd->nc[d])
5737 zones->nizone = nzonep;
5738 for (i = 0; i < zones->nizone; i++)
5740 if (dd_zp[i][0] != i)
5742 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5744 izone = &zones->izone[i];
5745 izone->j0 = dd_zp[i][1];
5746 izone->j1 = dd_zp[i][2];
5747 for (dim = 0; dim < DIM; dim++)
5749 if (dd->nc[dim] == 1)
5751 /* All shifts should be allowed */
5752 izone->shift0[dim] = -1;
5753 izone->shift1[dim] = 1;
5758 izone->shift0[d] = 0;
5759 izone->shift1[d] = 0;
5760 for(j=izone->j0; j<izone->j1; j++) {
5761 if (dd->shift[j][d] > dd->shift[i][d])
5762 izone->shift0[d] = -1;
5763 if (dd->shift[j][d] < dd->shift[i][d])
5764 izone->shift1[d] = 1;
5770 /* Assume the shift are not more than 1 cell */
5771 izone->shift0[dim] = 1;
5772 izone->shift1[dim] = -1;
5773 for (j = izone->j0; j < izone->j1; j++)
5775 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5776 if (shift_diff < izone->shift0[dim])
5778 izone->shift0[dim] = shift_diff;
5780 if (shift_diff > izone->shift1[dim])
5782 izone->shift1[dim] = shift_diff;
5789 if (dd->comm->eDLB != edlbNO)
5791 snew(dd->comm->root, dd->ndim);
5794 if (dd->comm->bRecordLoad)
5796 make_load_communicators(dd);
5800 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5803 gmx_domdec_comm_t *comm;
5814 if (comm->bCartesianPP)
5816 /* Set up cartesian communication for the particle-particle part */
5819 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5820 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5823 for (i = 0; i < DIM; i++)
5827 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5829 /* We overwrite the old communicator with the new cartesian one */
5830 cr->mpi_comm_mygroup = comm_cart;
5833 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5834 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5836 if (comm->bCartesianPP_PME)
5838 /* Since we want to use the original cartesian setup for sim,
5839 * and not the one after split, we need to make an index.
5841 snew(comm->ddindex2ddnodeid, dd->nnodes);
5842 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5843 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5844 /* Get the rank of the DD master,
5845 * above we made sure that the master node is a PP node.
5855 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5857 else if (comm->bCartesianPP)
5859 if (cr->npmenodes == 0)
5861 /* The PP communicator is also
5862 * the communicator for this simulation
5864 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5866 cr->nodeid = dd->rank;
5868 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5870 /* We need to make an index to go from the coordinates
5871 * to the nodeid of this simulation.
5873 snew(comm->ddindex2simnodeid, dd->nnodes);
5874 snew(buf, dd->nnodes);
5875 if (cr->duty & DUTY_PP)
5877 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5879 /* Communicate the ddindex to simulation nodeid index */
5880 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5881 cr->mpi_comm_mysim);
5884 /* Determine the master coordinates and rank.
5885 * The DD master should be the same node as the master of this sim.
5887 for (i = 0; i < dd->nnodes; i++)
5889 if (comm->ddindex2simnodeid[i] == 0)
5891 ddindex2xyz(dd->nc, i, dd->master_ci);
5892 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5897 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5902 /* No Cartesian communicators */
5903 /* We use the rank in dd->comm->all as DD index */
5904 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5905 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5907 clear_ivec(dd->master_ci);
5914 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5915 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5920 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5921 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5925 static void receive_ddindex2simnodeid(t_commrec *cr)
5929 gmx_domdec_comm_t *comm;
5936 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5938 snew(comm->ddindex2simnodeid, dd->nnodes);
5939 snew(buf, dd->nnodes);
5940 if (cr->duty & DUTY_PP)
5942 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5945 /* Communicate the ddindex to simulation nodeid index */
5946 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5947 cr->mpi_comm_mysim);
5954 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5955 int ncg, int natoms)
5957 gmx_domdec_master_t *ma;
5962 snew(ma->ncg, dd->nnodes);
5963 snew(ma->index, dd->nnodes+1);
5965 snew(ma->nat, dd->nnodes);
5966 snew(ma->ibuf, dd->nnodes*2);
5967 snew(ma->cell_x, DIM);
5968 for (i = 0; i < DIM; i++)
5970 snew(ma->cell_x[i], dd->nc[i]+1);
5973 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5979 snew(ma->vbuf, natoms);
5985 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5989 gmx_domdec_comm_t *comm;
6000 if (comm->bCartesianPP)
6002 for (i = 1; i < DIM; i++)
6004 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6006 if (bDiv[YY] || bDiv[ZZ])
6008 comm->bCartesianPP_PME = TRUE;
6009 /* If we have 2D PME decomposition, which is always in x+y,
6010 * we stack the PME only nodes in z.
6011 * Otherwise we choose the direction that provides the thinnest slab
6012 * of PME only nodes as this will have the least effect
6013 * on the PP communication.
6014 * But for the PME communication the opposite might be better.
6016 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6018 dd->nc[YY] > dd->nc[ZZ]))
6020 comm->cartpmedim = ZZ;
6024 comm->cartpmedim = YY;
6026 comm->ntot[comm->cartpmedim]
6027 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6031 fprintf(fplog, "#pmenodes (%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]);
6033 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6038 if (comm->bCartesianPP_PME)
6042 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]);
6045 for (i = 0; i < DIM; i++)
6049 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6052 MPI_Comm_rank(comm_cart, &rank);
6053 if (MASTERNODE(cr) && rank != 0)
6055 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6058 /* With this assigment we loose the link to the original communicator
6059 * which will usually be MPI_COMM_WORLD, unless have multisim.
6061 cr->mpi_comm_mysim = comm_cart;
6062 cr->sim_nodeid = rank;
6064 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6068 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6069 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6072 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6076 if (cr->npmenodes == 0 ||
6077 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6079 cr->duty = DUTY_PME;
6082 /* Split the sim communicator into PP and PME only nodes */
6083 MPI_Comm_split(cr->mpi_comm_mysim,
6085 dd_index(comm->ntot, dd->ci),
6086 &cr->mpi_comm_mygroup);
6090 switch (dd_node_order)
6095 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6098 case ddnoINTERLEAVE:
6099 /* Interleave the PP-only and PME-only nodes,
6100 * as on clusters with dual-core machines this will double
6101 * the communication bandwidth of the PME processes
6102 * and thus speed up the PP <-> PME and inter PME communication.
6106 fprintf(fplog, "Interleaving PP and PME nodes\n");
6108 comm->pmenodes = dd_pmenodes(cr);
6113 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6116 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6118 cr->duty = DUTY_PME;
6125 /* Split the sim communicator into PP and PME only nodes */
6126 MPI_Comm_split(cr->mpi_comm_mysim,
6129 &cr->mpi_comm_mygroup);
6130 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6136 fprintf(fplog, "This is a %s only node\n\n",
6137 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6141 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6144 gmx_domdec_comm_t *comm;
6150 copy_ivec(dd->nc, comm->ntot);
6152 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6153 comm->bCartesianPP_PME = FALSE;
6155 /* Reorder the nodes by default. This might change the MPI ranks.
6156 * Real reordering is only supported on very few architectures,
6157 * Blue Gene is one of them.
6159 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6161 if (cr->npmenodes > 0)
6163 /* Split the communicator into a PP and PME part */
6164 split_communicator(fplog, cr, dd_node_order, CartReorder);
6165 if (comm->bCartesianPP_PME)
6167 /* We (possibly) reordered the nodes in split_communicator,
6168 * so it is no longer required in make_pp_communicator.
6170 CartReorder = FALSE;
6175 /* All nodes do PP and PME */
6177 /* We do not require separate communicators */
6178 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6182 if (cr->duty & DUTY_PP)
6184 /* Copy or make a new PP communicator */
6185 make_pp_communicator(fplog, cr, CartReorder);
6189 receive_ddindex2simnodeid(cr);
6192 if (!(cr->duty & DUTY_PME))
6194 /* Set up the commnuication to our PME node */
6195 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6196 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6199 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6200 dd->pme_nodeid, dd->pme_receive_vir_ener);
6205 dd->pme_nodeid = -1;
6210 dd->ma = init_gmx_domdec_master_t(dd,
6212 comm->cgs_gl.index[comm->cgs_gl.nr]);
6216 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6218 real *slb_frac, tot;
6223 if (nc > 1 && size_string != NULL)
6227 fprintf(fplog, "Using static load balancing for the %s direction\n",
6232 for (i = 0; i < nc; i++)
6235 sscanf(size_string, "%lf%n", &dbl, &n);
6238 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6247 fprintf(fplog, "Relative cell sizes:");
6249 for (i = 0; i < nc; i++)
6254 fprintf(fplog, " %5.3f", slb_frac[i]);
6259 fprintf(fplog, "\n");
6266 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6269 gmx_mtop_ilistloop_t iloop;
6273 iloop = gmx_mtop_ilistloop_init(mtop);
6274 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6276 for (ftype = 0; ftype < F_NRE; ftype++)
6278 if ((interaction_function[ftype].flags & IF_BOND) &&
6281 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6289 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6295 val = getenv(env_var);
6298 if (sscanf(val, "%d", &nst) <= 0)
6304 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6312 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6316 fprintf(stderr, "\n%s\n", warn_string);
6320 fprintf(fplog, "\n%s\n", warn_string);
6324 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6325 t_inputrec *ir, FILE *fplog)
6327 if (ir->ePBC == epbcSCREW &&
6328 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6330 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6333 if (ir->ns_type == ensSIMPLE)
6335 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6338 if (ir->nstlist == 0)
6340 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6343 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6345 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6349 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6354 r = ddbox->box_size[XX];
6355 for (di = 0; di < dd->ndim; di++)
6358 /* Check using the initial average cell size */
6359 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6365 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6366 const char *dlb_opt, gmx_bool bRecordLoad,
6367 unsigned long Flags, t_inputrec *ir)
6375 case 'a': eDLB = edlbAUTO; break;
6376 case 'n': eDLB = edlbNO; break;
6377 case 'y': eDLB = edlbYES; break;
6378 default: gmx_incons("Unknown dlb_opt");
6381 if (Flags & MD_RERUN)
6386 if (!EI_DYNAMICS(ir->eI))
6388 if (eDLB == edlbYES)
6390 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6391 dd_warning(cr, fplog, buf);
6399 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6404 if (Flags & MD_REPRODUCIBLE)
6411 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6415 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6418 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6426 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6431 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6433 /* Decomposition order z,y,x */
6436 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6438 for (dim = DIM-1; dim >= 0; dim--)
6440 if (dd->nc[dim] > 1)
6442 dd->dim[dd->ndim++] = dim;
6448 /* Decomposition order x,y,z */
6449 for (dim = 0; dim < DIM; dim++)
6451 if (dd->nc[dim] > 1)
6453 dd->dim[dd->ndim++] = dim;
6459 static gmx_domdec_comm_t *init_dd_comm()
6461 gmx_domdec_comm_t *comm;
6465 snew(comm->cggl_flag, DIM*2);
6466 snew(comm->cgcm_state, DIM*2);
6467 for (i = 0; i < DIM*2; i++)
6469 comm->cggl_flag_nalloc[i] = 0;
6470 comm->cgcm_state_nalloc[i] = 0;
6473 comm->nalloc_int = 0;
6474 comm->buf_int = NULL;
6476 vec_rvec_init(&comm->vbuf);
6478 comm->n_load_have = 0;
6479 comm->n_load_collect = 0;
6481 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6483 comm->sum_nat[i] = 0;
6487 comm->load_step = 0;
6490 clear_ivec(comm->load_lim);
6497 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6498 unsigned long Flags,
6500 real comm_distance_min, real rconstr,
6501 const char *dlb_opt, real dlb_scale,
6502 const char *sizex, const char *sizey, const char *sizez,
6503 gmx_mtop_t *mtop, t_inputrec *ir,
6504 matrix box, rvec *x,
6506 int *npme_x, int *npme_y)
6509 gmx_domdec_comm_t *comm;
6512 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6519 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6524 dd->comm = init_dd_comm();
6526 snew(comm->cggl_flag, DIM*2);
6527 snew(comm->cgcm_state, DIM*2);
6529 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6530 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6532 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6533 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6534 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6535 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6536 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6537 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6538 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6539 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6541 dd->pme_recv_f_alloc = 0;
6542 dd->pme_recv_f_buf = NULL;
6544 if (dd->bSendRecv2 && fplog)
6546 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");
6552 fprintf(fplog, "Will load balance based on FLOP count\n");
6554 if (comm->eFlop > 1)
6556 srand(1+cr->nodeid);
6558 comm->bRecordLoad = TRUE;
6562 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6566 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6568 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6571 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6573 dd->bGridJump = comm->bDynLoadBal;
6574 comm->bPMELoadBalDLBLimits = FALSE;
6576 if (comm->nstSortCG)
6580 if (comm->nstSortCG == 1)
6582 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6586 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6590 snew(comm->sort, 1);
6596 fprintf(fplog, "Will not sort the charge groups\n");
6600 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6602 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6603 if (comm->bInterCGBondeds)
6605 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6609 comm->bInterCGMultiBody = FALSE;
6612 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6613 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6615 if (ir->rlistlong == 0)
6617 /* Set the cut-off to some very large value,
6618 * so we don't need if statements everywhere in the code.
6619 * We use sqrt, since the cut-off is squared in some places.
6621 comm->cutoff = GMX_CUTOFF_INF;
6625 comm->cutoff = ir->rlistlong;
6627 comm->cutoff_mbody = 0;
6629 comm->cellsize_limit = 0;
6630 comm->bBondComm = FALSE;
6632 if (comm->bInterCGBondeds)
6634 if (comm_distance_min > 0)
6636 comm->cutoff_mbody = comm_distance_min;
6637 if (Flags & MD_DDBONDCOMM)
6639 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6643 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6645 r_bonded_limit = comm->cutoff_mbody;
6647 else if (ir->bPeriodicMols)
6649 /* Can not easily determine the required cut-off */
6650 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");
6651 comm->cutoff_mbody = comm->cutoff/2;
6652 r_bonded_limit = comm->cutoff_mbody;
6658 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6659 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6661 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6662 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6664 /* We use an initial margin of 10% for the minimum cell size,
6665 * except when we are just below the non-bonded cut-off.
6667 if (Flags & MD_DDBONDCOMM)
6669 if (max(r_2b, r_mb) > comm->cutoff)
6671 r_bonded = max(r_2b, r_mb);
6672 r_bonded_limit = 1.1*r_bonded;
6673 comm->bBondComm = TRUE;
6678 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6680 /* We determine cutoff_mbody later */
6684 /* No special bonded communication,
6685 * simply increase the DD cut-off.
6687 r_bonded_limit = 1.1*max(r_2b, r_mb);
6688 comm->cutoff_mbody = r_bonded_limit;
6689 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6692 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6696 "Minimum cell size due to bonded interactions: %.3f nm\n",
6697 comm->cellsize_limit);
6701 if (dd->bInterCGcons && rconstr <= 0)
6703 /* There is a cell size limit due to the constraints (P-LINCS) */
6704 rconstr = constr_r_max(fplog, mtop, ir);
6708 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6710 if (rconstr > comm->cellsize_limit)
6712 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6716 else if (rconstr > 0 && fplog)
6718 /* Here we do not check for dd->bInterCGcons,
6719 * because one can also set a cell size limit for virtual sites only
6720 * and at this point we don't know yet if there are intercg v-sites.
6723 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6726 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6728 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6732 copy_ivec(nc, dd->nc);
6733 set_dd_dim(fplog, dd);
6734 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6736 if (cr->npmenodes == -1)
6740 acs = average_cellsize_min(dd, ddbox);
6741 if (acs < comm->cellsize_limit)
6745 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6747 gmx_fatal_collective(FARGS, cr, NULL,
6748 "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",
6749 acs, comm->cellsize_limit);
6754 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6756 /* We need to choose the optimal DD grid and possibly PME nodes */
6757 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6758 comm->eDLB != edlbNO, dlb_scale,
6759 comm->cellsize_limit, comm->cutoff,
6760 comm->bInterCGBondeds);
6762 if (dd->nc[XX] == 0)
6764 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6765 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6766 !bC ? "-rdd" : "-rcon",
6767 comm->eDLB != edlbNO ? " or -dds" : "",
6768 bC ? " or your LINCS settings" : "");
6770 gmx_fatal_collective(FARGS, cr, NULL,
6771 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6773 "Look in the log file for details on the domain decomposition",
6774 cr->nnodes-cr->npmenodes, limit, buf);
6776 set_dd_dim(fplog, dd);
6782 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6783 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6786 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6787 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6789 gmx_fatal_collective(FARGS, cr, NULL,
6790 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6791 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6793 if (cr->npmenodes > dd->nnodes)
6795 gmx_fatal_collective(FARGS, cr, NULL,
6796 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6798 if (cr->npmenodes > 0)
6800 comm->npmenodes = cr->npmenodes;
6804 comm->npmenodes = dd->nnodes;
6807 if (EEL_PME(ir->coulombtype))
6809 /* The following choices should match those
6810 * in comm_cost_est in domdec_setup.c.
6811 * Note that here the checks have to take into account
6812 * that the decomposition might occur in a different order than xyz
6813 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6814 * in which case they will not match those in comm_cost_est,
6815 * but since that is mainly for testing purposes that's fine.
6817 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6818 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6819 getenv("GMX_PMEONEDD") == NULL)
6821 comm->npmedecompdim = 2;
6822 comm->npmenodes_x = dd->nc[XX];
6823 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6827 /* In case nc is 1 in both x and y we could still choose to
6828 * decompose pme in y instead of x, but we use x for simplicity.
6830 comm->npmedecompdim = 1;
6831 if (dd->dim[0] == YY)
6833 comm->npmenodes_x = 1;
6834 comm->npmenodes_y = comm->npmenodes;
6838 comm->npmenodes_x = comm->npmenodes;
6839 comm->npmenodes_y = 1;
6844 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6845 comm->npmenodes_x, comm->npmenodes_y, 1);
6850 comm->npmedecompdim = 0;
6851 comm->npmenodes_x = 0;
6852 comm->npmenodes_y = 0;
6855 /* Technically we don't need both of these,
6856 * but it simplifies code not having to recalculate it.
6858 *npme_x = comm->npmenodes_x;
6859 *npme_y = comm->npmenodes_y;
6861 snew(comm->slb_frac, DIM);
6862 if (comm->eDLB == edlbNO)
6864 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6865 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6866 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6869 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6871 if (comm->bBondComm || comm->eDLB != edlbNO)
6873 /* Set the bonded communication distance to halfway
6874 * the minimum and the maximum,
6875 * since the extra communication cost is nearly zero.
6877 acs = average_cellsize_min(dd, ddbox);
6878 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6879 if (comm->eDLB != edlbNO)
6881 /* Check if this does not limit the scaling */
6882 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6884 if (!comm->bBondComm)
6886 /* Without bBondComm do not go beyond the n.b. cut-off */
6887 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6888 if (comm->cellsize_limit >= comm->cutoff)
6890 /* We don't loose a lot of efficieny
6891 * when increasing it to the n.b. cut-off.
6892 * It can even be slightly faster, because we need
6893 * less checks for the communication setup.
6895 comm->cutoff_mbody = comm->cutoff;
6898 /* Check if we did not end up below our original limit */
6899 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6901 if (comm->cutoff_mbody > comm->cellsize_limit)
6903 comm->cellsize_limit = comm->cutoff_mbody;
6906 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6911 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6912 "cellsize limit %f\n",
6913 comm->bBondComm, comm->cellsize_limit);
6918 check_dd_restrictions(cr, dd, ir, fplog);
6921 comm->partition_step = INT_MIN;
6924 clear_dd_cycle_counts(dd);
6929 static void set_dlb_limits(gmx_domdec_t *dd)
6934 for (d = 0; d < dd->ndim; d++)
6936 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6937 dd->comm->cellsize_min[dd->dim[d]] =
6938 dd->comm->cellsize_min_dlb[dd->dim[d]];
6943 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6946 gmx_domdec_comm_t *comm;
6956 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);
6959 cellsize_min = comm->cellsize_min[dd->dim[0]];
6960 for (d = 1; d < dd->ndim; d++)
6962 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6965 if (cellsize_min < comm->cellsize_limit*1.05)
6967 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");
6969 /* Change DLB from "auto" to "no". */
6970 comm->eDLB = edlbNO;
6975 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6976 comm->bDynLoadBal = TRUE;
6977 dd->bGridJump = TRUE;
6981 /* We can set the required cell size info here,
6982 * so we do not need to communicate this.
6983 * The grid is completely uniform.
6985 for (d = 0; d < dd->ndim; d++)
6989 comm->load[d].sum_m = comm->load[d].sum;
6991 nc = dd->nc[dd->dim[d]];
6992 for (i = 0; i < nc; i++)
6994 comm->root[d]->cell_f[i] = i/(real)nc;
6997 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6998 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7001 comm->root[d]->cell_f[nc] = 1.0;
7006 static char *init_bLocalCG(gmx_mtop_t *mtop)
7011 ncg = ncg_mtop(mtop);
7012 snew(bLocalCG, ncg);
7013 for (cg = 0; cg < ncg; cg++)
7015 bLocalCG[cg] = FALSE;
7021 void dd_init_bondeds(FILE *fplog,
7022 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7024 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7026 gmx_domdec_comm_t *comm;
7030 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7034 if (comm->bBondComm)
7036 /* Communicate atoms beyond the cut-off for bonded interactions */
7039 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7041 comm->bLocalCG = init_bLocalCG(mtop);
7045 /* Only communicate atoms based on cut-off */
7046 comm->cglink = NULL;
7047 comm->bLocalCG = NULL;
7051 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7053 gmx_bool bDynLoadBal, real dlb_scale,
7056 gmx_domdec_comm_t *comm;
7071 fprintf(fplog, "The maximum number of communication pulses is:");
7072 for (d = 0; d < dd->ndim; d++)
7074 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7076 fprintf(fplog, "\n");
7077 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7078 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7079 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7080 for (d = 0; d < DIM; d++)
7084 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7091 comm->cellsize_min_dlb[d]/
7092 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7094 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7097 fprintf(fplog, "\n");
7101 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7102 fprintf(fplog, "The initial number of communication pulses is:");
7103 for (d = 0; d < dd->ndim; d++)
7105 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7107 fprintf(fplog, "\n");
7108 fprintf(fplog, "The initial domain decomposition cell size is:");
7109 for (d = 0; d < DIM; d++)
7113 fprintf(fplog, " %c %.2f nm",
7114 dim2char(d), dd->comm->cellsize_min[d]);
7117 fprintf(fplog, "\n\n");
7120 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7122 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7123 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7124 "non-bonded interactions", "", comm->cutoff);
7128 limit = dd->comm->cellsize_limit;
7132 if (dynamic_dd_box(ddbox, ir))
7134 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7136 limit = dd->comm->cellsize_min[XX];
7137 for (d = 1; d < DIM; d++)
7139 limit = min(limit, dd->comm->cellsize_min[d]);
7143 if (comm->bInterCGBondeds)
7145 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7146 "two-body bonded interactions", "(-rdd)",
7147 max(comm->cutoff, comm->cutoff_mbody));
7148 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7149 "multi-body bonded interactions", "(-rdd)",
7150 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7154 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7155 "virtual site constructions", "(-rcon)", limit);
7157 if (dd->constraint_comm)
7159 sprintf(buf, "atoms separated by up to %d constraints",
7161 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7162 buf, "(-rcon)", limit);
7164 fprintf(fplog, "\n");
7170 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7172 const t_inputrec *ir,
7173 const gmx_ddbox_t *ddbox)
7175 gmx_domdec_comm_t *comm;
7176 int d, dim, npulse, npulse_d_max, npulse_d;
7181 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7183 /* Determine the maximum number of comm. pulses in one dimension */
7185 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7187 /* Determine the maximum required number of grid pulses */
7188 if (comm->cellsize_limit >= comm->cutoff)
7190 /* Only a single pulse is required */
7193 else if (!bNoCutOff && comm->cellsize_limit > 0)
7195 /* We round down slightly here to avoid overhead due to the latency
7196 * of extra communication calls when the cut-off
7197 * would be only slightly longer than the cell size.
7198 * Later cellsize_limit is redetermined,
7199 * so we can not miss interactions due to this rounding.
7201 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7205 /* There is no cell size limit */
7206 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7209 if (!bNoCutOff && npulse > 1)
7211 /* See if we can do with less pulses, based on dlb_scale */
7213 for (d = 0; d < dd->ndim; d++)
7216 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7217 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7218 npulse_d_max = max(npulse_d_max, npulse_d);
7220 npulse = min(npulse, npulse_d_max);
7223 /* This env var can override npulse */
7224 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7231 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7232 for (d = 0; d < dd->ndim; d++)
7234 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7235 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7236 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7237 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7238 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7240 comm->bVacDLBNoLimit = FALSE;
7244 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7245 if (!comm->bVacDLBNoLimit)
7247 comm->cellsize_limit = max(comm->cellsize_limit,
7248 comm->cutoff/comm->maxpulse);
7250 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7251 /* Set the minimum cell size for each DD dimension */
7252 for (d = 0; d < dd->ndim; d++)
7254 if (comm->bVacDLBNoLimit ||
7255 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7257 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7261 comm->cellsize_min_dlb[dd->dim[d]] =
7262 comm->cutoff/comm->cd[d].np_dlb;
7265 if (comm->cutoff_mbody <= 0)
7267 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7269 if (comm->bDynLoadBal)
7275 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7277 /* If each molecule is a single charge group
7278 * or we use domain decomposition for each periodic dimension,
7279 * we do not need to take pbc into account for the bonded interactions.
7281 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7284 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7287 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7288 t_inputrec *ir, gmx_ddbox_t *ddbox)
7290 gmx_domdec_comm_t *comm;
7296 /* Initialize the thread data.
7297 * This can not be done in init_domain_decomposition,
7298 * as the numbers of threads is determined later.
7300 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7303 snew(comm->dth, comm->nth);
7306 if (EEL_PME(ir->coulombtype))
7308 init_ddpme(dd, &comm->ddpme[0], 0);
7309 if (comm->npmedecompdim >= 2)
7311 init_ddpme(dd, &comm->ddpme[1], 1);
7316 comm->npmenodes = 0;
7317 if (dd->pme_nodeid >= 0)
7319 gmx_fatal_collective(FARGS, NULL, dd,
7320 "Can not have separate PME nodes without PME electrostatics");
7326 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7328 if (comm->eDLB != edlbNO)
7330 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7333 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7334 if (comm->eDLB == edlbAUTO)
7338 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7340 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7343 if (ir->ePBC == epbcNONE)
7345 vol_frac = 1 - 1/(double)dd->nnodes;
7350 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7354 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7356 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7358 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7361 static gmx_bool test_dd_cutoff(t_commrec *cr,
7362 t_state *state, t_inputrec *ir,
7373 set_ddbox(dd, FALSE, cr, ir, state->box,
7374 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7378 for (d = 0; d < dd->ndim; d++)
7382 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7383 if (dynamic_dd_box(&ddbox, ir))
7385 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7388 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7390 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7391 dd->comm->cd[d].np_dlb > 0)
7393 if (np > dd->comm->cd[d].np_dlb)
7398 /* If a current local cell size is smaller than the requested
7399 * cut-off, we could still fix it, but this gets very complicated.
7400 * Without fixing here, we might actually need more checks.
7402 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7409 if (dd->comm->eDLB != edlbNO)
7411 /* If DLB is not active yet, we don't need to check the grid jumps.
7412 * Actually we shouldn't, because then the grid jump data is not set.
7414 if (dd->comm->bDynLoadBal &&
7415 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7420 gmx_sumi(1, &LocallyLimited, cr);
7422 if (LocallyLimited > 0)
7431 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7434 gmx_bool bCutoffAllowed;
7436 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7440 cr->dd->comm->cutoff = cutoff_req;
7443 return bCutoffAllowed;
7446 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7448 gmx_domdec_comm_t *comm;
7450 comm = cr->dd->comm;
7452 /* Turn on the DLB limiting (might have been on already) */
7453 comm->bPMELoadBalDLBLimits = TRUE;
7455 /* Change the cut-off limit */
7456 comm->PMELoadBal_max_cutoff = comm->cutoff;
7459 static void merge_cg_buffers(int ncell,
7460 gmx_domdec_comm_dim_t *cd, int pulse,
7462 int *index_gl, int *recv_i,
7463 rvec *cg_cm, rvec *recv_vr,
7465 cginfo_mb_t *cginfo_mb, int *cginfo)
7467 gmx_domdec_ind_t *ind, *ind_p;
7468 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7469 int shift, shift_at;
7471 ind = &cd->ind[pulse];
7473 /* First correct the already stored data */
7474 shift = ind->nrecv[ncell];
7475 for (cell = ncell-1; cell >= 0; cell--)
7477 shift -= ind->nrecv[cell];
7480 /* Move the cg's present from previous grid pulses */
7481 cg0 = ncg_cell[ncell+cell];
7482 cg1 = ncg_cell[ncell+cell+1];
7483 cgindex[cg1+shift] = cgindex[cg1];
7484 for (cg = cg1-1; cg >= cg0; cg--)
7486 index_gl[cg+shift] = index_gl[cg];
7487 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7488 cgindex[cg+shift] = cgindex[cg];
7489 cginfo[cg+shift] = cginfo[cg];
7491 /* Correct the already stored send indices for the shift */
7492 for (p = 1; p <= pulse; p++)
7494 ind_p = &cd->ind[p];
7496 for (c = 0; c < cell; c++)
7498 cg0 += ind_p->nsend[c];
7500 cg1 = cg0 + ind_p->nsend[cell];
7501 for (cg = cg0; cg < cg1; cg++)
7503 ind_p->index[cg] += shift;
7509 /* Merge in the communicated buffers */
7513 for (cell = 0; cell < ncell; cell++)
7515 cg1 = ncg_cell[ncell+cell+1] + shift;
7518 /* Correct the old cg indices */
7519 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7521 cgindex[cg+1] += shift_at;
7524 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7526 /* Copy this charge group from the buffer */
7527 index_gl[cg1] = recv_i[cg0];
7528 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7529 /* Add it to the cgindex */
7530 cg_gl = index_gl[cg1];
7531 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7532 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7533 cgindex[cg1+1] = cgindex[cg1] + nat;
7538 shift += ind->nrecv[cell];
7539 ncg_cell[ncell+cell+1] = cg1;
7543 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7544 int nzone, int cg0, const int *cgindex)
7548 /* Store the atom block boundaries for easy copying of communication buffers
7551 for (zone = 0; zone < nzone; zone++)
7553 for (p = 0; p < cd->np; p++)
7555 cd->ind[p].cell2at0[zone] = cgindex[cg];
7556 cg += cd->ind[p].nrecv[zone];
7557 cd->ind[p].cell2at1[zone] = cgindex[cg];
7562 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7568 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7570 if (!bLocalCG[link->a[i]])
7579 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7581 real c[DIM][4]; /* the corners for the non-bonded communication */
7582 real cr0; /* corner for rounding */
7583 real cr1[4]; /* corners for rounding */
7584 real bc[DIM]; /* corners for bounded communication */
7585 real bcr1; /* corner for rounding for bonded communication */
7588 /* Determine the corners of the domain(s) we are communicating with */
7590 set_dd_corners(const gmx_domdec_t *dd,
7591 int dim0, int dim1, int dim2,
7595 const gmx_domdec_comm_t *comm;
7596 const gmx_domdec_zones_t *zones;
7601 zones = &comm->zones;
7603 /* Keep the compiler happy */
7607 /* The first dimension is equal for all cells */
7608 c->c[0][0] = comm->cell_x0[dim0];
7611 c->bc[0] = c->c[0][0];
7616 /* This cell row is only seen from the first row */
7617 c->c[1][0] = comm->cell_x0[dim1];
7618 /* All rows can see this row */
7619 c->c[1][1] = comm->cell_x0[dim1];
7622 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7625 /* For the multi-body distance we need the maximum */
7626 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7629 /* Set the upper-right corner for rounding */
7630 c->cr0 = comm->cell_x1[dim0];
7635 for (j = 0; j < 4; j++)
7637 c->c[2][j] = comm->cell_x0[dim2];
7641 /* Use the maximum of the i-cells that see a j-cell */
7642 for (i = 0; i < zones->nizone; i++)
7644 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7650 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7656 /* For the multi-body distance we need the maximum */
7657 c->bc[2] = comm->cell_x0[dim2];
7658 for (i = 0; i < 2; i++)
7660 for (j = 0; j < 2; j++)
7662 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7668 /* Set the upper-right corner for rounding */
7669 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7670 * Only cell (0,0,0) can see cell 7 (1,1,1)
7672 c->cr1[0] = comm->cell_x1[dim1];
7673 c->cr1[3] = comm->cell_x1[dim1];
7676 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7679 /* For the multi-body distance we need the maximum */
7680 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7687 /* Determine which cg's we need to send in this pulse from this zone */
7689 get_zone_pulse_cgs(gmx_domdec_t *dd,
7690 int zonei, int zone,
7692 const int *index_gl,
7694 int dim, int dim_ind,
7695 int dim0, int dim1, int dim2,
7696 real r_comm2, real r_bcomm2,
7700 real skew_fac2_d, real skew_fac_01,
7701 rvec *v_d, rvec *v_0, rvec *v_1,
7702 const dd_corners_t *c,
7704 gmx_bool bDistBonded,
7710 gmx_domdec_ind_t *ind,
7711 int **ibuf, int *ibuf_nalloc,
7717 gmx_domdec_comm_t *comm;
7719 gmx_bool bDistMB_pulse;
7721 real r2, rb2, r, tric_sh;
7724 int nsend_z, nsend, nat;
7728 bScrew = (dd->bScrewPBC && dim == XX);
7730 bDistMB_pulse = (bDistMB && bDistBonded);
7736 for (cg = cg0; cg < cg1; cg++)
7740 if (tric_dist[dim_ind] == 0)
7742 /* Rectangular direction, easy */
7743 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7750 r = cg_cm[cg][dim] - c->bc[dim_ind];
7756 /* Rounding gives at most a 16% reduction
7757 * in communicated atoms
7759 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7761 r = cg_cm[cg][dim0] - c->cr0;
7762 /* This is the first dimension, so always r >= 0 */
7769 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7771 r = cg_cm[cg][dim1] - c->cr1[zone];
7778 r = cg_cm[cg][dim1] - c->bcr1;
7788 /* Triclinic direction, more complicated */
7791 /* Rounding, conservative as the skew_fac multiplication
7792 * will slightly underestimate the distance.
7794 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7796 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7797 for (i = dim0+1; i < DIM; i++)
7799 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7801 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7804 rb[dim0] = rn[dim0];
7807 /* Take care that the cell planes along dim0 might not
7808 * be orthogonal to those along dim1 and dim2.
7810 for (i = 1; i <= dim_ind; i++)
7813 if (normal[dim0][dimd] > 0)
7815 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7818 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7823 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7825 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7827 for (i = dim1+1; i < DIM; i++)
7829 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7831 rn[dim1] += tric_sh;
7834 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7835 /* Take care of coupling of the distances
7836 * to the planes along dim0 and dim1 through dim2.
7838 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7839 /* Take care that the cell planes along dim1
7840 * might not be orthogonal to that along dim2.
7842 if (normal[dim1][dim2] > 0)
7844 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7850 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7853 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7854 /* Take care of coupling of the distances
7855 * to the planes along dim0 and dim1 through dim2.
7857 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7858 /* Take care that the cell planes along dim1
7859 * might not be orthogonal to that along dim2.
7861 if (normal[dim1][dim2] > 0)
7863 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7868 /* The distance along the communication direction */
7869 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7871 for (i = dim+1; i < DIM; i++)
7873 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7878 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7879 /* Take care of coupling of the distances
7880 * to the planes along dim0 and dim1 through dim2.
7882 if (dim_ind == 1 && zonei == 1)
7884 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7890 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7893 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7894 /* Take care of coupling of the distances
7895 * to the planes along dim0 and dim1 through dim2.
7897 if (dim_ind == 1 && zonei == 1)
7899 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7907 ((bDistMB && rb2 < r_bcomm2) ||
7908 (bDist2B && r2 < r_bcomm2)) &&
7910 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7911 missing_link(comm->cglink, index_gl[cg],
7914 /* Make an index to the local charge groups */
7915 if (nsend+1 > ind->nalloc)
7917 ind->nalloc = over_alloc_large(nsend+1);
7918 srenew(ind->index, ind->nalloc);
7920 if (nsend+1 > *ibuf_nalloc)
7922 *ibuf_nalloc = over_alloc_large(nsend+1);
7923 srenew(*ibuf, *ibuf_nalloc);
7925 ind->index[nsend] = cg;
7926 (*ibuf)[nsend] = index_gl[cg];
7928 vec_rvec_check_alloc(vbuf, nsend+1);
7930 if (dd->ci[dim] == 0)
7932 /* Correct cg_cm for pbc */
7933 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7936 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7937 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7942 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7945 nat += cgindex[cg+1] - cgindex[cg];
7951 *nsend_z_ptr = nsend_z;
7954 static void setup_dd_communication(gmx_domdec_t *dd,
7955 matrix box, gmx_ddbox_t *ddbox,
7956 t_forcerec *fr, t_state *state, rvec **f)
7958 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7959 int nzone, nzone_send, zone, zonei, cg0, cg1;
7960 int c, i, j, cg, cg_gl, nrcg;
7961 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7962 gmx_domdec_comm_t *comm;
7963 gmx_domdec_zones_t *zones;
7964 gmx_domdec_comm_dim_t *cd;
7965 gmx_domdec_ind_t *ind;
7966 cginfo_mb_t *cginfo_mb;
7967 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7968 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7969 dd_corners_t corners;
7971 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7972 real skew_fac2_d, skew_fac_01;
7979 fprintf(debug, "Setting up DD communication\n");
7984 switch (fr->cutoff_scheme)
7993 gmx_incons("unimplemented");
7997 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7999 dim = dd->dim[dim_ind];
8001 /* Check if we need to use triclinic distances */
8002 tric_dist[dim_ind] = 0;
8003 for (i = 0; i <= dim_ind; i++)
8005 if (ddbox->tric_dir[dd->dim[i]])
8007 tric_dist[dim_ind] = 1;
8012 bBondComm = comm->bBondComm;
8014 /* Do we need to determine extra distances for multi-body bondeds? */
8015 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8017 /* Do we need to determine extra distances for only two-body bondeds? */
8018 bDist2B = (bBondComm && !bDistMB);
8020 r_comm2 = sqr(comm->cutoff);
8021 r_bcomm2 = sqr(comm->cutoff_mbody);
8025 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8028 zones = &comm->zones;
8031 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8032 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8034 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8036 /* Triclinic stuff */
8037 normal = ddbox->normal;
8041 v_0 = ddbox->v[dim0];
8042 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8044 /* Determine the coupling coefficient for the distances
8045 * to the cell planes along dim0 and dim1 through dim2.
8046 * This is required for correct rounding.
8049 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8052 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8058 v_1 = ddbox->v[dim1];
8061 zone_cg_range = zones->cg_range;
8062 index_gl = dd->index_gl;
8063 cgindex = dd->cgindex;
8064 cginfo_mb = fr->cginfo_mb;
8066 zone_cg_range[0] = 0;
8067 zone_cg_range[1] = dd->ncg_home;
8068 comm->zone_ncg1[0] = dd->ncg_home;
8069 pos_cg = dd->ncg_home;
8071 nat_tot = dd->nat_home;
8073 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8075 dim = dd->dim[dim_ind];
8076 cd = &comm->cd[dim_ind];
8078 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8080 /* No pbc in this dimension, the first node should not comm. */
8088 v_d = ddbox->v[dim];
8089 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8091 cd->bInPlace = TRUE;
8092 for (p = 0; p < cd->np; p++)
8094 /* Only atoms communicated in the first pulse are used
8095 * for multi-body bonded interactions or for bBondComm.
8097 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8102 for (zone = 0; zone < nzone_send; zone++)
8104 if (tric_dist[dim_ind] && dim_ind > 0)
8106 /* Determine slightly more optimized skew_fac's
8108 * This reduces the number of communicated atoms
8109 * by about 10% for 3D DD of rhombic dodecahedra.
8111 for (dimd = 0; dimd < dim; dimd++)
8113 sf2_round[dimd] = 1;
8114 if (ddbox->tric_dir[dimd])
8116 for (i = dd->dim[dimd]+1; i < DIM; i++)
8118 /* If we are shifted in dimension i
8119 * and the cell plane is tilted forward
8120 * in dimension i, skip this coupling.
8122 if (!(zones->shift[nzone+zone][i] &&
8123 ddbox->v[dimd][i][dimd] >= 0))
8126 sqr(ddbox->v[dimd][i][dimd]);
8129 sf2_round[dimd] = 1/sf2_round[dimd];
8134 zonei = zone_perm[dim_ind][zone];
8137 /* Here we permutate the zones to obtain a convenient order
8138 * for neighbor searching
8140 cg0 = zone_cg_range[zonei];
8141 cg1 = zone_cg_range[zonei+1];
8145 /* Look only at the cg's received in the previous grid pulse
8147 cg1 = zone_cg_range[nzone+zone+1];
8148 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8151 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8152 for (th = 0; th < comm->nth; th++)
8154 gmx_domdec_ind_t *ind_p;
8155 int **ibuf_p, *ibuf_nalloc_p;
8157 int *nsend_p, *nat_p;
8163 /* Thread 0 writes in the comm buffers */
8165 ibuf_p = &comm->buf_int;
8166 ibuf_nalloc_p = &comm->nalloc_int;
8167 vbuf_p = &comm->vbuf;
8170 nsend_zone_p = &ind->nsend[zone];
8174 /* Other threads write into temp buffers */
8175 ind_p = &comm->dth[th].ind;
8176 ibuf_p = &comm->dth[th].ibuf;
8177 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8178 vbuf_p = &comm->dth[th].vbuf;
8179 nsend_p = &comm->dth[th].nsend;
8180 nat_p = &comm->dth[th].nat;
8181 nsend_zone_p = &comm->dth[th].nsend_zone;
8183 comm->dth[th].nsend = 0;
8184 comm->dth[th].nat = 0;
8185 comm->dth[th].nsend_zone = 0;
8195 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8196 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8199 /* Get the cg's for this pulse in this zone */
8200 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8202 dim, dim_ind, dim0, dim1, dim2,
8205 normal, skew_fac2_d, skew_fac_01,
8206 v_d, v_0, v_1, &corners, sf2_round,
8207 bDistBonded, bBondComm,
8211 ibuf_p, ibuf_nalloc_p,
8217 /* Append data of threads>=1 to the communication buffers */
8218 for (th = 1; th < comm->nth; th++)
8220 dd_comm_setup_work_t *dth;
8223 dth = &comm->dth[th];
8225 ns1 = nsend + dth->nsend_zone;
8226 if (ns1 > ind->nalloc)
8228 ind->nalloc = over_alloc_dd(ns1);
8229 srenew(ind->index, ind->nalloc);
8231 if (ns1 > comm->nalloc_int)
8233 comm->nalloc_int = over_alloc_dd(ns1);
8234 srenew(comm->buf_int, comm->nalloc_int);
8236 if (ns1 > comm->vbuf.nalloc)
8238 comm->vbuf.nalloc = over_alloc_dd(ns1);
8239 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8242 for (i = 0; i < dth->nsend_zone; i++)
8244 ind->index[nsend] = dth->ind.index[i];
8245 comm->buf_int[nsend] = dth->ibuf[i];
8246 copy_rvec(dth->vbuf.v[i],
8247 comm->vbuf.v[nsend]);
8251 ind->nsend[zone] += dth->nsend_zone;
8254 /* Clear the counts in case we do not have pbc */
8255 for (zone = nzone_send; zone < nzone; zone++)
8257 ind->nsend[zone] = 0;
8259 ind->nsend[nzone] = nsend;
8260 ind->nsend[nzone+1] = nat;
8261 /* Communicate the number of cg's and atoms to receive */
8262 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8263 ind->nsend, nzone+2,
8264 ind->nrecv, nzone+2);
8266 /* The rvec buffer is also required for atom buffers of size nsend
8267 * in dd_move_x and dd_move_f.
8269 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8273 /* We can receive in place if only the last zone is not empty */
8274 for (zone = 0; zone < nzone-1; zone++)
8276 if (ind->nrecv[zone] > 0)
8278 cd->bInPlace = FALSE;
8283 /* The int buffer is only required here for the cg indices */
8284 if (ind->nrecv[nzone] > comm->nalloc_int2)
8286 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8287 srenew(comm->buf_int2, comm->nalloc_int2);
8289 /* The rvec buffer is also required for atom buffers
8290 * of size nrecv in dd_move_x and dd_move_f.
8292 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8293 vec_rvec_check_alloc(&comm->vbuf2, i);
8297 /* Make space for the global cg indices */
8298 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8299 || dd->cg_nalloc == 0)
8301 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8302 srenew(index_gl, dd->cg_nalloc);
8303 srenew(cgindex, dd->cg_nalloc+1);
8305 /* Communicate the global cg indices */
8308 recv_i = index_gl + pos_cg;
8312 recv_i = comm->buf_int2;
8314 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8315 comm->buf_int, nsend,
8316 recv_i, ind->nrecv[nzone]);
8318 /* Make space for cg_cm */
8319 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8320 if (fr->cutoff_scheme == ecutsGROUP)
8328 /* Communicate cg_cm */
8331 recv_vr = cg_cm + pos_cg;
8335 recv_vr = comm->vbuf2.v;
8337 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8338 comm->vbuf.v, nsend,
8339 recv_vr, ind->nrecv[nzone]);
8341 /* Make the charge group index */
8344 zone = (p == 0 ? 0 : nzone - 1);
8345 while (zone < nzone)
8347 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8349 cg_gl = index_gl[pos_cg];
8350 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8351 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8352 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8355 /* Update the charge group presence,
8356 * so we can use it in the next pass of the loop.
8358 comm->bLocalCG[cg_gl] = TRUE;
8364 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8367 zone_cg_range[nzone+zone] = pos_cg;
8372 /* This part of the code is never executed with bBondComm. */
8373 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8374 index_gl, recv_i, cg_cm, recv_vr,
8375 cgindex, fr->cginfo_mb, fr->cginfo);
8376 pos_cg += ind->nrecv[nzone];
8378 nat_tot += ind->nrecv[nzone+1];
8382 /* Store the atom block for easy copying of communication buffers */
8383 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8387 dd->index_gl = index_gl;
8388 dd->cgindex = cgindex;
8390 dd->ncg_tot = zone_cg_range[zones->n];
8391 dd->nat_tot = nat_tot;
8392 comm->nat[ddnatHOME] = dd->nat_home;
8393 for (i = ddnatZONE; i < ddnatNR; i++)
8395 comm->nat[i] = dd->nat_tot;
8400 /* We don't need to update cginfo, since that was alrady done above.
8401 * So we pass NULL for the forcerec.
8403 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8404 NULL, comm->bLocalCG);
8409 fprintf(debug, "Finished setting up DD communication, zones:");
8410 for (c = 0; c < zones->n; c++)
8412 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8414 fprintf(debug, "\n");
8418 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8422 for (c = 0; c < zones->nizone; c++)
8424 zones->izone[c].cg1 = zones->cg_range[c+1];
8425 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8426 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8430 static void set_zones_size(gmx_domdec_t *dd,
8431 matrix box, const gmx_ddbox_t *ddbox,
8432 int zone_start, int zone_end)
8434 gmx_domdec_comm_t *comm;
8435 gmx_domdec_zones_t *zones;
8437 int z, zi, zj0, zj1, d, dim;
8440 real size_j, add_tric;
8445 zones = &comm->zones;
8447 /* Do we need to determine extra distances for multi-body bondeds? */
8448 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8450 for (z = zone_start; z < zone_end; z++)
8452 /* Copy cell limits to zone limits.
8453 * Valid for non-DD dims and non-shifted dims.
8455 copy_rvec(comm->cell_x0, zones->size[z].x0);
8456 copy_rvec(comm->cell_x1, zones->size[z].x1);
8459 for (d = 0; d < dd->ndim; d++)
8463 for (z = 0; z < zones->n; z++)
8465 /* With a staggered grid we have different sizes
8466 * for non-shifted dimensions.
8468 if (dd->bGridJump && zones->shift[z][dim] == 0)
8472 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8473 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8477 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8478 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8484 rcmbs = comm->cutoff_mbody;
8485 if (ddbox->tric_dir[dim])
8487 rcs /= ddbox->skew_fac[dim];
8488 rcmbs /= ddbox->skew_fac[dim];
8491 /* Set the lower limit for the shifted zone dimensions */
8492 for (z = zone_start; z < zone_end; z++)
8494 if (zones->shift[z][dim] > 0)
8497 if (!dd->bGridJump || d == 0)
8499 zones->size[z].x0[dim] = comm->cell_x1[dim];
8500 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8504 /* Here we take the lower limit of the zone from
8505 * the lowest domain of the zone below.
8509 zones->size[z].x0[dim] =
8510 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8516 zones->size[z].x0[dim] =
8517 zones->size[zone_perm[2][z-4]].x0[dim];
8521 zones->size[z].x0[dim] =
8522 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8525 /* A temporary limit, is updated below */
8526 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8530 for (zi = 0; zi < zones->nizone; zi++)
8532 if (zones->shift[zi][dim] == 0)
8534 /* This takes the whole zone into account.
8535 * With multiple pulses this will lead
8536 * to a larger zone then strictly necessary.
8538 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8539 zones->size[zi].x1[dim]+rcmbs);
8547 /* Loop over the i-zones to set the upper limit of each
8550 for (zi = 0; zi < zones->nizone; zi++)
8552 if (zones->shift[zi][dim] == 0)
8554 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8556 if (zones->shift[z][dim] > 0)
8558 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8559 zones->size[zi].x1[dim]+rcs);
8566 for (z = zone_start; z < zone_end; z++)
8568 /* Initialization only required to keep the compiler happy */
8569 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8572 /* To determine the bounding box for a zone we need to find
8573 * the extreme corners of 4, 2 or 1 corners.
8575 nc = 1 << (ddbox->npbcdim - 1);
8577 for (c = 0; c < nc; c++)
8579 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8583 corner[YY] = zones->size[z].x0[YY];
8587 corner[YY] = zones->size[z].x1[YY];
8591 corner[ZZ] = zones->size[z].x0[ZZ];
8595 corner[ZZ] = zones->size[z].x1[ZZ];
8597 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8599 /* With 1D domain decomposition the cg's are not in
8600 * the triclinic box, but triclinic x-y and rectangular y-z.
8601 * Shift y back, so it will later end up at 0.
8603 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8605 /* Apply the triclinic couplings */
8606 for (i = YY; i < ddbox->npbcdim; i++)
8608 for (j = XX; j < i; j++)
8610 corner[j] += corner[i]*box[i][j]/box[i][i];
8615 copy_rvec(corner, corner_min);
8616 copy_rvec(corner, corner_max);
8620 for (i = 0; i < DIM; i++)
8622 corner_min[i] = min(corner_min[i], corner[i]);
8623 corner_max[i] = max(corner_max[i], corner[i]);
8627 /* Copy the extreme cornes without offset along x */
8628 for (i = 0; i < DIM; i++)
8630 zones->size[z].bb_x0[i] = corner_min[i];
8631 zones->size[z].bb_x1[i] = corner_max[i];
8633 /* Add the offset along x */
8634 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8635 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8638 if (zone_start == 0)
8641 for (dim = 0; dim < DIM; dim++)
8643 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8645 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8650 for (z = zone_start; z < zone_end; z++)
8652 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8654 zones->size[z].x0[XX], zones->size[z].x1[XX],
8655 zones->size[z].x0[YY], zones->size[z].x1[YY],
8656 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8657 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8659 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8660 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8661 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8666 static int comp_cgsort(const void *a, const void *b)
8670 gmx_cgsort_t *cga, *cgb;
8671 cga = (gmx_cgsort_t *)a;
8672 cgb = (gmx_cgsort_t *)b;
8674 comp = cga->nsc - cgb->nsc;
8677 comp = cga->ind_gl - cgb->ind_gl;
8683 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8688 /* Order the data */
8689 for (i = 0; i < n; i++)
8691 buf[i] = a[sort[i].ind];
8694 /* Copy back to the original array */
8695 for (i = 0; i < n; i++)
8701 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8706 /* Order the data */
8707 for (i = 0; i < n; i++)
8709 copy_rvec(v[sort[i].ind], buf[i]);
8712 /* Copy back to the original array */
8713 for (i = 0; i < n; i++)
8715 copy_rvec(buf[i], v[i]);
8719 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8722 int a, atot, cg, cg0, cg1, i;
8724 if (cgindex == NULL)
8726 /* Avoid the useless loop of the atoms within a cg */
8727 order_vec_cg(ncg, sort, v, buf);
8732 /* Order the data */
8734 for (cg = 0; cg < ncg; cg++)
8736 cg0 = cgindex[sort[cg].ind];
8737 cg1 = cgindex[sort[cg].ind+1];
8738 for (i = cg0; i < cg1; i++)
8740 copy_rvec(v[i], buf[a]);
8746 /* Copy back to the original array */
8747 for (a = 0; a < atot; a++)
8749 copy_rvec(buf[a], v[a]);
8753 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8754 int nsort_new, gmx_cgsort_t *sort_new,
8755 gmx_cgsort_t *sort1)
8759 /* The new indices are not very ordered, so we qsort them */
8760 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8762 /* sort2 is already ordered, so now we can merge the two arrays */
8766 while (i2 < nsort2 || i_new < nsort_new)
8770 sort1[i1++] = sort_new[i_new++];
8772 else if (i_new == nsort_new)
8774 sort1[i1++] = sort2[i2++];
8776 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8777 (sort2[i2].nsc == sort_new[i_new].nsc &&
8778 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8780 sort1[i1++] = sort2[i2++];
8784 sort1[i1++] = sort_new[i_new++];
8789 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8791 gmx_domdec_sort_t *sort;
8792 gmx_cgsort_t *cgsort, *sort_i;
8793 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8794 int sort_last, sort_skip;
8796 sort = dd->comm->sort;
8798 a = fr->ns.grid->cell_index;
8800 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8802 if (ncg_home_old >= 0)
8804 /* The charge groups that remained in the same ns grid cell
8805 * are completely ordered. So we can sort efficiently by sorting
8806 * the charge groups that did move into the stationary list.
8811 for (i = 0; i < dd->ncg_home; i++)
8813 /* Check if this cg did not move to another node */
8816 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8818 /* This cg is new on this node or moved ns grid cell */
8819 if (nsort_new >= sort->sort_new_nalloc)
8821 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8822 srenew(sort->sort_new, sort->sort_new_nalloc);
8824 sort_i = &(sort->sort_new[nsort_new++]);
8828 /* This cg did not move */
8829 sort_i = &(sort->sort2[nsort2++]);
8831 /* Sort on the ns grid cell indices
8832 * and the global topology index.
8833 * index_gl is irrelevant with cell ns,
8834 * but we set it here anyhow to avoid a conditional.
8837 sort_i->ind_gl = dd->index_gl[i];
8844 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8847 /* Sort efficiently */
8848 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8853 cgsort = sort->sort;
8855 for (i = 0; i < dd->ncg_home; i++)
8857 /* Sort on the ns grid cell indices
8858 * and the global topology index
8860 cgsort[i].nsc = a[i];
8861 cgsort[i].ind_gl = dd->index_gl[i];
8863 if (cgsort[i].nsc < moved)
8870 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8872 /* Determine the order of the charge groups using qsort */
8873 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8879 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8882 int ncg_new, i, *a, na;
8884 sort = dd->comm->sort->sort;
8886 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8889 for (i = 0; i < na; i++)
8893 sort[ncg_new].ind = a[i];
8901 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
8904 gmx_domdec_sort_t *sort;
8905 gmx_cgsort_t *cgsort, *sort_i;
8907 int ncg_new, i, *ibuf, cgsize;
8910 sort = dd->comm->sort;
8912 if (dd->ncg_home > sort->sort_nalloc)
8914 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8915 srenew(sort->sort, sort->sort_nalloc);
8916 srenew(sort->sort2, sort->sort_nalloc);
8918 cgsort = sort->sort;
8920 switch (fr->cutoff_scheme)
8923 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8926 ncg_new = dd_sort_order_nbnxn(dd, fr);
8929 gmx_incons("unimplemented");
8933 /* We alloc with the old size, since cgindex is still old */
8934 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8935 vbuf = dd->comm->vbuf.v;
8939 cgindex = dd->cgindex;
8946 /* Remove the charge groups which are no longer at home here */
8947 dd->ncg_home = ncg_new;
8950 fprintf(debug, "Set the new home charge group count to %d\n",
8954 /* Reorder the state */
8955 for (i = 0; i < estNR; i++)
8957 if (EST_DISTR(i) && (state->flags & (1<<i)))
8962 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8965 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8968 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8971 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8975 case estDISRE_INITF:
8976 case estDISRE_RM3TAV:
8977 case estORIRE_INITF:
8979 /* No ordering required */
8982 gmx_incons("Unknown state entry encountered in dd_sort_state");
8987 if (fr->cutoff_scheme == ecutsGROUP)
8990 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8993 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8995 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8996 srenew(sort->ibuf, sort->ibuf_nalloc);
8999 /* Reorder the global cg index */
9000 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9001 /* Reorder the cginfo */
9002 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9003 /* Rebuild the local cg index */
9007 for (i = 0; i < dd->ncg_home; i++)
9009 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9010 ibuf[i+1] = ibuf[i] + cgsize;
9012 for (i = 0; i < dd->ncg_home+1; i++)
9014 dd->cgindex[i] = ibuf[i];
9019 for (i = 0; i < dd->ncg_home+1; i++)
9024 /* Set the home atom number */
9025 dd->nat_home = dd->cgindex[dd->ncg_home];
9027 if (fr->cutoff_scheme == ecutsVERLET)
9029 /* The atoms are now exactly in grid order, update the grid order */
9030 nbnxn_set_atomorder(fr->nbv->nbs);
9034 /* Copy the sorted ns cell indices back to the ns grid struct */
9035 for (i = 0; i < dd->ncg_home; i++)
9037 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9039 fr->ns.grid->nr = dd->ncg_home;
9043 static void add_dd_statistics(gmx_domdec_t *dd)
9045 gmx_domdec_comm_t *comm;
9050 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9052 comm->sum_nat[ddnat-ddnatZONE] +=
9053 comm->nat[ddnat] - comm->nat[ddnat-1];
9058 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9060 gmx_domdec_comm_t *comm;
9065 /* Reset all the statistics and counters for total run counting */
9066 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9068 comm->sum_nat[ddnat-ddnatZONE] = 0;
9072 comm->load_step = 0;
9075 clear_ivec(comm->load_lim);
9080 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9082 gmx_domdec_comm_t *comm;
9086 comm = cr->dd->comm;
9088 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9095 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");
9097 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9099 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9104 " av. #atoms communicated per step for force: %d x %.1f\n",
9108 if (cr->dd->vsite_comm)
9111 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9112 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9117 if (cr->dd->constraint_comm)
9120 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9121 1 + ir->nLincsIter, av);
9125 gmx_incons(" Unknown type for DD statistics");
9128 fprintf(fplog, "\n");
9130 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9132 print_dd_load_av(fplog, cr->dd);
9136 void dd_partition_system(FILE *fplog,
9137 gmx_large_int_t step,
9139 gmx_bool bMasterState,
9141 t_state *state_global,
9142 gmx_mtop_t *top_global,
9144 t_state *state_local,
9147 gmx_localtop_t *top_local,
9150 gmx_shellfc_t shellfc,
9151 gmx_constr_t constr,
9153 gmx_wallcycle_t wcycle,
9157 gmx_domdec_comm_t *comm;
9158 gmx_ddbox_t ddbox = {0};
9160 gmx_large_int_t step_pcoupl;
9161 rvec cell_ns_x0, cell_ns_x1;
9162 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9163 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9164 gmx_bool bRedist, bSortCG, bResortAll;
9165 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9172 bBoxChanged = (bMasterState || DEFORM(*ir));
9173 if (ir->epc != epcNO)
9175 /* With nstpcouple > 1 pressure coupling happens.
9176 * one step after calculating the pressure.
9177 * Box scaling happens at the end of the MD step,
9178 * after the DD partitioning.
9179 * We therefore have to do DLB in the first partitioning
9180 * after an MD step where P-coupling occured.
9181 * We need to determine the last step in which p-coupling occurred.
9182 * MRS -- need to validate this for vv?
9187 step_pcoupl = step - 1;
9191 step_pcoupl = ((step - 1)/n)*n + 1;
9193 if (step_pcoupl >= comm->partition_step)
9199 bNStGlobalComm = (step % nstglobalcomm == 0);
9201 if (!comm->bDynLoadBal)
9207 /* Should we do dynamic load balacing this step?
9208 * Since it requires (possibly expensive) global communication,
9209 * we might want to do DLB less frequently.
9211 if (bBoxChanged || ir->epc != epcNO)
9213 bDoDLB = bBoxChanged;
9217 bDoDLB = bNStGlobalComm;
9221 /* Check if we have recorded loads on the nodes */
9222 if (comm->bRecordLoad && dd_load_count(comm))
9224 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9226 /* Check if we should use DLB at the second partitioning
9227 * and every 100 partitionings,
9228 * so the extra communication cost is negligible.
9230 n = max(100, nstglobalcomm);
9231 bCheckDLB = (comm->n_load_collect == 0 ||
9232 comm->n_load_have % n == n-1);
9239 /* Print load every nstlog, first and last step to the log file */
9240 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9241 comm->n_load_collect == 0 ||
9243 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9245 /* Avoid extra communication due to verbose screen output
9246 * when nstglobalcomm is set.
9248 if (bDoDLB || bLogLoad || bCheckDLB ||
9249 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9251 get_load_distribution(dd, wcycle);
9256 dd_print_load(fplog, dd, step-1);
9260 dd_print_load_verbose(dd);
9263 comm->n_load_collect++;
9267 /* Since the timings are node dependent, the master decides */
9271 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9274 fprintf(debug, "step %s, imb loss %f\n",
9275 gmx_step_str(step, sbuf),
9276 dd_force_imb_perf_loss(dd));
9279 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9282 turn_on_dlb(fplog, cr, step);
9287 comm->n_load_have++;
9290 cgs_gl = &comm->cgs_gl;
9295 /* Clear the old state */
9296 clear_dd_indices(dd, 0, 0);
9299 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9300 TRUE, cgs_gl, state_global->x, &ddbox);
9302 get_cg_distribution(fplog, step, dd, cgs_gl,
9303 state_global->box, &ddbox, state_global->x);
9305 dd_distribute_state(dd, cgs_gl,
9306 state_global, state_local, f);
9308 dd_make_local_cgs(dd, &top_local->cgs);
9310 /* Ensure that we have space for the new distribution */
9311 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9313 if (fr->cutoff_scheme == ecutsGROUP)
9315 calc_cgcm(fplog, 0, dd->ncg_home,
9316 &top_local->cgs, state_local->x, fr->cg_cm);
9319 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9321 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9323 else if (state_local->ddp_count != dd->ddp_count)
9325 if (state_local->ddp_count > dd->ddp_count)
9327 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9330 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9332 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);
9335 /* Clear the old state */
9336 clear_dd_indices(dd, 0, 0);
9338 /* Build the new indices */
9339 rebuild_cgindex(dd, cgs_gl->index, state_local);
9340 make_dd_indices(dd, cgs_gl->index, 0);
9341 ncgindex_set = dd->ncg_home;
9343 if (fr->cutoff_scheme == ecutsGROUP)
9345 /* Redetermine the cg COMs */
9346 calc_cgcm(fplog, 0, dd->ncg_home,
9347 &top_local->cgs, state_local->x, fr->cg_cm);
9350 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9352 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9354 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9355 TRUE, &top_local->cgs, state_local->x, &ddbox);
9357 bRedist = comm->bDynLoadBal;
9361 /* We have the full state, only redistribute the cgs */
9363 /* Clear the non-home indices */
9364 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9367 /* Avoid global communication for dim's without pbc and -gcom */
9368 if (!bNStGlobalComm)
9370 copy_rvec(comm->box0, ddbox.box0 );
9371 copy_rvec(comm->box_size, ddbox.box_size);
9373 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9374 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9379 /* For dim's without pbc and -gcom */
9380 copy_rvec(ddbox.box0, comm->box0 );
9381 copy_rvec(ddbox.box_size, comm->box_size);
9383 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9386 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9388 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9391 /* Check if we should sort the charge groups */
9392 if (comm->nstSortCG > 0)
9394 bSortCG = (bMasterState ||
9395 (bRedist && (step % comm->nstSortCG == 0)));
9402 ncg_home_old = dd->ncg_home;
9407 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9409 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9411 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9413 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9416 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9418 &comm->cell_x0, &comm->cell_x1,
9419 dd->ncg_home, fr->cg_cm,
9420 cell_ns_x0, cell_ns_x1, &grid_density);
9424 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9427 switch (fr->cutoff_scheme)
9430 copy_ivec(fr->ns.grid->n, ncells_old);
9431 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9432 state_local->box, cell_ns_x0, cell_ns_x1,
9433 fr->rlistlong, grid_density);
9436 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9439 gmx_incons("unimplemented");
9441 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9442 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9446 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9448 /* Sort the state on charge group position.
9449 * This enables exact restarts from this step.
9450 * It also improves performance by about 15% with larger numbers
9451 * of atoms per node.
9454 /* Fill the ns grid with the home cell,
9455 * so we can sort with the indices.
9457 set_zones_ncg_home(dd);
9459 switch (fr->cutoff_scheme)
9462 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9464 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9466 comm->zones.size[0].bb_x0,
9467 comm->zones.size[0].bb_x1,
9469 comm->zones.dens_zone0,
9472 ncg_moved, bRedist ? comm->moved : NULL,
9473 fr->nbv->grp[eintLocal].kernel_type,
9474 fr->nbv->grp[eintLocal].nbat);
9476 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9479 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9480 0, dd->ncg_home, fr->cg_cm);
9482 copy_ivec(fr->ns.grid->n, ncells_new);
9485 gmx_incons("unimplemented");
9488 bResortAll = bMasterState;
9490 /* Check if we can user the old order and ns grid cell indices
9491 * of the charge groups to sort the charge groups efficiently.
9493 if (ncells_new[XX] != ncells_old[XX] ||
9494 ncells_new[YY] != ncells_old[YY] ||
9495 ncells_new[ZZ] != ncells_old[ZZ])
9502 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9503 gmx_step_str(step, sbuf), dd->ncg_home);
9505 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9506 bResortAll ? -1 : ncg_home_old);
9507 /* Rebuild all the indices */
9508 ga2la_clear(dd->ga2la);
9511 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9514 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9516 /* Setup up the communication and communicate the coordinates */
9517 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9519 /* Set the indices */
9520 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9522 /* Set the charge group boundaries for neighbor searching */
9523 set_cg_boundaries(&comm->zones);
9525 if (fr->cutoff_scheme == ecutsVERLET)
9527 set_zones_size(dd, state_local->box, &ddbox,
9528 bSortCG ? 1 : 0, comm->zones.n);
9531 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9534 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9535 -1,state_local->x,state_local->box);
9538 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9540 /* Extract a local topology from the global topology */
9541 for (i = 0; i < dd->ndim; i++)
9543 np[dd->dim[i]] = comm->cd[i].np;
9545 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9546 comm->cellsize_min, np,
9548 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9549 vsite, top_global, top_local);
9551 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9553 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9555 /* Set up the special atom communication */
9556 n = comm->nat[ddnatZONE];
9557 for (i = ddnatZONE+1; i < ddnatNR; i++)
9562 if (vsite && vsite->n_intercg_vsite)
9564 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9568 if (dd->bInterCGcons || dd->bInterCGsettles)
9570 /* Only for inter-cg constraints we need special code */
9571 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9572 constr, ir->nProjOrder,
9573 top_local->idef.il);
9577 gmx_incons("Unknown special atom type setup");
9582 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9584 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9586 /* Make space for the extra coordinates for virtual site
9587 * or constraint communication.
9589 state_local->natoms = comm->nat[ddnatNR-1];
9590 if (state_local->natoms > state_local->nalloc)
9592 dd_realloc_state(state_local, f, state_local->natoms);
9595 if (fr->bF_NoVirSum)
9597 if (vsite && vsite->n_intercg_vsite)
9599 nat_f_novirsum = comm->nat[ddnatVSITE];
9603 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9605 nat_f_novirsum = dd->nat_tot;
9609 nat_f_novirsum = dd->nat_home;
9618 /* Set the number of atoms required for the force calculation.
9619 * Forces need to be constrained when using a twin-range setup
9620 * or with energy minimization. For simple simulations we could
9621 * avoid some allocation, zeroing and copying, but this is
9622 * probably not worth the complications ande checking.
9624 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9625 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9627 /* We make the all mdatoms up to nat_tot_con.
9628 * We could save some work by only setting invmass
9629 * between nat_tot and nat_tot_con.
9631 /* This call also sets the new number of home particles to dd->nat_home */
9632 atoms2md(top_global, ir,
9633 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9635 /* Now we have the charges we can sort the FE interactions */
9636 dd_sort_local_top(dd, mdatoms, top_local);
9640 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9641 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9646 /* Make the local shell stuff, currently no communication is done */
9647 make_local_shells(cr, mdatoms, shellfc);
9650 if (ir->implicit_solvent)
9652 make_local_gb(cr, fr->born, ir->gb_algorithm);
9655 init_bonded_thread_force_reduction(fr, &top_local->idef);
9657 if (!(cr->duty & DUTY_PME))
9659 /* Send the charges to our PME only node */
9660 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9661 mdatoms->chargeA, mdatoms->chargeB,
9662 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9667 set_constraints(constr, top_local, ir, mdatoms, cr);
9670 if (ir->ePull != epullNO)
9672 /* Update the local pull groups */
9673 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9678 /* Update the local rotation groups */
9679 dd_make_local_rotation_groups(dd, ir->rot);
9683 add_dd_statistics(dd);
9685 /* Make sure we only count the cycles for this DD partitioning */
9686 clear_dd_cycle_counts(dd);
9688 /* Because the order of the atoms might have changed since
9689 * the last vsite construction, we need to communicate the constructing
9690 * atom coordinates again (for spreading the forces this MD step).
9692 dd_move_x_vsites(dd, state_local->box, state_local->x);
9694 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9696 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9698 dd_move_x(dd, state_local->box, state_local->x);
9699 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9700 -1, state_local->x, state_local->box);
9703 /* Store the partitioning step */
9704 comm->partition_step = step;
9706 /* Increase the DD partitioning counter */
9708 /* The state currently matches this DD partitioning count, store it */
9709 state_local->ddp_count = dd->ddp_count;
9712 /* The DD master node knows the complete cg distribution,
9713 * store the count so we can possibly skip the cg info communication.
9715 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9718 if (comm->DD_debug > 0)
9720 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9721 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9722 "after partitioning");