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 && comm->bPMELoadBalDLBLimits)
2726 cellsize_min = max(cellsize_min,
2727 comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2730 return cellsize_min;
2733 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2736 real grid_jump_limit;
2738 /* The distance between the boundaries of cells at distance
2739 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2740 * and by the fact that cells should not be shifted by more than
2741 * half their size, such that cg's only shift by one cell
2742 * at redecomposition.
2744 grid_jump_limit = comm->cellsize_limit;
2745 if (!comm->bVacDLBNoLimit)
2747 if (comm->bPMELoadBalDLBLimits)
2749 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2751 grid_jump_limit = max(grid_jump_limit,
2752 cutoff/comm->cd[dim_ind].np);
2755 return grid_jump_limit;
2758 static gmx_bool check_grid_jump(gmx_large_int_t step,
2764 gmx_domdec_comm_t *comm;
2773 for (d = 1; d < dd->ndim; d++)
2776 limit = grid_jump_limit(comm, cutoff, d);
2777 bfac = ddbox->box_size[dim];
2778 if (ddbox->tric_dir[dim])
2780 bfac *= ddbox->skew_fac[dim];
2782 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2783 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2791 /* This error should never be triggered under normal
2792 * circumstances, but you never know ...
2794 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.",
2795 gmx_step_str(step, buf),
2796 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2804 static int dd_load_count(gmx_domdec_comm_t *comm)
2806 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2809 static float dd_force_load(gmx_domdec_comm_t *comm)
2816 if (comm->eFlop > 1)
2818 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2823 load = comm->cycl[ddCyclF];
2824 if (comm->cycl_n[ddCyclF] > 1)
2826 /* Subtract the maximum of the last n cycle counts
2827 * to get rid of possible high counts due to other soures,
2828 * for instance system activity, that would otherwise
2829 * affect the dynamic load balancing.
2831 load -= comm->cycl_max[ddCyclF];
2838 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2840 gmx_domdec_comm_t *comm;
2845 snew(*dim_f, dd->nc[dim]+1);
2847 for (i = 1; i < dd->nc[dim]; i++)
2849 if (comm->slb_frac[dim])
2851 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2855 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2858 (*dim_f)[dd->nc[dim]] = 1;
2861 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2863 int pmeindex, slab, nso, i;
2866 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2872 ddpme->dim = dimind;
2874 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2876 ddpme->nslab = (ddpme->dim == 0 ?
2877 dd->comm->npmenodes_x :
2878 dd->comm->npmenodes_y);
2880 if (ddpme->nslab <= 1)
2885 nso = dd->comm->npmenodes/ddpme->nslab;
2886 /* Determine for each PME slab the PP location range for dimension dim */
2887 snew(ddpme->pp_min, ddpme->nslab);
2888 snew(ddpme->pp_max, ddpme->nslab);
2889 for (slab = 0; slab < ddpme->nslab; slab++)
2891 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2892 ddpme->pp_max[slab] = 0;
2894 for (i = 0; i < dd->nnodes; i++)
2896 ddindex2xyz(dd->nc, i, xyz);
2897 /* For y only use our y/z slab.
2898 * This assumes that the PME x grid size matches the DD grid size.
2900 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2902 pmeindex = ddindex2pmeindex(dd, i);
2905 slab = pmeindex/nso;
2909 slab = pmeindex % ddpme->nslab;
2911 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2912 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2916 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2919 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2921 if (dd->comm->ddpme[0].dim == XX)
2923 return dd->comm->ddpme[0].maxshift;
2931 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2933 if (dd->comm->ddpme[0].dim == YY)
2935 return dd->comm->ddpme[0].maxshift;
2937 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2939 return dd->comm->ddpme[1].maxshift;
2947 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2948 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2950 gmx_domdec_comm_t *comm;
2953 real range, pme_boundary;
2957 nc = dd->nc[ddpme->dim];
2960 if (!ddpme->dim_match)
2962 /* PP decomposition is not along dim: the worst situation */
2965 else if (ns <= 3 || (bUniform && ns == nc))
2967 /* The optimal situation */
2972 /* We need to check for all pme nodes which nodes they
2973 * could possibly need to communicate with.
2975 xmin = ddpme->pp_min;
2976 xmax = ddpme->pp_max;
2977 /* Allow for atoms to be maximally 2/3 times the cut-off
2978 * out of their DD cell. This is a reasonable balance between
2979 * between performance and support for most charge-group/cut-off
2982 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2983 /* Avoid extra communication when we are exactly at a boundary */
2987 for (s = 0; s < ns; s++)
2989 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2990 pme_boundary = (real)s/ns;
2993 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
2995 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
2999 pme_boundary = (real)(s+1)/ns;
3002 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3004 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3011 ddpme->maxshift = sh;
3015 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3016 ddpme->dim, ddpme->maxshift);
3020 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3024 for (d = 0; d < dd->ndim; d++)
3027 if (dim < ddbox->nboundeddim &&
3028 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3029 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3031 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",
3032 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3033 dd->nc[dim], dd->comm->cellsize_limit);
3038 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3039 gmx_bool bMaster, ivec npulse)
3041 gmx_domdec_comm_t *comm;
3044 real *cell_x, cell_dx, cellsize;
3048 for (d = 0; d < DIM; d++)
3050 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3052 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3055 cell_dx = ddbox->box_size[d]/dd->nc[d];
3058 for (j = 0; j < dd->nc[d]+1; j++)
3060 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3065 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3066 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3068 cellsize = cell_dx*ddbox->skew_fac[d];
3069 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3073 cellsize_min[d] = cellsize;
3077 /* Statically load balanced grid */
3078 /* Also when we are not doing a master distribution we determine
3079 * all cell borders in a loop to obtain identical values
3080 * to the master distribution case and to determine npulse.
3084 cell_x = dd->ma->cell_x[d];
3088 snew(cell_x, dd->nc[d]+1);
3090 cell_x[0] = ddbox->box0[d];
3091 for (j = 0; j < dd->nc[d]; j++)
3093 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3094 cell_x[j+1] = cell_x[j] + cell_dx;
3095 cellsize = cell_dx*ddbox->skew_fac[d];
3096 while (cellsize*npulse[d] < comm->cutoff &&
3097 npulse[d] < dd->nc[d]-1)
3101 cellsize_min[d] = min(cellsize_min[d], cellsize);
3105 comm->cell_x0[d] = cell_x[dd->ci[d]];
3106 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3110 /* The following limitation is to avoid that a cell would receive
3111 * some of its own home charge groups back over the periodic boundary.
3112 * Double charge groups cause trouble with the global indices.
3114 if (d < ddbox->npbcdim &&
3115 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3117 gmx_fatal_collective(FARGS, NULL, dd,
3118 "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",
3119 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3121 dd->nc[d], dd->nc[d],
3122 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3126 if (!comm->bDynLoadBal)
3128 copy_rvec(cellsize_min, comm->cellsize_min);
3131 for (d = 0; d < comm->npmedecompdim; d++)
3133 set_pme_maxshift(dd, &comm->ddpme[d],
3134 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3135 comm->ddpme[d].slb_dim_f);
3140 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3141 int d, int dim, gmx_domdec_root_t *root,
3143 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3145 gmx_domdec_comm_t *comm;
3146 int ncd, i, j, nmin, nmin_old;
3147 gmx_bool bLimLo, bLimHi;
3149 real fac, halfway, cellsize_limit_f_i, region_size;
3150 gmx_bool bPBC, bLastHi = FALSE;
3151 int nrange[] = {range[0], range[1]};
3153 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3159 bPBC = (dim < ddbox->npbcdim);
3161 cell_size = root->buf_ncd;
3165 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3168 /* First we need to check if the scaling does not make cells
3169 * smaller than the smallest allowed size.
3170 * We need to do this iteratively, since if a cell is too small,
3171 * it needs to be enlarged, which makes all the other cells smaller,
3172 * which could in turn make another cell smaller than allowed.
3174 for (i = range[0]; i < range[1]; i++)
3176 root->bCellMin[i] = FALSE;
3182 /* We need the total for normalization */
3184 for (i = range[0]; i < range[1]; i++)
3186 if (root->bCellMin[i] == FALSE)
3188 fac += cell_size[i];
3191 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3192 /* Determine the cell boundaries */
3193 for (i = range[0]; i < range[1]; i++)
3195 if (root->bCellMin[i] == FALSE)
3197 cell_size[i] *= fac;
3198 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3200 cellsize_limit_f_i = 0;
3204 cellsize_limit_f_i = cellsize_limit_f;
3206 if (cell_size[i] < cellsize_limit_f_i)
3208 root->bCellMin[i] = TRUE;
3209 cell_size[i] = cellsize_limit_f_i;
3213 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3216 while (nmin > nmin_old);
3219 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3220 /* For this check we should not use DD_CELL_MARGIN,
3221 * but a slightly smaller factor,
3222 * since rounding could get use below the limit.
3224 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3227 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",
3228 gmx_step_str(step, buf),
3229 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3230 ncd, comm->cellsize_min[dim]);
3233 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3237 /* Check if the boundary did not displace more than halfway
3238 * each of the cells it bounds, as this could cause problems,
3239 * especially when the differences between cell sizes are large.
3240 * If changes are applied, they will not make cells smaller
3241 * than the cut-off, as we check all the boundaries which
3242 * might be affected by a change and if the old state was ok,
3243 * the cells will at most be shrunk back to their old size.
3245 for (i = range[0]+1; i < range[1]; i++)
3247 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3248 if (root->cell_f[i] < halfway)
3250 root->cell_f[i] = halfway;
3251 /* Check if the change also causes shifts of the next boundaries */
3252 for (j = i+1; j < range[1]; j++)
3254 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3256 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3260 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3261 if (root->cell_f[i] > halfway)
3263 root->cell_f[i] = halfway;
3264 /* Check if the change also causes shifts of the next boundaries */
3265 for (j = i-1; j >= range[0]+1; j--)
3267 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3269 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3276 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3277 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3278 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3279 * for a and b nrange is used */
3282 /* Take care of the staggering of the cell boundaries */
3285 for (i = range[0]; i < range[1]; i++)
3287 root->cell_f_max0[i] = root->cell_f[i];
3288 root->cell_f_min1[i] = root->cell_f[i+1];
3293 for (i = range[0]+1; i < range[1]; i++)
3295 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3296 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3297 if (bLimLo && bLimHi)
3299 /* Both limits violated, try the best we can */
3300 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3301 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3302 nrange[0] = range[0];
3304 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3307 nrange[1] = range[1];
3308 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3314 /* root->cell_f[i] = root->bound_min[i]; */
3315 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3318 else if (bLimHi && !bLastHi)
3321 if (nrange[1] < range[1]) /* found a LimLo before */
3323 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3324 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3325 nrange[0] = nrange[1];
3327 root->cell_f[i] = root->bound_max[i];
3329 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3331 nrange[1] = range[1];
3334 if (nrange[1] < range[1]) /* found last a LimLo */
3336 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3337 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3338 nrange[0] = nrange[1];
3339 nrange[1] = range[1];
3340 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3342 else if (nrange[0] > range[0]) /* found at least one LimHi */
3344 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3351 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3352 int d, int dim, gmx_domdec_root_t *root,
3353 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3354 gmx_bool bUniform, gmx_large_int_t step)
3356 gmx_domdec_comm_t *comm;
3357 int ncd, d1, i, j, pos;
3359 real load_aver, load_i, imbalance, change, change_max, sc;
3360 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3364 int range[] = { 0, 0 };
3368 /* Convert the maximum change from the input percentage to a fraction */
3369 change_limit = comm->dlb_scale_lim*0.01;
3373 bPBC = (dim < ddbox->npbcdim);
3375 cell_size = root->buf_ncd;
3377 /* Store the original boundaries */
3378 for (i = 0; i < ncd+1; i++)
3380 root->old_cell_f[i] = root->cell_f[i];
3384 for (i = 0; i < ncd; i++)
3386 cell_size[i] = 1.0/ncd;
3389 else if (dd_load_count(comm))
3391 load_aver = comm->load[d].sum_m/ncd;
3393 for (i = 0; i < ncd; i++)
3395 /* Determine the relative imbalance of cell i */
3396 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3397 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3398 /* Determine the change of the cell size using underrelaxation */
3399 change = -relax*imbalance;
3400 change_max = max(change_max, max(change, -change));
3402 /* Limit the amount of scaling.
3403 * We need to use the same rescaling for all cells in one row,
3404 * otherwise the load balancing might not converge.
3407 if (change_max > change_limit)
3409 sc *= change_limit/change_max;
3411 for (i = 0; i < ncd; i++)
3413 /* Determine the relative imbalance of cell i */
3414 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3415 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3416 /* Determine the change of the cell size using underrelaxation */
3417 change = -sc*imbalance;
3418 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3422 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3423 cellsize_limit_f *= DD_CELL_MARGIN;
3424 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3425 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3426 if (ddbox->tric_dir[dim])
3428 cellsize_limit_f /= ddbox->skew_fac[dim];
3429 dist_min_f /= ddbox->skew_fac[dim];
3431 if (bDynamicBox && d > 0)
3433 dist_min_f *= DD_PRES_SCALE_MARGIN;
3435 if (d > 0 && !bUniform)
3437 /* Make sure that the grid is not shifted too much */
3438 for (i = 1; i < ncd; i++)
3440 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3442 gmx_incons("Inconsistent DD boundary staggering limits!");
3444 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3445 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3448 root->bound_min[i] += 0.5*space;
3450 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3451 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3454 root->bound_max[i] += 0.5*space;
3459 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3461 root->cell_f_max0[i-1] + dist_min_f,
3462 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3463 root->cell_f_min1[i] - dist_min_f);
3468 root->cell_f[0] = 0;
3469 root->cell_f[ncd] = 1;
3470 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3473 /* After the checks above, the cells should obey the cut-off
3474 * restrictions, but it does not hurt to check.
3476 for (i = 0; i < ncd; i++)
3480 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3481 dim, i, root->cell_f[i], root->cell_f[i+1]);
3484 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3485 root->cell_f[i+1] - root->cell_f[i] <
3486 cellsize_limit_f/DD_CELL_MARGIN)
3490 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3491 gmx_step_str(step, buf), dim2char(dim), i,
3492 (root->cell_f[i+1] - root->cell_f[i])
3493 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3498 /* Store the cell boundaries of the lower dimensions at the end */
3499 for (d1 = 0; d1 < d; d1++)
3501 root->cell_f[pos++] = comm->cell_f0[d1];
3502 root->cell_f[pos++] = comm->cell_f1[d1];
3505 if (d < comm->npmedecompdim)
3507 /* The master determines the maximum shift for
3508 * the coordinate communication between separate PME nodes.
3510 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3512 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3515 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3519 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3520 gmx_ddbox_t *ddbox, int dimind)
3522 gmx_domdec_comm_t *comm;
3527 /* Set the cell dimensions */
3528 dim = dd->dim[dimind];
3529 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3530 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3531 if (dim >= ddbox->nboundeddim)
3533 comm->cell_x0[dim] += ddbox->box0[dim];
3534 comm->cell_x1[dim] += ddbox->box0[dim];
3538 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3539 int d, int dim, real *cell_f_row,
3542 gmx_domdec_comm_t *comm;
3548 /* Each node would only need to know two fractions,
3549 * but it is probably cheaper to broadcast the whole array.
3551 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3552 0, comm->mpi_comm_load[d]);
3554 /* Copy the fractions for this dimension from the buffer */
3555 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3556 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3557 /* The whole array was communicated, so set the buffer position */
3558 pos = dd->nc[dim] + 1;
3559 for (d1 = 0; d1 <= d; d1++)
3563 /* Copy the cell fractions of the lower dimensions */
3564 comm->cell_f0[d1] = cell_f_row[pos++];
3565 comm->cell_f1[d1] = cell_f_row[pos++];
3567 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3569 /* Convert the communicated shift from float to int */
3570 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3573 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3577 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3578 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3579 gmx_bool bUniform, gmx_large_int_t step)
3581 gmx_domdec_comm_t *comm;
3583 gmx_bool bRowMember, bRowRoot;
3588 for (d = 0; d < dd->ndim; d++)
3593 for (d1 = d; d1 < dd->ndim; d1++)
3595 if (dd->ci[dd->dim[d1]] > 0)
3608 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3609 ddbox, bDynamicBox, bUniform, step);
3610 cell_f_row = comm->root[d]->cell_f;
3614 cell_f_row = comm->cell_f_row;
3616 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3621 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3625 /* This function assumes the box is static and should therefore
3626 * not be called when the box has changed since the last
3627 * call to dd_partition_system.
3629 for (d = 0; d < dd->ndim; d++)
3631 relative_to_absolute_cell_bounds(dd, ddbox, d);
3637 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3638 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3639 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3640 gmx_wallcycle_t wcycle)
3642 gmx_domdec_comm_t *comm;
3649 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3650 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3651 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3653 else if (bDynamicBox)
3655 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3658 /* Set the dimensions for which no DD is used */
3659 for (dim = 0; dim < DIM; dim++)
3661 if (dd->nc[dim] == 1)
3663 comm->cell_x0[dim] = 0;
3664 comm->cell_x1[dim] = ddbox->box_size[dim];
3665 if (dim >= ddbox->nboundeddim)
3667 comm->cell_x0[dim] += ddbox->box0[dim];
3668 comm->cell_x1[dim] += ddbox->box0[dim];
3674 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3677 gmx_domdec_comm_dim_t *cd;
3679 for (d = 0; d < dd->ndim; d++)
3681 cd = &dd->comm->cd[d];
3682 np = npulse[dd->dim[d]];
3683 if (np > cd->np_nalloc)
3687 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3688 dim2char(dd->dim[d]), np);
3690 if (DDMASTER(dd) && cd->np_nalloc > 0)
3692 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3694 srenew(cd->ind, np);
3695 for (i = cd->np_nalloc; i < np; i++)
3697 cd->ind[i].index = NULL;
3698 cd->ind[i].nalloc = 0;
3707 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3708 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3709 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3710 gmx_wallcycle_t wcycle)
3712 gmx_domdec_comm_t *comm;
3718 /* Copy the old cell boundaries for the cg displacement check */
3719 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3720 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3722 if (comm->bDynLoadBal)
3726 check_box_size(dd, ddbox);
3728 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3732 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3733 realloc_comm_ind(dd, npulse);
3738 for (d = 0; d < DIM; d++)
3740 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3741 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3746 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3748 rvec cell_ns_x0, rvec cell_ns_x1,
3749 gmx_large_int_t step)
3751 gmx_domdec_comm_t *comm;
3756 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3758 dim = dd->dim[dim_ind];
3760 /* Without PBC we don't have restrictions on the outer cells */
3761 if (!(dim >= ddbox->npbcdim &&
3762 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3763 comm->bDynLoadBal &&
3764 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3765 comm->cellsize_min[dim])
3768 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",
3769 gmx_step_str(step, buf), dim2char(dim),
3770 comm->cell_x1[dim] - comm->cell_x0[dim],
3771 ddbox->skew_fac[dim],
3772 dd->comm->cellsize_min[dim],
3773 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3777 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3779 /* Communicate the boundaries and update cell_ns_x0/1 */
3780 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3781 if (dd->bGridJump && dd->ndim > 1)
3783 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3788 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3792 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3800 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3801 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3810 static void check_screw_box(matrix box)
3812 /* Mathematical limitation */
3813 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3815 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3818 /* Limitation due to the asymmetry of the eighth shell method */
3819 if (box[ZZ][YY] != 0)
3821 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3825 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3826 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3829 gmx_domdec_master_t *ma;
3830 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3831 int i, icg, j, k, k0, k1, d, npbcdim;
3833 rvec box_size, cg_cm;
3835 real nrcg, inv_ncg, pos_d;
3837 gmx_bool bUnbounded, bScrew;
3841 if (tmp_ind == NULL)
3843 snew(tmp_nalloc, dd->nnodes);
3844 snew(tmp_ind, dd->nnodes);
3845 for (i = 0; i < dd->nnodes; i++)
3847 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3848 snew(tmp_ind[i], tmp_nalloc[i]);
3852 /* Clear the count */
3853 for (i = 0; i < dd->nnodes; i++)
3859 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3861 cgindex = cgs->index;
3863 /* Compute the center of geometry for all charge groups */
3864 for (icg = 0; icg < cgs->nr; icg++)
3867 k1 = cgindex[icg+1];
3871 copy_rvec(pos[k0], cg_cm);
3878 for (k = k0; (k < k1); k++)
3880 rvec_inc(cg_cm, pos[k]);
3882 for (d = 0; (d < DIM); d++)
3884 cg_cm[d] *= inv_ncg;
3887 /* Put the charge group in the box and determine the cell index */
3888 for (d = DIM-1; d >= 0; d--)
3891 if (d < dd->npbcdim)
3893 bScrew = (dd->bScrewPBC && d == XX);
3894 if (tric_dir[d] && dd->nc[d] > 1)
3896 /* Use triclinic coordintates for this dimension */
3897 for (j = d+1; j < DIM; j++)
3899 pos_d += cg_cm[j]*tcm[j][d];
3902 while (pos_d >= box[d][d])
3905 rvec_dec(cg_cm, box[d]);
3908 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3909 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3911 for (k = k0; (k < k1); k++)
3913 rvec_dec(pos[k], box[d]);
3916 pos[k][YY] = box[YY][YY] - pos[k][YY];
3917 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3924 rvec_inc(cg_cm, box[d]);
3927 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3928 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3930 for (k = k0; (k < k1); k++)
3932 rvec_inc(pos[k], box[d]);
3935 pos[k][YY] = box[YY][YY] - pos[k][YY];
3936 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3941 /* This could be done more efficiently */
3943 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3948 i = dd_index(dd->nc, ind);
3949 if (ma->ncg[i] == tmp_nalloc[i])
3951 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3952 srenew(tmp_ind[i], tmp_nalloc[i]);
3954 tmp_ind[i][ma->ncg[i]] = icg;
3956 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3960 for (i = 0; i < dd->nnodes; i++)
3963 for (k = 0; k < ma->ncg[i]; k++)
3965 ma->cg[k1++] = tmp_ind[i][k];
3968 ma->index[dd->nnodes] = k1;
3970 for (i = 0; i < dd->nnodes; i++)
3980 fprintf(fplog, "Charge group distribution at step %s:",
3981 gmx_step_str(step, buf));
3982 for (i = 0; i < dd->nnodes; i++)
3984 fprintf(fplog, " %d", ma->ncg[i]);
3986 fprintf(fplog, "\n");
3990 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
3991 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
3994 gmx_domdec_master_t *ma = NULL;
3997 int *ibuf, buf2[2] = { 0, 0 };
3998 gmx_bool bMaster = DDMASTER(dd);
4005 check_screw_box(box);
4008 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4010 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4011 for (i = 0; i < dd->nnodes; i++)
4013 ma->ibuf[2*i] = ma->ncg[i];
4014 ma->ibuf[2*i+1] = ma->nat[i];
4022 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4024 dd->ncg_home = buf2[0];
4025 dd->nat_home = buf2[1];
4026 dd->ncg_tot = dd->ncg_home;
4027 dd->nat_tot = dd->nat_home;
4028 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4030 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4031 srenew(dd->index_gl, dd->cg_nalloc);
4032 srenew(dd->cgindex, dd->cg_nalloc+1);
4036 for (i = 0; i < dd->nnodes; i++)
4038 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4039 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4044 DDMASTER(dd) ? ma->ibuf : NULL,
4045 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4046 DDMASTER(dd) ? ma->cg : NULL,
4047 dd->ncg_home*sizeof(int), dd->index_gl);
4049 /* Determine the home charge group sizes */
4051 for (i = 0; i < dd->ncg_home; i++)
4053 cg_gl = dd->index_gl[i];
4055 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4060 fprintf(debug, "Home charge groups:\n");
4061 for (i = 0; i < dd->ncg_home; i++)
4063 fprintf(debug, " %d", dd->index_gl[i]);
4066 fprintf(debug, "\n");
4069 fprintf(debug, "\n");
4073 static int compact_and_copy_vec_at(int ncg, int *move,
4076 rvec *src, gmx_domdec_comm_t *comm,
4079 int m, icg, i, i0, i1, nrcg;
4085 for (m = 0; m < DIM*2; m++)
4091 for (icg = 0; icg < ncg; icg++)
4093 i1 = cgindex[icg+1];
4099 /* Compact the home array in place */
4100 for (i = i0; i < i1; i++)
4102 copy_rvec(src[i], src[home_pos++]);
4108 /* Copy to the communication buffer */
4110 pos_vec[m] += 1 + vec*nrcg;
4111 for (i = i0; i < i1; i++)
4113 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4115 pos_vec[m] += (nvec - vec - 1)*nrcg;
4119 home_pos += i1 - i0;
4127 static int compact_and_copy_vec_cg(int ncg, int *move,
4129 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4132 int m, icg, i0, i1, nrcg;
4138 for (m = 0; m < DIM*2; m++)
4144 for (icg = 0; icg < ncg; icg++)
4146 i1 = cgindex[icg+1];
4152 /* Compact the home array in place */
4153 copy_rvec(src[icg], src[home_pos++]);
4159 /* Copy to the communication buffer */
4160 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4161 pos_vec[m] += 1 + nrcg*nvec;
4173 static int compact_ind(int ncg, int *move,
4174 int *index_gl, int *cgindex,
4176 gmx_ga2la_t ga2la, char *bLocalCG,
4179 int cg, nat, a0, a1, a, a_gl;
4184 for (cg = 0; cg < ncg; cg++)
4190 /* Compact the home arrays in place.
4191 * Anything that can be done here avoids access to global arrays.
4193 cgindex[home_pos] = nat;
4194 for (a = a0; a < a1; a++)
4197 gatindex[nat] = a_gl;
4198 /* The cell number stays 0, so we don't need to set it */
4199 ga2la_change_la(ga2la, a_gl, nat);
4202 index_gl[home_pos] = index_gl[cg];
4203 cginfo[home_pos] = cginfo[cg];
4204 /* The charge group remains local, so bLocalCG does not change */
4209 /* Clear the global indices */
4210 for (a = a0; a < a1; a++)
4212 ga2la_del(ga2la, gatindex[a]);
4216 bLocalCG[index_gl[cg]] = FALSE;
4220 cgindex[home_pos] = nat;
4225 static void clear_and_mark_ind(int ncg, int *move,
4226 int *index_gl, int *cgindex, int *gatindex,
4227 gmx_ga2la_t ga2la, char *bLocalCG,
4232 for (cg = 0; cg < ncg; cg++)
4238 /* Clear the global indices */
4239 for (a = a0; a < a1; a++)
4241 ga2la_del(ga2la, gatindex[a]);
4245 bLocalCG[index_gl[cg]] = FALSE;
4247 /* Signal that this cg has moved using the ns cell index.
4248 * Here we set it to -1. fill_grid will change it
4249 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4251 cell_index[cg] = -1;
4256 static void print_cg_move(FILE *fplog,
4258 gmx_large_int_t step, int cg, int dim, int dir,
4259 gmx_bool bHaveLimitdAndCMOld, real limitd,
4260 rvec cm_old, rvec cm_new, real pos_d)
4262 gmx_domdec_comm_t *comm;
4267 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4268 if (bHaveLimitdAndCMOld)
4270 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4271 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4275 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4276 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4278 fprintf(fplog, "distance out of cell %f\n",
4279 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4280 if (bHaveLimitdAndCMOld)
4282 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4283 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4285 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4286 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4287 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4289 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4290 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4292 comm->cell_x0[dim], comm->cell_x1[dim]);
4295 static void cg_move_error(FILE *fplog,
4297 gmx_large_int_t step, int cg, int dim, int dir,
4298 gmx_bool bHaveLimitdAndCMOld, real limitd,
4299 rvec cm_old, rvec cm_new, real pos_d)
4303 print_cg_move(fplog, dd, step, cg, dim, dir,
4304 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4306 print_cg_move(stderr, dd, step, cg, dim, dir,
4307 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4309 "A charge group moved too far between two domain decomposition steps\n"
4310 "This usually means that your system is not well equilibrated");
4313 static void rotate_state_atom(t_state *state, int a)
4317 for (est = 0; est < estNR; est++)
4319 if (EST_DISTR(est) && (state->flags & (1<<est)))
4324 /* Rotate the complete state; for a rectangular box only */
4325 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4326 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4329 state->v[a][YY] = -state->v[a][YY];
4330 state->v[a][ZZ] = -state->v[a][ZZ];
4333 state->sd_X[a][YY] = -state->sd_X[a][YY];
4334 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4337 state->cg_p[a][YY] = -state->cg_p[a][YY];
4338 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4340 case estDISRE_INITF:
4341 case estDISRE_RM3TAV:
4342 case estORIRE_INITF:
4344 /* These are distances, so not affected by rotation */
4347 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4353 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4355 if (natoms > comm->moved_nalloc)
4357 /* Contents should be preserved here */
4358 comm->moved_nalloc = over_alloc_dd(natoms);
4359 srenew(comm->moved, comm->moved_nalloc);
4365 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4368 ivec tric_dir, matrix tcm,
4369 rvec cell_x0, rvec cell_x1,
4370 rvec limitd, rvec limit0, rvec limit1,
4372 int cg_start, int cg_end,
4377 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4378 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4382 real inv_ncg, pos_d;
4385 npbcdim = dd->npbcdim;
4387 for (cg = cg_start; cg < cg_end; cg++)
4394 copy_rvec(state->x[k0], cm_new);
4401 for (k = k0; (k < k1); k++)
4403 rvec_inc(cm_new, state->x[k]);
4405 for (d = 0; (d < DIM); d++)
4407 cm_new[d] = inv_ncg*cm_new[d];
4412 /* Do pbc and check DD cell boundary crossings */
4413 for (d = DIM-1; d >= 0; d--)
4417 bScrew = (dd->bScrewPBC && d == XX);
4418 /* Determine the location of this cg in lattice coordinates */
4422 for (d2 = d+1; d2 < DIM; d2++)
4424 pos_d += cm_new[d2]*tcm[d2][d];
4427 /* Put the charge group in the triclinic unit-cell */
4428 if (pos_d >= cell_x1[d])
4430 if (pos_d >= limit1[d])
4432 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4433 cg_cm[cg], cm_new, pos_d);
4436 if (dd->ci[d] == dd->nc[d] - 1)
4438 rvec_dec(cm_new, state->box[d]);
4441 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4442 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4444 for (k = k0; (k < k1); k++)
4446 rvec_dec(state->x[k], state->box[d]);
4449 rotate_state_atom(state, k);
4454 else if (pos_d < cell_x0[d])
4456 if (pos_d < limit0[d])
4458 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4459 cg_cm[cg], cm_new, pos_d);
4464 rvec_inc(cm_new, state->box[d]);
4467 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4468 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4470 for (k = k0; (k < k1); k++)
4472 rvec_inc(state->x[k], state->box[d]);
4475 rotate_state_atom(state, k);
4481 else if (d < npbcdim)
4483 /* Put the charge group in the rectangular unit-cell */
4484 while (cm_new[d] >= state->box[d][d])
4486 rvec_dec(cm_new, state->box[d]);
4487 for (k = k0; (k < k1); k++)
4489 rvec_dec(state->x[k], state->box[d]);
4492 while (cm_new[d] < 0)
4494 rvec_inc(cm_new, state->box[d]);
4495 for (k = k0; (k < k1); k++)
4497 rvec_inc(state->x[k], state->box[d]);
4503 copy_rvec(cm_new, cg_cm[cg]);
4505 /* Determine where this cg should go */
4508 for (d = 0; d < dd->ndim; d++)
4513 flag |= DD_FLAG_FW(d);
4519 else if (dev[dim] == -1)
4521 flag |= DD_FLAG_BW(d);
4524 if (dd->nc[dim] > 2)
4535 /* Temporarily store the flag in move */
4536 move[cg] = mc + flag;
4540 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4541 gmx_domdec_t *dd, ivec tric_dir,
4542 t_state *state, rvec **f,
4543 t_forcerec *fr, t_mdatoms *md,
4551 int ncg[DIM*2], nat[DIM*2];
4552 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4553 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4554 int sbuf[2], rbuf[2];
4555 int home_pos_cg, home_pos_at, buf_pos;
4557 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4560 real inv_ncg, pos_d;
4562 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4564 cginfo_mb_t *cginfo_mb;
4565 gmx_domdec_comm_t *comm;
4567 int nthread, thread;
4571 check_screw_box(state->box);
4575 if (fr->cutoff_scheme == ecutsGROUP)
4580 for (i = 0; i < estNR; i++)
4586 case estX: /* Always present */ break;
4587 case estV: bV = (state->flags & (1<<i)); break;
4588 case estSDX: bSDX = (state->flags & (1<<i)); break;
4589 case estCGP: bCGP = (state->flags & (1<<i)); break;
4592 case estDISRE_INITF:
4593 case estDISRE_RM3TAV:
4594 case estORIRE_INITF:
4596 /* No processing required */
4599 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4604 if (dd->ncg_tot > comm->nalloc_int)
4606 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4607 srenew(comm->buf_int, comm->nalloc_int);
4609 move = comm->buf_int;
4611 /* Clear the count */
4612 for (c = 0; c < dd->ndim*2; c++)
4618 npbcdim = dd->npbcdim;
4620 for (d = 0; (d < DIM); d++)
4622 limitd[d] = dd->comm->cellsize_min[d];
4623 if (d >= npbcdim && dd->ci[d] == 0)
4625 cell_x0[d] = -GMX_FLOAT_MAX;
4629 cell_x0[d] = comm->cell_x0[d];
4631 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4633 cell_x1[d] = GMX_FLOAT_MAX;
4637 cell_x1[d] = comm->cell_x1[d];
4641 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4642 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4646 /* We check after communication if a charge group moved
4647 * more than one cell. Set the pre-comm check limit to float_max.
4649 limit0[d] = -GMX_FLOAT_MAX;
4650 limit1[d] = GMX_FLOAT_MAX;
4654 make_tric_corr_matrix(npbcdim, state->box, tcm);
4656 cgindex = dd->cgindex;
4658 nthread = gmx_omp_nthreads_get(emntDomdec);
4660 /* Compute the center of geometry for all home charge groups
4661 * and put them in the box and determine where they should go.
4663 #pragma omp parallel for num_threads(nthread) schedule(static)
4664 for (thread = 0; thread < nthread; thread++)
4666 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4667 cell_x0, cell_x1, limitd, limit0, limit1,
4669 ( thread *dd->ncg_home)/nthread,
4670 ((thread+1)*dd->ncg_home)/nthread,
4671 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4675 for (cg = 0; cg < dd->ncg_home; cg++)
4680 flag = mc & ~DD_FLAG_NRCG;
4681 mc = mc & DD_FLAG_NRCG;
4684 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4686 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4687 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4689 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4690 /* We store the cg size in the lower 16 bits
4691 * and the place where the charge group should go
4692 * in the next 6 bits. This saves some communication volume.
4694 nrcg = cgindex[cg+1] - cgindex[cg];
4695 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4701 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4702 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4705 for (i = 0; i < dd->ndim*2; i++)
4707 *ncg_moved += ncg[i];
4724 /* Make sure the communication buffers are large enough */
4725 for (mc = 0; mc < dd->ndim*2; mc++)
4727 nvr = ncg[mc] + nat[mc]*nvec;
4728 if (nvr > comm->cgcm_state_nalloc[mc])
4730 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4731 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4735 switch (fr->cutoff_scheme)
4738 /* Recalculating cg_cm might be cheaper than communicating,
4739 * but that could give rise to rounding issues.
4742 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4743 nvec, cg_cm, comm, bCompact);
4746 /* Without charge groups we send the moved atom coordinates
4747 * over twice. This is so the code below can be used without
4748 * many conditionals for both for with and without charge groups.
4751 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4752 nvec, state->x, comm, FALSE);
4755 home_pos_cg -= *ncg_moved;
4759 gmx_incons("unimplemented");
4765 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4766 nvec, vec++, state->x, comm, bCompact);
4769 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4770 nvec, vec++, state->v, comm, bCompact);
4774 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4775 nvec, vec++, state->sd_X, comm, bCompact);
4779 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4780 nvec, vec++, state->cg_p, comm, bCompact);
4785 compact_ind(dd->ncg_home, move,
4786 dd->index_gl, dd->cgindex, dd->gatindex,
4787 dd->ga2la, comm->bLocalCG,
4792 if (fr->cutoff_scheme == ecutsVERLET)
4794 moved = get_moved(comm, dd->ncg_home);
4796 for (k = 0; k < dd->ncg_home; k++)
4803 moved = fr->ns.grid->cell_index;
4806 clear_and_mark_ind(dd->ncg_home, move,
4807 dd->index_gl, dd->cgindex, dd->gatindex,
4808 dd->ga2la, comm->bLocalCG,
4812 cginfo_mb = fr->cginfo_mb;
4814 *ncg_stay_home = home_pos_cg;
4815 for (d = 0; d < dd->ndim; d++)
4821 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4824 /* Communicate the cg and atom counts */
4829 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4830 d, dir, sbuf[0], sbuf[1]);
4832 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4834 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4836 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4837 srenew(comm->buf_int, comm->nalloc_int);
4840 /* Communicate the charge group indices, sizes and flags */
4841 dd_sendrecv_int(dd, d, dir,
4842 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4843 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4845 nvs = ncg[cdd] + nat[cdd]*nvec;
4846 i = rbuf[0] + rbuf[1] *nvec;
4847 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4849 /* Communicate cgcm and state */
4850 dd_sendrecv_rvec(dd, d, dir,
4851 comm->cgcm_state[cdd], nvs,
4852 comm->vbuf.v+nvr, i);
4853 ncg_recv += rbuf[0];
4854 nat_recv += rbuf[1];
4858 /* Process the received charge groups */
4860 for (cg = 0; cg < ncg_recv; cg++)
4862 flag = comm->buf_int[cg*DD_CGIBS+1];
4864 if (dim >= npbcdim && dd->nc[dim] > 2)
4866 /* No pbc in this dim and more than one domain boundary.
4867 * We do a separate check if a charge group didn't move too far.
4869 if (((flag & DD_FLAG_FW(d)) &&
4870 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4871 ((flag & DD_FLAG_BW(d)) &&
4872 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4874 cg_move_error(fplog, dd, step, cg, dim,
4875 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4877 comm->vbuf.v[buf_pos],
4878 comm->vbuf.v[buf_pos],
4879 comm->vbuf.v[buf_pos][dim]);
4886 /* Check which direction this cg should go */
4887 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4891 /* The cell boundaries for dimension d2 are not equal
4892 * for each cell row of the lower dimension(s),
4893 * therefore we might need to redetermine where
4894 * this cg should go.
4897 /* If this cg crosses the box boundary in dimension d2
4898 * we can use the communicated flag, so we do not
4899 * have to worry about pbc.
4901 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4902 (flag & DD_FLAG_FW(d2))) ||
4903 (dd->ci[dim2] == 0 &&
4904 (flag & DD_FLAG_BW(d2)))))
4906 /* Clear the two flags for this dimension */
4907 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4908 /* Determine the location of this cg
4909 * in lattice coordinates
4911 pos_d = comm->vbuf.v[buf_pos][dim2];
4914 for (d3 = dim2+1; d3 < DIM; d3++)
4917 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4920 /* Check of we are not at the box edge.
4921 * pbc is only handled in the first step above,
4922 * but this check could move over pbc while
4923 * the first step did not due to different rounding.
4925 if (pos_d >= cell_x1[dim2] &&
4926 dd->ci[dim2] != dd->nc[dim2]-1)
4928 flag |= DD_FLAG_FW(d2);
4930 else if (pos_d < cell_x0[dim2] &&
4933 flag |= DD_FLAG_BW(d2);
4935 comm->buf_int[cg*DD_CGIBS+1] = flag;
4938 /* Set to which neighboring cell this cg should go */
4939 if (flag & DD_FLAG_FW(d2))
4943 else if (flag & DD_FLAG_BW(d2))
4945 if (dd->nc[dd->dim[d2]] > 2)
4957 nrcg = flag & DD_FLAG_NRCG;
4960 if (home_pos_cg+1 > dd->cg_nalloc)
4962 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4963 srenew(dd->index_gl, dd->cg_nalloc);
4964 srenew(dd->cgindex, dd->cg_nalloc+1);
4966 /* Set the global charge group index and size */
4967 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4968 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4969 /* Copy the state from the buffer */
4970 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4971 if (fr->cutoff_scheme == ecutsGROUP)
4974 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4978 /* Set the cginfo */
4979 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4980 dd->index_gl[home_pos_cg]);
4983 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4986 if (home_pos_at+nrcg > state->nalloc)
4988 dd_realloc_state(state, f, home_pos_at+nrcg);
4990 for (i = 0; i < nrcg; i++)
4992 copy_rvec(comm->vbuf.v[buf_pos++],
4993 state->x[home_pos_at+i]);
4997 for (i = 0; i < nrcg; i++)
4999 copy_rvec(comm->vbuf.v[buf_pos++],
5000 state->v[home_pos_at+i]);
5005 for (i = 0; i < nrcg; i++)
5007 copy_rvec(comm->vbuf.v[buf_pos++],
5008 state->sd_X[home_pos_at+i]);
5013 for (i = 0; i < nrcg; i++)
5015 copy_rvec(comm->vbuf.v[buf_pos++],
5016 state->cg_p[home_pos_at+i]);
5020 home_pos_at += nrcg;
5024 /* Reallocate the buffers if necessary */
5025 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5027 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5028 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5030 nvr = ncg[mc] + nat[mc]*nvec;
5031 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5033 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5034 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5036 /* Copy from the receive to the send buffers */
5037 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5038 comm->buf_int + cg*DD_CGIBS,
5039 DD_CGIBS*sizeof(int));
5040 memcpy(comm->cgcm_state[mc][nvr],
5041 comm->vbuf.v[buf_pos],
5042 (1+nrcg*nvec)*sizeof(rvec));
5043 buf_pos += 1 + nrcg*nvec;
5050 /* With sorting (!bCompact) the indices are now only partially up to date
5051 * and ncg_home and nat_home are not the real count, since there are
5052 * "holes" in the arrays for the charge groups that moved to neighbors.
5054 if (fr->cutoff_scheme == ecutsVERLET)
5056 moved = get_moved(comm, home_pos_cg);
5058 for (i = dd->ncg_home; i < home_pos_cg; i++)
5063 dd->ncg_home = home_pos_cg;
5064 dd->nat_home = home_pos_at;
5069 "Finished repartitioning: cgs moved out %d, new home %d\n",
5070 *ncg_moved, dd->ncg_home-*ncg_moved);
5075 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5077 dd->comm->cycl[ddCycl] += cycles;
5078 dd->comm->cycl_n[ddCycl]++;
5079 if (cycles > dd->comm->cycl_max[ddCycl])
5081 dd->comm->cycl_max[ddCycl] = cycles;
5085 static double force_flop_count(t_nrnb *nrnb)
5092 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5094 /* To get closer to the real timings, we half the count
5095 * for the normal loops and again half it for water loops.
5098 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5100 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5104 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5107 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5110 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5112 sum += nrnb->n[i]*cost_nrnb(i);
5115 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5117 sum += nrnb->n[i]*cost_nrnb(i);
5123 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5125 if (dd->comm->eFlop)
5127 dd->comm->flop -= force_flop_count(nrnb);
5130 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5132 if (dd->comm->eFlop)
5134 dd->comm->flop += force_flop_count(nrnb);
5139 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5143 for (i = 0; i < ddCyclNr; i++)
5145 dd->comm->cycl[i] = 0;
5146 dd->comm->cycl_n[i] = 0;
5147 dd->comm->cycl_max[i] = 0;
5150 dd->comm->flop_n = 0;
5153 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5155 gmx_domdec_comm_t *comm;
5156 gmx_domdec_load_t *load;
5157 gmx_domdec_root_t *root = NULL;
5158 int d, dim, cid, i, pos;
5159 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5164 fprintf(debug, "get_load_distribution start\n");
5167 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5171 bSepPME = (dd->pme_nodeid >= 0);
5173 for (d = dd->ndim-1; d >= 0; d--)
5176 /* Check if we participate in the communication in this dimension */
5177 if (d == dd->ndim-1 ||
5178 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5180 load = &comm->load[d];
5183 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5186 if (d == dd->ndim-1)
5188 sbuf[pos++] = dd_force_load(comm);
5189 sbuf[pos++] = sbuf[0];
5192 sbuf[pos++] = sbuf[0];
5193 sbuf[pos++] = cell_frac;
5196 sbuf[pos++] = comm->cell_f_max0[d];
5197 sbuf[pos++] = comm->cell_f_min1[d];
5202 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5203 sbuf[pos++] = comm->cycl[ddCyclPME];
5208 sbuf[pos++] = comm->load[d+1].sum;
5209 sbuf[pos++] = comm->load[d+1].max;
5212 sbuf[pos++] = comm->load[d+1].sum_m;
5213 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5214 sbuf[pos++] = comm->load[d+1].flags;
5217 sbuf[pos++] = comm->cell_f_max0[d];
5218 sbuf[pos++] = comm->cell_f_min1[d];
5223 sbuf[pos++] = comm->load[d+1].mdf;
5224 sbuf[pos++] = comm->load[d+1].pme;
5228 /* Communicate a row in DD direction d.
5229 * The communicators are setup such that the root always has rank 0.
5232 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5233 load->load, load->nload*sizeof(float), MPI_BYTE,
5234 0, comm->mpi_comm_load[d]);
5236 if (dd->ci[dim] == dd->master_ci[dim])
5238 /* We are the root, process this row */
5239 if (comm->bDynLoadBal)
5241 root = comm->root[d];
5251 for (i = 0; i < dd->nc[dim]; i++)
5253 load->sum += load->load[pos++];
5254 load->max = max(load->max, load->load[pos]);
5260 /* This direction could not be load balanced properly,
5261 * therefore we need to use the maximum iso the average load.
5263 load->sum_m = max(load->sum_m, load->load[pos]);
5267 load->sum_m += load->load[pos];
5270 load->cvol_min = min(load->cvol_min, load->load[pos]);
5274 load->flags = (int)(load->load[pos++] + 0.5);
5278 root->cell_f_max0[i] = load->load[pos++];
5279 root->cell_f_min1[i] = load->load[pos++];
5284 load->mdf = max(load->mdf, load->load[pos]);
5286 load->pme = max(load->pme, load->load[pos]);
5290 if (comm->bDynLoadBal && root->bLimited)
5292 load->sum_m *= dd->nc[dim];
5293 load->flags |= (1<<d);
5301 comm->nload += dd_load_count(comm);
5302 comm->load_step += comm->cycl[ddCyclStep];
5303 comm->load_sum += comm->load[0].sum;
5304 comm->load_max += comm->load[0].max;
5305 if (comm->bDynLoadBal)
5307 for (d = 0; d < dd->ndim; d++)
5309 if (comm->load[0].flags & (1<<d))
5311 comm->load_lim[d]++;
5317 comm->load_mdf += comm->load[0].mdf;
5318 comm->load_pme += comm->load[0].pme;
5322 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5326 fprintf(debug, "get_load_distribution finished\n");
5330 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5332 /* Return the relative performance loss on the total run time
5333 * due to the force calculation load imbalance.
5335 if (dd->comm->nload > 0)
5338 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5339 (dd->comm->load_step*dd->nnodes);
5347 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5350 int npp, npme, nnodes, d, limp;
5351 float imbal, pme_f_ratio, lossf, lossp = 0;
5353 gmx_domdec_comm_t *comm;
5356 if (DDMASTER(dd) && comm->nload > 0)
5359 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5360 nnodes = npp + npme;
5361 imbal = comm->load_max*npp/comm->load_sum - 1;
5362 lossf = dd_force_imb_perf_loss(dd);
5363 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5364 fprintf(fplog, "%s", buf);
5365 fprintf(stderr, "\n");
5366 fprintf(stderr, "%s", buf);
5367 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5368 fprintf(fplog, "%s", buf);
5369 fprintf(stderr, "%s", buf);
5371 if (comm->bDynLoadBal)
5373 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5374 for (d = 0; d < dd->ndim; d++)
5376 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5377 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5383 sprintf(buf+strlen(buf), "\n");
5384 fprintf(fplog, "%s", buf);
5385 fprintf(stderr, "%s", buf);
5389 pme_f_ratio = comm->load_pme/comm->load_mdf;
5390 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5393 lossp *= (float)npme/(float)nnodes;
5397 lossp *= (float)npp/(float)nnodes;
5399 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5400 fprintf(fplog, "%s", buf);
5401 fprintf(stderr, "%s", buf);
5402 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5403 fprintf(fplog, "%s", buf);
5404 fprintf(stderr, "%s", buf);
5406 fprintf(fplog, "\n");
5407 fprintf(stderr, "\n");
5409 if (lossf >= DD_PERF_LOSS)
5412 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5413 " in the domain decomposition.\n", lossf*100);
5414 if (!comm->bDynLoadBal)
5416 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5420 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5422 fprintf(fplog, "%s\n", buf);
5423 fprintf(stderr, "%s\n", buf);
5425 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5428 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5429 " had %s work to do than the PP nodes.\n"
5430 " You might want to %s the number of PME nodes\n"
5431 " or %s the cut-off and the grid spacing.\n",
5433 (lossp < 0) ? "less" : "more",
5434 (lossp < 0) ? "decrease" : "increase",
5435 (lossp < 0) ? "decrease" : "increase");
5436 fprintf(fplog, "%s\n", buf);
5437 fprintf(stderr, "%s\n", buf);
5442 static float dd_vol_min(gmx_domdec_t *dd)
5444 return dd->comm->load[0].cvol_min*dd->nnodes;
5447 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5449 return dd->comm->load[0].flags;
5452 static float dd_f_imbal(gmx_domdec_t *dd)
5454 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5457 float dd_pme_f_ratio(gmx_domdec_t *dd)
5459 if (dd->comm->cycl_n[ddCyclPME] > 0)
5461 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5469 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5474 flags = dd_load_flags(dd);
5478 "DD load balancing is limited by minimum cell size in dimension");
5479 for (d = 0; d < dd->ndim; d++)
5483 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5486 fprintf(fplog, "\n");
5488 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5489 if (dd->comm->bDynLoadBal)
5491 fprintf(fplog, " vol min/aver %5.3f%c",
5492 dd_vol_min(dd), flags ? '!' : ' ');
5494 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5495 if (dd->comm->cycl_n[ddCyclPME])
5497 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5499 fprintf(fplog, "\n\n");
5502 static void dd_print_load_verbose(gmx_domdec_t *dd)
5504 if (dd->comm->bDynLoadBal)
5506 fprintf(stderr, "vol %4.2f%c ",
5507 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5509 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5510 if (dd->comm->cycl_n[ddCyclPME])
5512 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5517 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5522 gmx_domdec_root_t *root;
5523 gmx_bool bPartOfGroup = FALSE;
5525 dim = dd->dim[dim_ind];
5526 copy_ivec(loc, loc_c);
5527 for (i = 0; i < dd->nc[dim]; i++)
5530 rank = dd_index(dd->nc, loc_c);
5531 if (rank == dd->rank)
5533 /* This process is part of the group */
5534 bPartOfGroup = TRUE;
5537 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5541 dd->comm->mpi_comm_load[dim_ind] = c_row;
5542 if (dd->comm->eDLB != edlbNO)
5544 if (dd->ci[dim] == dd->master_ci[dim])
5546 /* This is the root process of this row */
5547 snew(dd->comm->root[dim_ind], 1);
5548 root = dd->comm->root[dim_ind];
5549 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5550 snew(root->old_cell_f, dd->nc[dim]+1);
5551 snew(root->bCellMin, dd->nc[dim]);
5554 snew(root->cell_f_max0, dd->nc[dim]);
5555 snew(root->cell_f_min1, dd->nc[dim]);
5556 snew(root->bound_min, dd->nc[dim]);
5557 snew(root->bound_max, dd->nc[dim]);
5559 snew(root->buf_ncd, dd->nc[dim]);
5563 /* This is not a root process, we only need to receive cell_f */
5564 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5567 if (dd->ci[dim] == dd->master_ci[dim])
5569 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5575 static void make_load_communicators(gmx_domdec_t *dd)
5578 int dim0, dim1, i, j;
5583 fprintf(debug, "Making load communicators\n");
5586 snew(dd->comm->load, dd->ndim);
5587 snew(dd->comm->mpi_comm_load, dd->ndim);
5590 make_load_communicator(dd, 0, loc);
5594 for (i = 0; i < dd->nc[dim0]; i++)
5597 make_load_communicator(dd, 1, loc);
5603 for (i = 0; i < dd->nc[dim0]; i++)
5607 for (j = 0; j < dd->nc[dim1]; j++)
5610 make_load_communicator(dd, 2, loc);
5617 fprintf(debug, "Finished making load communicators\n");
5622 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5625 int d, dim, i, j, m;
5628 ivec dd_zp[DD_MAXIZONE];
5629 gmx_domdec_zones_t *zones;
5630 gmx_domdec_ns_ranges_t *izone;
5632 for (d = 0; d < dd->ndim; d++)
5635 copy_ivec(dd->ci, tmp);
5636 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5637 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5638 copy_ivec(dd->ci, tmp);
5639 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5640 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5643 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5646 dd->neighbor[d][1]);
5652 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5654 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5655 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5662 for (i = 0; i < nzonep; i++)
5664 copy_ivec(dd_zp3[i], dd_zp[i]);
5670 for (i = 0; i < nzonep; i++)
5672 copy_ivec(dd_zp2[i], dd_zp[i]);
5678 for (i = 0; i < nzonep; i++)
5680 copy_ivec(dd_zp1[i], dd_zp[i]);
5684 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5689 zones = &dd->comm->zones;
5691 for (i = 0; i < nzone; i++)
5694 clear_ivec(zones->shift[i]);
5695 for (d = 0; d < dd->ndim; d++)
5697 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5702 for (i = 0; i < nzone; i++)
5704 for (d = 0; d < DIM; d++)
5706 s[d] = dd->ci[d] - zones->shift[i][d];
5711 else if (s[d] >= dd->nc[d])
5717 zones->nizone = nzonep;
5718 for (i = 0; i < zones->nizone; i++)
5720 if (dd_zp[i][0] != i)
5722 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5724 izone = &zones->izone[i];
5725 izone->j0 = dd_zp[i][1];
5726 izone->j1 = dd_zp[i][2];
5727 for (dim = 0; dim < DIM; dim++)
5729 if (dd->nc[dim] == 1)
5731 /* All shifts should be allowed */
5732 izone->shift0[dim] = -1;
5733 izone->shift1[dim] = 1;
5738 izone->shift0[d] = 0;
5739 izone->shift1[d] = 0;
5740 for(j=izone->j0; j<izone->j1; j++) {
5741 if (dd->shift[j][d] > dd->shift[i][d])
5742 izone->shift0[d] = -1;
5743 if (dd->shift[j][d] < dd->shift[i][d])
5744 izone->shift1[d] = 1;
5750 /* Assume the shift are not more than 1 cell */
5751 izone->shift0[dim] = 1;
5752 izone->shift1[dim] = -1;
5753 for (j = izone->j0; j < izone->j1; j++)
5755 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5756 if (shift_diff < izone->shift0[dim])
5758 izone->shift0[dim] = shift_diff;
5760 if (shift_diff > izone->shift1[dim])
5762 izone->shift1[dim] = shift_diff;
5769 if (dd->comm->eDLB != edlbNO)
5771 snew(dd->comm->root, dd->ndim);
5774 if (dd->comm->bRecordLoad)
5776 make_load_communicators(dd);
5780 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5783 gmx_domdec_comm_t *comm;
5794 if (comm->bCartesianPP)
5796 /* Set up cartesian communication for the particle-particle part */
5799 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5800 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5803 for (i = 0; i < DIM; i++)
5807 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5809 /* We overwrite the old communicator with the new cartesian one */
5810 cr->mpi_comm_mygroup = comm_cart;
5813 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5814 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5816 if (comm->bCartesianPP_PME)
5818 /* Since we want to use the original cartesian setup for sim,
5819 * and not the one after split, we need to make an index.
5821 snew(comm->ddindex2ddnodeid, dd->nnodes);
5822 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5823 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5824 /* Get the rank of the DD master,
5825 * above we made sure that the master node is a PP node.
5835 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5837 else if (comm->bCartesianPP)
5839 if (cr->npmenodes == 0)
5841 /* The PP communicator is also
5842 * the communicator for this simulation
5844 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5846 cr->nodeid = dd->rank;
5848 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5850 /* We need to make an index to go from the coordinates
5851 * to the nodeid of this simulation.
5853 snew(comm->ddindex2simnodeid, dd->nnodes);
5854 snew(buf, dd->nnodes);
5855 if (cr->duty & DUTY_PP)
5857 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5859 /* Communicate the ddindex to simulation nodeid index */
5860 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5861 cr->mpi_comm_mysim);
5864 /* Determine the master coordinates and rank.
5865 * The DD master should be the same node as the master of this sim.
5867 for (i = 0; i < dd->nnodes; i++)
5869 if (comm->ddindex2simnodeid[i] == 0)
5871 ddindex2xyz(dd->nc, i, dd->master_ci);
5872 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5877 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5882 /* No Cartesian communicators */
5883 /* We use the rank in dd->comm->all as DD index */
5884 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5885 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5887 clear_ivec(dd->master_ci);
5894 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5895 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5900 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5901 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5905 static void receive_ddindex2simnodeid(t_commrec *cr)
5909 gmx_domdec_comm_t *comm;
5916 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5918 snew(comm->ddindex2simnodeid, dd->nnodes);
5919 snew(buf, dd->nnodes);
5920 if (cr->duty & DUTY_PP)
5922 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5925 /* Communicate the ddindex to simulation nodeid index */
5926 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5927 cr->mpi_comm_mysim);
5934 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5935 int ncg, int natoms)
5937 gmx_domdec_master_t *ma;
5942 snew(ma->ncg, dd->nnodes);
5943 snew(ma->index, dd->nnodes+1);
5945 snew(ma->nat, dd->nnodes);
5946 snew(ma->ibuf, dd->nnodes*2);
5947 snew(ma->cell_x, DIM);
5948 for (i = 0; i < DIM; i++)
5950 snew(ma->cell_x[i], dd->nc[i]+1);
5953 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5959 snew(ma->vbuf, natoms);
5965 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5969 gmx_domdec_comm_t *comm;
5980 if (comm->bCartesianPP)
5982 for (i = 1; i < DIM; i++)
5984 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5986 if (bDiv[YY] || bDiv[ZZ])
5988 comm->bCartesianPP_PME = TRUE;
5989 /* If we have 2D PME decomposition, which is always in x+y,
5990 * we stack the PME only nodes in z.
5991 * Otherwise we choose the direction that provides the thinnest slab
5992 * of PME only nodes as this will have the least effect
5993 * on the PP communication.
5994 * But for the PME communication the opposite might be better.
5996 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
5998 dd->nc[YY] > dd->nc[ZZ]))
6000 comm->cartpmedim = ZZ;
6004 comm->cartpmedim = YY;
6006 comm->ntot[comm->cartpmedim]
6007 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6011 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]);
6013 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6018 if (comm->bCartesianPP_PME)
6022 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]);
6025 for (i = 0; i < DIM; i++)
6029 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6032 MPI_Comm_rank(comm_cart, &rank);
6033 if (MASTERNODE(cr) && rank != 0)
6035 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6038 /* With this assigment we loose the link to the original communicator
6039 * which will usually be MPI_COMM_WORLD, unless have multisim.
6041 cr->mpi_comm_mysim = comm_cart;
6042 cr->sim_nodeid = rank;
6044 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6048 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6049 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6052 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6056 if (cr->npmenodes == 0 ||
6057 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6059 cr->duty = DUTY_PME;
6062 /* Split the sim communicator into PP and PME only nodes */
6063 MPI_Comm_split(cr->mpi_comm_mysim,
6065 dd_index(comm->ntot, dd->ci),
6066 &cr->mpi_comm_mygroup);
6070 switch (dd_node_order)
6075 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6078 case ddnoINTERLEAVE:
6079 /* Interleave the PP-only and PME-only nodes,
6080 * as on clusters with dual-core machines this will double
6081 * the communication bandwidth of the PME processes
6082 * and thus speed up the PP <-> PME and inter PME communication.
6086 fprintf(fplog, "Interleaving PP and PME nodes\n");
6088 comm->pmenodes = dd_pmenodes(cr);
6093 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6096 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6098 cr->duty = DUTY_PME;
6105 /* Split the sim communicator into PP and PME only nodes */
6106 MPI_Comm_split(cr->mpi_comm_mysim,
6109 &cr->mpi_comm_mygroup);
6110 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6116 fprintf(fplog, "This is a %s only node\n\n",
6117 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6121 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6124 gmx_domdec_comm_t *comm;
6130 copy_ivec(dd->nc, comm->ntot);
6132 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6133 comm->bCartesianPP_PME = FALSE;
6135 /* Reorder the nodes by default. This might change the MPI ranks.
6136 * Real reordering is only supported on very few architectures,
6137 * Blue Gene is one of them.
6139 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6141 if (cr->npmenodes > 0)
6143 /* Split the communicator into a PP and PME part */
6144 split_communicator(fplog, cr, dd_node_order, CartReorder);
6145 if (comm->bCartesianPP_PME)
6147 /* We (possibly) reordered the nodes in split_communicator,
6148 * so it is no longer required in make_pp_communicator.
6150 CartReorder = FALSE;
6155 /* All nodes do PP and PME */
6157 /* We do not require separate communicators */
6158 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6162 if (cr->duty & DUTY_PP)
6164 /* Copy or make a new PP communicator */
6165 make_pp_communicator(fplog, cr, CartReorder);
6169 receive_ddindex2simnodeid(cr);
6172 if (!(cr->duty & DUTY_PME))
6174 /* Set up the commnuication to our PME node */
6175 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6176 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6179 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6180 dd->pme_nodeid, dd->pme_receive_vir_ener);
6185 dd->pme_nodeid = -1;
6190 dd->ma = init_gmx_domdec_master_t(dd,
6192 comm->cgs_gl.index[comm->cgs_gl.nr]);
6196 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6198 real *slb_frac, tot;
6203 if (nc > 1 && size_string != NULL)
6207 fprintf(fplog, "Using static load balancing for the %s direction\n",
6212 for (i = 0; i < nc; i++)
6215 sscanf(size_string, "%lf%n", &dbl, &n);
6218 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6227 fprintf(fplog, "Relative cell sizes:");
6229 for (i = 0; i < nc; i++)
6234 fprintf(fplog, " %5.3f", slb_frac[i]);
6239 fprintf(fplog, "\n");
6246 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6249 gmx_mtop_ilistloop_t iloop;
6253 iloop = gmx_mtop_ilistloop_init(mtop);
6254 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6256 for (ftype = 0; ftype < F_NRE; ftype++)
6258 if ((interaction_function[ftype].flags & IF_BOND) &&
6261 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6269 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6275 val = getenv(env_var);
6278 if (sscanf(val, "%d", &nst) <= 0)
6284 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6292 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6296 fprintf(stderr, "\n%s\n", warn_string);
6300 fprintf(fplog, "\n%s\n", warn_string);
6304 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6305 t_inputrec *ir, FILE *fplog)
6307 if (ir->ePBC == epbcSCREW &&
6308 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6310 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6313 if (ir->ns_type == ensSIMPLE)
6315 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6318 if (ir->nstlist == 0)
6320 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6323 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6325 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6329 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6334 r = ddbox->box_size[XX];
6335 for (di = 0; di < dd->ndim; di++)
6338 /* Check using the initial average cell size */
6339 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6345 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6346 const char *dlb_opt, gmx_bool bRecordLoad,
6347 unsigned long Flags, t_inputrec *ir)
6355 case 'a': eDLB = edlbAUTO; break;
6356 case 'n': eDLB = edlbNO; break;
6357 case 'y': eDLB = edlbYES; break;
6358 default: gmx_incons("Unknown dlb_opt");
6361 if (Flags & MD_RERUN)
6366 if (!EI_DYNAMICS(ir->eI))
6368 if (eDLB == edlbYES)
6370 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6371 dd_warning(cr, fplog, buf);
6379 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6384 if (Flags & MD_REPRODUCIBLE)
6391 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6395 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6398 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6406 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6411 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6413 /* Decomposition order z,y,x */
6416 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6418 for (dim = DIM-1; dim >= 0; dim--)
6420 if (dd->nc[dim] > 1)
6422 dd->dim[dd->ndim++] = dim;
6428 /* Decomposition order x,y,z */
6429 for (dim = 0; dim < DIM; dim++)
6431 if (dd->nc[dim] > 1)
6433 dd->dim[dd->ndim++] = dim;
6439 static gmx_domdec_comm_t *init_dd_comm()
6441 gmx_domdec_comm_t *comm;
6445 snew(comm->cggl_flag, DIM*2);
6446 snew(comm->cgcm_state, DIM*2);
6447 for (i = 0; i < DIM*2; i++)
6449 comm->cggl_flag_nalloc[i] = 0;
6450 comm->cgcm_state_nalloc[i] = 0;
6453 comm->nalloc_int = 0;
6454 comm->buf_int = NULL;
6456 vec_rvec_init(&comm->vbuf);
6458 comm->n_load_have = 0;
6459 comm->n_load_collect = 0;
6461 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6463 comm->sum_nat[i] = 0;
6467 comm->load_step = 0;
6470 clear_ivec(comm->load_lim);
6477 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6478 unsigned long Flags,
6480 real comm_distance_min, real rconstr,
6481 const char *dlb_opt, real dlb_scale,
6482 const char *sizex, const char *sizey, const char *sizez,
6483 gmx_mtop_t *mtop, t_inputrec *ir,
6484 matrix box, rvec *x,
6486 int *npme_x, int *npme_y)
6489 gmx_domdec_comm_t *comm;
6492 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6499 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6504 dd->comm = init_dd_comm();
6506 snew(comm->cggl_flag, DIM*2);
6507 snew(comm->cgcm_state, DIM*2);
6509 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6510 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6512 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6513 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6514 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6515 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6516 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6517 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6518 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6519 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6521 dd->pme_recv_f_alloc = 0;
6522 dd->pme_recv_f_buf = NULL;
6524 if (dd->bSendRecv2 && fplog)
6526 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");
6532 fprintf(fplog, "Will load balance based on FLOP count\n");
6534 if (comm->eFlop > 1)
6536 srand(1+cr->nodeid);
6538 comm->bRecordLoad = TRUE;
6542 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6546 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6548 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6551 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6553 dd->bGridJump = comm->bDynLoadBal;
6554 comm->bPMELoadBalDLBLimits = FALSE;
6556 if (comm->nstSortCG)
6560 if (comm->nstSortCG == 1)
6562 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6566 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6570 snew(comm->sort, 1);
6576 fprintf(fplog, "Will not sort the charge groups\n");
6580 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6582 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6583 if (comm->bInterCGBondeds)
6585 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6589 comm->bInterCGMultiBody = FALSE;
6592 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6593 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6595 if (ir->rlistlong == 0)
6597 /* Set the cut-off to some very large value,
6598 * so we don't need if statements everywhere in the code.
6599 * We use sqrt, since the cut-off is squared in some places.
6601 comm->cutoff = GMX_CUTOFF_INF;
6605 comm->cutoff = ir->rlistlong;
6607 comm->cutoff_mbody = 0;
6609 comm->cellsize_limit = 0;
6610 comm->bBondComm = FALSE;
6612 if (comm->bInterCGBondeds)
6614 if (comm_distance_min > 0)
6616 comm->cutoff_mbody = comm_distance_min;
6617 if (Flags & MD_DDBONDCOMM)
6619 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6623 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6625 r_bonded_limit = comm->cutoff_mbody;
6627 else if (ir->bPeriodicMols)
6629 /* Can not easily determine the required cut-off */
6630 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");
6631 comm->cutoff_mbody = comm->cutoff/2;
6632 r_bonded_limit = comm->cutoff_mbody;
6638 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6639 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6641 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6642 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6644 /* We use an initial margin of 10% for the minimum cell size,
6645 * except when we are just below the non-bonded cut-off.
6647 if (Flags & MD_DDBONDCOMM)
6649 if (max(r_2b, r_mb) > comm->cutoff)
6651 r_bonded = max(r_2b, r_mb);
6652 r_bonded_limit = 1.1*r_bonded;
6653 comm->bBondComm = TRUE;
6658 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6660 /* We determine cutoff_mbody later */
6664 /* No special bonded communication,
6665 * simply increase the DD cut-off.
6667 r_bonded_limit = 1.1*max(r_2b, r_mb);
6668 comm->cutoff_mbody = r_bonded_limit;
6669 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6672 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6676 "Minimum cell size due to bonded interactions: %.3f nm\n",
6677 comm->cellsize_limit);
6681 if (dd->bInterCGcons && rconstr <= 0)
6683 /* There is a cell size limit due to the constraints (P-LINCS) */
6684 rconstr = constr_r_max(fplog, mtop, ir);
6688 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6690 if (rconstr > comm->cellsize_limit)
6692 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6696 else if (rconstr > 0 && fplog)
6698 /* Here we do not check for dd->bInterCGcons,
6699 * because one can also set a cell size limit for virtual sites only
6700 * and at this point we don't know yet if there are intercg v-sites.
6703 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6706 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6708 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6712 copy_ivec(nc, dd->nc);
6713 set_dd_dim(fplog, dd);
6714 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6716 if (cr->npmenodes == -1)
6720 acs = average_cellsize_min(dd, ddbox);
6721 if (acs < comm->cellsize_limit)
6725 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6727 gmx_fatal_collective(FARGS, cr, NULL,
6728 "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",
6729 acs, comm->cellsize_limit);
6734 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6736 /* We need to choose the optimal DD grid and possibly PME nodes */
6737 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6738 comm->eDLB != edlbNO, dlb_scale,
6739 comm->cellsize_limit, comm->cutoff,
6740 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6742 if (dd->nc[XX] == 0)
6744 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6745 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6746 !bC ? "-rdd" : "-rcon",
6747 comm->eDLB != edlbNO ? " or -dds" : "",
6748 bC ? " or your LINCS settings" : "");
6750 gmx_fatal_collective(FARGS, cr, NULL,
6751 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6753 "Look in the log file for details on the domain decomposition",
6754 cr->nnodes-cr->npmenodes, limit, buf);
6756 set_dd_dim(fplog, dd);
6762 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6763 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6766 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6767 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6769 gmx_fatal_collective(FARGS, cr, NULL,
6770 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6771 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6773 if (cr->npmenodes > dd->nnodes)
6775 gmx_fatal_collective(FARGS, cr, NULL,
6776 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6778 if (cr->npmenodes > 0)
6780 comm->npmenodes = cr->npmenodes;
6784 comm->npmenodes = dd->nnodes;
6787 if (EEL_PME(ir->coulombtype))
6789 /* The following choices should match those
6790 * in comm_cost_est in domdec_setup.c.
6791 * Note that here the checks have to take into account
6792 * that the decomposition might occur in a different order than xyz
6793 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6794 * in which case they will not match those in comm_cost_est,
6795 * but since that is mainly for testing purposes that's fine.
6797 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6798 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6799 getenv("GMX_PMEONEDD") == NULL)
6801 comm->npmedecompdim = 2;
6802 comm->npmenodes_x = dd->nc[XX];
6803 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6807 /* In case nc is 1 in both x and y we could still choose to
6808 * decompose pme in y instead of x, but we use x for simplicity.
6810 comm->npmedecompdim = 1;
6811 if (dd->dim[0] == YY)
6813 comm->npmenodes_x = 1;
6814 comm->npmenodes_y = comm->npmenodes;
6818 comm->npmenodes_x = comm->npmenodes;
6819 comm->npmenodes_y = 1;
6824 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6825 comm->npmenodes_x, comm->npmenodes_y, 1);
6830 comm->npmedecompdim = 0;
6831 comm->npmenodes_x = 0;
6832 comm->npmenodes_y = 0;
6835 /* Technically we don't need both of these,
6836 * but it simplifies code not having to recalculate it.
6838 *npme_x = comm->npmenodes_x;
6839 *npme_y = comm->npmenodes_y;
6841 snew(comm->slb_frac, DIM);
6842 if (comm->eDLB == edlbNO)
6844 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6845 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6846 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6849 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6851 if (comm->bBondComm || comm->eDLB != edlbNO)
6853 /* Set the bonded communication distance to halfway
6854 * the minimum and the maximum,
6855 * since the extra communication cost is nearly zero.
6857 acs = average_cellsize_min(dd, ddbox);
6858 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6859 if (comm->eDLB != edlbNO)
6861 /* Check if this does not limit the scaling */
6862 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6864 if (!comm->bBondComm)
6866 /* Without bBondComm do not go beyond the n.b. cut-off */
6867 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6868 if (comm->cellsize_limit >= comm->cutoff)
6870 /* We don't loose a lot of efficieny
6871 * when increasing it to the n.b. cut-off.
6872 * It can even be slightly faster, because we need
6873 * less checks for the communication setup.
6875 comm->cutoff_mbody = comm->cutoff;
6878 /* Check if we did not end up below our original limit */
6879 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6881 if (comm->cutoff_mbody > comm->cellsize_limit)
6883 comm->cellsize_limit = comm->cutoff_mbody;
6886 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6891 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6892 "cellsize limit %f\n",
6893 comm->bBondComm, comm->cellsize_limit);
6898 check_dd_restrictions(cr, dd, ir, fplog);
6901 comm->partition_step = INT_MIN;
6904 clear_dd_cycle_counts(dd);
6909 static void set_dlb_limits(gmx_domdec_t *dd)
6914 for (d = 0; d < dd->ndim; d++)
6916 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6917 dd->comm->cellsize_min[dd->dim[d]] =
6918 dd->comm->cellsize_min_dlb[dd->dim[d]];
6923 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6926 gmx_domdec_comm_t *comm;
6936 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);
6939 cellsize_min = comm->cellsize_min[dd->dim[0]];
6940 for (d = 1; d < dd->ndim; d++)
6942 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6945 if (cellsize_min < comm->cellsize_limit*1.05)
6947 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");
6949 /* Change DLB from "auto" to "no". */
6950 comm->eDLB = edlbNO;
6955 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6956 comm->bDynLoadBal = TRUE;
6957 dd->bGridJump = TRUE;
6961 /* We can set the required cell size info here,
6962 * so we do not need to communicate this.
6963 * The grid is completely uniform.
6965 for (d = 0; d < dd->ndim; d++)
6969 comm->load[d].sum_m = comm->load[d].sum;
6971 nc = dd->nc[dd->dim[d]];
6972 for (i = 0; i < nc; i++)
6974 comm->root[d]->cell_f[i] = i/(real)nc;
6977 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6978 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6981 comm->root[d]->cell_f[nc] = 1.0;
6986 static char *init_bLocalCG(gmx_mtop_t *mtop)
6991 ncg = ncg_mtop(mtop);
6992 snew(bLocalCG, ncg);
6993 for (cg = 0; cg < ncg; cg++)
6995 bLocalCG[cg] = FALSE;
7001 void dd_init_bondeds(FILE *fplog,
7002 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7003 gmx_vsite_t *vsite, gmx_constr_t constr,
7004 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7006 gmx_domdec_comm_t *comm;
7010 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7014 if (comm->bBondComm)
7016 /* Communicate atoms beyond the cut-off for bonded interactions */
7019 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7021 comm->bLocalCG = init_bLocalCG(mtop);
7025 /* Only communicate atoms based on cut-off */
7026 comm->cglink = NULL;
7027 comm->bLocalCG = NULL;
7031 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7033 gmx_bool bDynLoadBal, real dlb_scale,
7036 gmx_domdec_comm_t *comm;
7051 fprintf(fplog, "The maximum number of communication pulses is:");
7052 for (d = 0; d < dd->ndim; d++)
7054 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7056 fprintf(fplog, "\n");
7057 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7058 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7059 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7060 for (d = 0; d < DIM; d++)
7064 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7071 comm->cellsize_min_dlb[d]/
7072 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7074 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7077 fprintf(fplog, "\n");
7081 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7082 fprintf(fplog, "The initial number of communication pulses is:");
7083 for (d = 0; d < dd->ndim; d++)
7085 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7087 fprintf(fplog, "\n");
7088 fprintf(fplog, "The initial domain decomposition cell size is:");
7089 for (d = 0; d < DIM; d++)
7093 fprintf(fplog, " %c %.2f nm",
7094 dim2char(d), dd->comm->cellsize_min[d]);
7097 fprintf(fplog, "\n\n");
7100 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7102 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7103 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7104 "non-bonded interactions", "", comm->cutoff);
7108 limit = dd->comm->cellsize_limit;
7112 if (dynamic_dd_box(ddbox, ir))
7114 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7116 limit = dd->comm->cellsize_min[XX];
7117 for (d = 1; d < DIM; d++)
7119 limit = min(limit, dd->comm->cellsize_min[d]);
7123 if (comm->bInterCGBondeds)
7125 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7126 "two-body bonded interactions", "(-rdd)",
7127 max(comm->cutoff, comm->cutoff_mbody));
7128 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7129 "multi-body bonded interactions", "(-rdd)",
7130 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7134 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7135 "virtual site constructions", "(-rcon)", limit);
7137 if (dd->constraint_comm)
7139 sprintf(buf, "atoms separated by up to %d constraints",
7141 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7142 buf, "(-rcon)", limit);
7144 fprintf(fplog, "\n");
7150 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7152 const t_inputrec *ir,
7153 const gmx_ddbox_t *ddbox)
7155 gmx_domdec_comm_t *comm;
7156 int d, dim, npulse, npulse_d_max, npulse_d;
7161 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7163 /* Determine the maximum number of comm. pulses in one dimension */
7165 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7167 /* Determine the maximum required number of grid pulses */
7168 if (comm->cellsize_limit >= comm->cutoff)
7170 /* Only a single pulse is required */
7173 else if (!bNoCutOff && comm->cellsize_limit > 0)
7175 /* We round down slightly here to avoid overhead due to the latency
7176 * of extra communication calls when the cut-off
7177 * would be only slightly longer than the cell size.
7178 * Later cellsize_limit is redetermined,
7179 * so we can not miss interactions due to this rounding.
7181 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7185 /* There is no cell size limit */
7186 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7189 if (!bNoCutOff && npulse > 1)
7191 /* See if we can do with less pulses, based on dlb_scale */
7193 for (d = 0; d < dd->ndim; d++)
7196 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7197 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7198 npulse_d_max = max(npulse_d_max, npulse_d);
7200 npulse = min(npulse, npulse_d_max);
7203 /* This env var can override npulse */
7204 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7211 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7212 for (d = 0; d < dd->ndim; d++)
7214 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7215 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7216 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7217 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7218 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7220 comm->bVacDLBNoLimit = FALSE;
7224 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7225 if (!comm->bVacDLBNoLimit)
7227 comm->cellsize_limit = max(comm->cellsize_limit,
7228 comm->cutoff/comm->maxpulse);
7230 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7231 /* Set the minimum cell size for each DD dimension */
7232 for (d = 0; d < dd->ndim; d++)
7234 if (comm->bVacDLBNoLimit ||
7235 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7237 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7241 comm->cellsize_min_dlb[dd->dim[d]] =
7242 comm->cutoff/comm->cd[d].np_dlb;
7245 if (comm->cutoff_mbody <= 0)
7247 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7249 if (comm->bDynLoadBal)
7255 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7257 /* If each molecule is a single charge group
7258 * or we use domain decomposition for each periodic dimension,
7259 * we do not need to take pbc into account for the bonded interactions.
7261 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7264 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7267 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7268 t_inputrec *ir, t_forcerec *fr,
7271 gmx_domdec_comm_t *comm;
7277 /* Initialize the thread data.
7278 * This can not be done in init_domain_decomposition,
7279 * as the numbers of threads is determined later.
7281 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7284 snew(comm->dth, comm->nth);
7287 if (EEL_PME(ir->coulombtype))
7289 init_ddpme(dd, &comm->ddpme[0], 0);
7290 if (comm->npmedecompdim >= 2)
7292 init_ddpme(dd, &comm->ddpme[1], 1);
7297 comm->npmenodes = 0;
7298 if (dd->pme_nodeid >= 0)
7300 gmx_fatal_collective(FARGS, NULL, dd,
7301 "Can not have separate PME nodes without PME electrostatics");
7307 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7309 if (comm->eDLB != edlbNO)
7311 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7314 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7315 if (comm->eDLB == edlbAUTO)
7319 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7321 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7324 if (ir->ePBC == epbcNONE)
7326 vol_frac = 1 - 1/(double)dd->nnodes;
7331 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7335 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7337 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7339 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7342 static gmx_bool test_dd_cutoff(t_commrec *cr,
7343 t_state *state, t_inputrec *ir,
7354 set_ddbox(dd, FALSE, cr, ir, state->box,
7355 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7359 for (d = 0; d < dd->ndim; d++)
7363 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7364 if (dynamic_dd_box(&ddbox, ir))
7366 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7369 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7371 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7372 dd->comm->cd[d].np_dlb > 0)
7374 if (np > dd->comm->cd[d].np_dlb)
7379 /* If a current local cell size is smaller than the requested
7380 * cut-off, we could still fix it, but this gets very complicated.
7381 * Without fixing here, we might actually need more checks.
7383 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7390 if (dd->comm->eDLB != edlbNO)
7392 /* If DLB is not active yet, we don't need to check the grid jumps.
7393 * Actually we shouldn't, because then the grid jump data is not set.
7395 if (dd->comm->bDynLoadBal &&
7396 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7401 gmx_sumi(1, &LocallyLimited, cr);
7403 if (LocallyLimited > 0)
7412 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7415 gmx_bool bCutoffAllowed;
7417 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7421 cr->dd->comm->cutoff = cutoff_req;
7424 return bCutoffAllowed;
7427 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7429 gmx_domdec_comm_t *comm;
7431 comm = cr->dd->comm;
7433 /* Turn on the DLB limiting (might have been on already) */
7434 comm->bPMELoadBalDLBLimits = TRUE;
7436 /* Change the cut-off limit */
7437 comm->PMELoadBal_max_cutoff = comm->cutoff;
7440 static void merge_cg_buffers(int ncell,
7441 gmx_domdec_comm_dim_t *cd, int pulse,
7443 int *index_gl, int *recv_i,
7444 rvec *cg_cm, rvec *recv_vr,
7446 cginfo_mb_t *cginfo_mb, int *cginfo)
7448 gmx_domdec_ind_t *ind, *ind_p;
7449 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7450 int shift, shift_at;
7452 ind = &cd->ind[pulse];
7454 /* First correct the already stored data */
7455 shift = ind->nrecv[ncell];
7456 for (cell = ncell-1; cell >= 0; cell--)
7458 shift -= ind->nrecv[cell];
7461 /* Move the cg's present from previous grid pulses */
7462 cg0 = ncg_cell[ncell+cell];
7463 cg1 = ncg_cell[ncell+cell+1];
7464 cgindex[cg1+shift] = cgindex[cg1];
7465 for (cg = cg1-1; cg >= cg0; cg--)
7467 index_gl[cg+shift] = index_gl[cg];
7468 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7469 cgindex[cg+shift] = cgindex[cg];
7470 cginfo[cg+shift] = cginfo[cg];
7472 /* Correct the already stored send indices for the shift */
7473 for (p = 1; p <= pulse; p++)
7475 ind_p = &cd->ind[p];
7477 for (c = 0; c < cell; c++)
7479 cg0 += ind_p->nsend[c];
7481 cg1 = cg0 + ind_p->nsend[cell];
7482 for (cg = cg0; cg < cg1; cg++)
7484 ind_p->index[cg] += shift;
7490 /* Merge in the communicated buffers */
7494 for (cell = 0; cell < ncell; cell++)
7496 cg1 = ncg_cell[ncell+cell+1] + shift;
7499 /* Correct the old cg indices */
7500 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7502 cgindex[cg+1] += shift_at;
7505 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7507 /* Copy this charge group from the buffer */
7508 index_gl[cg1] = recv_i[cg0];
7509 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7510 /* Add it to the cgindex */
7511 cg_gl = index_gl[cg1];
7512 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7513 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7514 cgindex[cg1+1] = cgindex[cg1] + nat;
7519 shift += ind->nrecv[cell];
7520 ncg_cell[ncell+cell+1] = cg1;
7524 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7525 int nzone, int cg0, const int *cgindex)
7529 /* Store the atom block boundaries for easy copying of communication buffers
7532 for (zone = 0; zone < nzone; zone++)
7534 for (p = 0; p < cd->np; p++)
7536 cd->ind[p].cell2at0[zone] = cgindex[cg];
7537 cg += cd->ind[p].nrecv[zone];
7538 cd->ind[p].cell2at1[zone] = cgindex[cg];
7543 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7549 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7551 if (!bLocalCG[link->a[i]])
7560 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7562 real c[DIM][4]; /* the corners for the non-bonded communication */
7563 real cr0; /* corner for rounding */
7564 real cr1[4]; /* corners for rounding */
7565 real bc[DIM]; /* corners for bounded communication */
7566 real bcr1; /* corner for rounding for bonded communication */
7569 /* Determine the corners of the domain(s) we are communicating with */
7571 set_dd_corners(const gmx_domdec_t *dd,
7572 int dim0, int dim1, int dim2,
7576 const gmx_domdec_comm_t *comm;
7577 const gmx_domdec_zones_t *zones;
7582 zones = &comm->zones;
7584 /* Keep the compiler happy */
7588 /* The first dimension is equal for all cells */
7589 c->c[0][0] = comm->cell_x0[dim0];
7592 c->bc[0] = c->c[0][0];
7597 /* This cell row is only seen from the first row */
7598 c->c[1][0] = comm->cell_x0[dim1];
7599 /* All rows can see this row */
7600 c->c[1][1] = comm->cell_x0[dim1];
7603 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7606 /* For the multi-body distance we need the maximum */
7607 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7610 /* Set the upper-right corner for rounding */
7611 c->cr0 = comm->cell_x1[dim0];
7616 for (j = 0; j < 4; j++)
7618 c->c[2][j] = comm->cell_x0[dim2];
7622 /* Use the maximum of the i-cells that see a j-cell */
7623 for (i = 0; i < zones->nizone; i++)
7625 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7631 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7637 /* For the multi-body distance we need the maximum */
7638 c->bc[2] = comm->cell_x0[dim2];
7639 for (i = 0; i < 2; i++)
7641 for (j = 0; j < 2; j++)
7643 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7649 /* Set the upper-right corner for rounding */
7650 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7651 * Only cell (0,0,0) can see cell 7 (1,1,1)
7653 c->cr1[0] = comm->cell_x1[dim1];
7654 c->cr1[3] = comm->cell_x1[dim1];
7657 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7660 /* For the multi-body distance we need the maximum */
7661 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7668 /* Determine which cg's we need to send in this pulse from this zone */
7670 get_zone_pulse_cgs(gmx_domdec_t *dd,
7671 int zonei, int zone,
7673 const int *index_gl,
7675 int dim, int dim_ind,
7676 int dim0, int dim1, int dim2,
7677 real r_comm2, real r_bcomm2,
7681 real skew_fac2_d, real skew_fac_01,
7682 rvec *v_d, rvec *v_0, rvec *v_1,
7683 const dd_corners_t *c,
7685 gmx_bool bDistBonded,
7691 gmx_domdec_ind_t *ind,
7692 int **ibuf, int *ibuf_nalloc,
7698 gmx_domdec_comm_t *comm;
7700 gmx_bool bDistMB_pulse;
7702 real r2, rb2, r, tric_sh;
7705 int nsend_z, nsend, nat;
7709 bScrew = (dd->bScrewPBC && dim == XX);
7711 bDistMB_pulse = (bDistMB && bDistBonded);
7717 for (cg = cg0; cg < cg1; cg++)
7721 if (tric_dist[dim_ind] == 0)
7723 /* Rectangular direction, easy */
7724 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7731 r = cg_cm[cg][dim] - c->bc[dim_ind];
7737 /* Rounding gives at most a 16% reduction
7738 * in communicated atoms
7740 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7742 r = cg_cm[cg][dim0] - c->cr0;
7743 /* This is the first dimension, so always r >= 0 */
7750 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7752 r = cg_cm[cg][dim1] - c->cr1[zone];
7759 r = cg_cm[cg][dim1] - c->bcr1;
7769 /* Triclinic direction, more complicated */
7772 /* Rounding, conservative as the skew_fac multiplication
7773 * will slightly underestimate the distance.
7775 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7777 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7778 for (i = dim0+1; i < DIM; i++)
7780 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7782 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7785 rb[dim0] = rn[dim0];
7788 /* Take care that the cell planes along dim0 might not
7789 * be orthogonal to those along dim1 and dim2.
7791 for (i = 1; i <= dim_ind; i++)
7794 if (normal[dim0][dimd] > 0)
7796 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7799 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7804 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7806 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7808 for (i = dim1+1; i < DIM; i++)
7810 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7812 rn[dim1] += tric_sh;
7815 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7816 /* Take care of coupling of the distances
7817 * to the planes along dim0 and dim1 through dim2.
7819 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7820 /* Take care that the cell planes along dim1
7821 * might not be orthogonal to that along dim2.
7823 if (normal[dim1][dim2] > 0)
7825 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7831 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7834 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7835 /* Take care of coupling of the distances
7836 * to the planes along dim0 and dim1 through dim2.
7838 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7839 /* Take care that the cell planes along dim1
7840 * might not be orthogonal to that along dim2.
7842 if (normal[dim1][dim2] > 0)
7844 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7849 /* The distance along the communication direction */
7850 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7852 for (i = dim+1; i < DIM; i++)
7854 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7859 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7860 /* Take care of coupling of the distances
7861 * to the planes along dim0 and dim1 through dim2.
7863 if (dim_ind == 1 && zonei == 1)
7865 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7871 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7874 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7875 /* Take care of coupling of the distances
7876 * to the planes along dim0 and dim1 through dim2.
7878 if (dim_ind == 1 && zonei == 1)
7880 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7888 ((bDistMB && rb2 < r_bcomm2) ||
7889 (bDist2B && r2 < r_bcomm2)) &&
7891 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7892 missing_link(comm->cglink, index_gl[cg],
7895 /* Make an index to the local charge groups */
7896 if (nsend+1 > ind->nalloc)
7898 ind->nalloc = over_alloc_large(nsend+1);
7899 srenew(ind->index, ind->nalloc);
7901 if (nsend+1 > *ibuf_nalloc)
7903 *ibuf_nalloc = over_alloc_large(nsend+1);
7904 srenew(*ibuf, *ibuf_nalloc);
7906 ind->index[nsend] = cg;
7907 (*ibuf)[nsend] = index_gl[cg];
7909 vec_rvec_check_alloc(vbuf, nsend+1);
7911 if (dd->ci[dim] == 0)
7913 /* Correct cg_cm for pbc */
7914 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7917 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7918 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7923 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7926 nat += cgindex[cg+1] - cgindex[cg];
7932 *nsend_z_ptr = nsend_z;
7935 static void setup_dd_communication(gmx_domdec_t *dd,
7936 matrix box, gmx_ddbox_t *ddbox,
7937 t_forcerec *fr, t_state *state, rvec **f)
7939 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7940 int nzone, nzone_send, zone, zonei, cg0, cg1;
7941 int c, i, j, cg, cg_gl, nrcg;
7942 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7943 gmx_domdec_comm_t *comm;
7944 gmx_domdec_zones_t *zones;
7945 gmx_domdec_comm_dim_t *cd;
7946 gmx_domdec_ind_t *ind;
7947 cginfo_mb_t *cginfo_mb;
7948 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7949 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7950 dd_corners_t corners;
7952 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7953 real skew_fac2_d, skew_fac_01;
7960 fprintf(debug, "Setting up DD communication\n");
7965 switch (fr->cutoff_scheme)
7974 gmx_incons("unimplemented");
7978 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7980 dim = dd->dim[dim_ind];
7982 /* Check if we need to use triclinic distances */
7983 tric_dist[dim_ind] = 0;
7984 for (i = 0; i <= dim_ind; i++)
7986 if (ddbox->tric_dir[dd->dim[i]])
7988 tric_dist[dim_ind] = 1;
7993 bBondComm = comm->bBondComm;
7995 /* Do we need to determine extra distances for multi-body bondeds? */
7996 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
7998 /* Do we need to determine extra distances for only two-body bondeds? */
7999 bDist2B = (bBondComm && !bDistMB);
8001 r_comm2 = sqr(comm->cutoff);
8002 r_bcomm2 = sqr(comm->cutoff_mbody);
8006 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8009 zones = &comm->zones;
8012 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8013 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8015 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8017 /* Triclinic stuff */
8018 normal = ddbox->normal;
8022 v_0 = ddbox->v[dim0];
8023 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8025 /* Determine the coupling coefficient for the distances
8026 * to the cell planes along dim0 and dim1 through dim2.
8027 * This is required for correct rounding.
8030 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8033 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8039 v_1 = ddbox->v[dim1];
8042 zone_cg_range = zones->cg_range;
8043 index_gl = dd->index_gl;
8044 cgindex = dd->cgindex;
8045 cginfo_mb = fr->cginfo_mb;
8047 zone_cg_range[0] = 0;
8048 zone_cg_range[1] = dd->ncg_home;
8049 comm->zone_ncg1[0] = dd->ncg_home;
8050 pos_cg = dd->ncg_home;
8052 nat_tot = dd->nat_home;
8054 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8056 dim = dd->dim[dim_ind];
8057 cd = &comm->cd[dim_ind];
8059 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8061 /* No pbc in this dimension, the first node should not comm. */
8069 v_d = ddbox->v[dim];
8070 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8072 cd->bInPlace = TRUE;
8073 for (p = 0; p < cd->np; p++)
8075 /* Only atoms communicated in the first pulse are used
8076 * for multi-body bonded interactions or for bBondComm.
8078 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8083 for (zone = 0; zone < nzone_send; zone++)
8085 if (tric_dist[dim_ind] && dim_ind > 0)
8087 /* Determine slightly more optimized skew_fac's
8089 * This reduces the number of communicated atoms
8090 * by about 10% for 3D DD of rhombic dodecahedra.
8092 for (dimd = 0; dimd < dim; dimd++)
8094 sf2_round[dimd] = 1;
8095 if (ddbox->tric_dir[dimd])
8097 for (i = dd->dim[dimd]+1; i < DIM; i++)
8099 /* If we are shifted in dimension i
8100 * and the cell plane is tilted forward
8101 * in dimension i, skip this coupling.
8103 if (!(zones->shift[nzone+zone][i] &&
8104 ddbox->v[dimd][i][dimd] >= 0))
8107 sqr(ddbox->v[dimd][i][dimd]);
8110 sf2_round[dimd] = 1/sf2_round[dimd];
8115 zonei = zone_perm[dim_ind][zone];
8118 /* Here we permutate the zones to obtain a convenient order
8119 * for neighbor searching
8121 cg0 = zone_cg_range[zonei];
8122 cg1 = zone_cg_range[zonei+1];
8126 /* Look only at the cg's received in the previous grid pulse
8128 cg1 = zone_cg_range[nzone+zone+1];
8129 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8132 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8133 for (th = 0; th < comm->nth; th++)
8135 gmx_domdec_ind_t *ind_p;
8136 int **ibuf_p, *ibuf_nalloc_p;
8138 int *nsend_p, *nat_p;
8144 /* Thread 0 writes in the comm buffers */
8146 ibuf_p = &comm->buf_int;
8147 ibuf_nalloc_p = &comm->nalloc_int;
8148 vbuf_p = &comm->vbuf;
8151 nsend_zone_p = &ind->nsend[zone];
8155 /* Other threads write into temp buffers */
8156 ind_p = &comm->dth[th].ind;
8157 ibuf_p = &comm->dth[th].ibuf;
8158 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8159 vbuf_p = &comm->dth[th].vbuf;
8160 nsend_p = &comm->dth[th].nsend;
8161 nat_p = &comm->dth[th].nat;
8162 nsend_zone_p = &comm->dth[th].nsend_zone;
8164 comm->dth[th].nsend = 0;
8165 comm->dth[th].nat = 0;
8166 comm->dth[th].nsend_zone = 0;
8176 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8177 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8180 /* Get the cg's for this pulse in this zone */
8181 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8183 dim, dim_ind, dim0, dim1, dim2,
8186 normal, skew_fac2_d, skew_fac_01,
8187 v_d, v_0, v_1, &corners, sf2_round,
8188 bDistBonded, bBondComm,
8192 ibuf_p, ibuf_nalloc_p,
8198 /* Append data of threads>=1 to the communication buffers */
8199 for (th = 1; th < comm->nth; th++)
8201 dd_comm_setup_work_t *dth;
8204 dth = &comm->dth[th];
8206 ns1 = nsend + dth->nsend_zone;
8207 if (ns1 > ind->nalloc)
8209 ind->nalloc = over_alloc_dd(ns1);
8210 srenew(ind->index, ind->nalloc);
8212 if (ns1 > comm->nalloc_int)
8214 comm->nalloc_int = over_alloc_dd(ns1);
8215 srenew(comm->buf_int, comm->nalloc_int);
8217 if (ns1 > comm->vbuf.nalloc)
8219 comm->vbuf.nalloc = over_alloc_dd(ns1);
8220 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8223 for (i = 0; i < dth->nsend_zone; i++)
8225 ind->index[nsend] = dth->ind.index[i];
8226 comm->buf_int[nsend] = dth->ibuf[i];
8227 copy_rvec(dth->vbuf.v[i],
8228 comm->vbuf.v[nsend]);
8232 ind->nsend[zone] += dth->nsend_zone;
8235 /* Clear the counts in case we do not have pbc */
8236 for (zone = nzone_send; zone < nzone; zone++)
8238 ind->nsend[zone] = 0;
8240 ind->nsend[nzone] = nsend;
8241 ind->nsend[nzone+1] = nat;
8242 /* Communicate the number of cg's and atoms to receive */
8243 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8244 ind->nsend, nzone+2,
8245 ind->nrecv, nzone+2);
8247 /* The rvec buffer is also required for atom buffers of size nsend
8248 * in dd_move_x and dd_move_f.
8250 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8254 /* We can receive in place if only the last zone is not empty */
8255 for (zone = 0; zone < nzone-1; zone++)
8257 if (ind->nrecv[zone] > 0)
8259 cd->bInPlace = FALSE;
8264 /* The int buffer is only required here for the cg indices */
8265 if (ind->nrecv[nzone] > comm->nalloc_int2)
8267 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8268 srenew(comm->buf_int2, comm->nalloc_int2);
8270 /* The rvec buffer is also required for atom buffers
8271 * of size nrecv in dd_move_x and dd_move_f.
8273 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8274 vec_rvec_check_alloc(&comm->vbuf2, i);
8278 /* Make space for the global cg indices */
8279 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8280 || dd->cg_nalloc == 0)
8282 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8283 srenew(index_gl, dd->cg_nalloc);
8284 srenew(cgindex, dd->cg_nalloc+1);
8286 /* Communicate the global cg indices */
8289 recv_i = index_gl + pos_cg;
8293 recv_i = comm->buf_int2;
8295 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8296 comm->buf_int, nsend,
8297 recv_i, ind->nrecv[nzone]);
8299 /* Make space for cg_cm */
8300 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8301 if (fr->cutoff_scheme == ecutsGROUP)
8309 /* Communicate cg_cm */
8312 recv_vr = cg_cm + pos_cg;
8316 recv_vr = comm->vbuf2.v;
8318 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8319 comm->vbuf.v, nsend,
8320 recv_vr, ind->nrecv[nzone]);
8322 /* Make the charge group index */
8325 zone = (p == 0 ? 0 : nzone - 1);
8326 while (zone < nzone)
8328 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8330 cg_gl = index_gl[pos_cg];
8331 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8332 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8333 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8336 /* Update the charge group presence,
8337 * so we can use it in the next pass of the loop.
8339 comm->bLocalCG[cg_gl] = TRUE;
8345 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8348 zone_cg_range[nzone+zone] = pos_cg;
8353 /* This part of the code is never executed with bBondComm. */
8354 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8355 index_gl, recv_i, cg_cm, recv_vr,
8356 cgindex, fr->cginfo_mb, fr->cginfo);
8357 pos_cg += ind->nrecv[nzone];
8359 nat_tot += ind->nrecv[nzone+1];
8363 /* Store the atom block for easy copying of communication buffers */
8364 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8368 dd->index_gl = index_gl;
8369 dd->cgindex = cgindex;
8371 dd->ncg_tot = zone_cg_range[zones->n];
8372 dd->nat_tot = nat_tot;
8373 comm->nat[ddnatHOME] = dd->nat_home;
8374 for (i = ddnatZONE; i < ddnatNR; i++)
8376 comm->nat[i] = dd->nat_tot;
8381 /* We don't need to update cginfo, since that was alrady done above.
8382 * So we pass NULL for the forcerec.
8384 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8385 NULL, comm->bLocalCG);
8390 fprintf(debug, "Finished setting up DD communication, zones:");
8391 for (c = 0; c < zones->n; c++)
8393 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8395 fprintf(debug, "\n");
8399 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8403 for (c = 0; c < zones->nizone; c++)
8405 zones->izone[c].cg1 = zones->cg_range[c+1];
8406 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8407 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8411 static void set_zones_size(gmx_domdec_t *dd,
8412 matrix box, const gmx_ddbox_t *ddbox,
8413 int zone_start, int zone_end)
8415 gmx_domdec_comm_t *comm;
8416 gmx_domdec_zones_t *zones;
8418 int z, zi, zj0, zj1, d, dim;
8421 real size_j, add_tric;
8426 zones = &comm->zones;
8428 /* Do we need to determine extra distances for multi-body bondeds? */
8429 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8431 for (z = zone_start; z < zone_end; z++)
8433 /* Copy cell limits to zone limits.
8434 * Valid for non-DD dims and non-shifted dims.
8436 copy_rvec(comm->cell_x0, zones->size[z].x0);
8437 copy_rvec(comm->cell_x1, zones->size[z].x1);
8440 for (d = 0; d < dd->ndim; d++)
8444 for (z = 0; z < zones->n; z++)
8446 /* With a staggered grid we have different sizes
8447 * for non-shifted dimensions.
8449 if (dd->bGridJump && zones->shift[z][dim] == 0)
8453 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8454 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8458 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8459 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8465 rcmbs = comm->cutoff_mbody;
8466 if (ddbox->tric_dir[dim])
8468 rcs /= ddbox->skew_fac[dim];
8469 rcmbs /= ddbox->skew_fac[dim];
8472 /* Set the lower limit for the shifted zone dimensions */
8473 for (z = zone_start; z < zone_end; z++)
8475 if (zones->shift[z][dim] > 0)
8478 if (!dd->bGridJump || d == 0)
8480 zones->size[z].x0[dim] = comm->cell_x1[dim];
8481 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8485 /* Here we take the lower limit of the zone from
8486 * the lowest domain of the zone below.
8490 zones->size[z].x0[dim] =
8491 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8497 zones->size[z].x0[dim] =
8498 zones->size[zone_perm[2][z-4]].x0[dim];
8502 zones->size[z].x0[dim] =
8503 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8506 /* A temporary limit, is updated below */
8507 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8511 for (zi = 0; zi < zones->nizone; zi++)
8513 if (zones->shift[zi][dim] == 0)
8515 /* This takes the whole zone into account.
8516 * With multiple pulses this will lead
8517 * to a larger zone then strictly necessary.
8519 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8520 zones->size[zi].x1[dim]+rcmbs);
8528 /* Loop over the i-zones to set the upper limit of each
8531 for (zi = 0; zi < zones->nizone; zi++)
8533 if (zones->shift[zi][dim] == 0)
8535 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8537 if (zones->shift[z][dim] > 0)
8539 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8540 zones->size[zi].x1[dim]+rcs);
8547 for (z = zone_start; z < zone_end; z++)
8549 /* Initialization only required to keep the compiler happy */
8550 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8553 /* To determine the bounding box for a zone we need to find
8554 * the extreme corners of 4, 2 or 1 corners.
8556 nc = 1 << (ddbox->npbcdim - 1);
8558 for (c = 0; c < nc; c++)
8560 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8564 corner[YY] = zones->size[z].x0[YY];
8568 corner[YY] = zones->size[z].x1[YY];
8572 corner[ZZ] = zones->size[z].x0[ZZ];
8576 corner[ZZ] = zones->size[z].x1[ZZ];
8578 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8580 /* With 1D domain decomposition the cg's are not in
8581 * the triclinic box, but triclinic x-y and rectangular y-z.
8582 * Shift y back, so it will later end up at 0.
8584 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8586 /* Apply the triclinic couplings */
8587 for (i = YY; i < ddbox->npbcdim; i++)
8589 for (j = XX; j < i; j++)
8591 corner[j] += corner[i]*box[i][j]/box[i][i];
8596 copy_rvec(corner, corner_min);
8597 copy_rvec(corner, corner_max);
8601 for (i = 0; i < DIM; i++)
8603 corner_min[i] = min(corner_min[i], corner[i]);
8604 corner_max[i] = max(corner_max[i], corner[i]);
8608 /* Copy the extreme cornes without offset along x */
8609 for (i = 0; i < DIM; i++)
8611 zones->size[z].bb_x0[i] = corner_min[i];
8612 zones->size[z].bb_x1[i] = corner_max[i];
8614 /* Add the offset along x */
8615 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8616 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8619 if (zone_start == 0)
8622 for (dim = 0; dim < DIM; dim++)
8624 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8626 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8631 for (z = zone_start; z < zone_end; z++)
8633 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8635 zones->size[z].x0[XX], zones->size[z].x1[XX],
8636 zones->size[z].x0[YY], zones->size[z].x1[YY],
8637 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8638 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8640 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8641 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8642 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8647 static int comp_cgsort(const void *a, const void *b)
8651 gmx_cgsort_t *cga, *cgb;
8652 cga = (gmx_cgsort_t *)a;
8653 cgb = (gmx_cgsort_t *)b;
8655 comp = cga->nsc - cgb->nsc;
8658 comp = cga->ind_gl - cgb->ind_gl;
8664 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8669 /* Order the data */
8670 for (i = 0; i < n; i++)
8672 buf[i] = a[sort[i].ind];
8675 /* Copy back to the original array */
8676 for (i = 0; i < n; i++)
8682 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8687 /* Order the data */
8688 for (i = 0; i < n; i++)
8690 copy_rvec(v[sort[i].ind], buf[i]);
8693 /* Copy back to the original array */
8694 for (i = 0; i < n; i++)
8696 copy_rvec(buf[i], v[i]);
8700 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8703 int a, atot, cg, cg0, cg1, i;
8705 if (cgindex == NULL)
8707 /* Avoid the useless loop of the atoms within a cg */
8708 order_vec_cg(ncg, sort, v, buf);
8713 /* Order the data */
8715 for (cg = 0; cg < ncg; cg++)
8717 cg0 = cgindex[sort[cg].ind];
8718 cg1 = cgindex[sort[cg].ind+1];
8719 for (i = cg0; i < cg1; i++)
8721 copy_rvec(v[i], buf[a]);
8727 /* Copy back to the original array */
8728 for (a = 0; a < atot; a++)
8730 copy_rvec(buf[a], v[a]);
8734 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8735 int nsort_new, gmx_cgsort_t *sort_new,
8736 gmx_cgsort_t *sort1)
8740 /* The new indices are not very ordered, so we qsort them */
8741 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8743 /* sort2 is already ordered, so now we can merge the two arrays */
8747 while (i2 < nsort2 || i_new < nsort_new)
8751 sort1[i1++] = sort_new[i_new++];
8753 else if (i_new == nsort_new)
8755 sort1[i1++] = sort2[i2++];
8757 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8758 (sort2[i2].nsc == sort_new[i_new].nsc &&
8759 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8761 sort1[i1++] = sort2[i2++];
8765 sort1[i1++] = sort_new[i_new++];
8770 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8772 gmx_domdec_sort_t *sort;
8773 gmx_cgsort_t *cgsort, *sort_i;
8774 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8775 int sort_last, sort_skip;
8777 sort = dd->comm->sort;
8779 a = fr->ns.grid->cell_index;
8781 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8783 if (ncg_home_old >= 0)
8785 /* The charge groups that remained in the same ns grid cell
8786 * are completely ordered. So we can sort efficiently by sorting
8787 * the charge groups that did move into the stationary list.
8792 for (i = 0; i < dd->ncg_home; i++)
8794 /* Check if this cg did not move to another node */
8797 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8799 /* This cg is new on this node or moved ns grid cell */
8800 if (nsort_new >= sort->sort_new_nalloc)
8802 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8803 srenew(sort->sort_new, sort->sort_new_nalloc);
8805 sort_i = &(sort->sort_new[nsort_new++]);
8809 /* This cg did not move */
8810 sort_i = &(sort->sort2[nsort2++]);
8812 /* Sort on the ns grid cell indices
8813 * and the global topology index.
8814 * index_gl is irrelevant with cell ns,
8815 * but we set it here anyhow to avoid a conditional.
8818 sort_i->ind_gl = dd->index_gl[i];
8825 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8828 /* Sort efficiently */
8829 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8834 cgsort = sort->sort;
8836 for (i = 0; i < dd->ncg_home; i++)
8838 /* Sort on the ns grid cell indices
8839 * and the global topology index
8841 cgsort[i].nsc = a[i];
8842 cgsort[i].ind_gl = dd->index_gl[i];
8844 if (cgsort[i].nsc < moved)
8851 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8853 /* Determine the order of the charge groups using qsort */
8854 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8860 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8863 int ncg_new, i, *a, na;
8865 sort = dd->comm->sort->sort;
8867 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8870 for (i = 0; i < na; i++)
8874 sort[ncg_new].ind = a[i];
8882 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8883 rvec *cgcm, t_forcerec *fr, t_state *state,
8886 gmx_domdec_sort_t *sort;
8887 gmx_cgsort_t *cgsort, *sort_i;
8889 int ncg_new, i, *ibuf, cgsize;
8892 sort = dd->comm->sort;
8894 if (dd->ncg_home > sort->sort_nalloc)
8896 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8897 srenew(sort->sort, sort->sort_nalloc);
8898 srenew(sort->sort2, sort->sort_nalloc);
8900 cgsort = sort->sort;
8902 switch (fr->cutoff_scheme)
8905 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8908 ncg_new = dd_sort_order_nbnxn(dd, fr);
8911 gmx_incons("unimplemented");
8915 /* We alloc with the old size, since cgindex is still old */
8916 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8917 vbuf = dd->comm->vbuf.v;
8921 cgindex = dd->cgindex;
8928 /* Remove the charge groups which are no longer at home here */
8929 dd->ncg_home = ncg_new;
8932 fprintf(debug, "Set the new home charge group count to %d\n",
8936 /* Reorder the state */
8937 for (i = 0; i < estNR; i++)
8939 if (EST_DISTR(i) && (state->flags & (1<<i)))
8944 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8947 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8950 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8953 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8957 case estDISRE_INITF:
8958 case estDISRE_RM3TAV:
8959 case estORIRE_INITF:
8961 /* No ordering required */
8964 gmx_incons("Unknown state entry encountered in dd_sort_state");
8969 if (fr->cutoff_scheme == ecutsGROUP)
8972 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8975 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8977 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8978 srenew(sort->ibuf, sort->ibuf_nalloc);
8981 /* Reorder the global cg index */
8982 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
8983 /* Reorder the cginfo */
8984 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
8985 /* Rebuild the local cg index */
8989 for (i = 0; i < dd->ncg_home; i++)
8991 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8992 ibuf[i+1] = ibuf[i] + cgsize;
8994 for (i = 0; i < dd->ncg_home+1; i++)
8996 dd->cgindex[i] = ibuf[i];
9001 for (i = 0; i < dd->ncg_home+1; i++)
9006 /* Set the home atom number */
9007 dd->nat_home = dd->cgindex[dd->ncg_home];
9009 if (fr->cutoff_scheme == ecutsVERLET)
9011 /* The atoms are now exactly in grid order, update the grid order */
9012 nbnxn_set_atomorder(fr->nbv->nbs);
9016 /* Copy the sorted ns cell indices back to the ns grid struct */
9017 for (i = 0; i < dd->ncg_home; i++)
9019 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9021 fr->ns.grid->nr = dd->ncg_home;
9025 static void add_dd_statistics(gmx_domdec_t *dd)
9027 gmx_domdec_comm_t *comm;
9032 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9034 comm->sum_nat[ddnat-ddnatZONE] +=
9035 comm->nat[ddnat] - comm->nat[ddnat-1];
9040 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9042 gmx_domdec_comm_t *comm;
9047 /* Reset all the statistics and counters for total run counting */
9048 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9050 comm->sum_nat[ddnat-ddnatZONE] = 0;
9054 comm->load_step = 0;
9057 clear_ivec(comm->load_lim);
9062 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9064 gmx_domdec_comm_t *comm;
9068 comm = cr->dd->comm;
9070 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9077 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");
9079 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9081 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9086 " av. #atoms communicated per step for force: %d x %.1f\n",
9090 if (cr->dd->vsite_comm)
9093 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9094 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9099 if (cr->dd->constraint_comm)
9102 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9103 1 + ir->nLincsIter, av);
9107 gmx_incons(" Unknown type for DD statistics");
9110 fprintf(fplog, "\n");
9112 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9114 print_dd_load_av(fplog, cr->dd);
9118 void dd_partition_system(FILE *fplog,
9119 gmx_large_int_t step,
9121 gmx_bool bMasterState,
9123 t_state *state_global,
9124 gmx_mtop_t *top_global,
9126 t_state *state_local,
9129 gmx_localtop_t *top_local,
9132 gmx_shellfc_t shellfc,
9133 gmx_constr_t constr,
9135 gmx_wallcycle_t wcycle,
9139 gmx_domdec_comm_t *comm;
9140 gmx_ddbox_t ddbox = {0};
9142 gmx_large_int_t step_pcoupl;
9143 rvec cell_ns_x0, cell_ns_x1;
9144 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9145 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9146 gmx_bool bRedist, bSortCG, bResortAll;
9147 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9154 bBoxChanged = (bMasterState || DEFORM(*ir));
9155 if (ir->epc != epcNO)
9157 /* With nstpcouple > 1 pressure coupling happens.
9158 * one step after calculating the pressure.
9159 * Box scaling happens at the end of the MD step,
9160 * after the DD partitioning.
9161 * We therefore have to do DLB in the first partitioning
9162 * after an MD step where P-coupling occured.
9163 * We need to determine the last step in which p-coupling occurred.
9164 * MRS -- need to validate this for vv?
9169 step_pcoupl = step - 1;
9173 step_pcoupl = ((step - 1)/n)*n + 1;
9175 if (step_pcoupl >= comm->partition_step)
9181 bNStGlobalComm = (step % nstglobalcomm == 0);
9183 if (!comm->bDynLoadBal)
9189 /* Should we do dynamic load balacing this step?
9190 * Since it requires (possibly expensive) global communication,
9191 * we might want to do DLB less frequently.
9193 if (bBoxChanged || ir->epc != epcNO)
9195 bDoDLB = bBoxChanged;
9199 bDoDLB = bNStGlobalComm;
9203 /* Check if we have recorded loads on the nodes */
9204 if (comm->bRecordLoad && dd_load_count(comm))
9206 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9208 /* Check if we should use DLB at the second partitioning
9209 * and every 100 partitionings,
9210 * so the extra communication cost is negligible.
9212 n = max(100, nstglobalcomm);
9213 bCheckDLB = (comm->n_load_collect == 0 ||
9214 comm->n_load_have % n == n-1);
9221 /* Print load every nstlog, first and last step to the log file */
9222 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9223 comm->n_load_collect == 0 ||
9225 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9227 /* Avoid extra communication due to verbose screen output
9228 * when nstglobalcomm is set.
9230 if (bDoDLB || bLogLoad || bCheckDLB ||
9231 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9233 get_load_distribution(dd, wcycle);
9238 dd_print_load(fplog, dd, step-1);
9242 dd_print_load_verbose(dd);
9245 comm->n_load_collect++;
9249 /* Since the timings are node dependent, the master decides */
9253 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9256 fprintf(debug, "step %s, imb loss %f\n",
9257 gmx_step_str(step, sbuf),
9258 dd_force_imb_perf_loss(dd));
9261 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9264 turn_on_dlb(fplog, cr, step);
9269 comm->n_load_have++;
9272 cgs_gl = &comm->cgs_gl;
9277 /* Clear the old state */
9278 clear_dd_indices(dd, 0, 0);
9280 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9281 TRUE, cgs_gl, state_global->x, &ddbox);
9283 get_cg_distribution(fplog, step, dd, cgs_gl,
9284 state_global->box, &ddbox, state_global->x);
9286 dd_distribute_state(dd, cgs_gl,
9287 state_global, state_local, f);
9289 dd_make_local_cgs(dd, &top_local->cgs);
9291 /* Ensure that we have space for the new distribution */
9292 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9294 if (fr->cutoff_scheme == ecutsGROUP)
9296 calc_cgcm(fplog, 0, dd->ncg_home,
9297 &top_local->cgs, state_local->x, fr->cg_cm);
9300 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9302 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9306 else if (state_local->ddp_count != dd->ddp_count)
9308 if (state_local->ddp_count > dd->ddp_count)
9310 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9313 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9315 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);
9318 /* Clear the old state */
9319 clear_dd_indices(dd, 0, 0);
9321 /* Build the new indices */
9322 rebuild_cgindex(dd, cgs_gl->index, state_local);
9323 make_dd_indices(dd, cgs_gl->index, 0);
9325 if (fr->cutoff_scheme == ecutsGROUP)
9327 /* Redetermine the cg COMs */
9328 calc_cgcm(fplog, 0, dd->ncg_home,
9329 &top_local->cgs, state_local->x, fr->cg_cm);
9332 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9334 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9336 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9337 TRUE, &top_local->cgs, state_local->x, &ddbox);
9339 bRedist = comm->bDynLoadBal;
9343 /* We have the full state, only redistribute the cgs */
9345 /* Clear the non-home indices */
9346 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9348 /* Avoid global communication for dim's without pbc and -gcom */
9349 if (!bNStGlobalComm)
9351 copy_rvec(comm->box0, ddbox.box0 );
9352 copy_rvec(comm->box_size, ddbox.box_size);
9354 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9355 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9360 /* For dim's without pbc and -gcom */
9361 copy_rvec(ddbox.box0, comm->box0 );
9362 copy_rvec(ddbox.box_size, comm->box_size);
9364 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9367 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9369 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9372 /* Check if we should sort the charge groups */
9373 if (comm->nstSortCG > 0)
9375 bSortCG = (bMasterState ||
9376 (bRedist && (step % comm->nstSortCG == 0)));
9383 ncg_home_old = dd->ncg_home;
9388 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9390 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9391 state_local, f, fr, mdatoms,
9392 !bSortCG, nrnb, &cg0, &ncg_moved);
9394 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9397 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9399 &comm->cell_x0, &comm->cell_x1,
9400 dd->ncg_home, fr->cg_cm,
9401 cell_ns_x0, cell_ns_x1, &grid_density);
9405 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9408 switch (fr->cutoff_scheme)
9411 copy_ivec(fr->ns.grid->n, ncells_old);
9412 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9413 state_local->box, cell_ns_x0, cell_ns_x1,
9414 fr->rlistlong, grid_density);
9417 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9420 gmx_incons("unimplemented");
9422 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9423 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9427 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9429 /* Sort the state on charge group position.
9430 * This enables exact restarts from this step.
9431 * It also improves performance by about 15% with larger numbers
9432 * of atoms per node.
9435 /* Fill the ns grid with the home cell,
9436 * so we can sort with the indices.
9438 set_zones_ncg_home(dd);
9440 switch (fr->cutoff_scheme)
9443 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9445 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9447 comm->zones.size[0].bb_x0,
9448 comm->zones.size[0].bb_x1,
9450 comm->zones.dens_zone0,
9453 ncg_moved, bRedist ? comm->moved : NULL,
9454 fr->nbv->grp[eintLocal].kernel_type,
9455 fr->nbv->grp[eintLocal].nbat);
9457 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9460 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9461 0, dd->ncg_home, fr->cg_cm);
9463 copy_ivec(fr->ns.grid->n, ncells_new);
9466 gmx_incons("unimplemented");
9469 bResortAll = bMasterState;
9471 /* Check if we can user the old order and ns grid cell indices
9472 * of the charge groups to sort the charge groups efficiently.
9474 if (ncells_new[XX] != ncells_old[XX] ||
9475 ncells_new[YY] != ncells_old[YY] ||
9476 ncells_new[ZZ] != ncells_old[ZZ])
9483 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9484 gmx_step_str(step, sbuf), dd->ncg_home);
9486 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9487 bResortAll ? -1 : ncg_home_old);
9488 /* Rebuild all the indices */
9490 ga2la_clear(dd->ga2la);
9492 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9495 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9497 /* Setup up the communication and communicate the coordinates */
9498 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9500 /* Set the indices */
9501 make_dd_indices(dd, cgs_gl->index, cg0);
9503 /* Set the charge group boundaries for neighbor searching */
9504 set_cg_boundaries(&comm->zones);
9506 if (fr->cutoff_scheme == ecutsVERLET)
9508 set_zones_size(dd, state_local->box, &ddbox,
9509 bSortCG ? 1 : 0, comm->zones.n);
9512 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9515 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9516 -1,state_local->x,state_local->box);
9519 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9521 /* Extract a local topology from the global topology */
9522 for (i = 0; i < dd->ndim; i++)
9524 np[dd->dim[i]] = comm->cd[i].np;
9526 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9527 comm->cellsize_min, np,
9529 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9530 vsite, top_global, top_local);
9532 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9534 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9536 /* Set up the special atom communication */
9537 n = comm->nat[ddnatZONE];
9538 for (i = ddnatZONE+1; i < ddnatNR; i++)
9543 if (vsite && vsite->n_intercg_vsite)
9545 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9549 if (dd->bInterCGcons || dd->bInterCGsettles)
9551 /* Only for inter-cg constraints we need special code */
9552 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9553 constr, ir->nProjOrder,
9554 top_local->idef.il);
9558 gmx_incons("Unknown special atom type setup");
9563 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9565 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9567 /* Make space for the extra coordinates for virtual site
9568 * or constraint communication.
9570 state_local->natoms = comm->nat[ddnatNR-1];
9571 if (state_local->natoms > state_local->nalloc)
9573 dd_realloc_state(state_local, f, state_local->natoms);
9576 if (fr->bF_NoVirSum)
9578 if (vsite && vsite->n_intercg_vsite)
9580 nat_f_novirsum = comm->nat[ddnatVSITE];
9584 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9586 nat_f_novirsum = dd->nat_tot;
9590 nat_f_novirsum = dd->nat_home;
9599 /* Set the number of atoms required for the force calculation.
9600 * Forces need to be constrained when using a twin-range setup
9601 * or with energy minimization. For simple simulations we could
9602 * avoid some allocation, zeroing and copying, but this is
9603 * probably not worth the complications ande checking.
9605 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9606 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9608 /* We make the all mdatoms up to nat_tot_con.
9609 * We could save some work by only setting invmass
9610 * between nat_tot and nat_tot_con.
9612 /* This call also sets the new number of home particles to dd->nat_home */
9613 atoms2md(top_global, ir,
9614 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9616 /* Now we have the charges we can sort the FE interactions */
9617 dd_sort_local_top(dd, mdatoms, top_local);
9621 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9622 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9627 /* Make the local shell stuff, currently no communication is done */
9628 make_local_shells(cr, mdatoms, shellfc);
9631 if (ir->implicit_solvent)
9633 make_local_gb(cr, fr->born, ir->gb_algorithm);
9636 init_bonded_thread_force_reduction(fr, &top_local->idef);
9638 if (!(cr->duty & DUTY_PME))
9640 /* Send the charges to our PME only node */
9641 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9642 mdatoms->chargeA, mdatoms->chargeB,
9643 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9648 set_constraints(constr, top_local, ir, mdatoms, cr);
9651 if (ir->ePull != epullNO)
9653 /* Update the local pull groups */
9654 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9659 /* Update the local rotation groups */
9660 dd_make_local_rotation_groups(dd, ir->rot);
9664 add_dd_statistics(dd);
9666 /* Make sure we only count the cycles for this DD partitioning */
9667 clear_dd_cycle_counts(dd);
9669 /* Because the order of the atoms might have changed since
9670 * the last vsite construction, we need to communicate the constructing
9671 * atom coordinates again (for spreading the forces this MD step).
9673 dd_move_x_vsites(dd, state_local->box, state_local->x);
9675 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9677 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9679 dd_move_x(dd, state_local->box, state_local->x);
9680 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9681 -1, state_local->x, state_local->box);
9684 /* Store the partitioning step */
9685 comm->partition_step = step;
9687 /* Increase the DD partitioning counter */
9689 /* The state currently matches this DD partitioning count, store it */
9690 state_local->ddp_count = dd->ddp_count;
9693 /* The DD master node knows the complete cg distribution,
9694 * store the count so we can possibly skip the cg info communication.
9696 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9699 if (comm->DD_debug > 0)
9701 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9702 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9703 "after partitioning");