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"
67 #define DDRANK(dd, rank) (rank)
68 #define DDMASTERRANK(dd) (dd->masterrank)
70 typedef struct gmx_domdec_master
72 /* The cell boundaries */
74 /* The global charge group division */
75 int *ncg; /* Number of home charge groups for each node */
76 int *index; /* Index of nnodes+1 into cg */
77 int *cg; /* Global charge group index */
78 int *nat; /* Number of home atoms for each node. */
79 int *ibuf; /* Buffer for communication */
80 rvec *vbuf; /* Buffer for state scattering and gathering */
81 } gmx_domdec_master_t;
85 /* The numbers of charge groups to send and receive for each cell
86 * that requires communication, the last entry contains the total
87 * number of atoms that needs to be communicated.
89 int nsend[DD_MAXIZONE+2];
90 int nrecv[DD_MAXIZONE+2];
91 /* The charge groups to send */
94 /* The atom range for non-in-place communication */
95 int cell2at0[DD_MAXIZONE];
96 int cell2at1[DD_MAXIZONE];
101 int np; /* Number of grid pulses in this dimension */
102 int np_dlb; /* For dlb, for use with edlbAUTO */
103 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
105 gmx_bool bInPlace; /* Can we communicate in place? */
106 } gmx_domdec_comm_dim_t;
110 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
111 real *cell_f; /* State var.: cell boundaries, box relative */
112 real *old_cell_f; /* Temp. var.: old cell size */
113 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
114 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
115 real *bound_min; /* Temp. var.: lower limit for cell boundary */
116 real *bound_max; /* Temp. var.: upper limit for cell boundary */
117 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
118 real *buf_ncd; /* Temp. var. */
121 #define DD_NLOAD_MAX 9
123 /* Here floats are accurate enough, since these variables
124 * only influence the load balancing, not the actual MD results.
151 gmx_cgsort_t *sort_new;
163 /* This enum determines the order of the coordinates.
164 * ddnatHOME and ddnatZONE should be first and second,
165 * the others can be ordered as wanted.
168 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
172 edlbAUTO, edlbNO, edlbYES, edlbNR
174 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
178 int dim; /* The dimension */
179 gmx_bool dim_match; /* Tells if DD and PME dims match */
180 int nslab; /* The number of PME slabs in this dimension */
181 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
182 int *pp_min; /* The minimum pp node location, size nslab */
183 int *pp_max; /* The maximum pp node location,size nslab */
184 int maxshift; /* The maximum shift for coordinate redistribution in PME */
189 real min0; /* The minimum bottom of this zone */
190 real max1; /* The maximum top of this zone */
191 real min1; /* The minimum top of this zone */
192 real mch0; /* The maximum bottom communicaton height for this zone */
193 real mch1; /* The maximum top communicaton height for this zone */
194 real p1_0; /* The bottom value of the first cell in this zone */
195 real p1_1; /* The top value of the first cell in this zone */
200 gmx_domdec_ind_t ind;
207 } dd_comm_setup_work_t;
209 typedef struct gmx_domdec_comm
211 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
212 * unless stated otherwise.
215 /* The number of decomposition dimensions for PME, 0: no PME */
217 /* The number of nodes doing PME (PP/PME or only PME) */
221 /* The communication setup including the PME only nodes */
222 gmx_bool bCartesianPP_PME;
225 int *pmenodes; /* size npmenodes */
226 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
227 * but with bCartesianPP_PME */
228 gmx_ddpme_t ddpme[2];
230 /* The DD particle-particle nodes only */
231 gmx_bool bCartesianPP;
232 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
234 /* The global charge groups */
237 /* Should we sort the cgs */
239 gmx_domdec_sort_t *sort;
241 /* Are there charge groups? */
244 /* Are there bonded and multi-body interactions between charge groups? */
245 gmx_bool bInterCGBondeds;
246 gmx_bool bInterCGMultiBody;
248 /* Data for the optional bonded interaction atom communication range */
255 /* Are we actually using DLB? */
256 gmx_bool bDynLoadBal;
258 /* Cell sizes for static load balancing, first index cartesian */
261 /* The width of the communicated boundaries */
264 /* The minimum cell size (including triclinic correction) */
266 /* For dlb, for use with edlbAUTO */
267 rvec cellsize_min_dlb;
268 /* The lower limit for the DD cell size with DLB */
270 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
271 gmx_bool bVacDLBNoLimit;
273 /* With PME load balancing we set limits on DLB */
274 gmx_bool bPMELoadBalDLBLimits;
275 /* DLB needs to take into account that we want to allow this maximum
276 * cut-off (for PME load balancing), this could limit cell boundaries.
278 real PMELoadBal_max_cutoff;
280 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
282 /* box0 and box_size are required with dim's without pbc and -gcom */
286 /* The cell boundaries */
290 /* The old location of the cell boundaries, to check cg displacements */
294 /* The communication setup and charge group boundaries for the zones */
295 gmx_domdec_zones_t zones;
297 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
298 * cell boundaries of neighboring cells for dynamic load balancing.
300 gmx_ddzone_t zone_d1[2];
301 gmx_ddzone_t zone_d2[2][2];
303 /* The coordinate/force communication setup and indices */
304 gmx_domdec_comm_dim_t cd[DIM];
305 /* The maximum number of cells to communicate with in one dimension */
308 /* Which cg distribution is stored on the master node */
309 int master_cg_ddp_count;
311 /* The number of cg's received from the direct neighbors */
312 int zone_ncg1[DD_MAXZONE];
314 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
317 /* Array for signalling if atoms have moved to another domain */
321 /* Communication buffer for general use */
325 /* Communication buffer for general use */
328 /* Temporary storage for thread parallel communication setup */
330 dd_comm_setup_work_t *dth;
332 /* Communication buffers only used with multiple grid pulses */
337 /* Communication buffers for local redistribution */
339 int cggl_flag_nalloc[DIM*2];
341 int cgcm_state_nalloc[DIM*2];
343 /* Cell sizes for dynamic load balancing */
344 gmx_domdec_root_t **root;
348 real cell_f_max0[DIM];
349 real cell_f_min1[DIM];
351 /* Stuff for load communication */
352 gmx_bool bRecordLoad;
353 gmx_domdec_load_t *load;
355 MPI_Comm *mpi_comm_load;
358 /* Maximum DLB scaling per load balancing step in percent */
362 float cycl[ddCyclNr];
363 int cycl_n[ddCyclNr];
364 float cycl_max[ddCyclNr];
365 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
369 /* Have often have did we have load measurements */
371 /* Have often have we collected the load measurements */
375 double sum_nat[ddnatNR-ddnatZONE];
385 /* The last partition step */
386 gmx_large_int_t partition_step;
394 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
397 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
398 #define DD_FLAG_NRCG 65535
399 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
400 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
402 /* Zone permutation required to obtain consecutive charge groups
403 * for neighbor searching.
405 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
407 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
408 * components see only j zones with that component 0.
411 /* The DD zone order */
412 static const ivec dd_zo[DD_MAXZONE] =
413 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
418 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
423 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
428 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
430 /* Factors used to avoid problems due to rounding issues */
431 #define DD_CELL_MARGIN 1.0001
432 #define DD_CELL_MARGIN2 1.00005
433 /* Factor to account for pressure scaling during nstlist steps */
434 #define DD_PRES_SCALE_MARGIN 1.02
436 /* Allowed performance loss before we DLB or warn */
437 #define DD_PERF_LOSS 0.05
439 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
441 /* Use separate MPI send and receive commands
442 * when nnodes <= GMX_DD_NNODES_SENDRECV.
443 * This saves memory (and some copying for small nnodes).
444 * For high parallelization scatter and gather calls are used.
446 #define GMX_DD_NNODES_SENDRECV 4
450 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
452 static void index2xyz(ivec nc,int ind,ivec xyz)
454 xyz[XX] = ind % nc[XX];
455 xyz[YY] = (ind / nc[XX]) % nc[YY];
456 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
460 /* This order is required to minimize the coordinate communication in PME
461 * which uses decomposition in the x direction.
463 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
465 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
467 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
468 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
469 xyz[ZZ] = ind % nc[ZZ];
472 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
477 ddindex = dd_index(dd->nc, c);
478 if (dd->comm->bCartesianPP_PME)
480 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
482 else if (dd->comm->bCartesianPP)
485 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
496 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
498 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
501 int ddglatnr(gmx_domdec_t *dd, int i)
511 if (i >= dd->comm->nat[ddnatNR-1])
513 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
515 atnr = dd->gatindex[i] + 1;
521 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
523 return &dd->comm->cgs_gl;
526 static void vec_rvec_init(vec_rvec_t *v)
532 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
536 v->nalloc = over_alloc_dd(n);
537 srenew(v->v, v->nalloc);
541 void dd_store_state(gmx_domdec_t *dd, t_state *state)
545 if (state->ddp_count != dd->ddp_count)
547 gmx_incons("The state does not the domain decomposition state");
550 state->ncg_gl = dd->ncg_home;
551 if (state->ncg_gl > state->cg_gl_nalloc)
553 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
554 srenew(state->cg_gl, state->cg_gl_nalloc);
556 for (i = 0; i < state->ncg_gl; i++)
558 state->cg_gl[i] = dd->index_gl[i];
561 state->ddp_count_cg_gl = dd->ddp_count;
564 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
566 return &dd->comm->zones;
569 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
570 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
572 gmx_domdec_zones_t *zones;
575 zones = &dd->comm->zones;
578 while (icg >= zones->izone[izone].cg1)
587 else if (izone < zones->nizone)
589 *jcg0 = zones->izone[izone].jcg0;
593 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
594 icg, izone, zones->nizone);
597 *jcg1 = zones->izone[izone].jcg1;
599 for (d = 0; d < dd->ndim; d++)
602 shift0[dim] = zones->izone[izone].shift0[dim];
603 shift1[dim] = zones->izone[izone].shift1[dim];
604 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
606 /* A conservative approach, this can be optimized */
613 int dd_natoms_vsite(gmx_domdec_t *dd)
615 return dd->comm->nat[ddnatVSITE];
618 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
620 *at_start = dd->comm->nat[ddnatCON-1];
621 *at_end = dd->comm->nat[ddnatCON];
624 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
626 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
627 int *index, *cgindex;
628 gmx_domdec_comm_t *comm;
629 gmx_domdec_comm_dim_t *cd;
630 gmx_domdec_ind_t *ind;
631 rvec shift = {0, 0, 0}, *buf, *rbuf;
632 gmx_bool bPBC, bScrew;
636 cgindex = dd->cgindex;
641 nat_tot = dd->nat_home;
642 for (d = 0; d < dd->ndim; d++)
644 bPBC = (dd->ci[dd->dim[d]] == 0);
645 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
648 copy_rvec(box[dd->dim[d]], shift);
651 for (p = 0; p < cd->np; p++)
658 for (i = 0; i < ind->nsend[nzone]; i++)
660 at0 = cgindex[index[i]];
661 at1 = cgindex[index[i]+1];
662 for (j = at0; j < at1; j++)
664 copy_rvec(x[j], buf[n]);
671 for (i = 0; i < ind->nsend[nzone]; i++)
673 at0 = cgindex[index[i]];
674 at1 = cgindex[index[i]+1];
675 for (j = at0; j < at1; j++)
677 /* We need to shift the coordinates */
678 rvec_add(x[j], shift, buf[n]);
685 for (i = 0; i < ind->nsend[nzone]; i++)
687 at0 = cgindex[index[i]];
688 at1 = cgindex[index[i]+1];
689 for (j = at0; j < at1; j++)
692 buf[n][XX] = x[j][XX] + shift[XX];
694 * This operation requires a special shift force
695 * treatment, which is performed in calc_vir.
697 buf[n][YY] = box[YY][YY] - x[j][YY];
698 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
710 rbuf = comm->vbuf2.v;
712 /* Send and receive the coordinates */
713 dd_sendrecv_rvec(dd, d, dddirBackward,
714 buf, ind->nsend[nzone+1],
715 rbuf, ind->nrecv[nzone+1]);
719 for (zone = 0; zone < nzone; zone++)
721 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
723 copy_rvec(rbuf[j], x[i]);
728 nat_tot += ind->nrecv[nzone+1];
734 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
736 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
737 int *index, *cgindex;
738 gmx_domdec_comm_t *comm;
739 gmx_domdec_comm_dim_t *cd;
740 gmx_domdec_ind_t *ind;
744 gmx_bool bPBC, bScrew;
748 cgindex = dd->cgindex;
753 nzone = comm->zones.n/2;
754 nat_tot = dd->nat_tot;
755 for (d = dd->ndim-1; d >= 0; d--)
757 bPBC = (dd->ci[dd->dim[d]] == 0);
758 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
759 if (fshift == NULL && !bScrew)
763 /* Determine which shift vector we need */
769 for (p = cd->np-1; p >= 0; p--)
772 nat_tot -= ind->nrecv[nzone+1];
779 sbuf = comm->vbuf2.v;
781 for (zone = 0; zone < nzone; zone++)
783 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
785 copy_rvec(f[i], sbuf[j]);
790 /* Communicate the forces */
791 dd_sendrecv_rvec(dd, d, dddirForward,
792 sbuf, ind->nrecv[nzone+1],
793 buf, ind->nsend[nzone+1]);
795 /* Add the received forces */
799 for (i = 0; i < ind->nsend[nzone]; i++)
801 at0 = cgindex[index[i]];
802 at1 = cgindex[index[i]+1];
803 for (j = at0; j < at1; j++)
805 rvec_inc(f[j], buf[n]);
812 for (i = 0; i < ind->nsend[nzone]; i++)
814 at0 = cgindex[index[i]];
815 at1 = cgindex[index[i]+1];
816 for (j = at0; j < at1; j++)
818 rvec_inc(f[j], buf[n]);
819 /* Add this force to the shift force */
820 rvec_inc(fshift[is], buf[n]);
827 for (i = 0; i < ind->nsend[nzone]; i++)
829 at0 = cgindex[index[i]];
830 at1 = cgindex[index[i]+1];
831 for (j = at0; j < at1; j++)
833 /* Rotate the force */
834 f[j][XX] += buf[n][XX];
835 f[j][YY] -= buf[n][YY];
836 f[j][ZZ] -= buf[n][ZZ];
839 /* Add this force to the shift force */
840 rvec_inc(fshift[is], buf[n]);
851 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
853 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
854 int *index, *cgindex;
855 gmx_domdec_comm_t *comm;
856 gmx_domdec_comm_dim_t *cd;
857 gmx_domdec_ind_t *ind;
862 cgindex = dd->cgindex;
864 buf = &comm->vbuf.v[0][0];
867 nat_tot = dd->nat_home;
868 for (d = 0; d < dd->ndim; d++)
871 for (p = 0; p < cd->np; p++)
876 for (i = 0; i < ind->nsend[nzone]; i++)
878 at0 = cgindex[index[i]];
879 at1 = cgindex[index[i]+1];
880 for (j = at0; j < at1; j++)
893 rbuf = &comm->vbuf2.v[0][0];
895 /* Send and receive the coordinates */
896 dd_sendrecv_real(dd, d, dddirBackward,
897 buf, ind->nsend[nzone+1],
898 rbuf, ind->nrecv[nzone+1]);
902 for (zone = 0; zone < nzone; zone++)
904 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
911 nat_tot += ind->nrecv[nzone+1];
917 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
919 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
920 int *index, *cgindex;
921 gmx_domdec_comm_t *comm;
922 gmx_domdec_comm_dim_t *cd;
923 gmx_domdec_ind_t *ind;
928 cgindex = dd->cgindex;
930 buf = &comm->vbuf.v[0][0];
933 nzone = comm->zones.n/2;
934 nat_tot = dd->nat_tot;
935 for (d = dd->ndim-1; d >= 0; d--)
938 for (p = cd->np-1; p >= 0; p--)
941 nat_tot -= ind->nrecv[nzone+1];
948 sbuf = &comm->vbuf2.v[0][0];
950 for (zone = 0; zone < nzone; zone++)
952 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
959 /* Communicate the forces */
960 dd_sendrecv_real(dd, d, dddirForward,
961 sbuf, ind->nrecv[nzone+1],
962 buf, ind->nsend[nzone+1]);
964 /* Add the received forces */
966 for (i = 0; i < ind->nsend[nzone]; i++)
968 at0 = cgindex[index[i]];
969 at1 = cgindex[index[i]+1];
970 for (j = at0; j < at1; j++)
981 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
983 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",
985 zone->min0, zone->max1,
986 zone->mch0, zone->mch0,
987 zone->p1_0, zone->p1_1);
991 #define DDZONECOMM_MAXZONE 5
992 #define DDZONECOMM_BUFSIZE 3
994 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
995 int ddimind, int direction,
996 gmx_ddzone_t *buf_s, int n_s,
997 gmx_ddzone_t *buf_r, int n_r)
999 #define ZBS DDZONECOMM_BUFSIZE
1000 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1001 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1004 for (i = 0; i < n_s; i++)
1006 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1007 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1008 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1009 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1010 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1011 vbuf_s[i*ZBS+1][2] = 0;
1012 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1013 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1014 vbuf_s[i*ZBS+2][2] = 0;
1017 dd_sendrecv_rvec(dd, ddimind, direction,
1021 for (i = 0; i < n_r; i++)
1023 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1024 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1025 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1026 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1027 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1028 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1029 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1035 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1036 rvec cell_ns_x0, rvec cell_ns_x1)
1038 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1040 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1041 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1042 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1043 rvec extr_s[2], extr_r[2];
1045 real dist_d, c = 0, det;
1046 gmx_domdec_comm_t *comm;
1047 gmx_bool bPBC, bUse;
1051 for (d = 1; d < dd->ndim; d++)
1054 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1055 zp->min0 = cell_ns_x0[dim];
1056 zp->max1 = cell_ns_x1[dim];
1057 zp->min1 = cell_ns_x1[dim];
1058 zp->mch0 = cell_ns_x0[dim];
1059 zp->mch1 = cell_ns_x1[dim];
1060 zp->p1_0 = cell_ns_x0[dim];
1061 zp->p1_1 = cell_ns_x1[dim];
1064 for (d = dd->ndim-2; d >= 0; d--)
1067 bPBC = (dim < ddbox->npbcdim);
1069 /* Use an rvec to store two reals */
1070 extr_s[d][0] = comm->cell_f0[d+1];
1071 extr_s[d][1] = comm->cell_f1[d+1];
1072 extr_s[d][2] = comm->cell_f1[d+1];
1075 /* Store the extremes in the backward sending buffer,
1076 * so the get updated separately from the forward communication.
1078 for (d1 = d; d1 < dd->ndim-1; d1++)
1080 /* We invert the order to be able to use the same loop for buf_e */
1081 buf_s[pos].min0 = extr_s[d1][1];
1082 buf_s[pos].max1 = extr_s[d1][0];
1083 buf_s[pos].min1 = extr_s[d1][2];
1084 buf_s[pos].mch0 = 0;
1085 buf_s[pos].mch1 = 0;
1086 /* Store the cell corner of the dimension we communicate along */
1087 buf_s[pos].p1_0 = comm->cell_x0[dim];
1088 buf_s[pos].p1_1 = 0;
1092 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1095 if (dd->ndim == 3 && d == 0)
1097 buf_s[pos] = comm->zone_d2[0][1];
1099 buf_s[pos] = comm->zone_d1[0];
1103 /* We only need to communicate the extremes
1104 * in the forward direction
1106 npulse = comm->cd[d].np;
1109 /* Take the minimum to avoid double communication */
1110 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1114 /* Without PBC we should really not communicate over
1115 * the boundaries, but implementing that complicates
1116 * the communication setup and therefore we simply
1117 * do all communication, but ignore some data.
1119 npulse_min = npulse;
1121 for (p = 0; p < npulse_min; p++)
1123 /* Communicate the extremes forward */
1124 bUse = (bPBC || dd->ci[dim] > 0);
1126 dd_sendrecv_rvec(dd, d, dddirForward,
1127 extr_s+d, dd->ndim-d-1,
1128 extr_r+d, dd->ndim-d-1);
1132 for (d1 = d; d1 < dd->ndim-1; d1++)
1134 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1135 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1136 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1142 for (p = 0; p < npulse; p++)
1144 /* Communicate all the zone information backward */
1145 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1147 dd_sendrecv_ddzone(dd, d, dddirBackward,
1154 for (d1 = d+1; d1 < dd->ndim; d1++)
1156 /* Determine the decrease of maximum required
1157 * communication height along d1 due to the distance along d,
1158 * this avoids a lot of useless atom communication.
1160 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1162 if (ddbox->tric_dir[dim])
1164 /* c is the off-diagonal coupling between the cell planes
1165 * along directions d and d1.
1167 c = ddbox->v[dim][dd->dim[d1]][dim];
1173 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1176 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1180 /* A negative value signals out of range */
1186 /* Accumulate the extremes over all pulses */
1187 for (i = 0; i < buf_size; i++)
1191 buf_e[i] = buf_r[i];
1197 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1198 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1199 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1202 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1210 if (bUse && dh[d1] >= 0)
1212 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1213 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1216 /* Copy the received buffer to the send buffer,
1217 * to pass the data through with the next pulse.
1219 buf_s[i] = buf_r[i];
1221 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1222 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1224 /* Store the extremes */
1227 for (d1 = d; d1 < dd->ndim-1; d1++)
1229 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1230 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1231 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1235 if (d == 1 || (d == 0 && dd->ndim == 3))
1237 for (i = d; i < 2; i++)
1239 comm->zone_d2[1-d][i] = buf_e[pos];
1245 comm->zone_d1[1] = buf_e[pos];
1255 for (i = 0; i < 2; i++)
1259 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1261 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1262 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1268 for (i = 0; i < 2; i++)
1270 for (j = 0; j < 2; j++)
1274 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1276 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1277 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1281 for (d = 1; d < dd->ndim; d++)
1283 comm->cell_f_max0[d] = extr_s[d-1][0];
1284 comm->cell_f_min1[d] = extr_s[d-1][1];
1287 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1288 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1293 static void dd_collect_cg(gmx_domdec_t *dd,
1294 t_state *state_local)
1296 gmx_domdec_master_t *ma = NULL;
1297 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1300 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1302 /* The master has the correct distribution */
1306 if (state_local->ddp_count == dd->ddp_count)
1308 ncg_home = dd->ncg_home;
1310 nat_home = dd->nat_home;
1312 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1314 cgs_gl = &dd->comm->cgs_gl;
1316 ncg_home = state_local->ncg_gl;
1317 cg = state_local->cg_gl;
1319 for (i = 0; i < ncg_home; i++)
1321 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1326 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1329 buf2[0] = dd->ncg_home;
1330 buf2[1] = dd->nat_home;
1340 /* Collect the charge group and atom counts on the master */
1341 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1346 for (i = 0; i < dd->nnodes; i++)
1348 ma->ncg[i] = ma->ibuf[2*i];
1349 ma->nat[i] = ma->ibuf[2*i+1];
1350 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1353 /* Make byte counts and indices */
1354 for (i = 0; i < dd->nnodes; i++)
1356 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1357 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1361 fprintf(debug, "Initial charge group distribution: ");
1362 for (i = 0; i < dd->nnodes; i++)
1364 fprintf(debug, " %d", ma->ncg[i]);
1366 fprintf(debug, "\n");
1370 /* Collect the charge group indices on the master */
1372 dd->ncg_home*sizeof(int), dd->index_gl,
1373 DDMASTER(dd) ? ma->ibuf : NULL,
1374 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1375 DDMASTER(dd) ? ma->cg : NULL);
1377 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1380 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1383 gmx_domdec_master_t *ma;
1384 int n, i, c, a, nalloc = 0;
1393 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1394 dd->rank, dd->mpi_comm_all);
1399 /* Copy the master coordinates to the global array */
1400 cgs_gl = &dd->comm->cgs_gl;
1402 n = DDMASTERRANK(dd);
1404 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1406 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1408 copy_rvec(lv[a++], v[c]);
1412 for (n = 0; n < dd->nnodes; n++)
1416 if (ma->nat[n] > nalloc)
1418 nalloc = over_alloc_dd(ma->nat[n]);
1419 srenew(buf, nalloc);
1422 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1423 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1426 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1428 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1430 copy_rvec(buf[a++], v[c]);
1439 static void get_commbuffer_counts(gmx_domdec_t *dd,
1440 int **counts, int **disps)
1442 gmx_domdec_master_t *ma;
1447 /* Make the rvec count and displacment arrays */
1449 *disps = ma->ibuf + dd->nnodes;
1450 for (n = 0; n < dd->nnodes; n++)
1452 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1453 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1457 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1460 gmx_domdec_master_t *ma;
1461 int *rcounts = NULL, *disps = NULL;
1470 get_commbuffer_counts(dd, &rcounts, &disps);
1475 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1479 cgs_gl = &dd->comm->cgs_gl;
1482 for (n = 0; n < dd->nnodes; n++)
1484 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1486 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1488 copy_rvec(buf[a++], v[c]);
1495 void dd_collect_vec(gmx_domdec_t *dd,
1496 t_state *state_local, rvec *lv, rvec *v)
1498 gmx_domdec_master_t *ma;
1499 int n, i, c, a, nalloc = 0;
1502 dd_collect_cg(dd, state_local);
1504 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1506 dd_collect_vec_sendrecv(dd, lv, v);
1510 dd_collect_vec_gatherv(dd, lv, v);
1515 void dd_collect_state(gmx_domdec_t *dd,
1516 t_state *state_local, t_state *state)
1520 nh = state->nhchainlength;
1524 for (i = 0; i < efptNR; i++)
1526 state->lambda[i] = state_local->lambda[i];
1528 state->fep_state = state_local->fep_state;
1529 state->veta = state_local->veta;
1530 state->vol0 = state_local->vol0;
1531 copy_mat(state_local->box, state->box);
1532 copy_mat(state_local->boxv, state->boxv);
1533 copy_mat(state_local->svir_prev, state->svir_prev);
1534 copy_mat(state_local->fvir_prev, state->fvir_prev);
1535 copy_mat(state_local->pres_prev, state->pres_prev);
1538 for (i = 0; i < state_local->ngtc; i++)
1540 for (j = 0; j < nh; j++)
1542 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1543 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1545 state->therm_integral[i] = state_local->therm_integral[i];
1547 for (i = 0; i < state_local->nnhpres; i++)
1549 for (j = 0; j < nh; j++)
1551 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1552 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1556 for (est = 0; est < estNR; est++)
1558 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1563 dd_collect_vec(dd, state_local, state_local->x, state->x);
1566 dd_collect_vec(dd, state_local, state_local->v, state->v);
1569 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1572 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1575 if (state->nrngi == 1)
1579 for (i = 0; i < state_local->nrng; i++)
1581 state->ld_rng[i] = state_local->ld_rng[i];
1587 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1588 state_local->ld_rng, state->ld_rng);
1592 if (state->nrngi == 1)
1596 state->ld_rngi[0] = state_local->ld_rngi[0];
1601 dd_gather(dd, sizeof(state->ld_rngi[0]),
1602 state_local->ld_rngi, state->ld_rngi);
1605 case estDISRE_INITF:
1606 case estDISRE_RM3TAV:
1607 case estORIRE_INITF:
1611 gmx_incons("Unknown state entry encountered in dd_collect_state");
1617 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1623 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1626 state->nalloc = over_alloc_dd(nalloc);
1628 for (est = 0; est < estNR; est++)
1630 if (EST_DISTR(est) && (state->flags & (1<<est)))
1635 srenew(state->x, state->nalloc);
1638 srenew(state->v, state->nalloc);
1641 srenew(state->sd_X, state->nalloc);
1644 srenew(state->cg_p, state->nalloc);
1648 case estDISRE_INITF:
1649 case estDISRE_RM3TAV:
1650 case estORIRE_INITF:
1652 /* No reallocation required */
1655 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1662 srenew(*f, state->nalloc);
1666 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1669 if (nalloc > fr->cg_nalloc)
1673 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1675 fr->cg_nalloc = over_alloc_dd(nalloc);
1676 srenew(fr->cginfo, fr->cg_nalloc);
1677 if (fr->cutoff_scheme == ecutsGROUP)
1679 srenew(fr->cg_cm, fr->cg_nalloc);
1682 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1684 /* We don't use charge groups, we use x in state to set up
1685 * the atom communication.
1687 dd_realloc_state(state, f, nalloc);
1691 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1694 gmx_domdec_master_t *ma;
1695 int n, i, c, a, nalloc = 0;
1702 for (n = 0; n < dd->nnodes; n++)
1706 if (ma->nat[n] > nalloc)
1708 nalloc = over_alloc_dd(ma->nat[n]);
1709 srenew(buf, nalloc);
1711 /* Use lv as a temporary buffer */
1713 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1715 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1717 copy_rvec(v[c], buf[a++]);
1720 if (a != ma->nat[n])
1722 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1727 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1728 DDRANK(dd, n), n, dd->mpi_comm_all);
1733 n = DDMASTERRANK(dd);
1735 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1737 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1739 copy_rvec(v[c], lv[a++]);
1746 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1747 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1752 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1755 gmx_domdec_master_t *ma;
1756 int *scounts = NULL, *disps = NULL;
1757 int n, i, c, a, nalloc = 0;
1764 get_commbuffer_counts(dd, &scounts, &disps);
1768 for (n = 0; n < dd->nnodes; n++)
1770 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1772 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1774 copy_rvec(v[c], buf[a++]);
1780 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1783 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1785 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1787 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1791 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1795 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1796 t_state *state, t_state *state_local,
1801 nh = state->nhchainlength;
1805 for (i = 0; i < efptNR; i++)
1807 state_local->lambda[i] = state->lambda[i];
1809 state_local->fep_state = state->fep_state;
1810 state_local->veta = state->veta;
1811 state_local->vol0 = state->vol0;
1812 copy_mat(state->box, state_local->box);
1813 copy_mat(state->box_rel, state_local->box_rel);
1814 copy_mat(state->boxv, state_local->boxv);
1815 copy_mat(state->svir_prev, state_local->svir_prev);
1816 copy_mat(state->fvir_prev, state_local->fvir_prev);
1817 for (i = 0; i < state_local->ngtc; i++)
1819 for (j = 0; j < nh; j++)
1821 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1822 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1824 state_local->therm_integral[i] = state->therm_integral[i];
1826 for (i = 0; i < state_local->nnhpres; i++)
1828 for (j = 0; j < nh; j++)
1830 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1831 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1835 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1836 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1837 dd_bcast(dd, sizeof(real), &state_local->veta);
1838 dd_bcast(dd, sizeof(real), &state_local->vol0);
1839 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1840 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1841 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1842 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1843 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1844 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1845 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1846 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1847 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1848 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1850 if (dd->nat_home > state_local->nalloc)
1852 dd_realloc_state(state_local, f, dd->nat_home);
1854 for (i = 0; i < estNR; i++)
1856 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1861 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1864 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1867 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1870 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1873 if (state->nrngi == 1)
1876 state_local->nrng*sizeof(state_local->ld_rng[0]),
1877 state->ld_rng, state_local->ld_rng);
1882 state_local->nrng*sizeof(state_local->ld_rng[0]),
1883 state->ld_rng, state_local->ld_rng);
1887 if (state->nrngi == 1)
1889 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1890 state->ld_rngi, state_local->ld_rngi);
1894 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1895 state->ld_rngi, state_local->ld_rngi);
1898 case estDISRE_INITF:
1899 case estDISRE_RM3TAV:
1900 case estORIRE_INITF:
1902 /* Not implemented yet */
1905 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1911 static char dim2char(int dim)
1917 case XX: c = 'X'; break;
1918 case YY: c = 'Y'; break;
1919 case ZZ: c = 'Z'; break;
1920 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1926 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1927 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1929 rvec grid_s[2], *grid_r = NULL, cx, r;
1930 char fname[STRLEN], format[STRLEN], buf[22];
1932 int a, i, d, z, y, x;
1936 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1937 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1941 snew(grid_r, 2*dd->nnodes);
1944 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1948 for (d = 0; d < DIM; d++)
1950 for (i = 0; i < DIM; i++)
1958 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1960 tric[d][i] = box[i][d]/box[i][i];
1969 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1970 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
1971 out = gmx_fio_fopen(fname, "w");
1972 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1974 for (i = 0; i < dd->nnodes; i++)
1976 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1977 for (d = 0; d < DIM; d++)
1979 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1981 for (z = 0; z < 2; z++)
1983 for (y = 0; y < 2; y++)
1985 for (x = 0; x < 2; x++)
1987 cx[XX] = grid_r[i*2+x][XX];
1988 cx[YY] = grid_r[i*2+y][YY];
1989 cx[ZZ] = grid_r[i*2+z][ZZ];
1991 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
1992 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
1996 for (d = 0; d < DIM; d++)
1998 for (x = 0; x < 4; x++)
2002 case 0: y = 1 + i*8 + 2*x; break;
2003 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2004 case 2: y = 1 + i*8 + x; break;
2006 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2010 gmx_fio_fclose(out);
2015 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2016 gmx_mtop_t *mtop, t_commrec *cr,
2017 int natoms, rvec x[], matrix box)
2019 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2021 int i, ii, resnr, c;
2022 char *atomname, *resname;
2029 natoms = dd->comm->nat[ddnatVSITE];
2032 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2034 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
2035 sprintf(format4, "%s%s\n", get_pdbformat4(), "%6.2f%6.2f");
2037 out = gmx_fio_fopen(fname, "w");
2039 fprintf(out, "TITLE %s\n", title);
2040 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2041 for (i = 0; i < natoms; i++)
2043 ii = dd->gatindex[i];
2044 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2045 if (i < dd->comm->nat[ddnatZONE])
2048 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2054 else if (i < dd->comm->nat[ddnatVSITE])
2056 b = dd->comm->zones.n;
2060 b = dd->comm->zones.n + 1;
2062 fprintf(out, strlen(atomname) < 4 ? format : format4,
2063 "ATOM", (ii+1)%100000,
2064 atomname, resname, ' ', resnr%10000, ' ',
2065 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2067 fprintf(out, "TER\n");
2069 gmx_fio_fclose(out);
2072 real dd_cutoff_mbody(gmx_domdec_t *dd)
2074 gmx_domdec_comm_t *comm;
2081 if (comm->bInterCGBondeds)
2083 if (comm->cutoff_mbody > 0)
2085 r = comm->cutoff_mbody;
2089 /* cutoff_mbody=0 means we do not have DLB */
2090 r = comm->cellsize_min[dd->dim[0]];
2091 for (di = 1; di < dd->ndim; di++)
2093 r = min(r, comm->cellsize_min[dd->dim[di]]);
2095 if (comm->bBondComm)
2097 r = max(r, comm->cutoff_mbody);
2101 r = min(r, comm->cutoff);
2109 real dd_cutoff_twobody(gmx_domdec_t *dd)
2113 r_mb = dd_cutoff_mbody(dd);
2115 return max(dd->comm->cutoff, r_mb);
2119 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2123 nc = dd->nc[dd->comm->cartpmedim];
2124 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2125 copy_ivec(coord, coord_pme);
2126 coord_pme[dd->comm->cartpmedim] =
2127 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2130 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2132 /* Here we assign a PME node to communicate with this DD node
2133 * by assuming that the major index of both is x.
2134 * We add cr->npmenodes/2 to obtain an even distribution.
2136 return (ddindex*npme + npme/2)/ndd;
2139 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2141 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2144 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2146 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2149 static int *dd_pmenodes(t_commrec *cr)
2154 snew(pmenodes, cr->npmenodes);
2156 for (i = 0; i < cr->dd->nnodes; i++)
2158 p0 = cr_ddindex2pmeindex(cr, i);
2159 p1 = cr_ddindex2pmeindex(cr, i+1);
2160 if (i+1 == cr->dd->nnodes || p1 > p0)
2164 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2166 pmenodes[n] = i + 1 + n;
2174 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2177 ivec coords, coords_pme, nc;
2182 if (dd->comm->bCartesian) {
2183 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2184 dd_coords2pmecoords(dd,coords,coords_pme);
2185 copy_ivec(dd->ntot,nc);
2186 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2187 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2189 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2191 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2197 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2202 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2204 gmx_domdec_comm_t *comm;
2206 int ddindex, nodeid = -1;
2208 comm = cr->dd->comm;
2213 if (comm->bCartesianPP_PME)
2216 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2221 ddindex = dd_index(cr->dd->nc, coords);
2222 if (comm->bCartesianPP)
2224 nodeid = comm->ddindex2simnodeid[ddindex];
2230 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2242 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2245 gmx_domdec_comm_t *comm;
2246 ivec coord, coord_pme;
2253 /* This assumes a uniform x domain decomposition grid cell size */
2254 if (comm->bCartesianPP_PME)
2257 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2258 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2260 /* This is a PP node */
2261 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2262 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2266 else if (comm->bCartesianPP)
2268 if (sim_nodeid < dd->nnodes)
2270 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2275 /* This assumes DD cells with identical x coordinates
2276 * are numbered sequentially.
2278 if (dd->comm->pmenodes == NULL)
2280 if (sim_nodeid < dd->nnodes)
2282 /* The DD index equals the nodeid */
2283 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2289 while (sim_nodeid > dd->comm->pmenodes[i])
2293 if (sim_nodeid < dd->comm->pmenodes[i])
2295 pmenode = dd->comm->pmenodes[i];
2303 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2305 gmx_bool bPMEOnlyNode;
2307 if (DOMAINDECOMP(cr))
2309 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2313 bPMEOnlyNode = FALSE;
2316 return bPMEOnlyNode;
2319 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2320 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2324 ivec coord, coord_pme;
2328 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2331 for (x = 0; x < dd->nc[XX]; x++)
2333 for (y = 0; y < dd->nc[YY]; y++)
2335 for (z = 0; z < dd->nc[ZZ]; z++)
2337 if (dd->comm->bCartesianPP_PME)
2342 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2343 if (dd->ci[XX] == coord_pme[XX] &&
2344 dd->ci[YY] == coord_pme[YY] &&
2345 dd->ci[ZZ] == coord_pme[ZZ])
2347 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2352 /* The slab corresponds to the nodeid in the PME group */
2353 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2355 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2362 /* The last PP-only node is the peer node */
2363 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2367 fprintf(debug, "Receive coordinates from PP nodes:");
2368 for (x = 0; x < *nmy_ddnodes; x++)
2370 fprintf(debug, " %d", (*my_ddnodes)[x]);
2372 fprintf(debug, "\n");
2376 static gmx_bool receive_vir_ener(t_commrec *cr)
2378 gmx_domdec_comm_t *comm;
2379 int pmenode, coords[DIM], rank;
2383 if (cr->npmenodes < cr->dd->nnodes)
2385 comm = cr->dd->comm;
2386 if (comm->bCartesianPP_PME)
2388 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2390 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2391 coords[comm->cartpmedim]++;
2392 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2394 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2395 if (dd_simnode2pmenode(cr, rank) == pmenode)
2397 /* This is not the last PP node for pmenode */
2405 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2406 if (cr->sim_nodeid+1 < cr->nnodes &&
2407 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2409 /* This is not the last PP node for pmenode */
2418 static void set_zones_ncg_home(gmx_domdec_t *dd)
2420 gmx_domdec_zones_t *zones;
2423 zones = &dd->comm->zones;
2425 zones->cg_range[0] = 0;
2426 for (i = 1; i < zones->n+1; i++)
2428 zones->cg_range[i] = dd->ncg_home;
2432 static void rebuild_cgindex(gmx_domdec_t *dd,
2433 const int *gcgs_index, t_state *state)
2435 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2438 dd_cg_gl = dd->index_gl;
2439 cgindex = dd->cgindex;
2442 for (i = 0; i < state->ncg_gl; i++)
2446 dd_cg_gl[i] = cg_gl;
2447 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2451 dd->ncg_home = state->ncg_gl;
2454 set_zones_ncg_home(dd);
2457 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2459 while (cg >= cginfo_mb->cg_end)
2464 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2467 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2468 t_forcerec *fr, char *bLocalCG)
2470 cginfo_mb_t *cginfo_mb;
2476 cginfo_mb = fr->cginfo_mb;
2477 cginfo = fr->cginfo;
2479 for (cg = cg0; cg < cg1; cg++)
2481 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2485 if (bLocalCG != NULL)
2487 for (cg = cg0; cg < cg1; cg++)
2489 bLocalCG[index_gl[cg]] = TRUE;
2494 static void make_dd_indices(gmx_domdec_t *dd,
2495 const int *gcgs_index, int cg_start)
2497 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2498 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2503 bLocalCG = dd->comm->bLocalCG;
2505 if (dd->nat_tot > dd->gatindex_nalloc)
2507 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2508 srenew(dd->gatindex, dd->gatindex_nalloc);
2511 nzone = dd->comm->zones.n;
2512 zone2cg = dd->comm->zones.cg_range;
2513 zone_ncg1 = dd->comm->zone_ncg1;
2514 index_gl = dd->index_gl;
2515 gatindex = dd->gatindex;
2516 bCGs = dd->comm->bCGs;
2518 if (zone2cg[1] != dd->ncg_home)
2520 gmx_incons("dd->ncg_zone is not up to date");
2523 /* Make the local to global and global to local atom index */
2524 a = dd->cgindex[cg_start];
2525 for (zone = 0; zone < nzone; zone++)
2533 cg0 = zone2cg[zone];
2535 cg1 = zone2cg[zone+1];
2536 cg1_p1 = cg0 + zone_ncg1[zone];
2538 for (cg = cg0; cg < cg1; cg++)
2543 /* Signal that this cg is from more than one pulse away */
2546 cg_gl = index_gl[cg];
2549 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2552 ga2la_set(dd->ga2la, a_gl, a, zone1);
2558 gatindex[a] = cg_gl;
2559 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2566 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2569 int ncg, i, ngl, nerr;
2572 if (bLocalCG == NULL)
2576 for (i = 0; i < dd->ncg_tot; i++)
2578 if (!bLocalCG[dd->index_gl[i]])
2581 "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);
2586 for (i = 0; i < ncg_sys; i++)
2593 if (ngl != dd->ncg_tot)
2595 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);
2602 static void check_index_consistency(gmx_domdec_t *dd,
2603 int natoms_sys, int ncg_sys,
2606 int nerr, ngl, i, a, cell;
2611 if (dd->comm->DD_debug > 1)
2613 snew(have, natoms_sys);
2614 for (a = 0; a < dd->nat_tot; a++)
2616 if (have[dd->gatindex[a]] > 0)
2618 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);
2622 have[dd->gatindex[a]] = a + 1;
2628 snew(have, dd->nat_tot);
2631 for (i = 0; i < natoms_sys; i++)
2633 if (ga2la_get(dd->ga2la, i, &a, &cell))
2635 if (a >= dd->nat_tot)
2637 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);
2643 if (dd->gatindex[a] != i)
2645 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);
2652 if (ngl != dd->nat_tot)
2655 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2656 dd->rank, where, ngl, dd->nat_tot);
2658 for (a = 0; a < dd->nat_tot; a++)
2663 "DD node %d, %s: local atom %d, global %d has no global index\n",
2664 dd->rank, where, a+1, dd->gatindex[a]+1);
2669 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2673 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2674 dd->rank, where, nerr);
2678 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2685 /* Clear the whole list without searching */
2686 ga2la_clear(dd->ga2la);
2690 for (i = a_start; i < dd->nat_tot; i++)
2692 ga2la_del(dd->ga2la, dd->gatindex[i]);
2696 bLocalCG = dd->comm->bLocalCG;
2699 for (i = cg_start; i < dd->ncg_tot; i++)
2701 bLocalCG[dd->index_gl[i]] = FALSE;
2705 dd_clear_local_vsite_indices(dd);
2707 if (dd->constraints)
2709 dd_clear_local_constraint_indices(dd);
2713 /* This function should be used for moving the domain boudaries during DLB,
2714 * for obtaining the minimum cell size. It checks the initially set limit
2715 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2716 * and, possibly, a longer cut-off limit set for PME load balancing.
2718 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2722 cellsize_min = comm->cellsize_min[dim];
2724 if (!comm->bVacDLBNoLimit)
2726 /* The cut-off might have changed, e.g. by PME load balacning,
2727 * from the value used to set comm->cellsize_min, so check it.
2729 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2731 if (comm->bPMELoadBalDLBLimits)
2733 /* Check for the cut-off limit set by the PME load balancing */
2734 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2738 return cellsize_min;
2741 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2744 real grid_jump_limit;
2746 /* The distance between the boundaries of cells at distance
2747 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2748 * and by the fact that cells should not be shifted by more than
2749 * half their size, such that cg's only shift by one cell
2750 * at redecomposition.
2752 grid_jump_limit = comm->cellsize_limit;
2753 if (!comm->bVacDLBNoLimit)
2755 if (comm->bPMELoadBalDLBLimits)
2757 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2759 grid_jump_limit = max(grid_jump_limit,
2760 cutoff/comm->cd[dim_ind].np);
2763 return grid_jump_limit;
2766 static gmx_bool check_grid_jump(gmx_large_int_t step,
2772 gmx_domdec_comm_t *comm;
2781 for (d = 1; d < dd->ndim; d++)
2784 limit = grid_jump_limit(comm, cutoff, d);
2785 bfac = ddbox->box_size[dim];
2786 if (ddbox->tric_dir[dim])
2788 bfac *= ddbox->skew_fac[dim];
2790 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2791 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2799 /* This error should never be triggered under normal
2800 * circumstances, but you never know ...
2802 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.",
2803 gmx_step_str(step, buf),
2804 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2812 static int dd_load_count(gmx_domdec_comm_t *comm)
2814 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2817 static float dd_force_load(gmx_domdec_comm_t *comm)
2824 if (comm->eFlop > 1)
2826 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2831 load = comm->cycl[ddCyclF];
2832 if (comm->cycl_n[ddCyclF] > 1)
2834 /* Subtract the maximum of the last n cycle counts
2835 * to get rid of possible high counts due to other soures,
2836 * for instance system activity, that would otherwise
2837 * affect the dynamic load balancing.
2839 load -= comm->cycl_max[ddCyclF];
2846 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2848 gmx_domdec_comm_t *comm;
2853 snew(*dim_f, dd->nc[dim]+1);
2855 for (i = 1; i < dd->nc[dim]; i++)
2857 if (comm->slb_frac[dim])
2859 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2863 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2866 (*dim_f)[dd->nc[dim]] = 1;
2869 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2871 int pmeindex, slab, nso, i;
2874 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2880 ddpme->dim = dimind;
2882 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2884 ddpme->nslab = (ddpme->dim == 0 ?
2885 dd->comm->npmenodes_x :
2886 dd->comm->npmenodes_y);
2888 if (ddpme->nslab <= 1)
2893 nso = dd->comm->npmenodes/ddpme->nslab;
2894 /* Determine for each PME slab the PP location range for dimension dim */
2895 snew(ddpme->pp_min, ddpme->nslab);
2896 snew(ddpme->pp_max, ddpme->nslab);
2897 for (slab = 0; slab < ddpme->nslab; slab++)
2899 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2900 ddpme->pp_max[slab] = 0;
2902 for (i = 0; i < dd->nnodes; i++)
2904 ddindex2xyz(dd->nc, i, xyz);
2905 /* For y only use our y/z slab.
2906 * This assumes that the PME x grid size matches the DD grid size.
2908 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2910 pmeindex = ddindex2pmeindex(dd, i);
2913 slab = pmeindex/nso;
2917 slab = pmeindex % ddpme->nslab;
2919 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2920 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2924 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2927 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2929 if (dd->comm->ddpme[0].dim == XX)
2931 return dd->comm->ddpme[0].maxshift;
2939 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2941 if (dd->comm->ddpme[0].dim == YY)
2943 return dd->comm->ddpme[0].maxshift;
2945 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2947 return dd->comm->ddpme[1].maxshift;
2955 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2956 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2958 gmx_domdec_comm_t *comm;
2961 real range, pme_boundary;
2965 nc = dd->nc[ddpme->dim];
2968 if (!ddpme->dim_match)
2970 /* PP decomposition is not along dim: the worst situation */
2973 else if (ns <= 3 || (bUniform && ns == nc))
2975 /* The optimal situation */
2980 /* We need to check for all pme nodes which nodes they
2981 * could possibly need to communicate with.
2983 xmin = ddpme->pp_min;
2984 xmax = ddpme->pp_max;
2985 /* Allow for atoms to be maximally 2/3 times the cut-off
2986 * out of their DD cell. This is a reasonable balance between
2987 * between performance and support for most charge-group/cut-off
2990 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2991 /* Avoid extra communication when we are exactly at a boundary */
2995 for (s = 0; s < ns; s++)
2997 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2998 pme_boundary = (real)s/ns;
3001 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3003 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3007 pme_boundary = (real)(s+1)/ns;
3010 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3012 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3019 ddpme->maxshift = sh;
3023 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3024 ddpme->dim, ddpme->maxshift);
3028 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3032 for (d = 0; d < dd->ndim; d++)
3035 if (dim < ddbox->nboundeddim &&
3036 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3037 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3039 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",
3040 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3041 dd->nc[dim], dd->comm->cellsize_limit);
3046 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3047 gmx_bool bMaster, ivec npulse)
3049 gmx_domdec_comm_t *comm;
3052 real *cell_x, cell_dx, cellsize;
3056 for (d = 0; d < DIM; d++)
3058 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3060 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3063 cell_dx = ddbox->box_size[d]/dd->nc[d];
3066 for (j = 0; j < dd->nc[d]+1; j++)
3068 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3073 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3074 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3076 cellsize = cell_dx*ddbox->skew_fac[d];
3077 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3081 cellsize_min[d] = cellsize;
3085 /* Statically load balanced grid */
3086 /* Also when we are not doing a master distribution we determine
3087 * all cell borders in a loop to obtain identical values
3088 * to the master distribution case and to determine npulse.
3092 cell_x = dd->ma->cell_x[d];
3096 snew(cell_x, dd->nc[d]+1);
3098 cell_x[0] = ddbox->box0[d];
3099 for (j = 0; j < dd->nc[d]; j++)
3101 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3102 cell_x[j+1] = cell_x[j] + cell_dx;
3103 cellsize = cell_dx*ddbox->skew_fac[d];
3104 while (cellsize*npulse[d] < comm->cutoff &&
3105 npulse[d] < dd->nc[d]-1)
3109 cellsize_min[d] = min(cellsize_min[d], cellsize);
3113 comm->cell_x0[d] = cell_x[dd->ci[d]];
3114 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3118 /* The following limitation is to avoid that a cell would receive
3119 * some of its own home charge groups back over the periodic boundary.
3120 * Double charge groups cause trouble with the global indices.
3122 if (d < ddbox->npbcdim &&
3123 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3125 gmx_fatal_collective(FARGS, NULL, dd,
3126 "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",
3127 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3129 dd->nc[d], dd->nc[d],
3130 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3134 if (!comm->bDynLoadBal)
3136 copy_rvec(cellsize_min, comm->cellsize_min);
3139 for (d = 0; d < comm->npmedecompdim; d++)
3141 set_pme_maxshift(dd, &comm->ddpme[d],
3142 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3143 comm->ddpme[d].slb_dim_f);
3148 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3149 int d, int dim, gmx_domdec_root_t *root,
3151 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3153 gmx_domdec_comm_t *comm;
3154 int ncd, i, j, nmin, nmin_old;
3155 gmx_bool bLimLo, bLimHi;
3157 real fac, halfway, cellsize_limit_f_i, region_size;
3158 gmx_bool bPBC, bLastHi = FALSE;
3159 int nrange[] = {range[0], range[1]};
3161 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3167 bPBC = (dim < ddbox->npbcdim);
3169 cell_size = root->buf_ncd;
3173 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3176 /* First we need to check if the scaling does not make cells
3177 * smaller than the smallest allowed size.
3178 * We need to do this iteratively, since if a cell is too small,
3179 * it needs to be enlarged, which makes all the other cells smaller,
3180 * which could in turn make another cell smaller than allowed.
3182 for (i = range[0]; i < range[1]; i++)
3184 root->bCellMin[i] = FALSE;
3190 /* We need the total for normalization */
3192 for (i = range[0]; i < range[1]; i++)
3194 if (root->bCellMin[i] == FALSE)
3196 fac += cell_size[i];
3199 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3200 /* Determine the cell boundaries */
3201 for (i = range[0]; i < range[1]; i++)
3203 if (root->bCellMin[i] == FALSE)
3205 cell_size[i] *= fac;
3206 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3208 cellsize_limit_f_i = 0;
3212 cellsize_limit_f_i = cellsize_limit_f;
3214 if (cell_size[i] < cellsize_limit_f_i)
3216 root->bCellMin[i] = TRUE;
3217 cell_size[i] = cellsize_limit_f_i;
3221 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3224 while (nmin > nmin_old);
3227 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3228 /* For this check we should not use DD_CELL_MARGIN,
3229 * but a slightly smaller factor,
3230 * since rounding could get use below the limit.
3232 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3235 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",
3236 gmx_step_str(step, buf),
3237 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3238 ncd, comm->cellsize_min[dim]);
3241 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3245 /* Check if the boundary did not displace more than halfway
3246 * each of the cells it bounds, as this could cause problems,
3247 * especially when the differences between cell sizes are large.
3248 * If changes are applied, they will not make cells smaller
3249 * than the cut-off, as we check all the boundaries which
3250 * might be affected by a change and if the old state was ok,
3251 * the cells will at most be shrunk back to their old size.
3253 for (i = range[0]+1; i < range[1]; i++)
3255 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3256 if (root->cell_f[i] < halfway)
3258 root->cell_f[i] = halfway;
3259 /* Check if the change also causes shifts of the next boundaries */
3260 for (j = i+1; j < range[1]; j++)
3262 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3264 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3268 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3269 if (root->cell_f[i] > halfway)
3271 root->cell_f[i] = halfway;
3272 /* Check if the change also causes shifts of the next boundaries */
3273 for (j = i-1; j >= range[0]+1; j--)
3275 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3277 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3284 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3285 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3286 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3287 * for a and b nrange is used */
3290 /* Take care of the staggering of the cell boundaries */
3293 for (i = range[0]; i < range[1]; i++)
3295 root->cell_f_max0[i] = root->cell_f[i];
3296 root->cell_f_min1[i] = root->cell_f[i+1];
3301 for (i = range[0]+1; i < range[1]; i++)
3303 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3304 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3305 if (bLimLo && bLimHi)
3307 /* Both limits violated, try the best we can */
3308 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3309 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3310 nrange[0] = range[0];
3312 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3315 nrange[1] = range[1];
3316 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3322 /* root->cell_f[i] = root->bound_min[i]; */
3323 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3326 else if (bLimHi && !bLastHi)
3329 if (nrange[1] < range[1]) /* found a LimLo before */
3331 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3332 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3333 nrange[0] = nrange[1];
3335 root->cell_f[i] = root->bound_max[i];
3337 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3339 nrange[1] = range[1];
3342 if (nrange[1] < range[1]) /* found last a LimLo */
3344 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3345 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3346 nrange[0] = nrange[1];
3347 nrange[1] = range[1];
3348 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3350 else if (nrange[0] > range[0]) /* found at least one LimHi */
3352 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3359 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3360 int d, int dim, gmx_domdec_root_t *root,
3361 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3362 gmx_bool bUniform, gmx_large_int_t step)
3364 gmx_domdec_comm_t *comm;
3365 int ncd, d1, i, j, pos;
3367 real load_aver, load_i, imbalance, change, change_max, sc;
3368 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3372 int range[] = { 0, 0 };
3376 /* Convert the maximum change from the input percentage to a fraction */
3377 change_limit = comm->dlb_scale_lim*0.01;
3381 bPBC = (dim < ddbox->npbcdim);
3383 cell_size = root->buf_ncd;
3385 /* Store the original boundaries */
3386 for (i = 0; i < ncd+1; i++)
3388 root->old_cell_f[i] = root->cell_f[i];
3392 for (i = 0; i < ncd; i++)
3394 cell_size[i] = 1.0/ncd;
3397 else if (dd_load_count(comm))
3399 load_aver = comm->load[d].sum_m/ncd;
3401 for (i = 0; i < ncd; i++)
3403 /* Determine the relative imbalance of cell i */
3404 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3405 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3406 /* Determine the change of the cell size using underrelaxation */
3407 change = -relax*imbalance;
3408 change_max = max(change_max, max(change, -change));
3410 /* Limit the amount of scaling.
3411 * We need to use the same rescaling for all cells in one row,
3412 * otherwise the load balancing might not converge.
3415 if (change_max > change_limit)
3417 sc *= change_limit/change_max;
3419 for (i = 0; i < ncd; i++)
3421 /* Determine the relative imbalance of cell i */
3422 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3423 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3424 /* Determine the change of the cell size using underrelaxation */
3425 change = -sc*imbalance;
3426 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3430 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3431 cellsize_limit_f *= DD_CELL_MARGIN;
3432 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3433 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3434 if (ddbox->tric_dir[dim])
3436 cellsize_limit_f /= ddbox->skew_fac[dim];
3437 dist_min_f /= ddbox->skew_fac[dim];
3439 if (bDynamicBox && d > 0)
3441 dist_min_f *= DD_PRES_SCALE_MARGIN;
3443 if (d > 0 && !bUniform)
3445 /* Make sure that the grid is not shifted too much */
3446 for (i = 1; i < ncd; i++)
3448 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3450 gmx_incons("Inconsistent DD boundary staggering limits!");
3452 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3453 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3456 root->bound_min[i] += 0.5*space;
3458 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3459 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3462 root->bound_max[i] += 0.5*space;
3467 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3469 root->cell_f_max0[i-1] + dist_min_f,
3470 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3471 root->cell_f_min1[i] - dist_min_f);
3476 root->cell_f[0] = 0;
3477 root->cell_f[ncd] = 1;
3478 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3481 /* After the checks above, the cells should obey the cut-off
3482 * restrictions, but it does not hurt to check.
3484 for (i = 0; i < ncd; i++)
3488 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3489 dim, i, root->cell_f[i], root->cell_f[i+1]);
3492 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3493 root->cell_f[i+1] - root->cell_f[i] <
3494 cellsize_limit_f/DD_CELL_MARGIN)
3498 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3499 gmx_step_str(step, buf), dim2char(dim), i,
3500 (root->cell_f[i+1] - root->cell_f[i])
3501 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3506 /* Store the cell boundaries of the lower dimensions at the end */
3507 for (d1 = 0; d1 < d; d1++)
3509 root->cell_f[pos++] = comm->cell_f0[d1];
3510 root->cell_f[pos++] = comm->cell_f1[d1];
3513 if (d < comm->npmedecompdim)
3515 /* The master determines the maximum shift for
3516 * the coordinate communication between separate PME nodes.
3518 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3520 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3523 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3527 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3528 gmx_ddbox_t *ddbox, int dimind)
3530 gmx_domdec_comm_t *comm;
3535 /* Set the cell dimensions */
3536 dim = dd->dim[dimind];
3537 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3538 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3539 if (dim >= ddbox->nboundeddim)
3541 comm->cell_x0[dim] += ddbox->box0[dim];
3542 comm->cell_x1[dim] += ddbox->box0[dim];
3546 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3547 int d, int dim, real *cell_f_row,
3550 gmx_domdec_comm_t *comm;
3556 /* Each node would only need to know two fractions,
3557 * but it is probably cheaper to broadcast the whole array.
3559 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3560 0, comm->mpi_comm_load[d]);
3562 /* Copy the fractions for this dimension from the buffer */
3563 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3564 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3565 /* The whole array was communicated, so set the buffer position */
3566 pos = dd->nc[dim] + 1;
3567 for (d1 = 0; d1 <= d; d1++)
3571 /* Copy the cell fractions of the lower dimensions */
3572 comm->cell_f0[d1] = cell_f_row[pos++];
3573 comm->cell_f1[d1] = cell_f_row[pos++];
3575 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3577 /* Convert the communicated shift from float to int */
3578 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3581 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3585 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3586 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3587 gmx_bool bUniform, gmx_large_int_t step)
3589 gmx_domdec_comm_t *comm;
3591 gmx_bool bRowMember, bRowRoot;
3596 for (d = 0; d < dd->ndim; d++)
3601 for (d1 = d; d1 < dd->ndim; d1++)
3603 if (dd->ci[dd->dim[d1]] > 0)
3616 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3617 ddbox, bDynamicBox, bUniform, step);
3618 cell_f_row = comm->root[d]->cell_f;
3622 cell_f_row = comm->cell_f_row;
3624 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3629 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3633 /* This function assumes the box is static and should therefore
3634 * not be called when the box has changed since the last
3635 * call to dd_partition_system.
3637 for (d = 0; d < dd->ndim; d++)
3639 relative_to_absolute_cell_bounds(dd, ddbox, d);
3645 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3646 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3647 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3648 gmx_wallcycle_t wcycle)
3650 gmx_domdec_comm_t *comm;
3657 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3658 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3659 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3661 else if (bDynamicBox)
3663 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3666 /* Set the dimensions for which no DD is used */
3667 for (dim = 0; dim < DIM; dim++)
3669 if (dd->nc[dim] == 1)
3671 comm->cell_x0[dim] = 0;
3672 comm->cell_x1[dim] = ddbox->box_size[dim];
3673 if (dim >= ddbox->nboundeddim)
3675 comm->cell_x0[dim] += ddbox->box0[dim];
3676 comm->cell_x1[dim] += ddbox->box0[dim];
3682 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3685 gmx_domdec_comm_dim_t *cd;
3687 for (d = 0; d < dd->ndim; d++)
3689 cd = &dd->comm->cd[d];
3690 np = npulse[dd->dim[d]];
3691 if (np > cd->np_nalloc)
3695 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3696 dim2char(dd->dim[d]), np);
3698 if (DDMASTER(dd) && cd->np_nalloc > 0)
3700 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3702 srenew(cd->ind, np);
3703 for (i = cd->np_nalloc; i < np; i++)
3705 cd->ind[i].index = NULL;
3706 cd->ind[i].nalloc = 0;
3715 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3716 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3717 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3718 gmx_wallcycle_t wcycle)
3720 gmx_domdec_comm_t *comm;
3726 /* Copy the old cell boundaries for the cg displacement check */
3727 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3728 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3730 if (comm->bDynLoadBal)
3734 check_box_size(dd, ddbox);
3736 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3740 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3741 realloc_comm_ind(dd, npulse);
3746 for (d = 0; d < DIM; d++)
3748 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3749 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3754 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3756 rvec cell_ns_x0, rvec cell_ns_x1,
3757 gmx_large_int_t step)
3759 gmx_domdec_comm_t *comm;
3764 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3766 dim = dd->dim[dim_ind];
3768 /* Without PBC we don't have restrictions on the outer cells */
3769 if (!(dim >= ddbox->npbcdim &&
3770 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3771 comm->bDynLoadBal &&
3772 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3773 comm->cellsize_min[dim])
3776 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",
3777 gmx_step_str(step, buf), dim2char(dim),
3778 comm->cell_x1[dim] - comm->cell_x0[dim],
3779 ddbox->skew_fac[dim],
3780 dd->comm->cellsize_min[dim],
3781 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3785 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3787 /* Communicate the boundaries and update cell_ns_x0/1 */
3788 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3789 if (dd->bGridJump && dd->ndim > 1)
3791 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3796 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3800 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3808 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3809 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3818 static void check_screw_box(matrix box)
3820 /* Mathematical limitation */
3821 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3823 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3826 /* Limitation due to the asymmetry of the eighth shell method */
3827 if (box[ZZ][YY] != 0)
3829 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3833 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3834 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3837 gmx_domdec_master_t *ma;
3838 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3839 int i, icg, j, k, k0, k1, d, npbcdim;
3841 rvec box_size, cg_cm;
3843 real nrcg, inv_ncg, pos_d;
3845 gmx_bool bUnbounded, bScrew;
3849 if (tmp_ind == NULL)
3851 snew(tmp_nalloc, dd->nnodes);
3852 snew(tmp_ind, dd->nnodes);
3853 for (i = 0; i < dd->nnodes; i++)
3855 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3856 snew(tmp_ind[i], tmp_nalloc[i]);
3860 /* Clear the count */
3861 for (i = 0; i < dd->nnodes; i++)
3867 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3869 cgindex = cgs->index;
3871 /* Compute the center of geometry for all charge groups */
3872 for (icg = 0; icg < cgs->nr; icg++)
3875 k1 = cgindex[icg+1];
3879 copy_rvec(pos[k0], cg_cm);
3886 for (k = k0; (k < k1); k++)
3888 rvec_inc(cg_cm, pos[k]);
3890 for (d = 0; (d < DIM); d++)
3892 cg_cm[d] *= inv_ncg;
3895 /* Put the charge group in the box and determine the cell index */
3896 for (d = DIM-1; d >= 0; d--)
3899 if (d < dd->npbcdim)
3901 bScrew = (dd->bScrewPBC && d == XX);
3902 if (tric_dir[d] && dd->nc[d] > 1)
3904 /* Use triclinic coordintates for this dimension */
3905 for (j = d+1; j < DIM; j++)
3907 pos_d += cg_cm[j]*tcm[j][d];
3910 while (pos_d >= box[d][d])
3913 rvec_dec(cg_cm, box[d]);
3916 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3917 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3919 for (k = k0; (k < k1); k++)
3921 rvec_dec(pos[k], box[d]);
3924 pos[k][YY] = box[YY][YY] - pos[k][YY];
3925 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3932 rvec_inc(cg_cm, box[d]);
3935 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3936 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3938 for (k = k0; (k < k1); k++)
3940 rvec_inc(pos[k], box[d]);
3943 pos[k][YY] = box[YY][YY] - pos[k][YY];
3944 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3949 /* This could be done more efficiently */
3951 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3956 i = dd_index(dd->nc, ind);
3957 if (ma->ncg[i] == tmp_nalloc[i])
3959 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3960 srenew(tmp_ind[i], tmp_nalloc[i]);
3962 tmp_ind[i][ma->ncg[i]] = icg;
3964 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3968 for (i = 0; i < dd->nnodes; i++)
3971 for (k = 0; k < ma->ncg[i]; k++)
3973 ma->cg[k1++] = tmp_ind[i][k];
3976 ma->index[dd->nnodes] = k1;
3978 for (i = 0; i < dd->nnodes; i++)
3988 fprintf(fplog, "Charge group distribution at step %s:",
3989 gmx_step_str(step, buf));
3990 for (i = 0; i < dd->nnodes; i++)
3992 fprintf(fplog, " %d", ma->ncg[i]);
3994 fprintf(fplog, "\n");
3998 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
3999 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4002 gmx_domdec_master_t *ma = NULL;
4005 int *ibuf, buf2[2] = { 0, 0 };
4006 gmx_bool bMaster = DDMASTER(dd);
4013 check_screw_box(box);
4016 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4018 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4019 for (i = 0; i < dd->nnodes; i++)
4021 ma->ibuf[2*i] = ma->ncg[i];
4022 ma->ibuf[2*i+1] = ma->nat[i];
4030 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4032 dd->ncg_home = buf2[0];
4033 dd->nat_home = buf2[1];
4034 dd->ncg_tot = dd->ncg_home;
4035 dd->nat_tot = dd->nat_home;
4036 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4038 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4039 srenew(dd->index_gl, dd->cg_nalloc);
4040 srenew(dd->cgindex, dd->cg_nalloc+1);
4044 for (i = 0; i < dd->nnodes; i++)
4046 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4047 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4052 DDMASTER(dd) ? ma->ibuf : NULL,
4053 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4054 DDMASTER(dd) ? ma->cg : NULL,
4055 dd->ncg_home*sizeof(int), dd->index_gl);
4057 /* Determine the home charge group sizes */
4059 for (i = 0; i < dd->ncg_home; i++)
4061 cg_gl = dd->index_gl[i];
4063 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4068 fprintf(debug, "Home charge groups:\n");
4069 for (i = 0; i < dd->ncg_home; i++)
4071 fprintf(debug, " %d", dd->index_gl[i]);
4074 fprintf(debug, "\n");
4077 fprintf(debug, "\n");
4081 static int compact_and_copy_vec_at(int ncg, int *move,
4084 rvec *src, gmx_domdec_comm_t *comm,
4087 int m, icg, i, i0, i1, nrcg;
4093 for (m = 0; m < DIM*2; m++)
4099 for (icg = 0; icg < ncg; icg++)
4101 i1 = cgindex[icg+1];
4107 /* Compact the home array in place */
4108 for (i = i0; i < i1; i++)
4110 copy_rvec(src[i], src[home_pos++]);
4116 /* Copy to the communication buffer */
4118 pos_vec[m] += 1 + vec*nrcg;
4119 for (i = i0; i < i1; i++)
4121 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4123 pos_vec[m] += (nvec - vec - 1)*nrcg;
4127 home_pos += i1 - i0;
4135 static int compact_and_copy_vec_cg(int ncg, int *move,
4137 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4140 int m, icg, i0, i1, nrcg;
4146 for (m = 0; m < DIM*2; m++)
4152 for (icg = 0; icg < ncg; icg++)
4154 i1 = cgindex[icg+1];
4160 /* Compact the home array in place */
4161 copy_rvec(src[icg], src[home_pos++]);
4167 /* Copy to the communication buffer */
4168 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4169 pos_vec[m] += 1 + nrcg*nvec;
4181 static int compact_ind(int ncg, int *move,
4182 int *index_gl, int *cgindex,
4184 gmx_ga2la_t ga2la, char *bLocalCG,
4187 int cg, nat, a0, a1, a, a_gl;
4192 for (cg = 0; cg < ncg; cg++)
4198 /* Compact the home arrays in place.
4199 * Anything that can be done here avoids access to global arrays.
4201 cgindex[home_pos] = nat;
4202 for (a = a0; a < a1; a++)
4205 gatindex[nat] = a_gl;
4206 /* The cell number stays 0, so we don't need to set it */
4207 ga2la_change_la(ga2la, a_gl, nat);
4210 index_gl[home_pos] = index_gl[cg];
4211 cginfo[home_pos] = cginfo[cg];
4212 /* The charge group remains local, so bLocalCG does not change */
4217 /* Clear the global indices */
4218 for (a = a0; a < a1; a++)
4220 ga2la_del(ga2la, gatindex[a]);
4224 bLocalCG[index_gl[cg]] = FALSE;
4228 cgindex[home_pos] = nat;
4233 static void clear_and_mark_ind(int ncg, int *move,
4234 int *index_gl, int *cgindex, int *gatindex,
4235 gmx_ga2la_t ga2la, char *bLocalCG,
4240 for (cg = 0; cg < ncg; cg++)
4246 /* Clear the global indices */
4247 for (a = a0; a < a1; a++)
4249 ga2la_del(ga2la, gatindex[a]);
4253 bLocalCG[index_gl[cg]] = FALSE;
4255 /* Signal that this cg has moved using the ns cell index.
4256 * Here we set it to -1. fill_grid will change it
4257 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4259 cell_index[cg] = -1;
4264 static void print_cg_move(FILE *fplog,
4266 gmx_large_int_t step, int cg, int dim, int dir,
4267 gmx_bool bHaveLimitdAndCMOld, real limitd,
4268 rvec cm_old, rvec cm_new, real pos_d)
4270 gmx_domdec_comm_t *comm;
4275 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4276 if (bHaveLimitdAndCMOld)
4278 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4279 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4283 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4284 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4286 fprintf(fplog, "distance out of cell %f\n",
4287 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4288 if (bHaveLimitdAndCMOld)
4290 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4291 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4293 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4294 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4295 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4297 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4298 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4300 comm->cell_x0[dim], comm->cell_x1[dim]);
4303 static void cg_move_error(FILE *fplog,
4305 gmx_large_int_t step, int cg, int dim, int dir,
4306 gmx_bool bHaveLimitdAndCMOld, real limitd,
4307 rvec cm_old, rvec cm_new, real pos_d)
4311 print_cg_move(fplog, dd, step, cg, dim, dir,
4312 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4314 print_cg_move(stderr, dd, step, cg, dim, dir,
4315 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4317 "A charge group moved too far between two domain decomposition steps\n"
4318 "This usually means that your system is not well equilibrated");
4321 static void rotate_state_atom(t_state *state, int a)
4325 for (est = 0; est < estNR; est++)
4327 if (EST_DISTR(est) && (state->flags & (1<<est)))
4332 /* Rotate the complete state; for a rectangular box only */
4333 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4334 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4337 state->v[a][YY] = -state->v[a][YY];
4338 state->v[a][ZZ] = -state->v[a][ZZ];
4341 state->sd_X[a][YY] = -state->sd_X[a][YY];
4342 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4345 state->cg_p[a][YY] = -state->cg_p[a][YY];
4346 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4348 case estDISRE_INITF:
4349 case estDISRE_RM3TAV:
4350 case estORIRE_INITF:
4352 /* These are distances, so not affected by rotation */
4355 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4361 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4363 if (natoms > comm->moved_nalloc)
4365 /* Contents should be preserved here */
4366 comm->moved_nalloc = over_alloc_dd(natoms);
4367 srenew(comm->moved, comm->moved_nalloc);
4373 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4376 ivec tric_dir, matrix tcm,
4377 rvec cell_x0, rvec cell_x1,
4378 rvec limitd, rvec limit0, rvec limit1,
4380 int cg_start, int cg_end,
4385 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4386 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4390 real inv_ncg, pos_d;
4393 npbcdim = dd->npbcdim;
4395 for (cg = cg_start; cg < cg_end; cg++)
4402 copy_rvec(state->x[k0], cm_new);
4409 for (k = k0; (k < k1); k++)
4411 rvec_inc(cm_new, state->x[k]);
4413 for (d = 0; (d < DIM); d++)
4415 cm_new[d] = inv_ncg*cm_new[d];
4420 /* Do pbc and check DD cell boundary crossings */
4421 for (d = DIM-1; d >= 0; d--)
4425 bScrew = (dd->bScrewPBC && d == XX);
4426 /* Determine the location of this cg in lattice coordinates */
4430 for (d2 = d+1; d2 < DIM; d2++)
4432 pos_d += cm_new[d2]*tcm[d2][d];
4435 /* Put the charge group in the triclinic unit-cell */
4436 if (pos_d >= cell_x1[d])
4438 if (pos_d >= limit1[d])
4440 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4441 cg_cm[cg], cm_new, pos_d);
4444 if (dd->ci[d] == dd->nc[d] - 1)
4446 rvec_dec(cm_new, state->box[d]);
4449 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4450 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4452 for (k = k0; (k < k1); k++)
4454 rvec_dec(state->x[k], state->box[d]);
4457 rotate_state_atom(state, k);
4462 else if (pos_d < cell_x0[d])
4464 if (pos_d < limit0[d])
4466 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4467 cg_cm[cg], cm_new, pos_d);
4472 rvec_inc(cm_new, state->box[d]);
4475 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4476 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4478 for (k = k0; (k < k1); k++)
4480 rvec_inc(state->x[k], state->box[d]);
4483 rotate_state_atom(state, k);
4489 else if (d < npbcdim)
4491 /* Put the charge group in the rectangular unit-cell */
4492 while (cm_new[d] >= state->box[d][d])
4494 rvec_dec(cm_new, state->box[d]);
4495 for (k = k0; (k < k1); k++)
4497 rvec_dec(state->x[k], state->box[d]);
4500 while (cm_new[d] < 0)
4502 rvec_inc(cm_new, state->box[d]);
4503 for (k = k0; (k < k1); k++)
4505 rvec_inc(state->x[k], state->box[d]);
4511 copy_rvec(cm_new, cg_cm[cg]);
4513 /* Determine where this cg should go */
4516 for (d = 0; d < dd->ndim; d++)
4521 flag |= DD_FLAG_FW(d);
4527 else if (dev[dim] == -1)
4529 flag |= DD_FLAG_BW(d);
4532 if (dd->nc[dim] > 2)
4543 /* Temporarily store the flag in move */
4544 move[cg] = mc + flag;
4548 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4549 gmx_domdec_t *dd, ivec tric_dir,
4550 t_state *state, rvec **f,
4551 t_forcerec *fr, t_mdatoms *md,
4559 int ncg[DIM*2], nat[DIM*2];
4560 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4561 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4562 int sbuf[2], rbuf[2];
4563 int home_pos_cg, home_pos_at, buf_pos;
4565 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4568 real inv_ncg, pos_d;
4570 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4572 cginfo_mb_t *cginfo_mb;
4573 gmx_domdec_comm_t *comm;
4575 int nthread, thread;
4579 check_screw_box(state->box);
4583 if (fr->cutoff_scheme == ecutsGROUP)
4588 for (i = 0; i < estNR; i++)
4594 case estX: /* Always present */ break;
4595 case estV: bV = (state->flags & (1<<i)); break;
4596 case estSDX: bSDX = (state->flags & (1<<i)); break;
4597 case estCGP: bCGP = (state->flags & (1<<i)); break;
4600 case estDISRE_INITF:
4601 case estDISRE_RM3TAV:
4602 case estORIRE_INITF:
4604 /* No processing required */
4607 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4612 if (dd->ncg_tot > comm->nalloc_int)
4614 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4615 srenew(comm->buf_int, comm->nalloc_int);
4617 move = comm->buf_int;
4619 /* Clear the count */
4620 for (c = 0; c < dd->ndim*2; c++)
4626 npbcdim = dd->npbcdim;
4628 for (d = 0; (d < DIM); d++)
4630 limitd[d] = dd->comm->cellsize_min[d];
4631 if (d >= npbcdim && dd->ci[d] == 0)
4633 cell_x0[d] = -GMX_FLOAT_MAX;
4637 cell_x0[d] = comm->cell_x0[d];
4639 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4641 cell_x1[d] = GMX_FLOAT_MAX;
4645 cell_x1[d] = comm->cell_x1[d];
4649 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4650 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4654 /* We check after communication if a charge group moved
4655 * more than one cell. Set the pre-comm check limit to float_max.
4657 limit0[d] = -GMX_FLOAT_MAX;
4658 limit1[d] = GMX_FLOAT_MAX;
4662 make_tric_corr_matrix(npbcdim, state->box, tcm);
4664 cgindex = dd->cgindex;
4666 nthread = gmx_omp_nthreads_get(emntDomdec);
4668 /* Compute the center of geometry for all home charge groups
4669 * and put them in the box and determine where they should go.
4671 #pragma omp parallel for num_threads(nthread) schedule(static)
4672 for (thread = 0; thread < nthread; thread++)
4674 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4675 cell_x0, cell_x1, limitd, limit0, limit1,
4677 ( thread *dd->ncg_home)/nthread,
4678 ((thread+1)*dd->ncg_home)/nthread,
4679 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4683 for (cg = 0; cg < dd->ncg_home; cg++)
4688 flag = mc & ~DD_FLAG_NRCG;
4689 mc = mc & DD_FLAG_NRCG;
4692 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4694 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4695 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4697 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4698 /* We store the cg size in the lower 16 bits
4699 * and the place where the charge group should go
4700 * in the next 6 bits. This saves some communication volume.
4702 nrcg = cgindex[cg+1] - cgindex[cg];
4703 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4709 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4710 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4713 for (i = 0; i < dd->ndim*2; i++)
4715 *ncg_moved += ncg[i];
4732 /* Make sure the communication buffers are large enough */
4733 for (mc = 0; mc < dd->ndim*2; mc++)
4735 nvr = ncg[mc] + nat[mc]*nvec;
4736 if (nvr > comm->cgcm_state_nalloc[mc])
4738 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4739 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4743 switch (fr->cutoff_scheme)
4746 /* Recalculating cg_cm might be cheaper than communicating,
4747 * but that could give rise to rounding issues.
4750 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4751 nvec, cg_cm, comm, bCompact);
4754 /* Without charge groups we send the moved atom coordinates
4755 * over twice. This is so the code below can be used without
4756 * many conditionals for both for with and without charge groups.
4759 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4760 nvec, state->x, comm, FALSE);
4763 home_pos_cg -= *ncg_moved;
4767 gmx_incons("unimplemented");
4773 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4774 nvec, vec++, state->x, comm, bCompact);
4777 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4778 nvec, vec++, state->v, comm, bCompact);
4782 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4783 nvec, vec++, state->sd_X, comm, bCompact);
4787 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4788 nvec, vec++, state->cg_p, comm, bCompact);
4793 compact_ind(dd->ncg_home, move,
4794 dd->index_gl, dd->cgindex, dd->gatindex,
4795 dd->ga2la, comm->bLocalCG,
4800 if (fr->cutoff_scheme == ecutsVERLET)
4802 moved = get_moved(comm, dd->ncg_home);
4804 for (k = 0; k < dd->ncg_home; k++)
4811 moved = fr->ns.grid->cell_index;
4814 clear_and_mark_ind(dd->ncg_home, move,
4815 dd->index_gl, dd->cgindex, dd->gatindex,
4816 dd->ga2la, comm->bLocalCG,
4820 cginfo_mb = fr->cginfo_mb;
4822 *ncg_stay_home = home_pos_cg;
4823 for (d = 0; d < dd->ndim; d++)
4829 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4832 /* Communicate the cg and atom counts */
4837 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4838 d, dir, sbuf[0], sbuf[1]);
4840 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4842 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4844 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4845 srenew(comm->buf_int, comm->nalloc_int);
4848 /* Communicate the charge group indices, sizes and flags */
4849 dd_sendrecv_int(dd, d, dir,
4850 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4851 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4853 nvs = ncg[cdd] + nat[cdd]*nvec;
4854 i = rbuf[0] + rbuf[1] *nvec;
4855 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4857 /* Communicate cgcm and state */
4858 dd_sendrecv_rvec(dd, d, dir,
4859 comm->cgcm_state[cdd], nvs,
4860 comm->vbuf.v+nvr, i);
4861 ncg_recv += rbuf[0];
4862 nat_recv += rbuf[1];
4866 /* Process the received charge groups */
4868 for (cg = 0; cg < ncg_recv; cg++)
4870 flag = comm->buf_int[cg*DD_CGIBS+1];
4872 if (dim >= npbcdim && dd->nc[dim] > 2)
4874 /* No pbc in this dim and more than one domain boundary.
4875 * We do a separate check if a charge group didn't move too far.
4877 if (((flag & DD_FLAG_FW(d)) &&
4878 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4879 ((flag & DD_FLAG_BW(d)) &&
4880 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4882 cg_move_error(fplog, dd, step, cg, dim,
4883 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4885 comm->vbuf.v[buf_pos],
4886 comm->vbuf.v[buf_pos],
4887 comm->vbuf.v[buf_pos][dim]);
4894 /* Check which direction this cg should go */
4895 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4899 /* The cell boundaries for dimension d2 are not equal
4900 * for each cell row of the lower dimension(s),
4901 * therefore we might need to redetermine where
4902 * this cg should go.
4905 /* If this cg crosses the box boundary in dimension d2
4906 * we can use the communicated flag, so we do not
4907 * have to worry about pbc.
4909 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4910 (flag & DD_FLAG_FW(d2))) ||
4911 (dd->ci[dim2] == 0 &&
4912 (flag & DD_FLAG_BW(d2)))))
4914 /* Clear the two flags for this dimension */
4915 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4916 /* Determine the location of this cg
4917 * in lattice coordinates
4919 pos_d = comm->vbuf.v[buf_pos][dim2];
4922 for (d3 = dim2+1; d3 < DIM; d3++)
4925 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4928 /* Check of we are not at the box edge.
4929 * pbc is only handled in the first step above,
4930 * but this check could move over pbc while
4931 * the first step did not due to different rounding.
4933 if (pos_d >= cell_x1[dim2] &&
4934 dd->ci[dim2] != dd->nc[dim2]-1)
4936 flag |= DD_FLAG_FW(d2);
4938 else if (pos_d < cell_x0[dim2] &&
4941 flag |= DD_FLAG_BW(d2);
4943 comm->buf_int[cg*DD_CGIBS+1] = flag;
4946 /* Set to which neighboring cell this cg should go */
4947 if (flag & DD_FLAG_FW(d2))
4951 else if (flag & DD_FLAG_BW(d2))
4953 if (dd->nc[dd->dim[d2]] > 2)
4965 nrcg = flag & DD_FLAG_NRCG;
4968 if (home_pos_cg+1 > dd->cg_nalloc)
4970 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4971 srenew(dd->index_gl, dd->cg_nalloc);
4972 srenew(dd->cgindex, dd->cg_nalloc+1);
4974 /* Set the global charge group index and size */
4975 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4976 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4977 /* Copy the state from the buffer */
4978 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4979 if (fr->cutoff_scheme == ecutsGROUP)
4982 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4986 /* Set the cginfo */
4987 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4988 dd->index_gl[home_pos_cg]);
4991 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4994 if (home_pos_at+nrcg > state->nalloc)
4996 dd_realloc_state(state, f, home_pos_at+nrcg);
4998 for (i = 0; i < nrcg; i++)
5000 copy_rvec(comm->vbuf.v[buf_pos++],
5001 state->x[home_pos_at+i]);
5005 for (i = 0; i < nrcg; i++)
5007 copy_rvec(comm->vbuf.v[buf_pos++],
5008 state->v[home_pos_at+i]);
5013 for (i = 0; i < nrcg; i++)
5015 copy_rvec(comm->vbuf.v[buf_pos++],
5016 state->sd_X[home_pos_at+i]);
5021 for (i = 0; i < nrcg; i++)
5023 copy_rvec(comm->vbuf.v[buf_pos++],
5024 state->cg_p[home_pos_at+i]);
5028 home_pos_at += nrcg;
5032 /* Reallocate the buffers if necessary */
5033 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5035 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5036 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5038 nvr = ncg[mc] + nat[mc]*nvec;
5039 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5041 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5042 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5044 /* Copy from the receive to the send buffers */
5045 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5046 comm->buf_int + cg*DD_CGIBS,
5047 DD_CGIBS*sizeof(int));
5048 memcpy(comm->cgcm_state[mc][nvr],
5049 comm->vbuf.v[buf_pos],
5050 (1+nrcg*nvec)*sizeof(rvec));
5051 buf_pos += 1 + nrcg*nvec;
5058 /* With sorting (!bCompact) the indices are now only partially up to date
5059 * and ncg_home and nat_home are not the real count, since there are
5060 * "holes" in the arrays for the charge groups that moved to neighbors.
5062 if (fr->cutoff_scheme == ecutsVERLET)
5064 moved = get_moved(comm, home_pos_cg);
5066 for (i = dd->ncg_home; i < home_pos_cg; i++)
5071 dd->ncg_home = home_pos_cg;
5072 dd->nat_home = home_pos_at;
5077 "Finished repartitioning: cgs moved out %d, new home %d\n",
5078 *ncg_moved, dd->ncg_home-*ncg_moved);
5083 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5085 dd->comm->cycl[ddCycl] += cycles;
5086 dd->comm->cycl_n[ddCycl]++;
5087 if (cycles > dd->comm->cycl_max[ddCycl])
5089 dd->comm->cycl_max[ddCycl] = cycles;
5093 static double force_flop_count(t_nrnb *nrnb)
5100 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5102 /* To get closer to the real timings, we half the count
5103 * for the normal loops and again half it for water loops.
5106 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5108 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5112 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5115 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5118 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5120 sum += nrnb->n[i]*cost_nrnb(i);
5123 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5125 sum += nrnb->n[i]*cost_nrnb(i);
5131 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5133 if (dd->comm->eFlop)
5135 dd->comm->flop -= force_flop_count(nrnb);
5138 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5140 if (dd->comm->eFlop)
5142 dd->comm->flop += force_flop_count(nrnb);
5147 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5151 for (i = 0; i < ddCyclNr; i++)
5153 dd->comm->cycl[i] = 0;
5154 dd->comm->cycl_n[i] = 0;
5155 dd->comm->cycl_max[i] = 0;
5158 dd->comm->flop_n = 0;
5161 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5163 gmx_domdec_comm_t *comm;
5164 gmx_domdec_load_t *load;
5165 gmx_domdec_root_t *root = NULL;
5166 int d, dim, cid, i, pos;
5167 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5172 fprintf(debug, "get_load_distribution start\n");
5175 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5179 bSepPME = (dd->pme_nodeid >= 0);
5181 for (d = dd->ndim-1; d >= 0; d--)
5184 /* Check if we participate in the communication in this dimension */
5185 if (d == dd->ndim-1 ||
5186 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5188 load = &comm->load[d];
5191 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5194 if (d == dd->ndim-1)
5196 sbuf[pos++] = dd_force_load(comm);
5197 sbuf[pos++] = sbuf[0];
5200 sbuf[pos++] = sbuf[0];
5201 sbuf[pos++] = cell_frac;
5204 sbuf[pos++] = comm->cell_f_max0[d];
5205 sbuf[pos++] = comm->cell_f_min1[d];
5210 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5211 sbuf[pos++] = comm->cycl[ddCyclPME];
5216 sbuf[pos++] = comm->load[d+1].sum;
5217 sbuf[pos++] = comm->load[d+1].max;
5220 sbuf[pos++] = comm->load[d+1].sum_m;
5221 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5222 sbuf[pos++] = comm->load[d+1].flags;
5225 sbuf[pos++] = comm->cell_f_max0[d];
5226 sbuf[pos++] = comm->cell_f_min1[d];
5231 sbuf[pos++] = comm->load[d+1].mdf;
5232 sbuf[pos++] = comm->load[d+1].pme;
5236 /* Communicate a row in DD direction d.
5237 * The communicators are setup such that the root always has rank 0.
5240 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5241 load->load, load->nload*sizeof(float), MPI_BYTE,
5242 0, comm->mpi_comm_load[d]);
5244 if (dd->ci[dim] == dd->master_ci[dim])
5246 /* We are the root, process this row */
5247 if (comm->bDynLoadBal)
5249 root = comm->root[d];
5259 for (i = 0; i < dd->nc[dim]; i++)
5261 load->sum += load->load[pos++];
5262 load->max = max(load->max, load->load[pos]);
5268 /* This direction could not be load balanced properly,
5269 * therefore we need to use the maximum iso the average load.
5271 load->sum_m = max(load->sum_m, load->load[pos]);
5275 load->sum_m += load->load[pos];
5278 load->cvol_min = min(load->cvol_min, load->load[pos]);
5282 load->flags = (int)(load->load[pos++] + 0.5);
5286 root->cell_f_max0[i] = load->load[pos++];
5287 root->cell_f_min1[i] = load->load[pos++];
5292 load->mdf = max(load->mdf, load->load[pos]);
5294 load->pme = max(load->pme, load->load[pos]);
5298 if (comm->bDynLoadBal && root->bLimited)
5300 load->sum_m *= dd->nc[dim];
5301 load->flags |= (1<<d);
5309 comm->nload += dd_load_count(comm);
5310 comm->load_step += comm->cycl[ddCyclStep];
5311 comm->load_sum += comm->load[0].sum;
5312 comm->load_max += comm->load[0].max;
5313 if (comm->bDynLoadBal)
5315 for (d = 0; d < dd->ndim; d++)
5317 if (comm->load[0].flags & (1<<d))
5319 comm->load_lim[d]++;
5325 comm->load_mdf += comm->load[0].mdf;
5326 comm->load_pme += comm->load[0].pme;
5330 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5334 fprintf(debug, "get_load_distribution finished\n");
5338 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5340 /* Return the relative performance loss on the total run time
5341 * due to the force calculation load imbalance.
5343 if (dd->comm->nload > 0)
5346 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5347 (dd->comm->load_step*dd->nnodes);
5355 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5358 int npp, npme, nnodes, d, limp;
5359 float imbal, pme_f_ratio, lossf, lossp = 0;
5361 gmx_domdec_comm_t *comm;
5364 if (DDMASTER(dd) && comm->nload > 0)
5367 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5368 nnodes = npp + npme;
5369 imbal = comm->load_max*npp/comm->load_sum - 1;
5370 lossf = dd_force_imb_perf_loss(dd);
5371 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5372 fprintf(fplog, "%s", buf);
5373 fprintf(stderr, "\n");
5374 fprintf(stderr, "%s", buf);
5375 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5376 fprintf(fplog, "%s", buf);
5377 fprintf(stderr, "%s", buf);
5379 if (comm->bDynLoadBal)
5381 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5382 for (d = 0; d < dd->ndim; d++)
5384 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5385 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5391 sprintf(buf+strlen(buf), "\n");
5392 fprintf(fplog, "%s", buf);
5393 fprintf(stderr, "%s", buf);
5397 pme_f_ratio = comm->load_pme/comm->load_mdf;
5398 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5401 lossp *= (float)npme/(float)nnodes;
5405 lossp *= (float)npp/(float)nnodes;
5407 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5408 fprintf(fplog, "%s", buf);
5409 fprintf(stderr, "%s", buf);
5410 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5411 fprintf(fplog, "%s", buf);
5412 fprintf(stderr, "%s", buf);
5414 fprintf(fplog, "\n");
5415 fprintf(stderr, "\n");
5417 if (lossf >= DD_PERF_LOSS)
5420 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5421 " in the domain decomposition.\n", lossf*100);
5422 if (!comm->bDynLoadBal)
5424 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5428 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5430 fprintf(fplog, "%s\n", buf);
5431 fprintf(stderr, "%s\n", buf);
5433 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5436 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5437 " had %s work to do than the PP nodes.\n"
5438 " You might want to %s the number of PME nodes\n"
5439 " or %s the cut-off and the grid spacing.\n",
5441 (lossp < 0) ? "less" : "more",
5442 (lossp < 0) ? "decrease" : "increase",
5443 (lossp < 0) ? "decrease" : "increase");
5444 fprintf(fplog, "%s\n", buf);
5445 fprintf(stderr, "%s\n", buf);
5450 static float dd_vol_min(gmx_domdec_t *dd)
5452 return dd->comm->load[0].cvol_min*dd->nnodes;
5455 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5457 return dd->comm->load[0].flags;
5460 static float dd_f_imbal(gmx_domdec_t *dd)
5462 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5465 float dd_pme_f_ratio(gmx_domdec_t *dd)
5467 if (dd->comm->cycl_n[ddCyclPME] > 0)
5469 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5477 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5482 flags = dd_load_flags(dd);
5486 "DD load balancing is limited by minimum cell size in dimension");
5487 for (d = 0; d < dd->ndim; d++)
5491 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5494 fprintf(fplog, "\n");
5496 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5497 if (dd->comm->bDynLoadBal)
5499 fprintf(fplog, " vol min/aver %5.3f%c",
5500 dd_vol_min(dd), flags ? '!' : ' ');
5502 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5503 if (dd->comm->cycl_n[ddCyclPME])
5505 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5507 fprintf(fplog, "\n\n");
5510 static void dd_print_load_verbose(gmx_domdec_t *dd)
5512 if (dd->comm->bDynLoadBal)
5514 fprintf(stderr, "vol %4.2f%c ",
5515 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5517 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5518 if (dd->comm->cycl_n[ddCyclPME])
5520 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5525 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5530 gmx_domdec_root_t *root;
5531 gmx_bool bPartOfGroup = FALSE;
5533 dim = dd->dim[dim_ind];
5534 copy_ivec(loc, loc_c);
5535 for (i = 0; i < dd->nc[dim]; i++)
5538 rank = dd_index(dd->nc, loc_c);
5539 if (rank == dd->rank)
5541 /* This process is part of the group */
5542 bPartOfGroup = TRUE;
5545 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5549 dd->comm->mpi_comm_load[dim_ind] = c_row;
5550 if (dd->comm->eDLB != edlbNO)
5552 if (dd->ci[dim] == dd->master_ci[dim])
5554 /* This is the root process of this row */
5555 snew(dd->comm->root[dim_ind], 1);
5556 root = dd->comm->root[dim_ind];
5557 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5558 snew(root->old_cell_f, dd->nc[dim]+1);
5559 snew(root->bCellMin, dd->nc[dim]);
5562 snew(root->cell_f_max0, dd->nc[dim]);
5563 snew(root->cell_f_min1, dd->nc[dim]);
5564 snew(root->bound_min, dd->nc[dim]);
5565 snew(root->bound_max, dd->nc[dim]);
5567 snew(root->buf_ncd, dd->nc[dim]);
5571 /* This is not a root process, we only need to receive cell_f */
5572 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5575 if (dd->ci[dim] == dd->master_ci[dim])
5577 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5583 static void make_load_communicators(gmx_domdec_t *dd)
5586 int dim0, dim1, i, j;
5591 fprintf(debug, "Making load communicators\n");
5594 snew(dd->comm->load, dd->ndim);
5595 snew(dd->comm->mpi_comm_load, dd->ndim);
5598 make_load_communicator(dd, 0, loc);
5602 for (i = 0; i < dd->nc[dim0]; i++)
5605 make_load_communicator(dd, 1, loc);
5611 for (i = 0; i < dd->nc[dim0]; i++)
5615 for (j = 0; j < dd->nc[dim1]; j++)
5618 make_load_communicator(dd, 2, loc);
5625 fprintf(debug, "Finished making load communicators\n");
5630 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5633 int d, dim, i, j, m;
5636 ivec dd_zp[DD_MAXIZONE];
5637 gmx_domdec_zones_t *zones;
5638 gmx_domdec_ns_ranges_t *izone;
5640 for (d = 0; d < dd->ndim; d++)
5643 copy_ivec(dd->ci, tmp);
5644 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5645 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5646 copy_ivec(dd->ci, tmp);
5647 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5648 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5651 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5654 dd->neighbor[d][1]);
5660 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5662 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5663 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5670 for (i = 0; i < nzonep; i++)
5672 copy_ivec(dd_zp3[i], dd_zp[i]);
5678 for (i = 0; i < nzonep; i++)
5680 copy_ivec(dd_zp2[i], dd_zp[i]);
5686 for (i = 0; i < nzonep; i++)
5688 copy_ivec(dd_zp1[i], dd_zp[i]);
5692 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5697 zones = &dd->comm->zones;
5699 for (i = 0; i < nzone; i++)
5702 clear_ivec(zones->shift[i]);
5703 for (d = 0; d < dd->ndim; d++)
5705 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5710 for (i = 0; i < nzone; i++)
5712 for (d = 0; d < DIM; d++)
5714 s[d] = dd->ci[d] - zones->shift[i][d];
5719 else if (s[d] >= dd->nc[d])
5725 zones->nizone = nzonep;
5726 for (i = 0; i < zones->nizone; i++)
5728 if (dd_zp[i][0] != i)
5730 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5732 izone = &zones->izone[i];
5733 izone->j0 = dd_zp[i][1];
5734 izone->j1 = dd_zp[i][2];
5735 for (dim = 0; dim < DIM; dim++)
5737 if (dd->nc[dim] == 1)
5739 /* All shifts should be allowed */
5740 izone->shift0[dim] = -1;
5741 izone->shift1[dim] = 1;
5746 izone->shift0[d] = 0;
5747 izone->shift1[d] = 0;
5748 for(j=izone->j0; j<izone->j1; j++) {
5749 if (dd->shift[j][d] > dd->shift[i][d])
5750 izone->shift0[d] = -1;
5751 if (dd->shift[j][d] < dd->shift[i][d])
5752 izone->shift1[d] = 1;
5758 /* Assume the shift are not more than 1 cell */
5759 izone->shift0[dim] = 1;
5760 izone->shift1[dim] = -1;
5761 for (j = izone->j0; j < izone->j1; j++)
5763 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5764 if (shift_diff < izone->shift0[dim])
5766 izone->shift0[dim] = shift_diff;
5768 if (shift_diff > izone->shift1[dim])
5770 izone->shift1[dim] = shift_diff;
5777 if (dd->comm->eDLB != edlbNO)
5779 snew(dd->comm->root, dd->ndim);
5782 if (dd->comm->bRecordLoad)
5784 make_load_communicators(dd);
5788 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5791 gmx_domdec_comm_t *comm;
5802 if (comm->bCartesianPP)
5804 /* Set up cartesian communication for the particle-particle part */
5807 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5808 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5811 for (i = 0; i < DIM; i++)
5815 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5817 /* We overwrite the old communicator with the new cartesian one */
5818 cr->mpi_comm_mygroup = comm_cart;
5821 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5822 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5824 if (comm->bCartesianPP_PME)
5826 /* Since we want to use the original cartesian setup for sim,
5827 * and not the one after split, we need to make an index.
5829 snew(comm->ddindex2ddnodeid, dd->nnodes);
5830 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5831 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5832 /* Get the rank of the DD master,
5833 * above we made sure that the master node is a PP node.
5843 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5845 else if (comm->bCartesianPP)
5847 if (cr->npmenodes == 0)
5849 /* The PP communicator is also
5850 * the communicator for this simulation
5852 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5854 cr->nodeid = dd->rank;
5856 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5858 /* We need to make an index to go from the coordinates
5859 * to the nodeid of this simulation.
5861 snew(comm->ddindex2simnodeid, dd->nnodes);
5862 snew(buf, dd->nnodes);
5863 if (cr->duty & DUTY_PP)
5865 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5867 /* Communicate the ddindex to simulation nodeid index */
5868 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5869 cr->mpi_comm_mysim);
5872 /* Determine the master coordinates and rank.
5873 * The DD master should be the same node as the master of this sim.
5875 for (i = 0; i < dd->nnodes; i++)
5877 if (comm->ddindex2simnodeid[i] == 0)
5879 ddindex2xyz(dd->nc, i, dd->master_ci);
5880 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5885 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5890 /* No Cartesian communicators */
5891 /* We use the rank in dd->comm->all as DD index */
5892 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5893 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5895 clear_ivec(dd->master_ci);
5902 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5903 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5908 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5909 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5913 static void receive_ddindex2simnodeid(t_commrec *cr)
5917 gmx_domdec_comm_t *comm;
5924 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5926 snew(comm->ddindex2simnodeid, dd->nnodes);
5927 snew(buf, dd->nnodes);
5928 if (cr->duty & DUTY_PP)
5930 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5933 /* Communicate the ddindex to simulation nodeid index */
5934 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5935 cr->mpi_comm_mysim);
5942 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5943 int ncg, int natoms)
5945 gmx_domdec_master_t *ma;
5950 snew(ma->ncg, dd->nnodes);
5951 snew(ma->index, dd->nnodes+1);
5953 snew(ma->nat, dd->nnodes);
5954 snew(ma->ibuf, dd->nnodes*2);
5955 snew(ma->cell_x, DIM);
5956 for (i = 0; i < DIM; i++)
5958 snew(ma->cell_x[i], dd->nc[i]+1);
5961 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5967 snew(ma->vbuf, natoms);
5973 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5977 gmx_domdec_comm_t *comm;
5988 if (comm->bCartesianPP)
5990 for (i = 1; i < DIM; i++)
5992 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5994 if (bDiv[YY] || bDiv[ZZ])
5996 comm->bCartesianPP_PME = TRUE;
5997 /* If we have 2D PME decomposition, which is always in x+y,
5998 * we stack the PME only nodes in z.
5999 * Otherwise we choose the direction that provides the thinnest slab
6000 * of PME only nodes as this will have the least effect
6001 * on the PP communication.
6002 * But for the PME communication the opposite might be better.
6004 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6006 dd->nc[YY] > dd->nc[ZZ]))
6008 comm->cartpmedim = ZZ;
6012 comm->cartpmedim = YY;
6014 comm->ntot[comm->cartpmedim]
6015 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6019 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]);
6021 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6026 if (comm->bCartesianPP_PME)
6030 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]);
6033 for (i = 0; i < DIM; i++)
6037 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6040 MPI_Comm_rank(comm_cart, &rank);
6041 if (MASTERNODE(cr) && rank != 0)
6043 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6046 /* With this assigment we loose the link to the original communicator
6047 * which will usually be MPI_COMM_WORLD, unless have multisim.
6049 cr->mpi_comm_mysim = comm_cart;
6050 cr->sim_nodeid = rank;
6052 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6056 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6057 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6060 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6064 if (cr->npmenodes == 0 ||
6065 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6067 cr->duty = DUTY_PME;
6070 /* Split the sim communicator into PP and PME only nodes */
6071 MPI_Comm_split(cr->mpi_comm_mysim,
6073 dd_index(comm->ntot, dd->ci),
6074 &cr->mpi_comm_mygroup);
6078 switch (dd_node_order)
6083 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6086 case ddnoINTERLEAVE:
6087 /* Interleave the PP-only and PME-only nodes,
6088 * as on clusters with dual-core machines this will double
6089 * the communication bandwidth of the PME processes
6090 * and thus speed up the PP <-> PME and inter PME communication.
6094 fprintf(fplog, "Interleaving PP and PME nodes\n");
6096 comm->pmenodes = dd_pmenodes(cr);
6101 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6104 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6106 cr->duty = DUTY_PME;
6113 /* Split the sim communicator into PP and PME only nodes */
6114 MPI_Comm_split(cr->mpi_comm_mysim,
6117 &cr->mpi_comm_mygroup);
6118 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6124 fprintf(fplog, "This is a %s only node\n\n",
6125 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6129 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6132 gmx_domdec_comm_t *comm;
6138 copy_ivec(dd->nc, comm->ntot);
6140 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6141 comm->bCartesianPP_PME = FALSE;
6143 /* Reorder the nodes by default. This might change the MPI ranks.
6144 * Real reordering is only supported on very few architectures,
6145 * Blue Gene is one of them.
6147 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6149 if (cr->npmenodes > 0)
6151 /* Split the communicator into a PP and PME part */
6152 split_communicator(fplog, cr, dd_node_order, CartReorder);
6153 if (comm->bCartesianPP_PME)
6155 /* We (possibly) reordered the nodes in split_communicator,
6156 * so it is no longer required in make_pp_communicator.
6158 CartReorder = FALSE;
6163 /* All nodes do PP and PME */
6165 /* We do not require separate communicators */
6166 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6170 if (cr->duty & DUTY_PP)
6172 /* Copy or make a new PP communicator */
6173 make_pp_communicator(fplog, cr, CartReorder);
6177 receive_ddindex2simnodeid(cr);
6180 if (!(cr->duty & DUTY_PME))
6182 /* Set up the commnuication to our PME node */
6183 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6184 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6187 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6188 dd->pme_nodeid, dd->pme_receive_vir_ener);
6193 dd->pme_nodeid = -1;
6198 dd->ma = init_gmx_domdec_master_t(dd,
6200 comm->cgs_gl.index[comm->cgs_gl.nr]);
6204 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6206 real *slb_frac, tot;
6211 if (nc > 1 && size_string != NULL)
6215 fprintf(fplog, "Using static load balancing for the %s direction\n",
6220 for (i = 0; i < nc; i++)
6223 sscanf(size_string, "%lf%n", &dbl, &n);
6226 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6235 fprintf(fplog, "Relative cell sizes:");
6237 for (i = 0; i < nc; i++)
6242 fprintf(fplog, " %5.3f", slb_frac[i]);
6247 fprintf(fplog, "\n");
6254 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6257 gmx_mtop_ilistloop_t iloop;
6261 iloop = gmx_mtop_ilistloop_init(mtop);
6262 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6264 for (ftype = 0; ftype < F_NRE; ftype++)
6266 if ((interaction_function[ftype].flags & IF_BOND) &&
6269 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6277 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6283 val = getenv(env_var);
6286 if (sscanf(val, "%d", &nst) <= 0)
6292 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6300 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6304 fprintf(stderr, "\n%s\n", warn_string);
6308 fprintf(fplog, "\n%s\n", warn_string);
6312 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6313 t_inputrec *ir, FILE *fplog)
6315 if (ir->ePBC == epbcSCREW &&
6316 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6318 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6321 if (ir->ns_type == ensSIMPLE)
6323 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6326 if (ir->nstlist == 0)
6328 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6331 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6333 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6337 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6342 r = ddbox->box_size[XX];
6343 for (di = 0; di < dd->ndim; di++)
6346 /* Check using the initial average cell size */
6347 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6353 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6354 const char *dlb_opt, gmx_bool bRecordLoad,
6355 unsigned long Flags, t_inputrec *ir)
6363 case 'a': eDLB = edlbAUTO; break;
6364 case 'n': eDLB = edlbNO; break;
6365 case 'y': eDLB = edlbYES; break;
6366 default: gmx_incons("Unknown dlb_opt");
6369 if (Flags & MD_RERUN)
6374 if (!EI_DYNAMICS(ir->eI))
6376 if (eDLB == edlbYES)
6378 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6379 dd_warning(cr, fplog, buf);
6387 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6392 if (Flags & MD_REPRODUCIBLE)
6399 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6403 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6406 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6414 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6419 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6421 /* Decomposition order z,y,x */
6424 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6426 for (dim = DIM-1; dim >= 0; dim--)
6428 if (dd->nc[dim] > 1)
6430 dd->dim[dd->ndim++] = dim;
6436 /* Decomposition order x,y,z */
6437 for (dim = 0; dim < DIM; dim++)
6439 if (dd->nc[dim] > 1)
6441 dd->dim[dd->ndim++] = dim;
6447 static gmx_domdec_comm_t *init_dd_comm()
6449 gmx_domdec_comm_t *comm;
6453 snew(comm->cggl_flag, DIM*2);
6454 snew(comm->cgcm_state, DIM*2);
6455 for (i = 0; i < DIM*2; i++)
6457 comm->cggl_flag_nalloc[i] = 0;
6458 comm->cgcm_state_nalloc[i] = 0;
6461 comm->nalloc_int = 0;
6462 comm->buf_int = NULL;
6464 vec_rvec_init(&comm->vbuf);
6466 comm->n_load_have = 0;
6467 comm->n_load_collect = 0;
6469 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6471 comm->sum_nat[i] = 0;
6475 comm->load_step = 0;
6478 clear_ivec(comm->load_lim);
6485 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6486 unsigned long Flags,
6488 real comm_distance_min, real rconstr,
6489 const char *dlb_opt, real dlb_scale,
6490 const char *sizex, const char *sizey, const char *sizez,
6491 gmx_mtop_t *mtop, t_inputrec *ir,
6492 matrix box, rvec *x,
6494 int *npme_x, int *npme_y)
6497 gmx_domdec_comm_t *comm;
6500 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6507 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6512 dd->comm = init_dd_comm();
6514 snew(comm->cggl_flag, DIM*2);
6515 snew(comm->cgcm_state, DIM*2);
6517 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6518 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6520 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6521 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6522 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6523 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6524 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6525 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6526 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6527 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6529 dd->pme_recv_f_alloc = 0;
6530 dd->pme_recv_f_buf = NULL;
6532 if (dd->bSendRecv2 && fplog)
6534 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");
6540 fprintf(fplog, "Will load balance based on FLOP count\n");
6542 if (comm->eFlop > 1)
6544 srand(1+cr->nodeid);
6546 comm->bRecordLoad = TRUE;
6550 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6554 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6556 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6559 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6561 dd->bGridJump = comm->bDynLoadBal;
6562 comm->bPMELoadBalDLBLimits = FALSE;
6564 if (comm->nstSortCG)
6568 if (comm->nstSortCG == 1)
6570 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6574 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6578 snew(comm->sort, 1);
6584 fprintf(fplog, "Will not sort the charge groups\n");
6588 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6590 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6591 if (comm->bInterCGBondeds)
6593 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6597 comm->bInterCGMultiBody = FALSE;
6600 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6601 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6603 if (ir->rlistlong == 0)
6605 /* Set the cut-off to some very large value,
6606 * so we don't need if statements everywhere in the code.
6607 * We use sqrt, since the cut-off is squared in some places.
6609 comm->cutoff = GMX_CUTOFF_INF;
6613 comm->cutoff = ir->rlistlong;
6615 comm->cutoff_mbody = 0;
6617 comm->cellsize_limit = 0;
6618 comm->bBondComm = FALSE;
6620 if (comm->bInterCGBondeds)
6622 if (comm_distance_min > 0)
6624 comm->cutoff_mbody = comm_distance_min;
6625 if (Flags & MD_DDBONDCOMM)
6627 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6631 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6633 r_bonded_limit = comm->cutoff_mbody;
6635 else if (ir->bPeriodicMols)
6637 /* Can not easily determine the required cut-off */
6638 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");
6639 comm->cutoff_mbody = comm->cutoff/2;
6640 r_bonded_limit = comm->cutoff_mbody;
6646 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6647 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6649 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6650 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6652 /* We use an initial margin of 10% for the minimum cell size,
6653 * except when we are just below the non-bonded cut-off.
6655 if (Flags & MD_DDBONDCOMM)
6657 if (max(r_2b, r_mb) > comm->cutoff)
6659 r_bonded = max(r_2b, r_mb);
6660 r_bonded_limit = 1.1*r_bonded;
6661 comm->bBondComm = TRUE;
6666 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6668 /* We determine cutoff_mbody later */
6672 /* No special bonded communication,
6673 * simply increase the DD cut-off.
6675 r_bonded_limit = 1.1*max(r_2b, r_mb);
6676 comm->cutoff_mbody = r_bonded_limit;
6677 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6680 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6684 "Minimum cell size due to bonded interactions: %.3f nm\n",
6685 comm->cellsize_limit);
6689 if (dd->bInterCGcons && rconstr <= 0)
6691 /* There is a cell size limit due to the constraints (P-LINCS) */
6692 rconstr = constr_r_max(fplog, mtop, ir);
6696 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6698 if (rconstr > comm->cellsize_limit)
6700 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6704 else if (rconstr > 0 && fplog)
6706 /* Here we do not check for dd->bInterCGcons,
6707 * because one can also set a cell size limit for virtual sites only
6708 * and at this point we don't know yet if there are intercg v-sites.
6711 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6714 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6716 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6720 copy_ivec(nc, dd->nc);
6721 set_dd_dim(fplog, dd);
6722 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6724 if (cr->npmenodes == -1)
6728 acs = average_cellsize_min(dd, ddbox);
6729 if (acs < comm->cellsize_limit)
6733 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6735 gmx_fatal_collective(FARGS, cr, NULL,
6736 "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",
6737 acs, comm->cellsize_limit);
6742 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6744 /* We need to choose the optimal DD grid and possibly PME nodes */
6745 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6746 comm->eDLB != edlbNO, dlb_scale,
6747 comm->cellsize_limit, comm->cutoff,
6748 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6750 if (dd->nc[XX] == 0)
6752 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6753 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6754 !bC ? "-rdd" : "-rcon",
6755 comm->eDLB != edlbNO ? " or -dds" : "",
6756 bC ? " or your LINCS settings" : "");
6758 gmx_fatal_collective(FARGS, cr, NULL,
6759 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6761 "Look in the log file for details on the domain decomposition",
6762 cr->nnodes-cr->npmenodes, limit, buf);
6764 set_dd_dim(fplog, dd);
6770 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6771 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6774 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6775 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6777 gmx_fatal_collective(FARGS, cr, NULL,
6778 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6779 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6781 if (cr->npmenodes > dd->nnodes)
6783 gmx_fatal_collective(FARGS, cr, NULL,
6784 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6786 if (cr->npmenodes > 0)
6788 comm->npmenodes = cr->npmenodes;
6792 comm->npmenodes = dd->nnodes;
6795 if (EEL_PME(ir->coulombtype))
6797 /* The following choices should match those
6798 * in comm_cost_est in domdec_setup.c.
6799 * Note that here the checks have to take into account
6800 * that the decomposition might occur in a different order than xyz
6801 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6802 * in which case they will not match those in comm_cost_est,
6803 * but since that is mainly for testing purposes that's fine.
6805 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6806 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6807 getenv("GMX_PMEONEDD") == NULL)
6809 comm->npmedecompdim = 2;
6810 comm->npmenodes_x = dd->nc[XX];
6811 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6815 /* In case nc is 1 in both x and y we could still choose to
6816 * decompose pme in y instead of x, but we use x for simplicity.
6818 comm->npmedecompdim = 1;
6819 if (dd->dim[0] == YY)
6821 comm->npmenodes_x = 1;
6822 comm->npmenodes_y = comm->npmenodes;
6826 comm->npmenodes_x = comm->npmenodes;
6827 comm->npmenodes_y = 1;
6832 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6833 comm->npmenodes_x, comm->npmenodes_y, 1);
6838 comm->npmedecompdim = 0;
6839 comm->npmenodes_x = 0;
6840 comm->npmenodes_y = 0;
6843 /* Technically we don't need both of these,
6844 * but it simplifies code not having to recalculate it.
6846 *npme_x = comm->npmenodes_x;
6847 *npme_y = comm->npmenodes_y;
6849 snew(comm->slb_frac, DIM);
6850 if (comm->eDLB == edlbNO)
6852 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6853 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6854 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6857 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6859 if (comm->bBondComm || comm->eDLB != edlbNO)
6861 /* Set the bonded communication distance to halfway
6862 * the minimum and the maximum,
6863 * since the extra communication cost is nearly zero.
6865 acs = average_cellsize_min(dd, ddbox);
6866 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6867 if (comm->eDLB != edlbNO)
6869 /* Check if this does not limit the scaling */
6870 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6872 if (!comm->bBondComm)
6874 /* Without bBondComm do not go beyond the n.b. cut-off */
6875 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6876 if (comm->cellsize_limit >= comm->cutoff)
6878 /* We don't loose a lot of efficieny
6879 * when increasing it to the n.b. cut-off.
6880 * It can even be slightly faster, because we need
6881 * less checks for the communication setup.
6883 comm->cutoff_mbody = comm->cutoff;
6886 /* Check if we did not end up below our original limit */
6887 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6889 if (comm->cutoff_mbody > comm->cellsize_limit)
6891 comm->cellsize_limit = comm->cutoff_mbody;
6894 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6899 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6900 "cellsize limit %f\n",
6901 comm->bBondComm, comm->cellsize_limit);
6906 check_dd_restrictions(cr, dd, ir, fplog);
6909 comm->partition_step = INT_MIN;
6912 clear_dd_cycle_counts(dd);
6917 static void set_dlb_limits(gmx_domdec_t *dd)
6922 for (d = 0; d < dd->ndim; d++)
6924 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6925 dd->comm->cellsize_min[dd->dim[d]] =
6926 dd->comm->cellsize_min_dlb[dd->dim[d]];
6931 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6934 gmx_domdec_comm_t *comm;
6944 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);
6947 cellsize_min = comm->cellsize_min[dd->dim[0]];
6948 for (d = 1; d < dd->ndim; d++)
6950 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6953 if (cellsize_min < comm->cellsize_limit*1.05)
6955 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");
6957 /* Change DLB from "auto" to "no". */
6958 comm->eDLB = edlbNO;
6963 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6964 comm->bDynLoadBal = TRUE;
6965 dd->bGridJump = TRUE;
6969 /* We can set the required cell size info here,
6970 * so we do not need to communicate this.
6971 * The grid is completely uniform.
6973 for (d = 0; d < dd->ndim; d++)
6977 comm->load[d].sum_m = comm->load[d].sum;
6979 nc = dd->nc[dd->dim[d]];
6980 for (i = 0; i < nc; i++)
6982 comm->root[d]->cell_f[i] = i/(real)nc;
6985 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6986 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6989 comm->root[d]->cell_f[nc] = 1.0;
6994 static char *init_bLocalCG(gmx_mtop_t *mtop)
6999 ncg = ncg_mtop(mtop);
7000 snew(bLocalCG, ncg);
7001 for (cg = 0; cg < ncg; cg++)
7003 bLocalCG[cg] = FALSE;
7009 void dd_init_bondeds(FILE *fplog,
7010 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7011 gmx_vsite_t *vsite, gmx_constr_t constr,
7012 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7014 gmx_domdec_comm_t *comm;
7018 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7022 if (comm->bBondComm)
7024 /* Communicate atoms beyond the cut-off for bonded interactions */
7027 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7029 comm->bLocalCG = init_bLocalCG(mtop);
7033 /* Only communicate atoms based on cut-off */
7034 comm->cglink = NULL;
7035 comm->bLocalCG = NULL;
7039 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7041 gmx_bool bDynLoadBal, real dlb_scale,
7044 gmx_domdec_comm_t *comm;
7059 fprintf(fplog, "The maximum number of communication pulses is:");
7060 for (d = 0; d < dd->ndim; d++)
7062 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7064 fprintf(fplog, "\n");
7065 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7066 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7067 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7068 for (d = 0; d < DIM; d++)
7072 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7079 comm->cellsize_min_dlb[d]/
7080 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7082 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7085 fprintf(fplog, "\n");
7089 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7090 fprintf(fplog, "The initial number of communication pulses is:");
7091 for (d = 0; d < dd->ndim; d++)
7093 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7095 fprintf(fplog, "\n");
7096 fprintf(fplog, "The initial domain decomposition cell size is:");
7097 for (d = 0; d < DIM; d++)
7101 fprintf(fplog, " %c %.2f nm",
7102 dim2char(d), dd->comm->cellsize_min[d]);
7105 fprintf(fplog, "\n\n");
7108 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7110 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7111 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7112 "non-bonded interactions", "", comm->cutoff);
7116 limit = dd->comm->cellsize_limit;
7120 if (dynamic_dd_box(ddbox, ir))
7122 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7124 limit = dd->comm->cellsize_min[XX];
7125 for (d = 1; d < DIM; d++)
7127 limit = min(limit, dd->comm->cellsize_min[d]);
7131 if (comm->bInterCGBondeds)
7133 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7134 "two-body bonded interactions", "(-rdd)",
7135 max(comm->cutoff, comm->cutoff_mbody));
7136 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7137 "multi-body bonded interactions", "(-rdd)",
7138 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7142 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7143 "virtual site constructions", "(-rcon)", limit);
7145 if (dd->constraint_comm)
7147 sprintf(buf, "atoms separated by up to %d constraints",
7149 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7150 buf, "(-rcon)", limit);
7152 fprintf(fplog, "\n");
7158 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7160 const t_inputrec *ir,
7161 const gmx_ddbox_t *ddbox)
7163 gmx_domdec_comm_t *comm;
7164 int d, dim, npulse, npulse_d_max, npulse_d;
7169 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7171 /* Determine the maximum number of comm. pulses in one dimension */
7173 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7175 /* Determine the maximum required number of grid pulses */
7176 if (comm->cellsize_limit >= comm->cutoff)
7178 /* Only a single pulse is required */
7181 else if (!bNoCutOff && comm->cellsize_limit > 0)
7183 /* We round down slightly here to avoid overhead due to the latency
7184 * of extra communication calls when the cut-off
7185 * would be only slightly longer than the cell size.
7186 * Later cellsize_limit is redetermined,
7187 * so we can not miss interactions due to this rounding.
7189 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7193 /* There is no cell size limit */
7194 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7197 if (!bNoCutOff && npulse > 1)
7199 /* See if we can do with less pulses, based on dlb_scale */
7201 for (d = 0; d < dd->ndim; d++)
7204 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7205 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7206 npulse_d_max = max(npulse_d_max, npulse_d);
7208 npulse = min(npulse, npulse_d_max);
7211 /* This env var can override npulse */
7212 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7219 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7220 for (d = 0; d < dd->ndim; d++)
7222 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7223 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7224 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7225 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7226 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7228 comm->bVacDLBNoLimit = FALSE;
7232 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7233 if (!comm->bVacDLBNoLimit)
7235 comm->cellsize_limit = max(comm->cellsize_limit,
7236 comm->cutoff/comm->maxpulse);
7238 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7239 /* Set the minimum cell size for each DD dimension */
7240 for (d = 0; d < dd->ndim; d++)
7242 if (comm->bVacDLBNoLimit ||
7243 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7245 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7249 comm->cellsize_min_dlb[dd->dim[d]] =
7250 comm->cutoff/comm->cd[d].np_dlb;
7253 if (comm->cutoff_mbody <= 0)
7255 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7257 if (comm->bDynLoadBal)
7263 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7265 /* If each molecule is a single charge group
7266 * or we use domain decomposition for each periodic dimension,
7267 * we do not need to take pbc into account for the bonded interactions.
7269 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7272 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7275 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7276 t_inputrec *ir, t_forcerec *fr,
7279 gmx_domdec_comm_t *comm;
7285 /* Initialize the thread data.
7286 * This can not be done in init_domain_decomposition,
7287 * as the numbers of threads is determined later.
7289 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7292 snew(comm->dth, comm->nth);
7295 if (EEL_PME(ir->coulombtype))
7297 init_ddpme(dd, &comm->ddpme[0], 0);
7298 if (comm->npmedecompdim >= 2)
7300 init_ddpme(dd, &comm->ddpme[1], 1);
7305 comm->npmenodes = 0;
7306 if (dd->pme_nodeid >= 0)
7308 gmx_fatal_collective(FARGS, NULL, dd,
7309 "Can not have separate PME nodes without PME electrostatics");
7315 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7317 if (comm->eDLB != edlbNO)
7319 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7322 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7323 if (comm->eDLB == edlbAUTO)
7327 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7329 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7332 if (ir->ePBC == epbcNONE)
7334 vol_frac = 1 - 1/(double)dd->nnodes;
7339 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7343 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7345 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7347 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7350 static gmx_bool test_dd_cutoff(t_commrec *cr,
7351 t_state *state, t_inputrec *ir,
7362 set_ddbox(dd, FALSE, cr, ir, state->box,
7363 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7367 for (d = 0; d < dd->ndim; d++)
7371 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7372 if (dynamic_dd_box(&ddbox, ir))
7374 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7377 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7379 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7380 dd->comm->cd[d].np_dlb > 0)
7382 if (np > dd->comm->cd[d].np_dlb)
7387 /* If a current local cell size is smaller than the requested
7388 * cut-off, we could still fix it, but this gets very complicated.
7389 * Without fixing here, we might actually need more checks.
7391 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7398 if (dd->comm->eDLB != edlbNO)
7400 /* If DLB is not active yet, we don't need to check the grid jumps.
7401 * Actually we shouldn't, because then the grid jump data is not set.
7403 if (dd->comm->bDynLoadBal &&
7404 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7409 gmx_sumi(1, &LocallyLimited, cr);
7411 if (LocallyLimited > 0)
7420 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7423 gmx_bool bCutoffAllowed;
7425 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7429 cr->dd->comm->cutoff = cutoff_req;
7432 return bCutoffAllowed;
7435 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7437 gmx_domdec_comm_t *comm;
7439 comm = cr->dd->comm;
7441 /* Turn on the DLB limiting (might have been on already) */
7442 comm->bPMELoadBalDLBLimits = TRUE;
7444 /* Change the cut-off limit */
7445 comm->PMELoadBal_max_cutoff = comm->cutoff;
7448 static void merge_cg_buffers(int ncell,
7449 gmx_domdec_comm_dim_t *cd, int pulse,
7451 int *index_gl, int *recv_i,
7452 rvec *cg_cm, rvec *recv_vr,
7454 cginfo_mb_t *cginfo_mb, int *cginfo)
7456 gmx_domdec_ind_t *ind, *ind_p;
7457 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7458 int shift, shift_at;
7460 ind = &cd->ind[pulse];
7462 /* First correct the already stored data */
7463 shift = ind->nrecv[ncell];
7464 for (cell = ncell-1; cell >= 0; cell--)
7466 shift -= ind->nrecv[cell];
7469 /* Move the cg's present from previous grid pulses */
7470 cg0 = ncg_cell[ncell+cell];
7471 cg1 = ncg_cell[ncell+cell+1];
7472 cgindex[cg1+shift] = cgindex[cg1];
7473 for (cg = cg1-1; cg >= cg0; cg--)
7475 index_gl[cg+shift] = index_gl[cg];
7476 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7477 cgindex[cg+shift] = cgindex[cg];
7478 cginfo[cg+shift] = cginfo[cg];
7480 /* Correct the already stored send indices for the shift */
7481 for (p = 1; p <= pulse; p++)
7483 ind_p = &cd->ind[p];
7485 for (c = 0; c < cell; c++)
7487 cg0 += ind_p->nsend[c];
7489 cg1 = cg0 + ind_p->nsend[cell];
7490 for (cg = cg0; cg < cg1; cg++)
7492 ind_p->index[cg] += shift;
7498 /* Merge in the communicated buffers */
7502 for (cell = 0; cell < ncell; cell++)
7504 cg1 = ncg_cell[ncell+cell+1] + shift;
7507 /* Correct the old cg indices */
7508 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7510 cgindex[cg+1] += shift_at;
7513 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7515 /* Copy this charge group from the buffer */
7516 index_gl[cg1] = recv_i[cg0];
7517 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7518 /* Add it to the cgindex */
7519 cg_gl = index_gl[cg1];
7520 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7521 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7522 cgindex[cg1+1] = cgindex[cg1] + nat;
7527 shift += ind->nrecv[cell];
7528 ncg_cell[ncell+cell+1] = cg1;
7532 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7533 int nzone, int cg0, const int *cgindex)
7537 /* Store the atom block boundaries for easy copying of communication buffers
7540 for (zone = 0; zone < nzone; zone++)
7542 for (p = 0; p < cd->np; p++)
7544 cd->ind[p].cell2at0[zone] = cgindex[cg];
7545 cg += cd->ind[p].nrecv[zone];
7546 cd->ind[p].cell2at1[zone] = cgindex[cg];
7551 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7557 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7559 if (!bLocalCG[link->a[i]])
7568 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7570 real c[DIM][4]; /* the corners for the non-bonded communication */
7571 real cr0; /* corner for rounding */
7572 real cr1[4]; /* corners for rounding */
7573 real bc[DIM]; /* corners for bounded communication */
7574 real bcr1; /* corner for rounding for bonded communication */
7577 /* Determine the corners of the domain(s) we are communicating with */
7579 set_dd_corners(const gmx_domdec_t *dd,
7580 int dim0, int dim1, int dim2,
7584 const gmx_domdec_comm_t *comm;
7585 const gmx_domdec_zones_t *zones;
7590 zones = &comm->zones;
7592 /* Keep the compiler happy */
7596 /* The first dimension is equal for all cells */
7597 c->c[0][0] = comm->cell_x0[dim0];
7600 c->bc[0] = c->c[0][0];
7605 /* This cell row is only seen from the first row */
7606 c->c[1][0] = comm->cell_x0[dim1];
7607 /* All rows can see this row */
7608 c->c[1][1] = comm->cell_x0[dim1];
7611 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7614 /* For the multi-body distance we need the maximum */
7615 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7618 /* Set the upper-right corner for rounding */
7619 c->cr0 = comm->cell_x1[dim0];
7624 for (j = 0; j < 4; j++)
7626 c->c[2][j] = comm->cell_x0[dim2];
7630 /* Use the maximum of the i-cells that see a j-cell */
7631 for (i = 0; i < zones->nizone; i++)
7633 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7639 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7645 /* For the multi-body distance we need the maximum */
7646 c->bc[2] = comm->cell_x0[dim2];
7647 for (i = 0; i < 2; i++)
7649 for (j = 0; j < 2; j++)
7651 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7657 /* Set the upper-right corner for rounding */
7658 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7659 * Only cell (0,0,0) can see cell 7 (1,1,1)
7661 c->cr1[0] = comm->cell_x1[dim1];
7662 c->cr1[3] = comm->cell_x1[dim1];
7665 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7668 /* For the multi-body distance we need the maximum */
7669 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7676 /* Determine which cg's we need to send in this pulse from this zone */
7678 get_zone_pulse_cgs(gmx_domdec_t *dd,
7679 int zonei, int zone,
7681 const int *index_gl,
7683 int dim, int dim_ind,
7684 int dim0, int dim1, int dim2,
7685 real r_comm2, real r_bcomm2,
7689 real skew_fac2_d, real skew_fac_01,
7690 rvec *v_d, rvec *v_0, rvec *v_1,
7691 const dd_corners_t *c,
7693 gmx_bool bDistBonded,
7699 gmx_domdec_ind_t *ind,
7700 int **ibuf, int *ibuf_nalloc,
7706 gmx_domdec_comm_t *comm;
7708 gmx_bool bDistMB_pulse;
7710 real r2, rb2, r, tric_sh;
7713 int nsend_z, nsend, nat;
7717 bScrew = (dd->bScrewPBC && dim == XX);
7719 bDistMB_pulse = (bDistMB && bDistBonded);
7725 for (cg = cg0; cg < cg1; cg++)
7729 if (tric_dist[dim_ind] == 0)
7731 /* Rectangular direction, easy */
7732 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7739 r = cg_cm[cg][dim] - c->bc[dim_ind];
7745 /* Rounding gives at most a 16% reduction
7746 * in communicated atoms
7748 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7750 r = cg_cm[cg][dim0] - c->cr0;
7751 /* This is the first dimension, so always r >= 0 */
7758 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7760 r = cg_cm[cg][dim1] - c->cr1[zone];
7767 r = cg_cm[cg][dim1] - c->bcr1;
7777 /* Triclinic direction, more complicated */
7780 /* Rounding, conservative as the skew_fac multiplication
7781 * will slightly underestimate the distance.
7783 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7785 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7786 for (i = dim0+1; i < DIM; i++)
7788 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7790 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7793 rb[dim0] = rn[dim0];
7796 /* Take care that the cell planes along dim0 might not
7797 * be orthogonal to those along dim1 and dim2.
7799 for (i = 1; i <= dim_ind; i++)
7802 if (normal[dim0][dimd] > 0)
7804 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7807 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7812 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7814 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7816 for (i = dim1+1; i < DIM; i++)
7818 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7820 rn[dim1] += tric_sh;
7823 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7824 /* Take care of coupling of the distances
7825 * to the planes along dim0 and dim1 through dim2.
7827 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7828 /* Take care that the cell planes along dim1
7829 * might not be orthogonal to that along dim2.
7831 if (normal[dim1][dim2] > 0)
7833 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7839 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7842 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7843 /* Take care of coupling of the distances
7844 * to the planes along dim0 and dim1 through dim2.
7846 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7847 /* Take care that the cell planes along dim1
7848 * might not be orthogonal to that along dim2.
7850 if (normal[dim1][dim2] > 0)
7852 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7857 /* The distance along the communication direction */
7858 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7860 for (i = dim+1; i < DIM; i++)
7862 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7867 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7868 /* Take care of coupling of the distances
7869 * to the planes along dim0 and dim1 through dim2.
7871 if (dim_ind == 1 && zonei == 1)
7873 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7879 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7882 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7883 /* Take care of coupling of the distances
7884 * to the planes along dim0 and dim1 through dim2.
7886 if (dim_ind == 1 && zonei == 1)
7888 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7896 ((bDistMB && rb2 < r_bcomm2) ||
7897 (bDist2B && r2 < r_bcomm2)) &&
7899 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7900 missing_link(comm->cglink, index_gl[cg],
7903 /* Make an index to the local charge groups */
7904 if (nsend+1 > ind->nalloc)
7906 ind->nalloc = over_alloc_large(nsend+1);
7907 srenew(ind->index, ind->nalloc);
7909 if (nsend+1 > *ibuf_nalloc)
7911 *ibuf_nalloc = over_alloc_large(nsend+1);
7912 srenew(*ibuf, *ibuf_nalloc);
7914 ind->index[nsend] = cg;
7915 (*ibuf)[nsend] = index_gl[cg];
7917 vec_rvec_check_alloc(vbuf, nsend+1);
7919 if (dd->ci[dim] == 0)
7921 /* Correct cg_cm for pbc */
7922 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7925 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7926 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7931 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7934 nat += cgindex[cg+1] - cgindex[cg];
7940 *nsend_z_ptr = nsend_z;
7943 static void setup_dd_communication(gmx_domdec_t *dd,
7944 matrix box, gmx_ddbox_t *ddbox,
7945 t_forcerec *fr, t_state *state, rvec **f)
7947 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7948 int nzone, nzone_send, zone, zonei, cg0, cg1;
7949 int c, i, j, cg, cg_gl, nrcg;
7950 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7951 gmx_domdec_comm_t *comm;
7952 gmx_domdec_zones_t *zones;
7953 gmx_domdec_comm_dim_t *cd;
7954 gmx_domdec_ind_t *ind;
7955 cginfo_mb_t *cginfo_mb;
7956 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7957 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7958 dd_corners_t corners;
7960 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7961 real skew_fac2_d, skew_fac_01;
7968 fprintf(debug, "Setting up DD communication\n");
7973 switch (fr->cutoff_scheme)
7982 gmx_incons("unimplemented");
7986 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7988 dim = dd->dim[dim_ind];
7990 /* Check if we need to use triclinic distances */
7991 tric_dist[dim_ind] = 0;
7992 for (i = 0; i <= dim_ind; i++)
7994 if (ddbox->tric_dir[dd->dim[i]])
7996 tric_dist[dim_ind] = 1;
8001 bBondComm = comm->bBondComm;
8003 /* Do we need to determine extra distances for multi-body bondeds? */
8004 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8006 /* Do we need to determine extra distances for only two-body bondeds? */
8007 bDist2B = (bBondComm && !bDistMB);
8009 r_comm2 = sqr(comm->cutoff);
8010 r_bcomm2 = sqr(comm->cutoff_mbody);
8014 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8017 zones = &comm->zones;
8020 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8021 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8023 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8025 /* Triclinic stuff */
8026 normal = ddbox->normal;
8030 v_0 = ddbox->v[dim0];
8031 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8033 /* Determine the coupling coefficient for the distances
8034 * to the cell planes along dim0 and dim1 through dim2.
8035 * This is required for correct rounding.
8038 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8041 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8047 v_1 = ddbox->v[dim1];
8050 zone_cg_range = zones->cg_range;
8051 index_gl = dd->index_gl;
8052 cgindex = dd->cgindex;
8053 cginfo_mb = fr->cginfo_mb;
8055 zone_cg_range[0] = 0;
8056 zone_cg_range[1] = dd->ncg_home;
8057 comm->zone_ncg1[0] = dd->ncg_home;
8058 pos_cg = dd->ncg_home;
8060 nat_tot = dd->nat_home;
8062 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8064 dim = dd->dim[dim_ind];
8065 cd = &comm->cd[dim_ind];
8067 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8069 /* No pbc in this dimension, the first node should not comm. */
8077 v_d = ddbox->v[dim];
8078 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8080 cd->bInPlace = TRUE;
8081 for (p = 0; p < cd->np; p++)
8083 /* Only atoms communicated in the first pulse are used
8084 * for multi-body bonded interactions or for bBondComm.
8086 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8091 for (zone = 0; zone < nzone_send; zone++)
8093 if (tric_dist[dim_ind] && dim_ind > 0)
8095 /* Determine slightly more optimized skew_fac's
8097 * This reduces the number of communicated atoms
8098 * by about 10% for 3D DD of rhombic dodecahedra.
8100 for (dimd = 0; dimd < dim; dimd++)
8102 sf2_round[dimd] = 1;
8103 if (ddbox->tric_dir[dimd])
8105 for (i = dd->dim[dimd]+1; i < DIM; i++)
8107 /* If we are shifted in dimension i
8108 * and the cell plane is tilted forward
8109 * in dimension i, skip this coupling.
8111 if (!(zones->shift[nzone+zone][i] &&
8112 ddbox->v[dimd][i][dimd] >= 0))
8115 sqr(ddbox->v[dimd][i][dimd]);
8118 sf2_round[dimd] = 1/sf2_round[dimd];
8123 zonei = zone_perm[dim_ind][zone];
8126 /* Here we permutate the zones to obtain a convenient order
8127 * for neighbor searching
8129 cg0 = zone_cg_range[zonei];
8130 cg1 = zone_cg_range[zonei+1];
8134 /* Look only at the cg's received in the previous grid pulse
8136 cg1 = zone_cg_range[nzone+zone+1];
8137 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8140 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8141 for (th = 0; th < comm->nth; th++)
8143 gmx_domdec_ind_t *ind_p;
8144 int **ibuf_p, *ibuf_nalloc_p;
8146 int *nsend_p, *nat_p;
8152 /* Thread 0 writes in the comm buffers */
8154 ibuf_p = &comm->buf_int;
8155 ibuf_nalloc_p = &comm->nalloc_int;
8156 vbuf_p = &comm->vbuf;
8159 nsend_zone_p = &ind->nsend[zone];
8163 /* Other threads write into temp buffers */
8164 ind_p = &comm->dth[th].ind;
8165 ibuf_p = &comm->dth[th].ibuf;
8166 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8167 vbuf_p = &comm->dth[th].vbuf;
8168 nsend_p = &comm->dth[th].nsend;
8169 nat_p = &comm->dth[th].nat;
8170 nsend_zone_p = &comm->dth[th].nsend_zone;
8172 comm->dth[th].nsend = 0;
8173 comm->dth[th].nat = 0;
8174 comm->dth[th].nsend_zone = 0;
8184 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8185 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8188 /* Get the cg's for this pulse in this zone */
8189 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8191 dim, dim_ind, dim0, dim1, dim2,
8194 normal, skew_fac2_d, skew_fac_01,
8195 v_d, v_0, v_1, &corners, sf2_round,
8196 bDistBonded, bBondComm,
8200 ibuf_p, ibuf_nalloc_p,
8206 /* Append data of threads>=1 to the communication buffers */
8207 for (th = 1; th < comm->nth; th++)
8209 dd_comm_setup_work_t *dth;
8212 dth = &comm->dth[th];
8214 ns1 = nsend + dth->nsend_zone;
8215 if (ns1 > ind->nalloc)
8217 ind->nalloc = over_alloc_dd(ns1);
8218 srenew(ind->index, ind->nalloc);
8220 if (ns1 > comm->nalloc_int)
8222 comm->nalloc_int = over_alloc_dd(ns1);
8223 srenew(comm->buf_int, comm->nalloc_int);
8225 if (ns1 > comm->vbuf.nalloc)
8227 comm->vbuf.nalloc = over_alloc_dd(ns1);
8228 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8231 for (i = 0; i < dth->nsend_zone; i++)
8233 ind->index[nsend] = dth->ind.index[i];
8234 comm->buf_int[nsend] = dth->ibuf[i];
8235 copy_rvec(dth->vbuf.v[i],
8236 comm->vbuf.v[nsend]);
8240 ind->nsend[zone] += dth->nsend_zone;
8243 /* Clear the counts in case we do not have pbc */
8244 for (zone = nzone_send; zone < nzone; zone++)
8246 ind->nsend[zone] = 0;
8248 ind->nsend[nzone] = nsend;
8249 ind->nsend[nzone+1] = nat;
8250 /* Communicate the number of cg's and atoms to receive */
8251 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8252 ind->nsend, nzone+2,
8253 ind->nrecv, nzone+2);
8255 /* The rvec buffer is also required for atom buffers of size nsend
8256 * in dd_move_x and dd_move_f.
8258 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8262 /* We can receive in place if only the last zone is not empty */
8263 for (zone = 0; zone < nzone-1; zone++)
8265 if (ind->nrecv[zone] > 0)
8267 cd->bInPlace = FALSE;
8272 /* The int buffer is only required here for the cg indices */
8273 if (ind->nrecv[nzone] > comm->nalloc_int2)
8275 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8276 srenew(comm->buf_int2, comm->nalloc_int2);
8278 /* The rvec buffer is also required for atom buffers
8279 * of size nrecv in dd_move_x and dd_move_f.
8281 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8282 vec_rvec_check_alloc(&comm->vbuf2, i);
8286 /* Make space for the global cg indices */
8287 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8288 || dd->cg_nalloc == 0)
8290 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8291 srenew(index_gl, dd->cg_nalloc);
8292 srenew(cgindex, dd->cg_nalloc+1);
8294 /* Communicate the global cg indices */
8297 recv_i = index_gl + pos_cg;
8301 recv_i = comm->buf_int2;
8303 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8304 comm->buf_int, nsend,
8305 recv_i, ind->nrecv[nzone]);
8307 /* Make space for cg_cm */
8308 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8309 if (fr->cutoff_scheme == ecutsGROUP)
8317 /* Communicate cg_cm */
8320 recv_vr = cg_cm + pos_cg;
8324 recv_vr = comm->vbuf2.v;
8326 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8327 comm->vbuf.v, nsend,
8328 recv_vr, ind->nrecv[nzone]);
8330 /* Make the charge group index */
8333 zone = (p == 0 ? 0 : nzone - 1);
8334 while (zone < nzone)
8336 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8338 cg_gl = index_gl[pos_cg];
8339 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8340 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8341 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8344 /* Update the charge group presence,
8345 * so we can use it in the next pass of the loop.
8347 comm->bLocalCG[cg_gl] = TRUE;
8353 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8356 zone_cg_range[nzone+zone] = pos_cg;
8361 /* This part of the code is never executed with bBondComm. */
8362 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8363 index_gl, recv_i, cg_cm, recv_vr,
8364 cgindex, fr->cginfo_mb, fr->cginfo);
8365 pos_cg += ind->nrecv[nzone];
8367 nat_tot += ind->nrecv[nzone+1];
8371 /* Store the atom block for easy copying of communication buffers */
8372 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8376 dd->index_gl = index_gl;
8377 dd->cgindex = cgindex;
8379 dd->ncg_tot = zone_cg_range[zones->n];
8380 dd->nat_tot = nat_tot;
8381 comm->nat[ddnatHOME] = dd->nat_home;
8382 for (i = ddnatZONE; i < ddnatNR; i++)
8384 comm->nat[i] = dd->nat_tot;
8389 /* We don't need to update cginfo, since that was alrady done above.
8390 * So we pass NULL for the forcerec.
8392 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8393 NULL, comm->bLocalCG);
8398 fprintf(debug, "Finished setting up DD communication, zones:");
8399 for (c = 0; c < zones->n; c++)
8401 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8403 fprintf(debug, "\n");
8407 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8411 for (c = 0; c < zones->nizone; c++)
8413 zones->izone[c].cg1 = zones->cg_range[c+1];
8414 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8415 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8419 static void set_zones_size(gmx_domdec_t *dd,
8420 matrix box, const gmx_ddbox_t *ddbox,
8421 int zone_start, int zone_end)
8423 gmx_domdec_comm_t *comm;
8424 gmx_domdec_zones_t *zones;
8426 int z, zi, zj0, zj1, d, dim;
8429 real size_j, add_tric;
8434 zones = &comm->zones;
8436 /* Do we need to determine extra distances for multi-body bondeds? */
8437 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8439 for (z = zone_start; z < zone_end; z++)
8441 /* Copy cell limits to zone limits.
8442 * Valid for non-DD dims and non-shifted dims.
8444 copy_rvec(comm->cell_x0, zones->size[z].x0);
8445 copy_rvec(comm->cell_x1, zones->size[z].x1);
8448 for (d = 0; d < dd->ndim; d++)
8452 for (z = 0; z < zones->n; z++)
8454 /* With a staggered grid we have different sizes
8455 * for non-shifted dimensions.
8457 if (dd->bGridJump && zones->shift[z][dim] == 0)
8461 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8462 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8466 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8467 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8473 rcmbs = comm->cutoff_mbody;
8474 if (ddbox->tric_dir[dim])
8476 rcs /= ddbox->skew_fac[dim];
8477 rcmbs /= ddbox->skew_fac[dim];
8480 /* Set the lower limit for the shifted zone dimensions */
8481 for (z = zone_start; z < zone_end; z++)
8483 if (zones->shift[z][dim] > 0)
8486 if (!dd->bGridJump || d == 0)
8488 zones->size[z].x0[dim] = comm->cell_x1[dim];
8489 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8493 /* Here we take the lower limit of the zone from
8494 * the lowest domain of the zone below.
8498 zones->size[z].x0[dim] =
8499 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8505 zones->size[z].x0[dim] =
8506 zones->size[zone_perm[2][z-4]].x0[dim];
8510 zones->size[z].x0[dim] =
8511 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8514 /* A temporary limit, is updated below */
8515 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8519 for (zi = 0; zi < zones->nizone; zi++)
8521 if (zones->shift[zi][dim] == 0)
8523 /* This takes the whole zone into account.
8524 * With multiple pulses this will lead
8525 * to a larger zone then strictly necessary.
8527 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8528 zones->size[zi].x1[dim]+rcmbs);
8536 /* Loop over the i-zones to set the upper limit of each
8539 for (zi = 0; zi < zones->nizone; zi++)
8541 if (zones->shift[zi][dim] == 0)
8543 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8545 if (zones->shift[z][dim] > 0)
8547 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8548 zones->size[zi].x1[dim]+rcs);
8555 for (z = zone_start; z < zone_end; z++)
8557 /* Initialization only required to keep the compiler happy */
8558 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8561 /* To determine the bounding box for a zone we need to find
8562 * the extreme corners of 4, 2 or 1 corners.
8564 nc = 1 << (ddbox->npbcdim - 1);
8566 for (c = 0; c < nc; c++)
8568 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8572 corner[YY] = zones->size[z].x0[YY];
8576 corner[YY] = zones->size[z].x1[YY];
8580 corner[ZZ] = zones->size[z].x0[ZZ];
8584 corner[ZZ] = zones->size[z].x1[ZZ];
8586 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8588 /* With 1D domain decomposition the cg's are not in
8589 * the triclinic box, but triclinic x-y and rectangular y-z.
8590 * Shift y back, so it will later end up at 0.
8592 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8594 /* Apply the triclinic couplings */
8595 for (i = YY; i < ddbox->npbcdim; i++)
8597 for (j = XX; j < i; j++)
8599 corner[j] += corner[i]*box[i][j]/box[i][i];
8604 copy_rvec(corner, corner_min);
8605 copy_rvec(corner, corner_max);
8609 for (i = 0; i < DIM; i++)
8611 corner_min[i] = min(corner_min[i], corner[i]);
8612 corner_max[i] = max(corner_max[i], corner[i]);
8616 /* Copy the extreme cornes without offset along x */
8617 for (i = 0; i < DIM; i++)
8619 zones->size[z].bb_x0[i] = corner_min[i];
8620 zones->size[z].bb_x1[i] = corner_max[i];
8622 /* Add the offset along x */
8623 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8624 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8627 if (zone_start == 0)
8630 for (dim = 0; dim < DIM; dim++)
8632 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8634 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8639 for (z = zone_start; z < zone_end; z++)
8641 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8643 zones->size[z].x0[XX], zones->size[z].x1[XX],
8644 zones->size[z].x0[YY], zones->size[z].x1[YY],
8645 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8646 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8648 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8649 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8650 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8655 static int comp_cgsort(const void *a, const void *b)
8659 gmx_cgsort_t *cga, *cgb;
8660 cga = (gmx_cgsort_t *)a;
8661 cgb = (gmx_cgsort_t *)b;
8663 comp = cga->nsc - cgb->nsc;
8666 comp = cga->ind_gl - cgb->ind_gl;
8672 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8677 /* Order the data */
8678 for (i = 0; i < n; i++)
8680 buf[i] = a[sort[i].ind];
8683 /* Copy back to the original array */
8684 for (i = 0; i < n; i++)
8690 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8695 /* Order the data */
8696 for (i = 0; i < n; i++)
8698 copy_rvec(v[sort[i].ind], buf[i]);
8701 /* Copy back to the original array */
8702 for (i = 0; i < n; i++)
8704 copy_rvec(buf[i], v[i]);
8708 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8711 int a, atot, cg, cg0, cg1, i;
8713 if (cgindex == NULL)
8715 /* Avoid the useless loop of the atoms within a cg */
8716 order_vec_cg(ncg, sort, v, buf);
8721 /* Order the data */
8723 for (cg = 0; cg < ncg; cg++)
8725 cg0 = cgindex[sort[cg].ind];
8726 cg1 = cgindex[sort[cg].ind+1];
8727 for (i = cg0; i < cg1; i++)
8729 copy_rvec(v[i], buf[a]);
8735 /* Copy back to the original array */
8736 for (a = 0; a < atot; a++)
8738 copy_rvec(buf[a], v[a]);
8742 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8743 int nsort_new, gmx_cgsort_t *sort_new,
8744 gmx_cgsort_t *sort1)
8748 /* The new indices are not very ordered, so we qsort them */
8749 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8751 /* sort2 is already ordered, so now we can merge the two arrays */
8755 while (i2 < nsort2 || i_new < nsort_new)
8759 sort1[i1++] = sort_new[i_new++];
8761 else if (i_new == nsort_new)
8763 sort1[i1++] = sort2[i2++];
8765 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8766 (sort2[i2].nsc == sort_new[i_new].nsc &&
8767 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8769 sort1[i1++] = sort2[i2++];
8773 sort1[i1++] = sort_new[i_new++];
8778 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8780 gmx_domdec_sort_t *sort;
8781 gmx_cgsort_t *cgsort, *sort_i;
8782 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8783 int sort_last, sort_skip;
8785 sort = dd->comm->sort;
8787 a = fr->ns.grid->cell_index;
8789 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8791 if (ncg_home_old >= 0)
8793 /* The charge groups that remained in the same ns grid cell
8794 * are completely ordered. So we can sort efficiently by sorting
8795 * the charge groups that did move into the stationary list.
8800 for (i = 0; i < dd->ncg_home; i++)
8802 /* Check if this cg did not move to another node */
8805 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8807 /* This cg is new on this node or moved ns grid cell */
8808 if (nsort_new >= sort->sort_new_nalloc)
8810 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8811 srenew(sort->sort_new, sort->sort_new_nalloc);
8813 sort_i = &(sort->sort_new[nsort_new++]);
8817 /* This cg did not move */
8818 sort_i = &(sort->sort2[nsort2++]);
8820 /* Sort on the ns grid cell indices
8821 * and the global topology index.
8822 * index_gl is irrelevant with cell ns,
8823 * but we set it here anyhow to avoid a conditional.
8826 sort_i->ind_gl = dd->index_gl[i];
8833 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8836 /* Sort efficiently */
8837 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8842 cgsort = sort->sort;
8844 for (i = 0; i < dd->ncg_home; i++)
8846 /* Sort on the ns grid cell indices
8847 * and the global topology index
8849 cgsort[i].nsc = a[i];
8850 cgsort[i].ind_gl = dd->index_gl[i];
8852 if (cgsort[i].nsc < moved)
8859 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8861 /* Determine the order of the charge groups using qsort */
8862 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8868 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8871 int ncg_new, i, *a, na;
8873 sort = dd->comm->sort->sort;
8875 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8878 for (i = 0; i < na; i++)
8882 sort[ncg_new].ind = a[i];
8890 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8891 rvec *cgcm, t_forcerec *fr, t_state *state,
8894 gmx_domdec_sort_t *sort;
8895 gmx_cgsort_t *cgsort, *sort_i;
8897 int ncg_new, i, *ibuf, cgsize;
8900 sort = dd->comm->sort;
8902 if (dd->ncg_home > sort->sort_nalloc)
8904 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8905 srenew(sort->sort, sort->sort_nalloc);
8906 srenew(sort->sort2, sort->sort_nalloc);
8908 cgsort = sort->sort;
8910 switch (fr->cutoff_scheme)
8913 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8916 ncg_new = dd_sort_order_nbnxn(dd, fr);
8919 gmx_incons("unimplemented");
8923 /* We alloc with the old size, since cgindex is still old */
8924 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8925 vbuf = dd->comm->vbuf.v;
8929 cgindex = dd->cgindex;
8936 /* Remove the charge groups which are no longer at home here */
8937 dd->ncg_home = ncg_new;
8940 fprintf(debug, "Set the new home charge group count to %d\n",
8944 /* Reorder the state */
8945 for (i = 0; i < estNR; i++)
8947 if (EST_DISTR(i) && (state->flags & (1<<i)))
8952 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8955 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8958 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8961 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8965 case estDISRE_INITF:
8966 case estDISRE_RM3TAV:
8967 case estORIRE_INITF:
8969 /* No ordering required */
8972 gmx_incons("Unknown state entry encountered in dd_sort_state");
8977 if (fr->cutoff_scheme == ecutsGROUP)
8980 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8983 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8985 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8986 srenew(sort->ibuf, sort->ibuf_nalloc);
8989 /* Reorder the global cg index */
8990 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
8991 /* Reorder the cginfo */
8992 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
8993 /* Rebuild the local cg index */
8997 for (i = 0; i < dd->ncg_home; i++)
8999 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9000 ibuf[i+1] = ibuf[i] + cgsize;
9002 for (i = 0; i < dd->ncg_home+1; i++)
9004 dd->cgindex[i] = ibuf[i];
9009 for (i = 0; i < dd->ncg_home+1; i++)
9014 /* Set the home atom number */
9015 dd->nat_home = dd->cgindex[dd->ncg_home];
9017 if (fr->cutoff_scheme == ecutsVERLET)
9019 /* The atoms are now exactly in grid order, update the grid order */
9020 nbnxn_set_atomorder(fr->nbv->nbs);
9024 /* Copy the sorted ns cell indices back to the ns grid struct */
9025 for (i = 0; i < dd->ncg_home; i++)
9027 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9029 fr->ns.grid->nr = dd->ncg_home;
9033 static void add_dd_statistics(gmx_domdec_t *dd)
9035 gmx_domdec_comm_t *comm;
9040 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9042 comm->sum_nat[ddnat-ddnatZONE] +=
9043 comm->nat[ddnat] - comm->nat[ddnat-1];
9048 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9050 gmx_domdec_comm_t *comm;
9055 /* Reset all the statistics and counters for total run counting */
9056 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9058 comm->sum_nat[ddnat-ddnatZONE] = 0;
9062 comm->load_step = 0;
9065 clear_ivec(comm->load_lim);
9070 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9072 gmx_domdec_comm_t *comm;
9076 comm = cr->dd->comm;
9078 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9085 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");
9087 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9089 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9094 " av. #atoms communicated per step for force: %d x %.1f\n",
9098 if (cr->dd->vsite_comm)
9101 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9102 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9107 if (cr->dd->constraint_comm)
9110 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9111 1 + ir->nLincsIter, av);
9115 gmx_incons(" Unknown type for DD statistics");
9118 fprintf(fplog, "\n");
9120 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9122 print_dd_load_av(fplog, cr->dd);
9126 void dd_partition_system(FILE *fplog,
9127 gmx_large_int_t step,
9129 gmx_bool bMasterState,
9131 t_state *state_global,
9132 gmx_mtop_t *top_global,
9134 t_state *state_local,
9137 gmx_localtop_t *top_local,
9140 gmx_shellfc_t shellfc,
9141 gmx_constr_t constr,
9143 gmx_wallcycle_t wcycle,
9147 gmx_domdec_comm_t *comm;
9148 gmx_ddbox_t ddbox = {0};
9150 gmx_large_int_t step_pcoupl;
9151 rvec cell_ns_x0, cell_ns_x1;
9152 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9153 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9154 gmx_bool bRedist, bSortCG, bResortAll;
9155 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9162 bBoxChanged = (bMasterState || DEFORM(*ir));
9163 if (ir->epc != epcNO)
9165 /* With nstpcouple > 1 pressure coupling happens.
9166 * one step after calculating the pressure.
9167 * Box scaling happens at the end of the MD step,
9168 * after the DD partitioning.
9169 * We therefore have to do DLB in the first partitioning
9170 * after an MD step where P-coupling occured.
9171 * We need to determine the last step in which p-coupling occurred.
9172 * MRS -- need to validate this for vv?
9177 step_pcoupl = step - 1;
9181 step_pcoupl = ((step - 1)/n)*n + 1;
9183 if (step_pcoupl >= comm->partition_step)
9189 bNStGlobalComm = (step % nstglobalcomm == 0);
9191 if (!comm->bDynLoadBal)
9197 /* Should we do dynamic load balacing this step?
9198 * Since it requires (possibly expensive) global communication,
9199 * we might want to do DLB less frequently.
9201 if (bBoxChanged || ir->epc != epcNO)
9203 bDoDLB = bBoxChanged;
9207 bDoDLB = bNStGlobalComm;
9211 /* Check if we have recorded loads on the nodes */
9212 if (comm->bRecordLoad && dd_load_count(comm))
9214 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9216 /* Check if we should use DLB at the second partitioning
9217 * and every 100 partitionings,
9218 * so the extra communication cost is negligible.
9220 n = max(100, nstglobalcomm);
9221 bCheckDLB = (comm->n_load_collect == 0 ||
9222 comm->n_load_have % n == n-1);
9229 /* Print load every nstlog, first and last step to the log file */
9230 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9231 comm->n_load_collect == 0 ||
9233 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9235 /* Avoid extra communication due to verbose screen output
9236 * when nstglobalcomm is set.
9238 if (bDoDLB || bLogLoad || bCheckDLB ||
9239 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9241 get_load_distribution(dd, wcycle);
9246 dd_print_load(fplog, dd, step-1);
9250 dd_print_load_verbose(dd);
9253 comm->n_load_collect++;
9257 /* Since the timings are node dependent, the master decides */
9261 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9264 fprintf(debug, "step %s, imb loss %f\n",
9265 gmx_step_str(step, sbuf),
9266 dd_force_imb_perf_loss(dd));
9269 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9272 turn_on_dlb(fplog, cr, step);
9277 comm->n_load_have++;
9280 cgs_gl = &comm->cgs_gl;
9285 /* Clear the old state */
9286 clear_dd_indices(dd, 0, 0);
9288 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9289 TRUE, cgs_gl, state_global->x, &ddbox);
9291 get_cg_distribution(fplog, step, dd, cgs_gl,
9292 state_global->box, &ddbox, state_global->x);
9294 dd_distribute_state(dd, cgs_gl,
9295 state_global, state_local, f);
9297 dd_make_local_cgs(dd, &top_local->cgs);
9299 /* Ensure that we have space for the new distribution */
9300 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9302 if (fr->cutoff_scheme == ecutsGROUP)
9304 calc_cgcm(fplog, 0, dd->ncg_home,
9305 &top_local->cgs, state_local->x, fr->cg_cm);
9308 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9310 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9314 else if (state_local->ddp_count != dd->ddp_count)
9316 if (state_local->ddp_count > dd->ddp_count)
9318 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9321 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9323 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);
9326 /* Clear the old state */
9327 clear_dd_indices(dd, 0, 0);
9329 /* Build the new indices */
9330 rebuild_cgindex(dd, cgs_gl->index, state_local);
9331 make_dd_indices(dd, cgs_gl->index, 0);
9333 if (fr->cutoff_scheme == ecutsGROUP)
9335 /* Redetermine the cg COMs */
9336 calc_cgcm(fplog, 0, dd->ncg_home,
9337 &top_local->cgs, state_local->x, fr->cg_cm);
9340 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9342 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9344 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9345 TRUE, &top_local->cgs, state_local->x, &ddbox);
9347 bRedist = comm->bDynLoadBal;
9351 /* We have the full state, only redistribute the cgs */
9353 /* Clear the non-home indices */
9354 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9356 /* Avoid global communication for dim's without pbc and -gcom */
9357 if (!bNStGlobalComm)
9359 copy_rvec(comm->box0, ddbox.box0 );
9360 copy_rvec(comm->box_size, ddbox.box_size);
9362 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9363 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9368 /* For dim's without pbc and -gcom */
9369 copy_rvec(ddbox.box0, comm->box0 );
9370 copy_rvec(ddbox.box_size, comm->box_size);
9372 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9375 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9377 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9380 /* Check if we should sort the charge groups */
9381 if (comm->nstSortCG > 0)
9383 bSortCG = (bMasterState ||
9384 (bRedist && (step % comm->nstSortCG == 0)));
9391 ncg_home_old = dd->ncg_home;
9396 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9398 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9399 state_local, f, fr, mdatoms,
9400 !bSortCG, nrnb, &cg0, &ncg_moved);
9402 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9405 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9407 &comm->cell_x0, &comm->cell_x1,
9408 dd->ncg_home, fr->cg_cm,
9409 cell_ns_x0, cell_ns_x1, &grid_density);
9413 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9416 switch (fr->cutoff_scheme)
9419 copy_ivec(fr->ns.grid->n, ncells_old);
9420 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9421 state_local->box, cell_ns_x0, cell_ns_x1,
9422 fr->rlistlong, grid_density);
9425 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9428 gmx_incons("unimplemented");
9430 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9431 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9435 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9437 /* Sort the state on charge group position.
9438 * This enables exact restarts from this step.
9439 * It also improves performance by about 15% with larger numbers
9440 * of atoms per node.
9443 /* Fill the ns grid with the home cell,
9444 * so we can sort with the indices.
9446 set_zones_ncg_home(dd);
9448 switch (fr->cutoff_scheme)
9451 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9453 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9455 comm->zones.size[0].bb_x0,
9456 comm->zones.size[0].bb_x1,
9458 comm->zones.dens_zone0,
9461 ncg_moved, bRedist ? comm->moved : NULL,
9462 fr->nbv->grp[eintLocal].kernel_type,
9463 fr->nbv->grp[eintLocal].nbat);
9465 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9468 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9469 0, dd->ncg_home, fr->cg_cm);
9471 copy_ivec(fr->ns.grid->n, ncells_new);
9474 gmx_incons("unimplemented");
9477 bResortAll = bMasterState;
9479 /* Check if we can user the old order and ns grid cell indices
9480 * of the charge groups to sort the charge groups efficiently.
9482 if (ncells_new[XX] != ncells_old[XX] ||
9483 ncells_new[YY] != ncells_old[YY] ||
9484 ncells_new[ZZ] != ncells_old[ZZ])
9491 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9492 gmx_step_str(step, sbuf), dd->ncg_home);
9494 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9495 bResortAll ? -1 : ncg_home_old);
9496 /* Rebuild all the indices */
9498 ga2la_clear(dd->ga2la);
9500 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9503 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9505 /* Setup up the communication and communicate the coordinates */
9506 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9508 /* Set the indices */
9509 make_dd_indices(dd, cgs_gl->index, cg0);
9511 /* Set the charge group boundaries for neighbor searching */
9512 set_cg_boundaries(&comm->zones);
9514 if (fr->cutoff_scheme == ecutsVERLET)
9516 set_zones_size(dd, state_local->box, &ddbox,
9517 bSortCG ? 1 : 0, comm->zones.n);
9520 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9523 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9524 -1,state_local->x,state_local->box);
9527 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9529 /* Extract a local topology from the global topology */
9530 for (i = 0; i < dd->ndim; i++)
9532 np[dd->dim[i]] = comm->cd[i].np;
9534 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9535 comm->cellsize_min, np,
9537 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9538 vsite, top_global, top_local);
9540 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9542 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9544 /* Set up the special atom communication */
9545 n = comm->nat[ddnatZONE];
9546 for (i = ddnatZONE+1; i < ddnatNR; i++)
9551 if (vsite && vsite->n_intercg_vsite)
9553 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9557 if (dd->bInterCGcons || dd->bInterCGsettles)
9559 /* Only for inter-cg constraints we need special code */
9560 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9561 constr, ir->nProjOrder,
9562 top_local->idef.il);
9566 gmx_incons("Unknown special atom type setup");
9571 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9573 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9575 /* Make space for the extra coordinates for virtual site
9576 * or constraint communication.
9578 state_local->natoms = comm->nat[ddnatNR-1];
9579 if (state_local->natoms > state_local->nalloc)
9581 dd_realloc_state(state_local, f, state_local->natoms);
9584 if (fr->bF_NoVirSum)
9586 if (vsite && vsite->n_intercg_vsite)
9588 nat_f_novirsum = comm->nat[ddnatVSITE];
9592 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9594 nat_f_novirsum = dd->nat_tot;
9598 nat_f_novirsum = dd->nat_home;
9607 /* Set the number of atoms required for the force calculation.
9608 * Forces need to be constrained when using a twin-range setup
9609 * or with energy minimization. For simple simulations we could
9610 * avoid some allocation, zeroing and copying, but this is
9611 * probably not worth the complications ande checking.
9613 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9614 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9616 /* We make the all mdatoms up to nat_tot_con.
9617 * We could save some work by only setting invmass
9618 * between nat_tot and nat_tot_con.
9620 /* This call also sets the new number of home particles to dd->nat_home */
9621 atoms2md(top_global, ir,
9622 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9624 /* Now we have the charges we can sort the FE interactions */
9625 dd_sort_local_top(dd, mdatoms, top_local);
9629 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9630 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9635 /* Make the local shell stuff, currently no communication is done */
9636 make_local_shells(cr, mdatoms, shellfc);
9639 if (ir->implicit_solvent)
9641 make_local_gb(cr, fr->born, ir->gb_algorithm);
9644 init_bonded_thread_force_reduction(fr, &top_local->idef);
9646 if (!(cr->duty & DUTY_PME))
9648 /* Send the charges to our PME only node */
9649 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9650 mdatoms->chargeA, mdatoms->chargeB,
9651 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9656 set_constraints(constr, top_local, ir, mdatoms, cr);
9659 if (ir->ePull != epullNO)
9661 /* Update the local pull groups */
9662 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9667 /* Update the local rotation groups */
9668 dd_make_local_rotation_groups(dd, ir->rot);
9672 add_dd_statistics(dd);
9674 /* Make sure we only count the cycles for this DD partitioning */
9675 clear_dd_cycle_counts(dd);
9677 /* Because the order of the atoms might have changed since
9678 * the last vsite construction, we need to communicate the constructing
9679 * atom coordinates again (for spreading the forces this MD step).
9681 dd_move_x_vsites(dd, state_local->box, state_local->x);
9683 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9685 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9687 dd_move_x(dd, state_local->box, state_local->x);
9688 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9689 -1, state_local->x, state_local->box);
9692 /* Store the partitioning step */
9693 comm->partition_step = step;
9695 /* Increase the DD partitioning counter */
9697 /* The state currently matches this DD partitioning count, store it */
9698 state_local->ddp_count = dd->ddp_count;
9701 /* The DD master node knows the complete cg distribution,
9702 * store the count so we can possibly skip the cg info communication.
9704 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9707 if (comm->DD_debug > 0)
9709 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9710 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9711 "after partitioning");