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;
2442 static void rebuild_cgindex(gmx_domdec_t *dd,
2443 const int *gcgs_index, t_state *state)
2445 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2448 dd_cg_gl = dd->index_gl;
2449 cgindex = dd->cgindex;
2452 for (i = 0; i < state->ncg_gl; i++)
2456 dd_cg_gl[i] = cg_gl;
2457 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2461 dd->ncg_home = state->ncg_gl;
2464 set_zones_ncg_home(dd);
2467 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2469 while (cg >= cginfo_mb->cg_end)
2474 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2477 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2478 t_forcerec *fr, char *bLocalCG)
2480 cginfo_mb_t *cginfo_mb;
2486 cginfo_mb = fr->cginfo_mb;
2487 cginfo = fr->cginfo;
2489 for (cg = cg0; cg < cg1; cg++)
2491 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2495 if (bLocalCG != NULL)
2497 for (cg = cg0; cg < cg1; cg++)
2499 bLocalCG[index_gl[cg]] = TRUE;
2504 static void make_dd_indices(gmx_domdec_t *dd,
2505 const int *gcgs_index, int cg_start)
2507 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2508 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2513 bLocalCG = dd->comm->bLocalCG;
2515 if (dd->nat_tot > dd->gatindex_nalloc)
2517 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2518 srenew(dd->gatindex, dd->gatindex_nalloc);
2521 nzone = dd->comm->zones.n;
2522 zone2cg = dd->comm->zones.cg_range;
2523 zone_ncg1 = dd->comm->zone_ncg1;
2524 index_gl = dd->index_gl;
2525 gatindex = dd->gatindex;
2526 bCGs = dd->comm->bCGs;
2528 if (zone2cg[1] != dd->ncg_home)
2530 gmx_incons("dd->ncg_zone is not up to date");
2533 /* Make the local to global and global to local atom index */
2534 a = dd->cgindex[cg_start];
2535 for (zone = 0; zone < nzone; zone++)
2543 cg0 = zone2cg[zone];
2545 cg1 = zone2cg[zone+1];
2546 cg1_p1 = cg0 + zone_ncg1[zone];
2548 for (cg = cg0; cg < cg1; cg++)
2553 /* Signal that this cg is from more than one pulse away */
2556 cg_gl = index_gl[cg];
2559 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2562 ga2la_set(dd->ga2la, a_gl, a, zone1);
2568 gatindex[a] = cg_gl;
2569 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2576 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2579 int ncg, i, ngl, nerr;
2582 if (bLocalCG == NULL)
2586 for (i = 0; i < dd->ncg_tot; i++)
2588 if (!bLocalCG[dd->index_gl[i]])
2591 "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);
2596 for (i = 0; i < ncg_sys; i++)
2603 if (ngl != dd->ncg_tot)
2605 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);
2612 static void check_index_consistency(gmx_domdec_t *dd,
2613 int natoms_sys, int ncg_sys,
2616 int nerr, ngl, i, a, cell;
2621 if (dd->comm->DD_debug > 1)
2623 snew(have, natoms_sys);
2624 for (a = 0; a < dd->nat_tot; a++)
2626 if (have[dd->gatindex[a]] > 0)
2628 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);
2632 have[dd->gatindex[a]] = a + 1;
2638 snew(have, dd->nat_tot);
2641 for (i = 0; i < natoms_sys; i++)
2643 if (ga2la_get(dd->ga2la, i, &a, &cell))
2645 if (a >= dd->nat_tot)
2647 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);
2653 if (dd->gatindex[a] != i)
2655 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);
2662 if (ngl != dd->nat_tot)
2665 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2666 dd->rank, where, ngl, dd->nat_tot);
2668 for (a = 0; a < dd->nat_tot; a++)
2673 "DD node %d, %s: local atom %d, global %d has no global index\n",
2674 dd->rank, where, a+1, dd->gatindex[a]+1);
2679 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2683 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2684 dd->rank, where, nerr);
2688 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2695 /* Clear the whole list without searching */
2696 ga2la_clear(dd->ga2la);
2700 for (i = a_start; i < dd->nat_tot; i++)
2702 ga2la_del(dd->ga2la, dd->gatindex[i]);
2706 bLocalCG = dd->comm->bLocalCG;
2709 for (i = cg_start; i < dd->ncg_tot; i++)
2711 bLocalCG[dd->index_gl[i]] = FALSE;
2715 dd_clear_local_vsite_indices(dd);
2717 if (dd->constraints)
2719 dd_clear_local_constraint_indices(dd);
2723 /* This function should be used for moving the domain boudaries during DLB,
2724 * for obtaining the minimum cell size. It checks the initially set limit
2725 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2726 * and, possibly, a longer cut-off limit set for PME load balancing.
2728 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2732 cellsize_min = comm->cellsize_min[dim];
2734 if (!comm->bVacDLBNoLimit)
2736 /* The cut-off might have changed, e.g. by PME load balacning,
2737 * from the value used to set comm->cellsize_min, so check it.
2739 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2741 if (comm->bPMELoadBalDLBLimits)
2743 /* Check for the cut-off limit set by the PME load balancing */
2744 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2748 return cellsize_min;
2751 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2754 real grid_jump_limit;
2756 /* The distance between the boundaries of cells at distance
2757 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2758 * and by the fact that cells should not be shifted by more than
2759 * half their size, such that cg's only shift by one cell
2760 * at redecomposition.
2762 grid_jump_limit = comm->cellsize_limit;
2763 if (!comm->bVacDLBNoLimit)
2765 if (comm->bPMELoadBalDLBLimits)
2767 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2769 grid_jump_limit = max(grid_jump_limit,
2770 cutoff/comm->cd[dim_ind].np);
2773 return grid_jump_limit;
2776 static gmx_bool check_grid_jump(gmx_large_int_t step,
2782 gmx_domdec_comm_t *comm;
2791 for (d = 1; d < dd->ndim; d++)
2794 limit = grid_jump_limit(comm, cutoff, d);
2795 bfac = ddbox->box_size[dim];
2796 if (ddbox->tric_dir[dim])
2798 bfac *= ddbox->skew_fac[dim];
2800 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2801 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2809 /* This error should never be triggered under normal
2810 * circumstances, but you never know ...
2812 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.",
2813 gmx_step_str(step, buf),
2814 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2822 static int dd_load_count(gmx_domdec_comm_t *comm)
2824 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2827 static float dd_force_load(gmx_domdec_comm_t *comm)
2834 if (comm->eFlop > 1)
2836 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2841 load = comm->cycl[ddCyclF];
2842 if (comm->cycl_n[ddCyclF] > 1)
2844 /* Subtract the maximum of the last n cycle counts
2845 * to get rid of possible high counts due to other soures,
2846 * for instance system activity, that would otherwise
2847 * affect the dynamic load balancing.
2849 load -= comm->cycl_max[ddCyclF];
2856 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2858 gmx_domdec_comm_t *comm;
2863 snew(*dim_f, dd->nc[dim]+1);
2865 for (i = 1; i < dd->nc[dim]; i++)
2867 if (comm->slb_frac[dim])
2869 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2873 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2876 (*dim_f)[dd->nc[dim]] = 1;
2879 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2881 int pmeindex, slab, nso, i;
2884 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2890 ddpme->dim = dimind;
2892 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2894 ddpme->nslab = (ddpme->dim == 0 ?
2895 dd->comm->npmenodes_x :
2896 dd->comm->npmenodes_y);
2898 if (ddpme->nslab <= 1)
2903 nso = dd->comm->npmenodes/ddpme->nslab;
2904 /* Determine for each PME slab the PP location range for dimension dim */
2905 snew(ddpme->pp_min, ddpme->nslab);
2906 snew(ddpme->pp_max, ddpme->nslab);
2907 for (slab = 0; slab < ddpme->nslab; slab++)
2909 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2910 ddpme->pp_max[slab] = 0;
2912 for (i = 0; i < dd->nnodes; i++)
2914 ddindex2xyz(dd->nc, i, xyz);
2915 /* For y only use our y/z slab.
2916 * This assumes that the PME x grid size matches the DD grid size.
2918 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2920 pmeindex = ddindex2pmeindex(dd, i);
2923 slab = pmeindex/nso;
2927 slab = pmeindex % ddpme->nslab;
2929 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2930 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2934 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2937 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2939 if (dd->comm->ddpme[0].dim == XX)
2941 return dd->comm->ddpme[0].maxshift;
2949 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2951 if (dd->comm->ddpme[0].dim == YY)
2953 return dd->comm->ddpme[0].maxshift;
2955 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2957 return dd->comm->ddpme[1].maxshift;
2965 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2966 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2968 gmx_domdec_comm_t *comm;
2971 real range, pme_boundary;
2975 nc = dd->nc[ddpme->dim];
2978 if (!ddpme->dim_match)
2980 /* PP decomposition is not along dim: the worst situation */
2983 else if (ns <= 3 || (bUniform && ns == nc))
2985 /* The optimal situation */
2990 /* We need to check for all pme nodes which nodes they
2991 * could possibly need to communicate with.
2993 xmin = ddpme->pp_min;
2994 xmax = ddpme->pp_max;
2995 /* Allow for atoms to be maximally 2/3 times the cut-off
2996 * out of their DD cell. This is a reasonable balance between
2997 * between performance and support for most charge-group/cut-off
3000 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3001 /* Avoid extra communication when we are exactly at a boundary */
3005 for (s = 0; s < ns; s++)
3007 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3008 pme_boundary = (real)s/ns;
3011 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3013 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3017 pme_boundary = (real)(s+1)/ns;
3020 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3022 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3029 ddpme->maxshift = sh;
3033 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3034 ddpme->dim, ddpme->maxshift);
3038 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3042 for (d = 0; d < dd->ndim; d++)
3045 if (dim < ddbox->nboundeddim &&
3046 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3047 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3049 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",
3050 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3051 dd->nc[dim], dd->comm->cellsize_limit);
3056 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3057 gmx_bool bMaster, ivec npulse)
3059 gmx_domdec_comm_t *comm;
3062 real *cell_x, cell_dx, cellsize;
3066 for (d = 0; d < DIM; d++)
3068 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3070 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3073 cell_dx = ddbox->box_size[d]/dd->nc[d];
3076 for (j = 0; j < dd->nc[d]+1; j++)
3078 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3083 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3084 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3086 cellsize = cell_dx*ddbox->skew_fac[d];
3087 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3091 cellsize_min[d] = cellsize;
3095 /* Statically load balanced grid */
3096 /* Also when we are not doing a master distribution we determine
3097 * all cell borders in a loop to obtain identical values
3098 * to the master distribution case and to determine npulse.
3102 cell_x = dd->ma->cell_x[d];
3106 snew(cell_x, dd->nc[d]+1);
3108 cell_x[0] = ddbox->box0[d];
3109 for (j = 0; j < dd->nc[d]; j++)
3111 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3112 cell_x[j+1] = cell_x[j] + cell_dx;
3113 cellsize = cell_dx*ddbox->skew_fac[d];
3114 while (cellsize*npulse[d] < comm->cutoff &&
3115 npulse[d] < dd->nc[d]-1)
3119 cellsize_min[d] = min(cellsize_min[d], cellsize);
3123 comm->cell_x0[d] = cell_x[dd->ci[d]];
3124 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3128 /* The following limitation is to avoid that a cell would receive
3129 * some of its own home charge groups back over the periodic boundary.
3130 * Double charge groups cause trouble with the global indices.
3132 if (d < ddbox->npbcdim &&
3133 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3135 gmx_fatal_collective(FARGS, NULL, dd,
3136 "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",
3137 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3139 dd->nc[d], dd->nc[d],
3140 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3144 if (!comm->bDynLoadBal)
3146 copy_rvec(cellsize_min, comm->cellsize_min);
3149 for (d = 0; d < comm->npmedecompdim; d++)
3151 set_pme_maxshift(dd, &comm->ddpme[d],
3152 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3153 comm->ddpme[d].slb_dim_f);
3158 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3159 int d, int dim, gmx_domdec_root_t *root,
3161 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3163 gmx_domdec_comm_t *comm;
3164 int ncd, i, j, nmin, nmin_old;
3165 gmx_bool bLimLo, bLimHi;
3167 real fac, halfway, cellsize_limit_f_i, region_size;
3168 gmx_bool bPBC, bLastHi = FALSE;
3169 int nrange[] = {range[0], range[1]};
3171 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3177 bPBC = (dim < ddbox->npbcdim);
3179 cell_size = root->buf_ncd;
3183 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3186 /* First we need to check if the scaling does not make cells
3187 * smaller than the smallest allowed size.
3188 * We need to do this iteratively, since if a cell is too small,
3189 * it needs to be enlarged, which makes all the other cells smaller,
3190 * which could in turn make another cell smaller than allowed.
3192 for (i = range[0]; i < range[1]; i++)
3194 root->bCellMin[i] = FALSE;
3200 /* We need the total for normalization */
3202 for (i = range[0]; i < range[1]; i++)
3204 if (root->bCellMin[i] == FALSE)
3206 fac += cell_size[i];
3209 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3210 /* Determine the cell boundaries */
3211 for (i = range[0]; i < range[1]; i++)
3213 if (root->bCellMin[i] == FALSE)
3215 cell_size[i] *= fac;
3216 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3218 cellsize_limit_f_i = 0;
3222 cellsize_limit_f_i = cellsize_limit_f;
3224 if (cell_size[i] < cellsize_limit_f_i)
3226 root->bCellMin[i] = TRUE;
3227 cell_size[i] = cellsize_limit_f_i;
3231 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3234 while (nmin > nmin_old);
3237 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3238 /* For this check we should not use DD_CELL_MARGIN,
3239 * but a slightly smaller factor,
3240 * since rounding could get use below the limit.
3242 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3245 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",
3246 gmx_step_str(step, buf),
3247 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3248 ncd, comm->cellsize_min[dim]);
3251 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3255 /* Check if the boundary did not displace more than halfway
3256 * each of the cells it bounds, as this could cause problems,
3257 * especially when the differences between cell sizes are large.
3258 * If changes are applied, they will not make cells smaller
3259 * than the cut-off, as we check all the boundaries which
3260 * might be affected by a change and if the old state was ok,
3261 * the cells will at most be shrunk back to their old size.
3263 for (i = range[0]+1; i < range[1]; i++)
3265 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3266 if (root->cell_f[i] < halfway)
3268 root->cell_f[i] = halfway;
3269 /* Check if the change also causes shifts of the next boundaries */
3270 for (j = i+1; j < range[1]; j++)
3272 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3274 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3278 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3279 if (root->cell_f[i] > halfway)
3281 root->cell_f[i] = halfway;
3282 /* Check if the change also causes shifts of the next boundaries */
3283 for (j = i-1; j >= range[0]+1; j--)
3285 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3287 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3294 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3295 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3296 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3297 * for a and b nrange is used */
3300 /* Take care of the staggering of the cell boundaries */
3303 for (i = range[0]; i < range[1]; i++)
3305 root->cell_f_max0[i] = root->cell_f[i];
3306 root->cell_f_min1[i] = root->cell_f[i+1];
3311 for (i = range[0]+1; i < range[1]; i++)
3313 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3314 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3315 if (bLimLo && bLimHi)
3317 /* Both limits violated, try the best we can */
3318 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3319 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3320 nrange[0] = range[0];
3322 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3325 nrange[1] = range[1];
3326 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3332 /* root->cell_f[i] = root->bound_min[i]; */
3333 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3336 else if (bLimHi && !bLastHi)
3339 if (nrange[1] < range[1]) /* found a LimLo before */
3341 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3342 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3343 nrange[0] = nrange[1];
3345 root->cell_f[i] = root->bound_max[i];
3347 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3349 nrange[1] = range[1];
3352 if (nrange[1] < range[1]) /* found last a LimLo */
3354 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3355 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3356 nrange[0] = nrange[1];
3357 nrange[1] = range[1];
3358 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3360 else if (nrange[0] > range[0]) /* found at least one LimHi */
3362 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3369 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3370 int d, int dim, gmx_domdec_root_t *root,
3371 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3372 gmx_bool bUniform, gmx_large_int_t step)
3374 gmx_domdec_comm_t *comm;
3375 int ncd, d1, i, j, pos;
3377 real load_aver, load_i, imbalance, change, change_max, sc;
3378 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3382 int range[] = { 0, 0 };
3386 /* Convert the maximum change from the input percentage to a fraction */
3387 change_limit = comm->dlb_scale_lim*0.01;
3391 bPBC = (dim < ddbox->npbcdim);
3393 cell_size = root->buf_ncd;
3395 /* Store the original boundaries */
3396 for (i = 0; i < ncd+1; i++)
3398 root->old_cell_f[i] = root->cell_f[i];
3402 for (i = 0; i < ncd; i++)
3404 cell_size[i] = 1.0/ncd;
3407 else if (dd_load_count(comm))
3409 load_aver = comm->load[d].sum_m/ncd;
3411 for (i = 0; i < ncd; i++)
3413 /* Determine the relative imbalance of cell i */
3414 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3415 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3416 /* Determine the change of the cell size using underrelaxation */
3417 change = -relax*imbalance;
3418 change_max = max(change_max, max(change, -change));
3420 /* Limit the amount of scaling.
3421 * We need to use the same rescaling for all cells in one row,
3422 * otherwise the load balancing might not converge.
3425 if (change_max > change_limit)
3427 sc *= change_limit/change_max;
3429 for (i = 0; i < ncd; i++)
3431 /* Determine the relative imbalance of cell i */
3432 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3433 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3434 /* Determine the change of the cell size using underrelaxation */
3435 change = -sc*imbalance;
3436 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3440 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3441 cellsize_limit_f *= DD_CELL_MARGIN;
3442 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3443 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3444 if (ddbox->tric_dir[dim])
3446 cellsize_limit_f /= ddbox->skew_fac[dim];
3447 dist_min_f /= ddbox->skew_fac[dim];
3449 if (bDynamicBox && d > 0)
3451 dist_min_f *= DD_PRES_SCALE_MARGIN;
3453 if (d > 0 && !bUniform)
3455 /* Make sure that the grid is not shifted too much */
3456 for (i = 1; i < ncd; i++)
3458 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3460 gmx_incons("Inconsistent DD boundary staggering limits!");
3462 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3463 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3466 root->bound_min[i] += 0.5*space;
3468 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3469 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3472 root->bound_max[i] += 0.5*space;
3477 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3479 root->cell_f_max0[i-1] + dist_min_f,
3480 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3481 root->cell_f_min1[i] - dist_min_f);
3486 root->cell_f[0] = 0;
3487 root->cell_f[ncd] = 1;
3488 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3491 /* After the checks above, the cells should obey the cut-off
3492 * restrictions, but it does not hurt to check.
3494 for (i = 0; i < ncd; i++)
3498 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3499 dim, i, root->cell_f[i], root->cell_f[i+1]);
3502 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3503 root->cell_f[i+1] - root->cell_f[i] <
3504 cellsize_limit_f/DD_CELL_MARGIN)
3508 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3509 gmx_step_str(step, buf), dim2char(dim), i,
3510 (root->cell_f[i+1] - root->cell_f[i])
3511 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3516 /* Store the cell boundaries of the lower dimensions at the end */
3517 for (d1 = 0; d1 < d; d1++)
3519 root->cell_f[pos++] = comm->cell_f0[d1];
3520 root->cell_f[pos++] = comm->cell_f1[d1];
3523 if (d < comm->npmedecompdim)
3525 /* The master determines the maximum shift for
3526 * the coordinate communication between separate PME nodes.
3528 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3530 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3533 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3537 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3538 gmx_ddbox_t *ddbox, int dimind)
3540 gmx_domdec_comm_t *comm;
3545 /* Set the cell dimensions */
3546 dim = dd->dim[dimind];
3547 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3548 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3549 if (dim >= ddbox->nboundeddim)
3551 comm->cell_x0[dim] += ddbox->box0[dim];
3552 comm->cell_x1[dim] += ddbox->box0[dim];
3556 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3557 int d, int dim, real *cell_f_row,
3560 gmx_domdec_comm_t *comm;
3566 /* Each node would only need to know two fractions,
3567 * but it is probably cheaper to broadcast the whole array.
3569 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3570 0, comm->mpi_comm_load[d]);
3572 /* Copy the fractions for this dimension from the buffer */
3573 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3574 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3575 /* The whole array was communicated, so set the buffer position */
3576 pos = dd->nc[dim] + 1;
3577 for (d1 = 0; d1 <= d; d1++)
3581 /* Copy the cell fractions of the lower dimensions */
3582 comm->cell_f0[d1] = cell_f_row[pos++];
3583 comm->cell_f1[d1] = cell_f_row[pos++];
3585 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3587 /* Convert the communicated shift from float to int */
3588 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3591 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3595 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3596 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3597 gmx_bool bUniform, gmx_large_int_t step)
3599 gmx_domdec_comm_t *comm;
3601 gmx_bool bRowMember, bRowRoot;
3606 for (d = 0; d < dd->ndim; d++)
3611 for (d1 = d; d1 < dd->ndim; d1++)
3613 if (dd->ci[dd->dim[d1]] > 0)
3626 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3627 ddbox, bDynamicBox, bUniform, step);
3628 cell_f_row = comm->root[d]->cell_f;
3632 cell_f_row = comm->cell_f_row;
3634 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3639 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3643 /* This function assumes the box is static and should therefore
3644 * not be called when the box has changed since the last
3645 * call to dd_partition_system.
3647 for (d = 0; d < dd->ndim; d++)
3649 relative_to_absolute_cell_bounds(dd, ddbox, d);
3655 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3656 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3657 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3658 gmx_wallcycle_t wcycle)
3660 gmx_domdec_comm_t *comm;
3667 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3668 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3669 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3671 else if (bDynamicBox)
3673 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3676 /* Set the dimensions for which no DD is used */
3677 for (dim = 0; dim < DIM; dim++)
3679 if (dd->nc[dim] == 1)
3681 comm->cell_x0[dim] = 0;
3682 comm->cell_x1[dim] = ddbox->box_size[dim];
3683 if (dim >= ddbox->nboundeddim)
3685 comm->cell_x0[dim] += ddbox->box0[dim];
3686 comm->cell_x1[dim] += ddbox->box0[dim];
3692 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3695 gmx_domdec_comm_dim_t *cd;
3697 for (d = 0; d < dd->ndim; d++)
3699 cd = &dd->comm->cd[d];
3700 np = npulse[dd->dim[d]];
3701 if (np > cd->np_nalloc)
3705 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3706 dim2char(dd->dim[d]), np);
3708 if (DDMASTER(dd) && cd->np_nalloc > 0)
3710 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3712 srenew(cd->ind, np);
3713 for (i = cd->np_nalloc; i < np; i++)
3715 cd->ind[i].index = NULL;
3716 cd->ind[i].nalloc = 0;
3725 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3726 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3727 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3728 gmx_wallcycle_t wcycle)
3730 gmx_domdec_comm_t *comm;
3736 /* Copy the old cell boundaries for the cg displacement check */
3737 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3738 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3740 if (comm->bDynLoadBal)
3744 check_box_size(dd, ddbox);
3746 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3750 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3751 realloc_comm_ind(dd, npulse);
3756 for (d = 0; d < DIM; d++)
3758 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3759 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3764 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3766 rvec cell_ns_x0, rvec cell_ns_x1,
3767 gmx_large_int_t step)
3769 gmx_domdec_comm_t *comm;
3774 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3776 dim = dd->dim[dim_ind];
3778 /* Without PBC we don't have restrictions on the outer cells */
3779 if (!(dim >= ddbox->npbcdim &&
3780 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3781 comm->bDynLoadBal &&
3782 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3783 comm->cellsize_min[dim])
3786 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",
3787 gmx_step_str(step, buf), dim2char(dim),
3788 comm->cell_x1[dim] - comm->cell_x0[dim],
3789 ddbox->skew_fac[dim],
3790 dd->comm->cellsize_min[dim],
3791 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3795 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3797 /* Communicate the boundaries and update cell_ns_x0/1 */
3798 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3799 if (dd->bGridJump && dd->ndim > 1)
3801 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3806 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3810 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3818 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3819 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3828 static void check_screw_box(matrix box)
3830 /* Mathematical limitation */
3831 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3833 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3836 /* Limitation due to the asymmetry of the eighth shell method */
3837 if (box[ZZ][YY] != 0)
3839 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3843 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3844 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3847 gmx_domdec_master_t *ma;
3848 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3849 int i, icg, j, k, k0, k1, d, npbcdim;
3851 rvec box_size, cg_cm;
3853 real nrcg, inv_ncg, pos_d;
3855 gmx_bool bUnbounded, bScrew;
3859 if (tmp_ind == NULL)
3861 snew(tmp_nalloc, dd->nnodes);
3862 snew(tmp_ind, dd->nnodes);
3863 for (i = 0; i < dd->nnodes; i++)
3865 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3866 snew(tmp_ind[i], tmp_nalloc[i]);
3870 /* Clear the count */
3871 for (i = 0; i < dd->nnodes; i++)
3877 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3879 cgindex = cgs->index;
3881 /* Compute the center of geometry for all charge groups */
3882 for (icg = 0; icg < cgs->nr; icg++)
3885 k1 = cgindex[icg+1];
3889 copy_rvec(pos[k0], cg_cm);
3896 for (k = k0; (k < k1); k++)
3898 rvec_inc(cg_cm, pos[k]);
3900 for (d = 0; (d < DIM); d++)
3902 cg_cm[d] *= inv_ncg;
3905 /* Put the charge group in the box and determine the cell index */
3906 for (d = DIM-1; d >= 0; d--)
3909 if (d < dd->npbcdim)
3911 bScrew = (dd->bScrewPBC && d == XX);
3912 if (tric_dir[d] && dd->nc[d] > 1)
3914 /* Use triclinic coordintates for this dimension */
3915 for (j = d+1; j < DIM; j++)
3917 pos_d += cg_cm[j]*tcm[j][d];
3920 while (pos_d >= box[d][d])
3923 rvec_dec(cg_cm, box[d]);
3926 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3927 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3929 for (k = k0; (k < k1); k++)
3931 rvec_dec(pos[k], box[d]);
3934 pos[k][YY] = box[YY][YY] - pos[k][YY];
3935 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3942 rvec_inc(cg_cm, box[d]);
3945 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3946 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3948 for (k = k0; (k < k1); k++)
3950 rvec_inc(pos[k], box[d]);
3953 pos[k][YY] = box[YY][YY] - pos[k][YY];
3954 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3959 /* This could be done more efficiently */
3961 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3966 i = dd_index(dd->nc, ind);
3967 if (ma->ncg[i] == tmp_nalloc[i])
3969 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3970 srenew(tmp_ind[i], tmp_nalloc[i]);
3972 tmp_ind[i][ma->ncg[i]] = icg;
3974 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3978 for (i = 0; i < dd->nnodes; i++)
3981 for (k = 0; k < ma->ncg[i]; k++)
3983 ma->cg[k1++] = tmp_ind[i][k];
3986 ma->index[dd->nnodes] = k1;
3988 for (i = 0; i < dd->nnodes; i++)
3998 fprintf(fplog, "Charge group distribution at step %s:",
3999 gmx_step_str(step, buf));
4000 for (i = 0; i < dd->nnodes; i++)
4002 fprintf(fplog, " %d", ma->ncg[i]);
4004 fprintf(fplog, "\n");
4008 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4009 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4012 gmx_domdec_master_t *ma = NULL;
4015 int *ibuf, buf2[2] = { 0, 0 };
4016 gmx_bool bMaster = DDMASTER(dd);
4023 check_screw_box(box);
4026 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4028 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4029 for (i = 0; i < dd->nnodes; i++)
4031 ma->ibuf[2*i] = ma->ncg[i];
4032 ma->ibuf[2*i+1] = ma->nat[i];
4040 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4042 dd->ncg_home = buf2[0];
4043 dd->nat_home = buf2[1];
4044 dd->ncg_tot = dd->ncg_home;
4045 dd->nat_tot = dd->nat_home;
4046 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4048 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4049 srenew(dd->index_gl, dd->cg_nalloc);
4050 srenew(dd->cgindex, dd->cg_nalloc+1);
4054 for (i = 0; i < dd->nnodes; i++)
4056 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4057 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4062 DDMASTER(dd) ? ma->ibuf : NULL,
4063 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4064 DDMASTER(dd) ? ma->cg : NULL,
4065 dd->ncg_home*sizeof(int), dd->index_gl);
4067 /* Determine the home charge group sizes */
4069 for (i = 0; i < dd->ncg_home; i++)
4071 cg_gl = dd->index_gl[i];
4073 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4078 fprintf(debug, "Home charge groups:\n");
4079 for (i = 0; i < dd->ncg_home; i++)
4081 fprintf(debug, " %d", dd->index_gl[i]);
4084 fprintf(debug, "\n");
4087 fprintf(debug, "\n");
4091 static int compact_and_copy_vec_at(int ncg, int *move,
4094 rvec *src, gmx_domdec_comm_t *comm,
4097 int m, icg, i, i0, i1, nrcg;
4103 for (m = 0; m < DIM*2; m++)
4109 for (icg = 0; icg < ncg; icg++)
4111 i1 = cgindex[icg+1];
4117 /* Compact the home array in place */
4118 for (i = i0; i < i1; i++)
4120 copy_rvec(src[i], src[home_pos++]);
4126 /* Copy to the communication buffer */
4128 pos_vec[m] += 1 + vec*nrcg;
4129 for (i = i0; i < i1; i++)
4131 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4133 pos_vec[m] += (nvec - vec - 1)*nrcg;
4137 home_pos += i1 - i0;
4145 static int compact_and_copy_vec_cg(int ncg, int *move,
4147 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4150 int m, icg, i0, i1, nrcg;
4156 for (m = 0; m < DIM*2; m++)
4162 for (icg = 0; icg < ncg; icg++)
4164 i1 = cgindex[icg+1];
4170 /* Compact the home array in place */
4171 copy_rvec(src[icg], src[home_pos++]);
4177 /* Copy to the communication buffer */
4178 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4179 pos_vec[m] += 1 + nrcg*nvec;
4191 static int compact_ind(int ncg, int *move,
4192 int *index_gl, int *cgindex,
4194 gmx_ga2la_t ga2la, char *bLocalCG,
4197 int cg, nat, a0, a1, a, a_gl;
4202 for (cg = 0; cg < ncg; cg++)
4208 /* Compact the home arrays in place.
4209 * Anything that can be done here avoids access to global arrays.
4211 cgindex[home_pos] = nat;
4212 for (a = a0; a < a1; a++)
4215 gatindex[nat] = a_gl;
4216 /* The cell number stays 0, so we don't need to set it */
4217 ga2la_change_la(ga2la, a_gl, nat);
4220 index_gl[home_pos] = index_gl[cg];
4221 cginfo[home_pos] = cginfo[cg];
4222 /* The charge group remains local, so bLocalCG does not change */
4227 /* Clear the global indices */
4228 for (a = a0; a < a1; a++)
4230 ga2la_del(ga2la, gatindex[a]);
4234 bLocalCG[index_gl[cg]] = FALSE;
4238 cgindex[home_pos] = nat;
4243 static void clear_and_mark_ind(int ncg, int *move,
4244 int *index_gl, int *cgindex, int *gatindex,
4245 gmx_ga2la_t ga2la, char *bLocalCG,
4250 for (cg = 0; cg < ncg; cg++)
4256 /* Clear the global indices */
4257 for (a = a0; a < a1; a++)
4259 ga2la_del(ga2la, gatindex[a]);
4263 bLocalCG[index_gl[cg]] = FALSE;
4265 /* Signal that this cg has moved using the ns cell index.
4266 * Here we set it to -1. fill_grid will change it
4267 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4269 cell_index[cg] = -1;
4274 static void print_cg_move(FILE *fplog,
4276 gmx_large_int_t step, int cg, int dim, int dir,
4277 gmx_bool bHaveLimitdAndCMOld, real limitd,
4278 rvec cm_old, rvec cm_new, real pos_d)
4280 gmx_domdec_comm_t *comm;
4285 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4286 if (bHaveLimitdAndCMOld)
4288 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4289 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4293 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4294 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4296 fprintf(fplog, "distance out of cell %f\n",
4297 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4298 if (bHaveLimitdAndCMOld)
4300 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4301 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4303 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4304 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4305 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4307 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4308 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4310 comm->cell_x0[dim], comm->cell_x1[dim]);
4313 static void cg_move_error(FILE *fplog,
4315 gmx_large_int_t step, int cg, int dim, int dir,
4316 gmx_bool bHaveLimitdAndCMOld, real limitd,
4317 rvec cm_old, rvec cm_new, real pos_d)
4321 print_cg_move(fplog, dd, step, cg, dim, dir,
4322 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4324 print_cg_move(stderr, dd, step, cg, dim, dir,
4325 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4327 "A charge group moved too far between two domain decomposition steps\n"
4328 "This usually means that your system is not well equilibrated");
4331 static void rotate_state_atom(t_state *state, int a)
4335 for (est = 0; est < estNR; est++)
4337 if (EST_DISTR(est) && (state->flags & (1<<est)))
4342 /* Rotate the complete state; for a rectangular box only */
4343 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4344 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4347 state->v[a][YY] = -state->v[a][YY];
4348 state->v[a][ZZ] = -state->v[a][ZZ];
4351 state->sd_X[a][YY] = -state->sd_X[a][YY];
4352 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4355 state->cg_p[a][YY] = -state->cg_p[a][YY];
4356 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4358 case estDISRE_INITF:
4359 case estDISRE_RM3TAV:
4360 case estORIRE_INITF:
4362 /* These are distances, so not affected by rotation */
4365 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4371 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4373 if (natoms > comm->moved_nalloc)
4375 /* Contents should be preserved here */
4376 comm->moved_nalloc = over_alloc_dd(natoms);
4377 srenew(comm->moved, comm->moved_nalloc);
4383 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4386 ivec tric_dir, matrix tcm,
4387 rvec cell_x0, rvec cell_x1,
4388 rvec limitd, rvec limit0, rvec limit1,
4390 int cg_start, int cg_end,
4395 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4396 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4400 real inv_ncg, pos_d;
4403 npbcdim = dd->npbcdim;
4405 for (cg = cg_start; cg < cg_end; cg++)
4412 copy_rvec(state->x[k0], cm_new);
4419 for (k = k0; (k < k1); k++)
4421 rvec_inc(cm_new, state->x[k]);
4423 for (d = 0; (d < DIM); d++)
4425 cm_new[d] = inv_ncg*cm_new[d];
4430 /* Do pbc and check DD cell boundary crossings */
4431 for (d = DIM-1; d >= 0; d--)
4435 bScrew = (dd->bScrewPBC && d == XX);
4436 /* Determine the location of this cg in lattice coordinates */
4440 for (d2 = d+1; d2 < DIM; d2++)
4442 pos_d += cm_new[d2]*tcm[d2][d];
4445 /* Put the charge group in the triclinic unit-cell */
4446 if (pos_d >= cell_x1[d])
4448 if (pos_d >= limit1[d])
4450 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4451 cg_cm[cg], cm_new, pos_d);
4454 if (dd->ci[d] == dd->nc[d] - 1)
4456 rvec_dec(cm_new, state->box[d]);
4459 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4460 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4462 for (k = k0; (k < k1); k++)
4464 rvec_dec(state->x[k], state->box[d]);
4467 rotate_state_atom(state, k);
4472 else if (pos_d < cell_x0[d])
4474 if (pos_d < limit0[d])
4476 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4477 cg_cm[cg], cm_new, pos_d);
4482 rvec_inc(cm_new, state->box[d]);
4485 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4486 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4488 for (k = k0; (k < k1); k++)
4490 rvec_inc(state->x[k], state->box[d]);
4493 rotate_state_atom(state, k);
4499 else if (d < npbcdim)
4501 /* Put the charge group in the rectangular unit-cell */
4502 while (cm_new[d] >= state->box[d][d])
4504 rvec_dec(cm_new, state->box[d]);
4505 for (k = k0; (k < k1); k++)
4507 rvec_dec(state->x[k], state->box[d]);
4510 while (cm_new[d] < 0)
4512 rvec_inc(cm_new, state->box[d]);
4513 for (k = k0; (k < k1); k++)
4515 rvec_inc(state->x[k], state->box[d]);
4521 copy_rvec(cm_new, cg_cm[cg]);
4523 /* Determine where this cg should go */
4526 for (d = 0; d < dd->ndim; d++)
4531 flag |= DD_FLAG_FW(d);
4537 else if (dev[dim] == -1)
4539 flag |= DD_FLAG_BW(d);
4542 if (dd->nc[dim] > 2)
4553 /* Temporarily store the flag in move */
4554 move[cg] = mc + flag;
4558 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4559 gmx_domdec_t *dd, ivec tric_dir,
4560 t_state *state, rvec **f,
4561 t_forcerec *fr, t_mdatoms *md,
4569 int ncg[DIM*2], nat[DIM*2];
4570 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4571 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4572 int sbuf[2], rbuf[2];
4573 int home_pos_cg, home_pos_at, buf_pos;
4575 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4578 real inv_ncg, pos_d;
4580 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4582 cginfo_mb_t *cginfo_mb;
4583 gmx_domdec_comm_t *comm;
4585 int nthread, thread;
4589 check_screw_box(state->box);
4593 if (fr->cutoff_scheme == ecutsGROUP)
4598 for (i = 0; i < estNR; i++)
4604 case estX: /* Always present */ break;
4605 case estV: bV = (state->flags & (1<<i)); break;
4606 case estSDX: bSDX = (state->flags & (1<<i)); break;
4607 case estCGP: bCGP = (state->flags & (1<<i)); break;
4610 case estDISRE_INITF:
4611 case estDISRE_RM3TAV:
4612 case estORIRE_INITF:
4614 /* No processing required */
4617 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4622 if (dd->ncg_tot > comm->nalloc_int)
4624 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4625 srenew(comm->buf_int, comm->nalloc_int);
4627 move = comm->buf_int;
4629 /* Clear the count */
4630 for (c = 0; c < dd->ndim*2; c++)
4636 npbcdim = dd->npbcdim;
4638 for (d = 0; (d < DIM); d++)
4640 limitd[d] = dd->comm->cellsize_min[d];
4641 if (d >= npbcdim && dd->ci[d] == 0)
4643 cell_x0[d] = -GMX_FLOAT_MAX;
4647 cell_x0[d] = comm->cell_x0[d];
4649 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4651 cell_x1[d] = GMX_FLOAT_MAX;
4655 cell_x1[d] = comm->cell_x1[d];
4659 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4660 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4664 /* We check after communication if a charge group moved
4665 * more than one cell. Set the pre-comm check limit to float_max.
4667 limit0[d] = -GMX_FLOAT_MAX;
4668 limit1[d] = GMX_FLOAT_MAX;
4672 make_tric_corr_matrix(npbcdim, state->box, tcm);
4674 cgindex = dd->cgindex;
4676 nthread = gmx_omp_nthreads_get(emntDomdec);
4678 /* Compute the center of geometry for all home charge groups
4679 * and put them in the box and determine where they should go.
4681 #pragma omp parallel for num_threads(nthread) schedule(static)
4682 for (thread = 0; thread < nthread; thread++)
4684 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4685 cell_x0, cell_x1, limitd, limit0, limit1,
4687 ( thread *dd->ncg_home)/nthread,
4688 ((thread+1)*dd->ncg_home)/nthread,
4689 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4693 for (cg = 0; cg < dd->ncg_home; cg++)
4698 flag = mc & ~DD_FLAG_NRCG;
4699 mc = mc & DD_FLAG_NRCG;
4702 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4704 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4705 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4707 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4708 /* We store the cg size in the lower 16 bits
4709 * and the place where the charge group should go
4710 * in the next 6 bits. This saves some communication volume.
4712 nrcg = cgindex[cg+1] - cgindex[cg];
4713 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4719 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4720 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4723 for (i = 0; i < dd->ndim*2; i++)
4725 *ncg_moved += ncg[i];
4742 /* Make sure the communication buffers are large enough */
4743 for (mc = 0; mc < dd->ndim*2; mc++)
4745 nvr = ncg[mc] + nat[mc]*nvec;
4746 if (nvr > comm->cgcm_state_nalloc[mc])
4748 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4749 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4753 switch (fr->cutoff_scheme)
4756 /* Recalculating cg_cm might be cheaper than communicating,
4757 * but that could give rise to rounding issues.
4760 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4761 nvec, cg_cm, comm, bCompact);
4764 /* Without charge groups we send the moved atom coordinates
4765 * over twice. This is so the code below can be used without
4766 * many conditionals for both for with and without charge groups.
4769 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4770 nvec, state->x, comm, FALSE);
4773 home_pos_cg -= *ncg_moved;
4777 gmx_incons("unimplemented");
4783 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4784 nvec, vec++, state->x, comm, bCompact);
4787 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4788 nvec, vec++, state->v, comm, bCompact);
4792 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4793 nvec, vec++, state->sd_X, comm, bCompact);
4797 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4798 nvec, vec++, state->cg_p, comm, bCompact);
4803 compact_ind(dd->ncg_home, move,
4804 dd->index_gl, dd->cgindex, dd->gatindex,
4805 dd->ga2la, comm->bLocalCG,
4810 if (fr->cutoff_scheme == ecutsVERLET)
4812 moved = get_moved(comm, dd->ncg_home);
4814 for (k = 0; k < dd->ncg_home; k++)
4821 moved = fr->ns.grid->cell_index;
4824 clear_and_mark_ind(dd->ncg_home, move,
4825 dd->index_gl, dd->cgindex, dd->gatindex,
4826 dd->ga2la, comm->bLocalCG,
4830 cginfo_mb = fr->cginfo_mb;
4832 *ncg_stay_home = home_pos_cg;
4833 for (d = 0; d < dd->ndim; d++)
4839 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4842 /* Communicate the cg and atom counts */
4847 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4848 d, dir, sbuf[0], sbuf[1]);
4850 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4852 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4854 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4855 srenew(comm->buf_int, comm->nalloc_int);
4858 /* Communicate the charge group indices, sizes and flags */
4859 dd_sendrecv_int(dd, d, dir,
4860 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4861 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4863 nvs = ncg[cdd] + nat[cdd]*nvec;
4864 i = rbuf[0] + rbuf[1] *nvec;
4865 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4867 /* Communicate cgcm and state */
4868 dd_sendrecv_rvec(dd, d, dir,
4869 comm->cgcm_state[cdd], nvs,
4870 comm->vbuf.v+nvr, i);
4871 ncg_recv += rbuf[0];
4872 nat_recv += rbuf[1];
4876 /* Process the received charge groups */
4878 for (cg = 0; cg < ncg_recv; cg++)
4880 flag = comm->buf_int[cg*DD_CGIBS+1];
4882 if (dim >= npbcdim && dd->nc[dim] > 2)
4884 /* No pbc in this dim and more than one domain boundary.
4885 * We do a separate check if a charge group didn't move too far.
4887 if (((flag & DD_FLAG_FW(d)) &&
4888 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4889 ((flag & DD_FLAG_BW(d)) &&
4890 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4892 cg_move_error(fplog, dd, step, cg, dim,
4893 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4895 comm->vbuf.v[buf_pos],
4896 comm->vbuf.v[buf_pos],
4897 comm->vbuf.v[buf_pos][dim]);
4904 /* Check which direction this cg should go */
4905 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4909 /* The cell boundaries for dimension d2 are not equal
4910 * for each cell row of the lower dimension(s),
4911 * therefore we might need to redetermine where
4912 * this cg should go.
4915 /* If this cg crosses the box boundary in dimension d2
4916 * we can use the communicated flag, so we do not
4917 * have to worry about pbc.
4919 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4920 (flag & DD_FLAG_FW(d2))) ||
4921 (dd->ci[dim2] == 0 &&
4922 (flag & DD_FLAG_BW(d2)))))
4924 /* Clear the two flags for this dimension */
4925 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4926 /* Determine the location of this cg
4927 * in lattice coordinates
4929 pos_d = comm->vbuf.v[buf_pos][dim2];
4932 for (d3 = dim2+1; d3 < DIM; d3++)
4935 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4938 /* Check of we are not at the box edge.
4939 * pbc is only handled in the first step above,
4940 * but this check could move over pbc while
4941 * the first step did not due to different rounding.
4943 if (pos_d >= cell_x1[dim2] &&
4944 dd->ci[dim2] != dd->nc[dim2]-1)
4946 flag |= DD_FLAG_FW(d2);
4948 else if (pos_d < cell_x0[dim2] &&
4951 flag |= DD_FLAG_BW(d2);
4953 comm->buf_int[cg*DD_CGIBS+1] = flag;
4956 /* Set to which neighboring cell this cg should go */
4957 if (flag & DD_FLAG_FW(d2))
4961 else if (flag & DD_FLAG_BW(d2))
4963 if (dd->nc[dd->dim[d2]] > 2)
4975 nrcg = flag & DD_FLAG_NRCG;
4978 if (home_pos_cg+1 > dd->cg_nalloc)
4980 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4981 srenew(dd->index_gl, dd->cg_nalloc);
4982 srenew(dd->cgindex, dd->cg_nalloc+1);
4984 /* Set the global charge group index and size */
4985 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4986 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4987 /* Copy the state from the buffer */
4988 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4989 if (fr->cutoff_scheme == ecutsGROUP)
4992 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4996 /* Set the cginfo */
4997 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4998 dd->index_gl[home_pos_cg]);
5001 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5004 if (home_pos_at+nrcg > state->nalloc)
5006 dd_realloc_state(state, f, home_pos_at+nrcg);
5008 for (i = 0; i < nrcg; i++)
5010 copy_rvec(comm->vbuf.v[buf_pos++],
5011 state->x[home_pos_at+i]);
5015 for (i = 0; i < nrcg; i++)
5017 copy_rvec(comm->vbuf.v[buf_pos++],
5018 state->v[home_pos_at+i]);
5023 for (i = 0; i < nrcg; i++)
5025 copy_rvec(comm->vbuf.v[buf_pos++],
5026 state->sd_X[home_pos_at+i]);
5031 for (i = 0; i < nrcg; i++)
5033 copy_rvec(comm->vbuf.v[buf_pos++],
5034 state->cg_p[home_pos_at+i]);
5038 home_pos_at += nrcg;
5042 /* Reallocate the buffers if necessary */
5043 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5045 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5046 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5048 nvr = ncg[mc] + nat[mc]*nvec;
5049 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5051 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5052 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5054 /* Copy from the receive to the send buffers */
5055 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5056 comm->buf_int + cg*DD_CGIBS,
5057 DD_CGIBS*sizeof(int));
5058 memcpy(comm->cgcm_state[mc][nvr],
5059 comm->vbuf.v[buf_pos],
5060 (1+nrcg*nvec)*sizeof(rvec));
5061 buf_pos += 1 + nrcg*nvec;
5068 /* With sorting (!bCompact) the indices are now only partially up to date
5069 * and ncg_home and nat_home are not the real count, since there are
5070 * "holes" in the arrays for the charge groups that moved to neighbors.
5072 if (fr->cutoff_scheme == ecutsVERLET)
5074 moved = get_moved(comm, home_pos_cg);
5076 for (i = dd->ncg_home; i < home_pos_cg; i++)
5081 dd->ncg_home = home_pos_cg;
5082 dd->nat_home = home_pos_at;
5087 "Finished repartitioning: cgs moved out %d, new home %d\n",
5088 *ncg_moved, dd->ncg_home-*ncg_moved);
5093 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5095 dd->comm->cycl[ddCycl] += cycles;
5096 dd->comm->cycl_n[ddCycl]++;
5097 if (cycles > dd->comm->cycl_max[ddCycl])
5099 dd->comm->cycl_max[ddCycl] = cycles;
5103 static double force_flop_count(t_nrnb *nrnb)
5110 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5112 /* To get closer to the real timings, we half the count
5113 * for the normal loops and again half it for water loops.
5116 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5118 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5122 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5125 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5128 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5130 sum += nrnb->n[i]*cost_nrnb(i);
5133 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5135 sum += nrnb->n[i]*cost_nrnb(i);
5141 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5143 if (dd->comm->eFlop)
5145 dd->comm->flop -= force_flop_count(nrnb);
5148 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5150 if (dd->comm->eFlop)
5152 dd->comm->flop += force_flop_count(nrnb);
5157 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5161 for (i = 0; i < ddCyclNr; i++)
5163 dd->comm->cycl[i] = 0;
5164 dd->comm->cycl_n[i] = 0;
5165 dd->comm->cycl_max[i] = 0;
5168 dd->comm->flop_n = 0;
5171 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5173 gmx_domdec_comm_t *comm;
5174 gmx_domdec_load_t *load;
5175 gmx_domdec_root_t *root = NULL;
5176 int d, dim, cid, i, pos;
5177 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5182 fprintf(debug, "get_load_distribution start\n");
5185 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5189 bSepPME = (dd->pme_nodeid >= 0);
5191 for (d = dd->ndim-1; d >= 0; d--)
5194 /* Check if we participate in the communication in this dimension */
5195 if (d == dd->ndim-1 ||
5196 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5198 load = &comm->load[d];
5201 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5204 if (d == dd->ndim-1)
5206 sbuf[pos++] = dd_force_load(comm);
5207 sbuf[pos++] = sbuf[0];
5210 sbuf[pos++] = sbuf[0];
5211 sbuf[pos++] = cell_frac;
5214 sbuf[pos++] = comm->cell_f_max0[d];
5215 sbuf[pos++] = comm->cell_f_min1[d];
5220 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5221 sbuf[pos++] = comm->cycl[ddCyclPME];
5226 sbuf[pos++] = comm->load[d+1].sum;
5227 sbuf[pos++] = comm->load[d+1].max;
5230 sbuf[pos++] = comm->load[d+1].sum_m;
5231 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5232 sbuf[pos++] = comm->load[d+1].flags;
5235 sbuf[pos++] = comm->cell_f_max0[d];
5236 sbuf[pos++] = comm->cell_f_min1[d];
5241 sbuf[pos++] = comm->load[d+1].mdf;
5242 sbuf[pos++] = comm->load[d+1].pme;
5246 /* Communicate a row in DD direction d.
5247 * The communicators are setup such that the root always has rank 0.
5250 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5251 load->load, load->nload*sizeof(float), MPI_BYTE,
5252 0, comm->mpi_comm_load[d]);
5254 if (dd->ci[dim] == dd->master_ci[dim])
5256 /* We are the root, process this row */
5257 if (comm->bDynLoadBal)
5259 root = comm->root[d];
5269 for (i = 0; i < dd->nc[dim]; i++)
5271 load->sum += load->load[pos++];
5272 load->max = max(load->max, load->load[pos]);
5278 /* This direction could not be load balanced properly,
5279 * therefore we need to use the maximum iso the average load.
5281 load->sum_m = max(load->sum_m, load->load[pos]);
5285 load->sum_m += load->load[pos];
5288 load->cvol_min = min(load->cvol_min, load->load[pos]);
5292 load->flags = (int)(load->load[pos++] + 0.5);
5296 root->cell_f_max0[i] = load->load[pos++];
5297 root->cell_f_min1[i] = load->load[pos++];
5302 load->mdf = max(load->mdf, load->load[pos]);
5304 load->pme = max(load->pme, load->load[pos]);
5308 if (comm->bDynLoadBal && root->bLimited)
5310 load->sum_m *= dd->nc[dim];
5311 load->flags |= (1<<d);
5319 comm->nload += dd_load_count(comm);
5320 comm->load_step += comm->cycl[ddCyclStep];
5321 comm->load_sum += comm->load[0].sum;
5322 comm->load_max += comm->load[0].max;
5323 if (comm->bDynLoadBal)
5325 for (d = 0; d < dd->ndim; d++)
5327 if (comm->load[0].flags & (1<<d))
5329 comm->load_lim[d]++;
5335 comm->load_mdf += comm->load[0].mdf;
5336 comm->load_pme += comm->load[0].pme;
5340 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5344 fprintf(debug, "get_load_distribution finished\n");
5348 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5350 /* Return the relative performance loss on the total run time
5351 * due to the force calculation load imbalance.
5353 if (dd->comm->nload > 0)
5356 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5357 (dd->comm->load_step*dd->nnodes);
5365 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5368 int npp, npme, nnodes, d, limp;
5369 float imbal, pme_f_ratio, lossf, lossp = 0;
5371 gmx_domdec_comm_t *comm;
5374 if (DDMASTER(dd) && comm->nload > 0)
5377 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5378 nnodes = npp + npme;
5379 imbal = comm->load_max*npp/comm->load_sum - 1;
5380 lossf = dd_force_imb_perf_loss(dd);
5381 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5382 fprintf(fplog, "%s", buf);
5383 fprintf(stderr, "\n");
5384 fprintf(stderr, "%s", buf);
5385 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5386 fprintf(fplog, "%s", buf);
5387 fprintf(stderr, "%s", buf);
5389 if (comm->bDynLoadBal)
5391 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5392 for (d = 0; d < dd->ndim; d++)
5394 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5395 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5401 sprintf(buf+strlen(buf), "\n");
5402 fprintf(fplog, "%s", buf);
5403 fprintf(stderr, "%s", buf);
5407 pme_f_ratio = comm->load_pme/comm->load_mdf;
5408 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5411 lossp *= (float)npme/(float)nnodes;
5415 lossp *= (float)npp/(float)nnodes;
5417 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5418 fprintf(fplog, "%s", buf);
5419 fprintf(stderr, "%s", buf);
5420 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5421 fprintf(fplog, "%s", buf);
5422 fprintf(stderr, "%s", buf);
5424 fprintf(fplog, "\n");
5425 fprintf(stderr, "\n");
5427 if (lossf >= DD_PERF_LOSS)
5430 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5431 " in the domain decomposition.\n", lossf*100);
5432 if (!comm->bDynLoadBal)
5434 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5438 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5440 fprintf(fplog, "%s\n", buf);
5441 fprintf(stderr, "%s\n", buf);
5443 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5446 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5447 " had %s work to do than the PP nodes.\n"
5448 " You might want to %s the number of PME nodes\n"
5449 " or %s the cut-off and the grid spacing.\n",
5451 (lossp < 0) ? "less" : "more",
5452 (lossp < 0) ? "decrease" : "increase",
5453 (lossp < 0) ? "decrease" : "increase");
5454 fprintf(fplog, "%s\n", buf);
5455 fprintf(stderr, "%s\n", buf);
5460 static float dd_vol_min(gmx_domdec_t *dd)
5462 return dd->comm->load[0].cvol_min*dd->nnodes;
5465 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5467 return dd->comm->load[0].flags;
5470 static float dd_f_imbal(gmx_domdec_t *dd)
5472 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5475 float dd_pme_f_ratio(gmx_domdec_t *dd)
5477 if (dd->comm->cycl_n[ddCyclPME] > 0)
5479 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5487 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5492 flags = dd_load_flags(dd);
5496 "DD load balancing is limited by minimum cell size in dimension");
5497 for (d = 0; d < dd->ndim; d++)
5501 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5504 fprintf(fplog, "\n");
5506 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5507 if (dd->comm->bDynLoadBal)
5509 fprintf(fplog, " vol min/aver %5.3f%c",
5510 dd_vol_min(dd), flags ? '!' : ' ');
5512 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5513 if (dd->comm->cycl_n[ddCyclPME])
5515 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5517 fprintf(fplog, "\n\n");
5520 static void dd_print_load_verbose(gmx_domdec_t *dd)
5522 if (dd->comm->bDynLoadBal)
5524 fprintf(stderr, "vol %4.2f%c ",
5525 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5527 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5528 if (dd->comm->cycl_n[ddCyclPME])
5530 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5535 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5540 gmx_domdec_root_t *root;
5541 gmx_bool bPartOfGroup = FALSE;
5543 dim = dd->dim[dim_ind];
5544 copy_ivec(loc, loc_c);
5545 for (i = 0; i < dd->nc[dim]; i++)
5548 rank = dd_index(dd->nc, loc_c);
5549 if (rank == dd->rank)
5551 /* This process is part of the group */
5552 bPartOfGroup = TRUE;
5555 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5559 dd->comm->mpi_comm_load[dim_ind] = c_row;
5560 if (dd->comm->eDLB != edlbNO)
5562 if (dd->ci[dim] == dd->master_ci[dim])
5564 /* This is the root process of this row */
5565 snew(dd->comm->root[dim_ind], 1);
5566 root = dd->comm->root[dim_ind];
5567 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5568 snew(root->old_cell_f, dd->nc[dim]+1);
5569 snew(root->bCellMin, dd->nc[dim]);
5572 snew(root->cell_f_max0, dd->nc[dim]);
5573 snew(root->cell_f_min1, dd->nc[dim]);
5574 snew(root->bound_min, dd->nc[dim]);
5575 snew(root->bound_max, dd->nc[dim]);
5577 snew(root->buf_ncd, dd->nc[dim]);
5581 /* This is not a root process, we only need to receive cell_f */
5582 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5585 if (dd->ci[dim] == dd->master_ci[dim])
5587 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5593 static void make_load_communicators(gmx_domdec_t *dd)
5596 int dim0, dim1, i, j;
5601 fprintf(debug, "Making load communicators\n");
5604 snew(dd->comm->load, dd->ndim);
5605 snew(dd->comm->mpi_comm_load, dd->ndim);
5608 make_load_communicator(dd, 0, loc);
5612 for (i = 0; i < dd->nc[dim0]; i++)
5615 make_load_communicator(dd, 1, loc);
5621 for (i = 0; i < dd->nc[dim0]; i++)
5625 for (j = 0; j < dd->nc[dim1]; j++)
5628 make_load_communicator(dd, 2, loc);
5635 fprintf(debug, "Finished making load communicators\n");
5640 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5643 int d, dim, i, j, m;
5646 ivec dd_zp[DD_MAXIZONE];
5647 gmx_domdec_zones_t *zones;
5648 gmx_domdec_ns_ranges_t *izone;
5650 for (d = 0; d < dd->ndim; d++)
5653 copy_ivec(dd->ci, tmp);
5654 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5655 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5656 copy_ivec(dd->ci, tmp);
5657 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5658 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5661 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5664 dd->neighbor[d][1]);
5670 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5672 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5673 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5680 for (i = 0; i < nzonep; i++)
5682 copy_ivec(dd_zp3[i], dd_zp[i]);
5688 for (i = 0; i < nzonep; i++)
5690 copy_ivec(dd_zp2[i], dd_zp[i]);
5696 for (i = 0; i < nzonep; i++)
5698 copy_ivec(dd_zp1[i], dd_zp[i]);
5702 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5707 zones = &dd->comm->zones;
5709 for (i = 0; i < nzone; i++)
5712 clear_ivec(zones->shift[i]);
5713 for (d = 0; d < dd->ndim; d++)
5715 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5720 for (i = 0; i < nzone; i++)
5722 for (d = 0; d < DIM; d++)
5724 s[d] = dd->ci[d] - zones->shift[i][d];
5729 else if (s[d] >= dd->nc[d])
5735 zones->nizone = nzonep;
5736 for (i = 0; i < zones->nizone; i++)
5738 if (dd_zp[i][0] != i)
5740 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5742 izone = &zones->izone[i];
5743 izone->j0 = dd_zp[i][1];
5744 izone->j1 = dd_zp[i][2];
5745 for (dim = 0; dim < DIM; dim++)
5747 if (dd->nc[dim] == 1)
5749 /* All shifts should be allowed */
5750 izone->shift0[dim] = -1;
5751 izone->shift1[dim] = 1;
5756 izone->shift0[d] = 0;
5757 izone->shift1[d] = 0;
5758 for(j=izone->j0; j<izone->j1; j++) {
5759 if (dd->shift[j][d] > dd->shift[i][d])
5760 izone->shift0[d] = -1;
5761 if (dd->shift[j][d] < dd->shift[i][d])
5762 izone->shift1[d] = 1;
5768 /* Assume the shift are not more than 1 cell */
5769 izone->shift0[dim] = 1;
5770 izone->shift1[dim] = -1;
5771 for (j = izone->j0; j < izone->j1; j++)
5773 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5774 if (shift_diff < izone->shift0[dim])
5776 izone->shift0[dim] = shift_diff;
5778 if (shift_diff > izone->shift1[dim])
5780 izone->shift1[dim] = shift_diff;
5787 if (dd->comm->eDLB != edlbNO)
5789 snew(dd->comm->root, dd->ndim);
5792 if (dd->comm->bRecordLoad)
5794 make_load_communicators(dd);
5798 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5801 gmx_domdec_comm_t *comm;
5812 if (comm->bCartesianPP)
5814 /* Set up cartesian communication for the particle-particle part */
5817 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5818 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5821 for (i = 0; i < DIM; i++)
5825 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5827 /* We overwrite the old communicator with the new cartesian one */
5828 cr->mpi_comm_mygroup = comm_cart;
5831 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5832 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5834 if (comm->bCartesianPP_PME)
5836 /* Since we want to use the original cartesian setup for sim,
5837 * and not the one after split, we need to make an index.
5839 snew(comm->ddindex2ddnodeid, dd->nnodes);
5840 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5841 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5842 /* Get the rank of the DD master,
5843 * above we made sure that the master node is a PP node.
5853 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5855 else if (comm->bCartesianPP)
5857 if (cr->npmenodes == 0)
5859 /* The PP communicator is also
5860 * the communicator for this simulation
5862 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5864 cr->nodeid = dd->rank;
5866 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5868 /* We need to make an index to go from the coordinates
5869 * to the nodeid of this simulation.
5871 snew(comm->ddindex2simnodeid, dd->nnodes);
5872 snew(buf, dd->nnodes);
5873 if (cr->duty & DUTY_PP)
5875 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5877 /* Communicate the ddindex to simulation nodeid index */
5878 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5879 cr->mpi_comm_mysim);
5882 /* Determine the master coordinates and rank.
5883 * The DD master should be the same node as the master of this sim.
5885 for (i = 0; i < dd->nnodes; i++)
5887 if (comm->ddindex2simnodeid[i] == 0)
5889 ddindex2xyz(dd->nc, i, dd->master_ci);
5890 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5895 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5900 /* No Cartesian communicators */
5901 /* We use the rank in dd->comm->all as DD index */
5902 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5903 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5905 clear_ivec(dd->master_ci);
5912 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5913 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5918 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5919 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5923 static void receive_ddindex2simnodeid(t_commrec *cr)
5927 gmx_domdec_comm_t *comm;
5934 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5936 snew(comm->ddindex2simnodeid, dd->nnodes);
5937 snew(buf, dd->nnodes);
5938 if (cr->duty & DUTY_PP)
5940 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5943 /* Communicate the ddindex to simulation nodeid index */
5944 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5945 cr->mpi_comm_mysim);
5952 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5953 int ncg, int natoms)
5955 gmx_domdec_master_t *ma;
5960 snew(ma->ncg, dd->nnodes);
5961 snew(ma->index, dd->nnodes+1);
5963 snew(ma->nat, dd->nnodes);
5964 snew(ma->ibuf, dd->nnodes*2);
5965 snew(ma->cell_x, DIM);
5966 for (i = 0; i < DIM; i++)
5968 snew(ma->cell_x[i], dd->nc[i]+1);
5971 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5977 snew(ma->vbuf, natoms);
5983 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5987 gmx_domdec_comm_t *comm;
5998 if (comm->bCartesianPP)
6000 for (i = 1; i < DIM; i++)
6002 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6004 if (bDiv[YY] || bDiv[ZZ])
6006 comm->bCartesianPP_PME = TRUE;
6007 /* If we have 2D PME decomposition, which is always in x+y,
6008 * we stack the PME only nodes in z.
6009 * Otherwise we choose the direction that provides the thinnest slab
6010 * of PME only nodes as this will have the least effect
6011 * on the PP communication.
6012 * But for the PME communication the opposite might be better.
6014 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6016 dd->nc[YY] > dd->nc[ZZ]))
6018 comm->cartpmedim = ZZ;
6022 comm->cartpmedim = YY;
6024 comm->ntot[comm->cartpmedim]
6025 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6029 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]);
6031 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6036 if (comm->bCartesianPP_PME)
6040 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]);
6043 for (i = 0; i < DIM; i++)
6047 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6050 MPI_Comm_rank(comm_cart, &rank);
6051 if (MASTERNODE(cr) && rank != 0)
6053 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6056 /* With this assigment we loose the link to the original communicator
6057 * which will usually be MPI_COMM_WORLD, unless have multisim.
6059 cr->mpi_comm_mysim = comm_cart;
6060 cr->sim_nodeid = rank;
6062 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6066 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6067 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6070 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6074 if (cr->npmenodes == 0 ||
6075 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6077 cr->duty = DUTY_PME;
6080 /* Split the sim communicator into PP and PME only nodes */
6081 MPI_Comm_split(cr->mpi_comm_mysim,
6083 dd_index(comm->ntot, dd->ci),
6084 &cr->mpi_comm_mygroup);
6088 switch (dd_node_order)
6093 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6096 case ddnoINTERLEAVE:
6097 /* Interleave the PP-only and PME-only nodes,
6098 * as on clusters with dual-core machines this will double
6099 * the communication bandwidth of the PME processes
6100 * and thus speed up the PP <-> PME and inter PME communication.
6104 fprintf(fplog, "Interleaving PP and PME nodes\n");
6106 comm->pmenodes = dd_pmenodes(cr);
6111 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6114 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6116 cr->duty = DUTY_PME;
6123 /* Split the sim communicator into PP and PME only nodes */
6124 MPI_Comm_split(cr->mpi_comm_mysim,
6127 &cr->mpi_comm_mygroup);
6128 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6134 fprintf(fplog, "This is a %s only node\n\n",
6135 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6139 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6142 gmx_domdec_comm_t *comm;
6148 copy_ivec(dd->nc, comm->ntot);
6150 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6151 comm->bCartesianPP_PME = FALSE;
6153 /* Reorder the nodes by default. This might change the MPI ranks.
6154 * Real reordering is only supported on very few architectures,
6155 * Blue Gene is one of them.
6157 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6159 if (cr->npmenodes > 0)
6161 /* Split the communicator into a PP and PME part */
6162 split_communicator(fplog, cr, dd_node_order, CartReorder);
6163 if (comm->bCartesianPP_PME)
6165 /* We (possibly) reordered the nodes in split_communicator,
6166 * so it is no longer required in make_pp_communicator.
6168 CartReorder = FALSE;
6173 /* All nodes do PP and PME */
6175 /* We do not require separate communicators */
6176 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6180 if (cr->duty & DUTY_PP)
6182 /* Copy or make a new PP communicator */
6183 make_pp_communicator(fplog, cr, CartReorder);
6187 receive_ddindex2simnodeid(cr);
6190 if (!(cr->duty & DUTY_PME))
6192 /* Set up the commnuication to our PME node */
6193 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6194 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6197 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6198 dd->pme_nodeid, dd->pme_receive_vir_ener);
6203 dd->pme_nodeid = -1;
6208 dd->ma = init_gmx_domdec_master_t(dd,
6210 comm->cgs_gl.index[comm->cgs_gl.nr]);
6214 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6216 real *slb_frac, tot;
6221 if (nc > 1 && size_string != NULL)
6225 fprintf(fplog, "Using static load balancing for the %s direction\n",
6230 for (i = 0; i < nc; i++)
6233 sscanf(size_string, "%lf%n", &dbl, &n);
6236 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6245 fprintf(fplog, "Relative cell sizes:");
6247 for (i = 0; i < nc; i++)
6252 fprintf(fplog, " %5.3f", slb_frac[i]);
6257 fprintf(fplog, "\n");
6264 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6267 gmx_mtop_ilistloop_t iloop;
6271 iloop = gmx_mtop_ilistloop_init(mtop);
6272 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6274 for (ftype = 0; ftype < F_NRE; ftype++)
6276 if ((interaction_function[ftype].flags & IF_BOND) &&
6279 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6287 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6293 val = getenv(env_var);
6296 if (sscanf(val, "%d", &nst) <= 0)
6302 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6310 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6314 fprintf(stderr, "\n%s\n", warn_string);
6318 fprintf(fplog, "\n%s\n", warn_string);
6322 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6323 t_inputrec *ir, FILE *fplog)
6325 if (ir->ePBC == epbcSCREW &&
6326 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6328 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6331 if (ir->ns_type == ensSIMPLE)
6333 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6336 if (ir->nstlist == 0)
6338 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6341 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6343 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6347 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6352 r = ddbox->box_size[XX];
6353 for (di = 0; di < dd->ndim; di++)
6356 /* Check using the initial average cell size */
6357 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6363 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6364 const char *dlb_opt, gmx_bool bRecordLoad,
6365 unsigned long Flags, t_inputrec *ir)
6373 case 'a': eDLB = edlbAUTO; break;
6374 case 'n': eDLB = edlbNO; break;
6375 case 'y': eDLB = edlbYES; break;
6376 default: gmx_incons("Unknown dlb_opt");
6379 if (Flags & MD_RERUN)
6384 if (!EI_DYNAMICS(ir->eI))
6386 if (eDLB == edlbYES)
6388 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6389 dd_warning(cr, fplog, buf);
6397 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6402 if (Flags & MD_REPRODUCIBLE)
6409 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6413 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6416 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6424 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6429 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6431 /* Decomposition order z,y,x */
6434 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6436 for (dim = DIM-1; dim >= 0; dim--)
6438 if (dd->nc[dim] > 1)
6440 dd->dim[dd->ndim++] = dim;
6446 /* Decomposition order x,y,z */
6447 for (dim = 0; dim < DIM; dim++)
6449 if (dd->nc[dim] > 1)
6451 dd->dim[dd->ndim++] = dim;
6457 static gmx_domdec_comm_t *init_dd_comm()
6459 gmx_domdec_comm_t *comm;
6463 snew(comm->cggl_flag, DIM*2);
6464 snew(comm->cgcm_state, DIM*2);
6465 for (i = 0; i < DIM*2; i++)
6467 comm->cggl_flag_nalloc[i] = 0;
6468 comm->cgcm_state_nalloc[i] = 0;
6471 comm->nalloc_int = 0;
6472 comm->buf_int = NULL;
6474 vec_rvec_init(&comm->vbuf);
6476 comm->n_load_have = 0;
6477 comm->n_load_collect = 0;
6479 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6481 comm->sum_nat[i] = 0;
6485 comm->load_step = 0;
6488 clear_ivec(comm->load_lim);
6495 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6496 unsigned long Flags,
6498 real comm_distance_min, real rconstr,
6499 const char *dlb_opt, real dlb_scale,
6500 const char *sizex, const char *sizey, const char *sizez,
6501 gmx_mtop_t *mtop, t_inputrec *ir,
6502 matrix box, rvec *x,
6504 int *npme_x, int *npme_y)
6507 gmx_domdec_comm_t *comm;
6510 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6517 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6522 dd->comm = init_dd_comm();
6524 snew(comm->cggl_flag, DIM*2);
6525 snew(comm->cgcm_state, DIM*2);
6527 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6528 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6530 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6531 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6532 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6533 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6534 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6535 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6536 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6537 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6539 dd->pme_recv_f_alloc = 0;
6540 dd->pme_recv_f_buf = NULL;
6542 if (dd->bSendRecv2 && fplog)
6544 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");
6550 fprintf(fplog, "Will load balance based on FLOP count\n");
6552 if (comm->eFlop > 1)
6554 srand(1+cr->nodeid);
6556 comm->bRecordLoad = TRUE;
6560 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6564 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6566 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6569 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6571 dd->bGridJump = comm->bDynLoadBal;
6572 comm->bPMELoadBalDLBLimits = FALSE;
6574 if (comm->nstSortCG)
6578 if (comm->nstSortCG == 1)
6580 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6584 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6588 snew(comm->sort, 1);
6594 fprintf(fplog, "Will not sort the charge groups\n");
6598 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6600 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6601 if (comm->bInterCGBondeds)
6603 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6607 comm->bInterCGMultiBody = FALSE;
6610 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6611 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6613 if (ir->rlistlong == 0)
6615 /* Set the cut-off to some very large value,
6616 * so we don't need if statements everywhere in the code.
6617 * We use sqrt, since the cut-off is squared in some places.
6619 comm->cutoff = GMX_CUTOFF_INF;
6623 comm->cutoff = ir->rlistlong;
6625 comm->cutoff_mbody = 0;
6627 comm->cellsize_limit = 0;
6628 comm->bBondComm = FALSE;
6630 if (comm->bInterCGBondeds)
6632 if (comm_distance_min > 0)
6634 comm->cutoff_mbody = comm_distance_min;
6635 if (Flags & MD_DDBONDCOMM)
6637 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6641 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6643 r_bonded_limit = comm->cutoff_mbody;
6645 else if (ir->bPeriodicMols)
6647 /* Can not easily determine the required cut-off */
6648 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");
6649 comm->cutoff_mbody = comm->cutoff/2;
6650 r_bonded_limit = comm->cutoff_mbody;
6656 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6657 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6659 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6660 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6662 /* We use an initial margin of 10% for the minimum cell size,
6663 * except when we are just below the non-bonded cut-off.
6665 if (Flags & MD_DDBONDCOMM)
6667 if (max(r_2b, r_mb) > comm->cutoff)
6669 r_bonded = max(r_2b, r_mb);
6670 r_bonded_limit = 1.1*r_bonded;
6671 comm->bBondComm = TRUE;
6676 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6678 /* We determine cutoff_mbody later */
6682 /* No special bonded communication,
6683 * simply increase the DD cut-off.
6685 r_bonded_limit = 1.1*max(r_2b, r_mb);
6686 comm->cutoff_mbody = r_bonded_limit;
6687 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6690 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6694 "Minimum cell size due to bonded interactions: %.3f nm\n",
6695 comm->cellsize_limit);
6699 if (dd->bInterCGcons && rconstr <= 0)
6701 /* There is a cell size limit due to the constraints (P-LINCS) */
6702 rconstr = constr_r_max(fplog, mtop, ir);
6706 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6708 if (rconstr > comm->cellsize_limit)
6710 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6714 else if (rconstr > 0 && fplog)
6716 /* Here we do not check for dd->bInterCGcons,
6717 * because one can also set a cell size limit for virtual sites only
6718 * and at this point we don't know yet if there are intercg v-sites.
6721 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6724 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6726 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6730 copy_ivec(nc, dd->nc);
6731 set_dd_dim(fplog, dd);
6732 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6734 if (cr->npmenodes == -1)
6738 acs = average_cellsize_min(dd, ddbox);
6739 if (acs < comm->cellsize_limit)
6743 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6745 gmx_fatal_collective(FARGS, cr, NULL,
6746 "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",
6747 acs, comm->cellsize_limit);
6752 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6754 /* We need to choose the optimal DD grid and possibly PME nodes */
6755 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6756 comm->eDLB != edlbNO, dlb_scale,
6757 comm->cellsize_limit, comm->cutoff,
6758 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6760 if (dd->nc[XX] == 0)
6762 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6763 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6764 !bC ? "-rdd" : "-rcon",
6765 comm->eDLB != edlbNO ? " or -dds" : "",
6766 bC ? " or your LINCS settings" : "");
6768 gmx_fatal_collective(FARGS, cr, NULL,
6769 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6771 "Look in the log file for details on the domain decomposition",
6772 cr->nnodes-cr->npmenodes, limit, buf);
6774 set_dd_dim(fplog, dd);
6780 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6781 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6784 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6785 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6787 gmx_fatal_collective(FARGS, cr, NULL,
6788 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6789 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6791 if (cr->npmenodes > dd->nnodes)
6793 gmx_fatal_collective(FARGS, cr, NULL,
6794 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6796 if (cr->npmenodes > 0)
6798 comm->npmenodes = cr->npmenodes;
6802 comm->npmenodes = dd->nnodes;
6805 if (EEL_PME(ir->coulombtype))
6807 /* The following choices should match those
6808 * in comm_cost_est in domdec_setup.c.
6809 * Note that here the checks have to take into account
6810 * that the decomposition might occur in a different order than xyz
6811 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6812 * in which case they will not match those in comm_cost_est,
6813 * but since that is mainly for testing purposes that's fine.
6815 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6816 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6817 getenv("GMX_PMEONEDD") == NULL)
6819 comm->npmedecompdim = 2;
6820 comm->npmenodes_x = dd->nc[XX];
6821 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6825 /* In case nc is 1 in both x and y we could still choose to
6826 * decompose pme in y instead of x, but we use x for simplicity.
6828 comm->npmedecompdim = 1;
6829 if (dd->dim[0] == YY)
6831 comm->npmenodes_x = 1;
6832 comm->npmenodes_y = comm->npmenodes;
6836 comm->npmenodes_x = comm->npmenodes;
6837 comm->npmenodes_y = 1;
6842 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6843 comm->npmenodes_x, comm->npmenodes_y, 1);
6848 comm->npmedecompdim = 0;
6849 comm->npmenodes_x = 0;
6850 comm->npmenodes_y = 0;
6853 /* Technically we don't need both of these,
6854 * but it simplifies code not having to recalculate it.
6856 *npme_x = comm->npmenodes_x;
6857 *npme_y = comm->npmenodes_y;
6859 snew(comm->slb_frac, DIM);
6860 if (comm->eDLB == edlbNO)
6862 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6863 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6864 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6867 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6869 if (comm->bBondComm || comm->eDLB != edlbNO)
6871 /* Set the bonded communication distance to halfway
6872 * the minimum and the maximum,
6873 * since the extra communication cost is nearly zero.
6875 acs = average_cellsize_min(dd, ddbox);
6876 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6877 if (comm->eDLB != edlbNO)
6879 /* Check if this does not limit the scaling */
6880 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6882 if (!comm->bBondComm)
6884 /* Without bBondComm do not go beyond the n.b. cut-off */
6885 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6886 if (comm->cellsize_limit >= comm->cutoff)
6888 /* We don't loose a lot of efficieny
6889 * when increasing it to the n.b. cut-off.
6890 * It can even be slightly faster, because we need
6891 * less checks for the communication setup.
6893 comm->cutoff_mbody = comm->cutoff;
6896 /* Check if we did not end up below our original limit */
6897 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6899 if (comm->cutoff_mbody > comm->cellsize_limit)
6901 comm->cellsize_limit = comm->cutoff_mbody;
6904 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6909 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6910 "cellsize limit %f\n",
6911 comm->bBondComm, comm->cellsize_limit);
6916 check_dd_restrictions(cr, dd, ir, fplog);
6919 comm->partition_step = INT_MIN;
6922 clear_dd_cycle_counts(dd);
6927 static void set_dlb_limits(gmx_domdec_t *dd)
6932 for (d = 0; d < dd->ndim; d++)
6934 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6935 dd->comm->cellsize_min[dd->dim[d]] =
6936 dd->comm->cellsize_min_dlb[dd->dim[d]];
6941 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6944 gmx_domdec_comm_t *comm;
6954 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);
6957 cellsize_min = comm->cellsize_min[dd->dim[0]];
6958 for (d = 1; d < dd->ndim; d++)
6960 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6963 if (cellsize_min < comm->cellsize_limit*1.05)
6965 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");
6967 /* Change DLB from "auto" to "no". */
6968 comm->eDLB = edlbNO;
6973 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6974 comm->bDynLoadBal = TRUE;
6975 dd->bGridJump = TRUE;
6979 /* We can set the required cell size info here,
6980 * so we do not need to communicate this.
6981 * The grid is completely uniform.
6983 for (d = 0; d < dd->ndim; d++)
6987 comm->load[d].sum_m = comm->load[d].sum;
6989 nc = dd->nc[dd->dim[d]];
6990 for (i = 0; i < nc; i++)
6992 comm->root[d]->cell_f[i] = i/(real)nc;
6995 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6996 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6999 comm->root[d]->cell_f[nc] = 1.0;
7004 static char *init_bLocalCG(gmx_mtop_t *mtop)
7009 ncg = ncg_mtop(mtop);
7010 snew(bLocalCG, ncg);
7011 for (cg = 0; cg < ncg; cg++)
7013 bLocalCG[cg] = FALSE;
7019 void dd_init_bondeds(FILE *fplog,
7020 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7021 gmx_vsite_t *vsite, gmx_constr_t constr,
7022 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7024 gmx_domdec_comm_t *comm;
7028 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7032 if (comm->bBondComm)
7034 /* Communicate atoms beyond the cut-off for bonded interactions */
7037 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7039 comm->bLocalCG = init_bLocalCG(mtop);
7043 /* Only communicate atoms based on cut-off */
7044 comm->cglink = NULL;
7045 comm->bLocalCG = NULL;
7049 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7051 gmx_bool bDynLoadBal, real dlb_scale,
7054 gmx_domdec_comm_t *comm;
7069 fprintf(fplog, "The maximum number of communication pulses is:");
7070 for (d = 0; d < dd->ndim; d++)
7072 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7074 fprintf(fplog, "\n");
7075 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7076 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7077 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7078 for (d = 0; d < DIM; d++)
7082 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7089 comm->cellsize_min_dlb[d]/
7090 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7092 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7095 fprintf(fplog, "\n");
7099 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7100 fprintf(fplog, "The initial number of communication pulses is:");
7101 for (d = 0; d < dd->ndim; d++)
7103 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7105 fprintf(fplog, "\n");
7106 fprintf(fplog, "The initial domain decomposition cell size is:");
7107 for (d = 0; d < DIM; d++)
7111 fprintf(fplog, " %c %.2f nm",
7112 dim2char(d), dd->comm->cellsize_min[d]);
7115 fprintf(fplog, "\n\n");
7118 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7120 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7121 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7122 "non-bonded interactions", "", comm->cutoff);
7126 limit = dd->comm->cellsize_limit;
7130 if (dynamic_dd_box(ddbox, ir))
7132 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7134 limit = dd->comm->cellsize_min[XX];
7135 for (d = 1; d < DIM; d++)
7137 limit = min(limit, dd->comm->cellsize_min[d]);
7141 if (comm->bInterCGBondeds)
7143 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7144 "two-body bonded interactions", "(-rdd)",
7145 max(comm->cutoff, comm->cutoff_mbody));
7146 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7147 "multi-body bonded interactions", "(-rdd)",
7148 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7152 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7153 "virtual site constructions", "(-rcon)", limit);
7155 if (dd->constraint_comm)
7157 sprintf(buf, "atoms separated by up to %d constraints",
7159 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7160 buf, "(-rcon)", limit);
7162 fprintf(fplog, "\n");
7168 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7170 const t_inputrec *ir,
7171 const gmx_ddbox_t *ddbox)
7173 gmx_domdec_comm_t *comm;
7174 int d, dim, npulse, npulse_d_max, npulse_d;
7179 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7181 /* Determine the maximum number of comm. pulses in one dimension */
7183 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7185 /* Determine the maximum required number of grid pulses */
7186 if (comm->cellsize_limit >= comm->cutoff)
7188 /* Only a single pulse is required */
7191 else if (!bNoCutOff && comm->cellsize_limit > 0)
7193 /* We round down slightly here to avoid overhead due to the latency
7194 * of extra communication calls when the cut-off
7195 * would be only slightly longer than the cell size.
7196 * Later cellsize_limit is redetermined,
7197 * so we can not miss interactions due to this rounding.
7199 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7203 /* There is no cell size limit */
7204 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7207 if (!bNoCutOff && npulse > 1)
7209 /* See if we can do with less pulses, based on dlb_scale */
7211 for (d = 0; d < dd->ndim; d++)
7214 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7215 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7216 npulse_d_max = max(npulse_d_max, npulse_d);
7218 npulse = min(npulse, npulse_d_max);
7221 /* This env var can override npulse */
7222 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7229 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7230 for (d = 0; d < dd->ndim; d++)
7232 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7233 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7234 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7235 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7236 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7238 comm->bVacDLBNoLimit = FALSE;
7242 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7243 if (!comm->bVacDLBNoLimit)
7245 comm->cellsize_limit = max(comm->cellsize_limit,
7246 comm->cutoff/comm->maxpulse);
7248 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7249 /* Set the minimum cell size for each DD dimension */
7250 for (d = 0; d < dd->ndim; d++)
7252 if (comm->bVacDLBNoLimit ||
7253 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7255 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7259 comm->cellsize_min_dlb[dd->dim[d]] =
7260 comm->cutoff/comm->cd[d].np_dlb;
7263 if (comm->cutoff_mbody <= 0)
7265 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7267 if (comm->bDynLoadBal)
7273 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7275 /* If each molecule is a single charge group
7276 * or we use domain decomposition for each periodic dimension,
7277 * we do not need to take pbc into account for the bonded interactions.
7279 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7282 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7285 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7286 t_inputrec *ir, t_forcerec *fr,
7289 gmx_domdec_comm_t *comm;
7295 /* Initialize the thread data.
7296 * This can not be done in init_domain_decomposition,
7297 * as the numbers of threads is determined later.
7299 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7302 snew(comm->dth, comm->nth);
7305 if (EEL_PME(ir->coulombtype))
7307 init_ddpme(dd, &comm->ddpme[0], 0);
7308 if (comm->npmedecompdim >= 2)
7310 init_ddpme(dd, &comm->ddpme[1], 1);
7315 comm->npmenodes = 0;
7316 if (dd->pme_nodeid >= 0)
7318 gmx_fatal_collective(FARGS, NULL, dd,
7319 "Can not have separate PME nodes without PME electrostatics");
7325 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7327 if (comm->eDLB != edlbNO)
7329 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7332 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7333 if (comm->eDLB == edlbAUTO)
7337 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7339 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7342 if (ir->ePBC == epbcNONE)
7344 vol_frac = 1 - 1/(double)dd->nnodes;
7349 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7353 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7355 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7357 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7360 static gmx_bool test_dd_cutoff(t_commrec *cr,
7361 t_state *state, t_inputrec *ir,
7372 set_ddbox(dd, FALSE, cr, ir, state->box,
7373 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7377 for (d = 0; d < dd->ndim; d++)
7381 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7382 if (dynamic_dd_box(&ddbox, ir))
7384 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7387 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7389 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7390 dd->comm->cd[d].np_dlb > 0)
7392 if (np > dd->comm->cd[d].np_dlb)
7397 /* If a current local cell size is smaller than the requested
7398 * cut-off, we could still fix it, but this gets very complicated.
7399 * Without fixing here, we might actually need more checks.
7401 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7408 if (dd->comm->eDLB != edlbNO)
7410 /* If DLB is not active yet, we don't need to check the grid jumps.
7411 * Actually we shouldn't, because then the grid jump data is not set.
7413 if (dd->comm->bDynLoadBal &&
7414 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7419 gmx_sumi(1, &LocallyLimited, cr);
7421 if (LocallyLimited > 0)
7430 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7433 gmx_bool bCutoffAllowed;
7435 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7439 cr->dd->comm->cutoff = cutoff_req;
7442 return bCutoffAllowed;
7445 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7447 gmx_domdec_comm_t *comm;
7449 comm = cr->dd->comm;
7451 /* Turn on the DLB limiting (might have been on already) */
7452 comm->bPMELoadBalDLBLimits = TRUE;
7454 /* Change the cut-off limit */
7455 comm->PMELoadBal_max_cutoff = comm->cutoff;
7458 static void merge_cg_buffers(int ncell,
7459 gmx_domdec_comm_dim_t *cd, int pulse,
7461 int *index_gl, int *recv_i,
7462 rvec *cg_cm, rvec *recv_vr,
7464 cginfo_mb_t *cginfo_mb, int *cginfo)
7466 gmx_domdec_ind_t *ind, *ind_p;
7467 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7468 int shift, shift_at;
7470 ind = &cd->ind[pulse];
7472 /* First correct the already stored data */
7473 shift = ind->nrecv[ncell];
7474 for (cell = ncell-1; cell >= 0; cell--)
7476 shift -= ind->nrecv[cell];
7479 /* Move the cg's present from previous grid pulses */
7480 cg0 = ncg_cell[ncell+cell];
7481 cg1 = ncg_cell[ncell+cell+1];
7482 cgindex[cg1+shift] = cgindex[cg1];
7483 for (cg = cg1-1; cg >= cg0; cg--)
7485 index_gl[cg+shift] = index_gl[cg];
7486 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7487 cgindex[cg+shift] = cgindex[cg];
7488 cginfo[cg+shift] = cginfo[cg];
7490 /* Correct the already stored send indices for the shift */
7491 for (p = 1; p <= pulse; p++)
7493 ind_p = &cd->ind[p];
7495 for (c = 0; c < cell; c++)
7497 cg0 += ind_p->nsend[c];
7499 cg1 = cg0 + ind_p->nsend[cell];
7500 for (cg = cg0; cg < cg1; cg++)
7502 ind_p->index[cg] += shift;
7508 /* Merge in the communicated buffers */
7512 for (cell = 0; cell < ncell; cell++)
7514 cg1 = ncg_cell[ncell+cell+1] + shift;
7517 /* Correct the old cg indices */
7518 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7520 cgindex[cg+1] += shift_at;
7523 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7525 /* Copy this charge group from the buffer */
7526 index_gl[cg1] = recv_i[cg0];
7527 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7528 /* Add it to the cgindex */
7529 cg_gl = index_gl[cg1];
7530 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7531 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7532 cgindex[cg1+1] = cgindex[cg1] + nat;
7537 shift += ind->nrecv[cell];
7538 ncg_cell[ncell+cell+1] = cg1;
7542 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7543 int nzone, int cg0, const int *cgindex)
7547 /* Store the atom block boundaries for easy copying of communication buffers
7550 for (zone = 0; zone < nzone; zone++)
7552 for (p = 0; p < cd->np; p++)
7554 cd->ind[p].cell2at0[zone] = cgindex[cg];
7555 cg += cd->ind[p].nrecv[zone];
7556 cd->ind[p].cell2at1[zone] = cgindex[cg];
7561 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7567 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7569 if (!bLocalCG[link->a[i]])
7578 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7580 real c[DIM][4]; /* the corners for the non-bonded communication */
7581 real cr0; /* corner for rounding */
7582 real cr1[4]; /* corners for rounding */
7583 real bc[DIM]; /* corners for bounded communication */
7584 real bcr1; /* corner for rounding for bonded communication */
7587 /* Determine the corners of the domain(s) we are communicating with */
7589 set_dd_corners(const gmx_domdec_t *dd,
7590 int dim0, int dim1, int dim2,
7594 const gmx_domdec_comm_t *comm;
7595 const gmx_domdec_zones_t *zones;
7600 zones = &comm->zones;
7602 /* Keep the compiler happy */
7606 /* The first dimension is equal for all cells */
7607 c->c[0][0] = comm->cell_x0[dim0];
7610 c->bc[0] = c->c[0][0];
7615 /* This cell row is only seen from the first row */
7616 c->c[1][0] = comm->cell_x0[dim1];
7617 /* All rows can see this row */
7618 c->c[1][1] = comm->cell_x0[dim1];
7621 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7624 /* For the multi-body distance we need the maximum */
7625 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7628 /* Set the upper-right corner for rounding */
7629 c->cr0 = comm->cell_x1[dim0];
7634 for (j = 0; j < 4; j++)
7636 c->c[2][j] = comm->cell_x0[dim2];
7640 /* Use the maximum of the i-cells that see a j-cell */
7641 for (i = 0; i < zones->nizone; i++)
7643 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7649 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7655 /* For the multi-body distance we need the maximum */
7656 c->bc[2] = comm->cell_x0[dim2];
7657 for (i = 0; i < 2; i++)
7659 for (j = 0; j < 2; j++)
7661 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7667 /* Set the upper-right corner for rounding */
7668 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7669 * Only cell (0,0,0) can see cell 7 (1,1,1)
7671 c->cr1[0] = comm->cell_x1[dim1];
7672 c->cr1[3] = comm->cell_x1[dim1];
7675 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7678 /* For the multi-body distance we need the maximum */
7679 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7686 /* Determine which cg's we need to send in this pulse from this zone */
7688 get_zone_pulse_cgs(gmx_domdec_t *dd,
7689 int zonei, int zone,
7691 const int *index_gl,
7693 int dim, int dim_ind,
7694 int dim0, int dim1, int dim2,
7695 real r_comm2, real r_bcomm2,
7699 real skew_fac2_d, real skew_fac_01,
7700 rvec *v_d, rvec *v_0, rvec *v_1,
7701 const dd_corners_t *c,
7703 gmx_bool bDistBonded,
7709 gmx_domdec_ind_t *ind,
7710 int **ibuf, int *ibuf_nalloc,
7716 gmx_domdec_comm_t *comm;
7718 gmx_bool bDistMB_pulse;
7720 real r2, rb2, r, tric_sh;
7723 int nsend_z, nsend, nat;
7727 bScrew = (dd->bScrewPBC && dim == XX);
7729 bDistMB_pulse = (bDistMB && bDistBonded);
7735 for (cg = cg0; cg < cg1; cg++)
7739 if (tric_dist[dim_ind] == 0)
7741 /* Rectangular direction, easy */
7742 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7749 r = cg_cm[cg][dim] - c->bc[dim_ind];
7755 /* Rounding gives at most a 16% reduction
7756 * in communicated atoms
7758 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7760 r = cg_cm[cg][dim0] - c->cr0;
7761 /* This is the first dimension, so always r >= 0 */
7768 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7770 r = cg_cm[cg][dim1] - c->cr1[zone];
7777 r = cg_cm[cg][dim1] - c->bcr1;
7787 /* Triclinic direction, more complicated */
7790 /* Rounding, conservative as the skew_fac multiplication
7791 * will slightly underestimate the distance.
7793 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7795 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7796 for (i = dim0+1; i < DIM; i++)
7798 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7800 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7803 rb[dim0] = rn[dim0];
7806 /* Take care that the cell planes along dim0 might not
7807 * be orthogonal to those along dim1 and dim2.
7809 for (i = 1; i <= dim_ind; i++)
7812 if (normal[dim0][dimd] > 0)
7814 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7817 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7822 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7824 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7826 for (i = dim1+1; i < DIM; i++)
7828 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7830 rn[dim1] += tric_sh;
7833 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7834 /* Take care of coupling of the distances
7835 * to the planes along dim0 and dim1 through dim2.
7837 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7838 /* Take care that the cell planes along dim1
7839 * might not be orthogonal to that along dim2.
7841 if (normal[dim1][dim2] > 0)
7843 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7849 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7852 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7853 /* Take care of coupling of the distances
7854 * to the planes along dim0 and dim1 through dim2.
7856 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7857 /* Take care that the cell planes along dim1
7858 * might not be orthogonal to that along dim2.
7860 if (normal[dim1][dim2] > 0)
7862 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7867 /* The distance along the communication direction */
7868 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7870 for (i = dim+1; i < DIM; i++)
7872 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7877 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7878 /* Take care of coupling of the distances
7879 * to the planes along dim0 and dim1 through dim2.
7881 if (dim_ind == 1 && zonei == 1)
7883 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7889 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7892 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7893 /* Take care of coupling of the distances
7894 * to the planes along dim0 and dim1 through dim2.
7896 if (dim_ind == 1 && zonei == 1)
7898 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7906 ((bDistMB && rb2 < r_bcomm2) ||
7907 (bDist2B && r2 < r_bcomm2)) &&
7909 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7910 missing_link(comm->cglink, index_gl[cg],
7913 /* Make an index to the local charge groups */
7914 if (nsend+1 > ind->nalloc)
7916 ind->nalloc = over_alloc_large(nsend+1);
7917 srenew(ind->index, ind->nalloc);
7919 if (nsend+1 > *ibuf_nalloc)
7921 *ibuf_nalloc = over_alloc_large(nsend+1);
7922 srenew(*ibuf, *ibuf_nalloc);
7924 ind->index[nsend] = cg;
7925 (*ibuf)[nsend] = index_gl[cg];
7927 vec_rvec_check_alloc(vbuf, nsend+1);
7929 if (dd->ci[dim] == 0)
7931 /* Correct cg_cm for pbc */
7932 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7935 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7936 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7941 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7944 nat += cgindex[cg+1] - cgindex[cg];
7950 *nsend_z_ptr = nsend_z;
7953 static void setup_dd_communication(gmx_domdec_t *dd,
7954 matrix box, gmx_ddbox_t *ddbox,
7955 t_forcerec *fr, t_state *state, rvec **f)
7957 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7958 int nzone, nzone_send, zone, zonei, cg0, cg1;
7959 int c, i, j, cg, cg_gl, nrcg;
7960 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7961 gmx_domdec_comm_t *comm;
7962 gmx_domdec_zones_t *zones;
7963 gmx_domdec_comm_dim_t *cd;
7964 gmx_domdec_ind_t *ind;
7965 cginfo_mb_t *cginfo_mb;
7966 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7967 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7968 dd_corners_t corners;
7970 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7971 real skew_fac2_d, skew_fac_01;
7978 fprintf(debug, "Setting up DD communication\n");
7983 switch (fr->cutoff_scheme)
7992 gmx_incons("unimplemented");
7996 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7998 dim = dd->dim[dim_ind];
8000 /* Check if we need to use triclinic distances */
8001 tric_dist[dim_ind] = 0;
8002 for (i = 0; i <= dim_ind; i++)
8004 if (ddbox->tric_dir[dd->dim[i]])
8006 tric_dist[dim_ind] = 1;
8011 bBondComm = comm->bBondComm;
8013 /* Do we need to determine extra distances for multi-body bondeds? */
8014 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8016 /* Do we need to determine extra distances for only two-body bondeds? */
8017 bDist2B = (bBondComm && !bDistMB);
8019 r_comm2 = sqr(comm->cutoff);
8020 r_bcomm2 = sqr(comm->cutoff_mbody);
8024 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8027 zones = &comm->zones;
8030 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8031 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8033 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8035 /* Triclinic stuff */
8036 normal = ddbox->normal;
8040 v_0 = ddbox->v[dim0];
8041 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8043 /* Determine the coupling coefficient for the distances
8044 * to the cell planes along dim0 and dim1 through dim2.
8045 * This is required for correct rounding.
8048 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8051 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8057 v_1 = ddbox->v[dim1];
8060 zone_cg_range = zones->cg_range;
8061 index_gl = dd->index_gl;
8062 cgindex = dd->cgindex;
8063 cginfo_mb = fr->cginfo_mb;
8065 zone_cg_range[0] = 0;
8066 zone_cg_range[1] = dd->ncg_home;
8067 comm->zone_ncg1[0] = dd->ncg_home;
8068 pos_cg = dd->ncg_home;
8070 nat_tot = dd->nat_home;
8072 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8074 dim = dd->dim[dim_ind];
8075 cd = &comm->cd[dim_ind];
8077 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8079 /* No pbc in this dimension, the first node should not comm. */
8087 v_d = ddbox->v[dim];
8088 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8090 cd->bInPlace = TRUE;
8091 for (p = 0; p < cd->np; p++)
8093 /* Only atoms communicated in the first pulse are used
8094 * for multi-body bonded interactions or for bBondComm.
8096 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8101 for (zone = 0; zone < nzone_send; zone++)
8103 if (tric_dist[dim_ind] && dim_ind > 0)
8105 /* Determine slightly more optimized skew_fac's
8107 * This reduces the number of communicated atoms
8108 * by about 10% for 3D DD of rhombic dodecahedra.
8110 for (dimd = 0; dimd < dim; dimd++)
8112 sf2_round[dimd] = 1;
8113 if (ddbox->tric_dir[dimd])
8115 for (i = dd->dim[dimd]+1; i < DIM; i++)
8117 /* If we are shifted in dimension i
8118 * and the cell plane is tilted forward
8119 * in dimension i, skip this coupling.
8121 if (!(zones->shift[nzone+zone][i] &&
8122 ddbox->v[dimd][i][dimd] >= 0))
8125 sqr(ddbox->v[dimd][i][dimd]);
8128 sf2_round[dimd] = 1/sf2_round[dimd];
8133 zonei = zone_perm[dim_ind][zone];
8136 /* Here we permutate the zones to obtain a convenient order
8137 * for neighbor searching
8139 cg0 = zone_cg_range[zonei];
8140 cg1 = zone_cg_range[zonei+1];
8144 /* Look only at the cg's received in the previous grid pulse
8146 cg1 = zone_cg_range[nzone+zone+1];
8147 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8150 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8151 for (th = 0; th < comm->nth; th++)
8153 gmx_domdec_ind_t *ind_p;
8154 int **ibuf_p, *ibuf_nalloc_p;
8156 int *nsend_p, *nat_p;
8162 /* Thread 0 writes in the comm buffers */
8164 ibuf_p = &comm->buf_int;
8165 ibuf_nalloc_p = &comm->nalloc_int;
8166 vbuf_p = &comm->vbuf;
8169 nsend_zone_p = &ind->nsend[zone];
8173 /* Other threads write into temp buffers */
8174 ind_p = &comm->dth[th].ind;
8175 ibuf_p = &comm->dth[th].ibuf;
8176 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8177 vbuf_p = &comm->dth[th].vbuf;
8178 nsend_p = &comm->dth[th].nsend;
8179 nat_p = &comm->dth[th].nat;
8180 nsend_zone_p = &comm->dth[th].nsend_zone;
8182 comm->dth[th].nsend = 0;
8183 comm->dth[th].nat = 0;
8184 comm->dth[th].nsend_zone = 0;
8194 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8195 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8198 /* Get the cg's for this pulse in this zone */
8199 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8201 dim, dim_ind, dim0, dim1, dim2,
8204 normal, skew_fac2_d, skew_fac_01,
8205 v_d, v_0, v_1, &corners, sf2_round,
8206 bDistBonded, bBondComm,
8210 ibuf_p, ibuf_nalloc_p,
8216 /* Append data of threads>=1 to the communication buffers */
8217 for (th = 1; th < comm->nth; th++)
8219 dd_comm_setup_work_t *dth;
8222 dth = &comm->dth[th];
8224 ns1 = nsend + dth->nsend_zone;
8225 if (ns1 > ind->nalloc)
8227 ind->nalloc = over_alloc_dd(ns1);
8228 srenew(ind->index, ind->nalloc);
8230 if (ns1 > comm->nalloc_int)
8232 comm->nalloc_int = over_alloc_dd(ns1);
8233 srenew(comm->buf_int, comm->nalloc_int);
8235 if (ns1 > comm->vbuf.nalloc)
8237 comm->vbuf.nalloc = over_alloc_dd(ns1);
8238 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8241 for (i = 0; i < dth->nsend_zone; i++)
8243 ind->index[nsend] = dth->ind.index[i];
8244 comm->buf_int[nsend] = dth->ibuf[i];
8245 copy_rvec(dth->vbuf.v[i],
8246 comm->vbuf.v[nsend]);
8250 ind->nsend[zone] += dth->nsend_zone;
8253 /* Clear the counts in case we do not have pbc */
8254 for (zone = nzone_send; zone < nzone; zone++)
8256 ind->nsend[zone] = 0;
8258 ind->nsend[nzone] = nsend;
8259 ind->nsend[nzone+1] = nat;
8260 /* Communicate the number of cg's and atoms to receive */
8261 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8262 ind->nsend, nzone+2,
8263 ind->nrecv, nzone+2);
8265 /* The rvec buffer is also required for atom buffers of size nsend
8266 * in dd_move_x and dd_move_f.
8268 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8272 /* We can receive in place if only the last zone is not empty */
8273 for (zone = 0; zone < nzone-1; zone++)
8275 if (ind->nrecv[zone] > 0)
8277 cd->bInPlace = FALSE;
8282 /* The int buffer is only required here for the cg indices */
8283 if (ind->nrecv[nzone] > comm->nalloc_int2)
8285 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8286 srenew(comm->buf_int2, comm->nalloc_int2);
8288 /* The rvec buffer is also required for atom buffers
8289 * of size nrecv in dd_move_x and dd_move_f.
8291 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8292 vec_rvec_check_alloc(&comm->vbuf2, i);
8296 /* Make space for the global cg indices */
8297 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8298 || dd->cg_nalloc == 0)
8300 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8301 srenew(index_gl, dd->cg_nalloc);
8302 srenew(cgindex, dd->cg_nalloc+1);
8304 /* Communicate the global cg indices */
8307 recv_i = index_gl + pos_cg;
8311 recv_i = comm->buf_int2;
8313 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8314 comm->buf_int, nsend,
8315 recv_i, ind->nrecv[nzone]);
8317 /* Make space for cg_cm */
8318 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8319 if (fr->cutoff_scheme == ecutsGROUP)
8327 /* Communicate cg_cm */
8330 recv_vr = cg_cm + pos_cg;
8334 recv_vr = comm->vbuf2.v;
8336 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8337 comm->vbuf.v, nsend,
8338 recv_vr, ind->nrecv[nzone]);
8340 /* Make the charge group index */
8343 zone = (p == 0 ? 0 : nzone - 1);
8344 while (zone < nzone)
8346 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8348 cg_gl = index_gl[pos_cg];
8349 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8350 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8351 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8354 /* Update the charge group presence,
8355 * so we can use it in the next pass of the loop.
8357 comm->bLocalCG[cg_gl] = TRUE;
8363 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8366 zone_cg_range[nzone+zone] = pos_cg;
8371 /* This part of the code is never executed with bBondComm. */
8372 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8373 index_gl, recv_i, cg_cm, recv_vr,
8374 cgindex, fr->cginfo_mb, fr->cginfo);
8375 pos_cg += ind->nrecv[nzone];
8377 nat_tot += ind->nrecv[nzone+1];
8381 /* Store the atom block for easy copying of communication buffers */
8382 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8386 dd->index_gl = index_gl;
8387 dd->cgindex = cgindex;
8389 dd->ncg_tot = zone_cg_range[zones->n];
8390 dd->nat_tot = nat_tot;
8391 comm->nat[ddnatHOME] = dd->nat_home;
8392 for (i = ddnatZONE; i < ddnatNR; i++)
8394 comm->nat[i] = dd->nat_tot;
8399 /* We don't need to update cginfo, since that was alrady done above.
8400 * So we pass NULL for the forcerec.
8402 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8403 NULL, comm->bLocalCG);
8408 fprintf(debug, "Finished setting up DD communication, zones:");
8409 for (c = 0; c < zones->n; c++)
8411 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8413 fprintf(debug, "\n");
8417 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8421 for (c = 0; c < zones->nizone; c++)
8423 zones->izone[c].cg1 = zones->cg_range[c+1];
8424 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8425 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8429 static void set_zones_size(gmx_domdec_t *dd,
8430 matrix box, const gmx_ddbox_t *ddbox,
8431 int zone_start, int zone_end)
8433 gmx_domdec_comm_t *comm;
8434 gmx_domdec_zones_t *zones;
8436 int z, zi, zj0, zj1, d, dim;
8439 real size_j, add_tric;
8444 zones = &comm->zones;
8446 /* Do we need to determine extra distances for multi-body bondeds? */
8447 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8449 for (z = zone_start; z < zone_end; z++)
8451 /* Copy cell limits to zone limits.
8452 * Valid for non-DD dims and non-shifted dims.
8454 copy_rvec(comm->cell_x0, zones->size[z].x0);
8455 copy_rvec(comm->cell_x1, zones->size[z].x1);
8458 for (d = 0; d < dd->ndim; d++)
8462 for (z = 0; z < zones->n; z++)
8464 /* With a staggered grid we have different sizes
8465 * for non-shifted dimensions.
8467 if (dd->bGridJump && zones->shift[z][dim] == 0)
8471 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8472 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8476 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8477 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8483 rcmbs = comm->cutoff_mbody;
8484 if (ddbox->tric_dir[dim])
8486 rcs /= ddbox->skew_fac[dim];
8487 rcmbs /= ddbox->skew_fac[dim];
8490 /* Set the lower limit for the shifted zone dimensions */
8491 for (z = zone_start; z < zone_end; z++)
8493 if (zones->shift[z][dim] > 0)
8496 if (!dd->bGridJump || d == 0)
8498 zones->size[z].x0[dim] = comm->cell_x1[dim];
8499 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8503 /* Here we take the lower limit of the zone from
8504 * the lowest domain of the zone below.
8508 zones->size[z].x0[dim] =
8509 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8515 zones->size[z].x0[dim] =
8516 zones->size[zone_perm[2][z-4]].x0[dim];
8520 zones->size[z].x0[dim] =
8521 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8524 /* A temporary limit, is updated below */
8525 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8529 for (zi = 0; zi < zones->nizone; zi++)
8531 if (zones->shift[zi][dim] == 0)
8533 /* This takes the whole zone into account.
8534 * With multiple pulses this will lead
8535 * to a larger zone then strictly necessary.
8537 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8538 zones->size[zi].x1[dim]+rcmbs);
8546 /* Loop over the i-zones to set the upper limit of each
8549 for (zi = 0; zi < zones->nizone; zi++)
8551 if (zones->shift[zi][dim] == 0)
8553 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8555 if (zones->shift[z][dim] > 0)
8557 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8558 zones->size[zi].x1[dim]+rcs);
8565 for (z = zone_start; z < zone_end; z++)
8567 /* Initialization only required to keep the compiler happy */
8568 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8571 /* To determine the bounding box for a zone we need to find
8572 * the extreme corners of 4, 2 or 1 corners.
8574 nc = 1 << (ddbox->npbcdim - 1);
8576 for (c = 0; c < nc; c++)
8578 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8582 corner[YY] = zones->size[z].x0[YY];
8586 corner[YY] = zones->size[z].x1[YY];
8590 corner[ZZ] = zones->size[z].x0[ZZ];
8594 corner[ZZ] = zones->size[z].x1[ZZ];
8596 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8598 /* With 1D domain decomposition the cg's are not in
8599 * the triclinic box, but triclinic x-y and rectangular y-z.
8600 * Shift y back, so it will later end up at 0.
8602 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8604 /* Apply the triclinic couplings */
8605 for (i = YY; i < ddbox->npbcdim; i++)
8607 for (j = XX; j < i; j++)
8609 corner[j] += corner[i]*box[i][j]/box[i][i];
8614 copy_rvec(corner, corner_min);
8615 copy_rvec(corner, corner_max);
8619 for (i = 0; i < DIM; i++)
8621 corner_min[i] = min(corner_min[i], corner[i]);
8622 corner_max[i] = max(corner_max[i], corner[i]);
8626 /* Copy the extreme cornes without offset along x */
8627 for (i = 0; i < DIM; i++)
8629 zones->size[z].bb_x0[i] = corner_min[i];
8630 zones->size[z].bb_x1[i] = corner_max[i];
8632 /* Add the offset along x */
8633 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8634 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8637 if (zone_start == 0)
8640 for (dim = 0; dim < DIM; dim++)
8642 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8644 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8649 for (z = zone_start; z < zone_end; z++)
8651 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8653 zones->size[z].x0[XX], zones->size[z].x1[XX],
8654 zones->size[z].x0[YY], zones->size[z].x1[YY],
8655 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8656 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8658 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8659 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8660 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8665 static int comp_cgsort(const void *a, const void *b)
8669 gmx_cgsort_t *cga, *cgb;
8670 cga = (gmx_cgsort_t *)a;
8671 cgb = (gmx_cgsort_t *)b;
8673 comp = cga->nsc - cgb->nsc;
8676 comp = cga->ind_gl - cgb->ind_gl;
8682 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8687 /* Order the data */
8688 for (i = 0; i < n; i++)
8690 buf[i] = a[sort[i].ind];
8693 /* Copy back to the original array */
8694 for (i = 0; i < n; i++)
8700 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8705 /* Order the data */
8706 for (i = 0; i < n; i++)
8708 copy_rvec(v[sort[i].ind], buf[i]);
8711 /* Copy back to the original array */
8712 for (i = 0; i < n; i++)
8714 copy_rvec(buf[i], v[i]);
8718 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8721 int a, atot, cg, cg0, cg1, i;
8723 if (cgindex == NULL)
8725 /* Avoid the useless loop of the atoms within a cg */
8726 order_vec_cg(ncg, sort, v, buf);
8731 /* Order the data */
8733 for (cg = 0; cg < ncg; cg++)
8735 cg0 = cgindex[sort[cg].ind];
8736 cg1 = cgindex[sort[cg].ind+1];
8737 for (i = cg0; i < cg1; i++)
8739 copy_rvec(v[i], buf[a]);
8745 /* Copy back to the original array */
8746 for (a = 0; a < atot; a++)
8748 copy_rvec(buf[a], v[a]);
8752 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8753 int nsort_new, gmx_cgsort_t *sort_new,
8754 gmx_cgsort_t *sort1)
8758 /* The new indices are not very ordered, so we qsort them */
8759 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8761 /* sort2 is already ordered, so now we can merge the two arrays */
8765 while (i2 < nsort2 || i_new < nsort_new)
8769 sort1[i1++] = sort_new[i_new++];
8771 else if (i_new == nsort_new)
8773 sort1[i1++] = sort2[i2++];
8775 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8776 (sort2[i2].nsc == sort_new[i_new].nsc &&
8777 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8779 sort1[i1++] = sort2[i2++];
8783 sort1[i1++] = sort_new[i_new++];
8788 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8790 gmx_domdec_sort_t *sort;
8791 gmx_cgsort_t *cgsort, *sort_i;
8792 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8793 int sort_last, sort_skip;
8795 sort = dd->comm->sort;
8797 a = fr->ns.grid->cell_index;
8799 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8801 if (ncg_home_old >= 0)
8803 /* The charge groups that remained in the same ns grid cell
8804 * are completely ordered. So we can sort efficiently by sorting
8805 * the charge groups that did move into the stationary list.
8810 for (i = 0; i < dd->ncg_home; i++)
8812 /* Check if this cg did not move to another node */
8815 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8817 /* This cg is new on this node or moved ns grid cell */
8818 if (nsort_new >= sort->sort_new_nalloc)
8820 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8821 srenew(sort->sort_new, sort->sort_new_nalloc);
8823 sort_i = &(sort->sort_new[nsort_new++]);
8827 /* This cg did not move */
8828 sort_i = &(sort->sort2[nsort2++]);
8830 /* Sort on the ns grid cell indices
8831 * and the global topology index.
8832 * index_gl is irrelevant with cell ns,
8833 * but we set it here anyhow to avoid a conditional.
8836 sort_i->ind_gl = dd->index_gl[i];
8843 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8846 /* Sort efficiently */
8847 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8852 cgsort = sort->sort;
8854 for (i = 0; i < dd->ncg_home; i++)
8856 /* Sort on the ns grid cell indices
8857 * and the global topology index
8859 cgsort[i].nsc = a[i];
8860 cgsort[i].ind_gl = dd->index_gl[i];
8862 if (cgsort[i].nsc < moved)
8869 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8871 /* Determine the order of the charge groups using qsort */
8872 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8878 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8881 int ncg_new, i, *a, na;
8883 sort = dd->comm->sort->sort;
8885 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8888 for (i = 0; i < na; i++)
8892 sort[ncg_new].ind = a[i];
8900 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8901 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, cg0 = 0, 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);
9298 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9299 TRUE, cgs_gl, state_global->x, &ddbox);
9301 get_cg_distribution(fplog, step, dd, cgs_gl,
9302 state_global->box, &ddbox, state_global->x);
9304 dd_distribute_state(dd, cgs_gl,
9305 state_global, state_local, f);
9307 dd_make_local_cgs(dd, &top_local->cgs);
9309 /* Ensure that we have space for the new distribution */
9310 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9312 if (fr->cutoff_scheme == ecutsGROUP)
9314 calc_cgcm(fplog, 0, dd->ncg_home,
9315 &top_local->cgs, state_local->x, fr->cg_cm);
9318 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9320 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9324 else if (state_local->ddp_count != dd->ddp_count)
9326 if (state_local->ddp_count > dd->ddp_count)
9328 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9331 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9333 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);
9336 /* Clear the old state */
9337 clear_dd_indices(dd, 0, 0);
9339 /* Build the new indices */
9340 rebuild_cgindex(dd, cgs_gl->index, state_local);
9341 make_dd_indices(dd, cgs_gl->index, 0);
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);
9366 /* Avoid global communication for dim's without pbc and -gcom */
9367 if (!bNStGlobalComm)
9369 copy_rvec(comm->box0, ddbox.box0 );
9370 copy_rvec(comm->box_size, ddbox.box_size);
9372 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9373 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9378 /* For dim's without pbc and -gcom */
9379 copy_rvec(ddbox.box0, comm->box0 );
9380 copy_rvec(ddbox.box_size, comm->box_size);
9382 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9385 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9387 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9390 /* Check if we should sort the charge groups */
9391 if (comm->nstSortCG > 0)
9393 bSortCG = (bMasterState ||
9394 (bRedist && (step % comm->nstSortCG == 0)));
9401 ncg_home_old = dd->ncg_home;
9406 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9408 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9409 state_local, f, fr, mdatoms,
9410 !bSortCG, nrnb, &cg0, &ncg_moved);
9412 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9415 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9417 &comm->cell_x0, &comm->cell_x1,
9418 dd->ncg_home, fr->cg_cm,
9419 cell_ns_x0, cell_ns_x1, &grid_density);
9423 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9426 switch (fr->cutoff_scheme)
9429 copy_ivec(fr->ns.grid->n, ncells_old);
9430 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9431 state_local->box, cell_ns_x0, cell_ns_x1,
9432 fr->rlistlong, grid_density);
9435 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9438 gmx_incons("unimplemented");
9440 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9441 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9445 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9447 /* Sort the state on charge group position.
9448 * This enables exact restarts from this step.
9449 * It also improves performance by about 15% with larger numbers
9450 * of atoms per node.
9453 /* Fill the ns grid with the home cell,
9454 * so we can sort with the indices.
9456 set_zones_ncg_home(dd);
9458 switch (fr->cutoff_scheme)
9461 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9463 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9465 comm->zones.size[0].bb_x0,
9466 comm->zones.size[0].bb_x1,
9468 comm->zones.dens_zone0,
9471 ncg_moved, bRedist ? comm->moved : NULL,
9472 fr->nbv->grp[eintLocal].kernel_type,
9473 fr->nbv->grp[eintLocal].nbat);
9475 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9478 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9479 0, dd->ncg_home, fr->cg_cm);
9481 copy_ivec(fr->ns.grid->n, ncells_new);
9484 gmx_incons("unimplemented");
9487 bResortAll = bMasterState;
9489 /* Check if we can user the old order and ns grid cell indices
9490 * of the charge groups to sort the charge groups efficiently.
9492 if (ncells_new[XX] != ncells_old[XX] ||
9493 ncells_new[YY] != ncells_old[YY] ||
9494 ncells_new[ZZ] != ncells_old[ZZ])
9501 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9502 gmx_step_str(step, sbuf), dd->ncg_home);
9504 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9505 bResortAll ? -1 : ncg_home_old);
9506 /* Rebuild all the indices */
9508 ga2la_clear(dd->ga2la);
9510 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9513 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9515 /* Setup up the communication and communicate the coordinates */
9516 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9518 /* Set the indices */
9519 make_dd_indices(dd, cgs_gl->index, cg0);
9521 /* Set the charge group boundaries for neighbor searching */
9522 set_cg_boundaries(&comm->zones);
9524 if (fr->cutoff_scheme == ecutsVERLET)
9526 set_zones_size(dd, state_local->box, &ddbox,
9527 bSortCG ? 1 : 0, comm->zones.n);
9530 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9533 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9534 -1,state_local->x,state_local->box);
9537 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9539 /* Extract a local topology from the global topology */
9540 for (i = 0; i < dd->ndim; i++)
9542 np[dd->dim[i]] = comm->cd[i].np;
9544 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9545 comm->cellsize_min, np,
9547 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9548 vsite, top_global, top_local);
9550 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9552 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9554 /* Set up the special atom communication */
9555 n = comm->nat[ddnatZONE];
9556 for (i = ddnatZONE+1; i < ddnatNR; i++)
9561 if (vsite && vsite->n_intercg_vsite)
9563 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9567 if (dd->bInterCGcons || dd->bInterCGsettles)
9569 /* Only for inter-cg constraints we need special code */
9570 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9571 constr, ir->nProjOrder,
9572 top_local->idef.il);
9576 gmx_incons("Unknown special atom type setup");
9581 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9583 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9585 /* Make space for the extra coordinates for virtual site
9586 * or constraint communication.
9588 state_local->natoms = comm->nat[ddnatNR-1];
9589 if (state_local->natoms > state_local->nalloc)
9591 dd_realloc_state(state_local, f, state_local->natoms);
9594 if (fr->bF_NoVirSum)
9596 if (vsite && vsite->n_intercg_vsite)
9598 nat_f_novirsum = comm->nat[ddnatVSITE];
9602 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9604 nat_f_novirsum = dd->nat_tot;
9608 nat_f_novirsum = dd->nat_home;
9617 /* Set the number of atoms required for the force calculation.
9618 * Forces need to be constrained when using a twin-range setup
9619 * or with energy minimization. For simple simulations we could
9620 * avoid some allocation, zeroing and copying, but this is
9621 * probably not worth the complications ande checking.
9623 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9624 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9626 /* We make the all mdatoms up to nat_tot_con.
9627 * We could save some work by only setting invmass
9628 * between nat_tot and nat_tot_con.
9630 /* This call also sets the new number of home particles to dd->nat_home */
9631 atoms2md(top_global, ir,
9632 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9634 /* Now we have the charges we can sort the FE interactions */
9635 dd_sort_local_top(dd, mdatoms, top_local);
9639 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9640 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9645 /* Make the local shell stuff, currently no communication is done */
9646 make_local_shells(cr, mdatoms, shellfc);
9649 if (ir->implicit_solvent)
9651 make_local_gb(cr, fr->born, ir->gb_algorithm);
9654 init_bonded_thread_force_reduction(fr, &top_local->idef);
9656 if (!(cr->duty & DUTY_PME))
9658 /* Send the charges to our PME only node */
9659 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9660 mdatoms->chargeA, mdatoms->chargeB,
9661 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9666 set_constraints(constr, top_local, ir, mdatoms, cr);
9669 if (ir->ePull != epullNO)
9671 /* Update the local pull groups */
9672 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9677 /* Update the local rotation groups */
9678 dd_make_local_rotation_groups(dd, ir->rot);
9682 add_dd_statistics(dd);
9684 /* Make sure we only count the cycles for this DD partitioning */
9685 clear_dd_cycle_counts(dd);
9687 /* Because the order of the atoms might have changed since
9688 * the last vsite construction, we need to communicate the constructing
9689 * atom coordinates again (for spreading the forces this MD step).
9691 dd_move_x_vsites(dd, state_local->box, state_local->x);
9693 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9695 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9697 dd_move_x(dd, state_local->box, state_local->x);
9698 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9699 -1, state_local->x, state_local->box);
9702 /* Store the partitioning step */
9703 comm->partition_step = step;
9705 /* Increase the DD partitioning counter */
9707 /* The state currently matches this DD partitioning count, store it */
9708 state_local->ddp_count = dd->ddp_count;
9711 /* The DD master node knows the complete cg distribution,
9712 * store the count so we can possibly skip the cg info communication.
9714 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9717 if (comm->DD_debug > 0)
9719 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9720 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9721 "after partitioning");