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"
44 #include "pull_rotation.h"
48 #include "mtop_util.h"
49 #include "gmx_ga2la.h"
52 #include "nbnxn_search.h"
54 #include "gmx_omp_nthreads.h"
56 #include "gromacs/fileio/futil.h"
57 #include "gromacs/fileio/gmxfio.h"
58 #include "gromacs/fileio/pdbio.h"
59 #include "gromacs/timing/wallcycle.h"
60 #include "gromacs/utility/gmxmpi.h"
62 #define DDRANK(dd, rank) (rank)
63 #define DDMASTERRANK(dd) (dd->masterrank)
65 typedef struct gmx_domdec_master
67 /* The cell boundaries */
69 /* The global charge group division */
70 int *ncg; /* Number of home charge groups for each node */
71 int *index; /* Index of nnodes+1 into cg */
72 int *cg; /* Global charge group index */
73 int *nat; /* Number of home atoms for each node. */
74 int *ibuf; /* Buffer for communication */
75 rvec *vbuf; /* Buffer for state scattering and gathering */
76 } gmx_domdec_master_t;
80 /* The numbers of charge groups to send and receive for each cell
81 * that requires communication, the last entry contains the total
82 * number of atoms that needs to be communicated.
84 int nsend[DD_MAXIZONE+2];
85 int nrecv[DD_MAXIZONE+2];
86 /* The charge groups to send */
89 /* The atom range for non-in-place communication */
90 int cell2at0[DD_MAXIZONE];
91 int cell2at1[DD_MAXIZONE];
96 int np; /* Number of grid pulses in this dimension */
97 int np_dlb; /* For dlb, for use with edlbAUTO */
98 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
100 gmx_bool bInPlace; /* Can we communicate in place? */
101 } gmx_domdec_comm_dim_t;
105 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
106 real *cell_f; /* State var.: cell boundaries, box relative */
107 real *old_cell_f; /* Temp. var.: old cell size */
108 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
109 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
110 real *bound_min; /* Temp. var.: lower limit for cell boundary */
111 real *bound_max; /* Temp. var.: upper limit for cell boundary */
112 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
113 real *buf_ncd; /* Temp. var. */
116 #define DD_NLOAD_MAX 9
118 /* Here floats are accurate enough, since these variables
119 * only influence the load balancing, not the actual MD results.
146 gmx_cgsort_t *sort_new;
158 /* This enum determines the order of the coordinates.
159 * ddnatHOME and ddnatZONE should be first and second,
160 * the others can be ordered as wanted.
163 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
167 edlbAUTO, edlbNO, edlbYES, edlbNR
169 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
173 int dim; /* The dimension */
174 gmx_bool dim_match; /* Tells if DD and PME dims match */
175 int nslab; /* The number of PME slabs in this dimension */
176 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
177 int *pp_min; /* The minimum pp node location, size nslab */
178 int *pp_max; /* The maximum pp node location,size nslab */
179 int maxshift; /* The maximum shift for coordinate redistribution in PME */
184 real min0; /* The minimum bottom of this zone */
185 real max1; /* The maximum top of this zone */
186 real min1; /* The minimum top of this zone */
187 real mch0; /* The maximum bottom communicaton height for this zone */
188 real mch1; /* The maximum top communicaton height for this zone */
189 real p1_0; /* The bottom value of the first cell in this zone */
190 real p1_1; /* The top value of the first cell in this zone */
195 gmx_domdec_ind_t ind;
202 } dd_comm_setup_work_t;
204 typedef struct gmx_domdec_comm
206 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
207 * unless stated otherwise.
210 /* The number of decomposition dimensions for PME, 0: no PME */
212 /* The number of nodes doing PME (PP/PME or only PME) */
216 /* The communication setup including the PME only nodes */
217 gmx_bool bCartesianPP_PME;
220 int *pmenodes; /* size npmenodes */
221 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
222 * but with bCartesianPP_PME */
223 gmx_ddpme_t ddpme[2];
225 /* The DD particle-particle nodes only */
226 gmx_bool bCartesianPP;
227 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
229 /* The global charge groups */
232 /* Should we sort the cgs */
234 gmx_domdec_sort_t *sort;
236 /* Are there charge groups? */
239 /* Are there bonded and multi-body interactions between charge groups? */
240 gmx_bool bInterCGBondeds;
241 gmx_bool bInterCGMultiBody;
243 /* Data for the optional bonded interaction atom communication range */
250 /* Are we actually using DLB? */
251 gmx_bool bDynLoadBal;
253 /* Cell sizes for static load balancing, first index cartesian */
256 /* The width of the communicated boundaries */
259 /* The minimum cell size (including triclinic correction) */
261 /* For dlb, for use with edlbAUTO */
262 rvec cellsize_min_dlb;
263 /* The lower limit for the DD cell size with DLB */
265 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
266 gmx_bool bVacDLBNoLimit;
268 /* With PME load balancing we set limits on DLB */
269 gmx_bool bPMELoadBalDLBLimits;
270 /* DLB needs to take into account that we want to allow this maximum
271 * cut-off (for PME load balancing), this could limit cell boundaries.
273 real PMELoadBal_max_cutoff;
275 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
277 /* box0 and box_size are required with dim's without pbc and -gcom */
281 /* The cell boundaries */
285 /* The old location of the cell boundaries, to check cg displacements */
289 /* The communication setup and charge group boundaries for the zones */
290 gmx_domdec_zones_t zones;
292 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
293 * cell boundaries of neighboring cells for dynamic load balancing.
295 gmx_ddzone_t zone_d1[2];
296 gmx_ddzone_t zone_d2[2][2];
298 /* The coordinate/force communication setup and indices */
299 gmx_domdec_comm_dim_t cd[DIM];
300 /* The maximum number of cells to communicate with in one dimension */
303 /* Which cg distribution is stored on the master node */
304 int master_cg_ddp_count;
306 /* The number of cg's received from the direct neighbors */
307 int zone_ncg1[DD_MAXZONE];
309 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
312 /* Array for signalling if atoms have moved to another domain */
316 /* Communication buffer for general use */
320 /* Communication buffer for general use */
323 /* Temporary storage for thread parallel communication setup */
325 dd_comm_setup_work_t *dth;
327 /* Communication buffers only used with multiple grid pulses */
332 /* Communication buffers for local redistribution */
334 int cggl_flag_nalloc[DIM*2];
336 int cgcm_state_nalloc[DIM*2];
338 /* Cell sizes for dynamic load balancing */
339 gmx_domdec_root_t **root;
343 real cell_f_max0[DIM];
344 real cell_f_min1[DIM];
346 /* Stuff for load communication */
347 gmx_bool bRecordLoad;
348 gmx_domdec_load_t *load;
350 MPI_Comm *mpi_comm_load;
353 /* Maximum DLB scaling per load balancing step in percent */
357 float cycl[ddCyclNr];
358 int cycl_n[ddCyclNr];
359 float cycl_max[ddCyclNr];
360 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
364 /* Have often have did we have load measurements */
366 /* Have often have we collected the load measurements */
370 double sum_nat[ddnatNR-ddnatZONE];
380 /* The last partition step */
381 gmx_large_int_t partition_step;
389 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
392 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
393 #define DD_FLAG_NRCG 65535
394 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
395 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
397 /* Zone permutation required to obtain consecutive charge groups
398 * for neighbor searching.
400 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
402 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
403 * components see only j zones with that component 0.
406 /* The DD zone order */
407 static const ivec dd_zo[DD_MAXZONE] =
408 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
413 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
418 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
423 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
425 /* Factors used to avoid problems due to rounding issues */
426 #define DD_CELL_MARGIN 1.0001
427 #define DD_CELL_MARGIN2 1.00005
428 /* Factor to account for pressure scaling during nstlist steps */
429 #define DD_PRES_SCALE_MARGIN 1.02
431 /* Allowed performance loss before we DLB or warn */
432 #define DD_PERF_LOSS 0.05
434 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
436 /* Use separate MPI send and receive commands
437 * when nnodes <= GMX_DD_NNODES_SENDRECV.
438 * This saves memory (and some copying for small nnodes).
439 * For high parallelization scatter and gather calls are used.
441 #define GMX_DD_NNODES_SENDRECV 4
445 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
447 static void index2xyz(ivec nc,int ind,ivec xyz)
449 xyz[XX] = ind % nc[XX];
450 xyz[YY] = (ind / nc[XX]) % nc[YY];
451 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
455 /* This order is required to minimize the coordinate communication in PME
456 * which uses decomposition in the x direction.
458 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
460 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
462 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
463 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
464 xyz[ZZ] = ind % nc[ZZ];
467 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
472 ddindex = dd_index(dd->nc, c);
473 if (dd->comm->bCartesianPP_PME)
475 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
477 else if (dd->comm->bCartesianPP)
480 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
491 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
493 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
496 int ddglatnr(gmx_domdec_t *dd, int i)
506 if (i >= dd->comm->nat[ddnatNR-1])
508 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
510 atnr = dd->gatindex[i] + 1;
516 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
518 return &dd->comm->cgs_gl;
521 static void vec_rvec_init(vec_rvec_t *v)
527 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
531 v->nalloc = over_alloc_dd(n);
532 srenew(v->v, v->nalloc);
536 void dd_store_state(gmx_domdec_t *dd, t_state *state)
540 if (state->ddp_count != dd->ddp_count)
542 gmx_incons("The state does not the domain decomposition state");
545 state->ncg_gl = dd->ncg_home;
546 if (state->ncg_gl > state->cg_gl_nalloc)
548 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
549 srenew(state->cg_gl, state->cg_gl_nalloc);
551 for (i = 0; i < state->ncg_gl; i++)
553 state->cg_gl[i] = dd->index_gl[i];
556 state->ddp_count_cg_gl = dd->ddp_count;
559 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
561 return &dd->comm->zones;
564 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
565 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
567 gmx_domdec_zones_t *zones;
570 zones = &dd->comm->zones;
573 while (icg >= zones->izone[izone].cg1)
582 else if (izone < zones->nizone)
584 *jcg0 = zones->izone[izone].jcg0;
588 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
589 icg, izone, zones->nizone);
592 *jcg1 = zones->izone[izone].jcg1;
594 for (d = 0; d < dd->ndim; d++)
597 shift0[dim] = zones->izone[izone].shift0[dim];
598 shift1[dim] = zones->izone[izone].shift1[dim];
599 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
601 /* A conservative approach, this can be optimized */
608 int dd_natoms_vsite(gmx_domdec_t *dd)
610 return dd->comm->nat[ddnatVSITE];
613 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
615 *at_start = dd->comm->nat[ddnatCON-1];
616 *at_end = dd->comm->nat[ddnatCON];
619 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
621 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
622 int *index, *cgindex;
623 gmx_domdec_comm_t *comm;
624 gmx_domdec_comm_dim_t *cd;
625 gmx_domdec_ind_t *ind;
626 rvec shift = {0, 0, 0}, *buf, *rbuf;
627 gmx_bool bPBC, bScrew;
631 cgindex = dd->cgindex;
636 nat_tot = dd->nat_home;
637 for (d = 0; d < dd->ndim; d++)
639 bPBC = (dd->ci[dd->dim[d]] == 0);
640 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
643 copy_rvec(box[dd->dim[d]], shift);
646 for (p = 0; p < cd->np; p++)
653 for (i = 0; i < ind->nsend[nzone]; i++)
655 at0 = cgindex[index[i]];
656 at1 = cgindex[index[i]+1];
657 for (j = at0; j < at1; j++)
659 copy_rvec(x[j], buf[n]);
666 for (i = 0; i < ind->nsend[nzone]; i++)
668 at0 = cgindex[index[i]];
669 at1 = cgindex[index[i]+1];
670 for (j = at0; j < at1; j++)
672 /* We need to shift the coordinates */
673 rvec_add(x[j], shift, buf[n]);
680 for (i = 0; i < ind->nsend[nzone]; i++)
682 at0 = cgindex[index[i]];
683 at1 = cgindex[index[i]+1];
684 for (j = at0; j < at1; j++)
687 buf[n][XX] = x[j][XX] + shift[XX];
689 * This operation requires a special shift force
690 * treatment, which is performed in calc_vir.
692 buf[n][YY] = box[YY][YY] - x[j][YY];
693 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
705 rbuf = comm->vbuf2.v;
707 /* Send and receive the coordinates */
708 dd_sendrecv_rvec(dd, d, dddirBackward,
709 buf, ind->nsend[nzone+1],
710 rbuf, ind->nrecv[nzone+1]);
714 for (zone = 0; zone < nzone; zone++)
716 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
718 copy_rvec(rbuf[j], x[i]);
723 nat_tot += ind->nrecv[nzone+1];
729 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
731 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
732 int *index, *cgindex;
733 gmx_domdec_comm_t *comm;
734 gmx_domdec_comm_dim_t *cd;
735 gmx_domdec_ind_t *ind;
739 gmx_bool bPBC, bScrew;
743 cgindex = dd->cgindex;
748 nzone = comm->zones.n/2;
749 nat_tot = dd->nat_tot;
750 for (d = dd->ndim-1; d >= 0; d--)
752 bPBC = (dd->ci[dd->dim[d]] == 0);
753 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
754 if (fshift == NULL && !bScrew)
758 /* Determine which shift vector we need */
764 for (p = cd->np-1; p >= 0; p--)
767 nat_tot -= ind->nrecv[nzone+1];
774 sbuf = comm->vbuf2.v;
776 for (zone = 0; zone < nzone; zone++)
778 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
780 copy_rvec(f[i], sbuf[j]);
785 /* Communicate the forces */
786 dd_sendrecv_rvec(dd, d, dddirForward,
787 sbuf, ind->nrecv[nzone+1],
788 buf, ind->nsend[nzone+1]);
790 /* Add the received forces */
794 for (i = 0; i < ind->nsend[nzone]; i++)
796 at0 = cgindex[index[i]];
797 at1 = cgindex[index[i]+1];
798 for (j = at0; j < at1; j++)
800 rvec_inc(f[j], buf[n]);
807 for (i = 0; i < ind->nsend[nzone]; i++)
809 at0 = cgindex[index[i]];
810 at1 = cgindex[index[i]+1];
811 for (j = at0; j < at1; j++)
813 rvec_inc(f[j], buf[n]);
814 /* Add this force to the shift force */
815 rvec_inc(fshift[is], buf[n]);
822 for (i = 0; i < ind->nsend[nzone]; i++)
824 at0 = cgindex[index[i]];
825 at1 = cgindex[index[i]+1];
826 for (j = at0; j < at1; j++)
828 /* Rotate the force */
829 f[j][XX] += buf[n][XX];
830 f[j][YY] -= buf[n][YY];
831 f[j][ZZ] -= buf[n][ZZ];
834 /* Add this force to the shift force */
835 rvec_inc(fshift[is], buf[n]);
846 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
848 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
849 int *index, *cgindex;
850 gmx_domdec_comm_t *comm;
851 gmx_domdec_comm_dim_t *cd;
852 gmx_domdec_ind_t *ind;
857 cgindex = dd->cgindex;
859 buf = &comm->vbuf.v[0][0];
862 nat_tot = dd->nat_home;
863 for (d = 0; d < dd->ndim; d++)
866 for (p = 0; p < cd->np; p++)
871 for (i = 0; i < ind->nsend[nzone]; i++)
873 at0 = cgindex[index[i]];
874 at1 = cgindex[index[i]+1];
875 for (j = at0; j < at1; j++)
888 rbuf = &comm->vbuf2.v[0][0];
890 /* Send and receive the coordinates */
891 dd_sendrecv_real(dd, d, dddirBackward,
892 buf, ind->nsend[nzone+1],
893 rbuf, ind->nrecv[nzone+1]);
897 for (zone = 0; zone < nzone; zone++)
899 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
906 nat_tot += ind->nrecv[nzone+1];
912 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
914 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
915 int *index, *cgindex;
916 gmx_domdec_comm_t *comm;
917 gmx_domdec_comm_dim_t *cd;
918 gmx_domdec_ind_t *ind;
923 cgindex = dd->cgindex;
925 buf = &comm->vbuf.v[0][0];
928 nzone = comm->zones.n/2;
929 nat_tot = dd->nat_tot;
930 for (d = dd->ndim-1; d >= 0; d--)
933 for (p = cd->np-1; p >= 0; p--)
936 nat_tot -= ind->nrecv[nzone+1];
943 sbuf = &comm->vbuf2.v[0][0];
945 for (zone = 0; zone < nzone; zone++)
947 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
954 /* Communicate the forces */
955 dd_sendrecv_real(dd, d, dddirForward,
956 sbuf, ind->nrecv[nzone+1],
957 buf, ind->nsend[nzone+1]);
959 /* Add the received forces */
961 for (i = 0; i < ind->nsend[nzone]; i++)
963 at0 = cgindex[index[i]];
964 at1 = cgindex[index[i]+1];
965 for (j = at0; j < at1; j++)
976 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
978 fprintf(fp, "zone d0 %d d1 %d d2 %d min0 %6.3f max1 %6.3f mch0 %6.3f mch1 %6.3f p1_0 %6.3f p1_1 %6.3f\n",
980 zone->min0, zone->max1,
981 zone->mch0, zone->mch0,
982 zone->p1_0, zone->p1_1);
986 #define DDZONECOMM_MAXZONE 5
987 #define DDZONECOMM_BUFSIZE 3
989 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
990 int ddimind, int direction,
991 gmx_ddzone_t *buf_s, int n_s,
992 gmx_ddzone_t *buf_r, int n_r)
994 #define ZBS DDZONECOMM_BUFSIZE
995 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
996 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
999 for (i = 0; i < n_s; i++)
1001 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1002 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1003 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1004 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1005 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1006 vbuf_s[i*ZBS+1][2] = 0;
1007 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1008 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1009 vbuf_s[i*ZBS+2][2] = 0;
1012 dd_sendrecv_rvec(dd, ddimind, direction,
1016 for (i = 0; i < n_r; i++)
1018 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1019 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1020 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1021 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1022 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1023 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1024 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1030 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1031 rvec cell_ns_x0, rvec cell_ns_x1)
1033 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1035 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1036 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1037 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1038 rvec extr_s[2], extr_r[2];
1040 real dist_d, c = 0, det;
1041 gmx_domdec_comm_t *comm;
1042 gmx_bool bPBC, bUse;
1046 for (d = 1; d < dd->ndim; d++)
1049 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1050 zp->min0 = cell_ns_x0[dim];
1051 zp->max1 = cell_ns_x1[dim];
1052 zp->min1 = cell_ns_x1[dim];
1053 zp->mch0 = cell_ns_x0[dim];
1054 zp->mch1 = cell_ns_x1[dim];
1055 zp->p1_0 = cell_ns_x0[dim];
1056 zp->p1_1 = cell_ns_x1[dim];
1059 for (d = dd->ndim-2; d >= 0; d--)
1062 bPBC = (dim < ddbox->npbcdim);
1064 /* Use an rvec to store two reals */
1065 extr_s[d][0] = comm->cell_f0[d+1];
1066 extr_s[d][1] = comm->cell_f1[d+1];
1067 extr_s[d][2] = comm->cell_f1[d+1];
1070 /* Store the extremes in the backward sending buffer,
1071 * so the get updated separately from the forward communication.
1073 for (d1 = d; d1 < dd->ndim-1; d1++)
1075 /* We invert the order to be able to use the same loop for buf_e */
1076 buf_s[pos].min0 = extr_s[d1][1];
1077 buf_s[pos].max1 = extr_s[d1][0];
1078 buf_s[pos].min1 = extr_s[d1][2];
1079 buf_s[pos].mch0 = 0;
1080 buf_s[pos].mch1 = 0;
1081 /* Store the cell corner of the dimension we communicate along */
1082 buf_s[pos].p1_0 = comm->cell_x0[dim];
1083 buf_s[pos].p1_1 = 0;
1087 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1090 if (dd->ndim == 3 && d == 0)
1092 buf_s[pos] = comm->zone_d2[0][1];
1094 buf_s[pos] = comm->zone_d1[0];
1098 /* We only need to communicate the extremes
1099 * in the forward direction
1101 npulse = comm->cd[d].np;
1104 /* Take the minimum to avoid double communication */
1105 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1109 /* Without PBC we should really not communicate over
1110 * the boundaries, but implementing that complicates
1111 * the communication setup and therefore we simply
1112 * do all communication, but ignore some data.
1114 npulse_min = npulse;
1116 for (p = 0; p < npulse_min; p++)
1118 /* Communicate the extremes forward */
1119 bUse = (bPBC || dd->ci[dim] > 0);
1121 dd_sendrecv_rvec(dd, d, dddirForward,
1122 extr_s+d, dd->ndim-d-1,
1123 extr_r+d, dd->ndim-d-1);
1127 for (d1 = d; d1 < dd->ndim-1; d1++)
1129 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1130 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1131 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1137 for (p = 0; p < npulse; p++)
1139 /* Communicate all the zone information backward */
1140 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1142 dd_sendrecv_ddzone(dd, d, dddirBackward,
1149 for (d1 = d+1; d1 < dd->ndim; d1++)
1151 /* Determine the decrease of maximum required
1152 * communication height along d1 due to the distance along d,
1153 * this avoids a lot of useless atom communication.
1155 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1157 if (ddbox->tric_dir[dim])
1159 /* c is the off-diagonal coupling between the cell planes
1160 * along directions d and d1.
1162 c = ddbox->v[dim][dd->dim[d1]][dim];
1168 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1171 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1175 /* A negative value signals out of range */
1181 /* Accumulate the extremes over all pulses */
1182 for (i = 0; i < buf_size; i++)
1186 buf_e[i] = buf_r[i];
1192 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1193 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1194 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1197 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1205 if (bUse && dh[d1] >= 0)
1207 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1208 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1211 /* Copy the received buffer to the send buffer,
1212 * to pass the data through with the next pulse.
1214 buf_s[i] = buf_r[i];
1216 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1217 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1219 /* Store the extremes */
1222 for (d1 = d; d1 < dd->ndim-1; d1++)
1224 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1225 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1226 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1230 if (d == 1 || (d == 0 && dd->ndim == 3))
1232 for (i = d; i < 2; i++)
1234 comm->zone_d2[1-d][i] = buf_e[pos];
1240 comm->zone_d1[1] = buf_e[pos];
1250 for (i = 0; i < 2; i++)
1254 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1256 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1257 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1263 for (i = 0; i < 2; i++)
1265 for (j = 0; j < 2; j++)
1269 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1271 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1272 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1276 for (d = 1; d < dd->ndim; d++)
1278 comm->cell_f_max0[d] = extr_s[d-1][0];
1279 comm->cell_f_min1[d] = extr_s[d-1][1];
1282 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1283 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1288 static void dd_collect_cg(gmx_domdec_t *dd,
1289 t_state *state_local)
1291 gmx_domdec_master_t *ma = NULL;
1292 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1295 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1297 /* The master has the correct distribution */
1301 if (state_local->ddp_count == dd->ddp_count)
1303 ncg_home = dd->ncg_home;
1305 nat_home = dd->nat_home;
1307 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1309 cgs_gl = &dd->comm->cgs_gl;
1311 ncg_home = state_local->ncg_gl;
1312 cg = state_local->cg_gl;
1314 for (i = 0; i < ncg_home; i++)
1316 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1321 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1324 buf2[0] = dd->ncg_home;
1325 buf2[1] = dd->nat_home;
1335 /* Collect the charge group and atom counts on the master */
1336 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1341 for (i = 0; i < dd->nnodes; i++)
1343 ma->ncg[i] = ma->ibuf[2*i];
1344 ma->nat[i] = ma->ibuf[2*i+1];
1345 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1348 /* Make byte counts and indices */
1349 for (i = 0; i < dd->nnodes; i++)
1351 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1352 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1356 fprintf(debug, "Initial charge group distribution: ");
1357 for (i = 0; i < dd->nnodes; i++)
1359 fprintf(debug, " %d", ma->ncg[i]);
1361 fprintf(debug, "\n");
1365 /* Collect the charge group indices on the master */
1367 dd->ncg_home*sizeof(int), dd->index_gl,
1368 DDMASTER(dd) ? ma->ibuf : NULL,
1369 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1370 DDMASTER(dd) ? ma->cg : NULL);
1372 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1375 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1378 gmx_domdec_master_t *ma;
1379 int n, i, c, a, nalloc = 0;
1388 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1389 dd->rank, dd->mpi_comm_all);
1394 /* Copy the master coordinates to the global array */
1395 cgs_gl = &dd->comm->cgs_gl;
1397 n = DDMASTERRANK(dd);
1399 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1401 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1403 copy_rvec(lv[a++], v[c]);
1407 for (n = 0; n < dd->nnodes; n++)
1411 if (ma->nat[n] > nalloc)
1413 nalloc = over_alloc_dd(ma->nat[n]);
1414 srenew(buf, nalloc);
1417 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1418 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1421 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1423 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1425 copy_rvec(buf[a++], v[c]);
1434 static void get_commbuffer_counts(gmx_domdec_t *dd,
1435 int **counts, int **disps)
1437 gmx_domdec_master_t *ma;
1442 /* Make the rvec count and displacment arrays */
1444 *disps = ma->ibuf + dd->nnodes;
1445 for (n = 0; n < dd->nnodes; n++)
1447 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1448 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1452 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1455 gmx_domdec_master_t *ma;
1456 int *rcounts = NULL, *disps = NULL;
1465 get_commbuffer_counts(dd, &rcounts, &disps);
1470 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1474 cgs_gl = &dd->comm->cgs_gl;
1477 for (n = 0; n < dd->nnodes; n++)
1479 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1481 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1483 copy_rvec(buf[a++], v[c]);
1490 void dd_collect_vec(gmx_domdec_t *dd,
1491 t_state *state_local, rvec *lv, rvec *v)
1493 gmx_domdec_master_t *ma;
1494 int n, i, c, a, nalloc = 0;
1497 dd_collect_cg(dd, state_local);
1499 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1501 dd_collect_vec_sendrecv(dd, lv, v);
1505 dd_collect_vec_gatherv(dd, lv, v);
1510 void dd_collect_state(gmx_domdec_t *dd,
1511 t_state *state_local, t_state *state)
1515 nh = state->nhchainlength;
1519 for (i = 0; i < efptNR; i++)
1521 state->lambda[i] = state_local->lambda[i];
1523 state->fep_state = state_local->fep_state;
1524 state->veta = state_local->veta;
1525 state->vol0 = state_local->vol0;
1526 copy_mat(state_local->box, state->box);
1527 copy_mat(state_local->boxv, state->boxv);
1528 copy_mat(state_local->svir_prev, state->svir_prev);
1529 copy_mat(state_local->fvir_prev, state->fvir_prev);
1530 copy_mat(state_local->pres_prev, state->pres_prev);
1532 for (i = 0; i < state_local->ngtc; i++)
1534 for (j = 0; j < nh; j++)
1536 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1537 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1539 state->therm_integral[i] = state_local->therm_integral[i];
1541 for (i = 0; i < state_local->nnhpres; i++)
1543 for (j = 0; j < nh; j++)
1545 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1546 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1550 for (est = 0; est < estNR; est++)
1552 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1557 dd_collect_vec(dd, state_local, state_local->x, state->x);
1560 dd_collect_vec(dd, state_local, state_local->v, state->v);
1563 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1566 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1569 if (state->nrngi == 1)
1573 for (i = 0; i < state_local->nrng; i++)
1575 state->ld_rng[i] = state_local->ld_rng[i];
1581 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1582 state_local->ld_rng, state->ld_rng);
1586 if (state->nrngi == 1)
1590 state->ld_rngi[0] = state_local->ld_rngi[0];
1595 dd_gather(dd, sizeof(state->ld_rngi[0]),
1596 state_local->ld_rngi, state->ld_rngi);
1599 case estDISRE_INITF:
1600 case estDISRE_RM3TAV:
1601 case estORIRE_INITF:
1605 gmx_incons("Unknown state entry encountered in dd_collect_state");
1611 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1617 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1620 state->nalloc = over_alloc_dd(nalloc);
1622 for (est = 0; est < estNR; est++)
1624 if (EST_DISTR(est) && (state->flags & (1<<est)))
1629 srenew(state->x, state->nalloc);
1632 srenew(state->v, state->nalloc);
1635 srenew(state->sd_X, state->nalloc);
1638 srenew(state->cg_p, state->nalloc);
1642 case estDISRE_INITF:
1643 case estDISRE_RM3TAV:
1644 case estORIRE_INITF:
1646 /* No reallocation required */
1649 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1656 srenew(*f, state->nalloc);
1660 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1663 if (nalloc > fr->cg_nalloc)
1667 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1669 fr->cg_nalloc = over_alloc_dd(nalloc);
1670 srenew(fr->cginfo, fr->cg_nalloc);
1671 if (fr->cutoff_scheme == ecutsGROUP)
1673 srenew(fr->cg_cm, fr->cg_nalloc);
1676 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1678 /* We don't use charge groups, we use x in state to set up
1679 * the atom communication.
1681 dd_realloc_state(state, f, nalloc);
1685 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1688 gmx_domdec_master_t *ma;
1689 int n, i, c, a, nalloc = 0;
1696 for (n = 0; n < dd->nnodes; n++)
1700 if (ma->nat[n] > nalloc)
1702 nalloc = over_alloc_dd(ma->nat[n]);
1703 srenew(buf, nalloc);
1705 /* Use lv as a temporary buffer */
1707 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1709 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1711 copy_rvec(v[c], buf[a++]);
1714 if (a != ma->nat[n])
1716 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1721 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1722 DDRANK(dd, n), n, dd->mpi_comm_all);
1727 n = DDMASTERRANK(dd);
1729 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1731 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1733 copy_rvec(v[c], lv[a++]);
1740 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1741 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1746 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1749 gmx_domdec_master_t *ma;
1750 int *scounts = NULL, *disps = NULL;
1751 int n, i, c, a, nalloc = 0;
1758 get_commbuffer_counts(dd, &scounts, &disps);
1762 for (n = 0; n < dd->nnodes; n++)
1764 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1766 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1768 copy_rvec(v[c], buf[a++]);
1774 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1777 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1779 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1781 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1785 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1789 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1792 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1793 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1794 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1796 if (dfhist->nlambda > 0)
1798 int nlam = dfhist->nlambda;
1799 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1800 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1801 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1802 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1803 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1804 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1806 for (i = 0; i<nlam; i++) {
1807 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1808 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1809 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1810 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1811 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1812 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1817 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1818 t_state *state, t_state *state_local,
1823 nh = state->nhchainlength;
1827 for (i = 0; i < efptNR; i++)
1829 state_local->lambda[i] = state->lambda[i];
1831 state_local->fep_state = state->fep_state;
1832 state_local->veta = state->veta;
1833 state_local->vol0 = state->vol0;
1834 copy_mat(state->box, state_local->box);
1835 copy_mat(state->box_rel, state_local->box_rel);
1836 copy_mat(state->boxv, state_local->boxv);
1837 copy_mat(state->svir_prev, state_local->svir_prev);
1838 copy_mat(state->fvir_prev, state_local->fvir_prev);
1839 copy_df_history(&state_local->dfhist,&state->dfhist);
1840 for (i = 0; i < state_local->ngtc; i++)
1842 for (j = 0; j < nh; j++)
1844 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1845 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1847 state_local->therm_integral[i] = state->therm_integral[i];
1849 for (i = 0; i < state_local->nnhpres; i++)
1851 for (j = 0; j < nh; j++)
1853 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1854 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1858 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1859 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1860 dd_bcast(dd, sizeof(real), &state_local->veta);
1861 dd_bcast(dd, sizeof(real), &state_local->vol0);
1862 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1863 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1864 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1865 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1866 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1867 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1868 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1869 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1870 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1871 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1873 /* communicate df_history -- required for restarting from checkpoint */
1874 dd_distribute_dfhist(dd,&state_local->dfhist);
1876 if (dd->nat_home > state_local->nalloc)
1878 dd_realloc_state(state_local, f, dd->nat_home);
1880 for (i = 0; i < estNR; i++)
1882 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1887 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1890 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1893 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1896 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1899 if (state->nrngi == 1)
1902 state_local->nrng*sizeof(state_local->ld_rng[0]),
1903 state->ld_rng, state_local->ld_rng);
1908 state_local->nrng*sizeof(state_local->ld_rng[0]),
1909 state->ld_rng, state_local->ld_rng);
1913 if (state->nrngi == 1)
1915 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1916 state->ld_rngi, state_local->ld_rngi);
1920 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1921 state->ld_rngi, state_local->ld_rngi);
1924 case estDISRE_INITF:
1925 case estDISRE_RM3TAV:
1926 case estORIRE_INITF:
1928 /* Not implemented yet */
1931 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1937 static char dim2char(int dim)
1943 case XX: c = 'X'; break;
1944 case YY: c = 'Y'; break;
1945 case ZZ: c = 'Z'; break;
1946 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1952 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1953 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1955 rvec grid_s[2], *grid_r = NULL, cx, r;
1956 char fname[STRLEN], format[STRLEN], buf[22];
1958 int a, i, d, z, y, x;
1962 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1963 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1967 snew(grid_r, 2*dd->nnodes);
1970 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1974 for (d = 0; d < DIM; d++)
1976 for (i = 0; i < DIM; i++)
1984 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1986 tric[d][i] = box[i][d]/box[i][i];
1995 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1996 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
1997 out = gmx_fio_fopen(fname, "w");
1998 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2000 for (i = 0; i < dd->nnodes; i++)
2002 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2003 for (d = 0; d < DIM; d++)
2005 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2007 for (z = 0; z < 2; z++)
2009 for (y = 0; y < 2; y++)
2011 for (x = 0; x < 2; x++)
2013 cx[XX] = grid_r[i*2+x][XX];
2014 cx[YY] = grid_r[i*2+y][YY];
2015 cx[ZZ] = grid_r[i*2+z][ZZ];
2017 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
2018 ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
2022 for (d = 0; d < DIM; d++)
2024 for (x = 0; x < 4; x++)
2028 case 0: y = 1 + i*8 + 2*x; break;
2029 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2030 case 2: y = 1 + i*8 + x; break;
2032 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2036 gmx_fio_fclose(out);
2041 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2042 gmx_mtop_t *mtop, t_commrec *cr,
2043 int natoms, rvec x[], matrix box)
2045 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2047 int i, ii, resnr, c;
2048 char *atomname, *resname;
2055 natoms = dd->comm->nat[ddnatVSITE];
2058 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2060 sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
2061 sprintf(format4, "%s%s\n", get_pdbformat4(), "%6.2f%6.2f");
2063 out = gmx_fio_fopen(fname, "w");
2065 fprintf(out, "TITLE %s\n", title);
2066 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2067 for (i = 0; i < natoms; i++)
2069 ii = dd->gatindex[i];
2070 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2071 if (i < dd->comm->nat[ddnatZONE])
2074 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2080 else if (i < dd->comm->nat[ddnatVSITE])
2082 b = dd->comm->zones.n;
2086 b = dd->comm->zones.n + 1;
2088 fprintf(out, strlen(atomname) < 4 ? format : format4,
2089 "ATOM", (ii+1)%100000,
2090 atomname, resname, ' ', resnr%10000, ' ',
2091 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2093 fprintf(out, "TER\n");
2095 gmx_fio_fclose(out);
2098 real dd_cutoff_mbody(gmx_domdec_t *dd)
2100 gmx_domdec_comm_t *comm;
2107 if (comm->bInterCGBondeds)
2109 if (comm->cutoff_mbody > 0)
2111 r = comm->cutoff_mbody;
2115 /* cutoff_mbody=0 means we do not have DLB */
2116 r = comm->cellsize_min[dd->dim[0]];
2117 for (di = 1; di < dd->ndim; di++)
2119 r = min(r, comm->cellsize_min[dd->dim[di]]);
2121 if (comm->bBondComm)
2123 r = max(r, comm->cutoff_mbody);
2127 r = min(r, comm->cutoff);
2135 real dd_cutoff_twobody(gmx_domdec_t *dd)
2139 r_mb = dd_cutoff_mbody(dd);
2141 return max(dd->comm->cutoff, r_mb);
2145 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2149 nc = dd->nc[dd->comm->cartpmedim];
2150 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2151 copy_ivec(coord, coord_pme);
2152 coord_pme[dd->comm->cartpmedim] =
2153 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2156 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2158 /* Here we assign a PME node to communicate with this DD node
2159 * by assuming that the major index of both is x.
2160 * We add cr->npmenodes/2 to obtain an even distribution.
2162 return (ddindex*npme + npme/2)/ndd;
2165 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2167 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2170 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2172 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2175 static int *dd_pmenodes(t_commrec *cr)
2180 snew(pmenodes, cr->npmenodes);
2182 for (i = 0; i < cr->dd->nnodes; i++)
2184 p0 = cr_ddindex2pmeindex(cr, i);
2185 p1 = cr_ddindex2pmeindex(cr, i+1);
2186 if (i+1 == cr->dd->nnodes || p1 > p0)
2190 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2192 pmenodes[n] = i + 1 + n;
2200 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2203 ivec coords, coords_pme, nc;
2208 if (dd->comm->bCartesian) {
2209 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2210 dd_coords2pmecoords(dd,coords,coords_pme);
2211 copy_ivec(dd->ntot,nc);
2212 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2213 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2215 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2217 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2223 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2228 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2230 gmx_domdec_comm_t *comm;
2232 int ddindex, nodeid = -1;
2234 comm = cr->dd->comm;
2239 if (comm->bCartesianPP_PME)
2242 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2247 ddindex = dd_index(cr->dd->nc, coords);
2248 if (comm->bCartesianPP)
2250 nodeid = comm->ddindex2simnodeid[ddindex];
2256 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2268 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2271 gmx_domdec_comm_t *comm;
2272 ivec coord, coord_pme;
2279 /* This assumes a uniform x domain decomposition grid cell size */
2280 if (comm->bCartesianPP_PME)
2283 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2284 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2286 /* This is a PP node */
2287 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2288 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2292 else if (comm->bCartesianPP)
2294 if (sim_nodeid < dd->nnodes)
2296 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2301 /* This assumes DD cells with identical x coordinates
2302 * are numbered sequentially.
2304 if (dd->comm->pmenodes == NULL)
2306 if (sim_nodeid < dd->nnodes)
2308 /* The DD index equals the nodeid */
2309 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2315 while (sim_nodeid > dd->comm->pmenodes[i])
2319 if (sim_nodeid < dd->comm->pmenodes[i])
2321 pmenode = dd->comm->pmenodes[i];
2329 void get_pme_nnodes(const gmx_domdec_t *dd,
2330 int *npmenodes_x, int *npmenodes_y)
2334 *npmenodes_x = dd->comm->npmenodes_x;
2335 *npmenodes_y = dd->comm->npmenodes_y;
2344 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2346 gmx_bool bPMEOnlyNode;
2348 if (DOMAINDECOMP(cr))
2350 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2354 bPMEOnlyNode = FALSE;
2357 return bPMEOnlyNode;
2360 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2361 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2365 ivec coord, coord_pme;
2369 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2372 for (x = 0; x < dd->nc[XX]; x++)
2374 for (y = 0; y < dd->nc[YY]; y++)
2376 for (z = 0; z < dd->nc[ZZ]; z++)
2378 if (dd->comm->bCartesianPP_PME)
2383 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2384 if (dd->ci[XX] == coord_pme[XX] &&
2385 dd->ci[YY] == coord_pme[YY] &&
2386 dd->ci[ZZ] == coord_pme[ZZ])
2388 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2393 /* The slab corresponds to the nodeid in the PME group */
2394 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2396 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2403 /* The last PP-only node is the peer node */
2404 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2408 fprintf(debug, "Receive coordinates from PP nodes:");
2409 for (x = 0; x < *nmy_ddnodes; x++)
2411 fprintf(debug, " %d", (*my_ddnodes)[x]);
2413 fprintf(debug, "\n");
2417 static gmx_bool receive_vir_ener(t_commrec *cr)
2419 gmx_domdec_comm_t *comm;
2420 int pmenode, coords[DIM], rank;
2424 if (cr->npmenodes < cr->dd->nnodes)
2426 comm = cr->dd->comm;
2427 if (comm->bCartesianPP_PME)
2429 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2431 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2432 coords[comm->cartpmedim]++;
2433 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2435 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2436 if (dd_simnode2pmenode(cr, rank) == pmenode)
2438 /* This is not the last PP node for pmenode */
2446 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2447 if (cr->sim_nodeid+1 < cr->nnodes &&
2448 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2450 /* This is not the last PP node for pmenode */
2459 static void set_zones_ncg_home(gmx_domdec_t *dd)
2461 gmx_domdec_zones_t *zones;
2464 zones = &dd->comm->zones;
2466 zones->cg_range[0] = 0;
2467 for (i = 1; i < zones->n+1; i++)
2469 zones->cg_range[i] = dd->ncg_home;
2471 /* zone_ncg1[0] should always be equal to ncg_home */
2472 dd->comm->zone_ncg1[0] = dd->ncg_home;
2475 static void rebuild_cgindex(gmx_domdec_t *dd,
2476 const int *gcgs_index, t_state *state)
2478 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2481 dd_cg_gl = dd->index_gl;
2482 cgindex = dd->cgindex;
2485 for (i = 0; i < state->ncg_gl; i++)
2489 dd_cg_gl[i] = cg_gl;
2490 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2494 dd->ncg_home = state->ncg_gl;
2497 set_zones_ncg_home(dd);
2500 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2502 while (cg >= cginfo_mb->cg_end)
2507 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2510 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2511 t_forcerec *fr, char *bLocalCG)
2513 cginfo_mb_t *cginfo_mb;
2519 cginfo_mb = fr->cginfo_mb;
2520 cginfo = fr->cginfo;
2522 for (cg = cg0; cg < cg1; cg++)
2524 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2528 if (bLocalCG != NULL)
2530 for (cg = cg0; cg < cg1; cg++)
2532 bLocalCG[index_gl[cg]] = TRUE;
2537 static void make_dd_indices(gmx_domdec_t *dd,
2538 const int *gcgs_index, int cg_start)
2540 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2541 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2546 bLocalCG = dd->comm->bLocalCG;
2548 if (dd->nat_tot > dd->gatindex_nalloc)
2550 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2551 srenew(dd->gatindex, dd->gatindex_nalloc);
2554 nzone = dd->comm->zones.n;
2555 zone2cg = dd->comm->zones.cg_range;
2556 zone_ncg1 = dd->comm->zone_ncg1;
2557 index_gl = dd->index_gl;
2558 gatindex = dd->gatindex;
2559 bCGs = dd->comm->bCGs;
2561 if (zone2cg[1] != dd->ncg_home)
2563 gmx_incons("dd->ncg_zone is not up to date");
2566 /* Make the local to global and global to local atom index */
2567 a = dd->cgindex[cg_start];
2568 for (zone = 0; zone < nzone; zone++)
2576 cg0 = zone2cg[zone];
2578 cg1 = zone2cg[zone+1];
2579 cg1_p1 = cg0 + zone_ncg1[zone];
2581 for (cg = cg0; cg < cg1; cg++)
2586 /* Signal that this cg is from more than one pulse away */
2589 cg_gl = index_gl[cg];
2592 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2595 ga2la_set(dd->ga2la, a_gl, a, zone1);
2601 gatindex[a] = cg_gl;
2602 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2609 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2612 int ncg, i, ngl, nerr;
2615 if (bLocalCG == NULL)
2619 for (i = 0; i < dd->ncg_tot; i++)
2621 if (!bLocalCG[dd->index_gl[i]])
2624 "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);
2629 for (i = 0; i < ncg_sys; i++)
2636 if (ngl != dd->ncg_tot)
2638 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);
2645 static void check_index_consistency(gmx_domdec_t *dd,
2646 int natoms_sys, int ncg_sys,
2649 int nerr, ngl, i, a, cell;
2654 if (dd->comm->DD_debug > 1)
2656 snew(have, natoms_sys);
2657 for (a = 0; a < dd->nat_tot; a++)
2659 if (have[dd->gatindex[a]] > 0)
2661 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);
2665 have[dd->gatindex[a]] = a + 1;
2671 snew(have, dd->nat_tot);
2674 for (i = 0; i < natoms_sys; i++)
2676 if (ga2la_get(dd->ga2la, i, &a, &cell))
2678 if (a >= dd->nat_tot)
2680 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);
2686 if (dd->gatindex[a] != i)
2688 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);
2695 if (ngl != dd->nat_tot)
2698 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2699 dd->rank, where, ngl, dd->nat_tot);
2701 for (a = 0; a < dd->nat_tot; a++)
2706 "DD node %d, %s: local atom %d, global %d has no global index\n",
2707 dd->rank, where, a+1, dd->gatindex[a]+1);
2712 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2716 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2717 dd->rank, where, nerr);
2721 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2728 /* Clear the whole list without searching */
2729 ga2la_clear(dd->ga2la);
2733 for (i = a_start; i < dd->nat_tot; i++)
2735 ga2la_del(dd->ga2la, dd->gatindex[i]);
2739 bLocalCG = dd->comm->bLocalCG;
2742 for (i = cg_start; i < dd->ncg_tot; i++)
2744 bLocalCG[dd->index_gl[i]] = FALSE;
2748 dd_clear_local_vsite_indices(dd);
2750 if (dd->constraints)
2752 dd_clear_local_constraint_indices(dd);
2756 /* This function should be used for moving the domain boudaries during DLB,
2757 * for obtaining the minimum cell size. It checks the initially set limit
2758 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2759 * and, possibly, a longer cut-off limit set for PME load balancing.
2761 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2765 cellsize_min = comm->cellsize_min[dim];
2767 if (!comm->bVacDLBNoLimit)
2769 /* The cut-off might have changed, e.g. by PME load balacning,
2770 * from the value used to set comm->cellsize_min, so check it.
2772 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2774 if (comm->bPMELoadBalDLBLimits)
2776 /* Check for the cut-off limit set by the PME load balancing */
2777 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2781 return cellsize_min;
2784 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2787 real grid_jump_limit;
2789 /* The distance between the boundaries of cells at distance
2790 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2791 * and by the fact that cells should not be shifted by more than
2792 * half their size, such that cg's only shift by one cell
2793 * at redecomposition.
2795 grid_jump_limit = comm->cellsize_limit;
2796 if (!comm->bVacDLBNoLimit)
2798 if (comm->bPMELoadBalDLBLimits)
2800 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2802 grid_jump_limit = max(grid_jump_limit,
2803 cutoff/comm->cd[dim_ind].np);
2806 return grid_jump_limit;
2809 static gmx_bool check_grid_jump(gmx_large_int_t step,
2815 gmx_domdec_comm_t *comm;
2824 for (d = 1; d < dd->ndim; d++)
2827 limit = grid_jump_limit(comm, cutoff, d);
2828 bfac = ddbox->box_size[dim];
2829 if (ddbox->tric_dir[dim])
2831 bfac *= ddbox->skew_fac[dim];
2833 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2834 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2842 /* This error should never be triggered under normal
2843 * circumstances, but you never know ...
2845 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.",
2846 gmx_step_str(step, buf),
2847 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2855 static int dd_load_count(gmx_domdec_comm_t *comm)
2857 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2860 static float dd_force_load(gmx_domdec_comm_t *comm)
2867 if (comm->eFlop > 1)
2869 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2874 load = comm->cycl[ddCyclF];
2875 if (comm->cycl_n[ddCyclF] > 1)
2877 /* Subtract the maximum of the last n cycle counts
2878 * to get rid of possible high counts due to other soures,
2879 * for instance system activity, that would otherwise
2880 * affect the dynamic load balancing.
2882 load -= comm->cycl_max[ddCyclF];
2889 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2891 gmx_domdec_comm_t *comm;
2896 snew(*dim_f, dd->nc[dim]+1);
2898 for (i = 1; i < dd->nc[dim]; i++)
2900 if (comm->slb_frac[dim])
2902 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2906 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2909 (*dim_f)[dd->nc[dim]] = 1;
2912 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2914 int pmeindex, slab, nso, i;
2917 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2923 ddpme->dim = dimind;
2925 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2927 ddpme->nslab = (ddpme->dim == 0 ?
2928 dd->comm->npmenodes_x :
2929 dd->comm->npmenodes_y);
2931 if (ddpme->nslab <= 1)
2936 nso = dd->comm->npmenodes/ddpme->nslab;
2937 /* Determine for each PME slab the PP location range for dimension dim */
2938 snew(ddpme->pp_min, ddpme->nslab);
2939 snew(ddpme->pp_max, ddpme->nslab);
2940 for (slab = 0; slab < ddpme->nslab; slab++)
2942 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2943 ddpme->pp_max[slab] = 0;
2945 for (i = 0; i < dd->nnodes; i++)
2947 ddindex2xyz(dd->nc, i, xyz);
2948 /* For y only use our y/z slab.
2949 * This assumes that the PME x grid size matches the DD grid size.
2951 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2953 pmeindex = ddindex2pmeindex(dd, i);
2956 slab = pmeindex/nso;
2960 slab = pmeindex % ddpme->nslab;
2962 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2963 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2967 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2970 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2972 if (dd->comm->ddpme[0].dim == XX)
2974 return dd->comm->ddpme[0].maxshift;
2982 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2984 if (dd->comm->ddpme[0].dim == YY)
2986 return dd->comm->ddpme[0].maxshift;
2988 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2990 return dd->comm->ddpme[1].maxshift;
2998 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2999 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3001 gmx_domdec_comm_t *comm;
3004 real range, pme_boundary;
3008 nc = dd->nc[ddpme->dim];
3011 if (!ddpme->dim_match)
3013 /* PP decomposition is not along dim: the worst situation */
3016 else if (ns <= 3 || (bUniform && ns == nc))
3018 /* The optimal situation */
3023 /* We need to check for all pme nodes which nodes they
3024 * could possibly need to communicate with.
3026 xmin = ddpme->pp_min;
3027 xmax = ddpme->pp_max;
3028 /* Allow for atoms to be maximally 2/3 times the cut-off
3029 * out of their DD cell. This is a reasonable balance between
3030 * between performance and support for most charge-group/cut-off
3033 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3034 /* Avoid extra communication when we are exactly at a boundary */
3038 for (s = 0; s < ns; s++)
3040 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3041 pme_boundary = (real)s/ns;
3044 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3046 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3050 pme_boundary = (real)(s+1)/ns;
3053 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3055 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3062 ddpme->maxshift = sh;
3066 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3067 ddpme->dim, ddpme->maxshift);
3071 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3075 for (d = 0; d < dd->ndim; d++)
3078 if (dim < ddbox->nboundeddim &&
3079 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3080 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3082 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",
3083 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3084 dd->nc[dim], dd->comm->cellsize_limit);
3089 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3090 gmx_bool bMaster, ivec npulse)
3092 gmx_domdec_comm_t *comm;
3095 real *cell_x, cell_dx, cellsize;
3099 for (d = 0; d < DIM; d++)
3101 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3103 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3106 cell_dx = ddbox->box_size[d]/dd->nc[d];
3109 for (j = 0; j < dd->nc[d]+1; j++)
3111 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3116 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3117 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3119 cellsize = cell_dx*ddbox->skew_fac[d];
3120 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3124 cellsize_min[d] = cellsize;
3128 /* Statically load balanced grid */
3129 /* Also when we are not doing a master distribution we determine
3130 * all cell borders in a loop to obtain identical values
3131 * to the master distribution case and to determine npulse.
3135 cell_x = dd->ma->cell_x[d];
3139 snew(cell_x, dd->nc[d]+1);
3141 cell_x[0] = ddbox->box0[d];
3142 for (j = 0; j < dd->nc[d]; j++)
3144 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3145 cell_x[j+1] = cell_x[j] + cell_dx;
3146 cellsize = cell_dx*ddbox->skew_fac[d];
3147 while (cellsize*npulse[d] < comm->cutoff &&
3148 npulse[d] < dd->nc[d]-1)
3152 cellsize_min[d] = min(cellsize_min[d], cellsize);
3156 comm->cell_x0[d] = cell_x[dd->ci[d]];
3157 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3161 /* The following limitation is to avoid that a cell would receive
3162 * some of its own home charge groups back over the periodic boundary.
3163 * Double charge groups cause trouble with the global indices.
3165 if (d < ddbox->npbcdim &&
3166 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3168 gmx_fatal_collective(FARGS, NULL, dd,
3169 "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",
3170 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3172 dd->nc[d], dd->nc[d],
3173 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3177 if (!comm->bDynLoadBal)
3179 copy_rvec(cellsize_min, comm->cellsize_min);
3182 for (d = 0; d < comm->npmedecompdim; d++)
3184 set_pme_maxshift(dd, &comm->ddpme[d],
3185 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3186 comm->ddpme[d].slb_dim_f);
3191 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3192 int d, int dim, gmx_domdec_root_t *root,
3194 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3196 gmx_domdec_comm_t *comm;
3197 int ncd, i, j, nmin, nmin_old;
3198 gmx_bool bLimLo, bLimHi;
3200 real fac, halfway, cellsize_limit_f_i, region_size;
3201 gmx_bool bPBC, bLastHi = FALSE;
3202 int nrange[] = {range[0], range[1]};
3204 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3210 bPBC = (dim < ddbox->npbcdim);
3212 cell_size = root->buf_ncd;
3216 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3219 /* First we need to check if the scaling does not make cells
3220 * smaller than the smallest allowed size.
3221 * We need to do this iteratively, since if a cell is too small,
3222 * it needs to be enlarged, which makes all the other cells smaller,
3223 * which could in turn make another cell smaller than allowed.
3225 for (i = range[0]; i < range[1]; i++)
3227 root->bCellMin[i] = FALSE;
3233 /* We need the total for normalization */
3235 for (i = range[0]; i < range[1]; i++)
3237 if (root->bCellMin[i] == FALSE)
3239 fac += cell_size[i];
3242 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3243 /* Determine the cell boundaries */
3244 for (i = range[0]; i < range[1]; i++)
3246 if (root->bCellMin[i] == FALSE)
3248 cell_size[i] *= fac;
3249 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3251 cellsize_limit_f_i = 0;
3255 cellsize_limit_f_i = cellsize_limit_f;
3257 if (cell_size[i] < cellsize_limit_f_i)
3259 root->bCellMin[i] = TRUE;
3260 cell_size[i] = cellsize_limit_f_i;
3264 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3267 while (nmin > nmin_old);
3270 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3271 /* For this check we should not use DD_CELL_MARGIN,
3272 * but a slightly smaller factor,
3273 * since rounding could get use below the limit.
3275 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3278 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",
3279 gmx_step_str(step, buf),
3280 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3281 ncd, comm->cellsize_min[dim]);
3284 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3288 /* Check if the boundary did not displace more than halfway
3289 * each of the cells it bounds, as this could cause problems,
3290 * especially when the differences between cell sizes are large.
3291 * If changes are applied, they will not make cells smaller
3292 * than the cut-off, as we check all the boundaries which
3293 * might be affected by a change and if the old state was ok,
3294 * the cells will at most be shrunk back to their old size.
3296 for (i = range[0]+1; i < range[1]; i++)
3298 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3299 if (root->cell_f[i] < halfway)
3301 root->cell_f[i] = halfway;
3302 /* Check if the change also causes shifts of the next boundaries */
3303 for (j = i+1; j < range[1]; j++)
3305 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3307 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3311 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3312 if (root->cell_f[i] > halfway)
3314 root->cell_f[i] = halfway;
3315 /* Check if the change also causes shifts of the next boundaries */
3316 for (j = i-1; j >= range[0]+1; j--)
3318 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3320 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3327 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3328 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3329 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3330 * for a and b nrange is used */
3333 /* Take care of the staggering of the cell boundaries */
3336 for (i = range[0]; i < range[1]; i++)
3338 root->cell_f_max0[i] = root->cell_f[i];
3339 root->cell_f_min1[i] = root->cell_f[i+1];
3344 for (i = range[0]+1; i < range[1]; i++)
3346 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3347 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3348 if (bLimLo && bLimHi)
3350 /* Both limits violated, try the best we can */
3351 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3352 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3353 nrange[0] = range[0];
3355 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3358 nrange[1] = range[1];
3359 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3365 /* root->cell_f[i] = root->bound_min[i]; */
3366 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3369 else if (bLimHi && !bLastHi)
3372 if (nrange[1] < range[1]) /* found a LimLo before */
3374 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3375 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3376 nrange[0] = nrange[1];
3378 root->cell_f[i] = root->bound_max[i];
3380 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3382 nrange[1] = range[1];
3385 if (nrange[1] < range[1]) /* found last a LimLo */
3387 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3388 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3389 nrange[0] = nrange[1];
3390 nrange[1] = range[1];
3391 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3393 else if (nrange[0] > range[0]) /* found at least one LimHi */
3395 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3402 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3403 int d, int dim, gmx_domdec_root_t *root,
3404 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3405 gmx_bool bUniform, gmx_large_int_t step)
3407 gmx_domdec_comm_t *comm;
3408 int ncd, d1, i, j, pos;
3410 real load_aver, load_i, imbalance, change, change_max, sc;
3411 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3415 int range[] = { 0, 0 };
3419 /* Convert the maximum change from the input percentage to a fraction */
3420 change_limit = comm->dlb_scale_lim*0.01;
3424 bPBC = (dim < ddbox->npbcdim);
3426 cell_size = root->buf_ncd;
3428 /* Store the original boundaries */
3429 for (i = 0; i < ncd+1; i++)
3431 root->old_cell_f[i] = root->cell_f[i];
3435 for (i = 0; i < ncd; i++)
3437 cell_size[i] = 1.0/ncd;
3440 else if (dd_load_count(comm))
3442 load_aver = comm->load[d].sum_m/ncd;
3444 for (i = 0; i < ncd; i++)
3446 /* Determine the relative imbalance of cell i */
3447 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3448 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3449 /* Determine the change of the cell size using underrelaxation */
3450 change = -relax*imbalance;
3451 change_max = max(change_max, max(change, -change));
3453 /* Limit the amount of scaling.
3454 * We need to use the same rescaling for all cells in one row,
3455 * otherwise the load balancing might not converge.
3458 if (change_max > change_limit)
3460 sc *= change_limit/change_max;
3462 for (i = 0; i < ncd; i++)
3464 /* Determine the relative imbalance of cell i */
3465 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3466 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3467 /* Determine the change of the cell size using underrelaxation */
3468 change = -sc*imbalance;
3469 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3473 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3474 cellsize_limit_f *= DD_CELL_MARGIN;
3475 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3476 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3477 if (ddbox->tric_dir[dim])
3479 cellsize_limit_f /= ddbox->skew_fac[dim];
3480 dist_min_f /= ddbox->skew_fac[dim];
3482 if (bDynamicBox && d > 0)
3484 dist_min_f *= DD_PRES_SCALE_MARGIN;
3486 if (d > 0 && !bUniform)
3488 /* Make sure that the grid is not shifted too much */
3489 for (i = 1; i < ncd; i++)
3491 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3493 gmx_incons("Inconsistent DD boundary staggering limits!");
3495 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3496 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3499 root->bound_min[i] += 0.5*space;
3501 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3502 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3505 root->bound_max[i] += 0.5*space;
3510 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3512 root->cell_f_max0[i-1] + dist_min_f,
3513 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3514 root->cell_f_min1[i] - dist_min_f);
3519 root->cell_f[0] = 0;
3520 root->cell_f[ncd] = 1;
3521 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3524 /* After the checks above, the cells should obey the cut-off
3525 * restrictions, but it does not hurt to check.
3527 for (i = 0; i < ncd; i++)
3531 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3532 dim, i, root->cell_f[i], root->cell_f[i+1]);
3535 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3536 root->cell_f[i+1] - root->cell_f[i] <
3537 cellsize_limit_f/DD_CELL_MARGIN)
3541 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3542 gmx_step_str(step, buf), dim2char(dim), i,
3543 (root->cell_f[i+1] - root->cell_f[i])
3544 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3549 /* Store the cell boundaries of the lower dimensions at the end */
3550 for (d1 = 0; d1 < d; d1++)
3552 root->cell_f[pos++] = comm->cell_f0[d1];
3553 root->cell_f[pos++] = comm->cell_f1[d1];
3556 if (d < comm->npmedecompdim)
3558 /* The master determines the maximum shift for
3559 * the coordinate communication between separate PME nodes.
3561 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3563 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3566 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3570 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3571 gmx_ddbox_t *ddbox, int dimind)
3573 gmx_domdec_comm_t *comm;
3578 /* Set the cell dimensions */
3579 dim = dd->dim[dimind];
3580 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3581 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3582 if (dim >= ddbox->nboundeddim)
3584 comm->cell_x0[dim] += ddbox->box0[dim];
3585 comm->cell_x1[dim] += ddbox->box0[dim];
3589 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3590 int d, int dim, real *cell_f_row,
3593 gmx_domdec_comm_t *comm;
3599 /* Each node would only need to know two fractions,
3600 * but it is probably cheaper to broadcast the whole array.
3602 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3603 0, comm->mpi_comm_load[d]);
3605 /* Copy the fractions for this dimension from the buffer */
3606 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3607 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3608 /* The whole array was communicated, so set the buffer position */
3609 pos = dd->nc[dim] + 1;
3610 for (d1 = 0; d1 <= d; d1++)
3614 /* Copy the cell fractions of the lower dimensions */
3615 comm->cell_f0[d1] = cell_f_row[pos++];
3616 comm->cell_f1[d1] = cell_f_row[pos++];
3618 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3620 /* Convert the communicated shift from float to int */
3621 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3624 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3628 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3629 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3630 gmx_bool bUniform, gmx_large_int_t step)
3632 gmx_domdec_comm_t *comm;
3634 gmx_bool bRowMember, bRowRoot;
3639 for (d = 0; d < dd->ndim; d++)
3644 for (d1 = d; d1 < dd->ndim; d1++)
3646 if (dd->ci[dd->dim[d1]] > 0)
3659 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3660 ddbox, bDynamicBox, bUniform, step);
3661 cell_f_row = comm->root[d]->cell_f;
3665 cell_f_row = comm->cell_f_row;
3667 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3672 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3676 /* This function assumes the box is static and should therefore
3677 * not be called when the box has changed since the last
3678 * call to dd_partition_system.
3680 for (d = 0; d < dd->ndim; d++)
3682 relative_to_absolute_cell_bounds(dd, ddbox, d);
3688 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3689 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3690 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3691 gmx_wallcycle_t wcycle)
3693 gmx_domdec_comm_t *comm;
3700 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3701 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3702 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3704 else if (bDynamicBox)
3706 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3709 /* Set the dimensions for which no DD is used */
3710 for (dim = 0; dim < DIM; dim++)
3712 if (dd->nc[dim] == 1)
3714 comm->cell_x0[dim] = 0;
3715 comm->cell_x1[dim] = ddbox->box_size[dim];
3716 if (dim >= ddbox->nboundeddim)
3718 comm->cell_x0[dim] += ddbox->box0[dim];
3719 comm->cell_x1[dim] += ddbox->box0[dim];
3725 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3728 gmx_domdec_comm_dim_t *cd;
3730 for (d = 0; d < dd->ndim; d++)
3732 cd = &dd->comm->cd[d];
3733 np = npulse[dd->dim[d]];
3734 if (np > cd->np_nalloc)
3738 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3739 dim2char(dd->dim[d]), np);
3741 if (DDMASTER(dd) && cd->np_nalloc > 0)
3743 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3745 srenew(cd->ind, np);
3746 for (i = cd->np_nalloc; i < np; i++)
3748 cd->ind[i].index = NULL;
3749 cd->ind[i].nalloc = 0;
3758 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3759 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3760 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3761 gmx_wallcycle_t wcycle)
3763 gmx_domdec_comm_t *comm;
3769 /* Copy the old cell boundaries for the cg displacement check */
3770 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3771 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3773 if (comm->bDynLoadBal)
3777 check_box_size(dd, ddbox);
3779 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3783 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3784 realloc_comm_ind(dd, npulse);
3789 for (d = 0; d < DIM; d++)
3791 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3792 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3797 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3799 rvec cell_ns_x0, rvec cell_ns_x1,
3800 gmx_large_int_t step)
3802 gmx_domdec_comm_t *comm;
3807 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3809 dim = dd->dim[dim_ind];
3811 /* Without PBC we don't have restrictions on the outer cells */
3812 if (!(dim >= ddbox->npbcdim &&
3813 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3814 comm->bDynLoadBal &&
3815 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3816 comm->cellsize_min[dim])
3819 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",
3820 gmx_step_str(step, buf), dim2char(dim),
3821 comm->cell_x1[dim] - comm->cell_x0[dim],
3822 ddbox->skew_fac[dim],
3823 dd->comm->cellsize_min[dim],
3824 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3828 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3830 /* Communicate the boundaries and update cell_ns_x0/1 */
3831 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3832 if (dd->bGridJump && dd->ndim > 1)
3834 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3839 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3843 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3851 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3852 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3861 static void check_screw_box(matrix box)
3863 /* Mathematical limitation */
3864 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3866 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3869 /* Limitation due to the asymmetry of the eighth shell method */
3870 if (box[ZZ][YY] != 0)
3872 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3876 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3877 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3880 gmx_domdec_master_t *ma;
3881 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3882 int i, icg, j, k, k0, k1, d, npbcdim;
3884 rvec box_size, cg_cm;
3886 real nrcg, inv_ncg, pos_d;
3888 gmx_bool bUnbounded, bScrew;
3892 if (tmp_ind == NULL)
3894 snew(tmp_nalloc, dd->nnodes);
3895 snew(tmp_ind, dd->nnodes);
3896 for (i = 0; i < dd->nnodes; i++)
3898 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3899 snew(tmp_ind[i], tmp_nalloc[i]);
3903 /* Clear the count */
3904 for (i = 0; i < dd->nnodes; i++)
3910 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3912 cgindex = cgs->index;
3914 /* Compute the center of geometry for all charge groups */
3915 for (icg = 0; icg < cgs->nr; icg++)
3918 k1 = cgindex[icg+1];
3922 copy_rvec(pos[k0], cg_cm);
3929 for (k = k0; (k < k1); k++)
3931 rvec_inc(cg_cm, pos[k]);
3933 for (d = 0; (d < DIM); d++)
3935 cg_cm[d] *= inv_ncg;
3938 /* Put the charge group in the box and determine the cell index */
3939 for (d = DIM-1; d >= 0; d--)
3942 if (d < dd->npbcdim)
3944 bScrew = (dd->bScrewPBC && d == XX);
3945 if (tric_dir[d] && dd->nc[d] > 1)
3947 /* Use triclinic coordintates for this dimension */
3948 for (j = d+1; j < DIM; j++)
3950 pos_d += cg_cm[j]*tcm[j][d];
3953 while (pos_d >= box[d][d])
3956 rvec_dec(cg_cm, box[d]);
3959 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3960 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3962 for (k = k0; (k < k1); k++)
3964 rvec_dec(pos[k], box[d]);
3967 pos[k][YY] = box[YY][YY] - pos[k][YY];
3968 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3975 rvec_inc(cg_cm, box[d]);
3978 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3979 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3981 for (k = k0; (k < k1); k++)
3983 rvec_inc(pos[k], box[d]);
3986 pos[k][YY] = box[YY][YY] - pos[k][YY];
3987 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3992 /* This could be done more efficiently */
3994 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3999 i = dd_index(dd->nc, ind);
4000 if (ma->ncg[i] == tmp_nalloc[i])
4002 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4003 srenew(tmp_ind[i], tmp_nalloc[i]);
4005 tmp_ind[i][ma->ncg[i]] = icg;
4007 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4011 for (i = 0; i < dd->nnodes; i++)
4014 for (k = 0; k < ma->ncg[i]; k++)
4016 ma->cg[k1++] = tmp_ind[i][k];
4019 ma->index[dd->nnodes] = k1;
4021 for (i = 0; i < dd->nnodes; i++)
4031 fprintf(fplog, "Charge group distribution at step %s:",
4032 gmx_step_str(step, buf));
4033 for (i = 0; i < dd->nnodes; i++)
4035 fprintf(fplog, " %d", ma->ncg[i]);
4037 fprintf(fplog, "\n");
4041 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4042 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4045 gmx_domdec_master_t *ma = NULL;
4048 int *ibuf, buf2[2] = { 0, 0 };
4049 gmx_bool bMaster = DDMASTER(dd);
4056 check_screw_box(box);
4059 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4061 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4062 for (i = 0; i < dd->nnodes; i++)
4064 ma->ibuf[2*i] = ma->ncg[i];
4065 ma->ibuf[2*i+1] = ma->nat[i];
4073 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4075 dd->ncg_home = buf2[0];
4076 dd->nat_home = buf2[1];
4077 dd->ncg_tot = dd->ncg_home;
4078 dd->nat_tot = dd->nat_home;
4079 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4081 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4082 srenew(dd->index_gl, dd->cg_nalloc);
4083 srenew(dd->cgindex, dd->cg_nalloc+1);
4087 for (i = 0; i < dd->nnodes; i++)
4089 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4090 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4095 DDMASTER(dd) ? ma->ibuf : NULL,
4096 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4097 DDMASTER(dd) ? ma->cg : NULL,
4098 dd->ncg_home*sizeof(int), dd->index_gl);
4100 /* Determine the home charge group sizes */
4102 for (i = 0; i < dd->ncg_home; i++)
4104 cg_gl = dd->index_gl[i];
4106 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4111 fprintf(debug, "Home charge groups:\n");
4112 for (i = 0; i < dd->ncg_home; i++)
4114 fprintf(debug, " %d", dd->index_gl[i]);
4117 fprintf(debug, "\n");
4120 fprintf(debug, "\n");
4124 static int compact_and_copy_vec_at(int ncg, int *move,
4127 rvec *src, gmx_domdec_comm_t *comm,
4130 int m, icg, i, i0, i1, nrcg;
4136 for (m = 0; m < DIM*2; m++)
4142 for (icg = 0; icg < ncg; icg++)
4144 i1 = cgindex[icg+1];
4150 /* Compact the home array in place */
4151 for (i = i0; i < i1; i++)
4153 copy_rvec(src[i], src[home_pos++]);
4159 /* Copy to the communication buffer */
4161 pos_vec[m] += 1 + vec*nrcg;
4162 for (i = i0; i < i1; i++)
4164 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4166 pos_vec[m] += (nvec - vec - 1)*nrcg;
4170 home_pos += i1 - i0;
4178 static int compact_and_copy_vec_cg(int ncg, int *move,
4180 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4183 int m, icg, i0, i1, nrcg;
4189 for (m = 0; m < DIM*2; m++)
4195 for (icg = 0; icg < ncg; icg++)
4197 i1 = cgindex[icg+1];
4203 /* Compact the home array in place */
4204 copy_rvec(src[icg], src[home_pos++]);
4210 /* Copy to the communication buffer */
4211 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4212 pos_vec[m] += 1 + nrcg*nvec;
4224 static int compact_ind(int ncg, int *move,
4225 int *index_gl, int *cgindex,
4227 gmx_ga2la_t ga2la, char *bLocalCG,
4230 int cg, nat, a0, a1, a, a_gl;
4235 for (cg = 0; cg < ncg; cg++)
4241 /* Compact the home arrays in place.
4242 * Anything that can be done here avoids access to global arrays.
4244 cgindex[home_pos] = nat;
4245 for (a = a0; a < a1; a++)
4248 gatindex[nat] = a_gl;
4249 /* The cell number stays 0, so we don't need to set it */
4250 ga2la_change_la(ga2la, a_gl, nat);
4253 index_gl[home_pos] = index_gl[cg];
4254 cginfo[home_pos] = cginfo[cg];
4255 /* The charge group remains local, so bLocalCG does not change */
4260 /* Clear the global indices */
4261 for (a = a0; a < a1; a++)
4263 ga2la_del(ga2la, gatindex[a]);
4267 bLocalCG[index_gl[cg]] = FALSE;
4271 cgindex[home_pos] = nat;
4276 static void clear_and_mark_ind(int ncg, int *move,
4277 int *index_gl, int *cgindex, int *gatindex,
4278 gmx_ga2la_t ga2la, char *bLocalCG,
4283 for (cg = 0; cg < ncg; cg++)
4289 /* Clear the global indices */
4290 for (a = a0; a < a1; a++)
4292 ga2la_del(ga2la, gatindex[a]);
4296 bLocalCG[index_gl[cg]] = FALSE;
4298 /* Signal that this cg has moved using the ns cell index.
4299 * Here we set it to -1. fill_grid will change it
4300 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4302 cell_index[cg] = -1;
4307 static void print_cg_move(FILE *fplog,
4309 gmx_large_int_t step, int cg, int dim, int dir,
4310 gmx_bool bHaveLimitdAndCMOld, real limitd,
4311 rvec cm_old, rvec cm_new, real pos_d)
4313 gmx_domdec_comm_t *comm;
4318 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4319 if (bHaveLimitdAndCMOld)
4321 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4322 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4326 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4327 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4329 fprintf(fplog, "distance out of cell %f\n",
4330 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4331 if (bHaveLimitdAndCMOld)
4333 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4334 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4336 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4337 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4338 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4340 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4341 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4343 comm->cell_x0[dim], comm->cell_x1[dim]);
4346 static void cg_move_error(FILE *fplog,
4348 gmx_large_int_t step, int cg, int dim, int dir,
4349 gmx_bool bHaveLimitdAndCMOld, real limitd,
4350 rvec cm_old, rvec cm_new, real pos_d)
4354 print_cg_move(fplog, dd, step, cg, dim, dir,
4355 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4357 print_cg_move(stderr, dd, step, cg, dim, dir,
4358 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4360 "A charge group moved too far between two domain decomposition steps\n"
4361 "This usually means that your system is not well equilibrated");
4364 static void rotate_state_atom(t_state *state, int a)
4368 for (est = 0; est < estNR; est++)
4370 if (EST_DISTR(est) && (state->flags & (1<<est)))
4375 /* Rotate the complete state; for a rectangular box only */
4376 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4377 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4380 state->v[a][YY] = -state->v[a][YY];
4381 state->v[a][ZZ] = -state->v[a][ZZ];
4384 state->sd_X[a][YY] = -state->sd_X[a][YY];
4385 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4388 state->cg_p[a][YY] = -state->cg_p[a][YY];
4389 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4391 case estDISRE_INITF:
4392 case estDISRE_RM3TAV:
4393 case estORIRE_INITF:
4395 /* These are distances, so not affected by rotation */
4398 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4404 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4406 if (natoms > comm->moved_nalloc)
4408 /* Contents should be preserved here */
4409 comm->moved_nalloc = over_alloc_dd(natoms);
4410 srenew(comm->moved, comm->moved_nalloc);
4416 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4419 ivec tric_dir, matrix tcm,
4420 rvec cell_x0, rvec cell_x1,
4421 rvec limitd, rvec limit0, rvec limit1,
4423 int cg_start, int cg_end,
4428 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4429 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4433 real inv_ncg, pos_d;
4436 npbcdim = dd->npbcdim;
4438 for (cg = cg_start; cg < cg_end; cg++)
4445 copy_rvec(state->x[k0], cm_new);
4452 for (k = k0; (k < k1); k++)
4454 rvec_inc(cm_new, state->x[k]);
4456 for (d = 0; (d < DIM); d++)
4458 cm_new[d] = inv_ncg*cm_new[d];
4463 /* Do pbc and check DD cell boundary crossings */
4464 for (d = DIM-1; d >= 0; d--)
4468 bScrew = (dd->bScrewPBC && d == XX);
4469 /* Determine the location of this cg in lattice coordinates */
4473 for (d2 = d+1; d2 < DIM; d2++)
4475 pos_d += cm_new[d2]*tcm[d2][d];
4478 /* Put the charge group in the triclinic unit-cell */
4479 if (pos_d >= cell_x1[d])
4481 if (pos_d >= limit1[d])
4483 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4484 cg_cm[cg], cm_new, pos_d);
4487 if (dd->ci[d] == dd->nc[d] - 1)
4489 rvec_dec(cm_new, state->box[d]);
4492 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4493 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4495 for (k = k0; (k < k1); k++)
4497 rvec_dec(state->x[k], state->box[d]);
4500 rotate_state_atom(state, k);
4505 else if (pos_d < cell_x0[d])
4507 if (pos_d < limit0[d])
4509 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4510 cg_cm[cg], cm_new, pos_d);
4515 rvec_inc(cm_new, state->box[d]);
4518 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4519 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4521 for (k = k0; (k < k1); k++)
4523 rvec_inc(state->x[k], state->box[d]);
4526 rotate_state_atom(state, k);
4532 else if (d < npbcdim)
4534 /* Put the charge group in the rectangular unit-cell */
4535 while (cm_new[d] >= state->box[d][d])
4537 rvec_dec(cm_new, state->box[d]);
4538 for (k = k0; (k < k1); k++)
4540 rvec_dec(state->x[k], state->box[d]);
4543 while (cm_new[d] < 0)
4545 rvec_inc(cm_new, state->box[d]);
4546 for (k = k0; (k < k1); k++)
4548 rvec_inc(state->x[k], state->box[d]);
4554 copy_rvec(cm_new, cg_cm[cg]);
4556 /* Determine where this cg should go */
4559 for (d = 0; d < dd->ndim; d++)
4564 flag |= DD_FLAG_FW(d);
4570 else if (dev[dim] == -1)
4572 flag |= DD_FLAG_BW(d);
4575 if (dd->nc[dim] > 2)
4586 /* Temporarily store the flag in move */
4587 move[cg] = mc + flag;
4591 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4592 gmx_domdec_t *dd, ivec tric_dir,
4593 t_state *state, rvec **f,
4602 int ncg[DIM*2], nat[DIM*2];
4603 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4604 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4605 int sbuf[2], rbuf[2];
4606 int home_pos_cg, home_pos_at, buf_pos;
4608 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4611 real inv_ncg, pos_d;
4613 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4615 cginfo_mb_t *cginfo_mb;
4616 gmx_domdec_comm_t *comm;
4618 int nthread, thread;
4622 check_screw_box(state->box);
4626 if (fr->cutoff_scheme == ecutsGROUP)
4631 for (i = 0; i < estNR; i++)
4637 case estX: /* Always present */ break;
4638 case estV: bV = (state->flags & (1<<i)); break;
4639 case estSDX: bSDX = (state->flags & (1<<i)); break;
4640 case estCGP: bCGP = (state->flags & (1<<i)); break;
4643 case estDISRE_INITF:
4644 case estDISRE_RM3TAV:
4645 case estORIRE_INITF:
4647 /* No processing required */
4650 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4655 if (dd->ncg_tot > comm->nalloc_int)
4657 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4658 srenew(comm->buf_int, comm->nalloc_int);
4660 move = comm->buf_int;
4662 /* Clear the count */
4663 for (c = 0; c < dd->ndim*2; c++)
4669 npbcdim = dd->npbcdim;
4671 for (d = 0; (d < DIM); d++)
4673 limitd[d] = dd->comm->cellsize_min[d];
4674 if (d >= npbcdim && dd->ci[d] == 0)
4676 cell_x0[d] = -GMX_FLOAT_MAX;
4680 cell_x0[d] = comm->cell_x0[d];
4682 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4684 cell_x1[d] = GMX_FLOAT_MAX;
4688 cell_x1[d] = comm->cell_x1[d];
4692 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4693 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4697 /* We check after communication if a charge group moved
4698 * more than one cell. Set the pre-comm check limit to float_max.
4700 limit0[d] = -GMX_FLOAT_MAX;
4701 limit1[d] = GMX_FLOAT_MAX;
4705 make_tric_corr_matrix(npbcdim, state->box, tcm);
4707 cgindex = dd->cgindex;
4709 nthread = gmx_omp_nthreads_get(emntDomdec);
4711 /* Compute the center of geometry for all home charge groups
4712 * and put them in the box and determine where they should go.
4714 #pragma omp parallel for num_threads(nthread) schedule(static)
4715 for (thread = 0; thread < nthread; thread++)
4717 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4718 cell_x0, cell_x1, limitd, limit0, limit1,
4720 ( thread *dd->ncg_home)/nthread,
4721 ((thread+1)*dd->ncg_home)/nthread,
4722 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4726 for (cg = 0; cg < dd->ncg_home; cg++)
4731 flag = mc & ~DD_FLAG_NRCG;
4732 mc = mc & DD_FLAG_NRCG;
4735 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4737 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4738 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4740 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4741 /* We store the cg size in the lower 16 bits
4742 * and the place where the charge group should go
4743 * in the next 6 bits. This saves some communication volume.
4745 nrcg = cgindex[cg+1] - cgindex[cg];
4746 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4752 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4753 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4756 for (i = 0; i < dd->ndim*2; i++)
4758 *ncg_moved += ncg[i];
4775 /* Make sure the communication buffers are large enough */
4776 for (mc = 0; mc < dd->ndim*2; mc++)
4778 nvr = ncg[mc] + nat[mc]*nvec;
4779 if (nvr > comm->cgcm_state_nalloc[mc])
4781 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4782 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4786 switch (fr->cutoff_scheme)
4789 /* Recalculating cg_cm might be cheaper than communicating,
4790 * but that could give rise to rounding issues.
4793 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4794 nvec, cg_cm, comm, bCompact);
4797 /* Without charge groups we send the moved atom coordinates
4798 * over twice. This is so the code below can be used without
4799 * many conditionals for both for with and without charge groups.
4802 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4803 nvec, state->x, comm, FALSE);
4806 home_pos_cg -= *ncg_moved;
4810 gmx_incons("unimplemented");
4816 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4817 nvec, vec++, state->x, comm, bCompact);
4820 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4821 nvec, vec++, state->v, comm, bCompact);
4825 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4826 nvec, vec++, state->sd_X, comm, bCompact);
4830 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4831 nvec, vec++, state->cg_p, comm, bCompact);
4836 compact_ind(dd->ncg_home, move,
4837 dd->index_gl, dd->cgindex, dd->gatindex,
4838 dd->ga2la, comm->bLocalCG,
4843 if (fr->cutoff_scheme == ecutsVERLET)
4845 moved = get_moved(comm, dd->ncg_home);
4847 for (k = 0; k < dd->ncg_home; k++)
4854 moved = fr->ns.grid->cell_index;
4857 clear_and_mark_ind(dd->ncg_home, move,
4858 dd->index_gl, dd->cgindex, dd->gatindex,
4859 dd->ga2la, comm->bLocalCG,
4863 cginfo_mb = fr->cginfo_mb;
4865 *ncg_stay_home = home_pos_cg;
4866 for (d = 0; d < dd->ndim; d++)
4872 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4875 /* Communicate the cg and atom counts */
4880 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4881 d, dir, sbuf[0], sbuf[1]);
4883 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4885 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4887 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4888 srenew(comm->buf_int, comm->nalloc_int);
4891 /* Communicate the charge group indices, sizes and flags */
4892 dd_sendrecv_int(dd, d, dir,
4893 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4894 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4896 nvs = ncg[cdd] + nat[cdd]*nvec;
4897 i = rbuf[0] + rbuf[1] *nvec;
4898 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4900 /* Communicate cgcm and state */
4901 dd_sendrecv_rvec(dd, d, dir,
4902 comm->cgcm_state[cdd], nvs,
4903 comm->vbuf.v+nvr, i);
4904 ncg_recv += rbuf[0];
4905 nat_recv += rbuf[1];
4909 /* Process the received charge groups */
4911 for (cg = 0; cg < ncg_recv; cg++)
4913 flag = comm->buf_int[cg*DD_CGIBS+1];
4915 if (dim >= npbcdim && dd->nc[dim] > 2)
4917 /* No pbc in this dim and more than one domain boundary.
4918 * We do a separate check if a charge group didn't move too far.
4920 if (((flag & DD_FLAG_FW(d)) &&
4921 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4922 ((flag & DD_FLAG_BW(d)) &&
4923 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4925 cg_move_error(fplog, dd, step, cg, dim,
4926 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4928 comm->vbuf.v[buf_pos],
4929 comm->vbuf.v[buf_pos],
4930 comm->vbuf.v[buf_pos][dim]);
4937 /* Check which direction this cg should go */
4938 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4942 /* The cell boundaries for dimension d2 are not equal
4943 * for each cell row of the lower dimension(s),
4944 * therefore we might need to redetermine where
4945 * this cg should go.
4948 /* If this cg crosses the box boundary in dimension d2
4949 * we can use the communicated flag, so we do not
4950 * have to worry about pbc.
4952 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4953 (flag & DD_FLAG_FW(d2))) ||
4954 (dd->ci[dim2] == 0 &&
4955 (flag & DD_FLAG_BW(d2)))))
4957 /* Clear the two flags for this dimension */
4958 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4959 /* Determine the location of this cg
4960 * in lattice coordinates
4962 pos_d = comm->vbuf.v[buf_pos][dim2];
4965 for (d3 = dim2+1; d3 < DIM; d3++)
4968 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4971 /* Check of we are not at the box edge.
4972 * pbc is only handled in the first step above,
4973 * but this check could move over pbc while
4974 * the first step did not due to different rounding.
4976 if (pos_d >= cell_x1[dim2] &&
4977 dd->ci[dim2] != dd->nc[dim2]-1)
4979 flag |= DD_FLAG_FW(d2);
4981 else if (pos_d < cell_x0[dim2] &&
4984 flag |= DD_FLAG_BW(d2);
4986 comm->buf_int[cg*DD_CGIBS+1] = flag;
4989 /* Set to which neighboring cell this cg should go */
4990 if (flag & DD_FLAG_FW(d2))
4994 else if (flag & DD_FLAG_BW(d2))
4996 if (dd->nc[dd->dim[d2]] > 2)
5008 nrcg = flag & DD_FLAG_NRCG;
5011 if (home_pos_cg+1 > dd->cg_nalloc)
5013 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5014 srenew(dd->index_gl, dd->cg_nalloc);
5015 srenew(dd->cgindex, dd->cg_nalloc+1);
5017 /* Set the global charge group index and size */
5018 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5019 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5020 /* Copy the state from the buffer */
5021 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5022 if (fr->cutoff_scheme == ecutsGROUP)
5025 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5029 /* Set the cginfo */
5030 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5031 dd->index_gl[home_pos_cg]);
5034 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5037 if (home_pos_at+nrcg > state->nalloc)
5039 dd_realloc_state(state, f, home_pos_at+nrcg);
5041 for (i = 0; i < nrcg; i++)
5043 copy_rvec(comm->vbuf.v[buf_pos++],
5044 state->x[home_pos_at+i]);
5048 for (i = 0; i < nrcg; i++)
5050 copy_rvec(comm->vbuf.v[buf_pos++],
5051 state->v[home_pos_at+i]);
5056 for (i = 0; i < nrcg; i++)
5058 copy_rvec(comm->vbuf.v[buf_pos++],
5059 state->sd_X[home_pos_at+i]);
5064 for (i = 0; i < nrcg; i++)
5066 copy_rvec(comm->vbuf.v[buf_pos++],
5067 state->cg_p[home_pos_at+i]);
5071 home_pos_at += nrcg;
5075 /* Reallocate the buffers if necessary */
5076 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5078 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5079 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5081 nvr = ncg[mc] + nat[mc]*nvec;
5082 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5084 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5085 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5087 /* Copy from the receive to the send buffers */
5088 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5089 comm->buf_int + cg*DD_CGIBS,
5090 DD_CGIBS*sizeof(int));
5091 memcpy(comm->cgcm_state[mc][nvr],
5092 comm->vbuf.v[buf_pos],
5093 (1+nrcg*nvec)*sizeof(rvec));
5094 buf_pos += 1 + nrcg*nvec;
5101 /* With sorting (!bCompact) the indices are now only partially up to date
5102 * and ncg_home and nat_home are not the real count, since there are
5103 * "holes" in the arrays for the charge groups that moved to neighbors.
5105 if (fr->cutoff_scheme == ecutsVERLET)
5107 moved = get_moved(comm, home_pos_cg);
5109 for (i = dd->ncg_home; i < home_pos_cg; i++)
5114 dd->ncg_home = home_pos_cg;
5115 dd->nat_home = home_pos_at;
5120 "Finished repartitioning: cgs moved out %d, new home %d\n",
5121 *ncg_moved, dd->ncg_home-*ncg_moved);
5126 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5128 dd->comm->cycl[ddCycl] += cycles;
5129 dd->comm->cycl_n[ddCycl]++;
5130 if (cycles > dd->comm->cycl_max[ddCycl])
5132 dd->comm->cycl_max[ddCycl] = cycles;
5136 static double force_flop_count(t_nrnb *nrnb)
5143 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5145 /* To get closer to the real timings, we half the count
5146 * for the normal loops and again half it for water loops.
5149 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5151 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5155 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5158 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5161 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5163 sum += nrnb->n[i]*cost_nrnb(i);
5166 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5168 sum += nrnb->n[i]*cost_nrnb(i);
5174 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5176 if (dd->comm->eFlop)
5178 dd->comm->flop -= force_flop_count(nrnb);
5181 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5183 if (dd->comm->eFlop)
5185 dd->comm->flop += force_flop_count(nrnb);
5190 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5194 for (i = 0; i < ddCyclNr; i++)
5196 dd->comm->cycl[i] = 0;
5197 dd->comm->cycl_n[i] = 0;
5198 dd->comm->cycl_max[i] = 0;
5201 dd->comm->flop_n = 0;
5204 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5206 gmx_domdec_comm_t *comm;
5207 gmx_domdec_load_t *load;
5208 gmx_domdec_root_t *root = NULL;
5209 int d, dim, cid, i, pos;
5210 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5215 fprintf(debug, "get_load_distribution start\n");
5218 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5222 bSepPME = (dd->pme_nodeid >= 0);
5224 for (d = dd->ndim-1; d >= 0; d--)
5227 /* Check if we participate in the communication in this dimension */
5228 if (d == dd->ndim-1 ||
5229 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5231 load = &comm->load[d];
5234 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5237 if (d == dd->ndim-1)
5239 sbuf[pos++] = dd_force_load(comm);
5240 sbuf[pos++] = sbuf[0];
5243 sbuf[pos++] = sbuf[0];
5244 sbuf[pos++] = cell_frac;
5247 sbuf[pos++] = comm->cell_f_max0[d];
5248 sbuf[pos++] = comm->cell_f_min1[d];
5253 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5254 sbuf[pos++] = comm->cycl[ddCyclPME];
5259 sbuf[pos++] = comm->load[d+1].sum;
5260 sbuf[pos++] = comm->load[d+1].max;
5263 sbuf[pos++] = comm->load[d+1].sum_m;
5264 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5265 sbuf[pos++] = comm->load[d+1].flags;
5268 sbuf[pos++] = comm->cell_f_max0[d];
5269 sbuf[pos++] = comm->cell_f_min1[d];
5274 sbuf[pos++] = comm->load[d+1].mdf;
5275 sbuf[pos++] = comm->load[d+1].pme;
5279 /* Communicate a row in DD direction d.
5280 * The communicators are setup such that the root always has rank 0.
5283 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5284 load->load, load->nload*sizeof(float), MPI_BYTE,
5285 0, comm->mpi_comm_load[d]);
5287 if (dd->ci[dim] == dd->master_ci[dim])
5289 /* We are the root, process this row */
5290 if (comm->bDynLoadBal)
5292 root = comm->root[d];
5302 for (i = 0; i < dd->nc[dim]; i++)
5304 load->sum += load->load[pos++];
5305 load->max = max(load->max, load->load[pos]);
5311 /* This direction could not be load balanced properly,
5312 * therefore we need to use the maximum iso the average load.
5314 load->sum_m = max(load->sum_m, load->load[pos]);
5318 load->sum_m += load->load[pos];
5321 load->cvol_min = min(load->cvol_min, load->load[pos]);
5325 load->flags = (int)(load->load[pos++] + 0.5);
5329 root->cell_f_max0[i] = load->load[pos++];
5330 root->cell_f_min1[i] = load->load[pos++];
5335 load->mdf = max(load->mdf, load->load[pos]);
5337 load->pme = max(load->pme, load->load[pos]);
5341 if (comm->bDynLoadBal && root->bLimited)
5343 load->sum_m *= dd->nc[dim];
5344 load->flags |= (1<<d);
5352 comm->nload += dd_load_count(comm);
5353 comm->load_step += comm->cycl[ddCyclStep];
5354 comm->load_sum += comm->load[0].sum;
5355 comm->load_max += comm->load[0].max;
5356 if (comm->bDynLoadBal)
5358 for (d = 0; d < dd->ndim; d++)
5360 if (comm->load[0].flags & (1<<d))
5362 comm->load_lim[d]++;
5368 comm->load_mdf += comm->load[0].mdf;
5369 comm->load_pme += comm->load[0].pme;
5373 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5377 fprintf(debug, "get_load_distribution finished\n");
5381 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5383 /* Return the relative performance loss on the total run time
5384 * due to the force calculation load imbalance.
5386 if (dd->comm->nload > 0)
5389 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5390 (dd->comm->load_step*dd->nnodes);
5398 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5401 int npp, npme, nnodes, d, limp;
5402 float imbal, pme_f_ratio, lossf, lossp = 0;
5404 gmx_domdec_comm_t *comm;
5407 if (DDMASTER(dd) && comm->nload > 0)
5410 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5411 nnodes = npp + npme;
5412 imbal = comm->load_max*npp/comm->load_sum - 1;
5413 lossf = dd_force_imb_perf_loss(dd);
5414 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5415 fprintf(fplog, "%s", buf);
5416 fprintf(stderr, "\n");
5417 fprintf(stderr, "%s", buf);
5418 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5419 fprintf(fplog, "%s", buf);
5420 fprintf(stderr, "%s", buf);
5422 if (comm->bDynLoadBal)
5424 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5425 for (d = 0; d < dd->ndim; d++)
5427 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5428 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5434 sprintf(buf+strlen(buf), "\n");
5435 fprintf(fplog, "%s", buf);
5436 fprintf(stderr, "%s", buf);
5440 pme_f_ratio = comm->load_pme/comm->load_mdf;
5441 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5444 lossp *= (float)npme/(float)nnodes;
5448 lossp *= (float)npp/(float)nnodes;
5450 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5451 fprintf(fplog, "%s", buf);
5452 fprintf(stderr, "%s", buf);
5453 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5454 fprintf(fplog, "%s", buf);
5455 fprintf(stderr, "%s", buf);
5457 fprintf(fplog, "\n");
5458 fprintf(stderr, "\n");
5460 if (lossf >= DD_PERF_LOSS)
5463 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5464 " in the domain decomposition.\n", lossf*100);
5465 if (!comm->bDynLoadBal)
5467 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5471 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5473 fprintf(fplog, "%s\n", buf);
5474 fprintf(stderr, "%s\n", buf);
5476 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5479 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5480 " had %s work to do than the PP nodes.\n"
5481 " You might want to %s the number of PME nodes\n"
5482 " or %s the cut-off and the grid spacing.\n",
5484 (lossp < 0) ? "less" : "more",
5485 (lossp < 0) ? "decrease" : "increase",
5486 (lossp < 0) ? "decrease" : "increase");
5487 fprintf(fplog, "%s\n", buf);
5488 fprintf(stderr, "%s\n", buf);
5493 static float dd_vol_min(gmx_domdec_t *dd)
5495 return dd->comm->load[0].cvol_min*dd->nnodes;
5498 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5500 return dd->comm->load[0].flags;
5503 static float dd_f_imbal(gmx_domdec_t *dd)
5505 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5508 float dd_pme_f_ratio(gmx_domdec_t *dd)
5510 if (dd->comm->cycl_n[ddCyclPME] > 0)
5512 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5520 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5525 flags = dd_load_flags(dd);
5529 "DD load balancing is limited by minimum cell size in dimension");
5530 for (d = 0; d < dd->ndim; d++)
5534 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5537 fprintf(fplog, "\n");
5539 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5540 if (dd->comm->bDynLoadBal)
5542 fprintf(fplog, " vol min/aver %5.3f%c",
5543 dd_vol_min(dd), flags ? '!' : ' ');
5545 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5546 if (dd->comm->cycl_n[ddCyclPME])
5548 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5550 fprintf(fplog, "\n\n");
5553 static void dd_print_load_verbose(gmx_domdec_t *dd)
5555 if (dd->comm->bDynLoadBal)
5557 fprintf(stderr, "vol %4.2f%c ",
5558 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5560 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5561 if (dd->comm->cycl_n[ddCyclPME])
5563 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5568 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5573 gmx_domdec_root_t *root;
5574 gmx_bool bPartOfGroup = FALSE;
5576 dim = dd->dim[dim_ind];
5577 copy_ivec(loc, loc_c);
5578 for (i = 0; i < dd->nc[dim]; i++)
5581 rank = dd_index(dd->nc, loc_c);
5582 if (rank == dd->rank)
5584 /* This process is part of the group */
5585 bPartOfGroup = TRUE;
5588 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5592 dd->comm->mpi_comm_load[dim_ind] = c_row;
5593 if (dd->comm->eDLB != edlbNO)
5595 if (dd->ci[dim] == dd->master_ci[dim])
5597 /* This is the root process of this row */
5598 snew(dd->comm->root[dim_ind], 1);
5599 root = dd->comm->root[dim_ind];
5600 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5601 snew(root->old_cell_f, dd->nc[dim]+1);
5602 snew(root->bCellMin, dd->nc[dim]);
5605 snew(root->cell_f_max0, dd->nc[dim]);
5606 snew(root->cell_f_min1, dd->nc[dim]);
5607 snew(root->bound_min, dd->nc[dim]);
5608 snew(root->bound_max, dd->nc[dim]);
5610 snew(root->buf_ncd, dd->nc[dim]);
5614 /* This is not a root process, we only need to receive cell_f */
5615 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5618 if (dd->ci[dim] == dd->master_ci[dim])
5620 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5626 static void make_load_communicators(gmx_domdec_t *dd)
5629 int dim0, dim1, i, j;
5634 fprintf(debug, "Making load communicators\n");
5637 snew(dd->comm->load, dd->ndim);
5638 snew(dd->comm->mpi_comm_load, dd->ndim);
5641 make_load_communicator(dd, 0, loc);
5645 for (i = 0; i < dd->nc[dim0]; i++)
5648 make_load_communicator(dd, 1, loc);
5654 for (i = 0; i < dd->nc[dim0]; i++)
5658 for (j = 0; j < dd->nc[dim1]; j++)
5661 make_load_communicator(dd, 2, loc);
5668 fprintf(debug, "Finished making load communicators\n");
5673 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5676 int d, dim, i, j, m;
5679 ivec dd_zp[DD_MAXIZONE];
5680 gmx_domdec_zones_t *zones;
5681 gmx_domdec_ns_ranges_t *izone;
5683 for (d = 0; d < dd->ndim; d++)
5686 copy_ivec(dd->ci, tmp);
5687 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5688 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5689 copy_ivec(dd->ci, tmp);
5690 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5691 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5694 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5697 dd->neighbor[d][1]);
5703 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5705 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5706 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5713 for (i = 0; i < nzonep; i++)
5715 copy_ivec(dd_zp3[i], dd_zp[i]);
5721 for (i = 0; i < nzonep; i++)
5723 copy_ivec(dd_zp2[i], dd_zp[i]);
5729 for (i = 0; i < nzonep; i++)
5731 copy_ivec(dd_zp1[i], dd_zp[i]);
5735 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5740 zones = &dd->comm->zones;
5742 for (i = 0; i < nzone; i++)
5745 clear_ivec(zones->shift[i]);
5746 for (d = 0; d < dd->ndim; d++)
5748 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5753 for (i = 0; i < nzone; i++)
5755 for (d = 0; d < DIM; d++)
5757 s[d] = dd->ci[d] - zones->shift[i][d];
5762 else if (s[d] >= dd->nc[d])
5768 zones->nizone = nzonep;
5769 for (i = 0; i < zones->nizone; i++)
5771 if (dd_zp[i][0] != i)
5773 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5775 izone = &zones->izone[i];
5776 izone->j0 = dd_zp[i][1];
5777 izone->j1 = dd_zp[i][2];
5778 for (dim = 0; dim < DIM; dim++)
5780 if (dd->nc[dim] == 1)
5782 /* All shifts should be allowed */
5783 izone->shift0[dim] = -1;
5784 izone->shift1[dim] = 1;
5789 izone->shift0[d] = 0;
5790 izone->shift1[d] = 0;
5791 for(j=izone->j0; j<izone->j1; j++) {
5792 if (dd->shift[j][d] > dd->shift[i][d])
5793 izone->shift0[d] = -1;
5794 if (dd->shift[j][d] < dd->shift[i][d])
5795 izone->shift1[d] = 1;
5801 /* Assume the shift are not more than 1 cell */
5802 izone->shift0[dim] = 1;
5803 izone->shift1[dim] = -1;
5804 for (j = izone->j0; j < izone->j1; j++)
5806 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5807 if (shift_diff < izone->shift0[dim])
5809 izone->shift0[dim] = shift_diff;
5811 if (shift_diff > izone->shift1[dim])
5813 izone->shift1[dim] = shift_diff;
5820 if (dd->comm->eDLB != edlbNO)
5822 snew(dd->comm->root, dd->ndim);
5825 if (dd->comm->bRecordLoad)
5827 make_load_communicators(dd);
5831 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5834 gmx_domdec_comm_t *comm;
5845 if (comm->bCartesianPP)
5847 /* Set up cartesian communication for the particle-particle part */
5850 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5851 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5854 for (i = 0; i < DIM; i++)
5858 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5860 /* We overwrite the old communicator with the new cartesian one */
5861 cr->mpi_comm_mygroup = comm_cart;
5864 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5865 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5867 if (comm->bCartesianPP_PME)
5869 /* Since we want to use the original cartesian setup for sim,
5870 * and not the one after split, we need to make an index.
5872 snew(comm->ddindex2ddnodeid, dd->nnodes);
5873 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5874 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5875 /* Get the rank of the DD master,
5876 * above we made sure that the master node is a PP node.
5886 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5888 else if (comm->bCartesianPP)
5890 if (cr->npmenodes == 0)
5892 /* The PP communicator is also
5893 * the communicator for this simulation
5895 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5897 cr->nodeid = dd->rank;
5899 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5901 /* We need to make an index to go from the coordinates
5902 * to the nodeid of this simulation.
5904 snew(comm->ddindex2simnodeid, dd->nnodes);
5905 snew(buf, dd->nnodes);
5906 if (cr->duty & DUTY_PP)
5908 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5910 /* Communicate the ddindex to simulation nodeid index */
5911 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5912 cr->mpi_comm_mysim);
5915 /* Determine the master coordinates and rank.
5916 * The DD master should be the same node as the master of this sim.
5918 for (i = 0; i < dd->nnodes; i++)
5920 if (comm->ddindex2simnodeid[i] == 0)
5922 ddindex2xyz(dd->nc, i, dd->master_ci);
5923 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5928 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5933 /* No Cartesian communicators */
5934 /* We use the rank in dd->comm->all as DD index */
5935 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5936 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5938 clear_ivec(dd->master_ci);
5945 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5946 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5951 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5952 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5956 static void receive_ddindex2simnodeid(t_commrec *cr)
5960 gmx_domdec_comm_t *comm;
5967 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5969 snew(comm->ddindex2simnodeid, dd->nnodes);
5970 snew(buf, dd->nnodes);
5971 if (cr->duty & DUTY_PP)
5973 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5976 /* Communicate the ddindex to simulation nodeid index */
5977 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5978 cr->mpi_comm_mysim);
5985 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5986 int ncg, int natoms)
5988 gmx_domdec_master_t *ma;
5993 snew(ma->ncg, dd->nnodes);
5994 snew(ma->index, dd->nnodes+1);
5996 snew(ma->nat, dd->nnodes);
5997 snew(ma->ibuf, dd->nnodes*2);
5998 snew(ma->cell_x, DIM);
5999 for (i = 0; i < DIM; i++)
6001 snew(ma->cell_x[i], dd->nc[i]+1);
6004 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6010 snew(ma->vbuf, natoms);
6016 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
6020 gmx_domdec_comm_t *comm;
6031 if (comm->bCartesianPP)
6033 for (i = 1; i < DIM; i++)
6035 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6037 if (bDiv[YY] || bDiv[ZZ])
6039 comm->bCartesianPP_PME = TRUE;
6040 /* If we have 2D PME decomposition, which is always in x+y,
6041 * we stack the PME only nodes in z.
6042 * Otherwise we choose the direction that provides the thinnest slab
6043 * of PME only nodes as this will have the least effect
6044 * on the PP communication.
6045 * But for the PME communication the opposite might be better.
6047 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6049 dd->nc[YY] > dd->nc[ZZ]))
6051 comm->cartpmedim = ZZ;
6055 comm->cartpmedim = YY;
6057 comm->ntot[comm->cartpmedim]
6058 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6062 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]);
6064 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6069 if (comm->bCartesianPP_PME)
6073 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]);
6076 for (i = 0; i < DIM; i++)
6080 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6083 MPI_Comm_rank(comm_cart, &rank);
6084 if (MASTERNODE(cr) && rank != 0)
6086 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6089 /* With this assigment we loose the link to the original communicator
6090 * which will usually be MPI_COMM_WORLD, unless have multisim.
6092 cr->mpi_comm_mysim = comm_cart;
6093 cr->sim_nodeid = rank;
6095 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6099 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6100 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6103 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6107 if (cr->npmenodes == 0 ||
6108 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6110 cr->duty = DUTY_PME;
6113 /* Split the sim communicator into PP and PME only nodes */
6114 MPI_Comm_split(cr->mpi_comm_mysim,
6116 dd_index(comm->ntot, dd->ci),
6117 &cr->mpi_comm_mygroup);
6121 switch (dd_node_order)
6126 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6129 case ddnoINTERLEAVE:
6130 /* Interleave the PP-only and PME-only nodes,
6131 * as on clusters with dual-core machines this will double
6132 * the communication bandwidth of the PME processes
6133 * and thus speed up the PP <-> PME and inter PME communication.
6137 fprintf(fplog, "Interleaving PP and PME nodes\n");
6139 comm->pmenodes = dd_pmenodes(cr);
6144 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6147 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6149 cr->duty = DUTY_PME;
6156 /* Split the sim communicator into PP and PME only nodes */
6157 MPI_Comm_split(cr->mpi_comm_mysim,
6160 &cr->mpi_comm_mygroup);
6161 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6167 fprintf(fplog, "This is a %s only node\n\n",
6168 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6172 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6175 gmx_domdec_comm_t *comm;
6181 copy_ivec(dd->nc, comm->ntot);
6183 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6184 comm->bCartesianPP_PME = FALSE;
6186 /* Reorder the nodes by default. This might change the MPI ranks.
6187 * Real reordering is only supported on very few architectures,
6188 * Blue Gene is one of them.
6190 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6192 if (cr->npmenodes > 0)
6194 /* Split the communicator into a PP and PME part */
6195 split_communicator(fplog, cr, dd_node_order, CartReorder);
6196 if (comm->bCartesianPP_PME)
6198 /* We (possibly) reordered the nodes in split_communicator,
6199 * so it is no longer required in make_pp_communicator.
6201 CartReorder = FALSE;
6206 /* All nodes do PP and PME */
6208 /* We do not require separate communicators */
6209 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6213 if (cr->duty & DUTY_PP)
6215 /* Copy or make a new PP communicator */
6216 make_pp_communicator(fplog, cr, CartReorder);
6220 receive_ddindex2simnodeid(cr);
6223 if (!(cr->duty & DUTY_PME))
6225 /* Set up the commnuication to our PME node */
6226 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6227 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6230 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6231 dd->pme_nodeid, dd->pme_receive_vir_ener);
6236 dd->pme_nodeid = -1;
6241 dd->ma = init_gmx_domdec_master_t(dd,
6243 comm->cgs_gl.index[comm->cgs_gl.nr]);
6247 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6249 real *slb_frac, tot;
6254 if (nc > 1 && size_string != NULL)
6258 fprintf(fplog, "Using static load balancing for the %s direction\n",
6263 for (i = 0; i < nc; i++)
6266 sscanf(size_string, "%lf%n", &dbl, &n);
6269 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6278 fprintf(fplog, "Relative cell sizes:");
6280 for (i = 0; i < nc; i++)
6285 fprintf(fplog, " %5.3f", slb_frac[i]);
6290 fprintf(fplog, "\n");
6297 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6300 gmx_mtop_ilistloop_t iloop;
6304 iloop = gmx_mtop_ilistloop_init(mtop);
6305 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6307 for (ftype = 0; ftype < F_NRE; ftype++)
6309 if ((interaction_function[ftype].flags & IF_BOND) &&
6312 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6320 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6326 val = getenv(env_var);
6329 if (sscanf(val, "%d", &nst) <= 0)
6335 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6343 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6347 fprintf(stderr, "\n%s\n", warn_string);
6351 fprintf(fplog, "\n%s\n", warn_string);
6355 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6356 t_inputrec *ir, FILE *fplog)
6358 if (ir->ePBC == epbcSCREW &&
6359 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6361 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6364 if (ir->ns_type == ensSIMPLE)
6366 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6369 if (ir->nstlist == 0)
6371 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6374 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6376 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6380 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6385 r = ddbox->box_size[XX];
6386 for (di = 0; di < dd->ndim; di++)
6389 /* Check using the initial average cell size */
6390 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6396 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6397 const char *dlb_opt, gmx_bool bRecordLoad,
6398 unsigned long Flags, t_inputrec *ir)
6406 case 'a': eDLB = edlbAUTO; break;
6407 case 'n': eDLB = edlbNO; break;
6408 case 'y': eDLB = edlbYES; break;
6409 default: gmx_incons("Unknown dlb_opt");
6412 if (Flags & MD_RERUN)
6417 if (!EI_DYNAMICS(ir->eI))
6419 if (eDLB == edlbYES)
6421 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6422 dd_warning(cr, fplog, buf);
6430 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6435 if (Flags & MD_REPRODUCIBLE)
6442 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6446 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6449 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6457 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6462 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6464 /* Decomposition order z,y,x */
6467 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6469 for (dim = DIM-1; dim >= 0; dim--)
6471 if (dd->nc[dim] > 1)
6473 dd->dim[dd->ndim++] = dim;
6479 /* Decomposition order x,y,z */
6480 for (dim = 0; dim < DIM; dim++)
6482 if (dd->nc[dim] > 1)
6484 dd->dim[dd->ndim++] = dim;
6490 static gmx_domdec_comm_t *init_dd_comm()
6492 gmx_domdec_comm_t *comm;
6496 snew(comm->cggl_flag, DIM*2);
6497 snew(comm->cgcm_state, DIM*2);
6498 for (i = 0; i < DIM*2; i++)
6500 comm->cggl_flag_nalloc[i] = 0;
6501 comm->cgcm_state_nalloc[i] = 0;
6504 comm->nalloc_int = 0;
6505 comm->buf_int = NULL;
6507 vec_rvec_init(&comm->vbuf);
6509 comm->n_load_have = 0;
6510 comm->n_load_collect = 0;
6512 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6514 comm->sum_nat[i] = 0;
6518 comm->load_step = 0;
6521 clear_ivec(comm->load_lim);
6528 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6529 unsigned long Flags,
6531 real comm_distance_min, real rconstr,
6532 const char *dlb_opt, real dlb_scale,
6533 const char *sizex, const char *sizey, const char *sizez,
6534 gmx_mtop_t *mtop, t_inputrec *ir,
6535 matrix box, rvec *x,
6537 int *npme_x, int *npme_y)
6540 gmx_domdec_comm_t *comm;
6543 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6550 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6555 dd->comm = init_dd_comm();
6557 snew(comm->cggl_flag, DIM*2);
6558 snew(comm->cgcm_state, DIM*2);
6560 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6561 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6563 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6564 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6565 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6566 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6567 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6568 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6569 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6570 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6572 dd->pme_recv_f_alloc = 0;
6573 dd->pme_recv_f_buf = NULL;
6575 if (dd->bSendRecv2 && fplog)
6577 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");
6583 fprintf(fplog, "Will load balance based on FLOP count\n");
6585 if (comm->eFlop > 1)
6587 srand(1+cr->nodeid);
6589 comm->bRecordLoad = TRUE;
6593 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6597 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6599 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6602 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6604 dd->bGridJump = comm->bDynLoadBal;
6605 comm->bPMELoadBalDLBLimits = FALSE;
6607 if (comm->nstSortCG)
6611 if (comm->nstSortCG == 1)
6613 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6617 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6621 snew(comm->sort, 1);
6627 fprintf(fplog, "Will not sort the charge groups\n");
6631 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6633 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6634 if (comm->bInterCGBondeds)
6636 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6640 comm->bInterCGMultiBody = FALSE;
6643 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6644 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6646 if (ir->rlistlong == 0)
6648 /* Set the cut-off to some very large value,
6649 * so we don't need if statements everywhere in the code.
6650 * We use sqrt, since the cut-off is squared in some places.
6652 comm->cutoff = GMX_CUTOFF_INF;
6656 comm->cutoff = ir->rlistlong;
6658 comm->cutoff_mbody = 0;
6660 comm->cellsize_limit = 0;
6661 comm->bBondComm = FALSE;
6663 if (comm->bInterCGBondeds)
6665 if (comm_distance_min > 0)
6667 comm->cutoff_mbody = comm_distance_min;
6668 if (Flags & MD_DDBONDCOMM)
6670 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6674 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6676 r_bonded_limit = comm->cutoff_mbody;
6678 else if (ir->bPeriodicMols)
6680 /* Can not easily determine the required cut-off */
6681 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");
6682 comm->cutoff_mbody = comm->cutoff/2;
6683 r_bonded_limit = comm->cutoff_mbody;
6689 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6690 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6692 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6693 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6695 /* We use an initial margin of 10% for the minimum cell size,
6696 * except when we are just below the non-bonded cut-off.
6698 if (Flags & MD_DDBONDCOMM)
6700 if (max(r_2b, r_mb) > comm->cutoff)
6702 r_bonded = max(r_2b, r_mb);
6703 r_bonded_limit = 1.1*r_bonded;
6704 comm->bBondComm = TRUE;
6709 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6711 /* We determine cutoff_mbody later */
6715 /* No special bonded communication,
6716 * simply increase the DD cut-off.
6718 r_bonded_limit = 1.1*max(r_2b, r_mb);
6719 comm->cutoff_mbody = r_bonded_limit;
6720 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6723 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6727 "Minimum cell size due to bonded interactions: %.3f nm\n",
6728 comm->cellsize_limit);
6732 if (dd->bInterCGcons && rconstr <= 0)
6734 /* There is a cell size limit due to the constraints (P-LINCS) */
6735 rconstr = constr_r_max(fplog, mtop, ir);
6739 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6741 if (rconstr > comm->cellsize_limit)
6743 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6747 else if (rconstr > 0 && fplog)
6749 /* Here we do not check for dd->bInterCGcons,
6750 * because one can also set a cell size limit for virtual sites only
6751 * and at this point we don't know yet if there are intercg v-sites.
6754 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6757 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6759 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6763 copy_ivec(nc, dd->nc);
6764 set_dd_dim(fplog, dd);
6765 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6767 if (cr->npmenodes == -1)
6771 acs = average_cellsize_min(dd, ddbox);
6772 if (acs < comm->cellsize_limit)
6776 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6778 gmx_fatal_collective(FARGS, cr, NULL,
6779 "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",
6780 acs, comm->cellsize_limit);
6785 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6787 /* We need to choose the optimal DD grid and possibly PME nodes */
6788 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6789 comm->eDLB != edlbNO, dlb_scale,
6790 comm->cellsize_limit, comm->cutoff,
6791 comm->bInterCGBondeds);
6793 if (dd->nc[XX] == 0)
6795 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6796 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6797 !bC ? "-rdd" : "-rcon",
6798 comm->eDLB != edlbNO ? " or -dds" : "",
6799 bC ? " or your LINCS settings" : "");
6801 gmx_fatal_collective(FARGS, cr, NULL,
6802 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6804 "Look in the log file for details on the domain decomposition",
6805 cr->nnodes-cr->npmenodes, limit, buf);
6807 set_dd_dim(fplog, dd);
6813 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6814 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6817 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6818 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6820 gmx_fatal_collective(FARGS, cr, NULL,
6821 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6822 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6824 if (cr->npmenodes > dd->nnodes)
6826 gmx_fatal_collective(FARGS, cr, NULL,
6827 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6829 if (cr->npmenodes > 0)
6831 comm->npmenodes = cr->npmenodes;
6835 comm->npmenodes = dd->nnodes;
6838 if (EEL_PME(ir->coulombtype))
6840 /* The following choices should match those
6841 * in comm_cost_est in domdec_setup.c.
6842 * Note that here the checks have to take into account
6843 * that the decomposition might occur in a different order than xyz
6844 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6845 * in which case they will not match those in comm_cost_est,
6846 * but since that is mainly for testing purposes that's fine.
6848 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6849 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6850 getenv("GMX_PMEONEDD") == NULL)
6852 comm->npmedecompdim = 2;
6853 comm->npmenodes_x = dd->nc[XX];
6854 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6858 /* In case nc is 1 in both x and y we could still choose to
6859 * decompose pme in y instead of x, but we use x for simplicity.
6861 comm->npmedecompdim = 1;
6862 if (dd->dim[0] == YY)
6864 comm->npmenodes_x = 1;
6865 comm->npmenodes_y = comm->npmenodes;
6869 comm->npmenodes_x = comm->npmenodes;
6870 comm->npmenodes_y = 1;
6875 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6876 comm->npmenodes_x, comm->npmenodes_y, 1);
6881 comm->npmedecompdim = 0;
6882 comm->npmenodes_x = 0;
6883 comm->npmenodes_y = 0;
6886 /* Technically we don't need both of these,
6887 * but it simplifies code not having to recalculate it.
6889 *npme_x = comm->npmenodes_x;
6890 *npme_y = comm->npmenodes_y;
6892 snew(comm->slb_frac, DIM);
6893 if (comm->eDLB == edlbNO)
6895 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6896 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6897 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6900 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6902 if (comm->bBondComm || comm->eDLB != edlbNO)
6904 /* Set the bonded communication distance to halfway
6905 * the minimum and the maximum,
6906 * since the extra communication cost is nearly zero.
6908 acs = average_cellsize_min(dd, ddbox);
6909 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6910 if (comm->eDLB != edlbNO)
6912 /* Check if this does not limit the scaling */
6913 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6915 if (!comm->bBondComm)
6917 /* Without bBondComm do not go beyond the n.b. cut-off */
6918 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6919 if (comm->cellsize_limit >= comm->cutoff)
6921 /* We don't loose a lot of efficieny
6922 * when increasing it to the n.b. cut-off.
6923 * It can even be slightly faster, because we need
6924 * less checks for the communication setup.
6926 comm->cutoff_mbody = comm->cutoff;
6929 /* Check if we did not end up below our original limit */
6930 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6932 if (comm->cutoff_mbody > comm->cellsize_limit)
6934 comm->cellsize_limit = comm->cutoff_mbody;
6937 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6942 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6943 "cellsize limit %f\n",
6944 comm->bBondComm, comm->cellsize_limit);
6949 check_dd_restrictions(cr, dd, ir, fplog);
6952 comm->partition_step = INT_MIN;
6955 clear_dd_cycle_counts(dd);
6960 static void set_dlb_limits(gmx_domdec_t *dd)
6965 for (d = 0; d < dd->ndim; d++)
6967 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6968 dd->comm->cellsize_min[dd->dim[d]] =
6969 dd->comm->cellsize_min_dlb[dd->dim[d]];
6974 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6977 gmx_domdec_comm_t *comm;
6987 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);
6990 cellsize_min = comm->cellsize_min[dd->dim[0]];
6991 for (d = 1; d < dd->ndim; d++)
6993 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6996 if (cellsize_min < comm->cellsize_limit*1.05)
6998 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");
7000 /* Change DLB from "auto" to "no". */
7001 comm->eDLB = edlbNO;
7006 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7007 comm->bDynLoadBal = TRUE;
7008 dd->bGridJump = TRUE;
7012 /* We can set the required cell size info here,
7013 * so we do not need to communicate this.
7014 * The grid is completely uniform.
7016 for (d = 0; d < dd->ndim; d++)
7020 comm->load[d].sum_m = comm->load[d].sum;
7022 nc = dd->nc[dd->dim[d]];
7023 for (i = 0; i < nc; i++)
7025 comm->root[d]->cell_f[i] = i/(real)nc;
7028 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7029 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7032 comm->root[d]->cell_f[nc] = 1.0;
7037 static char *init_bLocalCG(gmx_mtop_t *mtop)
7042 ncg = ncg_mtop(mtop);
7043 snew(bLocalCG, ncg);
7044 for (cg = 0; cg < ncg; cg++)
7046 bLocalCG[cg] = FALSE;
7052 void dd_init_bondeds(FILE *fplog,
7053 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7055 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7057 gmx_domdec_comm_t *comm;
7061 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7065 if (comm->bBondComm)
7067 /* Communicate atoms beyond the cut-off for bonded interactions */
7070 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7072 comm->bLocalCG = init_bLocalCG(mtop);
7076 /* Only communicate atoms based on cut-off */
7077 comm->cglink = NULL;
7078 comm->bLocalCG = NULL;
7082 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7084 gmx_bool bDynLoadBal, real dlb_scale,
7087 gmx_domdec_comm_t *comm;
7102 fprintf(fplog, "The maximum number of communication pulses is:");
7103 for (d = 0; d < dd->ndim; d++)
7105 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7107 fprintf(fplog, "\n");
7108 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7109 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7110 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7111 for (d = 0; d < DIM; d++)
7115 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7122 comm->cellsize_min_dlb[d]/
7123 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7125 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7128 fprintf(fplog, "\n");
7132 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7133 fprintf(fplog, "The initial number of communication pulses is:");
7134 for (d = 0; d < dd->ndim; d++)
7136 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7138 fprintf(fplog, "\n");
7139 fprintf(fplog, "The initial domain decomposition cell size is:");
7140 for (d = 0; d < DIM; d++)
7144 fprintf(fplog, " %c %.2f nm",
7145 dim2char(d), dd->comm->cellsize_min[d]);
7148 fprintf(fplog, "\n\n");
7151 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7153 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7154 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7155 "non-bonded interactions", "", comm->cutoff);
7159 limit = dd->comm->cellsize_limit;
7163 if (dynamic_dd_box(ddbox, ir))
7165 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7167 limit = dd->comm->cellsize_min[XX];
7168 for (d = 1; d < DIM; d++)
7170 limit = min(limit, dd->comm->cellsize_min[d]);
7174 if (comm->bInterCGBondeds)
7176 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7177 "two-body bonded interactions", "(-rdd)",
7178 max(comm->cutoff, comm->cutoff_mbody));
7179 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7180 "multi-body bonded interactions", "(-rdd)",
7181 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7185 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7186 "virtual site constructions", "(-rcon)", limit);
7188 if (dd->constraint_comm)
7190 sprintf(buf, "atoms separated by up to %d constraints",
7192 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7193 buf, "(-rcon)", limit);
7195 fprintf(fplog, "\n");
7201 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7203 const t_inputrec *ir,
7204 const gmx_ddbox_t *ddbox)
7206 gmx_domdec_comm_t *comm;
7207 int d, dim, npulse, npulse_d_max, npulse_d;
7212 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7214 /* Determine the maximum number of comm. pulses in one dimension */
7216 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7218 /* Determine the maximum required number of grid pulses */
7219 if (comm->cellsize_limit >= comm->cutoff)
7221 /* Only a single pulse is required */
7224 else if (!bNoCutOff && comm->cellsize_limit > 0)
7226 /* We round down slightly here to avoid overhead due to the latency
7227 * of extra communication calls when the cut-off
7228 * would be only slightly longer than the cell size.
7229 * Later cellsize_limit is redetermined,
7230 * so we can not miss interactions due to this rounding.
7232 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7236 /* There is no cell size limit */
7237 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7240 if (!bNoCutOff && npulse > 1)
7242 /* See if we can do with less pulses, based on dlb_scale */
7244 for (d = 0; d < dd->ndim; d++)
7247 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7248 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7249 npulse_d_max = max(npulse_d_max, npulse_d);
7251 npulse = min(npulse, npulse_d_max);
7254 /* This env var can override npulse */
7255 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7262 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7263 for (d = 0; d < dd->ndim; d++)
7265 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7266 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7267 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7268 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7269 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7271 comm->bVacDLBNoLimit = FALSE;
7275 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7276 if (!comm->bVacDLBNoLimit)
7278 comm->cellsize_limit = max(comm->cellsize_limit,
7279 comm->cutoff/comm->maxpulse);
7281 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7282 /* Set the minimum cell size for each DD dimension */
7283 for (d = 0; d < dd->ndim; d++)
7285 if (comm->bVacDLBNoLimit ||
7286 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7288 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7292 comm->cellsize_min_dlb[dd->dim[d]] =
7293 comm->cutoff/comm->cd[d].np_dlb;
7296 if (comm->cutoff_mbody <= 0)
7298 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7300 if (comm->bDynLoadBal)
7306 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7308 /* If each molecule is a single charge group
7309 * or we use domain decomposition for each periodic dimension,
7310 * we do not need to take pbc into account for the bonded interactions.
7312 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7315 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7318 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7319 t_inputrec *ir, gmx_ddbox_t *ddbox)
7321 gmx_domdec_comm_t *comm;
7327 /* Initialize the thread data.
7328 * This can not be done in init_domain_decomposition,
7329 * as the numbers of threads is determined later.
7331 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7334 snew(comm->dth, comm->nth);
7337 if (EEL_PME(ir->coulombtype))
7339 init_ddpme(dd, &comm->ddpme[0], 0);
7340 if (comm->npmedecompdim >= 2)
7342 init_ddpme(dd, &comm->ddpme[1], 1);
7347 comm->npmenodes = 0;
7348 if (dd->pme_nodeid >= 0)
7350 gmx_fatal_collective(FARGS, NULL, dd,
7351 "Can not have separate PME nodes without PME electrostatics");
7357 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7359 if (comm->eDLB != edlbNO)
7361 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7364 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7365 if (comm->eDLB == edlbAUTO)
7369 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7371 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7374 if (ir->ePBC == epbcNONE)
7376 vol_frac = 1 - 1/(double)dd->nnodes;
7381 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7385 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7387 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7389 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7392 static gmx_bool test_dd_cutoff(t_commrec *cr,
7393 t_state *state, t_inputrec *ir,
7404 set_ddbox(dd, FALSE, cr, ir, state->box,
7405 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7409 for (d = 0; d < dd->ndim; d++)
7413 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7414 if (dynamic_dd_box(&ddbox, ir))
7416 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7419 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7421 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7422 dd->comm->cd[d].np_dlb > 0)
7424 if (np > dd->comm->cd[d].np_dlb)
7429 /* If a current local cell size is smaller than the requested
7430 * cut-off, we could still fix it, but this gets very complicated.
7431 * Without fixing here, we might actually need more checks.
7433 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7440 if (dd->comm->eDLB != edlbNO)
7442 /* If DLB is not active yet, we don't need to check the grid jumps.
7443 * Actually we shouldn't, because then the grid jump data is not set.
7445 if (dd->comm->bDynLoadBal &&
7446 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7451 gmx_sumi(1, &LocallyLimited, cr);
7453 if (LocallyLimited > 0)
7462 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7465 gmx_bool bCutoffAllowed;
7467 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7471 cr->dd->comm->cutoff = cutoff_req;
7474 return bCutoffAllowed;
7477 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7479 gmx_domdec_comm_t *comm;
7481 comm = cr->dd->comm;
7483 /* Turn on the DLB limiting (might have been on already) */
7484 comm->bPMELoadBalDLBLimits = TRUE;
7486 /* Change the cut-off limit */
7487 comm->PMELoadBal_max_cutoff = comm->cutoff;
7490 static void merge_cg_buffers(int ncell,
7491 gmx_domdec_comm_dim_t *cd, int pulse,
7493 int *index_gl, int *recv_i,
7494 rvec *cg_cm, rvec *recv_vr,
7496 cginfo_mb_t *cginfo_mb, int *cginfo)
7498 gmx_domdec_ind_t *ind, *ind_p;
7499 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7500 int shift, shift_at;
7502 ind = &cd->ind[pulse];
7504 /* First correct the already stored data */
7505 shift = ind->nrecv[ncell];
7506 for (cell = ncell-1; cell >= 0; cell--)
7508 shift -= ind->nrecv[cell];
7511 /* Move the cg's present from previous grid pulses */
7512 cg0 = ncg_cell[ncell+cell];
7513 cg1 = ncg_cell[ncell+cell+1];
7514 cgindex[cg1+shift] = cgindex[cg1];
7515 for (cg = cg1-1; cg >= cg0; cg--)
7517 index_gl[cg+shift] = index_gl[cg];
7518 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7519 cgindex[cg+shift] = cgindex[cg];
7520 cginfo[cg+shift] = cginfo[cg];
7522 /* Correct the already stored send indices for the shift */
7523 for (p = 1; p <= pulse; p++)
7525 ind_p = &cd->ind[p];
7527 for (c = 0; c < cell; c++)
7529 cg0 += ind_p->nsend[c];
7531 cg1 = cg0 + ind_p->nsend[cell];
7532 for (cg = cg0; cg < cg1; cg++)
7534 ind_p->index[cg] += shift;
7540 /* Merge in the communicated buffers */
7544 for (cell = 0; cell < ncell; cell++)
7546 cg1 = ncg_cell[ncell+cell+1] + shift;
7549 /* Correct the old cg indices */
7550 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7552 cgindex[cg+1] += shift_at;
7555 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7557 /* Copy this charge group from the buffer */
7558 index_gl[cg1] = recv_i[cg0];
7559 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7560 /* Add it to the cgindex */
7561 cg_gl = index_gl[cg1];
7562 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7563 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7564 cgindex[cg1+1] = cgindex[cg1] + nat;
7569 shift += ind->nrecv[cell];
7570 ncg_cell[ncell+cell+1] = cg1;
7574 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7575 int nzone, int cg0, const int *cgindex)
7579 /* Store the atom block boundaries for easy copying of communication buffers
7582 for (zone = 0; zone < nzone; zone++)
7584 for (p = 0; p < cd->np; p++)
7586 cd->ind[p].cell2at0[zone] = cgindex[cg];
7587 cg += cd->ind[p].nrecv[zone];
7588 cd->ind[p].cell2at1[zone] = cgindex[cg];
7593 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7599 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7601 if (!bLocalCG[link->a[i]])
7610 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7612 real c[DIM][4]; /* the corners for the non-bonded communication */
7613 real cr0; /* corner for rounding */
7614 real cr1[4]; /* corners for rounding */
7615 real bc[DIM]; /* corners for bounded communication */
7616 real bcr1; /* corner for rounding for bonded communication */
7619 /* Determine the corners of the domain(s) we are communicating with */
7621 set_dd_corners(const gmx_domdec_t *dd,
7622 int dim0, int dim1, int dim2,
7626 const gmx_domdec_comm_t *comm;
7627 const gmx_domdec_zones_t *zones;
7632 zones = &comm->zones;
7634 /* Keep the compiler happy */
7638 /* The first dimension is equal for all cells */
7639 c->c[0][0] = comm->cell_x0[dim0];
7642 c->bc[0] = c->c[0][0];
7647 /* This cell row is only seen from the first row */
7648 c->c[1][0] = comm->cell_x0[dim1];
7649 /* All rows can see this row */
7650 c->c[1][1] = comm->cell_x0[dim1];
7653 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7656 /* For the multi-body distance we need the maximum */
7657 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7660 /* Set the upper-right corner for rounding */
7661 c->cr0 = comm->cell_x1[dim0];
7666 for (j = 0; j < 4; j++)
7668 c->c[2][j] = comm->cell_x0[dim2];
7672 /* Use the maximum of the i-cells that see a j-cell */
7673 for (i = 0; i < zones->nizone; i++)
7675 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7681 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7687 /* For the multi-body distance we need the maximum */
7688 c->bc[2] = comm->cell_x0[dim2];
7689 for (i = 0; i < 2; i++)
7691 for (j = 0; j < 2; j++)
7693 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7699 /* Set the upper-right corner for rounding */
7700 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7701 * Only cell (0,0,0) can see cell 7 (1,1,1)
7703 c->cr1[0] = comm->cell_x1[dim1];
7704 c->cr1[3] = comm->cell_x1[dim1];
7707 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7710 /* For the multi-body distance we need the maximum */
7711 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7718 /* Determine which cg's we need to send in this pulse from this zone */
7720 get_zone_pulse_cgs(gmx_domdec_t *dd,
7721 int zonei, int zone,
7723 const int *index_gl,
7725 int dim, int dim_ind,
7726 int dim0, int dim1, int dim2,
7727 real r_comm2, real r_bcomm2,
7731 real skew_fac2_d, real skew_fac_01,
7732 rvec *v_d, rvec *v_0, rvec *v_1,
7733 const dd_corners_t *c,
7735 gmx_bool bDistBonded,
7741 gmx_domdec_ind_t *ind,
7742 int **ibuf, int *ibuf_nalloc,
7748 gmx_domdec_comm_t *comm;
7750 gmx_bool bDistMB_pulse;
7752 real r2, rb2, r, tric_sh;
7755 int nsend_z, nsend, nat;
7759 bScrew = (dd->bScrewPBC && dim == XX);
7761 bDistMB_pulse = (bDistMB && bDistBonded);
7767 for (cg = cg0; cg < cg1; cg++)
7771 if (tric_dist[dim_ind] == 0)
7773 /* Rectangular direction, easy */
7774 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7781 r = cg_cm[cg][dim] - c->bc[dim_ind];
7787 /* Rounding gives at most a 16% reduction
7788 * in communicated atoms
7790 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7792 r = cg_cm[cg][dim0] - c->cr0;
7793 /* This is the first dimension, so always r >= 0 */
7800 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7802 r = cg_cm[cg][dim1] - c->cr1[zone];
7809 r = cg_cm[cg][dim1] - c->bcr1;
7819 /* Triclinic direction, more complicated */
7822 /* Rounding, conservative as the skew_fac multiplication
7823 * will slightly underestimate the distance.
7825 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7827 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7828 for (i = dim0+1; i < DIM; i++)
7830 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7832 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7835 rb[dim0] = rn[dim0];
7838 /* Take care that the cell planes along dim0 might not
7839 * be orthogonal to those along dim1 and dim2.
7841 for (i = 1; i <= dim_ind; i++)
7844 if (normal[dim0][dimd] > 0)
7846 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7849 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7854 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7856 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7858 for (i = dim1+1; i < DIM; i++)
7860 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7862 rn[dim1] += tric_sh;
7865 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7866 /* Take care of coupling of the distances
7867 * to the planes along dim0 and dim1 through dim2.
7869 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7870 /* Take care that the cell planes along dim1
7871 * might not be orthogonal to that along dim2.
7873 if (normal[dim1][dim2] > 0)
7875 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7881 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7884 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7885 /* Take care of coupling of the distances
7886 * to the planes along dim0 and dim1 through dim2.
7888 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7889 /* Take care that the cell planes along dim1
7890 * might not be orthogonal to that along dim2.
7892 if (normal[dim1][dim2] > 0)
7894 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7899 /* The distance along the communication direction */
7900 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7902 for (i = dim+1; i < DIM; i++)
7904 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7909 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7910 /* Take care of coupling of the distances
7911 * to the planes along dim0 and dim1 through dim2.
7913 if (dim_ind == 1 && zonei == 1)
7915 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7921 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7924 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7925 /* Take care of coupling of the distances
7926 * to the planes along dim0 and dim1 through dim2.
7928 if (dim_ind == 1 && zonei == 1)
7930 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7938 ((bDistMB && rb2 < r_bcomm2) ||
7939 (bDist2B && r2 < r_bcomm2)) &&
7941 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7942 missing_link(comm->cglink, index_gl[cg],
7945 /* Make an index to the local charge groups */
7946 if (nsend+1 > ind->nalloc)
7948 ind->nalloc = over_alloc_large(nsend+1);
7949 srenew(ind->index, ind->nalloc);
7951 if (nsend+1 > *ibuf_nalloc)
7953 *ibuf_nalloc = over_alloc_large(nsend+1);
7954 srenew(*ibuf, *ibuf_nalloc);
7956 ind->index[nsend] = cg;
7957 (*ibuf)[nsend] = index_gl[cg];
7959 vec_rvec_check_alloc(vbuf, nsend+1);
7961 if (dd->ci[dim] == 0)
7963 /* Correct cg_cm for pbc */
7964 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7967 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7968 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7973 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7976 nat += cgindex[cg+1] - cgindex[cg];
7982 *nsend_z_ptr = nsend_z;
7985 static void setup_dd_communication(gmx_domdec_t *dd,
7986 matrix box, gmx_ddbox_t *ddbox,
7987 t_forcerec *fr, t_state *state, rvec **f)
7989 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7990 int nzone, nzone_send, zone, zonei, cg0, cg1;
7991 int c, i, j, cg, cg_gl, nrcg;
7992 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7993 gmx_domdec_comm_t *comm;
7994 gmx_domdec_zones_t *zones;
7995 gmx_domdec_comm_dim_t *cd;
7996 gmx_domdec_ind_t *ind;
7997 cginfo_mb_t *cginfo_mb;
7998 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7999 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
8000 dd_corners_t corners;
8002 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8003 real skew_fac2_d, skew_fac_01;
8010 fprintf(debug, "Setting up DD communication\n");
8015 switch (fr->cutoff_scheme)
8024 gmx_incons("unimplemented");
8028 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8030 dim = dd->dim[dim_ind];
8032 /* Check if we need to use triclinic distances */
8033 tric_dist[dim_ind] = 0;
8034 for (i = 0; i <= dim_ind; i++)
8036 if (ddbox->tric_dir[dd->dim[i]])
8038 tric_dist[dim_ind] = 1;
8043 bBondComm = comm->bBondComm;
8045 /* Do we need to determine extra distances for multi-body bondeds? */
8046 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8048 /* Do we need to determine extra distances for only two-body bondeds? */
8049 bDist2B = (bBondComm && !bDistMB);
8051 r_comm2 = sqr(comm->cutoff);
8052 r_bcomm2 = sqr(comm->cutoff_mbody);
8056 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8059 zones = &comm->zones;
8062 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8063 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8065 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8067 /* Triclinic stuff */
8068 normal = ddbox->normal;
8072 v_0 = ddbox->v[dim0];
8073 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8075 /* Determine the coupling coefficient for the distances
8076 * to the cell planes along dim0 and dim1 through dim2.
8077 * This is required for correct rounding.
8080 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8083 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8089 v_1 = ddbox->v[dim1];
8092 zone_cg_range = zones->cg_range;
8093 index_gl = dd->index_gl;
8094 cgindex = dd->cgindex;
8095 cginfo_mb = fr->cginfo_mb;
8097 zone_cg_range[0] = 0;
8098 zone_cg_range[1] = dd->ncg_home;
8099 comm->zone_ncg1[0] = dd->ncg_home;
8100 pos_cg = dd->ncg_home;
8102 nat_tot = dd->nat_home;
8104 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8106 dim = dd->dim[dim_ind];
8107 cd = &comm->cd[dim_ind];
8109 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8111 /* No pbc in this dimension, the first node should not comm. */
8119 v_d = ddbox->v[dim];
8120 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8122 cd->bInPlace = TRUE;
8123 for (p = 0; p < cd->np; p++)
8125 /* Only atoms communicated in the first pulse are used
8126 * for multi-body bonded interactions or for bBondComm.
8128 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8133 for (zone = 0; zone < nzone_send; zone++)
8135 if (tric_dist[dim_ind] && dim_ind > 0)
8137 /* Determine slightly more optimized skew_fac's
8139 * This reduces the number of communicated atoms
8140 * by about 10% for 3D DD of rhombic dodecahedra.
8142 for (dimd = 0; dimd < dim; dimd++)
8144 sf2_round[dimd] = 1;
8145 if (ddbox->tric_dir[dimd])
8147 for (i = dd->dim[dimd]+1; i < DIM; i++)
8149 /* If we are shifted in dimension i
8150 * and the cell plane is tilted forward
8151 * in dimension i, skip this coupling.
8153 if (!(zones->shift[nzone+zone][i] &&
8154 ddbox->v[dimd][i][dimd] >= 0))
8157 sqr(ddbox->v[dimd][i][dimd]);
8160 sf2_round[dimd] = 1/sf2_round[dimd];
8165 zonei = zone_perm[dim_ind][zone];
8168 /* Here we permutate the zones to obtain a convenient order
8169 * for neighbor searching
8171 cg0 = zone_cg_range[zonei];
8172 cg1 = zone_cg_range[zonei+1];
8176 /* Look only at the cg's received in the previous grid pulse
8178 cg1 = zone_cg_range[nzone+zone+1];
8179 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8182 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8183 for (th = 0; th < comm->nth; th++)
8185 gmx_domdec_ind_t *ind_p;
8186 int **ibuf_p, *ibuf_nalloc_p;
8188 int *nsend_p, *nat_p;
8194 /* Thread 0 writes in the comm buffers */
8196 ibuf_p = &comm->buf_int;
8197 ibuf_nalloc_p = &comm->nalloc_int;
8198 vbuf_p = &comm->vbuf;
8201 nsend_zone_p = &ind->nsend[zone];
8205 /* Other threads write into temp buffers */
8206 ind_p = &comm->dth[th].ind;
8207 ibuf_p = &comm->dth[th].ibuf;
8208 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8209 vbuf_p = &comm->dth[th].vbuf;
8210 nsend_p = &comm->dth[th].nsend;
8211 nat_p = &comm->dth[th].nat;
8212 nsend_zone_p = &comm->dth[th].nsend_zone;
8214 comm->dth[th].nsend = 0;
8215 comm->dth[th].nat = 0;
8216 comm->dth[th].nsend_zone = 0;
8226 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8227 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8230 /* Get the cg's for this pulse in this zone */
8231 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8233 dim, dim_ind, dim0, dim1, dim2,
8236 normal, skew_fac2_d, skew_fac_01,
8237 v_d, v_0, v_1, &corners, sf2_round,
8238 bDistBonded, bBondComm,
8242 ibuf_p, ibuf_nalloc_p,
8248 /* Append data of threads>=1 to the communication buffers */
8249 for (th = 1; th < comm->nth; th++)
8251 dd_comm_setup_work_t *dth;
8254 dth = &comm->dth[th];
8256 ns1 = nsend + dth->nsend_zone;
8257 if (ns1 > ind->nalloc)
8259 ind->nalloc = over_alloc_dd(ns1);
8260 srenew(ind->index, ind->nalloc);
8262 if (ns1 > comm->nalloc_int)
8264 comm->nalloc_int = over_alloc_dd(ns1);
8265 srenew(comm->buf_int, comm->nalloc_int);
8267 if (ns1 > comm->vbuf.nalloc)
8269 comm->vbuf.nalloc = over_alloc_dd(ns1);
8270 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8273 for (i = 0; i < dth->nsend_zone; i++)
8275 ind->index[nsend] = dth->ind.index[i];
8276 comm->buf_int[nsend] = dth->ibuf[i];
8277 copy_rvec(dth->vbuf.v[i],
8278 comm->vbuf.v[nsend]);
8282 ind->nsend[zone] += dth->nsend_zone;
8285 /* Clear the counts in case we do not have pbc */
8286 for (zone = nzone_send; zone < nzone; zone++)
8288 ind->nsend[zone] = 0;
8290 ind->nsend[nzone] = nsend;
8291 ind->nsend[nzone+1] = nat;
8292 /* Communicate the number of cg's and atoms to receive */
8293 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8294 ind->nsend, nzone+2,
8295 ind->nrecv, nzone+2);
8297 /* The rvec buffer is also required for atom buffers of size nsend
8298 * in dd_move_x and dd_move_f.
8300 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8304 /* We can receive in place if only the last zone is not empty */
8305 for (zone = 0; zone < nzone-1; zone++)
8307 if (ind->nrecv[zone] > 0)
8309 cd->bInPlace = FALSE;
8314 /* The int buffer is only required here for the cg indices */
8315 if (ind->nrecv[nzone] > comm->nalloc_int2)
8317 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8318 srenew(comm->buf_int2, comm->nalloc_int2);
8320 /* The rvec buffer is also required for atom buffers
8321 * of size nrecv in dd_move_x and dd_move_f.
8323 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8324 vec_rvec_check_alloc(&comm->vbuf2, i);
8328 /* Make space for the global cg indices */
8329 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8330 || dd->cg_nalloc == 0)
8332 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8333 srenew(index_gl, dd->cg_nalloc);
8334 srenew(cgindex, dd->cg_nalloc+1);
8336 /* Communicate the global cg indices */
8339 recv_i = index_gl + pos_cg;
8343 recv_i = comm->buf_int2;
8345 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8346 comm->buf_int, nsend,
8347 recv_i, ind->nrecv[nzone]);
8349 /* Make space for cg_cm */
8350 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8351 if (fr->cutoff_scheme == ecutsGROUP)
8359 /* Communicate cg_cm */
8362 recv_vr = cg_cm + pos_cg;
8366 recv_vr = comm->vbuf2.v;
8368 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8369 comm->vbuf.v, nsend,
8370 recv_vr, ind->nrecv[nzone]);
8372 /* Make the charge group index */
8375 zone = (p == 0 ? 0 : nzone - 1);
8376 while (zone < nzone)
8378 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8380 cg_gl = index_gl[pos_cg];
8381 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8382 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8383 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8386 /* Update the charge group presence,
8387 * so we can use it in the next pass of the loop.
8389 comm->bLocalCG[cg_gl] = TRUE;
8395 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8398 zone_cg_range[nzone+zone] = pos_cg;
8403 /* This part of the code is never executed with bBondComm. */
8404 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8405 index_gl, recv_i, cg_cm, recv_vr,
8406 cgindex, fr->cginfo_mb, fr->cginfo);
8407 pos_cg += ind->nrecv[nzone];
8409 nat_tot += ind->nrecv[nzone+1];
8413 /* Store the atom block for easy copying of communication buffers */
8414 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8418 dd->index_gl = index_gl;
8419 dd->cgindex = cgindex;
8421 dd->ncg_tot = zone_cg_range[zones->n];
8422 dd->nat_tot = nat_tot;
8423 comm->nat[ddnatHOME] = dd->nat_home;
8424 for (i = ddnatZONE; i < ddnatNR; i++)
8426 comm->nat[i] = dd->nat_tot;
8431 /* We don't need to update cginfo, since that was alrady done above.
8432 * So we pass NULL for the forcerec.
8434 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8435 NULL, comm->bLocalCG);
8440 fprintf(debug, "Finished setting up DD communication, zones:");
8441 for (c = 0; c < zones->n; c++)
8443 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8445 fprintf(debug, "\n");
8449 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8453 for (c = 0; c < zones->nizone; c++)
8455 zones->izone[c].cg1 = zones->cg_range[c+1];
8456 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8457 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8461 static void set_zones_size(gmx_domdec_t *dd,
8462 matrix box, const gmx_ddbox_t *ddbox,
8463 int zone_start, int zone_end)
8465 gmx_domdec_comm_t *comm;
8466 gmx_domdec_zones_t *zones;
8468 int z, zi, zj0, zj1, d, dim;
8471 real size_j, add_tric;
8476 zones = &comm->zones;
8478 /* Do we need to determine extra distances for multi-body bondeds? */
8479 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8481 for (z = zone_start; z < zone_end; z++)
8483 /* Copy cell limits to zone limits.
8484 * Valid for non-DD dims and non-shifted dims.
8486 copy_rvec(comm->cell_x0, zones->size[z].x0);
8487 copy_rvec(comm->cell_x1, zones->size[z].x1);
8490 for (d = 0; d < dd->ndim; d++)
8494 for (z = 0; z < zones->n; z++)
8496 /* With a staggered grid we have different sizes
8497 * for non-shifted dimensions.
8499 if (dd->bGridJump && zones->shift[z][dim] == 0)
8503 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8504 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8508 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8509 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8515 rcmbs = comm->cutoff_mbody;
8516 if (ddbox->tric_dir[dim])
8518 rcs /= ddbox->skew_fac[dim];
8519 rcmbs /= ddbox->skew_fac[dim];
8522 /* Set the lower limit for the shifted zone dimensions */
8523 for (z = zone_start; z < zone_end; z++)
8525 if (zones->shift[z][dim] > 0)
8528 if (!dd->bGridJump || d == 0)
8530 zones->size[z].x0[dim] = comm->cell_x1[dim];
8531 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8535 /* Here we take the lower limit of the zone from
8536 * the lowest domain of the zone below.
8540 zones->size[z].x0[dim] =
8541 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8547 zones->size[z].x0[dim] =
8548 zones->size[zone_perm[2][z-4]].x0[dim];
8552 zones->size[z].x0[dim] =
8553 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8556 /* A temporary limit, is updated below */
8557 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8561 for (zi = 0; zi < zones->nizone; zi++)
8563 if (zones->shift[zi][dim] == 0)
8565 /* This takes the whole zone into account.
8566 * With multiple pulses this will lead
8567 * to a larger zone then strictly necessary.
8569 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8570 zones->size[zi].x1[dim]+rcmbs);
8578 /* Loop over the i-zones to set the upper limit of each
8581 for (zi = 0; zi < zones->nizone; zi++)
8583 if (zones->shift[zi][dim] == 0)
8585 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8587 if (zones->shift[z][dim] > 0)
8589 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8590 zones->size[zi].x1[dim]+rcs);
8597 for (z = zone_start; z < zone_end; z++)
8599 /* Initialization only required to keep the compiler happy */
8600 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8603 /* To determine the bounding box for a zone we need to find
8604 * the extreme corners of 4, 2 or 1 corners.
8606 nc = 1 << (ddbox->npbcdim - 1);
8608 for (c = 0; c < nc; c++)
8610 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8614 corner[YY] = zones->size[z].x0[YY];
8618 corner[YY] = zones->size[z].x1[YY];
8622 corner[ZZ] = zones->size[z].x0[ZZ];
8626 corner[ZZ] = zones->size[z].x1[ZZ];
8628 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8630 /* With 1D domain decomposition the cg's are not in
8631 * the triclinic box, but triclinic x-y and rectangular y-z.
8632 * Shift y back, so it will later end up at 0.
8634 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8636 /* Apply the triclinic couplings */
8637 for (i = YY; i < ddbox->npbcdim; i++)
8639 for (j = XX; j < i; j++)
8641 corner[j] += corner[i]*box[i][j]/box[i][i];
8646 copy_rvec(corner, corner_min);
8647 copy_rvec(corner, corner_max);
8651 for (i = 0; i < DIM; i++)
8653 corner_min[i] = min(corner_min[i], corner[i]);
8654 corner_max[i] = max(corner_max[i], corner[i]);
8658 /* Copy the extreme cornes without offset along x */
8659 for (i = 0; i < DIM; i++)
8661 zones->size[z].bb_x0[i] = corner_min[i];
8662 zones->size[z].bb_x1[i] = corner_max[i];
8664 /* Add the offset along x */
8665 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8666 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8669 if (zone_start == 0)
8672 for (dim = 0; dim < DIM; dim++)
8674 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8676 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8681 for (z = zone_start; z < zone_end; z++)
8683 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8685 zones->size[z].x0[XX], zones->size[z].x1[XX],
8686 zones->size[z].x0[YY], zones->size[z].x1[YY],
8687 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8688 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8690 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8691 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8692 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8697 static int comp_cgsort(const void *a, const void *b)
8701 gmx_cgsort_t *cga, *cgb;
8702 cga = (gmx_cgsort_t *)a;
8703 cgb = (gmx_cgsort_t *)b;
8705 comp = cga->nsc - cgb->nsc;
8708 comp = cga->ind_gl - cgb->ind_gl;
8714 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8719 /* Order the data */
8720 for (i = 0; i < n; i++)
8722 buf[i] = a[sort[i].ind];
8725 /* Copy back to the original array */
8726 for (i = 0; i < n; i++)
8732 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8737 /* Order the data */
8738 for (i = 0; i < n; i++)
8740 copy_rvec(v[sort[i].ind], buf[i]);
8743 /* Copy back to the original array */
8744 for (i = 0; i < n; i++)
8746 copy_rvec(buf[i], v[i]);
8750 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8753 int a, atot, cg, cg0, cg1, i;
8755 if (cgindex == NULL)
8757 /* Avoid the useless loop of the atoms within a cg */
8758 order_vec_cg(ncg, sort, v, buf);
8763 /* Order the data */
8765 for (cg = 0; cg < ncg; cg++)
8767 cg0 = cgindex[sort[cg].ind];
8768 cg1 = cgindex[sort[cg].ind+1];
8769 for (i = cg0; i < cg1; i++)
8771 copy_rvec(v[i], buf[a]);
8777 /* Copy back to the original array */
8778 for (a = 0; a < atot; a++)
8780 copy_rvec(buf[a], v[a]);
8784 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8785 int nsort_new, gmx_cgsort_t *sort_new,
8786 gmx_cgsort_t *sort1)
8790 /* The new indices are not very ordered, so we qsort them */
8791 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8793 /* sort2 is already ordered, so now we can merge the two arrays */
8797 while (i2 < nsort2 || i_new < nsort_new)
8801 sort1[i1++] = sort_new[i_new++];
8803 else if (i_new == nsort_new)
8805 sort1[i1++] = sort2[i2++];
8807 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8808 (sort2[i2].nsc == sort_new[i_new].nsc &&
8809 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8811 sort1[i1++] = sort2[i2++];
8815 sort1[i1++] = sort_new[i_new++];
8820 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8822 gmx_domdec_sort_t *sort;
8823 gmx_cgsort_t *cgsort, *sort_i;
8824 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8825 int sort_last, sort_skip;
8827 sort = dd->comm->sort;
8829 a = fr->ns.grid->cell_index;
8831 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8833 if (ncg_home_old >= 0)
8835 /* The charge groups that remained in the same ns grid cell
8836 * are completely ordered. So we can sort efficiently by sorting
8837 * the charge groups that did move into the stationary list.
8842 for (i = 0; i < dd->ncg_home; i++)
8844 /* Check if this cg did not move to another node */
8847 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8849 /* This cg is new on this node or moved ns grid cell */
8850 if (nsort_new >= sort->sort_new_nalloc)
8852 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8853 srenew(sort->sort_new, sort->sort_new_nalloc);
8855 sort_i = &(sort->sort_new[nsort_new++]);
8859 /* This cg did not move */
8860 sort_i = &(sort->sort2[nsort2++]);
8862 /* Sort on the ns grid cell indices
8863 * and the global topology index.
8864 * index_gl is irrelevant with cell ns,
8865 * but we set it here anyhow to avoid a conditional.
8868 sort_i->ind_gl = dd->index_gl[i];
8875 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8878 /* Sort efficiently */
8879 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8884 cgsort = sort->sort;
8886 for (i = 0; i < dd->ncg_home; i++)
8888 /* Sort on the ns grid cell indices
8889 * and the global topology index
8891 cgsort[i].nsc = a[i];
8892 cgsort[i].ind_gl = dd->index_gl[i];
8894 if (cgsort[i].nsc < moved)
8901 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8903 /* Determine the order of the charge groups using qsort */
8904 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8910 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8913 int ncg_new, i, *a, na;
8915 sort = dd->comm->sort->sort;
8917 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8920 for (i = 0; i < na; i++)
8924 sort[ncg_new].ind = a[i];
8932 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
8935 gmx_domdec_sort_t *sort;
8936 gmx_cgsort_t *cgsort, *sort_i;
8938 int ncg_new, i, *ibuf, cgsize;
8941 sort = dd->comm->sort;
8943 if (dd->ncg_home > sort->sort_nalloc)
8945 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8946 srenew(sort->sort, sort->sort_nalloc);
8947 srenew(sort->sort2, sort->sort_nalloc);
8949 cgsort = sort->sort;
8951 switch (fr->cutoff_scheme)
8954 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8957 ncg_new = dd_sort_order_nbnxn(dd, fr);
8960 gmx_incons("unimplemented");
8964 /* We alloc with the old size, since cgindex is still old */
8965 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8966 vbuf = dd->comm->vbuf.v;
8970 cgindex = dd->cgindex;
8977 /* Remove the charge groups which are no longer at home here */
8978 dd->ncg_home = ncg_new;
8981 fprintf(debug, "Set the new home charge group count to %d\n",
8985 /* Reorder the state */
8986 for (i = 0; i < estNR; i++)
8988 if (EST_DISTR(i) && (state->flags & (1<<i)))
8993 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8996 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8999 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9002 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9006 case estDISRE_INITF:
9007 case estDISRE_RM3TAV:
9008 case estORIRE_INITF:
9010 /* No ordering required */
9013 gmx_incons("Unknown state entry encountered in dd_sort_state");
9018 if (fr->cutoff_scheme == ecutsGROUP)
9021 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9024 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9026 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9027 srenew(sort->ibuf, sort->ibuf_nalloc);
9030 /* Reorder the global cg index */
9031 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9032 /* Reorder the cginfo */
9033 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9034 /* Rebuild the local cg index */
9038 for (i = 0; i < dd->ncg_home; i++)
9040 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9041 ibuf[i+1] = ibuf[i] + cgsize;
9043 for (i = 0; i < dd->ncg_home+1; i++)
9045 dd->cgindex[i] = ibuf[i];
9050 for (i = 0; i < dd->ncg_home+1; i++)
9055 /* Set the home atom number */
9056 dd->nat_home = dd->cgindex[dd->ncg_home];
9058 if (fr->cutoff_scheme == ecutsVERLET)
9060 /* The atoms are now exactly in grid order, update the grid order */
9061 nbnxn_set_atomorder(fr->nbv->nbs);
9065 /* Copy the sorted ns cell indices back to the ns grid struct */
9066 for (i = 0; i < dd->ncg_home; i++)
9068 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9070 fr->ns.grid->nr = dd->ncg_home;
9074 static void add_dd_statistics(gmx_domdec_t *dd)
9076 gmx_domdec_comm_t *comm;
9081 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9083 comm->sum_nat[ddnat-ddnatZONE] +=
9084 comm->nat[ddnat] - comm->nat[ddnat-1];
9089 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9091 gmx_domdec_comm_t *comm;
9096 /* Reset all the statistics and counters for total run counting */
9097 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9099 comm->sum_nat[ddnat-ddnatZONE] = 0;
9103 comm->load_step = 0;
9106 clear_ivec(comm->load_lim);
9111 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9113 gmx_domdec_comm_t *comm;
9117 comm = cr->dd->comm;
9119 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9126 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");
9128 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9130 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9135 " av. #atoms communicated per step for force: %d x %.1f\n",
9139 if (cr->dd->vsite_comm)
9142 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9143 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9148 if (cr->dd->constraint_comm)
9151 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9152 1 + ir->nLincsIter, av);
9156 gmx_incons(" Unknown type for DD statistics");
9159 fprintf(fplog, "\n");
9161 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9163 print_dd_load_av(fplog, cr->dd);
9167 void dd_partition_system(FILE *fplog,
9168 gmx_large_int_t step,
9170 gmx_bool bMasterState,
9172 t_state *state_global,
9173 gmx_mtop_t *top_global,
9175 t_state *state_local,
9178 gmx_localtop_t *top_local,
9181 gmx_shellfc_t shellfc,
9182 gmx_constr_t constr,
9184 gmx_wallcycle_t wcycle,
9188 gmx_domdec_comm_t *comm;
9189 gmx_ddbox_t ddbox = {0};
9191 gmx_large_int_t step_pcoupl;
9192 rvec cell_ns_x0, cell_ns_x1;
9193 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9194 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9195 gmx_bool bRedist, bSortCG, bResortAll;
9196 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9203 bBoxChanged = (bMasterState || DEFORM(*ir));
9204 if (ir->epc != epcNO)
9206 /* With nstpcouple > 1 pressure coupling happens.
9207 * one step after calculating the pressure.
9208 * Box scaling happens at the end of the MD step,
9209 * after the DD partitioning.
9210 * We therefore have to do DLB in the first partitioning
9211 * after an MD step where P-coupling occured.
9212 * We need to determine the last step in which p-coupling occurred.
9213 * MRS -- need to validate this for vv?
9218 step_pcoupl = step - 1;
9222 step_pcoupl = ((step - 1)/n)*n + 1;
9224 if (step_pcoupl >= comm->partition_step)
9230 bNStGlobalComm = (step % nstglobalcomm == 0);
9232 if (!comm->bDynLoadBal)
9238 /* Should we do dynamic load balacing this step?
9239 * Since it requires (possibly expensive) global communication,
9240 * we might want to do DLB less frequently.
9242 if (bBoxChanged || ir->epc != epcNO)
9244 bDoDLB = bBoxChanged;
9248 bDoDLB = bNStGlobalComm;
9252 /* Check if we have recorded loads on the nodes */
9253 if (comm->bRecordLoad && dd_load_count(comm))
9255 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9257 /* Check if we should use DLB at the second partitioning
9258 * and every 100 partitionings,
9259 * so the extra communication cost is negligible.
9261 n = max(100, nstglobalcomm);
9262 bCheckDLB = (comm->n_load_collect == 0 ||
9263 comm->n_load_have % n == n-1);
9270 /* Print load every nstlog, first and last step to the log file */
9271 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9272 comm->n_load_collect == 0 ||
9274 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9276 /* Avoid extra communication due to verbose screen output
9277 * when nstglobalcomm is set.
9279 if (bDoDLB || bLogLoad || bCheckDLB ||
9280 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9282 get_load_distribution(dd, wcycle);
9287 dd_print_load(fplog, dd, step-1);
9291 dd_print_load_verbose(dd);
9294 comm->n_load_collect++;
9298 /* Since the timings are node dependent, the master decides */
9302 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9305 fprintf(debug, "step %s, imb loss %f\n",
9306 gmx_step_str(step, sbuf),
9307 dd_force_imb_perf_loss(dd));
9310 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9313 turn_on_dlb(fplog, cr, step);
9318 comm->n_load_have++;
9321 cgs_gl = &comm->cgs_gl;
9326 /* Clear the old state */
9327 clear_dd_indices(dd, 0, 0);
9330 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9331 TRUE, cgs_gl, state_global->x, &ddbox);
9333 get_cg_distribution(fplog, step, dd, cgs_gl,
9334 state_global->box, &ddbox, state_global->x);
9336 dd_distribute_state(dd, cgs_gl,
9337 state_global, state_local, f);
9339 dd_make_local_cgs(dd, &top_local->cgs);
9341 /* Ensure that we have space for the new distribution */
9342 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9344 if (fr->cutoff_scheme == ecutsGROUP)
9346 calc_cgcm(fplog, 0, dd->ncg_home,
9347 &top_local->cgs, state_local->x, fr->cg_cm);
9350 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9352 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9354 else if (state_local->ddp_count != dd->ddp_count)
9356 if (state_local->ddp_count > dd->ddp_count)
9358 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9361 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9363 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);
9366 /* Clear the old state */
9367 clear_dd_indices(dd, 0, 0);
9369 /* Build the new indices */
9370 rebuild_cgindex(dd, cgs_gl->index, state_local);
9371 make_dd_indices(dd, cgs_gl->index, 0);
9372 ncgindex_set = dd->ncg_home;
9374 if (fr->cutoff_scheme == ecutsGROUP)
9376 /* Redetermine the cg COMs */
9377 calc_cgcm(fplog, 0, dd->ncg_home,
9378 &top_local->cgs, state_local->x, fr->cg_cm);
9381 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9383 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9385 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9386 TRUE, &top_local->cgs, state_local->x, &ddbox);
9388 bRedist = comm->bDynLoadBal;
9392 /* We have the full state, only redistribute the cgs */
9394 /* Clear the non-home indices */
9395 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9398 /* Avoid global communication for dim's without pbc and -gcom */
9399 if (!bNStGlobalComm)
9401 copy_rvec(comm->box0, ddbox.box0 );
9402 copy_rvec(comm->box_size, ddbox.box_size);
9404 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9405 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9410 /* For dim's without pbc and -gcom */
9411 copy_rvec(ddbox.box0, comm->box0 );
9412 copy_rvec(ddbox.box_size, comm->box_size);
9414 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9417 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9419 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9422 /* Check if we should sort the charge groups */
9423 if (comm->nstSortCG > 0)
9425 bSortCG = (bMasterState ||
9426 (bRedist && (step % comm->nstSortCG == 0)));
9433 ncg_home_old = dd->ncg_home;
9438 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9440 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9442 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9444 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9447 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9449 &comm->cell_x0, &comm->cell_x1,
9450 dd->ncg_home, fr->cg_cm,
9451 cell_ns_x0, cell_ns_x1, &grid_density);
9455 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9458 switch (fr->cutoff_scheme)
9461 copy_ivec(fr->ns.grid->n, ncells_old);
9462 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9463 state_local->box, cell_ns_x0, cell_ns_x1,
9464 fr->rlistlong, grid_density);
9467 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9470 gmx_incons("unimplemented");
9472 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9473 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9477 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9479 /* Sort the state on charge group position.
9480 * This enables exact restarts from this step.
9481 * It also improves performance by about 15% with larger numbers
9482 * of atoms per node.
9485 /* Fill the ns grid with the home cell,
9486 * so we can sort with the indices.
9488 set_zones_ncg_home(dd);
9490 switch (fr->cutoff_scheme)
9493 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9495 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9497 comm->zones.size[0].bb_x0,
9498 comm->zones.size[0].bb_x1,
9500 comm->zones.dens_zone0,
9503 ncg_moved, bRedist ? comm->moved : NULL,
9504 fr->nbv->grp[eintLocal].kernel_type,
9505 fr->nbv->grp[eintLocal].nbat);
9507 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9510 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9511 0, dd->ncg_home, fr->cg_cm);
9513 copy_ivec(fr->ns.grid->n, ncells_new);
9516 gmx_incons("unimplemented");
9519 bResortAll = bMasterState;
9521 /* Check if we can user the old order and ns grid cell indices
9522 * of the charge groups to sort the charge groups efficiently.
9524 if (ncells_new[XX] != ncells_old[XX] ||
9525 ncells_new[YY] != ncells_old[YY] ||
9526 ncells_new[ZZ] != ncells_old[ZZ])
9533 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9534 gmx_step_str(step, sbuf), dd->ncg_home);
9536 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9537 bResortAll ? -1 : ncg_home_old);
9538 /* Rebuild all the indices */
9539 ga2la_clear(dd->ga2la);
9542 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9545 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9547 /* Setup up the communication and communicate the coordinates */
9548 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9550 /* Set the indices */
9551 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9553 /* Set the charge group boundaries for neighbor searching */
9554 set_cg_boundaries(&comm->zones);
9556 if (fr->cutoff_scheme == ecutsVERLET)
9558 set_zones_size(dd, state_local->box, &ddbox,
9559 bSortCG ? 1 : 0, comm->zones.n);
9562 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9565 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9566 -1,state_local->x,state_local->box);
9569 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9571 /* Extract a local topology from the global topology */
9572 for (i = 0; i < dd->ndim; i++)
9574 np[dd->dim[i]] = comm->cd[i].np;
9576 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9577 comm->cellsize_min, np,
9579 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9580 vsite, top_global, top_local);
9582 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9584 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9586 /* Set up the special atom communication */
9587 n = comm->nat[ddnatZONE];
9588 for (i = ddnatZONE+1; i < ddnatNR; i++)
9593 if (vsite && vsite->n_intercg_vsite)
9595 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9599 if (dd->bInterCGcons || dd->bInterCGsettles)
9601 /* Only for inter-cg constraints we need special code */
9602 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9603 constr, ir->nProjOrder,
9604 top_local->idef.il);
9608 gmx_incons("Unknown special atom type setup");
9613 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9615 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9617 /* Make space for the extra coordinates for virtual site
9618 * or constraint communication.
9620 state_local->natoms = comm->nat[ddnatNR-1];
9621 if (state_local->natoms > state_local->nalloc)
9623 dd_realloc_state(state_local, f, state_local->natoms);
9626 if (fr->bF_NoVirSum)
9628 if (vsite && vsite->n_intercg_vsite)
9630 nat_f_novirsum = comm->nat[ddnatVSITE];
9634 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9636 nat_f_novirsum = dd->nat_tot;
9640 nat_f_novirsum = dd->nat_home;
9649 /* Set the number of atoms required for the force calculation.
9650 * Forces need to be constrained when using a twin-range setup
9651 * or with energy minimization. For simple simulations we could
9652 * avoid some allocation, zeroing and copying, but this is
9653 * probably not worth the complications ande checking.
9655 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9656 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9658 /* We make the all mdatoms up to nat_tot_con.
9659 * We could save some work by only setting invmass
9660 * between nat_tot and nat_tot_con.
9662 /* This call also sets the new number of home particles to dd->nat_home */
9663 atoms2md(top_global, ir,
9664 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9666 /* Now we have the charges we can sort the FE interactions */
9667 dd_sort_local_top(dd, mdatoms, top_local);
9671 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9672 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9677 /* Make the local shell stuff, currently no communication is done */
9678 make_local_shells(cr, mdatoms, shellfc);
9681 if (ir->implicit_solvent)
9683 make_local_gb(cr, fr->born, ir->gb_algorithm);
9686 setup_bonded_threading(fr, &top_local->idef);
9688 if (!(cr->duty & DUTY_PME))
9690 /* Send the charges to our PME only node */
9691 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9692 mdatoms->chargeA, mdatoms->chargeB,
9693 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9698 set_constraints(constr, top_local, ir, mdatoms, cr);
9701 if (ir->ePull != epullNO)
9703 /* Update the local pull groups */
9704 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9709 /* Update the local rotation groups */
9710 dd_make_local_rotation_groups(dd, ir->rot);
9714 add_dd_statistics(dd);
9716 /* Make sure we only count the cycles for this DD partitioning */
9717 clear_dd_cycle_counts(dd);
9719 /* Because the order of the atoms might have changed since
9720 * the last vsite construction, we need to communicate the constructing
9721 * atom coordinates again (for spreading the forces this MD step).
9723 dd_move_x_vsites(dd, state_local->box, state_local->x);
9725 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9727 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9729 dd_move_x(dd, state_local->box, state_local->x);
9730 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9731 -1, state_local->x, state_local->box);
9734 /* Store the partitioning step */
9735 comm->partition_step = step;
9737 /* Increase the DD partitioning counter */
9739 /* The state currently matches this DD partitioning count, store it */
9740 state_local->ddp_count = dd->ddp_count;
9743 /* The DD master node knows the complete cg distribution,
9744 * store the count so we can possibly skip the cg info communication.
9746 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9749 if (comm->DD_debug > 0)
9751 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9752 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9753 "after partitioning");