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 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2300 gmx_bool bPMEOnlyNode;
2302 if (DOMAINDECOMP(cr))
2304 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2308 bPMEOnlyNode = FALSE;
2311 return bPMEOnlyNode;
2314 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2315 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2319 ivec coord, coord_pme;
2323 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2326 for (x = 0; x < dd->nc[XX]; x++)
2328 for (y = 0; y < dd->nc[YY]; y++)
2330 for (z = 0; z < dd->nc[ZZ]; z++)
2332 if (dd->comm->bCartesianPP_PME)
2337 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2338 if (dd->ci[XX] == coord_pme[XX] &&
2339 dd->ci[YY] == coord_pme[YY] &&
2340 dd->ci[ZZ] == coord_pme[ZZ])
2342 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2347 /* The slab corresponds to the nodeid in the PME group */
2348 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2350 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2357 /* The last PP-only node is the peer node */
2358 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2362 fprintf(debug, "Receive coordinates from PP nodes:");
2363 for (x = 0; x < *nmy_ddnodes; x++)
2365 fprintf(debug, " %d", (*my_ddnodes)[x]);
2367 fprintf(debug, "\n");
2371 static gmx_bool receive_vir_ener(t_commrec *cr)
2373 gmx_domdec_comm_t *comm;
2374 int pmenode, coords[DIM], rank;
2378 if (cr->npmenodes < cr->dd->nnodes)
2380 comm = cr->dd->comm;
2381 if (comm->bCartesianPP_PME)
2383 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2385 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2386 coords[comm->cartpmedim]++;
2387 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2389 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2390 if (dd_simnode2pmenode(cr, rank) == pmenode)
2392 /* This is not the last PP node for pmenode */
2400 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2401 if (cr->sim_nodeid+1 < cr->nnodes &&
2402 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2404 /* This is not the last PP node for pmenode */
2413 static void set_zones_ncg_home(gmx_domdec_t *dd)
2415 gmx_domdec_zones_t *zones;
2418 zones = &dd->comm->zones;
2420 zones->cg_range[0] = 0;
2421 for (i = 1; i < zones->n+1; i++)
2423 zones->cg_range[i] = dd->ncg_home;
2427 static void rebuild_cgindex(gmx_domdec_t *dd,
2428 const int *gcgs_index, t_state *state)
2430 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2433 dd_cg_gl = dd->index_gl;
2434 cgindex = dd->cgindex;
2437 for (i = 0; i < state->ncg_gl; i++)
2441 dd_cg_gl[i] = cg_gl;
2442 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2446 dd->ncg_home = state->ncg_gl;
2449 set_zones_ncg_home(dd);
2452 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2454 while (cg >= cginfo_mb->cg_end)
2459 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2462 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2463 t_forcerec *fr, char *bLocalCG)
2465 cginfo_mb_t *cginfo_mb;
2471 cginfo_mb = fr->cginfo_mb;
2472 cginfo = fr->cginfo;
2474 for (cg = cg0; cg < cg1; cg++)
2476 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2480 if (bLocalCG != NULL)
2482 for (cg = cg0; cg < cg1; cg++)
2484 bLocalCG[index_gl[cg]] = TRUE;
2489 static void make_dd_indices(gmx_domdec_t *dd,
2490 const int *gcgs_index, int cg_start)
2492 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2493 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2498 bLocalCG = dd->comm->bLocalCG;
2500 if (dd->nat_tot > dd->gatindex_nalloc)
2502 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2503 srenew(dd->gatindex, dd->gatindex_nalloc);
2506 nzone = dd->comm->zones.n;
2507 zone2cg = dd->comm->zones.cg_range;
2508 zone_ncg1 = dd->comm->zone_ncg1;
2509 index_gl = dd->index_gl;
2510 gatindex = dd->gatindex;
2511 bCGs = dd->comm->bCGs;
2513 if (zone2cg[1] != dd->ncg_home)
2515 gmx_incons("dd->ncg_zone is not up to date");
2518 /* Make the local to global and global to local atom index */
2519 a = dd->cgindex[cg_start];
2520 for (zone = 0; zone < nzone; zone++)
2528 cg0 = zone2cg[zone];
2530 cg1 = zone2cg[zone+1];
2531 cg1_p1 = cg0 + zone_ncg1[zone];
2533 for (cg = cg0; cg < cg1; cg++)
2538 /* Signal that this cg is from more than one pulse away */
2541 cg_gl = index_gl[cg];
2544 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2547 ga2la_set(dd->ga2la, a_gl, a, zone1);
2553 gatindex[a] = cg_gl;
2554 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2561 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2564 int ncg, i, ngl, nerr;
2567 if (bLocalCG == NULL)
2571 for (i = 0; i < dd->ncg_tot; i++)
2573 if (!bLocalCG[dd->index_gl[i]])
2576 "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);
2581 for (i = 0; i < ncg_sys; i++)
2588 if (ngl != dd->ncg_tot)
2590 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);
2597 static void check_index_consistency(gmx_domdec_t *dd,
2598 int natoms_sys, int ncg_sys,
2601 int nerr, ngl, i, a, cell;
2606 if (dd->comm->DD_debug > 1)
2608 snew(have, natoms_sys);
2609 for (a = 0; a < dd->nat_tot; a++)
2611 if (have[dd->gatindex[a]] > 0)
2613 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);
2617 have[dd->gatindex[a]] = a + 1;
2623 snew(have, dd->nat_tot);
2626 for (i = 0; i < natoms_sys; i++)
2628 if (ga2la_get(dd->ga2la, i, &a, &cell))
2630 if (a >= dd->nat_tot)
2632 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);
2638 if (dd->gatindex[a] != i)
2640 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);
2647 if (ngl != dd->nat_tot)
2650 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2651 dd->rank, where, ngl, dd->nat_tot);
2653 for (a = 0; a < dd->nat_tot; a++)
2658 "DD node %d, %s: local atom %d, global %d has no global index\n",
2659 dd->rank, where, a+1, dd->gatindex[a]+1);
2664 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2668 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2669 dd->rank, where, nerr);
2673 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2680 /* Clear the whole list without searching */
2681 ga2la_clear(dd->ga2la);
2685 for (i = a_start; i < dd->nat_tot; i++)
2687 ga2la_del(dd->ga2la, dd->gatindex[i]);
2691 bLocalCG = dd->comm->bLocalCG;
2694 for (i = cg_start; i < dd->ncg_tot; i++)
2696 bLocalCG[dd->index_gl[i]] = FALSE;
2700 dd_clear_local_vsite_indices(dd);
2702 if (dd->constraints)
2704 dd_clear_local_constraint_indices(dd);
2708 /* This function should be used for moving the domain boudaries during DLB,
2709 * for obtaining the minimum cell size. It checks the initially set limit
2710 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2711 * and, possibly, a longer cut-off limit set for PME load balancing.
2713 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2717 cellsize_min = comm->cellsize_min[dim];
2719 if (!comm->bVacDLBNoLimit)
2721 /* The cut-off might have changed, e.g. by PME load balacning,
2722 * from the value used to set comm->cellsize_min, so check it.
2724 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2726 if (comm->bPMELoadBalDLBLimits)
2728 /* Check for the cut-off limit set by the PME load balancing */
2729 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2733 return cellsize_min;
2736 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2739 real grid_jump_limit;
2741 /* The distance between the boundaries of cells at distance
2742 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2743 * and by the fact that cells should not be shifted by more than
2744 * half their size, such that cg's only shift by one cell
2745 * at redecomposition.
2747 grid_jump_limit = comm->cellsize_limit;
2748 if (!comm->bVacDLBNoLimit)
2750 if (comm->bPMELoadBalDLBLimits)
2752 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2754 grid_jump_limit = max(grid_jump_limit,
2755 cutoff/comm->cd[dim_ind].np);
2758 return grid_jump_limit;
2761 static gmx_bool check_grid_jump(gmx_large_int_t step,
2767 gmx_domdec_comm_t *comm;
2776 for (d = 1; d < dd->ndim; d++)
2779 limit = grid_jump_limit(comm, cutoff, d);
2780 bfac = ddbox->box_size[dim];
2781 if (ddbox->tric_dir[dim])
2783 bfac *= ddbox->skew_fac[dim];
2785 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2786 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2794 /* This error should never be triggered under normal
2795 * circumstances, but you never know ...
2797 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.",
2798 gmx_step_str(step, buf),
2799 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2807 static int dd_load_count(gmx_domdec_comm_t *comm)
2809 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2812 static float dd_force_load(gmx_domdec_comm_t *comm)
2819 if (comm->eFlop > 1)
2821 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2826 load = comm->cycl[ddCyclF];
2827 if (comm->cycl_n[ddCyclF] > 1)
2829 /* Subtract the maximum of the last n cycle counts
2830 * to get rid of possible high counts due to other soures,
2831 * for instance system activity, that would otherwise
2832 * affect the dynamic load balancing.
2834 load -= comm->cycl_max[ddCyclF];
2841 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2843 gmx_domdec_comm_t *comm;
2848 snew(*dim_f, dd->nc[dim]+1);
2850 for (i = 1; i < dd->nc[dim]; i++)
2852 if (comm->slb_frac[dim])
2854 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2858 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2861 (*dim_f)[dd->nc[dim]] = 1;
2864 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2866 int pmeindex, slab, nso, i;
2869 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2875 ddpme->dim = dimind;
2877 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2879 ddpme->nslab = (ddpme->dim == 0 ?
2880 dd->comm->npmenodes_x :
2881 dd->comm->npmenodes_y);
2883 if (ddpme->nslab <= 1)
2888 nso = dd->comm->npmenodes/ddpme->nslab;
2889 /* Determine for each PME slab the PP location range for dimension dim */
2890 snew(ddpme->pp_min, ddpme->nslab);
2891 snew(ddpme->pp_max, ddpme->nslab);
2892 for (slab = 0; slab < ddpme->nslab; slab++)
2894 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2895 ddpme->pp_max[slab] = 0;
2897 for (i = 0; i < dd->nnodes; i++)
2899 ddindex2xyz(dd->nc, i, xyz);
2900 /* For y only use our y/z slab.
2901 * This assumes that the PME x grid size matches the DD grid size.
2903 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2905 pmeindex = ddindex2pmeindex(dd, i);
2908 slab = pmeindex/nso;
2912 slab = pmeindex % ddpme->nslab;
2914 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2915 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2919 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2922 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2924 if (dd->comm->ddpme[0].dim == XX)
2926 return dd->comm->ddpme[0].maxshift;
2934 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2936 if (dd->comm->ddpme[0].dim == YY)
2938 return dd->comm->ddpme[0].maxshift;
2940 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2942 return dd->comm->ddpme[1].maxshift;
2950 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2951 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2953 gmx_domdec_comm_t *comm;
2956 real range, pme_boundary;
2960 nc = dd->nc[ddpme->dim];
2963 if (!ddpme->dim_match)
2965 /* PP decomposition is not along dim: the worst situation */
2968 else if (ns <= 3 || (bUniform && ns == nc))
2970 /* The optimal situation */
2975 /* We need to check for all pme nodes which nodes they
2976 * could possibly need to communicate with.
2978 xmin = ddpme->pp_min;
2979 xmax = ddpme->pp_max;
2980 /* Allow for atoms to be maximally 2/3 times the cut-off
2981 * out of their DD cell. This is a reasonable balance between
2982 * between performance and support for most charge-group/cut-off
2985 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2986 /* Avoid extra communication when we are exactly at a boundary */
2990 for (s = 0; s < ns; s++)
2992 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2993 pme_boundary = (real)s/ns;
2996 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
2998 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3002 pme_boundary = (real)(s+1)/ns;
3005 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3007 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3014 ddpme->maxshift = sh;
3018 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3019 ddpme->dim, ddpme->maxshift);
3023 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3027 for (d = 0; d < dd->ndim; d++)
3030 if (dim < ddbox->nboundeddim &&
3031 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3032 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3034 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",
3035 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3036 dd->nc[dim], dd->comm->cellsize_limit);
3041 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3042 gmx_bool bMaster, ivec npulse)
3044 gmx_domdec_comm_t *comm;
3047 real *cell_x, cell_dx, cellsize;
3051 for (d = 0; d < DIM; d++)
3053 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3055 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3058 cell_dx = ddbox->box_size[d]/dd->nc[d];
3061 for (j = 0; j < dd->nc[d]+1; j++)
3063 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3068 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3069 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3071 cellsize = cell_dx*ddbox->skew_fac[d];
3072 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3076 cellsize_min[d] = cellsize;
3080 /* Statically load balanced grid */
3081 /* Also when we are not doing a master distribution we determine
3082 * all cell borders in a loop to obtain identical values
3083 * to the master distribution case and to determine npulse.
3087 cell_x = dd->ma->cell_x[d];
3091 snew(cell_x, dd->nc[d]+1);
3093 cell_x[0] = ddbox->box0[d];
3094 for (j = 0; j < dd->nc[d]; j++)
3096 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3097 cell_x[j+1] = cell_x[j] + cell_dx;
3098 cellsize = cell_dx*ddbox->skew_fac[d];
3099 while (cellsize*npulse[d] < comm->cutoff &&
3100 npulse[d] < dd->nc[d]-1)
3104 cellsize_min[d] = min(cellsize_min[d], cellsize);
3108 comm->cell_x0[d] = cell_x[dd->ci[d]];
3109 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3113 /* The following limitation is to avoid that a cell would receive
3114 * some of its own home charge groups back over the periodic boundary.
3115 * Double charge groups cause trouble with the global indices.
3117 if (d < ddbox->npbcdim &&
3118 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3120 gmx_fatal_collective(FARGS, NULL, dd,
3121 "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",
3122 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3124 dd->nc[d], dd->nc[d],
3125 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3129 if (!comm->bDynLoadBal)
3131 copy_rvec(cellsize_min, comm->cellsize_min);
3134 for (d = 0; d < comm->npmedecompdim; d++)
3136 set_pme_maxshift(dd, &comm->ddpme[d],
3137 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3138 comm->ddpme[d].slb_dim_f);
3143 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3144 int d, int dim, gmx_domdec_root_t *root,
3146 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3148 gmx_domdec_comm_t *comm;
3149 int ncd, i, j, nmin, nmin_old;
3150 gmx_bool bLimLo, bLimHi;
3152 real fac, halfway, cellsize_limit_f_i, region_size;
3153 gmx_bool bPBC, bLastHi = FALSE;
3154 int nrange[] = {range[0], range[1]};
3156 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3162 bPBC = (dim < ddbox->npbcdim);
3164 cell_size = root->buf_ncd;
3168 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3171 /* First we need to check if the scaling does not make cells
3172 * smaller than the smallest allowed size.
3173 * We need to do this iteratively, since if a cell is too small,
3174 * it needs to be enlarged, which makes all the other cells smaller,
3175 * which could in turn make another cell smaller than allowed.
3177 for (i = range[0]; i < range[1]; i++)
3179 root->bCellMin[i] = FALSE;
3185 /* We need the total for normalization */
3187 for (i = range[0]; i < range[1]; i++)
3189 if (root->bCellMin[i] == FALSE)
3191 fac += cell_size[i];
3194 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3195 /* Determine the cell boundaries */
3196 for (i = range[0]; i < range[1]; i++)
3198 if (root->bCellMin[i] == FALSE)
3200 cell_size[i] *= fac;
3201 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3203 cellsize_limit_f_i = 0;
3207 cellsize_limit_f_i = cellsize_limit_f;
3209 if (cell_size[i] < cellsize_limit_f_i)
3211 root->bCellMin[i] = TRUE;
3212 cell_size[i] = cellsize_limit_f_i;
3216 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3219 while (nmin > nmin_old);
3222 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3223 /* For this check we should not use DD_CELL_MARGIN,
3224 * but a slightly smaller factor,
3225 * since rounding could get use below the limit.
3227 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3230 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",
3231 gmx_step_str(step, buf),
3232 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3233 ncd, comm->cellsize_min[dim]);
3236 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3240 /* Check if the boundary did not displace more than halfway
3241 * each of the cells it bounds, as this could cause problems,
3242 * especially when the differences between cell sizes are large.
3243 * If changes are applied, they will not make cells smaller
3244 * than the cut-off, as we check all the boundaries which
3245 * might be affected by a change and if the old state was ok,
3246 * the cells will at most be shrunk back to their old size.
3248 for (i = range[0]+1; i < range[1]; i++)
3250 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3251 if (root->cell_f[i] < halfway)
3253 root->cell_f[i] = halfway;
3254 /* Check if the change also causes shifts of the next boundaries */
3255 for (j = i+1; j < range[1]; j++)
3257 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3259 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3263 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3264 if (root->cell_f[i] > halfway)
3266 root->cell_f[i] = halfway;
3267 /* Check if the change also causes shifts of the next boundaries */
3268 for (j = i-1; j >= range[0]+1; j--)
3270 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3272 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3279 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3280 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3281 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3282 * for a and b nrange is used */
3285 /* Take care of the staggering of the cell boundaries */
3288 for (i = range[0]; i < range[1]; i++)
3290 root->cell_f_max0[i] = root->cell_f[i];
3291 root->cell_f_min1[i] = root->cell_f[i+1];
3296 for (i = range[0]+1; i < range[1]; i++)
3298 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3299 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3300 if (bLimLo && bLimHi)
3302 /* Both limits violated, try the best we can */
3303 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3304 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3305 nrange[0] = range[0];
3307 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3310 nrange[1] = range[1];
3311 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3317 /* root->cell_f[i] = root->bound_min[i]; */
3318 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3321 else if (bLimHi && !bLastHi)
3324 if (nrange[1] < range[1]) /* found a LimLo before */
3326 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3327 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3328 nrange[0] = nrange[1];
3330 root->cell_f[i] = root->bound_max[i];
3332 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3334 nrange[1] = range[1];
3337 if (nrange[1] < range[1]) /* found last a LimLo */
3339 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3340 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3341 nrange[0] = nrange[1];
3342 nrange[1] = range[1];
3343 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3345 else if (nrange[0] > range[0]) /* found at least one LimHi */
3347 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3354 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3355 int d, int dim, gmx_domdec_root_t *root,
3356 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3357 gmx_bool bUniform, gmx_large_int_t step)
3359 gmx_domdec_comm_t *comm;
3360 int ncd, d1, i, j, pos;
3362 real load_aver, load_i, imbalance, change, change_max, sc;
3363 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3367 int range[] = { 0, 0 };
3371 /* Convert the maximum change from the input percentage to a fraction */
3372 change_limit = comm->dlb_scale_lim*0.01;
3376 bPBC = (dim < ddbox->npbcdim);
3378 cell_size = root->buf_ncd;
3380 /* Store the original boundaries */
3381 for (i = 0; i < ncd+1; i++)
3383 root->old_cell_f[i] = root->cell_f[i];
3387 for (i = 0; i < ncd; i++)
3389 cell_size[i] = 1.0/ncd;
3392 else if (dd_load_count(comm))
3394 load_aver = comm->load[d].sum_m/ncd;
3396 for (i = 0; i < ncd; i++)
3398 /* Determine the relative imbalance of cell i */
3399 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3400 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3401 /* Determine the change of the cell size using underrelaxation */
3402 change = -relax*imbalance;
3403 change_max = max(change_max, max(change, -change));
3405 /* Limit the amount of scaling.
3406 * We need to use the same rescaling for all cells in one row,
3407 * otherwise the load balancing might not converge.
3410 if (change_max > change_limit)
3412 sc *= change_limit/change_max;
3414 for (i = 0; i < ncd; i++)
3416 /* Determine the relative imbalance of cell i */
3417 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3418 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3419 /* Determine the change of the cell size using underrelaxation */
3420 change = -sc*imbalance;
3421 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3425 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3426 cellsize_limit_f *= DD_CELL_MARGIN;
3427 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3428 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3429 if (ddbox->tric_dir[dim])
3431 cellsize_limit_f /= ddbox->skew_fac[dim];
3432 dist_min_f /= ddbox->skew_fac[dim];
3434 if (bDynamicBox && d > 0)
3436 dist_min_f *= DD_PRES_SCALE_MARGIN;
3438 if (d > 0 && !bUniform)
3440 /* Make sure that the grid is not shifted too much */
3441 for (i = 1; i < ncd; i++)
3443 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3445 gmx_incons("Inconsistent DD boundary staggering limits!");
3447 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3448 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3451 root->bound_min[i] += 0.5*space;
3453 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3454 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3457 root->bound_max[i] += 0.5*space;
3462 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3464 root->cell_f_max0[i-1] + dist_min_f,
3465 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3466 root->cell_f_min1[i] - dist_min_f);
3471 root->cell_f[0] = 0;
3472 root->cell_f[ncd] = 1;
3473 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3476 /* After the checks above, the cells should obey the cut-off
3477 * restrictions, but it does not hurt to check.
3479 for (i = 0; i < ncd; i++)
3483 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3484 dim, i, root->cell_f[i], root->cell_f[i+1]);
3487 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3488 root->cell_f[i+1] - root->cell_f[i] <
3489 cellsize_limit_f/DD_CELL_MARGIN)
3493 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3494 gmx_step_str(step, buf), dim2char(dim), i,
3495 (root->cell_f[i+1] - root->cell_f[i])
3496 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3501 /* Store the cell boundaries of the lower dimensions at the end */
3502 for (d1 = 0; d1 < d; d1++)
3504 root->cell_f[pos++] = comm->cell_f0[d1];
3505 root->cell_f[pos++] = comm->cell_f1[d1];
3508 if (d < comm->npmedecompdim)
3510 /* The master determines the maximum shift for
3511 * the coordinate communication between separate PME nodes.
3513 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3515 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3518 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3522 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3523 gmx_ddbox_t *ddbox, int dimind)
3525 gmx_domdec_comm_t *comm;
3530 /* Set the cell dimensions */
3531 dim = dd->dim[dimind];
3532 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3533 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3534 if (dim >= ddbox->nboundeddim)
3536 comm->cell_x0[dim] += ddbox->box0[dim];
3537 comm->cell_x1[dim] += ddbox->box0[dim];
3541 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3542 int d, int dim, real *cell_f_row,
3545 gmx_domdec_comm_t *comm;
3551 /* Each node would only need to know two fractions,
3552 * but it is probably cheaper to broadcast the whole array.
3554 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3555 0, comm->mpi_comm_load[d]);
3557 /* Copy the fractions for this dimension from the buffer */
3558 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3559 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3560 /* The whole array was communicated, so set the buffer position */
3561 pos = dd->nc[dim] + 1;
3562 for (d1 = 0; d1 <= d; d1++)
3566 /* Copy the cell fractions of the lower dimensions */
3567 comm->cell_f0[d1] = cell_f_row[pos++];
3568 comm->cell_f1[d1] = cell_f_row[pos++];
3570 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3572 /* Convert the communicated shift from float to int */
3573 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3576 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3580 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3581 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3582 gmx_bool bUniform, gmx_large_int_t step)
3584 gmx_domdec_comm_t *comm;
3586 gmx_bool bRowMember, bRowRoot;
3591 for (d = 0; d < dd->ndim; d++)
3596 for (d1 = d; d1 < dd->ndim; d1++)
3598 if (dd->ci[dd->dim[d1]] > 0)
3611 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3612 ddbox, bDynamicBox, bUniform, step);
3613 cell_f_row = comm->root[d]->cell_f;
3617 cell_f_row = comm->cell_f_row;
3619 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3624 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3628 /* This function assumes the box is static and should therefore
3629 * not be called when the box has changed since the last
3630 * call to dd_partition_system.
3632 for (d = 0; d < dd->ndim; d++)
3634 relative_to_absolute_cell_bounds(dd, ddbox, d);
3640 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3641 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3642 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3643 gmx_wallcycle_t wcycle)
3645 gmx_domdec_comm_t *comm;
3652 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3653 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3654 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3656 else if (bDynamicBox)
3658 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3661 /* Set the dimensions for which no DD is used */
3662 for (dim = 0; dim < DIM; dim++)
3664 if (dd->nc[dim] == 1)
3666 comm->cell_x0[dim] = 0;
3667 comm->cell_x1[dim] = ddbox->box_size[dim];
3668 if (dim >= ddbox->nboundeddim)
3670 comm->cell_x0[dim] += ddbox->box0[dim];
3671 comm->cell_x1[dim] += ddbox->box0[dim];
3677 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3680 gmx_domdec_comm_dim_t *cd;
3682 for (d = 0; d < dd->ndim; d++)
3684 cd = &dd->comm->cd[d];
3685 np = npulse[dd->dim[d]];
3686 if (np > cd->np_nalloc)
3690 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3691 dim2char(dd->dim[d]), np);
3693 if (DDMASTER(dd) && cd->np_nalloc > 0)
3695 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3697 srenew(cd->ind, np);
3698 for (i = cd->np_nalloc; i < np; i++)
3700 cd->ind[i].index = NULL;
3701 cd->ind[i].nalloc = 0;
3710 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3711 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3712 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3713 gmx_wallcycle_t wcycle)
3715 gmx_domdec_comm_t *comm;
3721 /* Copy the old cell boundaries for the cg displacement check */
3722 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3723 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3725 if (comm->bDynLoadBal)
3729 check_box_size(dd, ddbox);
3731 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3735 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3736 realloc_comm_ind(dd, npulse);
3741 for (d = 0; d < DIM; d++)
3743 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3744 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3749 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3751 rvec cell_ns_x0, rvec cell_ns_x1,
3752 gmx_large_int_t step)
3754 gmx_domdec_comm_t *comm;
3759 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3761 dim = dd->dim[dim_ind];
3763 /* Without PBC we don't have restrictions on the outer cells */
3764 if (!(dim >= ddbox->npbcdim &&
3765 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3766 comm->bDynLoadBal &&
3767 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3768 comm->cellsize_min[dim])
3771 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",
3772 gmx_step_str(step, buf), dim2char(dim),
3773 comm->cell_x1[dim] - comm->cell_x0[dim],
3774 ddbox->skew_fac[dim],
3775 dd->comm->cellsize_min[dim],
3776 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3780 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3782 /* Communicate the boundaries and update cell_ns_x0/1 */
3783 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3784 if (dd->bGridJump && dd->ndim > 1)
3786 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3791 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3795 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3803 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3804 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3813 static void check_screw_box(matrix box)
3815 /* Mathematical limitation */
3816 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3818 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3821 /* Limitation due to the asymmetry of the eighth shell method */
3822 if (box[ZZ][YY] != 0)
3824 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3828 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3829 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3832 gmx_domdec_master_t *ma;
3833 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3834 int i, icg, j, k, k0, k1, d, npbcdim;
3836 rvec box_size, cg_cm;
3838 real nrcg, inv_ncg, pos_d;
3840 gmx_bool bUnbounded, bScrew;
3844 if (tmp_ind == NULL)
3846 snew(tmp_nalloc, dd->nnodes);
3847 snew(tmp_ind, dd->nnodes);
3848 for (i = 0; i < dd->nnodes; i++)
3850 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3851 snew(tmp_ind[i], tmp_nalloc[i]);
3855 /* Clear the count */
3856 for (i = 0; i < dd->nnodes; i++)
3862 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3864 cgindex = cgs->index;
3866 /* Compute the center of geometry for all charge groups */
3867 for (icg = 0; icg < cgs->nr; icg++)
3870 k1 = cgindex[icg+1];
3874 copy_rvec(pos[k0], cg_cm);
3881 for (k = k0; (k < k1); k++)
3883 rvec_inc(cg_cm, pos[k]);
3885 for (d = 0; (d < DIM); d++)
3887 cg_cm[d] *= inv_ncg;
3890 /* Put the charge group in the box and determine the cell index */
3891 for (d = DIM-1; d >= 0; d--)
3894 if (d < dd->npbcdim)
3896 bScrew = (dd->bScrewPBC && d == XX);
3897 if (tric_dir[d] && dd->nc[d] > 1)
3899 /* Use triclinic coordintates for this dimension */
3900 for (j = d+1; j < DIM; j++)
3902 pos_d += cg_cm[j]*tcm[j][d];
3905 while (pos_d >= box[d][d])
3908 rvec_dec(cg_cm, box[d]);
3911 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3912 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3914 for (k = k0; (k < k1); k++)
3916 rvec_dec(pos[k], box[d]);
3919 pos[k][YY] = box[YY][YY] - pos[k][YY];
3920 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3927 rvec_inc(cg_cm, box[d]);
3930 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3931 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3933 for (k = k0; (k < k1); k++)
3935 rvec_inc(pos[k], box[d]);
3938 pos[k][YY] = box[YY][YY] - pos[k][YY];
3939 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3944 /* This could be done more efficiently */
3946 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3951 i = dd_index(dd->nc, ind);
3952 if (ma->ncg[i] == tmp_nalloc[i])
3954 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3955 srenew(tmp_ind[i], tmp_nalloc[i]);
3957 tmp_ind[i][ma->ncg[i]] = icg;
3959 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3963 for (i = 0; i < dd->nnodes; i++)
3966 for (k = 0; k < ma->ncg[i]; k++)
3968 ma->cg[k1++] = tmp_ind[i][k];
3971 ma->index[dd->nnodes] = k1;
3973 for (i = 0; i < dd->nnodes; i++)
3983 fprintf(fplog, "Charge group distribution at step %s:",
3984 gmx_step_str(step, buf));
3985 for (i = 0; i < dd->nnodes; i++)
3987 fprintf(fplog, " %d", ma->ncg[i]);
3989 fprintf(fplog, "\n");
3993 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
3994 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
3997 gmx_domdec_master_t *ma = NULL;
4000 int *ibuf, buf2[2] = { 0, 0 };
4001 gmx_bool bMaster = DDMASTER(dd);
4008 check_screw_box(box);
4011 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4013 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4014 for (i = 0; i < dd->nnodes; i++)
4016 ma->ibuf[2*i] = ma->ncg[i];
4017 ma->ibuf[2*i+1] = ma->nat[i];
4025 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4027 dd->ncg_home = buf2[0];
4028 dd->nat_home = buf2[1];
4029 dd->ncg_tot = dd->ncg_home;
4030 dd->nat_tot = dd->nat_home;
4031 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4033 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4034 srenew(dd->index_gl, dd->cg_nalloc);
4035 srenew(dd->cgindex, dd->cg_nalloc+1);
4039 for (i = 0; i < dd->nnodes; i++)
4041 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4042 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4047 DDMASTER(dd) ? ma->ibuf : NULL,
4048 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4049 DDMASTER(dd) ? ma->cg : NULL,
4050 dd->ncg_home*sizeof(int), dd->index_gl);
4052 /* Determine the home charge group sizes */
4054 for (i = 0; i < dd->ncg_home; i++)
4056 cg_gl = dd->index_gl[i];
4058 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4063 fprintf(debug, "Home charge groups:\n");
4064 for (i = 0; i < dd->ncg_home; i++)
4066 fprintf(debug, " %d", dd->index_gl[i]);
4069 fprintf(debug, "\n");
4072 fprintf(debug, "\n");
4076 static int compact_and_copy_vec_at(int ncg, int *move,
4079 rvec *src, gmx_domdec_comm_t *comm,
4082 int m, icg, i, i0, i1, nrcg;
4088 for (m = 0; m < DIM*2; m++)
4094 for (icg = 0; icg < ncg; icg++)
4096 i1 = cgindex[icg+1];
4102 /* Compact the home array in place */
4103 for (i = i0; i < i1; i++)
4105 copy_rvec(src[i], src[home_pos++]);
4111 /* Copy to the communication buffer */
4113 pos_vec[m] += 1 + vec*nrcg;
4114 for (i = i0; i < i1; i++)
4116 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4118 pos_vec[m] += (nvec - vec - 1)*nrcg;
4122 home_pos += i1 - i0;
4130 static int compact_and_copy_vec_cg(int ncg, int *move,
4132 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4135 int m, icg, i0, i1, nrcg;
4141 for (m = 0; m < DIM*2; m++)
4147 for (icg = 0; icg < ncg; icg++)
4149 i1 = cgindex[icg+1];
4155 /* Compact the home array in place */
4156 copy_rvec(src[icg], src[home_pos++]);
4162 /* Copy to the communication buffer */
4163 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4164 pos_vec[m] += 1 + nrcg*nvec;
4176 static int compact_ind(int ncg, int *move,
4177 int *index_gl, int *cgindex,
4179 gmx_ga2la_t ga2la, char *bLocalCG,
4182 int cg, nat, a0, a1, a, a_gl;
4187 for (cg = 0; cg < ncg; cg++)
4193 /* Compact the home arrays in place.
4194 * Anything that can be done here avoids access to global arrays.
4196 cgindex[home_pos] = nat;
4197 for (a = a0; a < a1; a++)
4200 gatindex[nat] = a_gl;
4201 /* The cell number stays 0, so we don't need to set it */
4202 ga2la_change_la(ga2la, a_gl, nat);
4205 index_gl[home_pos] = index_gl[cg];
4206 cginfo[home_pos] = cginfo[cg];
4207 /* The charge group remains local, so bLocalCG does not change */
4212 /* Clear the global indices */
4213 for (a = a0; a < a1; a++)
4215 ga2la_del(ga2la, gatindex[a]);
4219 bLocalCG[index_gl[cg]] = FALSE;
4223 cgindex[home_pos] = nat;
4228 static void clear_and_mark_ind(int ncg, int *move,
4229 int *index_gl, int *cgindex, int *gatindex,
4230 gmx_ga2la_t ga2la, char *bLocalCG,
4235 for (cg = 0; cg < ncg; cg++)
4241 /* Clear the global indices */
4242 for (a = a0; a < a1; a++)
4244 ga2la_del(ga2la, gatindex[a]);
4248 bLocalCG[index_gl[cg]] = FALSE;
4250 /* Signal that this cg has moved using the ns cell index.
4251 * Here we set it to -1. fill_grid will change it
4252 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4254 cell_index[cg] = -1;
4259 static void print_cg_move(FILE *fplog,
4261 gmx_large_int_t step, int cg, int dim, int dir,
4262 gmx_bool bHaveLimitdAndCMOld, real limitd,
4263 rvec cm_old, rvec cm_new, real pos_d)
4265 gmx_domdec_comm_t *comm;
4270 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4271 if (bHaveLimitdAndCMOld)
4273 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4274 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4278 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4279 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4281 fprintf(fplog, "distance out of cell %f\n",
4282 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4283 if (bHaveLimitdAndCMOld)
4285 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4286 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4288 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4289 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4290 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4292 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4293 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4295 comm->cell_x0[dim], comm->cell_x1[dim]);
4298 static void cg_move_error(FILE *fplog,
4300 gmx_large_int_t step, int cg, int dim, int dir,
4301 gmx_bool bHaveLimitdAndCMOld, real limitd,
4302 rvec cm_old, rvec cm_new, real pos_d)
4306 print_cg_move(fplog, dd, step, cg, dim, dir,
4307 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4309 print_cg_move(stderr, dd, step, cg, dim, dir,
4310 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4312 "A charge group moved too far between two domain decomposition steps\n"
4313 "This usually means that your system is not well equilibrated");
4316 static void rotate_state_atom(t_state *state, int a)
4320 for (est = 0; est < estNR; est++)
4322 if (EST_DISTR(est) && (state->flags & (1<<est)))
4327 /* Rotate the complete state; for a rectangular box only */
4328 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4329 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4332 state->v[a][YY] = -state->v[a][YY];
4333 state->v[a][ZZ] = -state->v[a][ZZ];
4336 state->sd_X[a][YY] = -state->sd_X[a][YY];
4337 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4340 state->cg_p[a][YY] = -state->cg_p[a][YY];
4341 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4343 case estDISRE_INITF:
4344 case estDISRE_RM3TAV:
4345 case estORIRE_INITF:
4347 /* These are distances, so not affected by rotation */
4350 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4356 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4358 if (natoms > comm->moved_nalloc)
4360 /* Contents should be preserved here */
4361 comm->moved_nalloc = over_alloc_dd(natoms);
4362 srenew(comm->moved, comm->moved_nalloc);
4368 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4371 ivec tric_dir, matrix tcm,
4372 rvec cell_x0, rvec cell_x1,
4373 rvec limitd, rvec limit0, rvec limit1,
4375 int cg_start, int cg_end,
4380 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4381 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4385 real inv_ncg, pos_d;
4388 npbcdim = dd->npbcdim;
4390 for (cg = cg_start; cg < cg_end; cg++)
4397 copy_rvec(state->x[k0], cm_new);
4404 for (k = k0; (k < k1); k++)
4406 rvec_inc(cm_new, state->x[k]);
4408 for (d = 0; (d < DIM); d++)
4410 cm_new[d] = inv_ncg*cm_new[d];
4415 /* Do pbc and check DD cell boundary crossings */
4416 for (d = DIM-1; d >= 0; d--)
4420 bScrew = (dd->bScrewPBC && d == XX);
4421 /* Determine the location of this cg in lattice coordinates */
4425 for (d2 = d+1; d2 < DIM; d2++)
4427 pos_d += cm_new[d2]*tcm[d2][d];
4430 /* Put the charge group in the triclinic unit-cell */
4431 if (pos_d >= cell_x1[d])
4433 if (pos_d >= limit1[d])
4435 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4436 cg_cm[cg], cm_new, pos_d);
4439 if (dd->ci[d] == dd->nc[d] - 1)
4441 rvec_dec(cm_new, state->box[d]);
4444 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4445 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4447 for (k = k0; (k < k1); k++)
4449 rvec_dec(state->x[k], state->box[d]);
4452 rotate_state_atom(state, k);
4457 else if (pos_d < cell_x0[d])
4459 if (pos_d < limit0[d])
4461 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4462 cg_cm[cg], cm_new, pos_d);
4467 rvec_inc(cm_new, state->box[d]);
4470 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4471 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4473 for (k = k0; (k < k1); k++)
4475 rvec_inc(state->x[k], state->box[d]);
4478 rotate_state_atom(state, k);
4484 else if (d < npbcdim)
4486 /* Put the charge group in the rectangular unit-cell */
4487 while (cm_new[d] >= state->box[d][d])
4489 rvec_dec(cm_new, state->box[d]);
4490 for (k = k0; (k < k1); k++)
4492 rvec_dec(state->x[k], state->box[d]);
4495 while (cm_new[d] < 0)
4497 rvec_inc(cm_new, state->box[d]);
4498 for (k = k0; (k < k1); k++)
4500 rvec_inc(state->x[k], state->box[d]);
4506 copy_rvec(cm_new, cg_cm[cg]);
4508 /* Determine where this cg should go */
4511 for (d = 0; d < dd->ndim; d++)
4516 flag |= DD_FLAG_FW(d);
4522 else if (dev[dim] == -1)
4524 flag |= DD_FLAG_BW(d);
4527 if (dd->nc[dim] > 2)
4538 /* Temporarily store the flag in move */
4539 move[cg] = mc + flag;
4543 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4544 gmx_domdec_t *dd, ivec tric_dir,
4545 t_state *state, rvec **f,
4546 t_forcerec *fr, t_mdatoms *md,
4554 int ncg[DIM*2], nat[DIM*2];
4555 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4556 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4557 int sbuf[2], rbuf[2];
4558 int home_pos_cg, home_pos_at, buf_pos;
4560 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4563 real inv_ncg, pos_d;
4565 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4567 cginfo_mb_t *cginfo_mb;
4568 gmx_domdec_comm_t *comm;
4570 int nthread, thread;
4574 check_screw_box(state->box);
4578 if (fr->cutoff_scheme == ecutsGROUP)
4583 for (i = 0; i < estNR; i++)
4589 case estX: /* Always present */ break;
4590 case estV: bV = (state->flags & (1<<i)); break;
4591 case estSDX: bSDX = (state->flags & (1<<i)); break;
4592 case estCGP: bCGP = (state->flags & (1<<i)); break;
4595 case estDISRE_INITF:
4596 case estDISRE_RM3TAV:
4597 case estORIRE_INITF:
4599 /* No processing required */
4602 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4607 if (dd->ncg_tot > comm->nalloc_int)
4609 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4610 srenew(comm->buf_int, comm->nalloc_int);
4612 move = comm->buf_int;
4614 /* Clear the count */
4615 for (c = 0; c < dd->ndim*2; c++)
4621 npbcdim = dd->npbcdim;
4623 for (d = 0; (d < DIM); d++)
4625 limitd[d] = dd->comm->cellsize_min[d];
4626 if (d >= npbcdim && dd->ci[d] == 0)
4628 cell_x0[d] = -GMX_FLOAT_MAX;
4632 cell_x0[d] = comm->cell_x0[d];
4634 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4636 cell_x1[d] = GMX_FLOAT_MAX;
4640 cell_x1[d] = comm->cell_x1[d];
4644 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4645 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4649 /* We check after communication if a charge group moved
4650 * more than one cell. Set the pre-comm check limit to float_max.
4652 limit0[d] = -GMX_FLOAT_MAX;
4653 limit1[d] = GMX_FLOAT_MAX;
4657 make_tric_corr_matrix(npbcdim, state->box, tcm);
4659 cgindex = dd->cgindex;
4661 nthread = gmx_omp_nthreads_get(emntDomdec);
4663 /* Compute the center of geometry for all home charge groups
4664 * and put them in the box and determine where they should go.
4666 #pragma omp parallel for num_threads(nthread) schedule(static)
4667 for (thread = 0; thread < nthread; thread++)
4669 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4670 cell_x0, cell_x1, limitd, limit0, limit1,
4672 ( thread *dd->ncg_home)/nthread,
4673 ((thread+1)*dd->ncg_home)/nthread,
4674 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4678 for (cg = 0; cg < dd->ncg_home; cg++)
4683 flag = mc & ~DD_FLAG_NRCG;
4684 mc = mc & DD_FLAG_NRCG;
4687 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4689 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4690 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4692 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4693 /* We store the cg size in the lower 16 bits
4694 * and the place where the charge group should go
4695 * in the next 6 bits. This saves some communication volume.
4697 nrcg = cgindex[cg+1] - cgindex[cg];
4698 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4704 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4705 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4708 for (i = 0; i < dd->ndim*2; i++)
4710 *ncg_moved += ncg[i];
4727 /* Make sure the communication buffers are large enough */
4728 for (mc = 0; mc < dd->ndim*2; mc++)
4730 nvr = ncg[mc] + nat[mc]*nvec;
4731 if (nvr > comm->cgcm_state_nalloc[mc])
4733 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4734 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4738 switch (fr->cutoff_scheme)
4741 /* Recalculating cg_cm might be cheaper than communicating,
4742 * but that could give rise to rounding issues.
4745 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4746 nvec, cg_cm, comm, bCompact);
4749 /* Without charge groups we send the moved atom coordinates
4750 * over twice. This is so the code below can be used without
4751 * many conditionals for both for with and without charge groups.
4754 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4755 nvec, state->x, comm, FALSE);
4758 home_pos_cg -= *ncg_moved;
4762 gmx_incons("unimplemented");
4768 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4769 nvec, vec++, state->x, comm, bCompact);
4772 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4773 nvec, vec++, state->v, comm, bCompact);
4777 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4778 nvec, vec++, state->sd_X, comm, bCompact);
4782 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4783 nvec, vec++, state->cg_p, comm, bCompact);
4788 compact_ind(dd->ncg_home, move,
4789 dd->index_gl, dd->cgindex, dd->gatindex,
4790 dd->ga2la, comm->bLocalCG,
4795 if (fr->cutoff_scheme == ecutsVERLET)
4797 moved = get_moved(comm, dd->ncg_home);
4799 for (k = 0; k < dd->ncg_home; k++)
4806 moved = fr->ns.grid->cell_index;
4809 clear_and_mark_ind(dd->ncg_home, move,
4810 dd->index_gl, dd->cgindex, dd->gatindex,
4811 dd->ga2la, comm->bLocalCG,
4815 cginfo_mb = fr->cginfo_mb;
4817 *ncg_stay_home = home_pos_cg;
4818 for (d = 0; d < dd->ndim; d++)
4824 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4827 /* Communicate the cg and atom counts */
4832 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4833 d, dir, sbuf[0], sbuf[1]);
4835 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4837 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4839 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4840 srenew(comm->buf_int, comm->nalloc_int);
4843 /* Communicate the charge group indices, sizes and flags */
4844 dd_sendrecv_int(dd, d, dir,
4845 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4846 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4848 nvs = ncg[cdd] + nat[cdd]*nvec;
4849 i = rbuf[0] + rbuf[1] *nvec;
4850 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4852 /* Communicate cgcm and state */
4853 dd_sendrecv_rvec(dd, d, dir,
4854 comm->cgcm_state[cdd], nvs,
4855 comm->vbuf.v+nvr, i);
4856 ncg_recv += rbuf[0];
4857 nat_recv += rbuf[1];
4861 /* Process the received charge groups */
4863 for (cg = 0; cg < ncg_recv; cg++)
4865 flag = comm->buf_int[cg*DD_CGIBS+1];
4867 if (dim >= npbcdim && dd->nc[dim] > 2)
4869 /* No pbc in this dim and more than one domain boundary.
4870 * We do a separate check if a charge group didn't move too far.
4872 if (((flag & DD_FLAG_FW(d)) &&
4873 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4874 ((flag & DD_FLAG_BW(d)) &&
4875 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4877 cg_move_error(fplog, dd, step, cg, dim,
4878 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4880 comm->vbuf.v[buf_pos],
4881 comm->vbuf.v[buf_pos],
4882 comm->vbuf.v[buf_pos][dim]);
4889 /* Check which direction this cg should go */
4890 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4894 /* The cell boundaries for dimension d2 are not equal
4895 * for each cell row of the lower dimension(s),
4896 * therefore we might need to redetermine where
4897 * this cg should go.
4900 /* If this cg crosses the box boundary in dimension d2
4901 * we can use the communicated flag, so we do not
4902 * have to worry about pbc.
4904 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4905 (flag & DD_FLAG_FW(d2))) ||
4906 (dd->ci[dim2] == 0 &&
4907 (flag & DD_FLAG_BW(d2)))))
4909 /* Clear the two flags for this dimension */
4910 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4911 /* Determine the location of this cg
4912 * in lattice coordinates
4914 pos_d = comm->vbuf.v[buf_pos][dim2];
4917 for (d3 = dim2+1; d3 < DIM; d3++)
4920 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4923 /* Check of we are not at the box edge.
4924 * pbc is only handled in the first step above,
4925 * but this check could move over pbc while
4926 * the first step did not due to different rounding.
4928 if (pos_d >= cell_x1[dim2] &&
4929 dd->ci[dim2] != dd->nc[dim2]-1)
4931 flag |= DD_FLAG_FW(d2);
4933 else if (pos_d < cell_x0[dim2] &&
4936 flag |= DD_FLAG_BW(d2);
4938 comm->buf_int[cg*DD_CGIBS+1] = flag;
4941 /* Set to which neighboring cell this cg should go */
4942 if (flag & DD_FLAG_FW(d2))
4946 else if (flag & DD_FLAG_BW(d2))
4948 if (dd->nc[dd->dim[d2]] > 2)
4960 nrcg = flag & DD_FLAG_NRCG;
4963 if (home_pos_cg+1 > dd->cg_nalloc)
4965 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4966 srenew(dd->index_gl, dd->cg_nalloc);
4967 srenew(dd->cgindex, dd->cg_nalloc+1);
4969 /* Set the global charge group index and size */
4970 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4971 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4972 /* Copy the state from the buffer */
4973 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4974 if (fr->cutoff_scheme == ecutsGROUP)
4977 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4981 /* Set the cginfo */
4982 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4983 dd->index_gl[home_pos_cg]);
4986 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4989 if (home_pos_at+nrcg > state->nalloc)
4991 dd_realloc_state(state, f, home_pos_at+nrcg);
4993 for (i = 0; i < nrcg; i++)
4995 copy_rvec(comm->vbuf.v[buf_pos++],
4996 state->x[home_pos_at+i]);
5000 for (i = 0; i < nrcg; i++)
5002 copy_rvec(comm->vbuf.v[buf_pos++],
5003 state->v[home_pos_at+i]);
5008 for (i = 0; i < nrcg; i++)
5010 copy_rvec(comm->vbuf.v[buf_pos++],
5011 state->sd_X[home_pos_at+i]);
5016 for (i = 0; i < nrcg; i++)
5018 copy_rvec(comm->vbuf.v[buf_pos++],
5019 state->cg_p[home_pos_at+i]);
5023 home_pos_at += nrcg;
5027 /* Reallocate the buffers if necessary */
5028 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5030 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5031 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5033 nvr = ncg[mc] + nat[mc]*nvec;
5034 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5036 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5037 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5039 /* Copy from the receive to the send buffers */
5040 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5041 comm->buf_int + cg*DD_CGIBS,
5042 DD_CGIBS*sizeof(int));
5043 memcpy(comm->cgcm_state[mc][nvr],
5044 comm->vbuf.v[buf_pos],
5045 (1+nrcg*nvec)*sizeof(rvec));
5046 buf_pos += 1 + nrcg*nvec;
5053 /* With sorting (!bCompact) the indices are now only partially up to date
5054 * and ncg_home and nat_home are not the real count, since there are
5055 * "holes" in the arrays for the charge groups that moved to neighbors.
5057 if (fr->cutoff_scheme == ecutsVERLET)
5059 moved = get_moved(comm, home_pos_cg);
5061 for (i = dd->ncg_home; i < home_pos_cg; i++)
5066 dd->ncg_home = home_pos_cg;
5067 dd->nat_home = home_pos_at;
5072 "Finished repartitioning: cgs moved out %d, new home %d\n",
5073 *ncg_moved, dd->ncg_home-*ncg_moved);
5078 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5080 dd->comm->cycl[ddCycl] += cycles;
5081 dd->comm->cycl_n[ddCycl]++;
5082 if (cycles > dd->comm->cycl_max[ddCycl])
5084 dd->comm->cycl_max[ddCycl] = cycles;
5088 static double force_flop_count(t_nrnb *nrnb)
5095 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5097 /* To get closer to the real timings, we half the count
5098 * for the normal loops and again half it for water loops.
5101 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5103 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5107 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5110 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5113 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5115 sum += nrnb->n[i]*cost_nrnb(i);
5118 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5120 sum += nrnb->n[i]*cost_nrnb(i);
5126 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5128 if (dd->comm->eFlop)
5130 dd->comm->flop -= force_flop_count(nrnb);
5133 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5135 if (dd->comm->eFlop)
5137 dd->comm->flop += force_flop_count(nrnb);
5142 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5146 for (i = 0; i < ddCyclNr; i++)
5148 dd->comm->cycl[i] = 0;
5149 dd->comm->cycl_n[i] = 0;
5150 dd->comm->cycl_max[i] = 0;
5153 dd->comm->flop_n = 0;
5156 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5158 gmx_domdec_comm_t *comm;
5159 gmx_domdec_load_t *load;
5160 gmx_domdec_root_t *root = NULL;
5161 int d, dim, cid, i, pos;
5162 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5167 fprintf(debug, "get_load_distribution start\n");
5170 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5174 bSepPME = (dd->pme_nodeid >= 0);
5176 for (d = dd->ndim-1; d >= 0; d--)
5179 /* Check if we participate in the communication in this dimension */
5180 if (d == dd->ndim-1 ||
5181 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5183 load = &comm->load[d];
5186 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5189 if (d == dd->ndim-1)
5191 sbuf[pos++] = dd_force_load(comm);
5192 sbuf[pos++] = sbuf[0];
5195 sbuf[pos++] = sbuf[0];
5196 sbuf[pos++] = cell_frac;
5199 sbuf[pos++] = comm->cell_f_max0[d];
5200 sbuf[pos++] = comm->cell_f_min1[d];
5205 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5206 sbuf[pos++] = comm->cycl[ddCyclPME];
5211 sbuf[pos++] = comm->load[d+1].sum;
5212 sbuf[pos++] = comm->load[d+1].max;
5215 sbuf[pos++] = comm->load[d+1].sum_m;
5216 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5217 sbuf[pos++] = comm->load[d+1].flags;
5220 sbuf[pos++] = comm->cell_f_max0[d];
5221 sbuf[pos++] = comm->cell_f_min1[d];
5226 sbuf[pos++] = comm->load[d+1].mdf;
5227 sbuf[pos++] = comm->load[d+1].pme;
5231 /* Communicate a row in DD direction d.
5232 * The communicators are setup such that the root always has rank 0.
5235 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5236 load->load, load->nload*sizeof(float), MPI_BYTE,
5237 0, comm->mpi_comm_load[d]);
5239 if (dd->ci[dim] == dd->master_ci[dim])
5241 /* We are the root, process this row */
5242 if (comm->bDynLoadBal)
5244 root = comm->root[d];
5254 for (i = 0; i < dd->nc[dim]; i++)
5256 load->sum += load->load[pos++];
5257 load->max = max(load->max, load->load[pos]);
5263 /* This direction could not be load balanced properly,
5264 * therefore we need to use the maximum iso the average load.
5266 load->sum_m = max(load->sum_m, load->load[pos]);
5270 load->sum_m += load->load[pos];
5273 load->cvol_min = min(load->cvol_min, load->load[pos]);
5277 load->flags = (int)(load->load[pos++] + 0.5);
5281 root->cell_f_max0[i] = load->load[pos++];
5282 root->cell_f_min1[i] = load->load[pos++];
5287 load->mdf = max(load->mdf, load->load[pos]);
5289 load->pme = max(load->pme, load->load[pos]);
5293 if (comm->bDynLoadBal && root->bLimited)
5295 load->sum_m *= dd->nc[dim];
5296 load->flags |= (1<<d);
5304 comm->nload += dd_load_count(comm);
5305 comm->load_step += comm->cycl[ddCyclStep];
5306 comm->load_sum += comm->load[0].sum;
5307 comm->load_max += comm->load[0].max;
5308 if (comm->bDynLoadBal)
5310 for (d = 0; d < dd->ndim; d++)
5312 if (comm->load[0].flags & (1<<d))
5314 comm->load_lim[d]++;
5320 comm->load_mdf += comm->load[0].mdf;
5321 comm->load_pme += comm->load[0].pme;
5325 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5329 fprintf(debug, "get_load_distribution finished\n");
5333 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5335 /* Return the relative performance loss on the total run time
5336 * due to the force calculation load imbalance.
5338 if (dd->comm->nload > 0)
5341 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5342 (dd->comm->load_step*dd->nnodes);
5350 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5353 int npp, npme, nnodes, d, limp;
5354 float imbal, pme_f_ratio, lossf, lossp = 0;
5356 gmx_domdec_comm_t *comm;
5359 if (DDMASTER(dd) && comm->nload > 0)
5362 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5363 nnodes = npp + npme;
5364 imbal = comm->load_max*npp/comm->load_sum - 1;
5365 lossf = dd_force_imb_perf_loss(dd);
5366 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5367 fprintf(fplog, "%s", buf);
5368 fprintf(stderr, "\n");
5369 fprintf(stderr, "%s", buf);
5370 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5371 fprintf(fplog, "%s", buf);
5372 fprintf(stderr, "%s", buf);
5374 if (comm->bDynLoadBal)
5376 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5377 for (d = 0; d < dd->ndim; d++)
5379 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5380 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5386 sprintf(buf+strlen(buf), "\n");
5387 fprintf(fplog, "%s", buf);
5388 fprintf(stderr, "%s", buf);
5392 pme_f_ratio = comm->load_pme/comm->load_mdf;
5393 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5396 lossp *= (float)npme/(float)nnodes;
5400 lossp *= (float)npp/(float)nnodes;
5402 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5403 fprintf(fplog, "%s", buf);
5404 fprintf(stderr, "%s", buf);
5405 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5406 fprintf(fplog, "%s", buf);
5407 fprintf(stderr, "%s", buf);
5409 fprintf(fplog, "\n");
5410 fprintf(stderr, "\n");
5412 if (lossf >= DD_PERF_LOSS)
5415 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5416 " in the domain decomposition.\n", lossf*100);
5417 if (!comm->bDynLoadBal)
5419 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5423 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5425 fprintf(fplog, "%s\n", buf);
5426 fprintf(stderr, "%s\n", buf);
5428 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5431 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5432 " had %s work to do than the PP nodes.\n"
5433 " You might want to %s the number of PME nodes\n"
5434 " or %s the cut-off and the grid spacing.\n",
5436 (lossp < 0) ? "less" : "more",
5437 (lossp < 0) ? "decrease" : "increase",
5438 (lossp < 0) ? "decrease" : "increase");
5439 fprintf(fplog, "%s\n", buf);
5440 fprintf(stderr, "%s\n", buf);
5445 static float dd_vol_min(gmx_domdec_t *dd)
5447 return dd->comm->load[0].cvol_min*dd->nnodes;
5450 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5452 return dd->comm->load[0].flags;
5455 static float dd_f_imbal(gmx_domdec_t *dd)
5457 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5460 float dd_pme_f_ratio(gmx_domdec_t *dd)
5462 if (dd->comm->cycl_n[ddCyclPME] > 0)
5464 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5472 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5477 flags = dd_load_flags(dd);
5481 "DD load balancing is limited by minimum cell size in dimension");
5482 for (d = 0; d < dd->ndim; d++)
5486 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5489 fprintf(fplog, "\n");
5491 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5492 if (dd->comm->bDynLoadBal)
5494 fprintf(fplog, " vol min/aver %5.3f%c",
5495 dd_vol_min(dd), flags ? '!' : ' ');
5497 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5498 if (dd->comm->cycl_n[ddCyclPME])
5500 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5502 fprintf(fplog, "\n\n");
5505 static void dd_print_load_verbose(gmx_domdec_t *dd)
5507 if (dd->comm->bDynLoadBal)
5509 fprintf(stderr, "vol %4.2f%c ",
5510 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5512 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5513 if (dd->comm->cycl_n[ddCyclPME])
5515 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5520 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5525 gmx_domdec_root_t *root;
5526 gmx_bool bPartOfGroup = FALSE;
5528 dim = dd->dim[dim_ind];
5529 copy_ivec(loc, loc_c);
5530 for (i = 0; i < dd->nc[dim]; i++)
5533 rank = dd_index(dd->nc, loc_c);
5534 if (rank == dd->rank)
5536 /* This process is part of the group */
5537 bPartOfGroup = TRUE;
5540 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5544 dd->comm->mpi_comm_load[dim_ind] = c_row;
5545 if (dd->comm->eDLB != edlbNO)
5547 if (dd->ci[dim] == dd->master_ci[dim])
5549 /* This is the root process of this row */
5550 snew(dd->comm->root[dim_ind], 1);
5551 root = dd->comm->root[dim_ind];
5552 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5553 snew(root->old_cell_f, dd->nc[dim]+1);
5554 snew(root->bCellMin, dd->nc[dim]);
5557 snew(root->cell_f_max0, dd->nc[dim]);
5558 snew(root->cell_f_min1, dd->nc[dim]);
5559 snew(root->bound_min, dd->nc[dim]);
5560 snew(root->bound_max, dd->nc[dim]);
5562 snew(root->buf_ncd, dd->nc[dim]);
5566 /* This is not a root process, we only need to receive cell_f */
5567 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5570 if (dd->ci[dim] == dd->master_ci[dim])
5572 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5578 static void make_load_communicators(gmx_domdec_t *dd)
5581 int dim0, dim1, i, j;
5586 fprintf(debug, "Making load communicators\n");
5589 snew(dd->comm->load, dd->ndim);
5590 snew(dd->comm->mpi_comm_load, dd->ndim);
5593 make_load_communicator(dd, 0, loc);
5597 for (i = 0; i < dd->nc[dim0]; i++)
5600 make_load_communicator(dd, 1, loc);
5606 for (i = 0; i < dd->nc[dim0]; i++)
5610 for (j = 0; j < dd->nc[dim1]; j++)
5613 make_load_communicator(dd, 2, loc);
5620 fprintf(debug, "Finished making load communicators\n");
5625 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5628 int d, dim, i, j, m;
5631 ivec dd_zp[DD_MAXIZONE];
5632 gmx_domdec_zones_t *zones;
5633 gmx_domdec_ns_ranges_t *izone;
5635 for (d = 0; d < dd->ndim; d++)
5638 copy_ivec(dd->ci, tmp);
5639 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5640 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5641 copy_ivec(dd->ci, tmp);
5642 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5643 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5646 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5649 dd->neighbor[d][1]);
5655 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5657 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5658 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5665 for (i = 0; i < nzonep; i++)
5667 copy_ivec(dd_zp3[i], dd_zp[i]);
5673 for (i = 0; i < nzonep; i++)
5675 copy_ivec(dd_zp2[i], dd_zp[i]);
5681 for (i = 0; i < nzonep; i++)
5683 copy_ivec(dd_zp1[i], dd_zp[i]);
5687 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5692 zones = &dd->comm->zones;
5694 for (i = 0; i < nzone; i++)
5697 clear_ivec(zones->shift[i]);
5698 for (d = 0; d < dd->ndim; d++)
5700 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5705 for (i = 0; i < nzone; i++)
5707 for (d = 0; d < DIM; d++)
5709 s[d] = dd->ci[d] - zones->shift[i][d];
5714 else if (s[d] >= dd->nc[d])
5720 zones->nizone = nzonep;
5721 for (i = 0; i < zones->nizone; i++)
5723 if (dd_zp[i][0] != i)
5725 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5727 izone = &zones->izone[i];
5728 izone->j0 = dd_zp[i][1];
5729 izone->j1 = dd_zp[i][2];
5730 for (dim = 0; dim < DIM; dim++)
5732 if (dd->nc[dim] == 1)
5734 /* All shifts should be allowed */
5735 izone->shift0[dim] = -1;
5736 izone->shift1[dim] = 1;
5741 izone->shift0[d] = 0;
5742 izone->shift1[d] = 0;
5743 for(j=izone->j0; j<izone->j1; j++) {
5744 if (dd->shift[j][d] > dd->shift[i][d])
5745 izone->shift0[d] = -1;
5746 if (dd->shift[j][d] < dd->shift[i][d])
5747 izone->shift1[d] = 1;
5753 /* Assume the shift are not more than 1 cell */
5754 izone->shift0[dim] = 1;
5755 izone->shift1[dim] = -1;
5756 for (j = izone->j0; j < izone->j1; j++)
5758 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5759 if (shift_diff < izone->shift0[dim])
5761 izone->shift0[dim] = shift_diff;
5763 if (shift_diff > izone->shift1[dim])
5765 izone->shift1[dim] = shift_diff;
5772 if (dd->comm->eDLB != edlbNO)
5774 snew(dd->comm->root, dd->ndim);
5777 if (dd->comm->bRecordLoad)
5779 make_load_communicators(dd);
5783 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5786 gmx_domdec_comm_t *comm;
5797 if (comm->bCartesianPP)
5799 /* Set up cartesian communication for the particle-particle part */
5802 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5803 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5806 for (i = 0; i < DIM; i++)
5810 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5812 /* We overwrite the old communicator with the new cartesian one */
5813 cr->mpi_comm_mygroup = comm_cart;
5816 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5817 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5819 if (comm->bCartesianPP_PME)
5821 /* Since we want to use the original cartesian setup for sim,
5822 * and not the one after split, we need to make an index.
5824 snew(comm->ddindex2ddnodeid, dd->nnodes);
5825 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5826 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5827 /* Get the rank of the DD master,
5828 * above we made sure that the master node is a PP node.
5838 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5840 else if (comm->bCartesianPP)
5842 if (cr->npmenodes == 0)
5844 /* The PP communicator is also
5845 * the communicator for this simulation
5847 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5849 cr->nodeid = dd->rank;
5851 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5853 /* We need to make an index to go from the coordinates
5854 * to the nodeid of this simulation.
5856 snew(comm->ddindex2simnodeid, dd->nnodes);
5857 snew(buf, dd->nnodes);
5858 if (cr->duty & DUTY_PP)
5860 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5862 /* Communicate the ddindex to simulation nodeid index */
5863 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5864 cr->mpi_comm_mysim);
5867 /* Determine the master coordinates and rank.
5868 * The DD master should be the same node as the master of this sim.
5870 for (i = 0; i < dd->nnodes; i++)
5872 if (comm->ddindex2simnodeid[i] == 0)
5874 ddindex2xyz(dd->nc, i, dd->master_ci);
5875 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5880 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5885 /* No Cartesian communicators */
5886 /* We use the rank in dd->comm->all as DD index */
5887 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5888 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5890 clear_ivec(dd->master_ci);
5897 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5898 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5903 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5904 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5908 static void receive_ddindex2simnodeid(t_commrec *cr)
5912 gmx_domdec_comm_t *comm;
5919 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5921 snew(comm->ddindex2simnodeid, dd->nnodes);
5922 snew(buf, dd->nnodes);
5923 if (cr->duty & DUTY_PP)
5925 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5928 /* Communicate the ddindex to simulation nodeid index */
5929 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5930 cr->mpi_comm_mysim);
5937 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5938 int ncg, int natoms)
5940 gmx_domdec_master_t *ma;
5945 snew(ma->ncg, dd->nnodes);
5946 snew(ma->index, dd->nnodes+1);
5948 snew(ma->nat, dd->nnodes);
5949 snew(ma->ibuf, dd->nnodes*2);
5950 snew(ma->cell_x, DIM);
5951 for (i = 0; i < DIM; i++)
5953 snew(ma->cell_x[i], dd->nc[i]+1);
5956 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5962 snew(ma->vbuf, natoms);
5968 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5972 gmx_domdec_comm_t *comm;
5983 if (comm->bCartesianPP)
5985 for (i = 1; i < DIM; i++)
5987 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5989 if (bDiv[YY] || bDiv[ZZ])
5991 comm->bCartesianPP_PME = TRUE;
5992 /* If we have 2D PME decomposition, which is always in x+y,
5993 * we stack the PME only nodes in z.
5994 * Otherwise we choose the direction that provides the thinnest slab
5995 * of PME only nodes as this will have the least effect
5996 * on the PP communication.
5997 * But for the PME communication the opposite might be better.
5999 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6001 dd->nc[YY] > dd->nc[ZZ]))
6003 comm->cartpmedim = ZZ;
6007 comm->cartpmedim = YY;
6009 comm->ntot[comm->cartpmedim]
6010 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6014 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]);
6016 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6021 if (comm->bCartesianPP_PME)
6025 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]);
6028 for (i = 0; i < DIM; i++)
6032 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6035 MPI_Comm_rank(comm_cart, &rank);
6036 if (MASTERNODE(cr) && rank != 0)
6038 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6041 /* With this assigment we loose the link to the original communicator
6042 * which will usually be MPI_COMM_WORLD, unless have multisim.
6044 cr->mpi_comm_mysim = comm_cart;
6045 cr->sim_nodeid = rank;
6047 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6051 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6052 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6055 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6059 if (cr->npmenodes == 0 ||
6060 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6062 cr->duty = DUTY_PME;
6065 /* Split the sim communicator into PP and PME only nodes */
6066 MPI_Comm_split(cr->mpi_comm_mysim,
6068 dd_index(comm->ntot, dd->ci),
6069 &cr->mpi_comm_mygroup);
6073 switch (dd_node_order)
6078 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6081 case ddnoINTERLEAVE:
6082 /* Interleave the PP-only and PME-only nodes,
6083 * as on clusters with dual-core machines this will double
6084 * the communication bandwidth of the PME processes
6085 * and thus speed up the PP <-> PME and inter PME communication.
6089 fprintf(fplog, "Interleaving PP and PME nodes\n");
6091 comm->pmenodes = dd_pmenodes(cr);
6096 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6099 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6101 cr->duty = DUTY_PME;
6108 /* Split the sim communicator into PP and PME only nodes */
6109 MPI_Comm_split(cr->mpi_comm_mysim,
6112 &cr->mpi_comm_mygroup);
6113 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6119 fprintf(fplog, "This is a %s only node\n\n",
6120 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6124 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6127 gmx_domdec_comm_t *comm;
6133 copy_ivec(dd->nc, comm->ntot);
6135 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6136 comm->bCartesianPP_PME = FALSE;
6138 /* Reorder the nodes by default. This might change the MPI ranks.
6139 * Real reordering is only supported on very few architectures,
6140 * Blue Gene is one of them.
6142 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6144 if (cr->npmenodes > 0)
6146 /* Split the communicator into a PP and PME part */
6147 split_communicator(fplog, cr, dd_node_order, CartReorder);
6148 if (comm->bCartesianPP_PME)
6150 /* We (possibly) reordered the nodes in split_communicator,
6151 * so it is no longer required in make_pp_communicator.
6153 CartReorder = FALSE;
6158 /* All nodes do PP and PME */
6160 /* We do not require separate communicators */
6161 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6165 if (cr->duty & DUTY_PP)
6167 /* Copy or make a new PP communicator */
6168 make_pp_communicator(fplog, cr, CartReorder);
6172 receive_ddindex2simnodeid(cr);
6175 if (!(cr->duty & DUTY_PME))
6177 /* Set up the commnuication to our PME node */
6178 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6179 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6182 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6183 dd->pme_nodeid, dd->pme_receive_vir_ener);
6188 dd->pme_nodeid = -1;
6193 dd->ma = init_gmx_domdec_master_t(dd,
6195 comm->cgs_gl.index[comm->cgs_gl.nr]);
6199 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6201 real *slb_frac, tot;
6206 if (nc > 1 && size_string != NULL)
6210 fprintf(fplog, "Using static load balancing for the %s direction\n",
6215 for (i = 0; i < nc; i++)
6218 sscanf(size_string, "%lf%n", &dbl, &n);
6221 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6230 fprintf(fplog, "Relative cell sizes:");
6232 for (i = 0; i < nc; i++)
6237 fprintf(fplog, " %5.3f", slb_frac[i]);
6242 fprintf(fplog, "\n");
6249 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6252 gmx_mtop_ilistloop_t iloop;
6256 iloop = gmx_mtop_ilistloop_init(mtop);
6257 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6259 for (ftype = 0; ftype < F_NRE; ftype++)
6261 if ((interaction_function[ftype].flags & IF_BOND) &&
6264 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6272 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6278 val = getenv(env_var);
6281 if (sscanf(val, "%d", &nst) <= 0)
6287 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6295 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6299 fprintf(stderr, "\n%s\n", warn_string);
6303 fprintf(fplog, "\n%s\n", warn_string);
6307 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6308 t_inputrec *ir, FILE *fplog)
6310 if (ir->ePBC == epbcSCREW &&
6311 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6313 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6316 if (ir->ns_type == ensSIMPLE)
6318 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6321 if (ir->nstlist == 0)
6323 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6326 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6328 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6332 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6337 r = ddbox->box_size[XX];
6338 for (di = 0; di < dd->ndim; di++)
6341 /* Check using the initial average cell size */
6342 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6348 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6349 const char *dlb_opt, gmx_bool bRecordLoad,
6350 unsigned long Flags, t_inputrec *ir)
6358 case 'a': eDLB = edlbAUTO; break;
6359 case 'n': eDLB = edlbNO; break;
6360 case 'y': eDLB = edlbYES; break;
6361 default: gmx_incons("Unknown dlb_opt");
6364 if (Flags & MD_RERUN)
6369 if (!EI_DYNAMICS(ir->eI))
6371 if (eDLB == edlbYES)
6373 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6374 dd_warning(cr, fplog, buf);
6382 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6387 if (Flags & MD_REPRODUCIBLE)
6394 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6398 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6401 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6409 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6414 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6416 /* Decomposition order z,y,x */
6419 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6421 for (dim = DIM-1; dim >= 0; dim--)
6423 if (dd->nc[dim] > 1)
6425 dd->dim[dd->ndim++] = dim;
6431 /* Decomposition order x,y,z */
6432 for (dim = 0; dim < DIM; dim++)
6434 if (dd->nc[dim] > 1)
6436 dd->dim[dd->ndim++] = dim;
6442 static gmx_domdec_comm_t *init_dd_comm()
6444 gmx_domdec_comm_t *comm;
6448 snew(comm->cggl_flag, DIM*2);
6449 snew(comm->cgcm_state, DIM*2);
6450 for (i = 0; i < DIM*2; i++)
6452 comm->cggl_flag_nalloc[i] = 0;
6453 comm->cgcm_state_nalloc[i] = 0;
6456 comm->nalloc_int = 0;
6457 comm->buf_int = NULL;
6459 vec_rvec_init(&comm->vbuf);
6461 comm->n_load_have = 0;
6462 comm->n_load_collect = 0;
6464 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6466 comm->sum_nat[i] = 0;
6470 comm->load_step = 0;
6473 clear_ivec(comm->load_lim);
6480 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6481 unsigned long Flags,
6483 real comm_distance_min, real rconstr,
6484 const char *dlb_opt, real dlb_scale,
6485 const char *sizex, const char *sizey, const char *sizez,
6486 gmx_mtop_t *mtop, t_inputrec *ir,
6487 matrix box, rvec *x,
6489 int *npme_x, int *npme_y)
6492 gmx_domdec_comm_t *comm;
6495 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6502 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6507 dd->comm = init_dd_comm();
6509 snew(comm->cggl_flag, DIM*2);
6510 snew(comm->cgcm_state, DIM*2);
6512 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6513 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6515 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6516 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6517 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6518 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6519 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6520 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6521 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6522 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6524 dd->pme_recv_f_alloc = 0;
6525 dd->pme_recv_f_buf = NULL;
6527 if (dd->bSendRecv2 && fplog)
6529 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");
6535 fprintf(fplog, "Will load balance based on FLOP count\n");
6537 if (comm->eFlop > 1)
6539 srand(1+cr->nodeid);
6541 comm->bRecordLoad = TRUE;
6545 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6549 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6551 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6554 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6556 dd->bGridJump = comm->bDynLoadBal;
6557 comm->bPMELoadBalDLBLimits = FALSE;
6559 if (comm->nstSortCG)
6563 if (comm->nstSortCG == 1)
6565 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6569 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6573 snew(comm->sort, 1);
6579 fprintf(fplog, "Will not sort the charge groups\n");
6583 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6585 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6586 if (comm->bInterCGBondeds)
6588 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6592 comm->bInterCGMultiBody = FALSE;
6595 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6596 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6598 if (ir->rlistlong == 0)
6600 /* Set the cut-off to some very large value,
6601 * so we don't need if statements everywhere in the code.
6602 * We use sqrt, since the cut-off is squared in some places.
6604 comm->cutoff = GMX_CUTOFF_INF;
6608 comm->cutoff = ir->rlistlong;
6610 comm->cutoff_mbody = 0;
6612 comm->cellsize_limit = 0;
6613 comm->bBondComm = FALSE;
6615 if (comm->bInterCGBondeds)
6617 if (comm_distance_min > 0)
6619 comm->cutoff_mbody = comm_distance_min;
6620 if (Flags & MD_DDBONDCOMM)
6622 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6626 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6628 r_bonded_limit = comm->cutoff_mbody;
6630 else if (ir->bPeriodicMols)
6632 /* Can not easily determine the required cut-off */
6633 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");
6634 comm->cutoff_mbody = comm->cutoff/2;
6635 r_bonded_limit = comm->cutoff_mbody;
6641 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6642 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6644 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6645 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6647 /* We use an initial margin of 10% for the minimum cell size,
6648 * except when we are just below the non-bonded cut-off.
6650 if (Flags & MD_DDBONDCOMM)
6652 if (max(r_2b, r_mb) > comm->cutoff)
6654 r_bonded = max(r_2b, r_mb);
6655 r_bonded_limit = 1.1*r_bonded;
6656 comm->bBondComm = TRUE;
6661 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6663 /* We determine cutoff_mbody later */
6667 /* No special bonded communication,
6668 * simply increase the DD cut-off.
6670 r_bonded_limit = 1.1*max(r_2b, r_mb);
6671 comm->cutoff_mbody = r_bonded_limit;
6672 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6675 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6679 "Minimum cell size due to bonded interactions: %.3f nm\n",
6680 comm->cellsize_limit);
6684 if (dd->bInterCGcons && rconstr <= 0)
6686 /* There is a cell size limit due to the constraints (P-LINCS) */
6687 rconstr = constr_r_max(fplog, mtop, ir);
6691 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6693 if (rconstr > comm->cellsize_limit)
6695 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6699 else if (rconstr > 0 && fplog)
6701 /* Here we do not check for dd->bInterCGcons,
6702 * because one can also set a cell size limit for virtual sites only
6703 * and at this point we don't know yet if there are intercg v-sites.
6706 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6709 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6711 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6715 copy_ivec(nc, dd->nc);
6716 set_dd_dim(fplog, dd);
6717 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6719 if (cr->npmenodes == -1)
6723 acs = average_cellsize_min(dd, ddbox);
6724 if (acs < comm->cellsize_limit)
6728 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6730 gmx_fatal_collective(FARGS, cr, NULL,
6731 "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",
6732 acs, comm->cellsize_limit);
6737 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6739 /* We need to choose the optimal DD grid and possibly PME nodes */
6740 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6741 comm->eDLB != edlbNO, dlb_scale,
6742 comm->cellsize_limit, comm->cutoff,
6743 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6745 if (dd->nc[XX] == 0)
6747 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6748 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6749 !bC ? "-rdd" : "-rcon",
6750 comm->eDLB != edlbNO ? " or -dds" : "",
6751 bC ? " or your LINCS settings" : "");
6753 gmx_fatal_collective(FARGS, cr, NULL,
6754 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6756 "Look in the log file for details on the domain decomposition",
6757 cr->nnodes-cr->npmenodes, limit, buf);
6759 set_dd_dim(fplog, dd);
6765 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6766 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6769 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6770 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6772 gmx_fatal_collective(FARGS, cr, NULL,
6773 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6774 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6776 if (cr->npmenodes > dd->nnodes)
6778 gmx_fatal_collective(FARGS, cr, NULL,
6779 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6781 if (cr->npmenodes > 0)
6783 comm->npmenodes = cr->npmenodes;
6787 comm->npmenodes = dd->nnodes;
6790 if (EEL_PME(ir->coulombtype))
6792 /* The following choices should match those
6793 * in comm_cost_est in domdec_setup.c.
6794 * Note that here the checks have to take into account
6795 * that the decomposition might occur in a different order than xyz
6796 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6797 * in which case they will not match those in comm_cost_est,
6798 * but since that is mainly for testing purposes that's fine.
6800 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6801 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6802 getenv("GMX_PMEONEDD") == NULL)
6804 comm->npmedecompdim = 2;
6805 comm->npmenodes_x = dd->nc[XX];
6806 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6810 /* In case nc is 1 in both x and y we could still choose to
6811 * decompose pme in y instead of x, but we use x for simplicity.
6813 comm->npmedecompdim = 1;
6814 if (dd->dim[0] == YY)
6816 comm->npmenodes_x = 1;
6817 comm->npmenodes_y = comm->npmenodes;
6821 comm->npmenodes_x = comm->npmenodes;
6822 comm->npmenodes_y = 1;
6827 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6828 comm->npmenodes_x, comm->npmenodes_y, 1);
6833 comm->npmedecompdim = 0;
6834 comm->npmenodes_x = 0;
6835 comm->npmenodes_y = 0;
6838 /* Technically we don't need both of these,
6839 * but it simplifies code not having to recalculate it.
6841 *npme_x = comm->npmenodes_x;
6842 *npme_y = comm->npmenodes_y;
6844 snew(comm->slb_frac, DIM);
6845 if (comm->eDLB == edlbNO)
6847 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6848 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6849 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6852 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6854 if (comm->bBondComm || comm->eDLB != edlbNO)
6856 /* Set the bonded communication distance to halfway
6857 * the minimum and the maximum,
6858 * since the extra communication cost is nearly zero.
6860 acs = average_cellsize_min(dd, ddbox);
6861 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6862 if (comm->eDLB != edlbNO)
6864 /* Check if this does not limit the scaling */
6865 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6867 if (!comm->bBondComm)
6869 /* Without bBondComm do not go beyond the n.b. cut-off */
6870 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6871 if (comm->cellsize_limit >= comm->cutoff)
6873 /* We don't loose a lot of efficieny
6874 * when increasing it to the n.b. cut-off.
6875 * It can even be slightly faster, because we need
6876 * less checks for the communication setup.
6878 comm->cutoff_mbody = comm->cutoff;
6881 /* Check if we did not end up below our original limit */
6882 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6884 if (comm->cutoff_mbody > comm->cellsize_limit)
6886 comm->cellsize_limit = comm->cutoff_mbody;
6889 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6894 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6895 "cellsize limit %f\n",
6896 comm->bBondComm, comm->cellsize_limit);
6901 check_dd_restrictions(cr, dd, ir, fplog);
6904 comm->partition_step = INT_MIN;
6907 clear_dd_cycle_counts(dd);
6912 static void set_dlb_limits(gmx_domdec_t *dd)
6917 for (d = 0; d < dd->ndim; d++)
6919 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6920 dd->comm->cellsize_min[dd->dim[d]] =
6921 dd->comm->cellsize_min_dlb[dd->dim[d]];
6926 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6929 gmx_domdec_comm_t *comm;
6939 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);
6942 cellsize_min = comm->cellsize_min[dd->dim[0]];
6943 for (d = 1; d < dd->ndim; d++)
6945 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6948 if (cellsize_min < comm->cellsize_limit*1.05)
6950 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");
6952 /* Change DLB from "auto" to "no". */
6953 comm->eDLB = edlbNO;
6958 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6959 comm->bDynLoadBal = TRUE;
6960 dd->bGridJump = TRUE;
6964 /* We can set the required cell size info here,
6965 * so we do not need to communicate this.
6966 * The grid is completely uniform.
6968 for (d = 0; d < dd->ndim; d++)
6972 comm->load[d].sum_m = comm->load[d].sum;
6974 nc = dd->nc[dd->dim[d]];
6975 for (i = 0; i < nc; i++)
6977 comm->root[d]->cell_f[i] = i/(real)nc;
6980 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6981 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6984 comm->root[d]->cell_f[nc] = 1.0;
6989 static char *init_bLocalCG(gmx_mtop_t *mtop)
6994 ncg = ncg_mtop(mtop);
6995 snew(bLocalCG, ncg);
6996 for (cg = 0; cg < ncg; cg++)
6998 bLocalCG[cg] = FALSE;
7004 void dd_init_bondeds(FILE *fplog,
7005 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7006 gmx_vsite_t *vsite, gmx_constr_t constr,
7007 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7009 gmx_domdec_comm_t *comm;
7013 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7017 if (comm->bBondComm)
7019 /* Communicate atoms beyond the cut-off for bonded interactions */
7022 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7024 comm->bLocalCG = init_bLocalCG(mtop);
7028 /* Only communicate atoms based on cut-off */
7029 comm->cglink = NULL;
7030 comm->bLocalCG = NULL;
7034 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7036 gmx_bool bDynLoadBal, real dlb_scale,
7039 gmx_domdec_comm_t *comm;
7054 fprintf(fplog, "The maximum number of communication pulses is:");
7055 for (d = 0; d < dd->ndim; d++)
7057 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7059 fprintf(fplog, "\n");
7060 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7061 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7062 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7063 for (d = 0; d < DIM; d++)
7067 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7074 comm->cellsize_min_dlb[d]/
7075 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7077 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7080 fprintf(fplog, "\n");
7084 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7085 fprintf(fplog, "The initial number of communication pulses is:");
7086 for (d = 0; d < dd->ndim; d++)
7088 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7090 fprintf(fplog, "\n");
7091 fprintf(fplog, "The initial domain decomposition cell size is:");
7092 for (d = 0; d < DIM; d++)
7096 fprintf(fplog, " %c %.2f nm",
7097 dim2char(d), dd->comm->cellsize_min[d]);
7100 fprintf(fplog, "\n\n");
7103 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7105 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7106 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7107 "non-bonded interactions", "", comm->cutoff);
7111 limit = dd->comm->cellsize_limit;
7115 if (dynamic_dd_box(ddbox, ir))
7117 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7119 limit = dd->comm->cellsize_min[XX];
7120 for (d = 1; d < DIM; d++)
7122 limit = min(limit, dd->comm->cellsize_min[d]);
7126 if (comm->bInterCGBondeds)
7128 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7129 "two-body bonded interactions", "(-rdd)",
7130 max(comm->cutoff, comm->cutoff_mbody));
7131 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7132 "multi-body bonded interactions", "(-rdd)",
7133 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7137 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7138 "virtual site constructions", "(-rcon)", limit);
7140 if (dd->constraint_comm)
7142 sprintf(buf, "atoms separated by up to %d constraints",
7144 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7145 buf, "(-rcon)", limit);
7147 fprintf(fplog, "\n");
7153 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7155 const t_inputrec *ir,
7156 const gmx_ddbox_t *ddbox)
7158 gmx_domdec_comm_t *comm;
7159 int d, dim, npulse, npulse_d_max, npulse_d;
7164 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7166 /* Determine the maximum number of comm. pulses in one dimension */
7168 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7170 /* Determine the maximum required number of grid pulses */
7171 if (comm->cellsize_limit >= comm->cutoff)
7173 /* Only a single pulse is required */
7176 else if (!bNoCutOff && comm->cellsize_limit > 0)
7178 /* We round down slightly here to avoid overhead due to the latency
7179 * of extra communication calls when the cut-off
7180 * would be only slightly longer than the cell size.
7181 * Later cellsize_limit is redetermined,
7182 * so we can not miss interactions due to this rounding.
7184 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7188 /* There is no cell size limit */
7189 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7192 if (!bNoCutOff && npulse > 1)
7194 /* See if we can do with less pulses, based on dlb_scale */
7196 for (d = 0; d < dd->ndim; d++)
7199 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7200 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7201 npulse_d_max = max(npulse_d_max, npulse_d);
7203 npulse = min(npulse, npulse_d_max);
7206 /* This env var can override npulse */
7207 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7214 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7215 for (d = 0; d < dd->ndim; d++)
7217 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7218 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7219 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7220 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7221 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7223 comm->bVacDLBNoLimit = FALSE;
7227 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7228 if (!comm->bVacDLBNoLimit)
7230 comm->cellsize_limit = max(comm->cellsize_limit,
7231 comm->cutoff/comm->maxpulse);
7233 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7234 /* Set the minimum cell size for each DD dimension */
7235 for (d = 0; d < dd->ndim; d++)
7237 if (comm->bVacDLBNoLimit ||
7238 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7240 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7244 comm->cellsize_min_dlb[dd->dim[d]] =
7245 comm->cutoff/comm->cd[d].np_dlb;
7248 if (comm->cutoff_mbody <= 0)
7250 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7252 if (comm->bDynLoadBal)
7258 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7260 /* If each molecule is a single charge group
7261 * or we use domain decomposition for each periodic dimension,
7262 * we do not need to take pbc into account for the bonded interactions.
7264 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7267 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7270 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7271 t_inputrec *ir, t_forcerec *fr,
7274 gmx_domdec_comm_t *comm;
7280 /* Initialize the thread data.
7281 * This can not be done in init_domain_decomposition,
7282 * as the numbers of threads is determined later.
7284 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7287 snew(comm->dth, comm->nth);
7290 if (EEL_PME(ir->coulombtype))
7292 init_ddpme(dd, &comm->ddpme[0], 0);
7293 if (comm->npmedecompdim >= 2)
7295 init_ddpme(dd, &comm->ddpme[1], 1);
7300 comm->npmenodes = 0;
7301 if (dd->pme_nodeid >= 0)
7303 gmx_fatal_collective(FARGS, NULL, dd,
7304 "Can not have separate PME nodes without PME electrostatics");
7310 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7312 if (comm->eDLB != edlbNO)
7314 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7317 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7318 if (comm->eDLB == edlbAUTO)
7322 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7324 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7327 if (ir->ePBC == epbcNONE)
7329 vol_frac = 1 - 1/(double)dd->nnodes;
7334 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7338 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7340 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7342 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7345 static gmx_bool test_dd_cutoff(t_commrec *cr,
7346 t_state *state, t_inputrec *ir,
7357 set_ddbox(dd, FALSE, cr, ir, state->box,
7358 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7362 for (d = 0; d < dd->ndim; d++)
7366 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7367 if (dynamic_dd_box(&ddbox, ir))
7369 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7372 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7374 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7375 dd->comm->cd[d].np_dlb > 0)
7377 if (np > dd->comm->cd[d].np_dlb)
7382 /* If a current local cell size is smaller than the requested
7383 * cut-off, we could still fix it, but this gets very complicated.
7384 * Without fixing here, we might actually need more checks.
7386 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7393 if (dd->comm->eDLB != edlbNO)
7395 /* If DLB is not active yet, we don't need to check the grid jumps.
7396 * Actually we shouldn't, because then the grid jump data is not set.
7398 if (dd->comm->bDynLoadBal &&
7399 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7404 gmx_sumi(1, &LocallyLimited, cr);
7406 if (LocallyLimited > 0)
7415 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7418 gmx_bool bCutoffAllowed;
7420 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7424 cr->dd->comm->cutoff = cutoff_req;
7427 return bCutoffAllowed;
7430 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7432 gmx_domdec_comm_t *comm;
7434 comm = cr->dd->comm;
7436 /* Turn on the DLB limiting (might have been on already) */
7437 comm->bPMELoadBalDLBLimits = TRUE;
7439 /* Change the cut-off limit */
7440 comm->PMELoadBal_max_cutoff = comm->cutoff;
7443 static void merge_cg_buffers(int ncell,
7444 gmx_domdec_comm_dim_t *cd, int pulse,
7446 int *index_gl, int *recv_i,
7447 rvec *cg_cm, rvec *recv_vr,
7449 cginfo_mb_t *cginfo_mb, int *cginfo)
7451 gmx_domdec_ind_t *ind, *ind_p;
7452 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7453 int shift, shift_at;
7455 ind = &cd->ind[pulse];
7457 /* First correct the already stored data */
7458 shift = ind->nrecv[ncell];
7459 for (cell = ncell-1; cell >= 0; cell--)
7461 shift -= ind->nrecv[cell];
7464 /* Move the cg's present from previous grid pulses */
7465 cg0 = ncg_cell[ncell+cell];
7466 cg1 = ncg_cell[ncell+cell+1];
7467 cgindex[cg1+shift] = cgindex[cg1];
7468 for (cg = cg1-1; cg >= cg0; cg--)
7470 index_gl[cg+shift] = index_gl[cg];
7471 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7472 cgindex[cg+shift] = cgindex[cg];
7473 cginfo[cg+shift] = cginfo[cg];
7475 /* Correct the already stored send indices for the shift */
7476 for (p = 1; p <= pulse; p++)
7478 ind_p = &cd->ind[p];
7480 for (c = 0; c < cell; c++)
7482 cg0 += ind_p->nsend[c];
7484 cg1 = cg0 + ind_p->nsend[cell];
7485 for (cg = cg0; cg < cg1; cg++)
7487 ind_p->index[cg] += shift;
7493 /* Merge in the communicated buffers */
7497 for (cell = 0; cell < ncell; cell++)
7499 cg1 = ncg_cell[ncell+cell+1] + shift;
7502 /* Correct the old cg indices */
7503 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7505 cgindex[cg+1] += shift_at;
7508 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7510 /* Copy this charge group from the buffer */
7511 index_gl[cg1] = recv_i[cg0];
7512 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7513 /* Add it to the cgindex */
7514 cg_gl = index_gl[cg1];
7515 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7516 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7517 cgindex[cg1+1] = cgindex[cg1] + nat;
7522 shift += ind->nrecv[cell];
7523 ncg_cell[ncell+cell+1] = cg1;
7527 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7528 int nzone, int cg0, const int *cgindex)
7532 /* Store the atom block boundaries for easy copying of communication buffers
7535 for (zone = 0; zone < nzone; zone++)
7537 for (p = 0; p < cd->np; p++)
7539 cd->ind[p].cell2at0[zone] = cgindex[cg];
7540 cg += cd->ind[p].nrecv[zone];
7541 cd->ind[p].cell2at1[zone] = cgindex[cg];
7546 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7552 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7554 if (!bLocalCG[link->a[i]])
7563 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7565 real c[DIM][4]; /* the corners for the non-bonded communication */
7566 real cr0; /* corner for rounding */
7567 real cr1[4]; /* corners for rounding */
7568 real bc[DIM]; /* corners for bounded communication */
7569 real bcr1; /* corner for rounding for bonded communication */
7572 /* Determine the corners of the domain(s) we are communicating with */
7574 set_dd_corners(const gmx_domdec_t *dd,
7575 int dim0, int dim1, int dim2,
7579 const gmx_domdec_comm_t *comm;
7580 const gmx_domdec_zones_t *zones;
7585 zones = &comm->zones;
7587 /* Keep the compiler happy */
7591 /* The first dimension is equal for all cells */
7592 c->c[0][0] = comm->cell_x0[dim0];
7595 c->bc[0] = c->c[0][0];
7600 /* This cell row is only seen from the first row */
7601 c->c[1][0] = comm->cell_x0[dim1];
7602 /* All rows can see this row */
7603 c->c[1][1] = comm->cell_x0[dim1];
7606 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7609 /* For the multi-body distance we need the maximum */
7610 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7613 /* Set the upper-right corner for rounding */
7614 c->cr0 = comm->cell_x1[dim0];
7619 for (j = 0; j < 4; j++)
7621 c->c[2][j] = comm->cell_x0[dim2];
7625 /* Use the maximum of the i-cells that see a j-cell */
7626 for (i = 0; i < zones->nizone; i++)
7628 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7634 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7640 /* For the multi-body distance we need the maximum */
7641 c->bc[2] = comm->cell_x0[dim2];
7642 for (i = 0; i < 2; i++)
7644 for (j = 0; j < 2; j++)
7646 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7652 /* Set the upper-right corner for rounding */
7653 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7654 * Only cell (0,0,0) can see cell 7 (1,1,1)
7656 c->cr1[0] = comm->cell_x1[dim1];
7657 c->cr1[3] = comm->cell_x1[dim1];
7660 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7663 /* For the multi-body distance we need the maximum */
7664 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7671 /* Determine which cg's we need to send in this pulse from this zone */
7673 get_zone_pulse_cgs(gmx_domdec_t *dd,
7674 int zonei, int zone,
7676 const int *index_gl,
7678 int dim, int dim_ind,
7679 int dim0, int dim1, int dim2,
7680 real r_comm2, real r_bcomm2,
7684 real skew_fac2_d, real skew_fac_01,
7685 rvec *v_d, rvec *v_0, rvec *v_1,
7686 const dd_corners_t *c,
7688 gmx_bool bDistBonded,
7694 gmx_domdec_ind_t *ind,
7695 int **ibuf, int *ibuf_nalloc,
7701 gmx_domdec_comm_t *comm;
7703 gmx_bool bDistMB_pulse;
7705 real r2, rb2, r, tric_sh;
7708 int nsend_z, nsend, nat;
7712 bScrew = (dd->bScrewPBC && dim == XX);
7714 bDistMB_pulse = (bDistMB && bDistBonded);
7720 for (cg = cg0; cg < cg1; cg++)
7724 if (tric_dist[dim_ind] == 0)
7726 /* Rectangular direction, easy */
7727 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7734 r = cg_cm[cg][dim] - c->bc[dim_ind];
7740 /* Rounding gives at most a 16% reduction
7741 * in communicated atoms
7743 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7745 r = cg_cm[cg][dim0] - c->cr0;
7746 /* This is the first dimension, so always r >= 0 */
7753 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7755 r = cg_cm[cg][dim1] - c->cr1[zone];
7762 r = cg_cm[cg][dim1] - c->bcr1;
7772 /* Triclinic direction, more complicated */
7775 /* Rounding, conservative as the skew_fac multiplication
7776 * will slightly underestimate the distance.
7778 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7780 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7781 for (i = dim0+1; i < DIM; i++)
7783 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7785 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7788 rb[dim0] = rn[dim0];
7791 /* Take care that the cell planes along dim0 might not
7792 * be orthogonal to those along dim1 and dim2.
7794 for (i = 1; i <= dim_ind; i++)
7797 if (normal[dim0][dimd] > 0)
7799 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7802 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7807 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7809 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7811 for (i = dim1+1; i < DIM; i++)
7813 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7815 rn[dim1] += tric_sh;
7818 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7819 /* Take care of coupling of the distances
7820 * to the planes along dim0 and dim1 through dim2.
7822 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7823 /* Take care that the cell planes along dim1
7824 * might not be orthogonal to that along dim2.
7826 if (normal[dim1][dim2] > 0)
7828 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7834 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7837 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7838 /* Take care of coupling of the distances
7839 * to the planes along dim0 and dim1 through dim2.
7841 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7842 /* Take care that the cell planes along dim1
7843 * might not be orthogonal to that along dim2.
7845 if (normal[dim1][dim2] > 0)
7847 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7852 /* The distance along the communication direction */
7853 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7855 for (i = dim+1; i < DIM; i++)
7857 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7862 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7863 /* Take care of coupling of the distances
7864 * to the planes along dim0 and dim1 through dim2.
7866 if (dim_ind == 1 && zonei == 1)
7868 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7874 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7877 rb2 += rb[dim]*rb[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 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7891 ((bDistMB && rb2 < r_bcomm2) ||
7892 (bDist2B && r2 < r_bcomm2)) &&
7894 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7895 missing_link(comm->cglink, index_gl[cg],
7898 /* Make an index to the local charge groups */
7899 if (nsend+1 > ind->nalloc)
7901 ind->nalloc = over_alloc_large(nsend+1);
7902 srenew(ind->index, ind->nalloc);
7904 if (nsend+1 > *ibuf_nalloc)
7906 *ibuf_nalloc = over_alloc_large(nsend+1);
7907 srenew(*ibuf, *ibuf_nalloc);
7909 ind->index[nsend] = cg;
7910 (*ibuf)[nsend] = index_gl[cg];
7912 vec_rvec_check_alloc(vbuf, nsend+1);
7914 if (dd->ci[dim] == 0)
7916 /* Correct cg_cm for pbc */
7917 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7920 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7921 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7926 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7929 nat += cgindex[cg+1] - cgindex[cg];
7935 *nsend_z_ptr = nsend_z;
7938 static void setup_dd_communication(gmx_domdec_t *dd,
7939 matrix box, gmx_ddbox_t *ddbox,
7940 t_forcerec *fr, t_state *state, rvec **f)
7942 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7943 int nzone, nzone_send, zone, zonei, cg0, cg1;
7944 int c, i, j, cg, cg_gl, nrcg;
7945 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7946 gmx_domdec_comm_t *comm;
7947 gmx_domdec_zones_t *zones;
7948 gmx_domdec_comm_dim_t *cd;
7949 gmx_domdec_ind_t *ind;
7950 cginfo_mb_t *cginfo_mb;
7951 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7952 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7953 dd_corners_t corners;
7955 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7956 real skew_fac2_d, skew_fac_01;
7963 fprintf(debug, "Setting up DD communication\n");
7968 switch (fr->cutoff_scheme)
7977 gmx_incons("unimplemented");
7981 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7983 dim = dd->dim[dim_ind];
7985 /* Check if we need to use triclinic distances */
7986 tric_dist[dim_ind] = 0;
7987 for (i = 0; i <= dim_ind; i++)
7989 if (ddbox->tric_dir[dd->dim[i]])
7991 tric_dist[dim_ind] = 1;
7996 bBondComm = comm->bBondComm;
7998 /* Do we need to determine extra distances for multi-body bondeds? */
7999 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8001 /* Do we need to determine extra distances for only two-body bondeds? */
8002 bDist2B = (bBondComm && !bDistMB);
8004 r_comm2 = sqr(comm->cutoff);
8005 r_bcomm2 = sqr(comm->cutoff_mbody);
8009 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8012 zones = &comm->zones;
8015 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8016 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8018 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8020 /* Triclinic stuff */
8021 normal = ddbox->normal;
8025 v_0 = ddbox->v[dim0];
8026 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8028 /* Determine the coupling coefficient for the distances
8029 * to the cell planes along dim0 and dim1 through dim2.
8030 * This is required for correct rounding.
8033 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8036 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8042 v_1 = ddbox->v[dim1];
8045 zone_cg_range = zones->cg_range;
8046 index_gl = dd->index_gl;
8047 cgindex = dd->cgindex;
8048 cginfo_mb = fr->cginfo_mb;
8050 zone_cg_range[0] = 0;
8051 zone_cg_range[1] = dd->ncg_home;
8052 comm->zone_ncg1[0] = dd->ncg_home;
8053 pos_cg = dd->ncg_home;
8055 nat_tot = dd->nat_home;
8057 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8059 dim = dd->dim[dim_ind];
8060 cd = &comm->cd[dim_ind];
8062 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8064 /* No pbc in this dimension, the first node should not comm. */
8072 v_d = ddbox->v[dim];
8073 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8075 cd->bInPlace = TRUE;
8076 for (p = 0; p < cd->np; p++)
8078 /* Only atoms communicated in the first pulse are used
8079 * for multi-body bonded interactions or for bBondComm.
8081 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8086 for (zone = 0; zone < nzone_send; zone++)
8088 if (tric_dist[dim_ind] && dim_ind > 0)
8090 /* Determine slightly more optimized skew_fac's
8092 * This reduces the number of communicated atoms
8093 * by about 10% for 3D DD of rhombic dodecahedra.
8095 for (dimd = 0; dimd < dim; dimd++)
8097 sf2_round[dimd] = 1;
8098 if (ddbox->tric_dir[dimd])
8100 for (i = dd->dim[dimd]+1; i < DIM; i++)
8102 /* If we are shifted in dimension i
8103 * and the cell plane is tilted forward
8104 * in dimension i, skip this coupling.
8106 if (!(zones->shift[nzone+zone][i] &&
8107 ddbox->v[dimd][i][dimd] >= 0))
8110 sqr(ddbox->v[dimd][i][dimd]);
8113 sf2_round[dimd] = 1/sf2_round[dimd];
8118 zonei = zone_perm[dim_ind][zone];
8121 /* Here we permutate the zones to obtain a convenient order
8122 * for neighbor searching
8124 cg0 = zone_cg_range[zonei];
8125 cg1 = zone_cg_range[zonei+1];
8129 /* Look only at the cg's received in the previous grid pulse
8131 cg1 = zone_cg_range[nzone+zone+1];
8132 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8135 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8136 for (th = 0; th < comm->nth; th++)
8138 gmx_domdec_ind_t *ind_p;
8139 int **ibuf_p, *ibuf_nalloc_p;
8141 int *nsend_p, *nat_p;
8147 /* Thread 0 writes in the comm buffers */
8149 ibuf_p = &comm->buf_int;
8150 ibuf_nalloc_p = &comm->nalloc_int;
8151 vbuf_p = &comm->vbuf;
8154 nsend_zone_p = &ind->nsend[zone];
8158 /* Other threads write into temp buffers */
8159 ind_p = &comm->dth[th].ind;
8160 ibuf_p = &comm->dth[th].ibuf;
8161 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8162 vbuf_p = &comm->dth[th].vbuf;
8163 nsend_p = &comm->dth[th].nsend;
8164 nat_p = &comm->dth[th].nat;
8165 nsend_zone_p = &comm->dth[th].nsend_zone;
8167 comm->dth[th].nsend = 0;
8168 comm->dth[th].nat = 0;
8169 comm->dth[th].nsend_zone = 0;
8179 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8180 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8183 /* Get the cg's for this pulse in this zone */
8184 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8186 dim, dim_ind, dim0, dim1, dim2,
8189 normal, skew_fac2_d, skew_fac_01,
8190 v_d, v_0, v_1, &corners, sf2_round,
8191 bDistBonded, bBondComm,
8195 ibuf_p, ibuf_nalloc_p,
8201 /* Append data of threads>=1 to the communication buffers */
8202 for (th = 1; th < comm->nth; th++)
8204 dd_comm_setup_work_t *dth;
8207 dth = &comm->dth[th];
8209 ns1 = nsend + dth->nsend_zone;
8210 if (ns1 > ind->nalloc)
8212 ind->nalloc = over_alloc_dd(ns1);
8213 srenew(ind->index, ind->nalloc);
8215 if (ns1 > comm->nalloc_int)
8217 comm->nalloc_int = over_alloc_dd(ns1);
8218 srenew(comm->buf_int, comm->nalloc_int);
8220 if (ns1 > comm->vbuf.nalloc)
8222 comm->vbuf.nalloc = over_alloc_dd(ns1);
8223 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8226 for (i = 0; i < dth->nsend_zone; i++)
8228 ind->index[nsend] = dth->ind.index[i];
8229 comm->buf_int[nsend] = dth->ibuf[i];
8230 copy_rvec(dth->vbuf.v[i],
8231 comm->vbuf.v[nsend]);
8235 ind->nsend[zone] += dth->nsend_zone;
8238 /* Clear the counts in case we do not have pbc */
8239 for (zone = nzone_send; zone < nzone; zone++)
8241 ind->nsend[zone] = 0;
8243 ind->nsend[nzone] = nsend;
8244 ind->nsend[nzone+1] = nat;
8245 /* Communicate the number of cg's and atoms to receive */
8246 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8247 ind->nsend, nzone+2,
8248 ind->nrecv, nzone+2);
8250 /* The rvec buffer is also required for atom buffers of size nsend
8251 * in dd_move_x and dd_move_f.
8253 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8257 /* We can receive in place if only the last zone is not empty */
8258 for (zone = 0; zone < nzone-1; zone++)
8260 if (ind->nrecv[zone] > 0)
8262 cd->bInPlace = FALSE;
8267 /* The int buffer is only required here for the cg indices */
8268 if (ind->nrecv[nzone] > comm->nalloc_int2)
8270 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8271 srenew(comm->buf_int2, comm->nalloc_int2);
8273 /* The rvec buffer is also required for atom buffers
8274 * of size nrecv in dd_move_x and dd_move_f.
8276 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8277 vec_rvec_check_alloc(&comm->vbuf2, i);
8281 /* Make space for the global cg indices */
8282 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8283 || dd->cg_nalloc == 0)
8285 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8286 srenew(index_gl, dd->cg_nalloc);
8287 srenew(cgindex, dd->cg_nalloc+1);
8289 /* Communicate the global cg indices */
8292 recv_i = index_gl + pos_cg;
8296 recv_i = comm->buf_int2;
8298 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8299 comm->buf_int, nsend,
8300 recv_i, ind->nrecv[nzone]);
8302 /* Make space for cg_cm */
8303 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8304 if (fr->cutoff_scheme == ecutsGROUP)
8312 /* Communicate cg_cm */
8315 recv_vr = cg_cm + pos_cg;
8319 recv_vr = comm->vbuf2.v;
8321 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8322 comm->vbuf.v, nsend,
8323 recv_vr, ind->nrecv[nzone]);
8325 /* Make the charge group index */
8328 zone = (p == 0 ? 0 : nzone - 1);
8329 while (zone < nzone)
8331 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8333 cg_gl = index_gl[pos_cg];
8334 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8335 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8336 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8339 /* Update the charge group presence,
8340 * so we can use it in the next pass of the loop.
8342 comm->bLocalCG[cg_gl] = TRUE;
8348 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8351 zone_cg_range[nzone+zone] = pos_cg;
8356 /* This part of the code is never executed with bBondComm. */
8357 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8358 index_gl, recv_i, cg_cm, recv_vr,
8359 cgindex, fr->cginfo_mb, fr->cginfo);
8360 pos_cg += ind->nrecv[nzone];
8362 nat_tot += ind->nrecv[nzone+1];
8366 /* Store the atom block for easy copying of communication buffers */
8367 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8371 dd->index_gl = index_gl;
8372 dd->cgindex = cgindex;
8374 dd->ncg_tot = zone_cg_range[zones->n];
8375 dd->nat_tot = nat_tot;
8376 comm->nat[ddnatHOME] = dd->nat_home;
8377 for (i = ddnatZONE; i < ddnatNR; i++)
8379 comm->nat[i] = dd->nat_tot;
8384 /* We don't need to update cginfo, since that was alrady done above.
8385 * So we pass NULL for the forcerec.
8387 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8388 NULL, comm->bLocalCG);
8393 fprintf(debug, "Finished setting up DD communication, zones:");
8394 for (c = 0; c < zones->n; c++)
8396 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8398 fprintf(debug, "\n");
8402 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8406 for (c = 0; c < zones->nizone; c++)
8408 zones->izone[c].cg1 = zones->cg_range[c+1];
8409 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8410 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8414 static void set_zones_size(gmx_domdec_t *dd,
8415 matrix box, const gmx_ddbox_t *ddbox,
8416 int zone_start, int zone_end)
8418 gmx_domdec_comm_t *comm;
8419 gmx_domdec_zones_t *zones;
8421 int z, zi, zj0, zj1, d, dim;
8424 real size_j, add_tric;
8429 zones = &comm->zones;
8431 /* Do we need to determine extra distances for multi-body bondeds? */
8432 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8434 for (z = zone_start; z < zone_end; z++)
8436 /* Copy cell limits to zone limits.
8437 * Valid for non-DD dims and non-shifted dims.
8439 copy_rvec(comm->cell_x0, zones->size[z].x0);
8440 copy_rvec(comm->cell_x1, zones->size[z].x1);
8443 for (d = 0; d < dd->ndim; d++)
8447 for (z = 0; z < zones->n; z++)
8449 /* With a staggered grid we have different sizes
8450 * for non-shifted dimensions.
8452 if (dd->bGridJump && zones->shift[z][dim] == 0)
8456 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8457 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8461 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8462 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8468 rcmbs = comm->cutoff_mbody;
8469 if (ddbox->tric_dir[dim])
8471 rcs /= ddbox->skew_fac[dim];
8472 rcmbs /= ddbox->skew_fac[dim];
8475 /* Set the lower limit for the shifted zone dimensions */
8476 for (z = zone_start; z < zone_end; z++)
8478 if (zones->shift[z][dim] > 0)
8481 if (!dd->bGridJump || d == 0)
8483 zones->size[z].x0[dim] = comm->cell_x1[dim];
8484 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8488 /* Here we take the lower limit of the zone from
8489 * the lowest domain of the zone below.
8493 zones->size[z].x0[dim] =
8494 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8500 zones->size[z].x0[dim] =
8501 zones->size[zone_perm[2][z-4]].x0[dim];
8505 zones->size[z].x0[dim] =
8506 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8509 /* A temporary limit, is updated below */
8510 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8514 for (zi = 0; zi < zones->nizone; zi++)
8516 if (zones->shift[zi][dim] == 0)
8518 /* This takes the whole zone into account.
8519 * With multiple pulses this will lead
8520 * to a larger zone then strictly necessary.
8522 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8523 zones->size[zi].x1[dim]+rcmbs);
8531 /* Loop over the i-zones to set the upper limit of each
8534 for (zi = 0; zi < zones->nizone; zi++)
8536 if (zones->shift[zi][dim] == 0)
8538 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8540 if (zones->shift[z][dim] > 0)
8542 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8543 zones->size[zi].x1[dim]+rcs);
8550 for (z = zone_start; z < zone_end; z++)
8552 /* Initialization only required to keep the compiler happy */
8553 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8556 /* To determine the bounding box for a zone we need to find
8557 * the extreme corners of 4, 2 or 1 corners.
8559 nc = 1 << (ddbox->npbcdim - 1);
8561 for (c = 0; c < nc; c++)
8563 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8567 corner[YY] = zones->size[z].x0[YY];
8571 corner[YY] = zones->size[z].x1[YY];
8575 corner[ZZ] = zones->size[z].x0[ZZ];
8579 corner[ZZ] = zones->size[z].x1[ZZ];
8581 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8583 /* With 1D domain decomposition the cg's are not in
8584 * the triclinic box, but triclinic x-y and rectangular y-z.
8585 * Shift y back, so it will later end up at 0.
8587 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8589 /* Apply the triclinic couplings */
8590 for (i = YY; i < ddbox->npbcdim; i++)
8592 for (j = XX; j < i; j++)
8594 corner[j] += corner[i]*box[i][j]/box[i][i];
8599 copy_rvec(corner, corner_min);
8600 copy_rvec(corner, corner_max);
8604 for (i = 0; i < DIM; i++)
8606 corner_min[i] = min(corner_min[i], corner[i]);
8607 corner_max[i] = max(corner_max[i], corner[i]);
8611 /* Copy the extreme cornes without offset along x */
8612 for (i = 0; i < DIM; i++)
8614 zones->size[z].bb_x0[i] = corner_min[i];
8615 zones->size[z].bb_x1[i] = corner_max[i];
8617 /* Add the offset along x */
8618 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8619 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8622 if (zone_start == 0)
8625 for (dim = 0; dim < DIM; dim++)
8627 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8629 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8634 for (z = zone_start; z < zone_end; z++)
8636 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8638 zones->size[z].x0[XX], zones->size[z].x1[XX],
8639 zones->size[z].x0[YY], zones->size[z].x1[YY],
8640 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8641 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8643 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8644 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8645 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8650 static int comp_cgsort(const void *a, const void *b)
8654 gmx_cgsort_t *cga, *cgb;
8655 cga = (gmx_cgsort_t *)a;
8656 cgb = (gmx_cgsort_t *)b;
8658 comp = cga->nsc - cgb->nsc;
8661 comp = cga->ind_gl - cgb->ind_gl;
8667 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8672 /* Order the data */
8673 for (i = 0; i < n; i++)
8675 buf[i] = a[sort[i].ind];
8678 /* Copy back to the original array */
8679 for (i = 0; i < n; i++)
8685 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8690 /* Order the data */
8691 for (i = 0; i < n; i++)
8693 copy_rvec(v[sort[i].ind], buf[i]);
8696 /* Copy back to the original array */
8697 for (i = 0; i < n; i++)
8699 copy_rvec(buf[i], v[i]);
8703 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8706 int a, atot, cg, cg0, cg1, i;
8708 if (cgindex == NULL)
8710 /* Avoid the useless loop of the atoms within a cg */
8711 order_vec_cg(ncg, sort, v, buf);
8716 /* Order the data */
8718 for (cg = 0; cg < ncg; cg++)
8720 cg0 = cgindex[sort[cg].ind];
8721 cg1 = cgindex[sort[cg].ind+1];
8722 for (i = cg0; i < cg1; i++)
8724 copy_rvec(v[i], buf[a]);
8730 /* Copy back to the original array */
8731 for (a = 0; a < atot; a++)
8733 copy_rvec(buf[a], v[a]);
8737 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8738 int nsort_new, gmx_cgsort_t *sort_new,
8739 gmx_cgsort_t *sort1)
8743 /* The new indices are not very ordered, so we qsort them */
8744 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8746 /* sort2 is already ordered, so now we can merge the two arrays */
8750 while (i2 < nsort2 || i_new < nsort_new)
8754 sort1[i1++] = sort_new[i_new++];
8756 else if (i_new == nsort_new)
8758 sort1[i1++] = sort2[i2++];
8760 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8761 (sort2[i2].nsc == sort_new[i_new].nsc &&
8762 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8764 sort1[i1++] = sort2[i2++];
8768 sort1[i1++] = sort_new[i_new++];
8773 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8775 gmx_domdec_sort_t *sort;
8776 gmx_cgsort_t *cgsort, *sort_i;
8777 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8778 int sort_last, sort_skip;
8780 sort = dd->comm->sort;
8782 a = fr->ns.grid->cell_index;
8784 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8786 if (ncg_home_old >= 0)
8788 /* The charge groups that remained in the same ns grid cell
8789 * are completely ordered. So we can sort efficiently by sorting
8790 * the charge groups that did move into the stationary list.
8795 for (i = 0; i < dd->ncg_home; i++)
8797 /* Check if this cg did not move to another node */
8800 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8802 /* This cg is new on this node or moved ns grid cell */
8803 if (nsort_new >= sort->sort_new_nalloc)
8805 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8806 srenew(sort->sort_new, sort->sort_new_nalloc);
8808 sort_i = &(sort->sort_new[nsort_new++]);
8812 /* This cg did not move */
8813 sort_i = &(sort->sort2[nsort2++]);
8815 /* Sort on the ns grid cell indices
8816 * and the global topology index.
8817 * index_gl is irrelevant with cell ns,
8818 * but we set it here anyhow to avoid a conditional.
8821 sort_i->ind_gl = dd->index_gl[i];
8828 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8831 /* Sort efficiently */
8832 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8837 cgsort = sort->sort;
8839 for (i = 0; i < dd->ncg_home; i++)
8841 /* Sort on the ns grid cell indices
8842 * and the global topology index
8844 cgsort[i].nsc = a[i];
8845 cgsort[i].ind_gl = dd->index_gl[i];
8847 if (cgsort[i].nsc < moved)
8854 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8856 /* Determine the order of the charge groups using qsort */
8857 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8863 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8866 int ncg_new, i, *a, na;
8868 sort = dd->comm->sort->sort;
8870 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8873 for (i = 0; i < na; i++)
8877 sort[ncg_new].ind = a[i];
8885 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8886 rvec *cgcm, t_forcerec *fr, t_state *state,
8889 gmx_domdec_sort_t *sort;
8890 gmx_cgsort_t *cgsort, *sort_i;
8892 int ncg_new, i, *ibuf, cgsize;
8895 sort = dd->comm->sort;
8897 if (dd->ncg_home > sort->sort_nalloc)
8899 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8900 srenew(sort->sort, sort->sort_nalloc);
8901 srenew(sort->sort2, sort->sort_nalloc);
8903 cgsort = sort->sort;
8905 switch (fr->cutoff_scheme)
8908 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8911 ncg_new = dd_sort_order_nbnxn(dd, fr);
8914 gmx_incons("unimplemented");
8918 /* We alloc with the old size, since cgindex is still old */
8919 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8920 vbuf = dd->comm->vbuf.v;
8924 cgindex = dd->cgindex;
8931 /* Remove the charge groups which are no longer at home here */
8932 dd->ncg_home = ncg_new;
8935 fprintf(debug, "Set the new home charge group count to %d\n",
8939 /* Reorder the state */
8940 for (i = 0; i < estNR; i++)
8942 if (EST_DISTR(i) && (state->flags & (1<<i)))
8947 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8950 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8953 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8956 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8960 case estDISRE_INITF:
8961 case estDISRE_RM3TAV:
8962 case estORIRE_INITF:
8964 /* No ordering required */
8967 gmx_incons("Unknown state entry encountered in dd_sort_state");
8972 if (fr->cutoff_scheme == ecutsGROUP)
8975 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8978 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8980 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8981 srenew(sort->ibuf, sort->ibuf_nalloc);
8984 /* Reorder the global cg index */
8985 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
8986 /* Reorder the cginfo */
8987 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
8988 /* Rebuild the local cg index */
8992 for (i = 0; i < dd->ncg_home; i++)
8994 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8995 ibuf[i+1] = ibuf[i] + cgsize;
8997 for (i = 0; i < dd->ncg_home+1; i++)
8999 dd->cgindex[i] = ibuf[i];
9004 for (i = 0; i < dd->ncg_home+1; i++)
9009 /* Set the home atom number */
9010 dd->nat_home = dd->cgindex[dd->ncg_home];
9012 if (fr->cutoff_scheme == ecutsVERLET)
9014 /* The atoms are now exactly in grid order, update the grid order */
9015 nbnxn_set_atomorder(fr->nbv->nbs);
9019 /* Copy the sorted ns cell indices back to the ns grid struct */
9020 for (i = 0; i < dd->ncg_home; i++)
9022 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9024 fr->ns.grid->nr = dd->ncg_home;
9028 static void add_dd_statistics(gmx_domdec_t *dd)
9030 gmx_domdec_comm_t *comm;
9035 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9037 comm->sum_nat[ddnat-ddnatZONE] +=
9038 comm->nat[ddnat] - comm->nat[ddnat-1];
9043 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9045 gmx_domdec_comm_t *comm;
9050 /* Reset all the statistics and counters for total run counting */
9051 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9053 comm->sum_nat[ddnat-ddnatZONE] = 0;
9057 comm->load_step = 0;
9060 clear_ivec(comm->load_lim);
9065 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9067 gmx_domdec_comm_t *comm;
9071 comm = cr->dd->comm;
9073 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9080 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");
9082 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9084 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9089 " av. #atoms communicated per step for force: %d x %.1f\n",
9093 if (cr->dd->vsite_comm)
9096 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9097 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9102 if (cr->dd->constraint_comm)
9105 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9106 1 + ir->nLincsIter, av);
9110 gmx_incons(" Unknown type for DD statistics");
9113 fprintf(fplog, "\n");
9115 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9117 print_dd_load_av(fplog, cr->dd);
9121 void dd_partition_system(FILE *fplog,
9122 gmx_large_int_t step,
9124 gmx_bool bMasterState,
9126 t_state *state_global,
9127 gmx_mtop_t *top_global,
9129 t_state *state_local,
9132 gmx_localtop_t *top_local,
9135 gmx_shellfc_t shellfc,
9136 gmx_constr_t constr,
9138 gmx_wallcycle_t wcycle,
9142 gmx_domdec_comm_t *comm;
9143 gmx_ddbox_t ddbox = {0};
9145 gmx_large_int_t step_pcoupl;
9146 rvec cell_ns_x0, cell_ns_x1;
9147 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9148 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9149 gmx_bool bRedist, bSortCG, bResortAll;
9150 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9157 bBoxChanged = (bMasterState || DEFORM(*ir));
9158 if (ir->epc != epcNO)
9160 /* With nstpcouple > 1 pressure coupling happens.
9161 * one step after calculating the pressure.
9162 * Box scaling happens at the end of the MD step,
9163 * after the DD partitioning.
9164 * We therefore have to do DLB in the first partitioning
9165 * after an MD step where P-coupling occured.
9166 * We need to determine the last step in which p-coupling occurred.
9167 * MRS -- need to validate this for vv?
9172 step_pcoupl = step - 1;
9176 step_pcoupl = ((step - 1)/n)*n + 1;
9178 if (step_pcoupl >= comm->partition_step)
9184 bNStGlobalComm = (step % nstglobalcomm == 0);
9186 if (!comm->bDynLoadBal)
9192 /* Should we do dynamic load balacing this step?
9193 * Since it requires (possibly expensive) global communication,
9194 * we might want to do DLB less frequently.
9196 if (bBoxChanged || ir->epc != epcNO)
9198 bDoDLB = bBoxChanged;
9202 bDoDLB = bNStGlobalComm;
9206 /* Check if we have recorded loads on the nodes */
9207 if (comm->bRecordLoad && dd_load_count(comm))
9209 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9211 /* Check if we should use DLB at the second partitioning
9212 * and every 100 partitionings,
9213 * so the extra communication cost is negligible.
9215 n = max(100, nstglobalcomm);
9216 bCheckDLB = (comm->n_load_collect == 0 ||
9217 comm->n_load_have % n == n-1);
9224 /* Print load every nstlog, first and last step to the log file */
9225 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9226 comm->n_load_collect == 0 ||
9228 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9230 /* Avoid extra communication due to verbose screen output
9231 * when nstglobalcomm is set.
9233 if (bDoDLB || bLogLoad || bCheckDLB ||
9234 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9236 get_load_distribution(dd, wcycle);
9241 dd_print_load(fplog, dd, step-1);
9245 dd_print_load_verbose(dd);
9248 comm->n_load_collect++;
9252 /* Since the timings are node dependent, the master decides */
9256 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9259 fprintf(debug, "step %s, imb loss %f\n",
9260 gmx_step_str(step, sbuf),
9261 dd_force_imb_perf_loss(dd));
9264 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9267 turn_on_dlb(fplog, cr, step);
9272 comm->n_load_have++;
9275 cgs_gl = &comm->cgs_gl;
9280 /* Clear the old state */
9281 clear_dd_indices(dd, 0, 0);
9283 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9284 TRUE, cgs_gl, state_global->x, &ddbox);
9286 get_cg_distribution(fplog, step, dd, cgs_gl,
9287 state_global->box, &ddbox, state_global->x);
9289 dd_distribute_state(dd, cgs_gl,
9290 state_global, state_local, f);
9292 dd_make_local_cgs(dd, &top_local->cgs);
9294 /* Ensure that we have space for the new distribution */
9295 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9297 if (fr->cutoff_scheme == ecutsGROUP)
9299 calc_cgcm(fplog, 0, dd->ncg_home,
9300 &top_local->cgs, state_local->x, fr->cg_cm);
9303 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9305 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9309 else if (state_local->ddp_count != dd->ddp_count)
9311 if (state_local->ddp_count > dd->ddp_count)
9313 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9316 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9318 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);
9321 /* Clear the old state */
9322 clear_dd_indices(dd, 0, 0);
9324 /* Build the new indices */
9325 rebuild_cgindex(dd, cgs_gl->index, state_local);
9326 make_dd_indices(dd, cgs_gl->index, 0);
9328 if (fr->cutoff_scheme == ecutsGROUP)
9330 /* Redetermine the cg COMs */
9331 calc_cgcm(fplog, 0, dd->ncg_home,
9332 &top_local->cgs, state_local->x, fr->cg_cm);
9335 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9337 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9339 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9340 TRUE, &top_local->cgs, state_local->x, &ddbox);
9342 bRedist = comm->bDynLoadBal;
9346 /* We have the full state, only redistribute the cgs */
9348 /* Clear the non-home indices */
9349 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9351 /* Avoid global communication for dim's without pbc and -gcom */
9352 if (!bNStGlobalComm)
9354 copy_rvec(comm->box0, ddbox.box0 );
9355 copy_rvec(comm->box_size, ddbox.box_size);
9357 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9358 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9363 /* For dim's without pbc and -gcom */
9364 copy_rvec(ddbox.box0, comm->box0 );
9365 copy_rvec(ddbox.box_size, comm->box_size);
9367 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9370 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9372 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9375 /* Check if we should sort the charge groups */
9376 if (comm->nstSortCG > 0)
9378 bSortCG = (bMasterState ||
9379 (bRedist && (step % comm->nstSortCG == 0)));
9386 ncg_home_old = dd->ncg_home;
9391 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9393 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9394 state_local, f, fr, mdatoms,
9395 !bSortCG, nrnb, &cg0, &ncg_moved);
9397 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9400 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9402 &comm->cell_x0, &comm->cell_x1,
9403 dd->ncg_home, fr->cg_cm,
9404 cell_ns_x0, cell_ns_x1, &grid_density);
9408 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9411 switch (fr->cutoff_scheme)
9414 copy_ivec(fr->ns.grid->n, ncells_old);
9415 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9416 state_local->box, cell_ns_x0, cell_ns_x1,
9417 fr->rlistlong, grid_density);
9420 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9423 gmx_incons("unimplemented");
9425 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9426 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9430 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9432 /* Sort the state on charge group position.
9433 * This enables exact restarts from this step.
9434 * It also improves performance by about 15% with larger numbers
9435 * of atoms per node.
9438 /* Fill the ns grid with the home cell,
9439 * so we can sort with the indices.
9441 set_zones_ncg_home(dd);
9443 switch (fr->cutoff_scheme)
9446 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9448 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9450 comm->zones.size[0].bb_x0,
9451 comm->zones.size[0].bb_x1,
9453 comm->zones.dens_zone0,
9456 ncg_moved, bRedist ? comm->moved : NULL,
9457 fr->nbv->grp[eintLocal].kernel_type,
9458 fr->nbv->grp[eintLocal].nbat);
9460 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9463 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9464 0, dd->ncg_home, fr->cg_cm);
9466 copy_ivec(fr->ns.grid->n, ncells_new);
9469 gmx_incons("unimplemented");
9472 bResortAll = bMasterState;
9474 /* Check if we can user the old order and ns grid cell indices
9475 * of the charge groups to sort the charge groups efficiently.
9477 if (ncells_new[XX] != ncells_old[XX] ||
9478 ncells_new[YY] != ncells_old[YY] ||
9479 ncells_new[ZZ] != ncells_old[ZZ])
9486 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9487 gmx_step_str(step, sbuf), dd->ncg_home);
9489 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9490 bResortAll ? -1 : ncg_home_old);
9491 /* Rebuild all the indices */
9493 ga2la_clear(dd->ga2la);
9495 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9498 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9500 /* Setup up the communication and communicate the coordinates */
9501 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9503 /* Set the indices */
9504 make_dd_indices(dd, cgs_gl->index, cg0);
9506 /* Set the charge group boundaries for neighbor searching */
9507 set_cg_boundaries(&comm->zones);
9509 if (fr->cutoff_scheme == ecutsVERLET)
9511 set_zones_size(dd, state_local->box, &ddbox,
9512 bSortCG ? 1 : 0, comm->zones.n);
9515 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9518 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9519 -1,state_local->x,state_local->box);
9522 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9524 /* Extract a local topology from the global topology */
9525 for (i = 0; i < dd->ndim; i++)
9527 np[dd->dim[i]] = comm->cd[i].np;
9529 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9530 comm->cellsize_min, np,
9532 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9533 vsite, top_global, top_local);
9535 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9537 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9539 /* Set up the special atom communication */
9540 n = comm->nat[ddnatZONE];
9541 for (i = ddnatZONE+1; i < ddnatNR; i++)
9546 if (vsite && vsite->n_intercg_vsite)
9548 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9552 if (dd->bInterCGcons || dd->bInterCGsettles)
9554 /* Only for inter-cg constraints we need special code */
9555 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9556 constr, ir->nProjOrder,
9557 top_local->idef.il);
9561 gmx_incons("Unknown special atom type setup");
9566 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9568 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9570 /* Make space for the extra coordinates for virtual site
9571 * or constraint communication.
9573 state_local->natoms = comm->nat[ddnatNR-1];
9574 if (state_local->natoms > state_local->nalloc)
9576 dd_realloc_state(state_local, f, state_local->natoms);
9579 if (fr->bF_NoVirSum)
9581 if (vsite && vsite->n_intercg_vsite)
9583 nat_f_novirsum = comm->nat[ddnatVSITE];
9587 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9589 nat_f_novirsum = dd->nat_tot;
9593 nat_f_novirsum = dd->nat_home;
9602 /* Set the number of atoms required for the force calculation.
9603 * Forces need to be constrained when using a twin-range setup
9604 * or with energy minimization. For simple simulations we could
9605 * avoid some allocation, zeroing and copying, but this is
9606 * probably not worth the complications ande checking.
9608 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9609 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9611 /* We make the all mdatoms up to nat_tot_con.
9612 * We could save some work by only setting invmass
9613 * between nat_tot and nat_tot_con.
9615 /* This call also sets the new number of home particles to dd->nat_home */
9616 atoms2md(top_global, ir,
9617 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9619 /* Now we have the charges we can sort the FE interactions */
9620 dd_sort_local_top(dd, mdatoms, top_local);
9624 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9625 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9630 /* Make the local shell stuff, currently no communication is done */
9631 make_local_shells(cr, mdatoms, shellfc);
9634 if (ir->implicit_solvent)
9636 make_local_gb(cr, fr->born, ir->gb_algorithm);
9639 init_bonded_thread_force_reduction(fr, &top_local->idef);
9641 if (!(cr->duty & DUTY_PME))
9643 /* Send the charges to our PME only node */
9644 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9645 mdatoms->chargeA, mdatoms->chargeB,
9646 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9651 set_constraints(constr, top_local, ir, mdatoms, cr);
9654 if (ir->ePull != epullNO)
9656 /* Update the local pull groups */
9657 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9662 /* Update the local rotation groups */
9663 dd_make_local_rotation_groups(dd, ir->rot);
9667 add_dd_statistics(dd);
9669 /* Make sure we only count the cycles for this DD partitioning */
9670 clear_dd_cycle_counts(dd);
9672 /* Because the order of the atoms might have changed since
9673 * the last vsite construction, we need to communicate the constructing
9674 * atom coordinates again (for spreading the forces this MD step).
9676 dd_move_x_vsites(dd, state_local->box, state_local->x);
9678 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9680 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9682 dd_move_x(dd, state_local->box, state_local->x);
9683 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9684 -1, state_local->x, state_local->box);
9687 /* Store the partitioning step */
9688 comm->partition_step = step;
9690 /* Increase the DD partitioning counter */
9692 /* The state currently matches this DD partitioning count, store it */
9693 state_local->ddp_count = dd->ddp_count;
9696 /* The DD master node knows the complete cg distribution,
9697 * store the count so we can possibly skip the cg info communication.
9699 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9702 if (comm->DD_debug > 0)
9704 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9705 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9706 "after partitioning");