2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2008
5 * Copyright (c) 2012,2013, by the GROMACS development team, led by
6 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
7 * others, as listed in the AUTHORS file in the top-level source
8 * directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gmx_fatal.h"
49 #include "gmx_fatal_collective.h"
52 #include "domdec_network.h"
55 #include "chargegroup.h"
64 #include "pull_rotation.h"
65 #include "gmx_wallcycle.h"
69 #include "mtop_util.h"
71 #include "gmx_ga2la.h"
73 #include "nbnxn_search.h"
75 #include "gmx_omp_nthreads.h"
76 #include "gpu_utils.h"
85 #define DDRANK(dd, rank) (rank)
86 #define DDMASTERRANK(dd) (dd->masterrank)
88 typedef struct gmx_domdec_master
90 /* The cell boundaries */
92 /* The global charge group division */
93 int *ncg; /* Number of home charge groups for each node */
94 int *index; /* Index of nnodes+1 into cg */
95 int *cg; /* Global charge group index */
96 int *nat; /* Number of home atoms for each node. */
97 int *ibuf; /* Buffer for communication */
98 rvec *vbuf; /* Buffer for state scattering and gathering */
99 } gmx_domdec_master_t;
103 /* The numbers of charge groups to send and receive for each cell
104 * that requires communication, the last entry contains the total
105 * number of atoms that needs to be communicated.
107 int nsend[DD_MAXIZONE+2];
108 int nrecv[DD_MAXIZONE+2];
109 /* The charge groups to send */
112 /* The atom range for non-in-place communication */
113 int cell2at0[DD_MAXIZONE];
114 int cell2at1[DD_MAXIZONE];
119 int np; /* Number of grid pulses in this dimension */
120 int np_dlb; /* For dlb, for use with edlbAUTO */
121 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
123 gmx_bool bInPlace; /* Can we communicate in place? */
124 } gmx_domdec_comm_dim_t;
128 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
129 real *cell_f; /* State var.: cell boundaries, box relative */
130 real *old_cell_f; /* Temp. var.: old cell size */
131 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
132 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
133 real *bound_min; /* Temp. var.: lower limit for cell boundary */
134 real *bound_max; /* Temp. var.: upper limit for cell boundary */
135 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
136 real *buf_ncd; /* Temp. var. */
139 #define DD_NLOAD_MAX 9
141 /* Here floats are accurate enough, since these variables
142 * only influence the load balancing, not the actual MD results.
169 gmx_cgsort_t *sort_new;
181 /* This enum determines the order of the coordinates.
182 * ddnatHOME and ddnatZONE should be first and second,
183 * the others can be ordered as wanted.
186 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
190 edlbAUTO, edlbNO, edlbYES, edlbNR
192 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
196 int dim; /* The dimension */
197 gmx_bool dim_match; /* Tells if DD and PME dims match */
198 int nslab; /* The number of PME slabs in this dimension */
199 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
200 int *pp_min; /* The minimum pp node location, size nslab */
201 int *pp_max; /* The maximum pp node location,size nslab */
202 int maxshift; /* The maximum shift for coordinate redistribution in PME */
207 real min0; /* The minimum bottom of this zone */
208 real max1; /* The maximum top of this zone */
209 real min1; /* The minimum top of this zone */
210 real mch0; /* The maximum bottom communicaton height for this zone */
211 real mch1; /* The maximum top communicaton height for this zone */
212 real p1_0; /* The bottom value of the first cell in this zone */
213 real p1_1; /* The top value of the first cell in this zone */
218 gmx_domdec_ind_t ind;
225 } dd_comm_setup_work_t;
227 typedef struct gmx_domdec_comm
229 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
230 * unless stated otherwise.
233 /* The number of decomposition dimensions for PME, 0: no PME */
235 /* The number of nodes doing PME (PP/PME or only PME) */
239 /* The communication setup including the PME only nodes */
240 gmx_bool bCartesianPP_PME;
243 int *pmenodes; /* size npmenodes */
244 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
245 * but with bCartesianPP_PME */
246 gmx_ddpme_t ddpme[2];
248 /* The DD particle-particle nodes only */
249 gmx_bool bCartesianPP;
250 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
252 /* The global charge groups */
255 /* Should we sort the cgs */
257 gmx_domdec_sort_t *sort;
259 /* Are there charge groups? */
262 /* Are there bonded and multi-body interactions between charge groups? */
263 gmx_bool bInterCGBondeds;
264 gmx_bool bInterCGMultiBody;
266 /* Data for the optional bonded interaction atom communication range */
273 /* Are we actually using DLB? */
274 gmx_bool bDynLoadBal;
276 /* Cell sizes for static load balancing, first index cartesian */
279 /* The width of the communicated boundaries */
282 /* The minimum cell size (including triclinic correction) */
284 /* For dlb, for use with edlbAUTO */
285 rvec cellsize_min_dlb;
286 /* The lower limit for the DD cell size with DLB */
288 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
289 gmx_bool bVacDLBNoLimit;
291 /* With PME load balancing we set limits on DLB */
292 gmx_bool bPMELoadBalDLBLimits;
293 /* DLB needs to take into account that we want to allow this maximum
294 * cut-off (for PME load balancing), this could limit cell boundaries.
296 real PMELoadBal_max_cutoff;
298 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
300 /* box0 and box_size are required with dim's without pbc and -gcom */
304 /* The cell boundaries */
308 /* The old location of the cell boundaries, to check cg displacements */
312 /* The communication setup and charge group boundaries for the zones */
313 gmx_domdec_zones_t zones;
315 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
316 * cell boundaries of neighboring cells for dynamic load balancing.
318 gmx_ddzone_t zone_d1[2];
319 gmx_ddzone_t zone_d2[2][2];
321 /* The coordinate/force communication setup and indices */
322 gmx_domdec_comm_dim_t cd[DIM];
323 /* The maximum number of cells to communicate with in one dimension */
326 /* Which cg distribution is stored on the master node */
327 int master_cg_ddp_count;
329 /* The number of cg's received from the direct neighbors */
330 int zone_ncg1[DD_MAXZONE];
332 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
335 /* Array for signalling if atoms have moved to another domain */
339 /* Communication buffer for general use */
343 /* Communication buffer for general use */
346 /* Temporary storage for thread parallel communication setup */
348 dd_comm_setup_work_t *dth;
350 /* Communication buffers only used with multiple grid pulses */
355 /* Communication buffers for local redistribution */
357 int cggl_flag_nalloc[DIM*2];
359 int cgcm_state_nalloc[DIM*2];
361 /* Cell sizes for dynamic load balancing */
362 gmx_domdec_root_t **root;
366 real cell_f_max0[DIM];
367 real cell_f_min1[DIM];
369 /* Stuff for load communication */
370 gmx_bool bRecordLoad;
371 gmx_domdec_load_t *load;
372 int nrank_gpu_shared;
374 MPI_Comm *mpi_comm_load;
375 MPI_Comm mpi_comm_gpu_shared;
378 /* Maximum DLB scaling per load balancing step in percent */
382 float cycl[ddCyclNr];
383 int cycl_n[ddCyclNr];
384 float cycl_max[ddCyclNr];
385 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
389 /* Have often have did we have load measurements */
391 /* Have often have we collected the load measurements */
395 double sum_nat[ddnatNR-ddnatZONE];
405 /* The last partition step */
406 gmx_large_int_t partition_step;
414 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
417 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
418 #define DD_FLAG_NRCG 65535
419 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
420 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
422 /* Zone permutation required to obtain consecutive charge groups
423 * for neighbor searching.
425 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
427 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
428 * components see only j zones with that component 0.
431 /* The DD zone order */
432 static const ivec dd_zo[DD_MAXZONE] =
433 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
438 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
443 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
448 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
450 /* Factors used to avoid problems due to rounding issues */
451 #define DD_CELL_MARGIN 1.0001
452 #define DD_CELL_MARGIN2 1.00005
453 /* Factor to account for pressure scaling during nstlist steps */
454 #define DD_PRES_SCALE_MARGIN 1.02
456 /* Allowed performance loss before we DLB or warn */
457 #define DD_PERF_LOSS 0.05
459 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
461 /* Use separate MPI send and receive commands
462 * when nnodes <= GMX_DD_NNODES_SENDRECV.
463 * This saves memory (and some copying for small nnodes).
464 * For high parallelization scatter and gather calls are used.
466 #define GMX_DD_NNODES_SENDRECV 4
470 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
472 static void index2xyz(ivec nc,int ind,ivec xyz)
474 xyz[XX] = ind % nc[XX];
475 xyz[YY] = (ind / nc[XX]) % nc[YY];
476 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
480 /* This order is required to minimize the coordinate communication in PME
481 * which uses decomposition in the x direction.
483 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
485 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
487 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
488 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
489 xyz[ZZ] = ind % nc[ZZ];
492 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
497 ddindex = dd_index(dd->nc, c);
498 if (dd->comm->bCartesianPP_PME)
500 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
502 else if (dd->comm->bCartesianPP)
505 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
516 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
518 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
521 int ddglatnr(gmx_domdec_t *dd, int i)
531 if (i >= dd->comm->nat[ddnatNR-1])
533 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
535 atnr = dd->gatindex[i] + 1;
541 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
543 return &dd->comm->cgs_gl;
546 static void vec_rvec_init(vec_rvec_t *v)
552 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
556 v->nalloc = over_alloc_dd(n);
557 srenew(v->v, v->nalloc);
561 void dd_store_state(gmx_domdec_t *dd, t_state *state)
565 if (state->ddp_count != dd->ddp_count)
567 gmx_incons("The state does not the domain decomposition state");
570 state->ncg_gl = dd->ncg_home;
571 if (state->ncg_gl > state->cg_gl_nalloc)
573 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
574 srenew(state->cg_gl, state->cg_gl_nalloc);
576 for (i = 0; i < state->ncg_gl; i++)
578 state->cg_gl[i] = dd->index_gl[i];
581 state->ddp_count_cg_gl = dd->ddp_count;
584 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
586 return &dd->comm->zones;
589 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
590 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
592 gmx_domdec_zones_t *zones;
595 zones = &dd->comm->zones;
598 while (icg >= zones->izone[izone].cg1)
607 else if (izone < zones->nizone)
609 *jcg0 = zones->izone[izone].jcg0;
613 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
614 icg, izone, zones->nizone);
617 *jcg1 = zones->izone[izone].jcg1;
619 for (d = 0; d < dd->ndim; d++)
622 shift0[dim] = zones->izone[izone].shift0[dim];
623 shift1[dim] = zones->izone[izone].shift1[dim];
624 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
626 /* A conservative approach, this can be optimized */
633 int dd_natoms_vsite(gmx_domdec_t *dd)
635 return dd->comm->nat[ddnatVSITE];
638 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
640 *at_start = dd->comm->nat[ddnatCON-1];
641 *at_end = dd->comm->nat[ddnatCON];
644 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
646 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
647 int *index, *cgindex;
648 gmx_domdec_comm_t *comm;
649 gmx_domdec_comm_dim_t *cd;
650 gmx_domdec_ind_t *ind;
651 rvec shift = {0, 0, 0}, *buf, *rbuf;
652 gmx_bool bPBC, bScrew;
656 cgindex = dd->cgindex;
661 nat_tot = dd->nat_home;
662 for (d = 0; d < dd->ndim; d++)
664 bPBC = (dd->ci[dd->dim[d]] == 0);
665 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
668 copy_rvec(box[dd->dim[d]], shift);
671 for (p = 0; p < cd->np; p++)
678 for (i = 0; i < ind->nsend[nzone]; i++)
680 at0 = cgindex[index[i]];
681 at1 = cgindex[index[i]+1];
682 for (j = at0; j < at1; j++)
684 copy_rvec(x[j], buf[n]);
691 for (i = 0; i < ind->nsend[nzone]; i++)
693 at0 = cgindex[index[i]];
694 at1 = cgindex[index[i]+1];
695 for (j = at0; j < at1; j++)
697 /* We need to shift the coordinates */
698 rvec_add(x[j], shift, buf[n]);
705 for (i = 0; i < ind->nsend[nzone]; i++)
707 at0 = cgindex[index[i]];
708 at1 = cgindex[index[i]+1];
709 for (j = at0; j < at1; j++)
712 buf[n][XX] = x[j][XX] + shift[XX];
714 * This operation requires a special shift force
715 * treatment, which is performed in calc_vir.
717 buf[n][YY] = box[YY][YY] - x[j][YY];
718 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
730 rbuf = comm->vbuf2.v;
732 /* Send and receive the coordinates */
733 dd_sendrecv_rvec(dd, d, dddirBackward,
734 buf, ind->nsend[nzone+1],
735 rbuf, ind->nrecv[nzone+1]);
739 for (zone = 0; zone < nzone; zone++)
741 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
743 copy_rvec(rbuf[j], x[i]);
748 nat_tot += ind->nrecv[nzone+1];
754 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
756 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
757 int *index, *cgindex;
758 gmx_domdec_comm_t *comm;
759 gmx_domdec_comm_dim_t *cd;
760 gmx_domdec_ind_t *ind;
764 gmx_bool bPBC, bScrew;
768 cgindex = dd->cgindex;
773 nzone = comm->zones.n/2;
774 nat_tot = dd->nat_tot;
775 for (d = dd->ndim-1; d >= 0; d--)
777 bPBC = (dd->ci[dd->dim[d]] == 0);
778 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
779 if (fshift == NULL && !bScrew)
783 /* Determine which shift vector we need */
789 for (p = cd->np-1; p >= 0; p--)
792 nat_tot -= ind->nrecv[nzone+1];
799 sbuf = comm->vbuf2.v;
801 for (zone = 0; zone < nzone; zone++)
803 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
805 copy_rvec(f[i], sbuf[j]);
810 /* Communicate the forces */
811 dd_sendrecv_rvec(dd, d, dddirForward,
812 sbuf, ind->nrecv[nzone+1],
813 buf, ind->nsend[nzone+1]);
815 /* Add the received forces */
819 for (i = 0; i < ind->nsend[nzone]; i++)
821 at0 = cgindex[index[i]];
822 at1 = cgindex[index[i]+1];
823 for (j = at0; j < at1; j++)
825 rvec_inc(f[j], buf[n]);
832 for (i = 0; i < ind->nsend[nzone]; i++)
834 at0 = cgindex[index[i]];
835 at1 = cgindex[index[i]+1];
836 for (j = at0; j < at1; j++)
838 rvec_inc(f[j], buf[n]);
839 /* Add this force to the shift force */
840 rvec_inc(fshift[is], buf[n]);
847 for (i = 0; i < ind->nsend[nzone]; i++)
849 at0 = cgindex[index[i]];
850 at1 = cgindex[index[i]+1];
851 for (j = at0; j < at1; j++)
853 /* Rotate the force */
854 f[j][XX] += buf[n][XX];
855 f[j][YY] -= buf[n][YY];
856 f[j][ZZ] -= buf[n][ZZ];
859 /* Add this force to the shift force */
860 rvec_inc(fshift[is], buf[n]);
871 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
873 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
874 int *index, *cgindex;
875 gmx_domdec_comm_t *comm;
876 gmx_domdec_comm_dim_t *cd;
877 gmx_domdec_ind_t *ind;
882 cgindex = dd->cgindex;
884 buf = &comm->vbuf.v[0][0];
887 nat_tot = dd->nat_home;
888 for (d = 0; d < dd->ndim; d++)
891 for (p = 0; p < cd->np; p++)
896 for (i = 0; i < ind->nsend[nzone]; i++)
898 at0 = cgindex[index[i]];
899 at1 = cgindex[index[i]+1];
900 for (j = at0; j < at1; j++)
913 rbuf = &comm->vbuf2.v[0][0];
915 /* Send and receive the coordinates */
916 dd_sendrecv_real(dd, d, dddirBackward,
917 buf, ind->nsend[nzone+1],
918 rbuf, ind->nrecv[nzone+1]);
922 for (zone = 0; zone < nzone; zone++)
924 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
931 nat_tot += ind->nrecv[nzone+1];
937 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
939 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
940 int *index, *cgindex;
941 gmx_domdec_comm_t *comm;
942 gmx_domdec_comm_dim_t *cd;
943 gmx_domdec_ind_t *ind;
948 cgindex = dd->cgindex;
950 buf = &comm->vbuf.v[0][0];
953 nzone = comm->zones.n/2;
954 nat_tot = dd->nat_tot;
955 for (d = dd->ndim-1; d >= 0; d--)
958 for (p = cd->np-1; p >= 0; p--)
961 nat_tot -= ind->nrecv[nzone+1];
968 sbuf = &comm->vbuf2.v[0][0];
970 for (zone = 0; zone < nzone; zone++)
972 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
979 /* Communicate the forces */
980 dd_sendrecv_real(dd, d, dddirForward,
981 sbuf, ind->nrecv[nzone+1],
982 buf, ind->nsend[nzone+1]);
984 /* Add the received forces */
986 for (i = 0; i < ind->nsend[nzone]; i++)
988 at0 = cgindex[index[i]];
989 at1 = cgindex[index[i]+1];
990 for (j = at0; j < at1; j++)
1001 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1003 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",
1005 zone->min0, zone->max1,
1006 zone->mch0, zone->mch0,
1007 zone->p1_0, zone->p1_1);
1011 #define DDZONECOMM_MAXZONE 5
1012 #define DDZONECOMM_BUFSIZE 3
1014 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1015 int ddimind, int direction,
1016 gmx_ddzone_t *buf_s, int n_s,
1017 gmx_ddzone_t *buf_r, int n_r)
1019 #define ZBS DDZONECOMM_BUFSIZE
1020 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1021 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1024 for (i = 0; i < n_s; i++)
1026 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1027 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1028 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1029 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1030 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1031 vbuf_s[i*ZBS+1][2] = 0;
1032 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1033 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1034 vbuf_s[i*ZBS+2][2] = 0;
1037 dd_sendrecv_rvec(dd, ddimind, direction,
1041 for (i = 0; i < n_r; i++)
1043 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1044 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1045 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1046 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1047 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1048 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1049 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1055 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1056 rvec cell_ns_x0, rvec cell_ns_x1)
1058 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1060 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1061 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1062 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1063 rvec extr_s[2], extr_r[2];
1065 real dist_d, c = 0, det;
1066 gmx_domdec_comm_t *comm;
1067 gmx_bool bPBC, bUse;
1071 for (d = 1; d < dd->ndim; d++)
1074 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1075 zp->min0 = cell_ns_x0[dim];
1076 zp->max1 = cell_ns_x1[dim];
1077 zp->min1 = cell_ns_x1[dim];
1078 zp->mch0 = cell_ns_x0[dim];
1079 zp->mch1 = cell_ns_x1[dim];
1080 zp->p1_0 = cell_ns_x0[dim];
1081 zp->p1_1 = cell_ns_x1[dim];
1084 for (d = dd->ndim-2; d >= 0; d--)
1087 bPBC = (dim < ddbox->npbcdim);
1089 /* Use an rvec to store two reals */
1090 extr_s[d][0] = comm->cell_f0[d+1];
1091 extr_s[d][1] = comm->cell_f1[d+1];
1092 extr_s[d][2] = comm->cell_f1[d+1];
1095 /* Store the extremes in the backward sending buffer,
1096 * so the get updated separately from the forward communication.
1098 for (d1 = d; d1 < dd->ndim-1; d1++)
1100 /* We invert the order to be able to use the same loop for buf_e */
1101 buf_s[pos].min0 = extr_s[d1][1];
1102 buf_s[pos].max1 = extr_s[d1][0];
1103 buf_s[pos].min1 = extr_s[d1][2];
1104 buf_s[pos].mch0 = 0;
1105 buf_s[pos].mch1 = 0;
1106 /* Store the cell corner of the dimension we communicate along */
1107 buf_s[pos].p1_0 = comm->cell_x0[dim];
1108 buf_s[pos].p1_1 = 0;
1112 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1115 if (dd->ndim == 3 && d == 0)
1117 buf_s[pos] = comm->zone_d2[0][1];
1119 buf_s[pos] = comm->zone_d1[0];
1123 /* We only need to communicate the extremes
1124 * in the forward direction
1126 npulse = comm->cd[d].np;
1129 /* Take the minimum to avoid double communication */
1130 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1134 /* Without PBC we should really not communicate over
1135 * the boundaries, but implementing that complicates
1136 * the communication setup and therefore we simply
1137 * do all communication, but ignore some data.
1139 npulse_min = npulse;
1141 for (p = 0; p < npulse_min; p++)
1143 /* Communicate the extremes forward */
1144 bUse = (bPBC || dd->ci[dim] > 0);
1146 dd_sendrecv_rvec(dd, d, dddirForward,
1147 extr_s+d, dd->ndim-d-1,
1148 extr_r+d, dd->ndim-d-1);
1152 for (d1 = d; d1 < dd->ndim-1; d1++)
1154 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1155 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1156 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1162 for (p = 0; p < npulse; p++)
1164 /* Communicate all the zone information backward */
1165 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1167 dd_sendrecv_ddzone(dd, d, dddirBackward,
1174 for (d1 = d+1; d1 < dd->ndim; d1++)
1176 /* Determine the decrease of maximum required
1177 * communication height along d1 due to the distance along d,
1178 * this avoids a lot of useless atom communication.
1180 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1182 if (ddbox->tric_dir[dim])
1184 /* c is the off-diagonal coupling between the cell planes
1185 * along directions d and d1.
1187 c = ddbox->v[dim][dd->dim[d1]][dim];
1193 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1196 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1200 /* A negative value signals out of range */
1206 /* Accumulate the extremes over all pulses */
1207 for (i = 0; i < buf_size; i++)
1211 buf_e[i] = buf_r[i];
1217 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1218 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1219 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1222 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1230 if (bUse && dh[d1] >= 0)
1232 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1233 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1236 /* Copy the received buffer to the send buffer,
1237 * to pass the data through with the next pulse.
1239 buf_s[i] = buf_r[i];
1241 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1242 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1244 /* Store the extremes */
1247 for (d1 = d; d1 < dd->ndim-1; d1++)
1249 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1250 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1251 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1255 if (d == 1 || (d == 0 && dd->ndim == 3))
1257 for (i = d; i < 2; i++)
1259 comm->zone_d2[1-d][i] = buf_e[pos];
1265 comm->zone_d1[1] = buf_e[pos];
1275 for (i = 0; i < 2; i++)
1279 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1281 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1282 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1288 for (i = 0; i < 2; i++)
1290 for (j = 0; j < 2; j++)
1294 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1296 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1297 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1301 for (d = 1; d < dd->ndim; d++)
1303 comm->cell_f_max0[d] = extr_s[d-1][0];
1304 comm->cell_f_min1[d] = extr_s[d-1][1];
1307 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1308 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1313 static void dd_collect_cg(gmx_domdec_t *dd,
1314 t_state *state_local)
1316 gmx_domdec_master_t *ma = NULL;
1317 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1320 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1322 /* The master has the correct distribution */
1326 if (state_local->ddp_count == dd->ddp_count)
1328 ncg_home = dd->ncg_home;
1330 nat_home = dd->nat_home;
1332 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1334 cgs_gl = &dd->comm->cgs_gl;
1336 ncg_home = state_local->ncg_gl;
1337 cg = state_local->cg_gl;
1339 for (i = 0; i < ncg_home; i++)
1341 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1346 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1349 buf2[0] = dd->ncg_home;
1350 buf2[1] = dd->nat_home;
1360 /* Collect the charge group and atom counts on the master */
1361 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1366 for (i = 0; i < dd->nnodes; i++)
1368 ma->ncg[i] = ma->ibuf[2*i];
1369 ma->nat[i] = ma->ibuf[2*i+1];
1370 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1373 /* Make byte counts and indices */
1374 for (i = 0; i < dd->nnodes; i++)
1376 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1377 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1381 fprintf(debug, "Initial charge group distribution: ");
1382 for (i = 0; i < dd->nnodes; i++)
1384 fprintf(debug, " %d", ma->ncg[i]);
1386 fprintf(debug, "\n");
1390 /* Collect the charge group indices on the master */
1392 dd->ncg_home*sizeof(int), dd->index_gl,
1393 DDMASTER(dd) ? ma->ibuf : NULL,
1394 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1395 DDMASTER(dd) ? ma->cg : NULL);
1397 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1400 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1403 gmx_domdec_master_t *ma;
1404 int n, i, c, a, nalloc = 0;
1413 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1414 dd->rank, dd->mpi_comm_all);
1419 /* Copy the master coordinates to the global array */
1420 cgs_gl = &dd->comm->cgs_gl;
1422 n = DDMASTERRANK(dd);
1424 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1426 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1428 copy_rvec(lv[a++], v[c]);
1432 for (n = 0; n < dd->nnodes; n++)
1436 if (ma->nat[n] > nalloc)
1438 nalloc = over_alloc_dd(ma->nat[n]);
1439 srenew(buf, nalloc);
1442 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1443 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1446 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1448 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1450 copy_rvec(buf[a++], v[c]);
1459 static void get_commbuffer_counts(gmx_domdec_t *dd,
1460 int **counts, int **disps)
1462 gmx_domdec_master_t *ma;
1467 /* Make the rvec count and displacment arrays */
1469 *disps = ma->ibuf + dd->nnodes;
1470 for (n = 0; n < dd->nnodes; n++)
1472 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1473 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1477 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1480 gmx_domdec_master_t *ma;
1481 int *rcounts = NULL, *disps = NULL;
1490 get_commbuffer_counts(dd, &rcounts, &disps);
1495 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1499 cgs_gl = &dd->comm->cgs_gl;
1502 for (n = 0; n < dd->nnodes; n++)
1504 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1506 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1508 copy_rvec(buf[a++], v[c]);
1515 void dd_collect_vec(gmx_domdec_t *dd,
1516 t_state *state_local, rvec *lv, rvec *v)
1518 gmx_domdec_master_t *ma;
1519 int n, i, c, a, nalloc = 0;
1522 dd_collect_cg(dd, state_local);
1524 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1526 dd_collect_vec_sendrecv(dd, lv, v);
1530 dd_collect_vec_gatherv(dd, lv, v);
1535 void dd_collect_state(gmx_domdec_t *dd,
1536 t_state *state_local, t_state *state)
1540 nh = state->nhchainlength;
1544 for (i = 0; i < efptNR; i++)
1546 state->lambda[i] = state_local->lambda[i];
1548 state->fep_state = state_local->fep_state;
1549 state->veta = state_local->veta;
1550 state->vol0 = state_local->vol0;
1551 copy_mat(state_local->box, state->box);
1552 copy_mat(state_local->boxv, state->boxv);
1553 copy_mat(state_local->svir_prev, state->svir_prev);
1554 copy_mat(state_local->fvir_prev, state->fvir_prev);
1555 copy_mat(state_local->pres_prev, state->pres_prev);
1557 for (i = 0; i < state_local->ngtc; i++)
1559 for (j = 0; j < nh; j++)
1561 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1562 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1564 state->therm_integral[i] = state_local->therm_integral[i];
1566 for (i = 0; i < state_local->nnhpres; i++)
1568 for (j = 0; j < nh; j++)
1570 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1571 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1575 for (est = 0; est < estNR; est++)
1577 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1582 dd_collect_vec(dd, state_local, state_local->x, state->x);
1585 dd_collect_vec(dd, state_local, state_local->v, state->v);
1588 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1591 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1594 if (state->nrngi == 1)
1598 for (i = 0; i < state_local->nrng; i++)
1600 state->ld_rng[i] = state_local->ld_rng[i];
1606 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1607 state_local->ld_rng, state->ld_rng);
1611 if (state->nrngi == 1)
1615 state->ld_rngi[0] = state_local->ld_rngi[0];
1620 dd_gather(dd, sizeof(state->ld_rngi[0]),
1621 state_local->ld_rngi, state->ld_rngi);
1624 case estDISRE_INITF:
1625 case estDISRE_RM3TAV:
1626 case estORIRE_INITF:
1630 gmx_incons("Unknown state entry encountered in dd_collect_state");
1636 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1642 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1645 state->nalloc = over_alloc_dd(nalloc);
1647 for (est = 0; est < estNR; est++)
1649 if (EST_DISTR(est) && (state->flags & (1<<est)))
1654 srenew(state->x, state->nalloc);
1657 srenew(state->v, state->nalloc);
1660 srenew(state->sd_X, state->nalloc);
1663 srenew(state->cg_p, state->nalloc);
1667 case estDISRE_INITF:
1668 case estDISRE_RM3TAV:
1669 case estORIRE_INITF:
1671 /* No reallocation required */
1674 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1681 srenew(*f, state->nalloc);
1685 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1688 if (nalloc > fr->cg_nalloc)
1692 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1694 fr->cg_nalloc = over_alloc_dd(nalloc);
1695 srenew(fr->cginfo, fr->cg_nalloc);
1696 if (fr->cutoff_scheme == ecutsGROUP)
1698 srenew(fr->cg_cm, fr->cg_nalloc);
1701 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1703 /* We don't use charge groups, we use x in state to set up
1704 * the atom communication.
1706 dd_realloc_state(state, f, nalloc);
1710 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1713 gmx_domdec_master_t *ma;
1714 int n, i, c, a, nalloc = 0;
1721 for (n = 0; n < dd->nnodes; n++)
1725 if (ma->nat[n] > nalloc)
1727 nalloc = over_alloc_dd(ma->nat[n]);
1728 srenew(buf, nalloc);
1730 /* Use lv as a temporary buffer */
1732 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1734 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1736 copy_rvec(v[c], buf[a++]);
1739 if (a != ma->nat[n])
1741 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1746 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1747 DDRANK(dd, n), n, dd->mpi_comm_all);
1752 n = DDMASTERRANK(dd);
1754 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1756 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1758 copy_rvec(v[c], lv[a++]);
1765 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1766 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1771 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1774 gmx_domdec_master_t *ma;
1775 int *scounts = NULL, *disps = NULL;
1776 int n, i, c, a, nalloc = 0;
1783 get_commbuffer_counts(dd, &scounts, &disps);
1787 for (n = 0; n < dd->nnodes; n++)
1789 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1791 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1793 copy_rvec(v[c], buf[a++]);
1799 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1802 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1804 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1806 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1810 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1814 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1817 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1818 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1819 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1821 if (dfhist->nlambda > 0)
1823 int nlam = dfhist->nlambda;
1824 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1825 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1826 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1827 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1828 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1829 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1831 for (i = 0; i<nlam; i++) {
1832 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1833 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1834 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1835 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1836 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1837 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1842 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1843 t_state *state, t_state *state_local,
1848 nh = state->nhchainlength;
1852 for (i = 0; i < efptNR; i++)
1854 state_local->lambda[i] = state->lambda[i];
1856 state_local->fep_state = state->fep_state;
1857 state_local->veta = state->veta;
1858 state_local->vol0 = state->vol0;
1859 copy_mat(state->box, state_local->box);
1860 copy_mat(state->box_rel, state_local->box_rel);
1861 copy_mat(state->boxv, state_local->boxv);
1862 copy_mat(state->svir_prev, state_local->svir_prev);
1863 copy_mat(state->fvir_prev, state_local->fvir_prev);
1864 copy_df_history(&state_local->dfhist,&state->dfhist);
1865 for (i = 0; i < state_local->ngtc; i++)
1867 for (j = 0; j < nh; j++)
1869 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1870 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1872 state_local->therm_integral[i] = state->therm_integral[i];
1874 for (i = 0; i < state_local->nnhpres; i++)
1876 for (j = 0; j < nh; j++)
1878 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1879 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1883 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1884 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1885 dd_bcast(dd, sizeof(real), &state_local->veta);
1886 dd_bcast(dd, sizeof(real), &state_local->vol0);
1887 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1888 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1889 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1890 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1891 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1892 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1893 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1894 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1895 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1896 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1898 /* communicate df_history -- required for restarting from checkpoint */
1899 dd_distribute_dfhist(dd,&state_local->dfhist);
1901 if (dd->nat_home > state_local->nalloc)
1903 dd_realloc_state(state_local, f, dd->nat_home);
1905 for (i = 0; i < estNR; i++)
1907 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1912 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1915 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1918 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1921 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1924 if (state->nrngi == 1)
1927 state_local->nrng*sizeof(state_local->ld_rng[0]),
1928 state->ld_rng, state_local->ld_rng);
1933 state_local->nrng*sizeof(state_local->ld_rng[0]),
1934 state->ld_rng, state_local->ld_rng);
1938 if (state->nrngi == 1)
1940 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1941 state->ld_rngi, state_local->ld_rngi);
1945 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1946 state->ld_rngi, state_local->ld_rngi);
1949 case estDISRE_INITF:
1950 case estDISRE_RM3TAV:
1951 case estORIRE_INITF:
1953 /* Not implemented yet */
1956 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1962 static char dim2char(int dim)
1968 case XX: c = 'X'; break;
1969 case YY: c = 'Y'; break;
1970 case ZZ: c = 'Z'; break;
1971 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1977 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1978 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1980 rvec grid_s[2], *grid_r = NULL, cx, r;
1981 char fname[STRLEN], format[STRLEN], buf[22];
1983 int a, i, d, z, y, x;
1987 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1988 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1992 snew(grid_r, 2*dd->nnodes);
1995 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1999 for (d = 0; d < DIM; d++)
2001 for (i = 0; i < DIM; i++)
2009 if (d < ddbox->npbcdim && dd->nc[d] > 1)
2011 tric[d][i] = box[i][d]/box[i][i];
2020 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2021 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2022 out = gmx_fio_fopen(fname, "w");
2023 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2025 for (i = 0; i < dd->nnodes; i++)
2027 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2028 for (d = 0; d < DIM; d++)
2030 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2032 for (z = 0; z < 2; z++)
2034 for (y = 0; y < 2; y++)
2036 for (x = 0; x < 2; x++)
2038 cx[XX] = grid_r[i*2+x][XX];
2039 cx[YY] = grid_r[i*2+y][YY];
2040 cx[ZZ] = grid_r[i*2+z][ZZ];
2042 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
2043 ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
2047 for (d = 0; d < DIM; d++)
2049 for (x = 0; x < 4; x++)
2053 case 0: y = 1 + i*8 + 2*x; break;
2054 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2055 case 2: y = 1 + i*8 + x; break;
2057 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2061 gmx_fio_fclose(out);
2066 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2067 gmx_mtop_t *mtop, t_commrec *cr,
2068 int natoms, rvec x[], matrix box)
2070 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2072 int i, ii, resnr, c;
2073 char *atomname, *resname;
2080 natoms = dd->comm->nat[ddnatVSITE];
2083 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2085 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2086 sprintf(format4, "%s%s\n", pdbformat4, "%6.2f%6.2f");
2088 out = gmx_fio_fopen(fname, "w");
2090 fprintf(out, "TITLE %s\n", title);
2091 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2092 for (i = 0; i < natoms; i++)
2094 ii = dd->gatindex[i];
2095 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2096 if (i < dd->comm->nat[ddnatZONE])
2099 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2105 else if (i < dd->comm->nat[ddnatVSITE])
2107 b = dd->comm->zones.n;
2111 b = dd->comm->zones.n + 1;
2113 fprintf(out, strlen(atomname) < 4 ? format : format4,
2114 "ATOM", (ii+1)%100000,
2115 atomname, resname, ' ', resnr%10000, ' ',
2116 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2118 fprintf(out, "TER\n");
2120 gmx_fio_fclose(out);
2123 real dd_cutoff_mbody(gmx_domdec_t *dd)
2125 gmx_domdec_comm_t *comm;
2132 if (comm->bInterCGBondeds)
2134 if (comm->cutoff_mbody > 0)
2136 r = comm->cutoff_mbody;
2140 /* cutoff_mbody=0 means we do not have DLB */
2141 r = comm->cellsize_min[dd->dim[0]];
2142 for (di = 1; di < dd->ndim; di++)
2144 r = min(r, comm->cellsize_min[dd->dim[di]]);
2146 if (comm->bBondComm)
2148 r = max(r, comm->cutoff_mbody);
2152 r = min(r, comm->cutoff);
2160 real dd_cutoff_twobody(gmx_domdec_t *dd)
2164 r_mb = dd_cutoff_mbody(dd);
2166 return max(dd->comm->cutoff, r_mb);
2170 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2174 nc = dd->nc[dd->comm->cartpmedim];
2175 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2176 copy_ivec(coord, coord_pme);
2177 coord_pme[dd->comm->cartpmedim] =
2178 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2181 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2183 /* Here we assign a PME node to communicate with this DD node
2184 * by assuming that the major index of both is x.
2185 * We add cr->npmenodes/2 to obtain an even distribution.
2187 return (ddindex*npme + npme/2)/ndd;
2190 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2192 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2195 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2197 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2200 static int *dd_pmenodes(t_commrec *cr)
2205 snew(pmenodes, cr->npmenodes);
2207 for (i = 0; i < cr->dd->nnodes; i++)
2209 p0 = cr_ddindex2pmeindex(cr, i);
2210 p1 = cr_ddindex2pmeindex(cr, i+1);
2211 if (i+1 == cr->dd->nnodes || p1 > p0)
2215 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2217 pmenodes[n] = i + 1 + n;
2225 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2228 ivec coords, coords_pme, nc;
2233 if (dd->comm->bCartesian) {
2234 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2235 dd_coords2pmecoords(dd,coords,coords_pme);
2236 copy_ivec(dd->ntot,nc);
2237 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2238 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2240 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2242 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2248 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2253 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2255 gmx_domdec_comm_t *comm;
2257 int ddindex, nodeid = -1;
2259 comm = cr->dd->comm;
2264 if (comm->bCartesianPP_PME)
2267 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2272 ddindex = dd_index(cr->dd->nc, coords);
2273 if (comm->bCartesianPP)
2275 nodeid = comm->ddindex2simnodeid[ddindex];
2281 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2293 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2296 gmx_domdec_comm_t *comm;
2297 ivec coord, coord_pme;
2304 /* This assumes a uniform x domain decomposition grid cell size */
2305 if (comm->bCartesianPP_PME)
2308 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2309 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2311 /* This is a PP node */
2312 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2313 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2317 else if (comm->bCartesianPP)
2319 if (sim_nodeid < dd->nnodes)
2321 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2326 /* This assumes DD cells with identical x coordinates
2327 * are numbered sequentially.
2329 if (dd->comm->pmenodes == NULL)
2331 if (sim_nodeid < dd->nnodes)
2333 /* The DD index equals the nodeid */
2334 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2340 while (sim_nodeid > dd->comm->pmenodes[i])
2344 if (sim_nodeid < dd->comm->pmenodes[i])
2346 pmenode = dd->comm->pmenodes[i];
2354 void get_pme_nnodes(const gmx_domdec_t *dd,
2355 int *npmenodes_x, int *npmenodes_y)
2359 *npmenodes_x = dd->comm->npmenodes_x;
2360 *npmenodes_y = dd->comm->npmenodes_y;
2369 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2371 gmx_bool bPMEOnlyNode;
2373 if (DOMAINDECOMP(cr))
2375 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2379 bPMEOnlyNode = FALSE;
2382 return bPMEOnlyNode;
2385 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2386 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2390 ivec coord, coord_pme;
2394 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2397 for (x = 0; x < dd->nc[XX]; x++)
2399 for (y = 0; y < dd->nc[YY]; y++)
2401 for (z = 0; z < dd->nc[ZZ]; z++)
2403 if (dd->comm->bCartesianPP_PME)
2408 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2409 if (dd->ci[XX] == coord_pme[XX] &&
2410 dd->ci[YY] == coord_pme[YY] &&
2411 dd->ci[ZZ] == coord_pme[ZZ])
2413 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2418 /* The slab corresponds to the nodeid in the PME group */
2419 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2421 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2428 /* The last PP-only node is the peer node */
2429 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2433 fprintf(debug, "Receive coordinates from PP nodes:");
2434 for (x = 0; x < *nmy_ddnodes; x++)
2436 fprintf(debug, " %d", (*my_ddnodes)[x]);
2438 fprintf(debug, "\n");
2442 static gmx_bool receive_vir_ener(t_commrec *cr)
2444 gmx_domdec_comm_t *comm;
2445 int pmenode, coords[DIM], rank;
2449 if (cr->npmenodes < cr->dd->nnodes)
2451 comm = cr->dd->comm;
2452 if (comm->bCartesianPP_PME)
2454 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2456 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2457 coords[comm->cartpmedim]++;
2458 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2460 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2461 if (dd_simnode2pmenode(cr, rank) == pmenode)
2463 /* This is not the last PP node for pmenode */
2471 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2472 if (cr->sim_nodeid+1 < cr->nnodes &&
2473 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2475 /* This is not the last PP node for pmenode */
2484 static void set_zones_ncg_home(gmx_domdec_t *dd)
2486 gmx_domdec_zones_t *zones;
2489 zones = &dd->comm->zones;
2491 zones->cg_range[0] = 0;
2492 for (i = 1; i < zones->n+1; i++)
2494 zones->cg_range[i] = dd->ncg_home;
2496 /* zone_ncg1[0] should always be equal to ncg_home */
2497 dd->comm->zone_ncg1[0] = dd->ncg_home;
2500 static void rebuild_cgindex(gmx_domdec_t *dd,
2501 const int *gcgs_index, t_state *state)
2503 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2506 dd_cg_gl = dd->index_gl;
2507 cgindex = dd->cgindex;
2510 for (i = 0; i < state->ncg_gl; i++)
2514 dd_cg_gl[i] = cg_gl;
2515 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2519 dd->ncg_home = state->ncg_gl;
2522 set_zones_ncg_home(dd);
2525 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2527 while (cg >= cginfo_mb->cg_end)
2532 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2535 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2536 t_forcerec *fr, char *bLocalCG)
2538 cginfo_mb_t *cginfo_mb;
2544 cginfo_mb = fr->cginfo_mb;
2545 cginfo = fr->cginfo;
2547 for (cg = cg0; cg < cg1; cg++)
2549 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2553 if (bLocalCG != NULL)
2555 for (cg = cg0; cg < cg1; cg++)
2557 bLocalCG[index_gl[cg]] = TRUE;
2562 static void make_dd_indices(gmx_domdec_t *dd,
2563 const int *gcgs_index, int cg_start)
2565 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2566 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2571 bLocalCG = dd->comm->bLocalCG;
2573 if (dd->nat_tot > dd->gatindex_nalloc)
2575 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2576 srenew(dd->gatindex, dd->gatindex_nalloc);
2579 nzone = dd->comm->zones.n;
2580 zone2cg = dd->comm->zones.cg_range;
2581 zone_ncg1 = dd->comm->zone_ncg1;
2582 index_gl = dd->index_gl;
2583 gatindex = dd->gatindex;
2584 bCGs = dd->comm->bCGs;
2586 if (zone2cg[1] != dd->ncg_home)
2588 gmx_incons("dd->ncg_zone is not up to date");
2591 /* Make the local to global and global to local atom index */
2592 a = dd->cgindex[cg_start];
2593 for (zone = 0; zone < nzone; zone++)
2601 cg0 = zone2cg[zone];
2603 cg1 = zone2cg[zone+1];
2604 cg1_p1 = cg0 + zone_ncg1[zone];
2606 for (cg = cg0; cg < cg1; cg++)
2611 /* Signal that this cg is from more than one pulse away */
2614 cg_gl = index_gl[cg];
2617 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2620 ga2la_set(dd->ga2la, a_gl, a, zone1);
2626 gatindex[a] = cg_gl;
2627 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2634 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2637 int ncg, i, ngl, nerr;
2640 if (bLocalCG == NULL)
2644 for (i = 0; i < dd->ncg_tot; i++)
2646 if (!bLocalCG[dd->index_gl[i]])
2649 "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);
2654 for (i = 0; i < ncg_sys; i++)
2661 if (ngl != dd->ncg_tot)
2663 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);
2670 static void check_index_consistency(gmx_domdec_t *dd,
2671 int natoms_sys, int ncg_sys,
2674 int nerr, ngl, i, a, cell;
2679 if (dd->comm->DD_debug > 1)
2681 snew(have, natoms_sys);
2682 for (a = 0; a < dd->nat_tot; a++)
2684 if (have[dd->gatindex[a]] > 0)
2686 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);
2690 have[dd->gatindex[a]] = a + 1;
2696 snew(have, dd->nat_tot);
2699 for (i = 0; i < natoms_sys; i++)
2701 if (ga2la_get(dd->ga2la, i, &a, &cell))
2703 if (a >= dd->nat_tot)
2705 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);
2711 if (dd->gatindex[a] != i)
2713 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);
2720 if (ngl != dd->nat_tot)
2723 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2724 dd->rank, where, ngl, dd->nat_tot);
2726 for (a = 0; a < dd->nat_tot; a++)
2731 "DD node %d, %s: local atom %d, global %d has no global index\n",
2732 dd->rank, where, a+1, dd->gatindex[a]+1);
2737 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2741 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2742 dd->rank, where, nerr);
2746 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2753 /* Clear the whole list without searching */
2754 ga2la_clear(dd->ga2la);
2758 for (i = a_start; i < dd->nat_tot; i++)
2760 ga2la_del(dd->ga2la, dd->gatindex[i]);
2764 bLocalCG = dd->comm->bLocalCG;
2767 for (i = cg_start; i < dd->ncg_tot; i++)
2769 bLocalCG[dd->index_gl[i]] = FALSE;
2773 dd_clear_local_vsite_indices(dd);
2775 if (dd->constraints)
2777 dd_clear_local_constraint_indices(dd);
2781 /* This function should be used for moving the domain boudaries during DLB,
2782 * for obtaining the minimum cell size. It checks the initially set limit
2783 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2784 * and, possibly, a longer cut-off limit set for PME load balancing.
2786 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2790 cellsize_min = comm->cellsize_min[dim];
2792 if (!comm->bVacDLBNoLimit)
2794 /* The cut-off might have changed, e.g. by PME load balacning,
2795 * from the value used to set comm->cellsize_min, so check it.
2797 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2799 if (comm->bPMELoadBalDLBLimits)
2801 /* Check for the cut-off limit set by the PME load balancing */
2802 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2806 return cellsize_min;
2809 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2812 real grid_jump_limit;
2814 /* The distance between the boundaries of cells at distance
2815 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2816 * and by the fact that cells should not be shifted by more than
2817 * half their size, such that cg's only shift by one cell
2818 * at redecomposition.
2820 grid_jump_limit = comm->cellsize_limit;
2821 if (!comm->bVacDLBNoLimit)
2823 if (comm->bPMELoadBalDLBLimits)
2825 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2827 grid_jump_limit = max(grid_jump_limit,
2828 cutoff/comm->cd[dim_ind].np);
2831 return grid_jump_limit;
2834 static gmx_bool check_grid_jump(gmx_large_int_t step,
2840 gmx_domdec_comm_t *comm;
2849 for (d = 1; d < dd->ndim; d++)
2852 limit = grid_jump_limit(comm, cutoff, d);
2853 bfac = ddbox->box_size[dim];
2854 if (ddbox->tric_dir[dim])
2856 bfac *= ddbox->skew_fac[dim];
2858 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2859 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2867 /* This error should never be triggered under normal
2868 * circumstances, but you never know ...
2870 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.",
2871 gmx_step_str(step, buf),
2872 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2880 static int dd_load_count(gmx_domdec_comm_t *comm)
2882 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2885 static float dd_force_load(gmx_domdec_comm_t *comm)
2892 if (comm->eFlop > 1)
2894 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2899 load = comm->cycl[ddCyclF];
2900 if (comm->cycl_n[ddCyclF] > 1)
2902 /* Subtract the maximum of the last n cycle counts
2903 * to get rid of possible high counts due to other sources,
2904 * for instance system activity, that would otherwise
2905 * affect the dynamic load balancing.
2907 load -= comm->cycl_max[ddCyclF];
2911 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2913 float gpu_wait, gpu_wait_sum;
2915 gpu_wait = comm->cycl[ddCyclWaitGPU];
2916 if (comm->cycl_n[ddCyclF] > 1)
2918 /* We should remove the WaitGPU time of the same MD step
2919 * as the one with the maximum F time, since the F time
2920 * and the wait time are not independent.
2921 * Furthermore, the step for the max F time should be chosen
2922 * the same on all ranks that share the same GPU.
2923 * But to keep the code simple, we remove the average instead.
2924 * The main reason for artificially long times at some steps
2925 * is spurious CPU activity or MPI time, so we don't expect
2926 * that changes in the GPU wait time matter a lot here.
2928 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2930 /* Sum the wait times over the ranks that share the same GPU */
2931 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2932 comm->mpi_comm_gpu_shared);
2933 /* Replace the wait time by the average over the ranks */
2934 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2942 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2944 gmx_domdec_comm_t *comm;
2949 snew(*dim_f, dd->nc[dim]+1);
2951 for (i = 1; i < dd->nc[dim]; i++)
2953 if (comm->slb_frac[dim])
2955 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2959 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2962 (*dim_f)[dd->nc[dim]] = 1;
2965 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2967 int pmeindex, slab, nso, i;
2970 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2976 ddpme->dim = dimind;
2978 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2980 ddpme->nslab = (ddpme->dim == 0 ?
2981 dd->comm->npmenodes_x :
2982 dd->comm->npmenodes_y);
2984 if (ddpme->nslab <= 1)
2989 nso = dd->comm->npmenodes/ddpme->nslab;
2990 /* Determine for each PME slab the PP location range for dimension dim */
2991 snew(ddpme->pp_min, ddpme->nslab);
2992 snew(ddpme->pp_max, ddpme->nslab);
2993 for (slab = 0; slab < ddpme->nslab; slab++)
2995 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2996 ddpme->pp_max[slab] = 0;
2998 for (i = 0; i < dd->nnodes; i++)
3000 ddindex2xyz(dd->nc, i, xyz);
3001 /* For y only use our y/z slab.
3002 * This assumes that the PME x grid size matches the DD grid size.
3004 if (dimind == 0 || xyz[XX] == dd->ci[XX])
3006 pmeindex = ddindex2pmeindex(dd, i);
3009 slab = pmeindex/nso;
3013 slab = pmeindex % ddpme->nslab;
3015 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
3016 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
3020 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
3023 int dd_pme_maxshift_x(gmx_domdec_t *dd)
3025 if (dd->comm->ddpme[0].dim == XX)
3027 return dd->comm->ddpme[0].maxshift;
3035 int dd_pme_maxshift_y(gmx_domdec_t *dd)
3037 if (dd->comm->ddpme[0].dim == YY)
3039 return dd->comm->ddpme[0].maxshift;
3041 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3043 return dd->comm->ddpme[1].maxshift;
3051 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3052 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3054 gmx_domdec_comm_t *comm;
3057 real range, pme_boundary;
3061 nc = dd->nc[ddpme->dim];
3064 if (!ddpme->dim_match)
3066 /* PP decomposition is not along dim: the worst situation */
3069 else if (ns <= 3 || (bUniform && ns == nc))
3071 /* The optimal situation */
3076 /* We need to check for all pme nodes which nodes they
3077 * could possibly need to communicate with.
3079 xmin = ddpme->pp_min;
3080 xmax = ddpme->pp_max;
3081 /* Allow for atoms to be maximally 2/3 times the cut-off
3082 * out of their DD cell. This is a reasonable balance between
3083 * between performance and support for most charge-group/cut-off
3086 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3087 /* Avoid extra communication when we are exactly at a boundary */
3091 for (s = 0; s < ns; s++)
3093 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3094 pme_boundary = (real)s/ns;
3097 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3099 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3103 pme_boundary = (real)(s+1)/ns;
3106 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3108 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3115 ddpme->maxshift = sh;
3119 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3120 ddpme->dim, ddpme->maxshift);
3124 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3128 for (d = 0; d < dd->ndim; d++)
3131 if (dim < ddbox->nboundeddim &&
3132 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3133 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3135 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",
3136 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3137 dd->nc[dim], dd->comm->cellsize_limit);
3142 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3143 gmx_bool bMaster, ivec npulse)
3145 gmx_domdec_comm_t *comm;
3148 real *cell_x, cell_dx, cellsize;
3152 for (d = 0; d < DIM; d++)
3154 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3156 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3159 cell_dx = ddbox->box_size[d]/dd->nc[d];
3162 for (j = 0; j < dd->nc[d]+1; j++)
3164 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3169 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3170 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3172 cellsize = cell_dx*ddbox->skew_fac[d];
3173 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3177 cellsize_min[d] = cellsize;
3181 /* Statically load balanced grid */
3182 /* Also when we are not doing a master distribution we determine
3183 * all cell borders in a loop to obtain identical values
3184 * to the master distribution case and to determine npulse.
3188 cell_x = dd->ma->cell_x[d];
3192 snew(cell_x, dd->nc[d]+1);
3194 cell_x[0] = ddbox->box0[d];
3195 for (j = 0; j < dd->nc[d]; j++)
3197 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3198 cell_x[j+1] = cell_x[j] + cell_dx;
3199 cellsize = cell_dx*ddbox->skew_fac[d];
3200 while (cellsize*npulse[d] < comm->cutoff &&
3201 npulse[d] < dd->nc[d]-1)
3205 cellsize_min[d] = min(cellsize_min[d], cellsize);
3209 comm->cell_x0[d] = cell_x[dd->ci[d]];
3210 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3214 /* The following limitation is to avoid that a cell would receive
3215 * some of its own home charge groups back over the periodic boundary.
3216 * Double charge groups cause trouble with the global indices.
3218 if (d < ddbox->npbcdim &&
3219 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3221 gmx_fatal_collective(FARGS, NULL, dd,
3222 "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",
3223 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3225 dd->nc[d], dd->nc[d],
3226 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3230 if (!comm->bDynLoadBal)
3232 copy_rvec(cellsize_min, comm->cellsize_min);
3235 for (d = 0; d < comm->npmedecompdim; d++)
3237 set_pme_maxshift(dd, &comm->ddpme[d],
3238 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3239 comm->ddpme[d].slb_dim_f);
3244 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3245 int d, int dim, gmx_domdec_root_t *root,
3247 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3249 gmx_domdec_comm_t *comm;
3250 int ncd, i, j, nmin, nmin_old;
3251 gmx_bool bLimLo, bLimHi;
3253 real fac, halfway, cellsize_limit_f_i, region_size;
3254 gmx_bool bPBC, bLastHi = FALSE;
3255 int nrange[] = {range[0], range[1]};
3257 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3263 bPBC = (dim < ddbox->npbcdim);
3265 cell_size = root->buf_ncd;
3269 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3272 /* First we need to check if the scaling does not make cells
3273 * smaller than the smallest allowed size.
3274 * We need to do this iteratively, since if a cell is too small,
3275 * it needs to be enlarged, which makes all the other cells smaller,
3276 * which could in turn make another cell smaller than allowed.
3278 for (i = range[0]; i < range[1]; i++)
3280 root->bCellMin[i] = FALSE;
3286 /* We need the total for normalization */
3288 for (i = range[0]; i < range[1]; i++)
3290 if (root->bCellMin[i] == FALSE)
3292 fac += cell_size[i];
3295 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3296 /* Determine the cell boundaries */
3297 for (i = range[0]; i < range[1]; i++)
3299 if (root->bCellMin[i] == FALSE)
3301 cell_size[i] *= fac;
3302 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3304 cellsize_limit_f_i = 0;
3308 cellsize_limit_f_i = cellsize_limit_f;
3310 if (cell_size[i] < cellsize_limit_f_i)
3312 root->bCellMin[i] = TRUE;
3313 cell_size[i] = cellsize_limit_f_i;
3317 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3320 while (nmin > nmin_old);
3323 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3324 /* For this check we should not use DD_CELL_MARGIN,
3325 * but a slightly smaller factor,
3326 * since rounding could get use below the limit.
3328 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3331 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",
3332 gmx_step_str(step, buf),
3333 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3334 ncd, comm->cellsize_min[dim]);
3337 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3341 /* Check if the boundary did not displace more than halfway
3342 * each of the cells it bounds, as this could cause problems,
3343 * especially when the differences between cell sizes are large.
3344 * If changes are applied, they will not make cells smaller
3345 * than the cut-off, as we check all the boundaries which
3346 * might be affected by a change and if the old state was ok,
3347 * the cells will at most be shrunk back to their old size.
3349 for (i = range[0]+1; i < range[1]; i++)
3351 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3352 if (root->cell_f[i] < halfway)
3354 root->cell_f[i] = halfway;
3355 /* Check if the change also causes shifts of the next boundaries */
3356 for (j = i+1; j < range[1]; j++)
3358 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3360 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3364 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3365 if (root->cell_f[i] > halfway)
3367 root->cell_f[i] = halfway;
3368 /* Check if the change also causes shifts of the next boundaries */
3369 for (j = i-1; j >= range[0]+1; j--)
3371 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3373 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3380 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3381 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3382 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3383 * for a and b nrange is used */
3386 /* Take care of the staggering of the cell boundaries */
3389 for (i = range[0]; i < range[1]; i++)
3391 root->cell_f_max0[i] = root->cell_f[i];
3392 root->cell_f_min1[i] = root->cell_f[i+1];
3397 for (i = range[0]+1; i < range[1]; i++)
3399 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3400 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3401 if (bLimLo && bLimHi)
3403 /* Both limits violated, try the best we can */
3404 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3405 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3406 nrange[0] = range[0];
3408 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3411 nrange[1] = range[1];
3412 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3418 /* root->cell_f[i] = root->bound_min[i]; */
3419 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3422 else if (bLimHi && !bLastHi)
3425 if (nrange[1] < range[1]) /* found a LimLo before */
3427 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3428 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3429 nrange[0] = nrange[1];
3431 root->cell_f[i] = root->bound_max[i];
3433 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3435 nrange[1] = range[1];
3438 if (nrange[1] < range[1]) /* found last a LimLo */
3440 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3441 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3442 nrange[0] = nrange[1];
3443 nrange[1] = range[1];
3444 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3446 else if (nrange[0] > range[0]) /* found at least one LimHi */
3448 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3455 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3456 int d, int dim, gmx_domdec_root_t *root,
3457 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3458 gmx_bool bUniform, gmx_large_int_t step)
3460 gmx_domdec_comm_t *comm;
3461 int ncd, d1, i, j, pos;
3463 real load_aver, load_i, imbalance, change, change_max, sc;
3464 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3468 int range[] = { 0, 0 };
3472 /* Convert the maximum change from the input percentage to a fraction */
3473 change_limit = comm->dlb_scale_lim*0.01;
3477 bPBC = (dim < ddbox->npbcdim);
3479 cell_size = root->buf_ncd;
3481 /* Store the original boundaries */
3482 for (i = 0; i < ncd+1; i++)
3484 root->old_cell_f[i] = root->cell_f[i];
3488 for (i = 0; i < ncd; i++)
3490 cell_size[i] = 1.0/ncd;
3493 else if (dd_load_count(comm))
3495 load_aver = comm->load[d].sum_m/ncd;
3497 for (i = 0; i < ncd; i++)
3499 /* Determine the relative imbalance of cell i */
3500 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3501 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3502 /* Determine the change of the cell size using underrelaxation */
3503 change = -relax*imbalance;
3504 change_max = max(change_max, max(change, -change));
3506 /* Limit the amount of scaling.
3507 * We need to use the same rescaling for all cells in one row,
3508 * otherwise the load balancing might not converge.
3511 if (change_max > change_limit)
3513 sc *= change_limit/change_max;
3515 for (i = 0; i < ncd; i++)
3517 /* Determine the relative imbalance of cell i */
3518 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3519 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3520 /* Determine the change of the cell size using underrelaxation */
3521 change = -sc*imbalance;
3522 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3526 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3527 cellsize_limit_f *= DD_CELL_MARGIN;
3528 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3529 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3530 if (ddbox->tric_dir[dim])
3532 cellsize_limit_f /= ddbox->skew_fac[dim];
3533 dist_min_f /= ddbox->skew_fac[dim];
3535 if (bDynamicBox && d > 0)
3537 dist_min_f *= DD_PRES_SCALE_MARGIN;
3539 if (d > 0 && !bUniform)
3541 /* Make sure that the grid is not shifted too much */
3542 for (i = 1; i < ncd; i++)
3544 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3546 gmx_incons("Inconsistent DD boundary staggering limits!");
3548 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3549 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3552 root->bound_min[i] += 0.5*space;
3554 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3555 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3558 root->bound_max[i] += 0.5*space;
3563 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3565 root->cell_f_max0[i-1] + dist_min_f,
3566 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3567 root->cell_f_min1[i] - dist_min_f);
3572 root->cell_f[0] = 0;
3573 root->cell_f[ncd] = 1;
3574 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3577 /* After the checks above, the cells should obey the cut-off
3578 * restrictions, but it does not hurt to check.
3580 for (i = 0; i < ncd; i++)
3584 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3585 dim, i, root->cell_f[i], root->cell_f[i+1]);
3588 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3589 root->cell_f[i+1] - root->cell_f[i] <
3590 cellsize_limit_f/DD_CELL_MARGIN)
3594 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3595 gmx_step_str(step, buf), dim2char(dim), i,
3596 (root->cell_f[i+1] - root->cell_f[i])
3597 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3602 /* Store the cell boundaries of the lower dimensions at the end */
3603 for (d1 = 0; d1 < d; d1++)
3605 root->cell_f[pos++] = comm->cell_f0[d1];
3606 root->cell_f[pos++] = comm->cell_f1[d1];
3609 if (d < comm->npmedecompdim)
3611 /* The master determines the maximum shift for
3612 * the coordinate communication between separate PME nodes.
3614 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3616 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3619 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3623 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3624 gmx_ddbox_t *ddbox, int dimind)
3626 gmx_domdec_comm_t *comm;
3631 /* Set the cell dimensions */
3632 dim = dd->dim[dimind];
3633 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3634 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3635 if (dim >= ddbox->nboundeddim)
3637 comm->cell_x0[dim] += ddbox->box0[dim];
3638 comm->cell_x1[dim] += ddbox->box0[dim];
3642 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3643 int d, int dim, real *cell_f_row,
3646 gmx_domdec_comm_t *comm;
3652 /* Each node would only need to know two fractions,
3653 * but it is probably cheaper to broadcast the whole array.
3655 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3656 0, comm->mpi_comm_load[d]);
3658 /* Copy the fractions for this dimension from the buffer */
3659 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3660 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3661 /* The whole array was communicated, so set the buffer position */
3662 pos = dd->nc[dim] + 1;
3663 for (d1 = 0; d1 <= d; d1++)
3667 /* Copy the cell fractions of the lower dimensions */
3668 comm->cell_f0[d1] = cell_f_row[pos++];
3669 comm->cell_f1[d1] = cell_f_row[pos++];
3671 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3673 /* Convert the communicated shift from float to int */
3674 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3677 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3681 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3682 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3683 gmx_bool bUniform, gmx_large_int_t step)
3685 gmx_domdec_comm_t *comm;
3687 gmx_bool bRowMember, bRowRoot;
3692 for (d = 0; d < dd->ndim; d++)
3697 for (d1 = d; d1 < dd->ndim; d1++)
3699 if (dd->ci[dd->dim[d1]] > 0)
3712 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3713 ddbox, bDynamicBox, bUniform, step);
3714 cell_f_row = comm->root[d]->cell_f;
3718 cell_f_row = comm->cell_f_row;
3720 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3725 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3729 /* This function assumes the box is static and should therefore
3730 * not be called when the box has changed since the last
3731 * call to dd_partition_system.
3733 for (d = 0; d < dd->ndim; d++)
3735 relative_to_absolute_cell_bounds(dd, ddbox, d);
3741 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3742 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3743 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3744 gmx_wallcycle_t wcycle)
3746 gmx_domdec_comm_t *comm;
3753 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3754 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3755 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3757 else if (bDynamicBox)
3759 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3762 /* Set the dimensions for which no DD is used */
3763 for (dim = 0; dim < DIM; dim++)
3765 if (dd->nc[dim] == 1)
3767 comm->cell_x0[dim] = 0;
3768 comm->cell_x1[dim] = ddbox->box_size[dim];
3769 if (dim >= ddbox->nboundeddim)
3771 comm->cell_x0[dim] += ddbox->box0[dim];
3772 comm->cell_x1[dim] += ddbox->box0[dim];
3778 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3781 gmx_domdec_comm_dim_t *cd;
3783 for (d = 0; d < dd->ndim; d++)
3785 cd = &dd->comm->cd[d];
3786 np = npulse[dd->dim[d]];
3787 if (np > cd->np_nalloc)
3791 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3792 dim2char(dd->dim[d]), np);
3794 if (DDMASTER(dd) && cd->np_nalloc > 0)
3796 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3798 srenew(cd->ind, np);
3799 for (i = cd->np_nalloc; i < np; i++)
3801 cd->ind[i].index = NULL;
3802 cd->ind[i].nalloc = 0;
3811 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3812 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3813 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3814 gmx_wallcycle_t wcycle)
3816 gmx_domdec_comm_t *comm;
3822 /* Copy the old cell boundaries for the cg displacement check */
3823 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3824 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3826 if (comm->bDynLoadBal)
3830 check_box_size(dd, ddbox);
3832 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3836 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3837 realloc_comm_ind(dd, npulse);
3842 for (d = 0; d < DIM; d++)
3844 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3845 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3850 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3852 rvec cell_ns_x0, rvec cell_ns_x1,
3853 gmx_large_int_t step)
3855 gmx_domdec_comm_t *comm;
3860 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3862 dim = dd->dim[dim_ind];
3864 /* Without PBC we don't have restrictions on the outer cells */
3865 if (!(dim >= ddbox->npbcdim &&
3866 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3867 comm->bDynLoadBal &&
3868 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3869 comm->cellsize_min[dim])
3872 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",
3873 gmx_step_str(step, buf), dim2char(dim),
3874 comm->cell_x1[dim] - comm->cell_x0[dim],
3875 ddbox->skew_fac[dim],
3876 dd->comm->cellsize_min[dim],
3877 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3881 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3883 /* Communicate the boundaries and update cell_ns_x0/1 */
3884 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3885 if (dd->bGridJump && dd->ndim > 1)
3887 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3892 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3896 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3904 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3905 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3914 static void check_screw_box(matrix box)
3916 /* Mathematical limitation */
3917 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3919 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3922 /* Limitation due to the asymmetry of the eighth shell method */
3923 if (box[ZZ][YY] != 0)
3925 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3929 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3930 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3933 gmx_domdec_master_t *ma;
3934 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3935 int i, icg, j, k, k0, k1, d, npbcdim;
3937 rvec box_size, cg_cm;
3939 real nrcg, inv_ncg, pos_d;
3941 gmx_bool bUnbounded, bScrew;
3945 if (tmp_ind == NULL)
3947 snew(tmp_nalloc, dd->nnodes);
3948 snew(tmp_ind, dd->nnodes);
3949 for (i = 0; i < dd->nnodes; i++)
3951 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3952 snew(tmp_ind[i], tmp_nalloc[i]);
3956 /* Clear the count */
3957 for (i = 0; i < dd->nnodes; i++)
3963 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3965 cgindex = cgs->index;
3967 /* Compute the center of geometry for all charge groups */
3968 for (icg = 0; icg < cgs->nr; icg++)
3971 k1 = cgindex[icg+1];
3975 copy_rvec(pos[k0], cg_cm);
3982 for (k = k0; (k < k1); k++)
3984 rvec_inc(cg_cm, pos[k]);
3986 for (d = 0; (d < DIM); d++)
3988 cg_cm[d] *= inv_ncg;
3991 /* Put the charge group in the box and determine the cell index */
3992 for (d = DIM-1; d >= 0; d--)
3995 if (d < dd->npbcdim)
3997 bScrew = (dd->bScrewPBC && d == XX);
3998 if (tric_dir[d] && dd->nc[d] > 1)
4000 /* Use triclinic coordintates for this dimension */
4001 for (j = d+1; j < DIM; j++)
4003 pos_d += cg_cm[j]*tcm[j][d];
4006 while (pos_d >= box[d][d])
4009 rvec_dec(cg_cm, box[d]);
4012 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4013 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4015 for (k = k0; (k < k1); k++)
4017 rvec_dec(pos[k], box[d]);
4020 pos[k][YY] = box[YY][YY] - pos[k][YY];
4021 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4028 rvec_inc(cg_cm, box[d]);
4031 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4032 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4034 for (k = k0; (k < k1); k++)
4036 rvec_inc(pos[k], box[d]);
4039 pos[k][YY] = box[YY][YY] - pos[k][YY];
4040 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4045 /* This could be done more efficiently */
4047 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4052 i = dd_index(dd->nc, ind);
4053 if (ma->ncg[i] == tmp_nalloc[i])
4055 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4056 srenew(tmp_ind[i], tmp_nalloc[i]);
4058 tmp_ind[i][ma->ncg[i]] = icg;
4060 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4064 for (i = 0; i < dd->nnodes; i++)
4067 for (k = 0; k < ma->ncg[i]; k++)
4069 ma->cg[k1++] = tmp_ind[i][k];
4072 ma->index[dd->nnodes] = k1;
4074 for (i = 0; i < dd->nnodes; i++)
4084 fprintf(fplog, "Charge group distribution at step %s:",
4085 gmx_step_str(step, buf));
4086 for (i = 0; i < dd->nnodes; i++)
4088 fprintf(fplog, " %d", ma->ncg[i]);
4090 fprintf(fplog, "\n");
4094 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4095 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4098 gmx_domdec_master_t *ma = NULL;
4101 int *ibuf, buf2[2] = { 0, 0 };
4102 gmx_bool bMaster = DDMASTER(dd);
4109 check_screw_box(box);
4112 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4114 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4115 for (i = 0; i < dd->nnodes; i++)
4117 ma->ibuf[2*i] = ma->ncg[i];
4118 ma->ibuf[2*i+1] = ma->nat[i];
4126 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4128 dd->ncg_home = buf2[0];
4129 dd->nat_home = buf2[1];
4130 dd->ncg_tot = dd->ncg_home;
4131 dd->nat_tot = dd->nat_home;
4132 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4134 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4135 srenew(dd->index_gl, dd->cg_nalloc);
4136 srenew(dd->cgindex, dd->cg_nalloc+1);
4140 for (i = 0; i < dd->nnodes; i++)
4142 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4143 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4148 DDMASTER(dd) ? ma->ibuf : NULL,
4149 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4150 DDMASTER(dd) ? ma->cg : NULL,
4151 dd->ncg_home*sizeof(int), dd->index_gl);
4153 /* Determine the home charge group sizes */
4155 for (i = 0; i < dd->ncg_home; i++)
4157 cg_gl = dd->index_gl[i];
4159 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4164 fprintf(debug, "Home charge groups:\n");
4165 for (i = 0; i < dd->ncg_home; i++)
4167 fprintf(debug, " %d", dd->index_gl[i]);
4170 fprintf(debug, "\n");
4173 fprintf(debug, "\n");
4177 static int compact_and_copy_vec_at(int ncg, int *move,
4180 rvec *src, gmx_domdec_comm_t *comm,
4183 int m, icg, i, 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 for (i = i0; i < i1; i++)
4206 copy_rvec(src[i], src[home_pos++]);
4212 /* Copy to the communication buffer */
4214 pos_vec[m] += 1 + vec*nrcg;
4215 for (i = i0; i < i1; i++)
4217 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4219 pos_vec[m] += (nvec - vec - 1)*nrcg;
4223 home_pos += i1 - i0;
4231 static int compact_and_copy_vec_cg(int ncg, int *move,
4233 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4236 int m, icg, i0, i1, nrcg;
4242 for (m = 0; m < DIM*2; m++)
4248 for (icg = 0; icg < ncg; icg++)
4250 i1 = cgindex[icg+1];
4256 /* Compact the home array in place */
4257 copy_rvec(src[icg], src[home_pos++]);
4263 /* Copy to the communication buffer */
4264 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4265 pos_vec[m] += 1 + nrcg*nvec;
4277 static int compact_ind(int ncg, int *move,
4278 int *index_gl, int *cgindex,
4280 gmx_ga2la_t ga2la, char *bLocalCG,
4283 int cg, nat, a0, a1, a, a_gl;
4288 for (cg = 0; cg < ncg; cg++)
4294 /* Compact the home arrays in place.
4295 * Anything that can be done here avoids access to global arrays.
4297 cgindex[home_pos] = nat;
4298 for (a = a0; a < a1; a++)
4301 gatindex[nat] = a_gl;
4302 /* The cell number stays 0, so we don't need to set it */
4303 ga2la_change_la(ga2la, a_gl, nat);
4306 index_gl[home_pos] = index_gl[cg];
4307 cginfo[home_pos] = cginfo[cg];
4308 /* The charge group remains local, so bLocalCG does not change */
4313 /* Clear the global indices */
4314 for (a = a0; a < a1; a++)
4316 ga2la_del(ga2la, gatindex[a]);
4320 bLocalCG[index_gl[cg]] = FALSE;
4324 cgindex[home_pos] = nat;
4329 static void clear_and_mark_ind(int ncg, int *move,
4330 int *index_gl, int *cgindex, int *gatindex,
4331 gmx_ga2la_t ga2la, char *bLocalCG,
4336 for (cg = 0; cg < ncg; cg++)
4342 /* Clear the global indices */
4343 for (a = a0; a < a1; a++)
4345 ga2la_del(ga2la, gatindex[a]);
4349 bLocalCG[index_gl[cg]] = FALSE;
4351 /* Signal that this cg has moved using the ns cell index.
4352 * Here we set it to -1. fill_grid will change it
4353 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4355 cell_index[cg] = -1;
4360 static void print_cg_move(FILE *fplog,
4362 gmx_large_int_t step, int cg, int dim, int dir,
4363 gmx_bool bHaveLimitdAndCMOld, real limitd,
4364 rvec cm_old, rvec cm_new, real pos_d)
4366 gmx_domdec_comm_t *comm;
4371 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4372 if (bHaveLimitdAndCMOld)
4374 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4375 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4379 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4380 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4382 fprintf(fplog, "distance out of cell %f\n",
4383 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4384 if (bHaveLimitdAndCMOld)
4386 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4387 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4389 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4390 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4391 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4393 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4394 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4396 comm->cell_x0[dim], comm->cell_x1[dim]);
4399 static void cg_move_error(FILE *fplog,
4401 gmx_large_int_t step, int cg, int dim, int dir,
4402 gmx_bool bHaveLimitdAndCMOld, real limitd,
4403 rvec cm_old, rvec cm_new, real pos_d)
4407 print_cg_move(fplog, dd, step, cg, dim, dir,
4408 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4410 print_cg_move(stderr, dd, step, cg, dim, dir,
4411 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4413 "A charge group moved too far between two domain decomposition steps\n"
4414 "This usually means that your system is not well equilibrated");
4417 static void rotate_state_atom(t_state *state, int a)
4421 for (est = 0; est < estNR; est++)
4423 if (EST_DISTR(est) && (state->flags & (1<<est)))
4428 /* Rotate the complete state; for a rectangular box only */
4429 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4430 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4433 state->v[a][YY] = -state->v[a][YY];
4434 state->v[a][ZZ] = -state->v[a][ZZ];
4437 state->sd_X[a][YY] = -state->sd_X[a][YY];
4438 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4441 state->cg_p[a][YY] = -state->cg_p[a][YY];
4442 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4444 case estDISRE_INITF:
4445 case estDISRE_RM3TAV:
4446 case estORIRE_INITF:
4448 /* These are distances, so not affected by rotation */
4451 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4457 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4459 if (natoms > comm->moved_nalloc)
4461 /* Contents should be preserved here */
4462 comm->moved_nalloc = over_alloc_dd(natoms);
4463 srenew(comm->moved, comm->moved_nalloc);
4469 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4472 ivec tric_dir, matrix tcm,
4473 rvec cell_x0, rvec cell_x1,
4474 rvec limitd, rvec limit0, rvec limit1,
4476 int cg_start, int cg_end,
4481 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4482 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4486 real inv_ncg, pos_d;
4489 npbcdim = dd->npbcdim;
4491 for (cg = cg_start; cg < cg_end; cg++)
4498 copy_rvec(state->x[k0], cm_new);
4505 for (k = k0; (k < k1); k++)
4507 rvec_inc(cm_new, state->x[k]);
4509 for (d = 0; (d < DIM); d++)
4511 cm_new[d] = inv_ncg*cm_new[d];
4516 /* Do pbc and check DD cell boundary crossings */
4517 for (d = DIM-1; d >= 0; d--)
4521 bScrew = (dd->bScrewPBC && d == XX);
4522 /* Determine the location of this cg in lattice coordinates */
4526 for (d2 = d+1; d2 < DIM; d2++)
4528 pos_d += cm_new[d2]*tcm[d2][d];
4531 /* Put the charge group in the triclinic unit-cell */
4532 if (pos_d >= cell_x1[d])
4534 if (pos_d >= limit1[d])
4536 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4537 cg_cm[cg], cm_new, pos_d);
4540 if (dd->ci[d] == dd->nc[d] - 1)
4542 rvec_dec(cm_new, state->box[d]);
4545 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4546 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4548 for (k = k0; (k < k1); k++)
4550 rvec_dec(state->x[k], state->box[d]);
4553 rotate_state_atom(state, k);
4558 else if (pos_d < cell_x0[d])
4560 if (pos_d < limit0[d])
4562 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4563 cg_cm[cg], cm_new, pos_d);
4568 rvec_inc(cm_new, state->box[d]);
4571 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4572 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4574 for (k = k0; (k < k1); k++)
4576 rvec_inc(state->x[k], state->box[d]);
4579 rotate_state_atom(state, k);
4585 else if (d < npbcdim)
4587 /* Put the charge group in the rectangular unit-cell */
4588 while (cm_new[d] >= state->box[d][d])
4590 rvec_dec(cm_new, state->box[d]);
4591 for (k = k0; (k < k1); k++)
4593 rvec_dec(state->x[k], state->box[d]);
4596 while (cm_new[d] < 0)
4598 rvec_inc(cm_new, state->box[d]);
4599 for (k = k0; (k < k1); k++)
4601 rvec_inc(state->x[k], state->box[d]);
4607 copy_rvec(cm_new, cg_cm[cg]);
4609 /* Determine where this cg should go */
4612 for (d = 0; d < dd->ndim; d++)
4617 flag |= DD_FLAG_FW(d);
4623 else if (dev[dim] == -1)
4625 flag |= DD_FLAG_BW(d);
4628 if (dd->nc[dim] > 2)
4639 /* Temporarily store the flag in move */
4640 move[cg] = mc + flag;
4644 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4645 gmx_domdec_t *dd, ivec tric_dir,
4646 t_state *state, rvec **f,
4647 t_forcerec *fr, t_mdatoms *md,
4655 int ncg[DIM*2], nat[DIM*2];
4656 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4657 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4658 int sbuf[2], rbuf[2];
4659 int home_pos_cg, home_pos_at, buf_pos;
4661 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4664 real inv_ncg, pos_d;
4666 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4668 cginfo_mb_t *cginfo_mb;
4669 gmx_domdec_comm_t *comm;
4671 int nthread, thread;
4675 check_screw_box(state->box);
4679 if (fr->cutoff_scheme == ecutsGROUP)
4684 for (i = 0; i < estNR; i++)
4690 case estX: /* Always present */ break;
4691 case estV: bV = (state->flags & (1<<i)); break;
4692 case estSDX: bSDX = (state->flags & (1<<i)); break;
4693 case estCGP: bCGP = (state->flags & (1<<i)); break;
4696 case estDISRE_INITF:
4697 case estDISRE_RM3TAV:
4698 case estORIRE_INITF:
4700 /* No processing required */
4703 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4708 if (dd->ncg_tot > comm->nalloc_int)
4710 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4711 srenew(comm->buf_int, comm->nalloc_int);
4713 move = comm->buf_int;
4715 /* Clear the count */
4716 for (c = 0; c < dd->ndim*2; c++)
4722 npbcdim = dd->npbcdim;
4724 for (d = 0; (d < DIM); d++)
4726 limitd[d] = dd->comm->cellsize_min[d];
4727 if (d >= npbcdim && dd->ci[d] == 0)
4729 cell_x0[d] = -GMX_FLOAT_MAX;
4733 cell_x0[d] = comm->cell_x0[d];
4735 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4737 cell_x1[d] = GMX_FLOAT_MAX;
4741 cell_x1[d] = comm->cell_x1[d];
4745 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4746 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4750 /* We check after communication if a charge group moved
4751 * more than one cell. Set the pre-comm check limit to float_max.
4753 limit0[d] = -GMX_FLOAT_MAX;
4754 limit1[d] = GMX_FLOAT_MAX;
4758 make_tric_corr_matrix(npbcdim, state->box, tcm);
4760 cgindex = dd->cgindex;
4762 nthread = gmx_omp_nthreads_get(emntDomdec);
4764 /* Compute the center of geometry for all home charge groups
4765 * and put them in the box and determine where they should go.
4767 #pragma omp parallel for num_threads(nthread) schedule(static)
4768 for (thread = 0; thread < nthread; thread++)
4770 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4771 cell_x0, cell_x1, limitd, limit0, limit1,
4773 ( thread *dd->ncg_home)/nthread,
4774 ((thread+1)*dd->ncg_home)/nthread,
4775 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4779 for (cg = 0; cg < dd->ncg_home; cg++)
4784 flag = mc & ~DD_FLAG_NRCG;
4785 mc = mc & DD_FLAG_NRCG;
4788 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4790 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4791 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4793 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4794 /* We store the cg size in the lower 16 bits
4795 * and the place where the charge group should go
4796 * in the next 6 bits. This saves some communication volume.
4798 nrcg = cgindex[cg+1] - cgindex[cg];
4799 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4805 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4806 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4809 for (i = 0; i < dd->ndim*2; i++)
4811 *ncg_moved += ncg[i];
4828 /* Make sure the communication buffers are large enough */
4829 for (mc = 0; mc < dd->ndim*2; mc++)
4831 nvr = ncg[mc] + nat[mc]*nvec;
4832 if (nvr > comm->cgcm_state_nalloc[mc])
4834 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4835 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4839 switch (fr->cutoff_scheme)
4842 /* Recalculating cg_cm might be cheaper than communicating,
4843 * but that could give rise to rounding issues.
4846 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4847 nvec, cg_cm, comm, bCompact);
4850 /* Without charge groups we send the moved atom coordinates
4851 * over twice. This is so the code below can be used without
4852 * many conditionals for both for with and without charge groups.
4855 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4856 nvec, state->x, comm, FALSE);
4859 home_pos_cg -= *ncg_moved;
4863 gmx_incons("unimplemented");
4869 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4870 nvec, vec++, state->x, comm, bCompact);
4873 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4874 nvec, vec++, state->v, comm, bCompact);
4878 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4879 nvec, vec++, state->sd_X, comm, bCompact);
4883 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4884 nvec, vec++, state->cg_p, comm, bCompact);
4889 compact_ind(dd->ncg_home, move,
4890 dd->index_gl, dd->cgindex, dd->gatindex,
4891 dd->ga2la, comm->bLocalCG,
4896 if (fr->cutoff_scheme == ecutsVERLET)
4898 moved = get_moved(comm, dd->ncg_home);
4900 for (k = 0; k < dd->ncg_home; k++)
4907 moved = fr->ns.grid->cell_index;
4910 clear_and_mark_ind(dd->ncg_home, move,
4911 dd->index_gl, dd->cgindex, dd->gatindex,
4912 dd->ga2la, comm->bLocalCG,
4916 cginfo_mb = fr->cginfo_mb;
4918 *ncg_stay_home = home_pos_cg;
4919 for (d = 0; d < dd->ndim; d++)
4925 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4928 /* Communicate the cg and atom counts */
4933 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4934 d, dir, sbuf[0], sbuf[1]);
4936 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4938 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4940 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4941 srenew(comm->buf_int, comm->nalloc_int);
4944 /* Communicate the charge group indices, sizes and flags */
4945 dd_sendrecv_int(dd, d, dir,
4946 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4947 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4949 nvs = ncg[cdd] + nat[cdd]*nvec;
4950 i = rbuf[0] + rbuf[1] *nvec;
4951 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4953 /* Communicate cgcm and state */
4954 dd_sendrecv_rvec(dd, d, dir,
4955 comm->cgcm_state[cdd], nvs,
4956 comm->vbuf.v+nvr, i);
4957 ncg_recv += rbuf[0];
4958 nat_recv += rbuf[1];
4962 /* Process the received charge groups */
4964 for (cg = 0; cg < ncg_recv; cg++)
4966 flag = comm->buf_int[cg*DD_CGIBS+1];
4968 if (dim >= npbcdim && dd->nc[dim] > 2)
4970 /* No pbc in this dim and more than one domain boundary.
4971 * We do a separate check if a charge group didn't move too far.
4973 if (((flag & DD_FLAG_FW(d)) &&
4974 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4975 ((flag & DD_FLAG_BW(d)) &&
4976 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4978 cg_move_error(fplog, dd, step, cg, dim,
4979 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4981 comm->vbuf.v[buf_pos],
4982 comm->vbuf.v[buf_pos],
4983 comm->vbuf.v[buf_pos][dim]);
4990 /* Check which direction this cg should go */
4991 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4995 /* The cell boundaries for dimension d2 are not equal
4996 * for each cell row of the lower dimension(s),
4997 * therefore we might need to redetermine where
4998 * this cg should go.
5001 /* If this cg crosses the box boundary in dimension d2
5002 * we can use the communicated flag, so we do not
5003 * have to worry about pbc.
5005 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5006 (flag & DD_FLAG_FW(d2))) ||
5007 (dd->ci[dim2] == 0 &&
5008 (flag & DD_FLAG_BW(d2)))))
5010 /* Clear the two flags for this dimension */
5011 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5012 /* Determine the location of this cg
5013 * in lattice coordinates
5015 pos_d = comm->vbuf.v[buf_pos][dim2];
5018 for (d3 = dim2+1; d3 < DIM; d3++)
5021 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5024 /* Check of we are not at the box edge.
5025 * pbc is only handled in the first step above,
5026 * but this check could move over pbc while
5027 * the first step did not due to different rounding.
5029 if (pos_d >= cell_x1[dim2] &&
5030 dd->ci[dim2] != dd->nc[dim2]-1)
5032 flag |= DD_FLAG_FW(d2);
5034 else if (pos_d < cell_x0[dim2] &&
5037 flag |= DD_FLAG_BW(d2);
5039 comm->buf_int[cg*DD_CGIBS+1] = flag;
5042 /* Set to which neighboring cell this cg should go */
5043 if (flag & DD_FLAG_FW(d2))
5047 else if (flag & DD_FLAG_BW(d2))
5049 if (dd->nc[dd->dim[d2]] > 2)
5061 nrcg = flag & DD_FLAG_NRCG;
5064 if (home_pos_cg+1 > dd->cg_nalloc)
5066 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5067 srenew(dd->index_gl, dd->cg_nalloc);
5068 srenew(dd->cgindex, dd->cg_nalloc+1);
5070 /* Set the global charge group index and size */
5071 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5072 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5073 /* Copy the state from the buffer */
5074 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5075 if (fr->cutoff_scheme == ecutsGROUP)
5078 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5082 /* Set the cginfo */
5083 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5084 dd->index_gl[home_pos_cg]);
5087 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5090 if (home_pos_at+nrcg > state->nalloc)
5092 dd_realloc_state(state, f, home_pos_at+nrcg);
5094 for (i = 0; i < nrcg; i++)
5096 copy_rvec(comm->vbuf.v[buf_pos++],
5097 state->x[home_pos_at+i]);
5101 for (i = 0; i < nrcg; i++)
5103 copy_rvec(comm->vbuf.v[buf_pos++],
5104 state->v[home_pos_at+i]);
5109 for (i = 0; i < nrcg; i++)
5111 copy_rvec(comm->vbuf.v[buf_pos++],
5112 state->sd_X[home_pos_at+i]);
5117 for (i = 0; i < nrcg; i++)
5119 copy_rvec(comm->vbuf.v[buf_pos++],
5120 state->cg_p[home_pos_at+i]);
5124 home_pos_at += nrcg;
5128 /* Reallocate the buffers if necessary */
5129 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5131 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5132 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5134 nvr = ncg[mc] + nat[mc]*nvec;
5135 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5137 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5138 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5140 /* Copy from the receive to the send buffers */
5141 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5142 comm->buf_int + cg*DD_CGIBS,
5143 DD_CGIBS*sizeof(int));
5144 memcpy(comm->cgcm_state[mc][nvr],
5145 comm->vbuf.v[buf_pos],
5146 (1+nrcg*nvec)*sizeof(rvec));
5147 buf_pos += 1 + nrcg*nvec;
5154 /* With sorting (!bCompact) the indices are now only partially up to date
5155 * and ncg_home and nat_home are not the real count, since there are
5156 * "holes" in the arrays for the charge groups that moved to neighbors.
5158 if (fr->cutoff_scheme == ecutsVERLET)
5160 moved = get_moved(comm, home_pos_cg);
5162 for (i = dd->ncg_home; i < home_pos_cg; i++)
5167 dd->ncg_home = home_pos_cg;
5168 dd->nat_home = home_pos_at;
5173 "Finished repartitioning: cgs moved out %d, new home %d\n",
5174 *ncg_moved, dd->ncg_home-*ncg_moved);
5179 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5181 dd->comm->cycl[ddCycl] += cycles;
5182 dd->comm->cycl_n[ddCycl]++;
5183 if (cycles > dd->comm->cycl_max[ddCycl])
5185 dd->comm->cycl_max[ddCycl] = cycles;
5189 static double force_flop_count(t_nrnb *nrnb)
5196 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5198 /* To get closer to the real timings, we half the count
5199 * for the normal loops and again half it for water loops.
5202 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5204 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5208 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5211 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5214 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5216 sum += nrnb->n[i]*cost_nrnb(i);
5219 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5221 sum += nrnb->n[i]*cost_nrnb(i);
5227 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5229 if (dd->comm->eFlop)
5231 dd->comm->flop -= force_flop_count(nrnb);
5234 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5236 if (dd->comm->eFlop)
5238 dd->comm->flop += force_flop_count(nrnb);
5243 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5247 for (i = 0; i < ddCyclNr; i++)
5249 dd->comm->cycl[i] = 0;
5250 dd->comm->cycl_n[i] = 0;
5251 dd->comm->cycl_max[i] = 0;
5254 dd->comm->flop_n = 0;
5257 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5259 gmx_domdec_comm_t *comm;
5260 gmx_domdec_load_t *load;
5261 gmx_domdec_root_t *root = NULL;
5262 int d, dim, cid, i, pos;
5263 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5268 fprintf(debug, "get_load_distribution start\n");
5271 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5275 bSepPME = (dd->pme_nodeid >= 0);
5277 for (d = dd->ndim-1; d >= 0; d--)
5280 /* Check if we participate in the communication in this dimension */
5281 if (d == dd->ndim-1 ||
5282 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5284 load = &comm->load[d];
5287 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5290 if (d == dd->ndim-1)
5292 sbuf[pos++] = dd_force_load(comm);
5293 sbuf[pos++] = sbuf[0];
5296 sbuf[pos++] = sbuf[0];
5297 sbuf[pos++] = cell_frac;
5300 sbuf[pos++] = comm->cell_f_max0[d];
5301 sbuf[pos++] = comm->cell_f_min1[d];
5306 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5307 sbuf[pos++] = comm->cycl[ddCyclPME];
5312 sbuf[pos++] = comm->load[d+1].sum;
5313 sbuf[pos++] = comm->load[d+1].max;
5316 sbuf[pos++] = comm->load[d+1].sum_m;
5317 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5318 sbuf[pos++] = comm->load[d+1].flags;
5321 sbuf[pos++] = comm->cell_f_max0[d];
5322 sbuf[pos++] = comm->cell_f_min1[d];
5327 sbuf[pos++] = comm->load[d+1].mdf;
5328 sbuf[pos++] = comm->load[d+1].pme;
5332 /* Communicate a row in DD direction d.
5333 * The communicators are setup such that the root always has rank 0.
5336 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5337 load->load, load->nload*sizeof(float), MPI_BYTE,
5338 0, comm->mpi_comm_load[d]);
5340 if (dd->ci[dim] == dd->master_ci[dim])
5342 /* We are the root, process this row */
5343 if (comm->bDynLoadBal)
5345 root = comm->root[d];
5355 for (i = 0; i < dd->nc[dim]; i++)
5357 load->sum += load->load[pos++];
5358 load->max = max(load->max, load->load[pos]);
5364 /* This direction could not be load balanced properly,
5365 * therefore we need to use the maximum iso the average load.
5367 load->sum_m = max(load->sum_m, load->load[pos]);
5371 load->sum_m += load->load[pos];
5374 load->cvol_min = min(load->cvol_min, load->load[pos]);
5378 load->flags = (int)(load->load[pos++] + 0.5);
5382 root->cell_f_max0[i] = load->load[pos++];
5383 root->cell_f_min1[i] = load->load[pos++];
5388 load->mdf = max(load->mdf, load->load[pos]);
5390 load->pme = max(load->pme, load->load[pos]);
5394 if (comm->bDynLoadBal && root->bLimited)
5396 load->sum_m *= dd->nc[dim];
5397 load->flags |= (1<<d);
5405 comm->nload += dd_load_count(comm);
5406 comm->load_step += comm->cycl[ddCyclStep];
5407 comm->load_sum += comm->load[0].sum;
5408 comm->load_max += comm->load[0].max;
5409 if (comm->bDynLoadBal)
5411 for (d = 0; d < dd->ndim; d++)
5413 if (comm->load[0].flags & (1<<d))
5415 comm->load_lim[d]++;
5421 comm->load_mdf += comm->load[0].mdf;
5422 comm->load_pme += comm->load[0].pme;
5426 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5430 fprintf(debug, "get_load_distribution finished\n");
5434 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5436 /* Return the relative performance loss on the total run time
5437 * due to the force calculation load imbalance.
5439 if (dd->comm->nload > 0)
5442 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5443 (dd->comm->load_step*dd->nnodes);
5451 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5454 int npp, npme, nnodes, d, limp;
5455 float imbal, pme_f_ratio, lossf, lossp = 0;
5457 gmx_domdec_comm_t *comm;
5460 if (DDMASTER(dd) && comm->nload > 0)
5463 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5464 nnodes = npp + npme;
5465 imbal = comm->load_max*npp/comm->load_sum - 1;
5466 lossf = dd_force_imb_perf_loss(dd);
5467 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5468 fprintf(fplog, "%s", buf);
5469 fprintf(stderr, "\n");
5470 fprintf(stderr, "%s", buf);
5471 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5472 fprintf(fplog, "%s", buf);
5473 fprintf(stderr, "%s", buf);
5475 if (comm->bDynLoadBal)
5477 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5478 for (d = 0; d < dd->ndim; d++)
5480 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5481 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5487 sprintf(buf+strlen(buf), "\n");
5488 fprintf(fplog, "%s", buf);
5489 fprintf(stderr, "%s", buf);
5493 pme_f_ratio = comm->load_pme/comm->load_mdf;
5494 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5497 lossp *= (float)npme/(float)nnodes;
5501 lossp *= (float)npp/(float)nnodes;
5503 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5504 fprintf(fplog, "%s", buf);
5505 fprintf(stderr, "%s", buf);
5506 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5507 fprintf(fplog, "%s", buf);
5508 fprintf(stderr, "%s", buf);
5510 fprintf(fplog, "\n");
5511 fprintf(stderr, "\n");
5513 if (lossf >= DD_PERF_LOSS)
5516 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5517 " in the domain decomposition.\n", lossf*100);
5518 if (!comm->bDynLoadBal)
5520 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5524 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5526 fprintf(fplog, "%s\n", buf);
5527 fprintf(stderr, "%s\n", buf);
5529 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5532 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5533 " had %s work to do than the PP nodes.\n"
5534 " You might want to %s the number of PME nodes\n"
5535 " or %s the cut-off and the grid spacing.\n",
5537 (lossp < 0) ? "less" : "more",
5538 (lossp < 0) ? "decrease" : "increase",
5539 (lossp < 0) ? "decrease" : "increase");
5540 fprintf(fplog, "%s\n", buf);
5541 fprintf(stderr, "%s\n", buf);
5546 static float dd_vol_min(gmx_domdec_t *dd)
5548 return dd->comm->load[0].cvol_min*dd->nnodes;
5551 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5553 return dd->comm->load[0].flags;
5556 static float dd_f_imbal(gmx_domdec_t *dd)
5558 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5561 float dd_pme_f_ratio(gmx_domdec_t *dd)
5563 if (dd->comm->cycl_n[ddCyclPME] > 0)
5565 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5573 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5578 flags = dd_load_flags(dd);
5582 "DD load balancing is limited by minimum cell size in dimension");
5583 for (d = 0; d < dd->ndim; d++)
5587 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5590 fprintf(fplog, "\n");
5592 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5593 if (dd->comm->bDynLoadBal)
5595 fprintf(fplog, " vol min/aver %5.3f%c",
5596 dd_vol_min(dd), flags ? '!' : ' ');
5598 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5599 if (dd->comm->cycl_n[ddCyclPME])
5601 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5603 fprintf(fplog, "\n\n");
5606 static void dd_print_load_verbose(gmx_domdec_t *dd)
5608 if (dd->comm->bDynLoadBal)
5610 fprintf(stderr, "vol %4.2f%c ",
5611 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5613 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5614 if (dd->comm->cycl_n[ddCyclPME])
5616 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5621 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5626 gmx_domdec_root_t *root;
5627 gmx_bool bPartOfGroup = FALSE;
5629 dim = dd->dim[dim_ind];
5630 copy_ivec(loc, loc_c);
5631 for (i = 0; i < dd->nc[dim]; i++)
5634 rank = dd_index(dd->nc, loc_c);
5635 if (rank == dd->rank)
5637 /* This process is part of the group */
5638 bPartOfGroup = TRUE;
5641 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5645 dd->comm->mpi_comm_load[dim_ind] = c_row;
5646 if (dd->comm->eDLB != edlbNO)
5648 if (dd->ci[dim] == dd->master_ci[dim])
5650 /* This is the root process of this row */
5651 snew(dd->comm->root[dim_ind], 1);
5652 root = dd->comm->root[dim_ind];
5653 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5654 snew(root->old_cell_f, dd->nc[dim]+1);
5655 snew(root->bCellMin, dd->nc[dim]);
5658 snew(root->cell_f_max0, dd->nc[dim]);
5659 snew(root->cell_f_min1, dd->nc[dim]);
5660 snew(root->bound_min, dd->nc[dim]);
5661 snew(root->bound_max, dd->nc[dim]);
5663 snew(root->buf_ncd, dd->nc[dim]);
5667 /* This is not a root process, we only need to receive cell_f */
5668 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5671 if (dd->ci[dim] == dd->master_ci[dim])
5673 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5679 void dd_setup_dlb_resource_sharing(t_commrec *cr,
5680 const gmx_hw_info_t *hwinfo,
5681 const gmx_hw_opt_t *hw_opt)
5684 int physicalnode_id_hash;
5687 MPI_Comm mpi_comm_pp_physicalnode;
5689 if (!(cr->duty & DUTY_PP) ||
5690 hw_opt->gpu_opt.ncuda_dev_use == 0)
5692 /* Only PP nodes (currently) use GPUs.
5693 * If we don't have GPUs, there are no resources to share.
5698 physicalnode_id_hash = gmx_physicalnode_id_hash();
5700 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5706 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5707 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5708 dd->rank, physicalnode_id_hash, gpu_id);
5710 /* Split the PP communicator over the physical nodes */
5711 /* TODO: See if we should store this (before), as it's also used for
5712 * for the nodecomm summution.
5714 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5715 &mpi_comm_pp_physicalnode);
5716 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5717 &dd->comm->mpi_comm_gpu_shared);
5718 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5719 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5723 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5726 /* Note that some ranks could share a GPU, while others don't */
5728 if (dd->comm->nrank_gpu_shared == 1)
5730 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5735 static void make_load_communicators(gmx_domdec_t *dd)
5738 int dim0, dim1, i, j;
5743 fprintf(debug, "Making load communicators\n");
5746 snew(dd->comm->load, dd->ndim);
5747 snew(dd->comm->mpi_comm_load, dd->ndim);
5750 make_load_communicator(dd, 0, loc);
5754 for (i = 0; i < dd->nc[dim0]; i++)
5757 make_load_communicator(dd, 1, loc);
5763 for (i = 0; i < dd->nc[dim0]; i++)
5767 for (j = 0; j < dd->nc[dim1]; j++)
5770 make_load_communicator(dd, 2, loc);
5777 fprintf(debug, "Finished making load communicators\n");
5782 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5785 int d, dim, i, j, m;
5788 ivec dd_zp[DD_MAXIZONE];
5789 gmx_domdec_zones_t *zones;
5790 gmx_domdec_ns_ranges_t *izone;
5792 for (d = 0; d < dd->ndim; d++)
5795 copy_ivec(dd->ci, tmp);
5796 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5797 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5798 copy_ivec(dd->ci, tmp);
5799 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5800 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5803 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5806 dd->neighbor[d][1]);
5812 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5814 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5815 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5822 for (i = 0; i < nzonep; i++)
5824 copy_ivec(dd_zp3[i], dd_zp[i]);
5830 for (i = 0; i < nzonep; i++)
5832 copy_ivec(dd_zp2[i], dd_zp[i]);
5838 for (i = 0; i < nzonep; i++)
5840 copy_ivec(dd_zp1[i], dd_zp[i]);
5844 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5849 zones = &dd->comm->zones;
5851 for (i = 0; i < nzone; i++)
5854 clear_ivec(zones->shift[i]);
5855 for (d = 0; d < dd->ndim; d++)
5857 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5862 for (i = 0; i < nzone; i++)
5864 for (d = 0; d < DIM; d++)
5866 s[d] = dd->ci[d] - zones->shift[i][d];
5871 else if (s[d] >= dd->nc[d])
5877 zones->nizone = nzonep;
5878 for (i = 0; i < zones->nizone; i++)
5880 if (dd_zp[i][0] != i)
5882 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5884 izone = &zones->izone[i];
5885 izone->j0 = dd_zp[i][1];
5886 izone->j1 = dd_zp[i][2];
5887 for (dim = 0; dim < DIM; dim++)
5889 if (dd->nc[dim] == 1)
5891 /* All shifts should be allowed */
5892 izone->shift0[dim] = -1;
5893 izone->shift1[dim] = 1;
5898 izone->shift0[d] = 0;
5899 izone->shift1[d] = 0;
5900 for(j=izone->j0; j<izone->j1; j++) {
5901 if (dd->shift[j][d] > dd->shift[i][d])
5902 izone->shift0[d] = -1;
5903 if (dd->shift[j][d] < dd->shift[i][d])
5904 izone->shift1[d] = 1;
5910 /* Assume the shift are not more than 1 cell */
5911 izone->shift0[dim] = 1;
5912 izone->shift1[dim] = -1;
5913 for (j = izone->j0; j < izone->j1; j++)
5915 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5916 if (shift_diff < izone->shift0[dim])
5918 izone->shift0[dim] = shift_diff;
5920 if (shift_diff > izone->shift1[dim])
5922 izone->shift1[dim] = shift_diff;
5929 if (dd->comm->eDLB != edlbNO)
5931 snew(dd->comm->root, dd->ndim);
5934 if (dd->comm->bRecordLoad)
5936 make_load_communicators(dd);
5940 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5943 gmx_domdec_comm_t *comm;
5954 if (comm->bCartesianPP)
5956 /* Set up cartesian communication for the particle-particle part */
5959 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5960 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5963 for (i = 0; i < DIM; i++)
5967 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5969 /* We overwrite the old communicator with the new cartesian one */
5970 cr->mpi_comm_mygroup = comm_cart;
5973 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5974 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5976 if (comm->bCartesianPP_PME)
5978 /* Since we want to use the original cartesian setup for sim,
5979 * and not the one after split, we need to make an index.
5981 snew(comm->ddindex2ddnodeid, dd->nnodes);
5982 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5983 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5984 /* Get the rank of the DD master,
5985 * above we made sure that the master node is a PP node.
5995 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5997 else if (comm->bCartesianPP)
5999 if (cr->npmenodes == 0)
6001 /* The PP communicator is also
6002 * the communicator for this simulation
6004 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6006 cr->nodeid = dd->rank;
6008 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6010 /* We need to make an index to go from the coordinates
6011 * to the nodeid of this simulation.
6013 snew(comm->ddindex2simnodeid, dd->nnodes);
6014 snew(buf, dd->nnodes);
6015 if (cr->duty & DUTY_PP)
6017 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6019 /* Communicate the ddindex to simulation nodeid index */
6020 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6021 cr->mpi_comm_mysim);
6024 /* Determine the master coordinates and rank.
6025 * The DD master should be the same node as the master of this sim.
6027 for (i = 0; i < dd->nnodes; i++)
6029 if (comm->ddindex2simnodeid[i] == 0)
6031 ddindex2xyz(dd->nc, i, dd->master_ci);
6032 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6037 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6042 /* No Cartesian communicators */
6043 /* We use the rank in dd->comm->all as DD index */
6044 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6045 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6047 clear_ivec(dd->master_ci);
6054 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
6055 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6060 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
6061 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6065 static void receive_ddindex2simnodeid(t_commrec *cr)
6069 gmx_domdec_comm_t *comm;
6076 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6078 snew(comm->ddindex2simnodeid, dd->nnodes);
6079 snew(buf, dd->nnodes);
6080 if (cr->duty & DUTY_PP)
6082 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6085 /* Communicate the ddindex to simulation nodeid index */
6086 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6087 cr->mpi_comm_mysim);
6094 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6095 int ncg, int natoms)
6097 gmx_domdec_master_t *ma;
6102 snew(ma->ncg, dd->nnodes);
6103 snew(ma->index, dd->nnodes+1);
6105 snew(ma->nat, dd->nnodes);
6106 snew(ma->ibuf, dd->nnodes*2);
6107 snew(ma->cell_x, DIM);
6108 for (i = 0; i < DIM; i++)
6110 snew(ma->cell_x[i], dd->nc[i]+1);
6113 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6119 snew(ma->vbuf, natoms);
6125 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
6129 gmx_domdec_comm_t *comm;
6140 if (comm->bCartesianPP)
6142 for (i = 1; i < DIM; i++)
6144 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6146 if (bDiv[YY] || bDiv[ZZ])
6148 comm->bCartesianPP_PME = TRUE;
6149 /* If we have 2D PME decomposition, which is always in x+y,
6150 * we stack the PME only nodes in z.
6151 * Otherwise we choose the direction that provides the thinnest slab
6152 * of PME only nodes as this will have the least effect
6153 * on the PP communication.
6154 * But for the PME communication the opposite might be better.
6156 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6158 dd->nc[YY] > dd->nc[ZZ]))
6160 comm->cartpmedim = ZZ;
6164 comm->cartpmedim = YY;
6166 comm->ntot[comm->cartpmedim]
6167 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6171 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]);
6173 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6178 if (comm->bCartesianPP_PME)
6182 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]);
6185 for (i = 0; i < DIM; i++)
6189 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6192 MPI_Comm_rank(comm_cart, &rank);
6193 if (MASTERNODE(cr) && rank != 0)
6195 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6198 /* With this assigment we loose the link to the original communicator
6199 * which will usually be MPI_COMM_WORLD, unless have multisim.
6201 cr->mpi_comm_mysim = comm_cart;
6202 cr->sim_nodeid = rank;
6204 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6208 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6209 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6212 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6216 if (cr->npmenodes == 0 ||
6217 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6219 cr->duty = DUTY_PME;
6222 /* Split the sim communicator into PP and PME only nodes */
6223 MPI_Comm_split(cr->mpi_comm_mysim,
6225 dd_index(comm->ntot, dd->ci),
6226 &cr->mpi_comm_mygroup);
6230 switch (dd_node_order)
6235 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6238 case ddnoINTERLEAVE:
6239 /* Interleave the PP-only and PME-only nodes,
6240 * as on clusters with dual-core machines this will double
6241 * the communication bandwidth of the PME processes
6242 * and thus speed up the PP <-> PME and inter PME communication.
6246 fprintf(fplog, "Interleaving PP and PME nodes\n");
6248 comm->pmenodes = dd_pmenodes(cr);
6253 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6256 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6258 cr->duty = DUTY_PME;
6265 /* Split the sim communicator into PP and PME only nodes */
6266 MPI_Comm_split(cr->mpi_comm_mysim,
6269 &cr->mpi_comm_mygroup);
6270 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6276 fprintf(fplog, "This is a %s only node\n\n",
6277 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6281 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6284 gmx_domdec_comm_t *comm;
6290 copy_ivec(dd->nc, comm->ntot);
6292 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6293 comm->bCartesianPP_PME = FALSE;
6295 /* Reorder the nodes by default. This might change the MPI ranks.
6296 * Real reordering is only supported on very few architectures,
6297 * Blue Gene is one of them.
6299 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6301 if (cr->npmenodes > 0)
6303 /* Split the communicator into a PP and PME part */
6304 split_communicator(fplog, cr, dd_node_order, CartReorder);
6305 if (comm->bCartesianPP_PME)
6307 /* We (possibly) reordered the nodes in split_communicator,
6308 * so it is no longer required in make_pp_communicator.
6310 CartReorder = FALSE;
6315 /* All nodes do PP and PME */
6317 /* We do not require separate communicators */
6318 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6322 if (cr->duty & DUTY_PP)
6324 /* Copy or make a new PP communicator */
6325 make_pp_communicator(fplog, cr, CartReorder);
6329 receive_ddindex2simnodeid(cr);
6332 if (!(cr->duty & DUTY_PME))
6334 /* Set up the commnuication to our PME node */
6335 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6336 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6339 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6340 dd->pme_nodeid, dd->pme_receive_vir_ener);
6345 dd->pme_nodeid = -1;
6350 dd->ma = init_gmx_domdec_master_t(dd,
6352 comm->cgs_gl.index[comm->cgs_gl.nr]);
6356 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6358 real *slb_frac, tot;
6363 if (nc > 1 && size_string != NULL)
6367 fprintf(fplog, "Using static load balancing for the %s direction\n",
6372 for (i = 0; i < nc; i++)
6375 sscanf(size_string, "%lf%n", &dbl, &n);
6378 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6387 fprintf(fplog, "Relative cell sizes:");
6389 for (i = 0; i < nc; i++)
6394 fprintf(fplog, " %5.3f", slb_frac[i]);
6399 fprintf(fplog, "\n");
6406 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6409 gmx_mtop_ilistloop_t iloop;
6413 iloop = gmx_mtop_ilistloop_init(mtop);
6414 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6416 for (ftype = 0; ftype < F_NRE; ftype++)
6418 if ((interaction_function[ftype].flags & IF_BOND) &&
6421 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6429 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6435 val = getenv(env_var);
6438 if (sscanf(val, "%d", &nst) <= 0)
6444 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6452 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6456 fprintf(stderr, "\n%s\n", warn_string);
6460 fprintf(fplog, "\n%s\n", warn_string);
6464 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6465 t_inputrec *ir, FILE *fplog)
6467 if (ir->ePBC == epbcSCREW &&
6468 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6470 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6473 if (ir->ns_type == ensSIMPLE)
6475 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6478 if (ir->nstlist == 0)
6480 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6483 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6485 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6489 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6494 r = ddbox->box_size[XX];
6495 for (di = 0; di < dd->ndim; di++)
6498 /* Check using the initial average cell size */
6499 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6505 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6506 const char *dlb_opt, gmx_bool bRecordLoad,
6507 unsigned long Flags, t_inputrec *ir)
6515 case 'a': eDLB = edlbAUTO; break;
6516 case 'n': eDLB = edlbNO; break;
6517 case 'y': eDLB = edlbYES; break;
6518 default: gmx_incons("Unknown dlb_opt");
6521 if (Flags & MD_RERUN)
6526 if (!EI_DYNAMICS(ir->eI))
6528 if (eDLB == edlbYES)
6530 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6531 dd_warning(cr, fplog, buf);
6539 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6544 if (Flags & MD_REPRODUCIBLE)
6551 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6555 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6558 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6566 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6571 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6573 /* Decomposition order z,y,x */
6576 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6578 for (dim = DIM-1; dim >= 0; dim--)
6580 if (dd->nc[dim] > 1)
6582 dd->dim[dd->ndim++] = dim;
6588 /* Decomposition order x,y,z */
6589 for (dim = 0; dim < DIM; dim++)
6591 if (dd->nc[dim] > 1)
6593 dd->dim[dd->ndim++] = dim;
6599 static gmx_domdec_comm_t *init_dd_comm()
6601 gmx_domdec_comm_t *comm;
6605 snew(comm->cggl_flag, DIM*2);
6606 snew(comm->cgcm_state, DIM*2);
6607 for (i = 0; i < DIM*2; i++)
6609 comm->cggl_flag_nalloc[i] = 0;
6610 comm->cgcm_state_nalloc[i] = 0;
6613 comm->nalloc_int = 0;
6614 comm->buf_int = NULL;
6616 vec_rvec_init(&comm->vbuf);
6618 comm->n_load_have = 0;
6619 comm->n_load_collect = 0;
6621 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6623 comm->sum_nat[i] = 0;
6627 comm->load_step = 0;
6630 clear_ivec(comm->load_lim);
6637 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6638 unsigned long Flags,
6640 real comm_distance_min, real rconstr,
6641 const char *dlb_opt, real dlb_scale,
6642 const char *sizex, const char *sizey, const char *sizez,
6643 gmx_mtop_t *mtop, t_inputrec *ir,
6644 matrix box, rvec *x,
6646 int *npme_x, int *npme_y)
6649 gmx_domdec_comm_t *comm;
6652 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6659 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6664 dd->comm = init_dd_comm();
6666 snew(comm->cggl_flag, DIM*2);
6667 snew(comm->cgcm_state, DIM*2);
6669 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6670 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6672 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6673 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6674 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6675 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6676 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6677 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6678 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6679 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6681 dd->pme_recv_f_alloc = 0;
6682 dd->pme_recv_f_buf = NULL;
6684 if (dd->bSendRecv2 && fplog)
6686 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");
6692 fprintf(fplog, "Will load balance based on FLOP count\n");
6694 if (comm->eFlop > 1)
6696 srand(1+cr->nodeid);
6698 comm->bRecordLoad = TRUE;
6702 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6706 /* Initialize to GPU share count to 0, might change later */
6707 comm->nrank_gpu_shared = 0;
6709 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6711 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6714 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6716 dd->bGridJump = comm->bDynLoadBal;
6717 comm->bPMELoadBalDLBLimits = FALSE;
6719 if (comm->nstSortCG)
6723 if (comm->nstSortCG == 1)
6725 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6729 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6733 snew(comm->sort, 1);
6739 fprintf(fplog, "Will not sort the charge groups\n");
6743 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6745 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6746 if (comm->bInterCGBondeds)
6748 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6752 comm->bInterCGMultiBody = FALSE;
6755 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6756 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6758 if (ir->rlistlong == 0)
6760 /* Set the cut-off to some very large value,
6761 * so we don't need if statements everywhere in the code.
6762 * We use sqrt, since the cut-off is squared in some places.
6764 comm->cutoff = GMX_CUTOFF_INF;
6768 comm->cutoff = ir->rlistlong;
6770 comm->cutoff_mbody = 0;
6772 comm->cellsize_limit = 0;
6773 comm->bBondComm = FALSE;
6775 if (comm->bInterCGBondeds)
6777 if (comm_distance_min > 0)
6779 comm->cutoff_mbody = comm_distance_min;
6780 if (Flags & MD_DDBONDCOMM)
6782 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6786 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6788 r_bonded_limit = comm->cutoff_mbody;
6790 else if (ir->bPeriodicMols)
6792 /* Can not easily determine the required cut-off */
6793 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");
6794 comm->cutoff_mbody = comm->cutoff/2;
6795 r_bonded_limit = comm->cutoff_mbody;
6801 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6802 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6804 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6805 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6807 /* We use an initial margin of 10% for the minimum cell size,
6808 * except when we are just below the non-bonded cut-off.
6810 if (Flags & MD_DDBONDCOMM)
6812 if (max(r_2b, r_mb) > comm->cutoff)
6814 r_bonded = max(r_2b, r_mb);
6815 r_bonded_limit = 1.1*r_bonded;
6816 comm->bBondComm = TRUE;
6821 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6823 /* We determine cutoff_mbody later */
6827 /* No special bonded communication,
6828 * simply increase the DD cut-off.
6830 r_bonded_limit = 1.1*max(r_2b, r_mb);
6831 comm->cutoff_mbody = r_bonded_limit;
6832 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6835 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6839 "Minimum cell size due to bonded interactions: %.3f nm\n",
6840 comm->cellsize_limit);
6844 if (dd->bInterCGcons && rconstr <= 0)
6846 /* There is a cell size limit due to the constraints (P-LINCS) */
6847 rconstr = constr_r_max(fplog, mtop, ir);
6851 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6853 if (rconstr > comm->cellsize_limit)
6855 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6859 else if (rconstr > 0 && fplog)
6861 /* Here we do not check for dd->bInterCGcons,
6862 * because one can also set a cell size limit for virtual sites only
6863 * and at this point we don't know yet if there are intercg v-sites.
6866 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6869 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6871 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6875 copy_ivec(nc, dd->nc);
6876 set_dd_dim(fplog, dd);
6877 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6879 if (cr->npmenodes == -1)
6883 acs = average_cellsize_min(dd, ddbox);
6884 if (acs < comm->cellsize_limit)
6888 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6890 gmx_fatal_collective(FARGS, cr, NULL,
6891 "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",
6892 acs, comm->cellsize_limit);
6897 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6899 /* We need to choose the optimal DD grid and possibly PME nodes */
6900 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6901 comm->eDLB != edlbNO, dlb_scale,
6902 comm->cellsize_limit, comm->cutoff,
6903 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6905 if (dd->nc[XX] == 0)
6907 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6908 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6909 !bC ? "-rdd" : "-rcon",
6910 comm->eDLB != edlbNO ? " or -dds" : "",
6911 bC ? " or your LINCS settings" : "");
6913 gmx_fatal_collective(FARGS, cr, NULL,
6914 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6916 "Look in the log file for details on the domain decomposition",
6917 cr->nnodes-cr->npmenodes, limit, buf);
6919 set_dd_dim(fplog, dd);
6925 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6926 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6929 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6930 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6932 gmx_fatal_collective(FARGS, cr, NULL,
6933 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6934 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6936 if (cr->npmenodes > dd->nnodes)
6938 gmx_fatal_collective(FARGS, cr, NULL,
6939 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6941 if (cr->npmenodes > 0)
6943 comm->npmenodes = cr->npmenodes;
6947 comm->npmenodes = dd->nnodes;
6950 if (EEL_PME(ir->coulombtype))
6952 /* The following choices should match those
6953 * in comm_cost_est in domdec_setup.c.
6954 * Note that here the checks have to take into account
6955 * that the decomposition might occur in a different order than xyz
6956 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6957 * in which case they will not match those in comm_cost_est,
6958 * but since that is mainly for testing purposes that's fine.
6960 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6961 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6962 getenv("GMX_PMEONEDD") == NULL)
6964 comm->npmedecompdim = 2;
6965 comm->npmenodes_x = dd->nc[XX];
6966 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6970 /* In case nc is 1 in both x and y we could still choose to
6971 * decompose pme in y instead of x, but we use x for simplicity.
6973 comm->npmedecompdim = 1;
6974 if (dd->dim[0] == YY)
6976 comm->npmenodes_x = 1;
6977 comm->npmenodes_y = comm->npmenodes;
6981 comm->npmenodes_x = comm->npmenodes;
6982 comm->npmenodes_y = 1;
6987 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6988 comm->npmenodes_x, comm->npmenodes_y, 1);
6993 comm->npmedecompdim = 0;
6994 comm->npmenodes_x = 0;
6995 comm->npmenodes_y = 0;
6998 /* Technically we don't need both of these,
6999 * but it simplifies code not having to recalculate it.
7001 *npme_x = comm->npmenodes_x;
7002 *npme_y = comm->npmenodes_y;
7004 snew(comm->slb_frac, DIM);
7005 if (comm->eDLB == edlbNO)
7007 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7008 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7009 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7012 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7014 if (comm->bBondComm || comm->eDLB != edlbNO)
7016 /* Set the bonded communication distance to halfway
7017 * the minimum and the maximum,
7018 * since the extra communication cost is nearly zero.
7020 acs = average_cellsize_min(dd, ddbox);
7021 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7022 if (comm->eDLB != edlbNO)
7024 /* Check if this does not limit the scaling */
7025 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
7027 if (!comm->bBondComm)
7029 /* Without bBondComm do not go beyond the n.b. cut-off */
7030 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
7031 if (comm->cellsize_limit >= comm->cutoff)
7033 /* We don't loose a lot of efficieny
7034 * when increasing it to the n.b. cut-off.
7035 * It can even be slightly faster, because we need
7036 * less checks for the communication setup.
7038 comm->cutoff_mbody = comm->cutoff;
7041 /* Check if we did not end up below our original limit */
7042 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
7044 if (comm->cutoff_mbody > comm->cellsize_limit)
7046 comm->cellsize_limit = comm->cutoff_mbody;
7049 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7054 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7055 "cellsize limit %f\n",
7056 comm->bBondComm, comm->cellsize_limit);
7061 check_dd_restrictions(cr, dd, ir, fplog);
7064 comm->partition_step = INT_MIN;
7067 clear_dd_cycle_counts(dd);
7072 static void set_dlb_limits(gmx_domdec_t *dd)
7077 for (d = 0; d < dd->ndim; d++)
7079 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7080 dd->comm->cellsize_min[dd->dim[d]] =
7081 dd->comm->cellsize_min_dlb[dd->dim[d]];
7086 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
7089 gmx_domdec_comm_t *comm;
7099 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);
7102 cellsize_min = comm->cellsize_min[dd->dim[0]];
7103 for (d = 1; d < dd->ndim; d++)
7105 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7108 if (cellsize_min < comm->cellsize_limit*1.05)
7110 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");
7112 /* Change DLB from "auto" to "no". */
7113 comm->eDLB = edlbNO;
7118 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7119 comm->bDynLoadBal = TRUE;
7120 dd->bGridJump = TRUE;
7124 /* We can set the required cell size info here,
7125 * so we do not need to communicate this.
7126 * The grid is completely uniform.
7128 for (d = 0; d < dd->ndim; d++)
7132 comm->load[d].sum_m = comm->load[d].sum;
7134 nc = dd->nc[dd->dim[d]];
7135 for (i = 0; i < nc; i++)
7137 comm->root[d]->cell_f[i] = i/(real)nc;
7140 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7141 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7144 comm->root[d]->cell_f[nc] = 1.0;
7149 static char *init_bLocalCG(gmx_mtop_t *mtop)
7154 ncg = ncg_mtop(mtop);
7155 snew(bLocalCG, ncg);
7156 for (cg = 0; cg < ncg; cg++)
7158 bLocalCG[cg] = FALSE;
7164 void dd_init_bondeds(FILE *fplog,
7165 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7166 gmx_vsite_t *vsite, gmx_constr_t constr,
7167 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7169 gmx_domdec_comm_t *comm;
7173 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7177 if (comm->bBondComm)
7179 /* Communicate atoms beyond the cut-off for bonded interactions */
7182 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7184 comm->bLocalCG = init_bLocalCG(mtop);
7188 /* Only communicate atoms based on cut-off */
7189 comm->cglink = NULL;
7190 comm->bLocalCG = NULL;
7194 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7196 gmx_bool bDynLoadBal, real dlb_scale,
7199 gmx_domdec_comm_t *comm;
7214 fprintf(fplog, "The maximum number of communication pulses is:");
7215 for (d = 0; d < dd->ndim; d++)
7217 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7219 fprintf(fplog, "\n");
7220 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7221 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7222 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7223 for (d = 0; d < DIM; d++)
7227 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7234 comm->cellsize_min_dlb[d]/
7235 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7237 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7240 fprintf(fplog, "\n");
7244 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7245 fprintf(fplog, "The initial number of communication pulses is:");
7246 for (d = 0; d < dd->ndim; d++)
7248 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7250 fprintf(fplog, "\n");
7251 fprintf(fplog, "The initial domain decomposition cell size is:");
7252 for (d = 0; d < DIM; d++)
7256 fprintf(fplog, " %c %.2f nm",
7257 dim2char(d), dd->comm->cellsize_min[d]);
7260 fprintf(fplog, "\n\n");
7263 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7265 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7266 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7267 "non-bonded interactions", "", comm->cutoff);
7271 limit = dd->comm->cellsize_limit;
7275 if (dynamic_dd_box(ddbox, ir))
7277 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7279 limit = dd->comm->cellsize_min[XX];
7280 for (d = 1; d < DIM; d++)
7282 limit = min(limit, dd->comm->cellsize_min[d]);
7286 if (comm->bInterCGBondeds)
7288 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7289 "two-body bonded interactions", "(-rdd)",
7290 max(comm->cutoff, comm->cutoff_mbody));
7291 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7292 "multi-body bonded interactions", "(-rdd)",
7293 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7297 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7298 "virtual site constructions", "(-rcon)", limit);
7300 if (dd->constraint_comm)
7302 sprintf(buf, "atoms separated by up to %d constraints",
7304 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7305 buf, "(-rcon)", limit);
7307 fprintf(fplog, "\n");
7313 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7315 const t_inputrec *ir,
7316 const gmx_ddbox_t *ddbox)
7318 gmx_domdec_comm_t *comm;
7319 int d, dim, npulse, npulse_d_max, npulse_d;
7324 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7326 /* Determine the maximum number of comm. pulses in one dimension */
7328 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7330 /* Determine the maximum required number of grid pulses */
7331 if (comm->cellsize_limit >= comm->cutoff)
7333 /* Only a single pulse is required */
7336 else if (!bNoCutOff && comm->cellsize_limit > 0)
7338 /* We round down slightly here to avoid overhead due to the latency
7339 * of extra communication calls when the cut-off
7340 * would be only slightly longer than the cell size.
7341 * Later cellsize_limit is redetermined,
7342 * so we can not miss interactions due to this rounding.
7344 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7348 /* There is no cell size limit */
7349 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7352 if (!bNoCutOff && npulse > 1)
7354 /* See if we can do with less pulses, based on dlb_scale */
7356 for (d = 0; d < dd->ndim; d++)
7359 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7360 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7361 npulse_d_max = max(npulse_d_max, npulse_d);
7363 npulse = min(npulse, npulse_d_max);
7366 /* This env var can override npulse */
7367 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7374 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7375 for (d = 0; d < dd->ndim; d++)
7377 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7378 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7379 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7380 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7381 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7383 comm->bVacDLBNoLimit = FALSE;
7387 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7388 if (!comm->bVacDLBNoLimit)
7390 comm->cellsize_limit = max(comm->cellsize_limit,
7391 comm->cutoff/comm->maxpulse);
7393 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7394 /* Set the minimum cell size for each DD dimension */
7395 for (d = 0; d < dd->ndim; d++)
7397 if (comm->bVacDLBNoLimit ||
7398 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7400 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7404 comm->cellsize_min_dlb[dd->dim[d]] =
7405 comm->cutoff/comm->cd[d].np_dlb;
7408 if (comm->cutoff_mbody <= 0)
7410 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7412 if (comm->bDynLoadBal)
7418 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7420 /* If each molecule is a single charge group
7421 * or we use domain decomposition for each periodic dimension,
7422 * we do not need to take pbc into account for the bonded interactions.
7424 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7427 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7430 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7431 t_inputrec *ir, t_forcerec *fr,
7434 gmx_domdec_comm_t *comm;
7440 /* Initialize the thread data.
7441 * This can not be done in init_domain_decomposition,
7442 * as the numbers of threads is determined later.
7444 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7447 snew(comm->dth, comm->nth);
7450 if (EEL_PME(ir->coulombtype))
7452 init_ddpme(dd, &comm->ddpme[0], 0);
7453 if (comm->npmedecompdim >= 2)
7455 init_ddpme(dd, &comm->ddpme[1], 1);
7460 comm->npmenodes = 0;
7461 if (dd->pme_nodeid >= 0)
7463 gmx_fatal_collective(FARGS, NULL, dd,
7464 "Can not have separate PME nodes without PME electrostatics");
7470 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7472 if (comm->eDLB != edlbNO)
7474 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7477 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7478 if (comm->eDLB == edlbAUTO)
7482 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7484 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7487 if (ir->ePBC == epbcNONE)
7489 vol_frac = 1 - 1/(double)dd->nnodes;
7494 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7498 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7500 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7502 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7505 static gmx_bool test_dd_cutoff(t_commrec *cr,
7506 t_state *state, t_inputrec *ir,
7517 set_ddbox(dd, FALSE, cr, ir, state->box,
7518 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7522 for (d = 0; d < dd->ndim; d++)
7526 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7527 if (dynamic_dd_box(&ddbox, ir))
7529 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7532 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7534 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7535 dd->comm->cd[d].np_dlb > 0)
7537 if (np > dd->comm->cd[d].np_dlb)
7542 /* If a current local cell size is smaller than the requested
7543 * cut-off, we could still fix it, but this gets very complicated.
7544 * Without fixing here, we might actually need more checks.
7546 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7553 if (dd->comm->eDLB != edlbNO)
7555 /* If DLB is not active yet, we don't need to check the grid jumps.
7556 * Actually we shouldn't, because then the grid jump data is not set.
7558 if (dd->comm->bDynLoadBal &&
7559 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7564 gmx_sumi(1, &LocallyLimited, cr);
7566 if (LocallyLimited > 0)
7575 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7578 gmx_bool bCutoffAllowed;
7580 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7584 cr->dd->comm->cutoff = cutoff_req;
7587 return bCutoffAllowed;
7590 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7592 gmx_domdec_comm_t *comm;
7594 comm = cr->dd->comm;
7596 /* Turn on the DLB limiting (might have been on already) */
7597 comm->bPMELoadBalDLBLimits = TRUE;
7599 /* Change the cut-off limit */
7600 comm->PMELoadBal_max_cutoff = comm->cutoff;
7603 static void merge_cg_buffers(int ncell,
7604 gmx_domdec_comm_dim_t *cd, int pulse,
7606 int *index_gl, int *recv_i,
7607 rvec *cg_cm, rvec *recv_vr,
7609 cginfo_mb_t *cginfo_mb, int *cginfo)
7611 gmx_domdec_ind_t *ind, *ind_p;
7612 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7613 int shift, shift_at;
7615 ind = &cd->ind[pulse];
7617 /* First correct the already stored data */
7618 shift = ind->nrecv[ncell];
7619 for (cell = ncell-1; cell >= 0; cell--)
7621 shift -= ind->nrecv[cell];
7624 /* Move the cg's present from previous grid pulses */
7625 cg0 = ncg_cell[ncell+cell];
7626 cg1 = ncg_cell[ncell+cell+1];
7627 cgindex[cg1+shift] = cgindex[cg1];
7628 for (cg = cg1-1; cg >= cg0; cg--)
7630 index_gl[cg+shift] = index_gl[cg];
7631 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7632 cgindex[cg+shift] = cgindex[cg];
7633 cginfo[cg+shift] = cginfo[cg];
7635 /* Correct the already stored send indices for the shift */
7636 for (p = 1; p <= pulse; p++)
7638 ind_p = &cd->ind[p];
7640 for (c = 0; c < cell; c++)
7642 cg0 += ind_p->nsend[c];
7644 cg1 = cg0 + ind_p->nsend[cell];
7645 for (cg = cg0; cg < cg1; cg++)
7647 ind_p->index[cg] += shift;
7653 /* Merge in the communicated buffers */
7657 for (cell = 0; cell < ncell; cell++)
7659 cg1 = ncg_cell[ncell+cell+1] + shift;
7662 /* Correct the old cg indices */
7663 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7665 cgindex[cg+1] += shift_at;
7668 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7670 /* Copy this charge group from the buffer */
7671 index_gl[cg1] = recv_i[cg0];
7672 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7673 /* Add it to the cgindex */
7674 cg_gl = index_gl[cg1];
7675 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7676 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7677 cgindex[cg1+1] = cgindex[cg1] + nat;
7682 shift += ind->nrecv[cell];
7683 ncg_cell[ncell+cell+1] = cg1;
7687 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7688 int nzone, int cg0, const int *cgindex)
7692 /* Store the atom block boundaries for easy copying of communication buffers
7695 for (zone = 0; zone < nzone; zone++)
7697 for (p = 0; p < cd->np; p++)
7699 cd->ind[p].cell2at0[zone] = cgindex[cg];
7700 cg += cd->ind[p].nrecv[zone];
7701 cd->ind[p].cell2at1[zone] = cgindex[cg];
7706 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7712 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7714 if (!bLocalCG[link->a[i]])
7723 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7725 real c[DIM][4]; /* the corners for the non-bonded communication */
7726 real cr0; /* corner for rounding */
7727 real cr1[4]; /* corners for rounding */
7728 real bc[DIM]; /* corners for bounded communication */
7729 real bcr1; /* corner for rounding for bonded communication */
7732 /* Determine the corners of the domain(s) we are communicating with */
7734 set_dd_corners(const gmx_domdec_t *dd,
7735 int dim0, int dim1, int dim2,
7739 const gmx_domdec_comm_t *comm;
7740 const gmx_domdec_zones_t *zones;
7745 zones = &comm->zones;
7747 /* Keep the compiler happy */
7751 /* The first dimension is equal for all cells */
7752 c->c[0][0] = comm->cell_x0[dim0];
7755 c->bc[0] = c->c[0][0];
7760 /* This cell row is only seen from the first row */
7761 c->c[1][0] = comm->cell_x0[dim1];
7762 /* All rows can see this row */
7763 c->c[1][1] = comm->cell_x0[dim1];
7766 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7769 /* For the multi-body distance we need the maximum */
7770 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7773 /* Set the upper-right corner for rounding */
7774 c->cr0 = comm->cell_x1[dim0];
7779 for (j = 0; j < 4; j++)
7781 c->c[2][j] = comm->cell_x0[dim2];
7785 /* Use the maximum of the i-cells that see a j-cell */
7786 for (i = 0; i < zones->nizone; i++)
7788 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7794 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7800 /* For the multi-body distance we need the maximum */
7801 c->bc[2] = comm->cell_x0[dim2];
7802 for (i = 0; i < 2; i++)
7804 for (j = 0; j < 2; j++)
7806 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7812 /* Set the upper-right corner for rounding */
7813 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7814 * Only cell (0,0,0) can see cell 7 (1,1,1)
7816 c->cr1[0] = comm->cell_x1[dim1];
7817 c->cr1[3] = comm->cell_x1[dim1];
7820 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7823 /* For the multi-body distance we need the maximum */
7824 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7831 /* Determine which cg's we need to send in this pulse from this zone */
7833 get_zone_pulse_cgs(gmx_domdec_t *dd,
7834 int zonei, int zone,
7836 const int *index_gl,
7838 int dim, int dim_ind,
7839 int dim0, int dim1, int dim2,
7840 real r_comm2, real r_bcomm2,
7844 real skew_fac2_d, real skew_fac_01,
7845 rvec *v_d, rvec *v_0, rvec *v_1,
7846 const dd_corners_t *c,
7848 gmx_bool bDistBonded,
7854 gmx_domdec_ind_t *ind,
7855 int **ibuf, int *ibuf_nalloc,
7861 gmx_domdec_comm_t *comm;
7863 gmx_bool bDistMB_pulse;
7865 real r2, rb2, r, tric_sh;
7868 int nsend_z, nsend, nat;
7872 bScrew = (dd->bScrewPBC && dim == XX);
7874 bDistMB_pulse = (bDistMB && bDistBonded);
7880 for (cg = cg0; cg < cg1; cg++)
7884 if (tric_dist[dim_ind] == 0)
7886 /* Rectangular direction, easy */
7887 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7894 r = cg_cm[cg][dim] - c->bc[dim_ind];
7900 /* Rounding gives at most a 16% reduction
7901 * in communicated atoms
7903 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7905 r = cg_cm[cg][dim0] - c->cr0;
7906 /* This is the first dimension, so always r >= 0 */
7913 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7915 r = cg_cm[cg][dim1] - c->cr1[zone];
7922 r = cg_cm[cg][dim1] - c->bcr1;
7932 /* Triclinic direction, more complicated */
7935 /* Rounding, conservative as the skew_fac multiplication
7936 * will slightly underestimate the distance.
7938 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7940 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7941 for (i = dim0+1; i < DIM; i++)
7943 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7945 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7948 rb[dim0] = rn[dim0];
7951 /* Take care that the cell planes along dim0 might not
7952 * be orthogonal to those along dim1 and dim2.
7954 for (i = 1; i <= dim_ind; i++)
7957 if (normal[dim0][dimd] > 0)
7959 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7962 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7967 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7969 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7971 for (i = dim1+1; i < DIM; i++)
7973 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7975 rn[dim1] += tric_sh;
7978 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7979 /* Take care of coupling of the distances
7980 * to the planes along dim0 and dim1 through dim2.
7982 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7983 /* Take care that the cell planes along dim1
7984 * might not be orthogonal to that along dim2.
7986 if (normal[dim1][dim2] > 0)
7988 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7994 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7997 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7998 /* Take care of coupling of the distances
7999 * to the planes along dim0 and dim1 through dim2.
8001 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8002 /* Take care that the cell planes along dim1
8003 * might not be orthogonal to that along dim2.
8005 if (normal[dim1][dim2] > 0)
8007 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8012 /* The distance along the communication direction */
8013 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8015 for (i = dim+1; i < DIM; i++)
8017 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8022 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8023 /* Take care of coupling of the distances
8024 * to the planes along dim0 and dim1 through dim2.
8026 if (dim_ind == 1 && zonei == 1)
8028 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8034 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8037 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8038 /* Take care of coupling of the distances
8039 * to the planes along dim0 and dim1 through dim2.
8041 if (dim_ind == 1 && zonei == 1)
8043 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8051 ((bDistMB && rb2 < r_bcomm2) ||
8052 (bDist2B && r2 < r_bcomm2)) &&
8054 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8055 missing_link(comm->cglink, index_gl[cg],
8058 /* Make an index to the local charge groups */
8059 if (nsend+1 > ind->nalloc)
8061 ind->nalloc = over_alloc_large(nsend+1);
8062 srenew(ind->index, ind->nalloc);
8064 if (nsend+1 > *ibuf_nalloc)
8066 *ibuf_nalloc = over_alloc_large(nsend+1);
8067 srenew(*ibuf, *ibuf_nalloc);
8069 ind->index[nsend] = cg;
8070 (*ibuf)[nsend] = index_gl[cg];
8072 vec_rvec_check_alloc(vbuf, nsend+1);
8074 if (dd->ci[dim] == 0)
8076 /* Correct cg_cm for pbc */
8077 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8080 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8081 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8086 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8089 nat += cgindex[cg+1] - cgindex[cg];
8095 *nsend_z_ptr = nsend_z;
8098 static void setup_dd_communication(gmx_domdec_t *dd,
8099 matrix box, gmx_ddbox_t *ddbox,
8100 t_forcerec *fr, t_state *state, rvec **f)
8102 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8103 int nzone, nzone_send, zone, zonei, cg0, cg1;
8104 int c, i, j, cg, cg_gl, nrcg;
8105 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8106 gmx_domdec_comm_t *comm;
8107 gmx_domdec_zones_t *zones;
8108 gmx_domdec_comm_dim_t *cd;
8109 gmx_domdec_ind_t *ind;
8110 cginfo_mb_t *cginfo_mb;
8111 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8112 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
8113 dd_corners_t corners;
8115 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8116 real skew_fac2_d, skew_fac_01;
8123 fprintf(debug, "Setting up DD communication\n");
8128 switch (fr->cutoff_scheme)
8137 gmx_incons("unimplemented");
8141 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8143 dim = dd->dim[dim_ind];
8145 /* Check if we need to use triclinic distances */
8146 tric_dist[dim_ind] = 0;
8147 for (i = 0; i <= dim_ind; i++)
8149 if (ddbox->tric_dir[dd->dim[i]])
8151 tric_dist[dim_ind] = 1;
8156 bBondComm = comm->bBondComm;
8158 /* Do we need to determine extra distances for multi-body bondeds? */
8159 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8161 /* Do we need to determine extra distances for only two-body bondeds? */
8162 bDist2B = (bBondComm && !bDistMB);
8164 r_comm2 = sqr(comm->cutoff);
8165 r_bcomm2 = sqr(comm->cutoff_mbody);
8169 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8172 zones = &comm->zones;
8175 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8176 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8178 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8180 /* Triclinic stuff */
8181 normal = ddbox->normal;
8185 v_0 = ddbox->v[dim0];
8186 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8188 /* Determine the coupling coefficient for the distances
8189 * to the cell planes along dim0 and dim1 through dim2.
8190 * This is required for correct rounding.
8193 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8196 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8202 v_1 = ddbox->v[dim1];
8205 zone_cg_range = zones->cg_range;
8206 index_gl = dd->index_gl;
8207 cgindex = dd->cgindex;
8208 cginfo_mb = fr->cginfo_mb;
8210 zone_cg_range[0] = 0;
8211 zone_cg_range[1] = dd->ncg_home;
8212 comm->zone_ncg1[0] = dd->ncg_home;
8213 pos_cg = dd->ncg_home;
8215 nat_tot = dd->nat_home;
8217 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8219 dim = dd->dim[dim_ind];
8220 cd = &comm->cd[dim_ind];
8222 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8224 /* No pbc in this dimension, the first node should not comm. */
8232 v_d = ddbox->v[dim];
8233 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8235 cd->bInPlace = TRUE;
8236 for (p = 0; p < cd->np; p++)
8238 /* Only atoms communicated in the first pulse are used
8239 * for multi-body bonded interactions or for bBondComm.
8241 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8246 for (zone = 0; zone < nzone_send; zone++)
8248 if (tric_dist[dim_ind] && dim_ind > 0)
8250 /* Determine slightly more optimized skew_fac's
8252 * This reduces the number of communicated atoms
8253 * by about 10% for 3D DD of rhombic dodecahedra.
8255 for (dimd = 0; dimd < dim; dimd++)
8257 sf2_round[dimd] = 1;
8258 if (ddbox->tric_dir[dimd])
8260 for (i = dd->dim[dimd]+1; i < DIM; i++)
8262 /* If we are shifted in dimension i
8263 * and the cell plane is tilted forward
8264 * in dimension i, skip this coupling.
8266 if (!(zones->shift[nzone+zone][i] &&
8267 ddbox->v[dimd][i][dimd] >= 0))
8270 sqr(ddbox->v[dimd][i][dimd]);
8273 sf2_round[dimd] = 1/sf2_round[dimd];
8278 zonei = zone_perm[dim_ind][zone];
8281 /* Here we permutate the zones to obtain a convenient order
8282 * for neighbor searching
8284 cg0 = zone_cg_range[zonei];
8285 cg1 = zone_cg_range[zonei+1];
8289 /* Look only at the cg's received in the previous grid pulse
8291 cg1 = zone_cg_range[nzone+zone+1];
8292 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8295 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8296 for (th = 0; th < comm->nth; th++)
8298 gmx_domdec_ind_t *ind_p;
8299 int **ibuf_p, *ibuf_nalloc_p;
8301 int *nsend_p, *nat_p;
8307 /* Thread 0 writes in the comm buffers */
8309 ibuf_p = &comm->buf_int;
8310 ibuf_nalloc_p = &comm->nalloc_int;
8311 vbuf_p = &comm->vbuf;
8314 nsend_zone_p = &ind->nsend[zone];
8318 /* Other threads write into temp buffers */
8319 ind_p = &comm->dth[th].ind;
8320 ibuf_p = &comm->dth[th].ibuf;
8321 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8322 vbuf_p = &comm->dth[th].vbuf;
8323 nsend_p = &comm->dth[th].nsend;
8324 nat_p = &comm->dth[th].nat;
8325 nsend_zone_p = &comm->dth[th].nsend_zone;
8327 comm->dth[th].nsend = 0;
8328 comm->dth[th].nat = 0;
8329 comm->dth[th].nsend_zone = 0;
8339 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8340 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8343 /* Get the cg's for this pulse in this zone */
8344 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8346 dim, dim_ind, dim0, dim1, dim2,
8349 normal, skew_fac2_d, skew_fac_01,
8350 v_d, v_0, v_1, &corners, sf2_round,
8351 bDistBonded, bBondComm,
8355 ibuf_p, ibuf_nalloc_p,
8361 /* Append data of threads>=1 to the communication buffers */
8362 for (th = 1; th < comm->nth; th++)
8364 dd_comm_setup_work_t *dth;
8367 dth = &comm->dth[th];
8369 ns1 = nsend + dth->nsend_zone;
8370 if (ns1 > ind->nalloc)
8372 ind->nalloc = over_alloc_dd(ns1);
8373 srenew(ind->index, ind->nalloc);
8375 if (ns1 > comm->nalloc_int)
8377 comm->nalloc_int = over_alloc_dd(ns1);
8378 srenew(comm->buf_int, comm->nalloc_int);
8380 if (ns1 > comm->vbuf.nalloc)
8382 comm->vbuf.nalloc = over_alloc_dd(ns1);
8383 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8386 for (i = 0; i < dth->nsend_zone; i++)
8388 ind->index[nsend] = dth->ind.index[i];
8389 comm->buf_int[nsend] = dth->ibuf[i];
8390 copy_rvec(dth->vbuf.v[i],
8391 comm->vbuf.v[nsend]);
8395 ind->nsend[zone] += dth->nsend_zone;
8398 /* Clear the counts in case we do not have pbc */
8399 for (zone = nzone_send; zone < nzone; zone++)
8401 ind->nsend[zone] = 0;
8403 ind->nsend[nzone] = nsend;
8404 ind->nsend[nzone+1] = nat;
8405 /* Communicate the number of cg's and atoms to receive */
8406 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8407 ind->nsend, nzone+2,
8408 ind->nrecv, nzone+2);
8410 /* The rvec buffer is also required for atom buffers of size nsend
8411 * in dd_move_x and dd_move_f.
8413 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8417 /* We can receive in place if only the last zone is not empty */
8418 for (zone = 0; zone < nzone-1; zone++)
8420 if (ind->nrecv[zone] > 0)
8422 cd->bInPlace = FALSE;
8427 /* The int buffer is only required here for the cg indices */
8428 if (ind->nrecv[nzone] > comm->nalloc_int2)
8430 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8431 srenew(comm->buf_int2, comm->nalloc_int2);
8433 /* The rvec buffer is also required for atom buffers
8434 * of size nrecv in dd_move_x and dd_move_f.
8436 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8437 vec_rvec_check_alloc(&comm->vbuf2, i);
8441 /* Make space for the global cg indices */
8442 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8443 || dd->cg_nalloc == 0)
8445 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8446 srenew(index_gl, dd->cg_nalloc);
8447 srenew(cgindex, dd->cg_nalloc+1);
8449 /* Communicate the global cg indices */
8452 recv_i = index_gl + pos_cg;
8456 recv_i = comm->buf_int2;
8458 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8459 comm->buf_int, nsend,
8460 recv_i, ind->nrecv[nzone]);
8462 /* Make space for cg_cm */
8463 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8464 if (fr->cutoff_scheme == ecutsGROUP)
8472 /* Communicate cg_cm */
8475 recv_vr = cg_cm + pos_cg;
8479 recv_vr = comm->vbuf2.v;
8481 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8482 comm->vbuf.v, nsend,
8483 recv_vr, ind->nrecv[nzone]);
8485 /* Make the charge group index */
8488 zone = (p == 0 ? 0 : nzone - 1);
8489 while (zone < nzone)
8491 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8493 cg_gl = index_gl[pos_cg];
8494 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8495 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8496 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8499 /* Update the charge group presence,
8500 * so we can use it in the next pass of the loop.
8502 comm->bLocalCG[cg_gl] = TRUE;
8508 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8511 zone_cg_range[nzone+zone] = pos_cg;
8516 /* This part of the code is never executed with bBondComm. */
8517 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8518 index_gl, recv_i, cg_cm, recv_vr,
8519 cgindex, fr->cginfo_mb, fr->cginfo);
8520 pos_cg += ind->nrecv[nzone];
8522 nat_tot += ind->nrecv[nzone+1];
8526 /* Store the atom block for easy copying of communication buffers */
8527 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8531 dd->index_gl = index_gl;
8532 dd->cgindex = cgindex;
8534 dd->ncg_tot = zone_cg_range[zones->n];
8535 dd->nat_tot = nat_tot;
8536 comm->nat[ddnatHOME] = dd->nat_home;
8537 for (i = ddnatZONE; i < ddnatNR; i++)
8539 comm->nat[i] = dd->nat_tot;
8544 /* We don't need to update cginfo, since that was alrady done above.
8545 * So we pass NULL for the forcerec.
8547 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8548 NULL, comm->bLocalCG);
8553 fprintf(debug, "Finished setting up DD communication, zones:");
8554 for (c = 0; c < zones->n; c++)
8556 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8558 fprintf(debug, "\n");
8562 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8566 for (c = 0; c < zones->nizone; c++)
8568 zones->izone[c].cg1 = zones->cg_range[c+1];
8569 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8570 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8574 static void set_zones_size(gmx_domdec_t *dd,
8575 matrix box, const gmx_ddbox_t *ddbox,
8576 int zone_start, int zone_end)
8578 gmx_domdec_comm_t *comm;
8579 gmx_domdec_zones_t *zones;
8581 int z, zi, zj0, zj1, d, dim;
8584 real size_j, add_tric;
8589 zones = &comm->zones;
8591 /* Do we need to determine extra distances for multi-body bondeds? */
8592 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8594 for (z = zone_start; z < zone_end; z++)
8596 /* Copy cell limits to zone limits.
8597 * Valid for non-DD dims and non-shifted dims.
8599 copy_rvec(comm->cell_x0, zones->size[z].x0);
8600 copy_rvec(comm->cell_x1, zones->size[z].x1);
8603 for (d = 0; d < dd->ndim; d++)
8607 for (z = 0; z < zones->n; z++)
8609 /* With a staggered grid we have different sizes
8610 * for non-shifted dimensions.
8612 if (dd->bGridJump && zones->shift[z][dim] == 0)
8616 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8617 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8621 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8622 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8628 rcmbs = comm->cutoff_mbody;
8629 if (ddbox->tric_dir[dim])
8631 rcs /= ddbox->skew_fac[dim];
8632 rcmbs /= ddbox->skew_fac[dim];
8635 /* Set the lower limit for the shifted zone dimensions */
8636 for (z = zone_start; z < zone_end; z++)
8638 if (zones->shift[z][dim] > 0)
8641 if (!dd->bGridJump || d == 0)
8643 zones->size[z].x0[dim] = comm->cell_x1[dim];
8644 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8648 /* Here we take the lower limit of the zone from
8649 * the lowest domain of the zone below.
8653 zones->size[z].x0[dim] =
8654 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8660 zones->size[z].x0[dim] =
8661 zones->size[zone_perm[2][z-4]].x0[dim];
8665 zones->size[z].x0[dim] =
8666 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8669 /* A temporary limit, is updated below */
8670 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8674 for (zi = 0; zi < zones->nizone; zi++)
8676 if (zones->shift[zi][dim] == 0)
8678 /* This takes the whole zone into account.
8679 * With multiple pulses this will lead
8680 * to a larger zone then strictly necessary.
8682 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8683 zones->size[zi].x1[dim]+rcmbs);
8691 /* Loop over the i-zones to set the upper limit of each
8694 for (zi = 0; zi < zones->nizone; zi++)
8696 if (zones->shift[zi][dim] == 0)
8698 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8700 if (zones->shift[z][dim] > 0)
8702 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8703 zones->size[zi].x1[dim]+rcs);
8710 for (z = zone_start; z < zone_end; z++)
8712 /* Initialization only required to keep the compiler happy */
8713 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8716 /* To determine the bounding box for a zone we need to find
8717 * the extreme corners of 4, 2 or 1 corners.
8719 nc = 1 << (ddbox->npbcdim - 1);
8721 for (c = 0; c < nc; c++)
8723 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8727 corner[YY] = zones->size[z].x0[YY];
8731 corner[YY] = zones->size[z].x1[YY];
8735 corner[ZZ] = zones->size[z].x0[ZZ];
8739 corner[ZZ] = zones->size[z].x1[ZZ];
8741 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8743 /* With 1D domain decomposition the cg's are not in
8744 * the triclinic box, but triclinic x-y and rectangular y-z.
8745 * Shift y back, so it will later end up at 0.
8747 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8749 /* Apply the triclinic couplings */
8750 for (i = YY; i < ddbox->npbcdim; i++)
8752 for (j = XX; j < i; j++)
8754 corner[j] += corner[i]*box[i][j]/box[i][i];
8759 copy_rvec(corner, corner_min);
8760 copy_rvec(corner, corner_max);
8764 for (i = 0; i < DIM; i++)
8766 corner_min[i] = min(corner_min[i], corner[i]);
8767 corner_max[i] = max(corner_max[i], corner[i]);
8771 /* Copy the extreme cornes without offset along x */
8772 for (i = 0; i < DIM; i++)
8774 zones->size[z].bb_x0[i] = corner_min[i];
8775 zones->size[z].bb_x1[i] = corner_max[i];
8777 /* Add the offset along x */
8778 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8779 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8782 if (zone_start == 0)
8785 for (dim = 0; dim < DIM; dim++)
8787 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8789 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8794 for (z = zone_start; z < zone_end; z++)
8796 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8798 zones->size[z].x0[XX], zones->size[z].x1[XX],
8799 zones->size[z].x0[YY], zones->size[z].x1[YY],
8800 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8801 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8803 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8804 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8805 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8810 static int comp_cgsort(const void *a, const void *b)
8814 gmx_cgsort_t *cga, *cgb;
8815 cga = (gmx_cgsort_t *)a;
8816 cgb = (gmx_cgsort_t *)b;
8818 comp = cga->nsc - cgb->nsc;
8821 comp = cga->ind_gl - cgb->ind_gl;
8827 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8832 /* Order the data */
8833 for (i = 0; i < n; i++)
8835 buf[i] = a[sort[i].ind];
8838 /* Copy back to the original array */
8839 for (i = 0; i < n; i++)
8845 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8850 /* Order the data */
8851 for (i = 0; i < n; i++)
8853 copy_rvec(v[sort[i].ind], buf[i]);
8856 /* Copy back to the original array */
8857 for (i = 0; i < n; i++)
8859 copy_rvec(buf[i], v[i]);
8863 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8866 int a, atot, cg, cg0, cg1, i;
8868 if (cgindex == NULL)
8870 /* Avoid the useless loop of the atoms within a cg */
8871 order_vec_cg(ncg, sort, v, buf);
8876 /* Order the data */
8878 for (cg = 0; cg < ncg; cg++)
8880 cg0 = cgindex[sort[cg].ind];
8881 cg1 = cgindex[sort[cg].ind+1];
8882 for (i = cg0; i < cg1; i++)
8884 copy_rvec(v[i], buf[a]);
8890 /* Copy back to the original array */
8891 for (a = 0; a < atot; a++)
8893 copy_rvec(buf[a], v[a]);
8897 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8898 int nsort_new, gmx_cgsort_t *sort_new,
8899 gmx_cgsort_t *sort1)
8903 /* The new indices are not very ordered, so we qsort them */
8904 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8906 /* sort2 is already ordered, so now we can merge the two arrays */
8910 while (i2 < nsort2 || i_new < nsort_new)
8914 sort1[i1++] = sort_new[i_new++];
8916 else if (i_new == nsort_new)
8918 sort1[i1++] = sort2[i2++];
8920 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8921 (sort2[i2].nsc == sort_new[i_new].nsc &&
8922 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8924 sort1[i1++] = sort2[i2++];
8928 sort1[i1++] = sort_new[i_new++];
8933 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8935 gmx_domdec_sort_t *sort;
8936 gmx_cgsort_t *cgsort, *sort_i;
8937 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8938 int sort_last, sort_skip;
8940 sort = dd->comm->sort;
8942 a = fr->ns.grid->cell_index;
8944 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8946 if (ncg_home_old >= 0)
8948 /* The charge groups that remained in the same ns grid cell
8949 * are completely ordered. So we can sort efficiently by sorting
8950 * the charge groups that did move into the stationary list.
8955 for (i = 0; i < dd->ncg_home; i++)
8957 /* Check if this cg did not move to another node */
8960 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8962 /* This cg is new on this node or moved ns grid cell */
8963 if (nsort_new >= sort->sort_new_nalloc)
8965 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8966 srenew(sort->sort_new, sort->sort_new_nalloc);
8968 sort_i = &(sort->sort_new[nsort_new++]);
8972 /* This cg did not move */
8973 sort_i = &(sort->sort2[nsort2++]);
8975 /* Sort on the ns grid cell indices
8976 * and the global topology index.
8977 * index_gl is irrelevant with cell ns,
8978 * but we set it here anyhow to avoid a conditional.
8981 sort_i->ind_gl = dd->index_gl[i];
8988 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8991 /* Sort efficiently */
8992 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8997 cgsort = sort->sort;
8999 for (i = 0; i < dd->ncg_home; i++)
9001 /* Sort on the ns grid cell indices
9002 * and the global topology index
9004 cgsort[i].nsc = a[i];
9005 cgsort[i].ind_gl = dd->index_gl[i];
9007 if (cgsort[i].nsc < moved)
9014 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9016 /* Determine the order of the charge groups using qsort */
9017 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9023 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9026 int ncg_new, i, *a, na;
9028 sort = dd->comm->sort->sort;
9030 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9033 for (i = 0; i < na; i++)
9037 sort[ncg_new].ind = a[i];
9045 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
9046 rvec *cgcm, t_forcerec *fr, t_state *state,
9049 gmx_domdec_sort_t *sort;
9050 gmx_cgsort_t *cgsort, *sort_i;
9052 int ncg_new, i, *ibuf, cgsize;
9055 sort = dd->comm->sort;
9057 if (dd->ncg_home > sort->sort_nalloc)
9059 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9060 srenew(sort->sort, sort->sort_nalloc);
9061 srenew(sort->sort2, sort->sort_nalloc);
9063 cgsort = sort->sort;
9065 switch (fr->cutoff_scheme)
9068 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9071 ncg_new = dd_sort_order_nbnxn(dd, fr);
9074 gmx_incons("unimplemented");
9078 /* We alloc with the old size, since cgindex is still old */
9079 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9080 vbuf = dd->comm->vbuf.v;
9084 cgindex = dd->cgindex;
9091 /* Remove the charge groups which are no longer at home here */
9092 dd->ncg_home = ncg_new;
9095 fprintf(debug, "Set the new home charge group count to %d\n",
9099 /* Reorder the state */
9100 for (i = 0; i < estNR; i++)
9102 if (EST_DISTR(i) && (state->flags & (1<<i)))
9107 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9110 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9113 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9116 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9120 case estDISRE_INITF:
9121 case estDISRE_RM3TAV:
9122 case estORIRE_INITF:
9124 /* No ordering required */
9127 gmx_incons("Unknown state entry encountered in dd_sort_state");
9132 if (fr->cutoff_scheme == ecutsGROUP)
9135 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9138 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9140 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9141 srenew(sort->ibuf, sort->ibuf_nalloc);
9144 /* Reorder the global cg index */
9145 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9146 /* Reorder the cginfo */
9147 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9148 /* Rebuild the local cg index */
9152 for (i = 0; i < dd->ncg_home; i++)
9154 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9155 ibuf[i+1] = ibuf[i] + cgsize;
9157 for (i = 0; i < dd->ncg_home+1; i++)
9159 dd->cgindex[i] = ibuf[i];
9164 for (i = 0; i < dd->ncg_home+1; i++)
9169 /* Set the home atom number */
9170 dd->nat_home = dd->cgindex[dd->ncg_home];
9172 if (fr->cutoff_scheme == ecutsVERLET)
9174 /* The atoms are now exactly in grid order, update the grid order */
9175 nbnxn_set_atomorder(fr->nbv->nbs);
9179 /* Copy the sorted ns cell indices back to the ns grid struct */
9180 for (i = 0; i < dd->ncg_home; i++)
9182 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9184 fr->ns.grid->nr = dd->ncg_home;
9188 static void add_dd_statistics(gmx_domdec_t *dd)
9190 gmx_domdec_comm_t *comm;
9195 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9197 comm->sum_nat[ddnat-ddnatZONE] +=
9198 comm->nat[ddnat] - comm->nat[ddnat-1];
9203 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9205 gmx_domdec_comm_t *comm;
9210 /* Reset all the statistics and counters for total run counting */
9211 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9213 comm->sum_nat[ddnat-ddnatZONE] = 0;
9217 comm->load_step = 0;
9220 clear_ivec(comm->load_lim);
9225 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9227 gmx_domdec_comm_t *comm;
9231 comm = cr->dd->comm;
9233 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9240 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");
9242 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9244 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9249 " av. #atoms communicated per step for force: %d x %.1f\n",
9253 if (cr->dd->vsite_comm)
9256 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9257 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9262 if (cr->dd->constraint_comm)
9265 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9266 1 + ir->nLincsIter, av);
9270 gmx_incons(" Unknown type for DD statistics");
9273 fprintf(fplog, "\n");
9275 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9277 print_dd_load_av(fplog, cr->dd);
9281 void dd_partition_system(FILE *fplog,
9282 gmx_large_int_t step,
9284 gmx_bool bMasterState,
9286 t_state *state_global,
9287 gmx_mtop_t *top_global,
9289 t_state *state_local,
9292 gmx_localtop_t *top_local,
9295 gmx_shellfc_t shellfc,
9296 gmx_constr_t constr,
9298 gmx_wallcycle_t wcycle,
9302 gmx_domdec_comm_t *comm;
9303 gmx_ddbox_t ddbox = {0};
9305 gmx_large_int_t step_pcoupl;
9306 rvec cell_ns_x0, cell_ns_x1;
9307 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9308 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9309 gmx_bool bRedist, bSortCG, bResortAll;
9310 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9317 bBoxChanged = (bMasterState || DEFORM(*ir));
9318 if (ir->epc != epcNO)
9320 /* With nstpcouple > 1 pressure coupling happens.
9321 * one step after calculating the pressure.
9322 * Box scaling happens at the end of the MD step,
9323 * after the DD partitioning.
9324 * We therefore have to do DLB in the first partitioning
9325 * after an MD step where P-coupling occured.
9326 * We need to determine the last step in which p-coupling occurred.
9327 * MRS -- need to validate this for vv?
9332 step_pcoupl = step - 1;
9336 step_pcoupl = ((step - 1)/n)*n + 1;
9338 if (step_pcoupl >= comm->partition_step)
9344 bNStGlobalComm = (step % nstglobalcomm == 0);
9346 if (!comm->bDynLoadBal)
9352 /* Should we do dynamic load balacing this step?
9353 * Since it requires (possibly expensive) global communication,
9354 * we might want to do DLB less frequently.
9356 if (bBoxChanged || ir->epc != epcNO)
9358 bDoDLB = bBoxChanged;
9362 bDoDLB = bNStGlobalComm;
9366 /* Check if we have recorded loads on the nodes */
9367 if (comm->bRecordLoad && dd_load_count(comm))
9369 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9371 /* Check if we should use DLB at the second partitioning
9372 * and every 100 partitionings,
9373 * so the extra communication cost is negligible.
9375 n = max(100, nstglobalcomm);
9376 bCheckDLB = (comm->n_load_collect == 0 ||
9377 comm->n_load_have % n == n-1);
9384 /* Print load every nstlog, first and last step to the log file */
9385 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9386 comm->n_load_collect == 0 ||
9388 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9390 /* Avoid extra communication due to verbose screen output
9391 * when nstglobalcomm is set.
9393 if (bDoDLB || bLogLoad || bCheckDLB ||
9394 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9396 get_load_distribution(dd, wcycle);
9401 dd_print_load(fplog, dd, step-1);
9405 dd_print_load_verbose(dd);
9408 comm->n_load_collect++;
9412 /* Since the timings are node dependent, the master decides */
9416 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9419 fprintf(debug, "step %s, imb loss %f\n",
9420 gmx_step_str(step, sbuf),
9421 dd_force_imb_perf_loss(dd));
9424 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9427 turn_on_dlb(fplog, cr, step);
9432 comm->n_load_have++;
9435 cgs_gl = &comm->cgs_gl;
9440 /* Clear the old state */
9441 clear_dd_indices(dd, 0, 0);
9444 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9445 TRUE, cgs_gl, state_global->x, &ddbox);
9447 get_cg_distribution(fplog, step, dd, cgs_gl,
9448 state_global->box, &ddbox, state_global->x);
9450 dd_distribute_state(dd, cgs_gl,
9451 state_global, state_local, f);
9453 dd_make_local_cgs(dd, &top_local->cgs);
9455 /* Ensure that we have space for the new distribution */
9456 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9458 if (fr->cutoff_scheme == ecutsGROUP)
9460 calc_cgcm(fplog, 0, dd->ncg_home,
9461 &top_local->cgs, state_local->x, fr->cg_cm);
9464 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9466 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9468 else if (state_local->ddp_count != dd->ddp_count)
9470 if (state_local->ddp_count > dd->ddp_count)
9472 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9475 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9477 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);
9480 /* Clear the old state */
9481 clear_dd_indices(dd, 0, 0);
9483 /* Build the new indices */
9484 rebuild_cgindex(dd, cgs_gl->index, state_local);
9485 make_dd_indices(dd, cgs_gl->index, 0);
9486 ncgindex_set = dd->ncg_home;
9488 if (fr->cutoff_scheme == ecutsGROUP)
9490 /* Redetermine the cg COMs */
9491 calc_cgcm(fplog, 0, dd->ncg_home,
9492 &top_local->cgs, state_local->x, fr->cg_cm);
9495 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9497 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9499 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9500 TRUE, &top_local->cgs, state_local->x, &ddbox);
9502 bRedist = comm->bDynLoadBal;
9506 /* We have the full state, only redistribute the cgs */
9508 /* Clear the non-home indices */
9509 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9512 /* Avoid global communication for dim's without pbc and -gcom */
9513 if (!bNStGlobalComm)
9515 copy_rvec(comm->box0, ddbox.box0 );
9516 copy_rvec(comm->box_size, ddbox.box_size);
9518 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9519 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9524 /* For dim's without pbc and -gcom */
9525 copy_rvec(ddbox.box0, comm->box0 );
9526 copy_rvec(ddbox.box_size, comm->box_size);
9528 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9531 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9533 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9536 /* Check if we should sort the charge groups */
9537 if (comm->nstSortCG > 0)
9539 bSortCG = (bMasterState ||
9540 (bRedist && (step % comm->nstSortCG == 0)));
9547 ncg_home_old = dd->ncg_home;
9552 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9554 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9555 state_local, f, fr, mdatoms,
9556 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9558 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9561 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9563 &comm->cell_x0, &comm->cell_x1,
9564 dd->ncg_home, fr->cg_cm,
9565 cell_ns_x0, cell_ns_x1, &grid_density);
9569 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9572 switch (fr->cutoff_scheme)
9575 copy_ivec(fr->ns.grid->n, ncells_old);
9576 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9577 state_local->box, cell_ns_x0, cell_ns_x1,
9578 fr->rlistlong, grid_density);
9581 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9584 gmx_incons("unimplemented");
9586 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9587 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9591 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9593 /* Sort the state on charge group position.
9594 * This enables exact restarts from this step.
9595 * It also improves performance by about 15% with larger numbers
9596 * of atoms per node.
9599 /* Fill the ns grid with the home cell,
9600 * so we can sort with the indices.
9602 set_zones_ncg_home(dd);
9604 switch (fr->cutoff_scheme)
9607 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9609 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9611 comm->zones.size[0].bb_x0,
9612 comm->zones.size[0].bb_x1,
9614 comm->zones.dens_zone0,
9617 ncg_moved, bRedist ? comm->moved : NULL,
9618 fr->nbv->grp[eintLocal].kernel_type,
9619 fr->nbv->grp[eintLocal].nbat);
9621 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9624 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9625 0, dd->ncg_home, fr->cg_cm);
9627 copy_ivec(fr->ns.grid->n, ncells_new);
9630 gmx_incons("unimplemented");
9633 bResortAll = bMasterState;
9635 /* Check if we can user the old order and ns grid cell indices
9636 * of the charge groups to sort the charge groups efficiently.
9638 if (ncells_new[XX] != ncells_old[XX] ||
9639 ncells_new[YY] != ncells_old[YY] ||
9640 ncells_new[ZZ] != ncells_old[ZZ])
9647 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9648 gmx_step_str(step, sbuf), dd->ncg_home);
9650 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9651 bResortAll ? -1 : ncg_home_old);
9652 /* Rebuild all the indices */
9653 ga2la_clear(dd->ga2la);
9656 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9659 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9661 /* Setup up the communication and communicate the coordinates */
9662 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9664 /* Set the indices */
9665 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9667 /* Set the charge group boundaries for neighbor searching */
9668 set_cg_boundaries(&comm->zones);
9670 if (fr->cutoff_scheme == ecutsVERLET)
9672 set_zones_size(dd, state_local->box, &ddbox,
9673 bSortCG ? 1 : 0, comm->zones.n);
9676 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9679 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9680 -1,state_local->x,state_local->box);
9683 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9685 /* Extract a local topology from the global topology */
9686 for (i = 0; i < dd->ndim; i++)
9688 np[dd->dim[i]] = comm->cd[i].np;
9690 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9691 comm->cellsize_min, np,
9693 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9694 vsite, top_global, top_local);
9696 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9698 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9700 /* Set up the special atom communication */
9701 n = comm->nat[ddnatZONE];
9702 for (i = ddnatZONE+1; i < ddnatNR; i++)
9707 if (vsite && vsite->n_intercg_vsite)
9709 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9713 if (dd->bInterCGcons || dd->bInterCGsettles)
9715 /* Only for inter-cg constraints we need special code */
9716 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9717 constr, ir->nProjOrder,
9718 top_local->idef.il);
9722 gmx_incons("Unknown special atom type setup");
9727 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9729 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9731 /* Make space for the extra coordinates for virtual site
9732 * or constraint communication.
9734 state_local->natoms = comm->nat[ddnatNR-1];
9735 if (state_local->natoms > state_local->nalloc)
9737 dd_realloc_state(state_local, f, state_local->natoms);
9740 if (fr->bF_NoVirSum)
9742 if (vsite && vsite->n_intercg_vsite)
9744 nat_f_novirsum = comm->nat[ddnatVSITE];
9748 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9750 nat_f_novirsum = dd->nat_tot;
9754 nat_f_novirsum = dd->nat_home;
9763 /* Set the number of atoms required for the force calculation.
9764 * Forces need to be constrained when using a twin-range setup
9765 * or with energy minimization. For simple simulations we could
9766 * avoid some allocation, zeroing and copying, but this is
9767 * probably not worth the complications ande checking.
9769 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9770 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9772 /* We make the all mdatoms up to nat_tot_con.
9773 * We could save some work by only setting invmass
9774 * between nat_tot and nat_tot_con.
9776 /* This call also sets the new number of home particles to dd->nat_home */
9777 atoms2md(top_global, ir,
9778 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9780 /* Now we have the charges we can sort the FE interactions */
9781 dd_sort_local_top(dd, mdatoms, top_local);
9785 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9786 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9791 /* Make the local shell stuff, currently no communication is done */
9792 make_local_shells(cr, mdatoms, shellfc);
9795 if (ir->implicit_solvent)
9797 make_local_gb(cr, fr->born, ir->gb_algorithm);
9800 setup_bonded_threading(fr, &top_local->idef);
9802 if (!(cr->duty & DUTY_PME))
9804 /* Send the charges to our PME only node */
9805 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9806 mdatoms->chargeA, mdatoms->chargeB,
9807 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9812 set_constraints(constr, top_local, ir, mdatoms, cr);
9815 if (ir->ePull != epullNO)
9817 /* Update the local pull groups */
9818 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9823 /* Update the local rotation groups */
9824 dd_make_local_rotation_groups(dd, ir->rot);
9828 add_dd_statistics(dd);
9830 /* Make sure we only count the cycles for this DD partitioning */
9831 clear_dd_cycle_counts(dd);
9833 /* Because the order of the atoms might have changed since
9834 * the last vsite construction, we need to communicate the constructing
9835 * atom coordinates again (for spreading the forces this MD step).
9837 dd_move_x_vsites(dd, state_local->box, state_local->x);
9839 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9841 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9843 dd_move_x(dd, state_local->box, state_local->x);
9844 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9845 -1, state_local->x, state_local->box);
9848 /* Store the partitioning step */
9849 comm->partition_step = step;
9851 /* Increase the DD partitioning counter */
9853 /* The state currently matches this DD partitioning count, store it */
9854 state_local->ddp_count = dd->ddp_count;
9857 /* The DD master node knows the complete cg distribution,
9858 * store the count so we can possibly skip the cg info communication.
9860 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9863 if (comm->DD_debug > 0)
9865 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9866 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9867 "after partitioning");