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"
84 #define DDRANK(dd, rank) (rank)
85 #define DDMASTERRANK(dd) (dd->masterrank)
87 typedef struct gmx_domdec_master
89 /* The cell boundaries */
91 /* The global charge group division */
92 int *ncg; /* Number of home charge groups for each node */
93 int *index; /* Index of nnodes+1 into cg */
94 int *cg; /* Global charge group index */
95 int *nat; /* Number of home atoms for each node. */
96 int *ibuf; /* Buffer for communication */
97 rvec *vbuf; /* Buffer for state scattering and gathering */
98 } gmx_domdec_master_t;
102 /* The numbers of charge groups to send and receive for each cell
103 * that requires communication, the last entry contains the total
104 * number of atoms that needs to be communicated.
106 int nsend[DD_MAXIZONE+2];
107 int nrecv[DD_MAXIZONE+2];
108 /* The charge groups to send */
111 /* The atom range for non-in-place communication */
112 int cell2at0[DD_MAXIZONE];
113 int cell2at1[DD_MAXIZONE];
118 int np; /* Number of grid pulses in this dimension */
119 int np_dlb; /* For dlb, for use with edlbAUTO */
120 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
122 gmx_bool bInPlace; /* Can we communicate in place? */
123 } gmx_domdec_comm_dim_t;
127 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
128 real *cell_f; /* State var.: cell boundaries, box relative */
129 real *old_cell_f; /* Temp. var.: old cell size */
130 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
131 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
132 real *bound_min; /* Temp. var.: lower limit for cell boundary */
133 real *bound_max; /* Temp. var.: upper limit for cell boundary */
134 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
135 real *buf_ncd; /* Temp. var. */
138 #define DD_NLOAD_MAX 9
140 /* Here floats are accurate enough, since these variables
141 * only influence the load balancing, not the actual MD results.
168 gmx_cgsort_t *sort_new;
180 /* This enum determines the order of the coordinates.
181 * ddnatHOME and ddnatZONE should be first and second,
182 * the others can be ordered as wanted.
185 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
189 edlbAUTO, edlbNO, edlbYES, edlbNR
191 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
195 int dim; /* The dimension */
196 gmx_bool dim_match; /* Tells if DD and PME dims match */
197 int nslab; /* The number of PME slabs in this dimension */
198 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
199 int *pp_min; /* The minimum pp node location, size nslab */
200 int *pp_max; /* The maximum pp node location,size nslab */
201 int maxshift; /* The maximum shift for coordinate redistribution in PME */
206 real min0; /* The minimum bottom of this zone */
207 real max1; /* The maximum top of this zone */
208 real min1; /* The minimum top of this zone */
209 real mch0; /* The maximum bottom communicaton height for this zone */
210 real mch1; /* The maximum top communicaton height for this zone */
211 real p1_0; /* The bottom value of the first cell in this zone */
212 real p1_1; /* The top value of the first cell in this zone */
217 gmx_domdec_ind_t ind;
224 } dd_comm_setup_work_t;
226 typedef struct gmx_domdec_comm
228 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
229 * unless stated otherwise.
232 /* The number of decomposition dimensions for PME, 0: no PME */
234 /* The number of nodes doing PME (PP/PME or only PME) */
238 /* The communication setup including the PME only nodes */
239 gmx_bool bCartesianPP_PME;
242 int *pmenodes; /* size npmenodes */
243 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
244 * but with bCartesianPP_PME */
245 gmx_ddpme_t ddpme[2];
247 /* The DD particle-particle nodes only */
248 gmx_bool bCartesianPP;
249 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
251 /* The global charge groups */
254 /* Should we sort the cgs */
256 gmx_domdec_sort_t *sort;
258 /* Are there charge groups? */
261 /* Are there bonded and multi-body interactions between charge groups? */
262 gmx_bool bInterCGBondeds;
263 gmx_bool bInterCGMultiBody;
265 /* Data for the optional bonded interaction atom communication range */
272 /* Are we actually using DLB? */
273 gmx_bool bDynLoadBal;
275 /* Cell sizes for static load balancing, first index cartesian */
278 /* The width of the communicated boundaries */
281 /* The minimum cell size (including triclinic correction) */
283 /* For dlb, for use with edlbAUTO */
284 rvec cellsize_min_dlb;
285 /* The lower limit for the DD cell size with DLB */
287 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
288 gmx_bool bVacDLBNoLimit;
290 /* With PME load balancing we set limits on DLB */
291 gmx_bool bPMELoadBalDLBLimits;
292 /* DLB needs to take into account that we want to allow this maximum
293 * cut-off (for PME load balancing), this could limit cell boundaries.
295 real PMELoadBal_max_cutoff;
297 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
299 /* box0 and box_size are required with dim's without pbc and -gcom */
303 /* The cell boundaries */
307 /* The old location of the cell boundaries, to check cg displacements */
311 /* The communication setup and charge group boundaries for the zones */
312 gmx_domdec_zones_t zones;
314 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
315 * cell boundaries of neighboring cells for dynamic load balancing.
317 gmx_ddzone_t zone_d1[2];
318 gmx_ddzone_t zone_d2[2][2];
320 /* The coordinate/force communication setup and indices */
321 gmx_domdec_comm_dim_t cd[DIM];
322 /* The maximum number of cells to communicate with in one dimension */
325 /* Which cg distribution is stored on the master node */
326 int master_cg_ddp_count;
328 /* The number of cg's received from the direct neighbors */
329 int zone_ncg1[DD_MAXZONE];
331 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
334 /* Array for signalling if atoms have moved to another domain */
338 /* Communication buffer for general use */
342 /* Communication buffer for general use */
345 /* Temporary storage for thread parallel communication setup */
347 dd_comm_setup_work_t *dth;
349 /* Communication buffers only used with multiple grid pulses */
354 /* Communication buffers for local redistribution */
356 int cggl_flag_nalloc[DIM*2];
358 int cgcm_state_nalloc[DIM*2];
360 /* Cell sizes for dynamic load balancing */
361 gmx_domdec_root_t **root;
365 real cell_f_max0[DIM];
366 real cell_f_min1[DIM];
368 /* Stuff for load communication */
369 gmx_bool bRecordLoad;
370 gmx_domdec_load_t *load;
372 MPI_Comm *mpi_comm_load;
375 /* Maximum DLB scaling per load balancing step in percent */
379 float cycl[ddCyclNr];
380 int cycl_n[ddCyclNr];
381 float cycl_max[ddCyclNr];
382 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
386 /* Have often have did we have load measurements */
388 /* Have often have we collected the load measurements */
392 double sum_nat[ddnatNR-ddnatZONE];
402 /* The last partition step */
403 gmx_large_int_t partition_step;
411 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
414 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
415 #define DD_FLAG_NRCG 65535
416 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
417 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
419 /* Zone permutation required to obtain consecutive charge groups
420 * for neighbor searching.
422 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
424 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
425 * components see only j zones with that component 0.
428 /* The DD zone order */
429 static const ivec dd_zo[DD_MAXZONE] =
430 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
435 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
440 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
445 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
447 /* Factors used to avoid problems due to rounding issues */
448 #define DD_CELL_MARGIN 1.0001
449 #define DD_CELL_MARGIN2 1.00005
450 /* Factor to account for pressure scaling during nstlist steps */
451 #define DD_PRES_SCALE_MARGIN 1.02
453 /* Allowed performance loss before we DLB or warn */
454 #define DD_PERF_LOSS 0.05
456 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
458 /* Use separate MPI send and receive commands
459 * when nnodes <= GMX_DD_NNODES_SENDRECV.
460 * This saves memory (and some copying for small nnodes).
461 * For high parallelization scatter and gather calls are used.
463 #define GMX_DD_NNODES_SENDRECV 4
467 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
469 static void index2xyz(ivec nc,int ind,ivec xyz)
471 xyz[XX] = ind % nc[XX];
472 xyz[YY] = (ind / nc[XX]) % nc[YY];
473 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
477 /* This order is required to minimize the coordinate communication in PME
478 * which uses decomposition in the x direction.
480 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
482 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
484 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
485 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
486 xyz[ZZ] = ind % nc[ZZ];
489 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
494 ddindex = dd_index(dd->nc, c);
495 if (dd->comm->bCartesianPP_PME)
497 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
499 else if (dd->comm->bCartesianPP)
502 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
513 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
515 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
518 int ddglatnr(gmx_domdec_t *dd, int i)
528 if (i >= dd->comm->nat[ddnatNR-1])
530 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
532 atnr = dd->gatindex[i] + 1;
538 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
540 return &dd->comm->cgs_gl;
543 static void vec_rvec_init(vec_rvec_t *v)
549 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
553 v->nalloc = over_alloc_dd(n);
554 srenew(v->v, v->nalloc);
558 void dd_store_state(gmx_domdec_t *dd, t_state *state)
562 if (state->ddp_count != dd->ddp_count)
564 gmx_incons("The state does not the domain decomposition state");
567 state->ncg_gl = dd->ncg_home;
568 if (state->ncg_gl > state->cg_gl_nalloc)
570 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
571 srenew(state->cg_gl, state->cg_gl_nalloc);
573 for (i = 0; i < state->ncg_gl; i++)
575 state->cg_gl[i] = dd->index_gl[i];
578 state->ddp_count_cg_gl = dd->ddp_count;
581 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
583 return &dd->comm->zones;
586 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
587 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
589 gmx_domdec_zones_t *zones;
592 zones = &dd->comm->zones;
595 while (icg >= zones->izone[izone].cg1)
604 else if (izone < zones->nizone)
606 *jcg0 = zones->izone[izone].jcg0;
610 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
611 icg, izone, zones->nizone);
614 *jcg1 = zones->izone[izone].jcg1;
616 for (d = 0; d < dd->ndim; d++)
619 shift0[dim] = zones->izone[izone].shift0[dim];
620 shift1[dim] = zones->izone[izone].shift1[dim];
621 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
623 /* A conservative approach, this can be optimized */
630 int dd_natoms_vsite(gmx_domdec_t *dd)
632 return dd->comm->nat[ddnatVSITE];
635 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
637 *at_start = dd->comm->nat[ddnatCON-1];
638 *at_end = dd->comm->nat[ddnatCON];
641 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
643 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
644 int *index, *cgindex;
645 gmx_domdec_comm_t *comm;
646 gmx_domdec_comm_dim_t *cd;
647 gmx_domdec_ind_t *ind;
648 rvec shift = {0, 0, 0}, *buf, *rbuf;
649 gmx_bool bPBC, bScrew;
653 cgindex = dd->cgindex;
658 nat_tot = dd->nat_home;
659 for (d = 0; d < dd->ndim; d++)
661 bPBC = (dd->ci[dd->dim[d]] == 0);
662 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
665 copy_rvec(box[dd->dim[d]], shift);
668 for (p = 0; p < cd->np; p++)
675 for (i = 0; i < ind->nsend[nzone]; i++)
677 at0 = cgindex[index[i]];
678 at1 = cgindex[index[i]+1];
679 for (j = at0; j < at1; j++)
681 copy_rvec(x[j], buf[n]);
688 for (i = 0; i < ind->nsend[nzone]; i++)
690 at0 = cgindex[index[i]];
691 at1 = cgindex[index[i]+1];
692 for (j = at0; j < at1; j++)
694 /* We need to shift the coordinates */
695 rvec_add(x[j], shift, buf[n]);
702 for (i = 0; i < ind->nsend[nzone]; i++)
704 at0 = cgindex[index[i]];
705 at1 = cgindex[index[i]+1];
706 for (j = at0; j < at1; j++)
709 buf[n][XX] = x[j][XX] + shift[XX];
711 * This operation requires a special shift force
712 * treatment, which is performed in calc_vir.
714 buf[n][YY] = box[YY][YY] - x[j][YY];
715 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
727 rbuf = comm->vbuf2.v;
729 /* Send and receive the coordinates */
730 dd_sendrecv_rvec(dd, d, dddirBackward,
731 buf, ind->nsend[nzone+1],
732 rbuf, ind->nrecv[nzone+1]);
736 for (zone = 0; zone < nzone; zone++)
738 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
740 copy_rvec(rbuf[j], x[i]);
745 nat_tot += ind->nrecv[nzone+1];
751 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
753 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
754 int *index, *cgindex;
755 gmx_domdec_comm_t *comm;
756 gmx_domdec_comm_dim_t *cd;
757 gmx_domdec_ind_t *ind;
761 gmx_bool bPBC, bScrew;
765 cgindex = dd->cgindex;
770 nzone = comm->zones.n/2;
771 nat_tot = dd->nat_tot;
772 for (d = dd->ndim-1; d >= 0; d--)
774 bPBC = (dd->ci[dd->dim[d]] == 0);
775 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
776 if (fshift == NULL && !bScrew)
780 /* Determine which shift vector we need */
786 for (p = cd->np-1; p >= 0; p--)
789 nat_tot -= ind->nrecv[nzone+1];
796 sbuf = comm->vbuf2.v;
798 for (zone = 0; zone < nzone; zone++)
800 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
802 copy_rvec(f[i], sbuf[j]);
807 /* Communicate the forces */
808 dd_sendrecv_rvec(dd, d, dddirForward,
809 sbuf, ind->nrecv[nzone+1],
810 buf, ind->nsend[nzone+1]);
812 /* Add the received forces */
816 for (i = 0; i < ind->nsend[nzone]; i++)
818 at0 = cgindex[index[i]];
819 at1 = cgindex[index[i]+1];
820 for (j = at0; j < at1; j++)
822 rvec_inc(f[j], buf[n]);
829 for (i = 0; i < ind->nsend[nzone]; i++)
831 at0 = cgindex[index[i]];
832 at1 = cgindex[index[i]+1];
833 for (j = at0; j < at1; j++)
835 rvec_inc(f[j], buf[n]);
836 /* Add this force to the shift force */
837 rvec_inc(fshift[is], buf[n]);
844 for (i = 0; i < ind->nsend[nzone]; i++)
846 at0 = cgindex[index[i]];
847 at1 = cgindex[index[i]+1];
848 for (j = at0; j < at1; j++)
850 /* Rotate the force */
851 f[j][XX] += buf[n][XX];
852 f[j][YY] -= buf[n][YY];
853 f[j][ZZ] -= buf[n][ZZ];
856 /* Add this force to the shift force */
857 rvec_inc(fshift[is], buf[n]);
868 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
870 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
871 int *index, *cgindex;
872 gmx_domdec_comm_t *comm;
873 gmx_domdec_comm_dim_t *cd;
874 gmx_domdec_ind_t *ind;
879 cgindex = dd->cgindex;
881 buf = &comm->vbuf.v[0][0];
884 nat_tot = dd->nat_home;
885 for (d = 0; d < dd->ndim; d++)
888 for (p = 0; p < cd->np; p++)
893 for (i = 0; i < ind->nsend[nzone]; i++)
895 at0 = cgindex[index[i]];
896 at1 = cgindex[index[i]+1];
897 for (j = at0; j < at1; j++)
910 rbuf = &comm->vbuf2.v[0][0];
912 /* Send and receive the coordinates */
913 dd_sendrecv_real(dd, d, dddirBackward,
914 buf, ind->nsend[nzone+1],
915 rbuf, ind->nrecv[nzone+1]);
919 for (zone = 0; zone < nzone; zone++)
921 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
928 nat_tot += ind->nrecv[nzone+1];
934 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
936 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
937 int *index, *cgindex;
938 gmx_domdec_comm_t *comm;
939 gmx_domdec_comm_dim_t *cd;
940 gmx_domdec_ind_t *ind;
945 cgindex = dd->cgindex;
947 buf = &comm->vbuf.v[0][0];
950 nzone = comm->zones.n/2;
951 nat_tot = dd->nat_tot;
952 for (d = dd->ndim-1; d >= 0; d--)
955 for (p = cd->np-1; p >= 0; p--)
958 nat_tot -= ind->nrecv[nzone+1];
965 sbuf = &comm->vbuf2.v[0][0];
967 for (zone = 0; zone < nzone; zone++)
969 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
976 /* Communicate the forces */
977 dd_sendrecv_real(dd, d, dddirForward,
978 sbuf, ind->nrecv[nzone+1],
979 buf, ind->nsend[nzone+1]);
981 /* Add the received forces */
983 for (i = 0; i < ind->nsend[nzone]; i++)
985 at0 = cgindex[index[i]];
986 at1 = cgindex[index[i]+1];
987 for (j = at0; j < at1; j++)
998 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1000 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",
1002 zone->min0, zone->max1,
1003 zone->mch0, zone->mch0,
1004 zone->p1_0, zone->p1_1);
1008 #define DDZONECOMM_MAXZONE 5
1009 #define DDZONECOMM_BUFSIZE 3
1011 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1012 int ddimind, int direction,
1013 gmx_ddzone_t *buf_s, int n_s,
1014 gmx_ddzone_t *buf_r, int n_r)
1016 #define ZBS DDZONECOMM_BUFSIZE
1017 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1018 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1021 for (i = 0; i < n_s; i++)
1023 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1024 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1025 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1026 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1027 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1028 vbuf_s[i*ZBS+1][2] = 0;
1029 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1030 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1031 vbuf_s[i*ZBS+2][2] = 0;
1034 dd_sendrecv_rvec(dd, ddimind, direction,
1038 for (i = 0; i < n_r; i++)
1040 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1041 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1042 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1043 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1044 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1045 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1046 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1052 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1053 rvec cell_ns_x0, rvec cell_ns_x1)
1055 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1057 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1058 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1059 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1060 rvec extr_s[2], extr_r[2];
1062 real dist_d, c = 0, det;
1063 gmx_domdec_comm_t *comm;
1064 gmx_bool bPBC, bUse;
1068 for (d = 1; d < dd->ndim; d++)
1071 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1072 zp->min0 = cell_ns_x0[dim];
1073 zp->max1 = cell_ns_x1[dim];
1074 zp->min1 = cell_ns_x1[dim];
1075 zp->mch0 = cell_ns_x0[dim];
1076 zp->mch1 = cell_ns_x1[dim];
1077 zp->p1_0 = cell_ns_x0[dim];
1078 zp->p1_1 = cell_ns_x1[dim];
1081 for (d = dd->ndim-2; d >= 0; d--)
1084 bPBC = (dim < ddbox->npbcdim);
1086 /* Use an rvec to store two reals */
1087 extr_s[d][0] = comm->cell_f0[d+1];
1088 extr_s[d][1] = comm->cell_f1[d+1];
1089 extr_s[d][2] = comm->cell_f1[d+1];
1092 /* Store the extremes in the backward sending buffer,
1093 * so the get updated separately from the forward communication.
1095 for (d1 = d; d1 < dd->ndim-1; d1++)
1097 /* We invert the order to be able to use the same loop for buf_e */
1098 buf_s[pos].min0 = extr_s[d1][1];
1099 buf_s[pos].max1 = extr_s[d1][0];
1100 buf_s[pos].min1 = extr_s[d1][2];
1101 buf_s[pos].mch0 = 0;
1102 buf_s[pos].mch1 = 0;
1103 /* Store the cell corner of the dimension we communicate along */
1104 buf_s[pos].p1_0 = comm->cell_x0[dim];
1105 buf_s[pos].p1_1 = 0;
1109 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1112 if (dd->ndim == 3 && d == 0)
1114 buf_s[pos] = comm->zone_d2[0][1];
1116 buf_s[pos] = comm->zone_d1[0];
1120 /* We only need to communicate the extremes
1121 * in the forward direction
1123 npulse = comm->cd[d].np;
1126 /* Take the minimum to avoid double communication */
1127 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1131 /* Without PBC we should really not communicate over
1132 * the boundaries, but implementing that complicates
1133 * the communication setup and therefore we simply
1134 * do all communication, but ignore some data.
1136 npulse_min = npulse;
1138 for (p = 0; p < npulse_min; p++)
1140 /* Communicate the extremes forward */
1141 bUse = (bPBC || dd->ci[dim] > 0);
1143 dd_sendrecv_rvec(dd, d, dddirForward,
1144 extr_s+d, dd->ndim-d-1,
1145 extr_r+d, dd->ndim-d-1);
1149 for (d1 = d; d1 < dd->ndim-1; d1++)
1151 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1152 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1153 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1159 for (p = 0; p < npulse; p++)
1161 /* Communicate all the zone information backward */
1162 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1164 dd_sendrecv_ddzone(dd, d, dddirBackward,
1171 for (d1 = d+1; d1 < dd->ndim; d1++)
1173 /* Determine the decrease of maximum required
1174 * communication height along d1 due to the distance along d,
1175 * this avoids a lot of useless atom communication.
1177 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1179 if (ddbox->tric_dir[dim])
1181 /* c is the off-diagonal coupling between the cell planes
1182 * along directions d and d1.
1184 c = ddbox->v[dim][dd->dim[d1]][dim];
1190 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1193 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1197 /* A negative value signals out of range */
1203 /* Accumulate the extremes over all pulses */
1204 for (i = 0; i < buf_size; i++)
1208 buf_e[i] = buf_r[i];
1214 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1215 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1216 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1219 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1227 if (bUse && dh[d1] >= 0)
1229 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1230 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1233 /* Copy the received buffer to the send buffer,
1234 * to pass the data through with the next pulse.
1236 buf_s[i] = buf_r[i];
1238 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1239 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1241 /* Store the extremes */
1244 for (d1 = d; d1 < dd->ndim-1; d1++)
1246 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1247 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1248 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1252 if (d == 1 || (d == 0 && dd->ndim == 3))
1254 for (i = d; i < 2; i++)
1256 comm->zone_d2[1-d][i] = buf_e[pos];
1262 comm->zone_d1[1] = buf_e[pos];
1272 for (i = 0; i < 2; i++)
1276 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1278 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1279 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1285 for (i = 0; i < 2; i++)
1287 for (j = 0; j < 2; j++)
1291 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1293 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1294 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1298 for (d = 1; d < dd->ndim; d++)
1300 comm->cell_f_max0[d] = extr_s[d-1][0];
1301 comm->cell_f_min1[d] = extr_s[d-1][1];
1304 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1305 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1310 static void dd_collect_cg(gmx_domdec_t *dd,
1311 t_state *state_local)
1313 gmx_domdec_master_t *ma = NULL;
1314 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1317 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1319 /* The master has the correct distribution */
1323 if (state_local->ddp_count == dd->ddp_count)
1325 ncg_home = dd->ncg_home;
1327 nat_home = dd->nat_home;
1329 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1331 cgs_gl = &dd->comm->cgs_gl;
1333 ncg_home = state_local->ncg_gl;
1334 cg = state_local->cg_gl;
1336 for (i = 0; i < ncg_home; i++)
1338 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1343 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1346 buf2[0] = dd->ncg_home;
1347 buf2[1] = dd->nat_home;
1357 /* Collect the charge group and atom counts on the master */
1358 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1363 for (i = 0; i < dd->nnodes; i++)
1365 ma->ncg[i] = ma->ibuf[2*i];
1366 ma->nat[i] = ma->ibuf[2*i+1];
1367 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1370 /* Make byte counts and indices */
1371 for (i = 0; i < dd->nnodes; i++)
1373 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1374 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1378 fprintf(debug, "Initial charge group distribution: ");
1379 for (i = 0; i < dd->nnodes; i++)
1381 fprintf(debug, " %d", ma->ncg[i]);
1383 fprintf(debug, "\n");
1387 /* Collect the charge group indices on the master */
1389 dd->ncg_home*sizeof(int), dd->index_gl,
1390 DDMASTER(dd) ? ma->ibuf : NULL,
1391 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1392 DDMASTER(dd) ? ma->cg : NULL);
1394 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1397 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1400 gmx_domdec_master_t *ma;
1401 int n, i, c, a, nalloc = 0;
1410 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1411 dd->rank, dd->mpi_comm_all);
1416 /* Copy the master coordinates to the global array */
1417 cgs_gl = &dd->comm->cgs_gl;
1419 n = DDMASTERRANK(dd);
1421 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1423 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1425 copy_rvec(lv[a++], v[c]);
1429 for (n = 0; n < dd->nnodes; n++)
1433 if (ma->nat[n] > nalloc)
1435 nalloc = over_alloc_dd(ma->nat[n]);
1436 srenew(buf, nalloc);
1439 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1440 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1443 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1445 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1447 copy_rvec(buf[a++], v[c]);
1456 static void get_commbuffer_counts(gmx_domdec_t *dd,
1457 int **counts, int **disps)
1459 gmx_domdec_master_t *ma;
1464 /* Make the rvec count and displacment arrays */
1466 *disps = ma->ibuf + dd->nnodes;
1467 for (n = 0; n < dd->nnodes; n++)
1469 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1470 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1474 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1477 gmx_domdec_master_t *ma;
1478 int *rcounts = NULL, *disps = NULL;
1487 get_commbuffer_counts(dd, &rcounts, &disps);
1492 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1496 cgs_gl = &dd->comm->cgs_gl;
1499 for (n = 0; n < dd->nnodes; n++)
1501 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1503 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1505 copy_rvec(buf[a++], v[c]);
1512 void dd_collect_vec(gmx_domdec_t *dd,
1513 t_state *state_local, rvec *lv, rvec *v)
1515 gmx_domdec_master_t *ma;
1516 int n, i, c, a, nalloc = 0;
1519 dd_collect_cg(dd, state_local);
1521 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1523 dd_collect_vec_sendrecv(dd, lv, v);
1527 dd_collect_vec_gatherv(dd, lv, v);
1532 void dd_collect_state(gmx_domdec_t *dd,
1533 t_state *state_local, t_state *state)
1537 nh = state->nhchainlength;
1541 for (i = 0; i < efptNR; i++)
1543 state->lambda[i] = state_local->lambda[i];
1545 state->fep_state = state_local->fep_state;
1546 state->veta = state_local->veta;
1547 state->vol0 = state_local->vol0;
1548 copy_mat(state_local->box, state->box);
1549 copy_mat(state_local->boxv, state->boxv);
1550 copy_mat(state_local->svir_prev, state->svir_prev);
1551 copy_mat(state_local->fvir_prev, state->fvir_prev);
1552 copy_mat(state_local->pres_prev, state->pres_prev);
1555 for (i = 0; i < state_local->ngtc; i++)
1557 for (j = 0; j < nh; j++)
1559 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1560 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1562 state->therm_integral[i] = state_local->therm_integral[i];
1564 for (i = 0; i < state_local->nnhpres; i++)
1566 for (j = 0; j < nh; j++)
1568 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1569 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1573 for (est = 0; est < estNR; est++)
1575 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1580 dd_collect_vec(dd, state_local, state_local->x, state->x);
1583 dd_collect_vec(dd, state_local, state_local->v, state->v);
1586 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1589 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1592 if (state->nrngi == 1)
1596 for (i = 0; i < state_local->nrng; i++)
1598 state->ld_rng[i] = state_local->ld_rng[i];
1604 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1605 state_local->ld_rng, state->ld_rng);
1609 if (state->nrngi == 1)
1613 state->ld_rngi[0] = state_local->ld_rngi[0];
1618 dd_gather(dd, sizeof(state->ld_rngi[0]),
1619 state_local->ld_rngi, state->ld_rngi);
1622 case estDISRE_INITF:
1623 case estDISRE_RM3TAV:
1624 case estORIRE_INITF:
1628 gmx_incons("Unknown state entry encountered in dd_collect_state");
1634 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1640 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1643 state->nalloc = over_alloc_dd(nalloc);
1645 for (est = 0; est < estNR; est++)
1647 if (EST_DISTR(est) && (state->flags & (1<<est)))
1652 srenew(state->x, state->nalloc);
1655 srenew(state->v, state->nalloc);
1658 srenew(state->sd_X, state->nalloc);
1661 srenew(state->cg_p, state->nalloc);
1665 case estDISRE_INITF:
1666 case estDISRE_RM3TAV:
1667 case estORIRE_INITF:
1669 /* No reallocation required */
1672 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1679 srenew(*f, state->nalloc);
1683 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1686 if (nalloc > fr->cg_nalloc)
1690 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1692 fr->cg_nalloc = over_alloc_dd(nalloc);
1693 srenew(fr->cginfo, fr->cg_nalloc);
1694 if (fr->cutoff_scheme == ecutsGROUP)
1696 srenew(fr->cg_cm, fr->cg_nalloc);
1699 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1701 /* We don't use charge groups, we use x in state to set up
1702 * the atom communication.
1704 dd_realloc_state(state, f, nalloc);
1708 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1711 gmx_domdec_master_t *ma;
1712 int n, i, c, a, nalloc = 0;
1719 for (n = 0; n < dd->nnodes; n++)
1723 if (ma->nat[n] > nalloc)
1725 nalloc = over_alloc_dd(ma->nat[n]);
1726 srenew(buf, nalloc);
1728 /* Use lv as a temporary buffer */
1730 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1732 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1734 copy_rvec(v[c], buf[a++]);
1737 if (a != ma->nat[n])
1739 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1744 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1745 DDRANK(dd, n), n, dd->mpi_comm_all);
1750 n = DDMASTERRANK(dd);
1752 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1754 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1756 copy_rvec(v[c], lv[a++]);
1763 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1764 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1769 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1772 gmx_domdec_master_t *ma;
1773 int *scounts = NULL, *disps = NULL;
1774 int n, i, c, a, nalloc = 0;
1781 get_commbuffer_counts(dd, &scounts, &disps);
1785 for (n = 0; n < dd->nnodes; n++)
1787 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1789 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1791 copy_rvec(v[c], buf[a++]);
1797 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1800 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1802 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1804 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1808 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1812 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1813 t_state *state, t_state *state_local,
1818 nh = state->nhchainlength;
1822 for (i = 0; i < efptNR; i++)
1824 state_local->lambda[i] = state->lambda[i];
1826 state_local->fep_state = state->fep_state;
1827 state_local->veta = state->veta;
1828 state_local->vol0 = state->vol0;
1829 copy_mat(state->box, state_local->box);
1830 copy_mat(state->box_rel, state_local->box_rel);
1831 copy_mat(state->boxv, state_local->boxv);
1832 copy_mat(state->svir_prev, state_local->svir_prev);
1833 copy_mat(state->fvir_prev, state_local->fvir_prev);
1834 for (i = 0; i < state_local->ngtc; i++)
1836 for (j = 0; j < nh; j++)
1838 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1839 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1841 state_local->therm_integral[i] = state->therm_integral[i];
1843 for (i = 0; i < state_local->nnhpres; i++)
1845 for (j = 0; j < nh; j++)
1847 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1848 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1852 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1853 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1854 dd_bcast(dd, sizeof(real), &state_local->veta);
1855 dd_bcast(dd, sizeof(real), &state_local->vol0);
1856 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1857 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1858 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1859 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1860 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1861 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1862 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1863 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1864 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1865 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1867 if (dd->nat_home > state_local->nalloc)
1869 dd_realloc_state(state_local, f, dd->nat_home);
1871 for (i = 0; i < estNR; i++)
1873 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1878 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1881 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1884 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1887 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1890 if (state->nrngi == 1)
1893 state_local->nrng*sizeof(state_local->ld_rng[0]),
1894 state->ld_rng, state_local->ld_rng);
1899 state_local->nrng*sizeof(state_local->ld_rng[0]),
1900 state->ld_rng, state_local->ld_rng);
1904 if (state->nrngi == 1)
1906 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1907 state->ld_rngi, state_local->ld_rngi);
1911 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1912 state->ld_rngi, state_local->ld_rngi);
1915 case estDISRE_INITF:
1916 case estDISRE_RM3TAV:
1917 case estORIRE_INITF:
1919 /* Not implemented yet */
1922 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1928 static char dim2char(int dim)
1934 case XX: c = 'X'; break;
1935 case YY: c = 'Y'; break;
1936 case ZZ: c = 'Z'; break;
1937 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1943 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1944 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1946 rvec grid_s[2], *grid_r = NULL, cx, r;
1947 char fname[STRLEN], format[STRLEN], buf[22];
1949 int a, i, d, z, y, x;
1953 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1954 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1958 snew(grid_r, 2*dd->nnodes);
1961 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1965 for (d = 0; d < DIM; d++)
1967 for (i = 0; i < DIM; i++)
1975 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1977 tric[d][i] = box[i][d]/box[i][i];
1986 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1987 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
1988 out = gmx_fio_fopen(fname, "w");
1989 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1991 for (i = 0; i < dd->nnodes; i++)
1993 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1994 for (d = 0; d < DIM; d++)
1996 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1998 for (z = 0; z < 2; z++)
2000 for (y = 0; y < 2; y++)
2002 for (x = 0; x < 2; x++)
2004 cx[XX] = grid_r[i*2+x][XX];
2005 cx[YY] = grid_r[i*2+y][YY];
2006 cx[ZZ] = grid_r[i*2+z][ZZ];
2008 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
2009 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
2013 for (d = 0; d < DIM; d++)
2015 for (x = 0; x < 4; x++)
2019 case 0: y = 1 + i*8 + 2*x; break;
2020 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2021 case 2: y = 1 + i*8 + x; break;
2023 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2027 gmx_fio_fclose(out);
2032 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2033 gmx_mtop_t *mtop, t_commrec *cr,
2034 int natoms, rvec x[], matrix box)
2036 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2038 int i, ii, resnr, c;
2039 char *atomname, *resname;
2046 natoms = dd->comm->nat[ddnatVSITE];
2049 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2051 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2052 sprintf(format4, "%s%s\n", pdbformat4, "%6.2f%6.2f");
2054 out = gmx_fio_fopen(fname, "w");
2056 fprintf(out, "TITLE %s\n", title);
2057 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2058 for (i = 0; i < natoms; i++)
2060 ii = dd->gatindex[i];
2061 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2062 if (i < dd->comm->nat[ddnatZONE])
2065 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2071 else if (i < dd->comm->nat[ddnatVSITE])
2073 b = dd->comm->zones.n;
2077 b = dd->comm->zones.n + 1;
2079 fprintf(out, strlen(atomname) < 4 ? format : format4,
2080 "ATOM", (ii+1)%100000,
2081 atomname, resname, ' ', resnr%10000, ' ',
2082 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2084 fprintf(out, "TER\n");
2086 gmx_fio_fclose(out);
2089 real dd_cutoff_mbody(gmx_domdec_t *dd)
2091 gmx_domdec_comm_t *comm;
2098 if (comm->bInterCGBondeds)
2100 if (comm->cutoff_mbody > 0)
2102 r = comm->cutoff_mbody;
2106 /* cutoff_mbody=0 means we do not have DLB */
2107 r = comm->cellsize_min[dd->dim[0]];
2108 for (di = 1; di < dd->ndim; di++)
2110 r = min(r, comm->cellsize_min[dd->dim[di]]);
2112 if (comm->bBondComm)
2114 r = max(r, comm->cutoff_mbody);
2118 r = min(r, comm->cutoff);
2126 real dd_cutoff_twobody(gmx_domdec_t *dd)
2130 r_mb = dd_cutoff_mbody(dd);
2132 return max(dd->comm->cutoff, r_mb);
2136 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2140 nc = dd->nc[dd->comm->cartpmedim];
2141 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2142 copy_ivec(coord, coord_pme);
2143 coord_pme[dd->comm->cartpmedim] =
2144 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2147 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2149 /* Here we assign a PME node to communicate with this DD node
2150 * by assuming that the major index of both is x.
2151 * We add cr->npmenodes/2 to obtain an even distribution.
2153 return (ddindex*npme + npme/2)/ndd;
2156 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2158 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2161 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2163 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2166 static int *dd_pmenodes(t_commrec *cr)
2171 snew(pmenodes, cr->npmenodes);
2173 for (i = 0; i < cr->dd->nnodes; i++)
2175 p0 = cr_ddindex2pmeindex(cr, i);
2176 p1 = cr_ddindex2pmeindex(cr, i+1);
2177 if (i+1 == cr->dd->nnodes || p1 > p0)
2181 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2183 pmenodes[n] = i + 1 + n;
2191 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2194 ivec coords, coords_pme, nc;
2199 if (dd->comm->bCartesian) {
2200 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2201 dd_coords2pmecoords(dd,coords,coords_pme);
2202 copy_ivec(dd->ntot,nc);
2203 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2204 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2206 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2208 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2214 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2219 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2221 gmx_domdec_comm_t *comm;
2223 int ddindex, nodeid = -1;
2225 comm = cr->dd->comm;
2230 if (comm->bCartesianPP_PME)
2233 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2238 ddindex = dd_index(cr->dd->nc, coords);
2239 if (comm->bCartesianPP)
2241 nodeid = comm->ddindex2simnodeid[ddindex];
2247 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2259 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2262 gmx_domdec_comm_t *comm;
2263 ivec coord, coord_pme;
2270 /* This assumes a uniform x domain decomposition grid cell size */
2271 if (comm->bCartesianPP_PME)
2274 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2275 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2277 /* This is a PP node */
2278 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2279 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2283 else if (comm->bCartesianPP)
2285 if (sim_nodeid < dd->nnodes)
2287 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2292 /* This assumes DD cells with identical x coordinates
2293 * are numbered sequentially.
2295 if (dd->comm->pmenodes == NULL)
2297 if (sim_nodeid < dd->nnodes)
2299 /* The DD index equals the nodeid */
2300 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2306 while (sim_nodeid > dd->comm->pmenodes[i])
2310 if (sim_nodeid < dd->comm->pmenodes[i])
2312 pmenode = dd->comm->pmenodes[i];
2320 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2322 gmx_bool bPMEOnlyNode;
2324 if (DOMAINDECOMP(cr))
2326 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2330 bPMEOnlyNode = FALSE;
2333 return bPMEOnlyNode;
2336 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2337 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2341 ivec coord, coord_pme;
2345 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2348 for (x = 0; x < dd->nc[XX]; x++)
2350 for (y = 0; y < dd->nc[YY]; y++)
2352 for (z = 0; z < dd->nc[ZZ]; z++)
2354 if (dd->comm->bCartesianPP_PME)
2359 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2360 if (dd->ci[XX] == coord_pme[XX] &&
2361 dd->ci[YY] == coord_pme[YY] &&
2362 dd->ci[ZZ] == coord_pme[ZZ])
2364 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2369 /* The slab corresponds to the nodeid in the PME group */
2370 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2372 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2379 /* The last PP-only node is the peer node */
2380 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2384 fprintf(debug, "Receive coordinates from PP nodes:");
2385 for (x = 0; x < *nmy_ddnodes; x++)
2387 fprintf(debug, " %d", (*my_ddnodes)[x]);
2389 fprintf(debug, "\n");
2393 static gmx_bool receive_vir_ener(t_commrec *cr)
2395 gmx_domdec_comm_t *comm;
2396 int pmenode, coords[DIM], rank;
2400 if (cr->npmenodes < cr->dd->nnodes)
2402 comm = cr->dd->comm;
2403 if (comm->bCartesianPP_PME)
2405 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2407 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2408 coords[comm->cartpmedim]++;
2409 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2411 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2412 if (dd_simnode2pmenode(cr, rank) == pmenode)
2414 /* This is not the last PP node for pmenode */
2422 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2423 if (cr->sim_nodeid+1 < cr->nnodes &&
2424 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2426 /* This is not the last PP node for pmenode */
2435 static void set_zones_ncg_home(gmx_domdec_t *dd)
2437 gmx_domdec_zones_t *zones;
2440 zones = &dd->comm->zones;
2442 zones->cg_range[0] = 0;
2443 for (i = 1; i < zones->n+1; i++)
2445 zones->cg_range[i] = dd->ncg_home;
2449 static void rebuild_cgindex(gmx_domdec_t *dd,
2450 const int *gcgs_index, t_state *state)
2452 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2455 dd_cg_gl = dd->index_gl;
2456 cgindex = dd->cgindex;
2459 for (i = 0; i < state->ncg_gl; i++)
2463 dd_cg_gl[i] = cg_gl;
2464 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2468 dd->ncg_home = state->ncg_gl;
2471 set_zones_ncg_home(dd);
2474 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2476 while (cg >= cginfo_mb->cg_end)
2481 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2484 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2485 t_forcerec *fr, char *bLocalCG)
2487 cginfo_mb_t *cginfo_mb;
2493 cginfo_mb = fr->cginfo_mb;
2494 cginfo = fr->cginfo;
2496 for (cg = cg0; cg < cg1; cg++)
2498 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2502 if (bLocalCG != NULL)
2504 for (cg = cg0; cg < cg1; cg++)
2506 bLocalCG[index_gl[cg]] = TRUE;
2511 static void make_dd_indices(gmx_domdec_t *dd,
2512 const int *gcgs_index, int cg_start)
2514 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2515 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2520 bLocalCG = dd->comm->bLocalCG;
2522 if (dd->nat_tot > dd->gatindex_nalloc)
2524 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2525 srenew(dd->gatindex, dd->gatindex_nalloc);
2528 nzone = dd->comm->zones.n;
2529 zone2cg = dd->comm->zones.cg_range;
2530 zone_ncg1 = dd->comm->zone_ncg1;
2531 index_gl = dd->index_gl;
2532 gatindex = dd->gatindex;
2533 bCGs = dd->comm->bCGs;
2535 if (zone2cg[1] != dd->ncg_home)
2537 gmx_incons("dd->ncg_zone is not up to date");
2540 /* Make the local to global and global to local atom index */
2541 a = dd->cgindex[cg_start];
2542 for (zone = 0; zone < nzone; zone++)
2550 cg0 = zone2cg[zone];
2552 cg1 = zone2cg[zone+1];
2553 cg1_p1 = cg0 + zone_ncg1[zone];
2555 for (cg = cg0; cg < cg1; cg++)
2560 /* Signal that this cg is from more than one pulse away */
2563 cg_gl = index_gl[cg];
2566 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2569 ga2la_set(dd->ga2la, a_gl, a, zone1);
2575 gatindex[a] = cg_gl;
2576 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2583 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2586 int ncg, i, ngl, nerr;
2589 if (bLocalCG == NULL)
2593 for (i = 0; i < dd->ncg_tot; i++)
2595 if (!bLocalCG[dd->index_gl[i]])
2598 "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);
2603 for (i = 0; i < ncg_sys; i++)
2610 if (ngl != dd->ncg_tot)
2612 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);
2619 static void check_index_consistency(gmx_domdec_t *dd,
2620 int natoms_sys, int ncg_sys,
2623 int nerr, ngl, i, a, cell;
2628 if (dd->comm->DD_debug > 1)
2630 snew(have, natoms_sys);
2631 for (a = 0; a < dd->nat_tot; a++)
2633 if (have[dd->gatindex[a]] > 0)
2635 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);
2639 have[dd->gatindex[a]] = a + 1;
2645 snew(have, dd->nat_tot);
2648 for (i = 0; i < natoms_sys; i++)
2650 if (ga2la_get(dd->ga2la, i, &a, &cell))
2652 if (a >= dd->nat_tot)
2654 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);
2660 if (dd->gatindex[a] != i)
2662 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);
2669 if (ngl != dd->nat_tot)
2672 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2673 dd->rank, where, ngl, dd->nat_tot);
2675 for (a = 0; a < dd->nat_tot; a++)
2680 "DD node %d, %s: local atom %d, global %d has no global index\n",
2681 dd->rank, where, a+1, dd->gatindex[a]+1);
2686 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2690 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2691 dd->rank, where, nerr);
2695 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2702 /* Clear the whole list without searching */
2703 ga2la_clear(dd->ga2la);
2707 for (i = a_start; i < dd->nat_tot; i++)
2709 ga2la_del(dd->ga2la, dd->gatindex[i]);
2713 bLocalCG = dd->comm->bLocalCG;
2716 for (i = cg_start; i < dd->ncg_tot; i++)
2718 bLocalCG[dd->index_gl[i]] = FALSE;
2722 dd_clear_local_vsite_indices(dd);
2724 if (dd->constraints)
2726 dd_clear_local_constraint_indices(dd);
2730 /* This function should be used for moving the domain boudaries during DLB,
2731 * for obtaining the minimum cell size. It checks the initially set limit
2732 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2733 * and, possibly, a longer cut-off limit set for PME load balancing.
2735 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2739 cellsize_min = comm->cellsize_min[dim];
2741 if (!comm->bVacDLBNoLimit && comm->bPMELoadBalDLBLimits)
2743 cellsize_min = max(cellsize_min,
2744 comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2747 return cellsize_min;
2750 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2753 real grid_jump_limit;
2755 /* The distance between the boundaries of cells at distance
2756 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2757 * and by the fact that cells should not be shifted by more than
2758 * half their size, such that cg's only shift by one cell
2759 * at redecomposition.
2761 grid_jump_limit = comm->cellsize_limit;
2762 if (!comm->bVacDLBNoLimit)
2764 if (comm->bPMELoadBalDLBLimits)
2766 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2768 grid_jump_limit = max(grid_jump_limit,
2769 cutoff/comm->cd[dim_ind].np);
2772 return grid_jump_limit;
2775 static gmx_bool check_grid_jump(gmx_large_int_t step,
2781 gmx_domdec_comm_t *comm;
2790 for (d = 1; d < dd->ndim; d++)
2793 limit = grid_jump_limit(comm, cutoff, d);
2794 bfac = ddbox->box_size[dim];
2795 if (ddbox->tric_dir[dim])
2797 bfac *= ddbox->skew_fac[dim];
2799 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2800 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2808 /* This error should never be triggered under normal
2809 * circumstances, but you never know ...
2811 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.",
2812 gmx_step_str(step, buf),
2813 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2821 static int dd_load_count(gmx_domdec_comm_t *comm)
2823 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2826 static float dd_force_load(gmx_domdec_comm_t *comm)
2833 if (comm->eFlop > 1)
2835 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2840 load = comm->cycl[ddCyclF];
2841 if (comm->cycl_n[ddCyclF] > 1)
2843 /* Subtract the maximum of the last n cycle counts
2844 * to get rid of possible high counts due to other soures,
2845 * for instance system activity, that would otherwise
2846 * affect the dynamic load balancing.
2848 load -= comm->cycl_max[ddCyclF];
2855 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2857 gmx_domdec_comm_t *comm;
2862 snew(*dim_f, dd->nc[dim]+1);
2864 for (i = 1; i < dd->nc[dim]; i++)
2866 if (comm->slb_frac[dim])
2868 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2872 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2875 (*dim_f)[dd->nc[dim]] = 1;
2878 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2880 int pmeindex, slab, nso, i;
2883 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2889 ddpme->dim = dimind;
2891 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2893 ddpme->nslab = (ddpme->dim == 0 ?
2894 dd->comm->npmenodes_x :
2895 dd->comm->npmenodes_y);
2897 if (ddpme->nslab <= 1)
2902 nso = dd->comm->npmenodes/ddpme->nslab;
2903 /* Determine for each PME slab the PP location range for dimension dim */
2904 snew(ddpme->pp_min, ddpme->nslab);
2905 snew(ddpme->pp_max, ddpme->nslab);
2906 for (slab = 0; slab < ddpme->nslab; slab++)
2908 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2909 ddpme->pp_max[slab] = 0;
2911 for (i = 0; i < dd->nnodes; i++)
2913 ddindex2xyz(dd->nc, i, xyz);
2914 /* For y only use our y/z slab.
2915 * This assumes that the PME x grid size matches the DD grid size.
2917 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2919 pmeindex = ddindex2pmeindex(dd, i);
2922 slab = pmeindex/nso;
2926 slab = pmeindex % ddpme->nslab;
2928 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2929 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2933 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2936 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2938 if (dd->comm->ddpme[0].dim == XX)
2940 return dd->comm->ddpme[0].maxshift;
2948 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2950 if (dd->comm->ddpme[0].dim == YY)
2952 return dd->comm->ddpme[0].maxshift;
2954 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2956 return dd->comm->ddpme[1].maxshift;
2964 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2965 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2967 gmx_domdec_comm_t *comm;
2970 real range, pme_boundary;
2974 nc = dd->nc[ddpme->dim];
2977 if (!ddpme->dim_match)
2979 /* PP decomposition is not along dim: the worst situation */
2982 else if (ns <= 3 || (bUniform && ns == nc))
2984 /* The optimal situation */
2989 /* We need to check for all pme nodes which nodes they
2990 * could possibly need to communicate with.
2992 xmin = ddpme->pp_min;
2993 xmax = ddpme->pp_max;
2994 /* Allow for atoms to be maximally 2/3 times the cut-off
2995 * out of their DD cell. This is a reasonable balance between
2996 * between performance and support for most charge-group/cut-off
2999 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3000 /* Avoid extra communication when we are exactly at a boundary */
3004 for (s = 0; s < ns; s++)
3006 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3007 pme_boundary = (real)s/ns;
3010 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3012 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3016 pme_boundary = (real)(s+1)/ns;
3019 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3021 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3028 ddpme->maxshift = sh;
3032 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3033 ddpme->dim, ddpme->maxshift);
3037 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3041 for (d = 0; d < dd->ndim; d++)
3044 if (dim < ddbox->nboundeddim &&
3045 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3046 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3048 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",
3049 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3050 dd->nc[dim], dd->comm->cellsize_limit);
3055 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3056 gmx_bool bMaster, ivec npulse)
3058 gmx_domdec_comm_t *comm;
3061 real *cell_x, cell_dx, cellsize;
3065 for (d = 0; d < DIM; d++)
3067 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3069 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3072 cell_dx = ddbox->box_size[d]/dd->nc[d];
3075 for (j = 0; j < dd->nc[d]+1; j++)
3077 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3082 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3083 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3085 cellsize = cell_dx*ddbox->skew_fac[d];
3086 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3090 cellsize_min[d] = cellsize;
3094 /* Statically load balanced grid */
3095 /* Also when we are not doing a master distribution we determine
3096 * all cell borders in a loop to obtain identical values
3097 * to the master distribution case and to determine npulse.
3101 cell_x = dd->ma->cell_x[d];
3105 snew(cell_x, dd->nc[d]+1);
3107 cell_x[0] = ddbox->box0[d];
3108 for (j = 0; j < dd->nc[d]; j++)
3110 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3111 cell_x[j+1] = cell_x[j] + cell_dx;
3112 cellsize = cell_dx*ddbox->skew_fac[d];
3113 while (cellsize*npulse[d] < comm->cutoff &&
3114 npulse[d] < dd->nc[d]-1)
3118 cellsize_min[d] = min(cellsize_min[d], cellsize);
3122 comm->cell_x0[d] = cell_x[dd->ci[d]];
3123 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3127 /* The following limitation is to avoid that a cell would receive
3128 * some of its own home charge groups back over the periodic boundary.
3129 * Double charge groups cause trouble with the global indices.
3131 if (d < ddbox->npbcdim &&
3132 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3134 gmx_fatal_collective(FARGS, NULL, dd,
3135 "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",
3136 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3138 dd->nc[d], dd->nc[d],
3139 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3143 if (!comm->bDynLoadBal)
3145 copy_rvec(cellsize_min, comm->cellsize_min);
3148 for (d = 0; d < comm->npmedecompdim; d++)
3150 set_pme_maxshift(dd, &comm->ddpme[d],
3151 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3152 comm->ddpme[d].slb_dim_f);
3157 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3158 int d, int dim, gmx_domdec_root_t *root,
3160 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3162 gmx_domdec_comm_t *comm;
3163 int ncd, i, j, nmin, nmin_old;
3164 gmx_bool bLimLo, bLimHi;
3166 real fac, halfway, cellsize_limit_f_i, region_size;
3167 gmx_bool bPBC, bLastHi = FALSE;
3168 int nrange[] = {range[0], range[1]};
3170 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3176 bPBC = (dim < ddbox->npbcdim);
3178 cell_size = root->buf_ncd;
3182 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3185 /* First we need to check if the scaling does not make cells
3186 * smaller than the smallest allowed size.
3187 * We need to do this iteratively, since if a cell is too small,
3188 * it needs to be enlarged, which makes all the other cells smaller,
3189 * which could in turn make another cell smaller than allowed.
3191 for (i = range[0]; i < range[1]; i++)
3193 root->bCellMin[i] = FALSE;
3199 /* We need the total for normalization */
3201 for (i = range[0]; i < range[1]; i++)
3203 if (root->bCellMin[i] == FALSE)
3205 fac += cell_size[i];
3208 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3209 /* Determine the cell boundaries */
3210 for (i = range[0]; i < range[1]; i++)
3212 if (root->bCellMin[i] == FALSE)
3214 cell_size[i] *= fac;
3215 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3217 cellsize_limit_f_i = 0;
3221 cellsize_limit_f_i = cellsize_limit_f;
3223 if (cell_size[i] < cellsize_limit_f_i)
3225 root->bCellMin[i] = TRUE;
3226 cell_size[i] = cellsize_limit_f_i;
3230 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3233 while (nmin > nmin_old);
3236 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3237 /* For this check we should not use DD_CELL_MARGIN,
3238 * but a slightly smaller factor,
3239 * since rounding could get use below the limit.
3241 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3244 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",
3245 gmx_step_str(step, buf),
3246 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3247 ncd, comm->cellsize_min[dim]);
3250 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3254 /* Check if the boundary did not displace more than halfway
3255 * each of the cells it bounds, as this could cause problems,
3256 * especially when the differences between cell sizes are large.
3257 * If changes are applied, they will not make cells smaller
3258 * than the cut-off, as we check all the boundaries which
3259 * might be affected by a change and if the old state was ok,
3260 * the cells will at most be shrunk back to their old size.
3262 for (i = range[0]+1; i < range[1]; i++)
3264 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3265 if (root->cell_f[i] < halfway)
3267 root->cell_f[i] = halfway;
3268 /* Check if the change also causes shifts of the next boundaries */
3269 for (j = i+1; j < range[1]; j++)
3271 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3273 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3277 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3278 if (root->cell_f[i] > halfway)
3280 root->cell_f[i] = halfway;
3281 /* Check if the change also causes shifts of the next boundaries */
3282 for (j = i-1; j >= range[0]+1; j--)
3284 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3286 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3293 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3294 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3295 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3296 * for a and b nrange is used */
3299 /* Take care of the staggering of the cell boundaries */
3302 for (i = range[0]; i < range[1]; i++)
3304 root->cell_f_max0[i] = root->cell_f[i];
3305 root->cell_f_min1[i] = root->cell_f[i+1];
3310 for (i = range[0]+1; i < range[1]; i++)
3312 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3313 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3314 if (bLimLo && bLimHi)
3316 /* Both limits violated, try the best we can */
3317 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3318 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3319 nrange[0] = range[0];
3321 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3324 nrange[1] = range[1];
3325 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3331 /* root->cell_f[i] = root->bound_min[i]; */
3332 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3335 else if (bLimHi && !bLastHi)
3338 if (nrange[1] < range[1]) /* found a LimLo before */
3340 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3341 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3342 nrange[0] = nrange[1];
3344 root->cell_f[i] = root->bound_max[i];
3346 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3348 nrange[1] = range[1];
3351 if (nrange[1] < range[1]) /* found last a LimLo */
3353 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3354 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3355 nrange[0] = nrange[1];
3356 nrange[1] = range[1];
3357 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3359 else if (nrange[0] > range[0]) /* found at least one LimHi */
3361 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3368 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3369 int d, int dim, gmx_domdec_root_t *root,
3370 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3371 gmx_bool bUniform, gmx_large_int_t step)
3373 gmx_domdec_comm_t *comm;
3374 int ncd, d1, i, j, pos;
3376 real load_aver, load_i, imbalance, change, change_max, sc;
3377 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3381 int range[] = { 0, 0 };
3385 /* Convert the maximum change from the input percentage to a fraction */
3386 change_limit = comm->dlb_scale_lim*0.01;
3390 bPBC = (dim < ddbox->npbcdim);
3392 cell_size = root->buf_ncd;
3394 /* Store the original boundaries */
3395 for (i = 0; i < ncd+1; i++)
3397 root->old_cell_f[i] = root->cell_f[i];
3401 for (i = 0; i < ncd; i++)
3403 cell_size[i] = 1.0/ncd;
3406 else if (dd_load_count(comm))
3408 load_aver = comm->load[d].sum_m/ncd;
3410 for (i = 0; i < ncd; i++)
3412 /* Determine the relative imbalance of cell i */
3413 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3414 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3415 /* Determine the change of the cell size using underrelaxation */
3416 change = -relax*imbalance;
3417 change_max = max(change_max, max(change, -change));
3419 /* Limit the amount of scaling.
3420 * We need to use the same rescaling for all cells in one row,
3421 * otherwise the load balancing might not converge.
3424 if (change_max > change_limit)
3426 sc *= change_limit/change_max;
3428 for (i = 0; i < ncd; i++)
3430 /* Determine the relative imbalance of cell i */
3431 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3432 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3433 /* Determine the change of the cell size using underrelaxation */
3434 change = -sc*imbalance;
3435 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3439 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3440 cellsize_limit_f *= DD_CELL_MARGIN;
3441 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3442 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3443 if (ddbox->tric_dir[dim])
3445 cellsize_limit_f /= ddbox->skew_fac[dim];
3446 dist_min_f /= ddbox->skew_fac[dim];
3448 if (bDynamicBox && d > 0)
3450 dist_min_f *= DD_PRES_SCALE_MARGIN;
3452 if (d > 0 && !bUniform)
3454 /* Make sure that the grid is not shifted too much */
3455 for (i = 1; i < ncd; i++)
3457 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3459 gmx_incons("Inconsistent DD boundary staggering limits!");
3461 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3462 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3465 root->bound_min[i] += 0.5*space;
3467 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3468 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3471 root->bound_max[i] += 0.5*space;
3476 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3478 root->cell_f_max0[i-1] + dist_min_f,
3479 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3480 root->cell_f_min1[i] - dist_min_f);
3485 root->cell_f[0] = 0;
3486 root->cell_f[ncd] = 1;
3487 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3490 /* After the checks above, the cells should obey the cut-off
3491 * restrictions, but it does not hurt to check.
3493 for (i = 0; i < ncd; i++)
3497 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3498 dim, i, root->cell_f[i], root->cell_f[i+1]);
3501 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3502 root->cell_f[i+1] - root->cell_f[i] <
3503 cellsize_limit_f/DD_CELL_MARGIN)
3507 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3508 gmx_step_str(step, buf), dim2char(dim), i,
3509 (root->cell_f[i+1] - root->cell_f[i])
3510 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3515 /* Store the cell boundaries of the lower dimensions at the end */
3516 for (d1 = 0; d1 < d; d1++)
3518 root->cell_f[pos++] = comm->cell_f0[d1];
3519 root->cell_f[pos++] = comm->cell_f1[d1];
3522 if (d < comm->npmedecompdim)
3524 /* The master determines the maximum shift for
3525 * the coordinate communication between separate PME nodes.
3527 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3529 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3532 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3536 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3537 gmx_ddbox_t *ddbox, int dimind)
3539 gmx_domdec_comm_t *comm;
3544 /* Set the cell dimensions */
3545 dim = dd->dim[dimind];
3546 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3547 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3548 if (dim >= ddbox->nboundeddim)
3550 comm->cell_x0[dim] += ddbox->box0[dim];
3551 comm->cell_x1[dim] += ddbox->box0[dim];
3555 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3556 int d, int dim, real *cell_f_row,
3559 gmx_domdec_comm_t *comm;
3565 /* Each node would only need to know two fractions,
3566 * but it is probably cheaper to broadcast the whole array.
3568 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3569 0, comm->mpi_comm_load[d]);
3571 /* Copy the fractions for this dimension from the buffer */
3572 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3573 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3574 /* The whole array was communicated, so set the buffer position */
3575 pos = dd->nc[dim] + 1;
3576 for (d1 = 0; d1 <= d; d1++)
3580 /* Copy the cell fractions of the lower dimensions */
3581 comm->cell_f0[d1] = cell_f_row[pos++];
3582 comm->cell_f1[d1] = cell_f_row[pos++];
3584 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3586 /* Convert the communicated shift from float to int */
3587 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3590 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3594 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3595 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3596 gmx_bool bUniform, gmx_large_int_t step)
3598 gmx_domdec_comm_t *comm;
3600 gmx_bool bRowMember, bRowRoot;
3605 for (d = 0; d < dd->ndim; d++)
3610 for (d1 = d; d1 < dd->ndim; d1++)
3612 if (dd->ci[dd->dim[d1]] > 0)
3625 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3626 ddbox, bDynamicBox, bUniform, step);
3627 cell_f_row = comm->root[d]->cell_f;
3631 cell_f_row = comm->cell_f_row;
3633 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3638 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3642 /* This function assumes the box is static and should therefore
3643 * not be called when the box has changed since the last
3644 * call to dd_partition_system.
3646 for (d = 0; d < dd->ndim; d++)
3648 relative_to_absolute_cell_bounds(dd, ddbox, d);
3654 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3655 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3656 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3657 gmx_wallcycle_t wcycle)
3659 gmx_domdec_comm_t *comm;
3666 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3667 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3668 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3670 else if (bDynamicBox)
3672 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3675 /* Set the dimensions for which no DD is used */
3676 for (dim = 0; dim < DIM; dim++)
3678 if (dd->nc[dim] == 1)
3680 comm->cell_x0[dim] = 0;
3681 comm->cell_x1[dim] = ddbox->box_size[dim];
3682 if (dim >= ddbox->nboundeddim)
3684 comm->cell_x0[dim] += ddbox->box0[dim];
3685 comm->cell_x1[dim] += ddbox->box0[dim];
3691 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3694 gmx_domdec_comm_dim_t *cd;
3696 for (d = 0; d < dd->ndim; d++)
3698 cd = &dd->comm->cd[d];
3699 np = npulse[dd->dim[d]];
3700 if (np > cd->np_nalloc)
3704 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3705 dim2char(dd->dim[d]), np);
3707 if (DDMASTER(dd) && cd->np_nalloc > 0)
3709 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3711 srenew(cd->ind, np);
3712 for (i = cd->np_nalloc; i < np; i++)
3714 cd->ind[i].index = NULL;
3715 cd->ind[i].nalloc = 0;
3724 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3725 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3726 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3727 gmx_wallcycle_t wcycle)
3729 gmx_domdec_comm_t *comm;
3735 /* Copy the old cell boundaries for the cg displacement check */
3736 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3737 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3739 if (comm->bDynLoadBal)
3743 check_box_size(dd, ddbox);
3745 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3749 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3750 realloc_comm_ind(dd, npulse);
3755 for (d = 0; d < DIM; d++)
3757 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3758 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3763 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3765 rvec cell_ns_x0, rvec cell_ns_x1,
3766 gmx_large_int_t step)
3768 gmx_domdec_comm_t *comm;
3773 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3775 dim = dd->dim[dim_ind];
3777 /* Without PBC we don't have restrictions on the outer cells */
3778 if (!(dim >= ddbox->npbcdim &&
3779 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3780 comm->bDynLoadBal &&
3781 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3782 comm->cellsize_min[dim])
3785 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",
3786 gmx_step_str(step, buf), dim2char(dim),
3787 comm->cell_x1[dim] - comm->cell_x0[dim],
3788 ddbox->skew_fac[dim],
3789 dd->comm->cellsize_min[dim],
3790 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3794 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3796 /* Communicate the boundaries and update cell_ns_x0/1 */
3797 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3798 if (dd->bGridJump && dd->ndim > 1)
3800 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3805 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3809 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3817 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3818 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3827 static void check_screw_box(matrix box)
3829 /* Mathematical limitation */
3830 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3832 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3835 /* Limitation due to the asymmetry of the eighth shell method */
3836 if (box[ZZ][YY] != 0)
3838 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3842 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3843 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3846 gmx_domdec_master_t *ma;
3847 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3848 int i, icg, j, k, k0, k1, d, npbcdim;
3850 rvec box_size, cg_cm;
3852 real nrcg, inv_ncg, pos_d;
3854 gmx_bool bUnbounded, bScrew;
3858 if (tmp_ind == NULL)
3860 snew(tmp_nalloc, dd->nnodes);
3861 snew(tmp_ind, dd->nnodes);
3862 for (i = 0; i < dd->nnodes; i++)
3864 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3865 snew(tmp_ind[i], tmp_nalloc[i]);
3869 /* Clear the count */
3870 for (i = 0; i < dd->nnodes; i++)
3876 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3878 cgindex = cgs->index;
3880 /* Compute the center of geometry for all charge groups */
3881 for (icg = 0; icg < cgs->nr; icg++)
3884 k1 = cgindex[icg+1];
3888 copy_rvec(pos[k0], cg_cm);
3895 for (k = k0; (k < k1); k++)
3897 rvec_inc(cg_cm, pos[k]);
3899 for (d = 0; (d < DIM); d++)
3901 cg_cm[d] *= inv_ncg;
3904 /* Put the charge group in the box and determine the cell index */
3905 for (d = DIM-1; d >= 0; d--)
3908 if (d < dd->npbcdim)
3910 bScrew = (dd->bScrewPBC && d == XX);
3911 if (tric_dir[d] && dd->nc[d] > 1)
3913 /* Use triclinic coordintates for this dimension */
3914 for (j = d+1; j < DIM; j++)
3916 pos_d += cg_cm[j]*tcm[j][d];
3919 while (pos_d >= box[d][d])
3922 rvec_dec(cg_cm, box[d]);
3925 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3926 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3928 for (k = k0; (k < k1); k++)
3930 rvec_dec(pos[k], box[d]);
3933 pos[k][YY] = box[YY][YY] - pos[k][YY];
3934 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3941 rvec_inc(cg_cm, box[d]);
3944 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3945 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3947 for (k = k0; (k < k1); k++)
3949 rvec_inc(pos[k], box[d]);
3952 pos[k][YY] = box[YY][YY] - pos[k][YY];
3953 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3958 /* This could be done more efficiently */
3960 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3965 i = dd_index(dd->nc, ind);
3966 if (ma->ncg[i] == tmp_nalloc[i])
3968 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3969 srenew(tmp_ind[i], tmp_nalloc[i]);
3971 tmp_ind[i][ma->ncg[i]] = icg;
3973 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3977 for (i = 0; i < dd->nnodes; i++)
3980 for (k = 0; k < ma->ncg[i]; k++)
3982 ma->cg[k1++] = tmp_ind[i][k];
3985 ma->index[dd->nnodes] = k1;
3987 for (i = 0; i < dd->nnodes; i++)
3997 fprintf(fplog, "Charge group distribution at step %s:",
3998 gmx_step_str(step, buf));
3999 for (i = 0; i < dd->nnodes; i++)
4001 fprintf(fplog, " %d", ma->ncg[i]);
4003 fprintf(fplog, "\n");
4007 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4008 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4011 gmx_domdec_master_t *ma = NULL;
4014 int *ibuf, buf2[2] = { 0, 0 };
4015 gmx_bool bMaster = DDMASTER(dd);
4022 check_screw_box(box);
4025 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4027 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4028 for (i = 0; i < dd->nnodes; i++)
4030 ma->ibuf[2*i] = ma->ncg[i];
4031 ma->ibuf[2*i+1] = ma->nat[i];
4039 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4041 dd->ncg_home = buf2[0];
4042 dd->nat_home = buf2[1];
4043 dd->ncg_tot = dd->ncg_home;
4044 dd->nat_tot = dd->nat_home;
4045 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4047 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4048 srenew(dd->index_gl, dd->cg_nalloc);
4049 srenew(dd->cgindex, dd->cg_nalloc+1);
4053 for (i = 0; i < dd->nnodes; i++)
4055 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4056 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4061 DDMASTER(dd) ? ma->ibuf : NULL,
4062 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4063 DDMASTER(dd) ? ma->cg : NULL,
4064 dd->ncg_home*sizeof(int), dd->index_gl);
4066 /* Determine the home charge group sizes */
4068 for (i = 0; i < dd->ncg_home; i++)
4070 cg_gl = dd->index_gl[i];
4072 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4077 fprintf(debug, "Home charge groups:\n");
4078 for (i = 0; i < dd->ncg_home; i++)
4080 fprintf(debug, " %d", dd->index_gl[i]);
4083 fprintf(debug, "\n");
4086 fprintf(debug, "\n");
4090 static int compact_and_copy_vec_at(int ncg, int *move,
4093 rvec *src, gmx_domdec_comm_t *comm,
4096 int m, icg, i, i0, i1, nrcg;
4102 for (m = 0; m < DIM*2; m++)
4108 for (icg = 0; icg < ncg; icg++)
4110 i1 = cgindex[icg+1];
4116 /* Compact the home array in place */
4117 for (i = i0; i < i1; i++)
4119 copy_rvec(src[i], src[home_pos++]);
4125 /* Copy to the communication buffer */
4127 pos_vec[m] += 1 + vec*nrcg;
4128 for (i = i0; i < i1; i++)
4130 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4132 pos_vec[m] += (nvec - vec - 1)*nrcg;
4136 home_pos += i1 - i0;
4144 static int compact_and_copy_vec_cg(int ncg, int *move,
4146 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4149 int m, icg, i0, i1, nrcg;
4155 for (m = 0; m < DIM*2; m++)
4161 for (icg = 0; icg < ncg; icg++)
4163 i1 = cgindex[icg+1];
4169 /* Compact the home array in place */
4170 copy_rvec(src[icg], src[home_pos++]);
4176 /* Copy to the communication buffer */
4177 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4178 pos_vec[m] += 1 + nrcg*nvec;
4190 static int compact_ind(int ncg, int *move,
4191 int *index_gl, int *cgindex,
4193 gmx_ga2la_t ga2la, char *bLocalCG,
4196 int cg, nat, a0, a1, a, a_gl;
4201 for (cg = 0; cg < ncg; cg++)
4207 /* Compact the home arrays in place.
4208 * Anything that can be done here avoids access to global arrays.
4210 cgindex[home_pos] = nat;
4211 for (a = a0; a < a1; a++)
4214 gatindex[nat] = a_gl;
4215 /* The cell number stays 0, so we don't need to set it */
4216 ga2la_change_la(ga2la, a_gl, nat);
4219 index_gl[home_pos] = index_gl[cg];
4220 cginfo[home_pos] = cginfo[cg];
4221 /* The charge group remains local, so bLocalCG does not change */
4226 /* Clear the global indices */
4227 for (a = a0; a < a1; a++)
4229 ga2la_del(ga2la, gatindex[a]);
4233 bLocalCG[index_gl[cg]] = FALSE;
4237 cgindex[home_pos] = nat;
4242 static void clear_and_mark_ind(int ncg, int *move,
4243 int *index_gl, int *cgindex, int *gatindex,
4244 gmx_ga2la_t ga2la, char *bLocalCG,
4249 for (cg = 0; cg < ncg; cg++)
4255 /* Clear the global indices */
4256 for (a = a0; a < a1; a++)
4258 ga2la_del(ga2la, gatindex[a]);
4262 bLocalCG[index_gl[cg]] = FALSE;
4264 /* Signal that this cg has moved using the ns cell index.
4265 * Here we set it to -1. fill_grid will change it
4266 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4268 cell_index[cg] = -1;
4273 static void print_cg_move(FILE *fplog,
4275 gmx_large_int_t step, int cg, int dim, int dir,
4276 gmx_bool bHaveLimitdAndCMOld, real limitd,
4277 rvec cm_old, rvec cm_new, real pos_d)
4279 gmx_domdec_comm_t *comm;
4284 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4285 if (bHaveLimitdAndCMOld)
4287 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4288 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4292 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4293 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4295 fprintf(fplog, "distance out of cell %f\n",
4296 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4297 if (bHaveLimitdAndCMOld)
4299 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4300 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4302 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4303 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4304 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4306 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4307 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4309 comm->cell_x0[dim], comm->cell_x1[dim]);
4312 static void cg_move_error(FILE *fplog,
4314 gmx_large_int_t step, int cg, int dim, int dir,
4315 gmx_bool bHaveLimitdAndCMOld, real limitd,
4316 rvec cm_old, rvec cm_new, real pos_d)
4320 print_cg_move(fplog, dd, step, cg, dim, dir,
4321 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4323 print_cg_move(stderr, dd, step, cg, dim, dir,
4324 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4326 "A charge group moved too far between two domain decomposition steps\n"
4327 "This usually means that your system is not well equilibrated");
4330 static void rotate_state_atom(t_state *state, int a)
4334 for (est = 0; est < estNR; est++)
4336 if (EST_DISTR(est) && (state->flags & (1<<est)))
4341 /* Rotate the complete state; for a rectangular box only */
4342 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4343 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4346 state->v[a][YY] = -state->v[a][YY];
4347 state->v[a][ZZ] = -state->v[a][ZZ];
4350 state->sd_X[a][YY] = -state->sd_X[a][YY];
4351 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4354 state->cg_p[a][YY] = -state->cg_p[a][YY];
4355 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4357 case estDISRE_INITF:
4358 case estDISRE_RM3TAV:
4359 case estORIRE_INITF:
4361 /* These are distances, so not affected by rotation */
4364 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4370 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4372 if (natoms > comm->moved_nalloc)
4374 /* Contents should be preserved here */
4375 comm->moved_nalloc = over_alloc_dd(natoms);
4376 srenew(comm->moved, comm->moved_nalloc);
4382 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4385 ivec tric_dir, matrix tcm,
4386 rvec cell_x0, rvec cell_x1,
4387 rvec limitd, rvec limit0, rvec limit1,
4389 int cg_start, int cg_end,
4394 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4395 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4399 real inv_ncg, pos_d;
4402 npbcdim = dd->npbcdim;
4404 for (cg = cg_start; cg < cg_end; cg++)
4411 copy_rvec(state->x[k0], cm_new);
4418 for (k = k0; (k < k1); k++)
4420 rvec_inc(cm_new, state->x[k]);
4422 for (d = 0; (d < DIM); d++)
4424 cm_new[d] = inv_ncg*cm_new[d];
4429 /* Do pbc and check DD cell boundary crossings */
4430 for (d = DIM-1; d >= 0; d--)
4434 bScrew = (dd->bScrewPBC && d == XX);
4435 /* Determine the location of this cg in lattice coordinates */
4439 for (d2 = d+1; d2 < DIM; d2++)
4441 pos_d += cm_new[d2]*tcm[d2][d];
4444 /* Put the charge group in the triclinic unit-cell */
4445 if (pos_d >= cell_x1[d])
4447 if (pos_d >= limit1[d])
4449 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4450 cg_cm[cg], cm_new, pos_d);
4453 if (dd->ci[d] == dd->nc[d] - 1)
4455 rvec_dec(cm_new, state->box[d]);
4458 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4459 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4461 for (k = k0; (k < k1); k++)
4463 rvec_dec(state->x[k], state->box[d]);
4466 rotate_state_atom(state, k);
4471 else if (pos_d < cell_x0[d])
4473 if (pos_d < limit0[d])
4475 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4476 cg_cm[cg], cm_new, pos_d);
4481 rvec_inc(cm_new, state->box[d]);
4484 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4485 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4487 for (k = k0; (k < k1); k++)
4489 rvec_inc(state->x[k], state->box[d]);
4492 rotate_state_atom(state, k);
4498 else if (d < npbcdim)
4500 /* Put the charge group in the rectangular unit-cell */
4501 while (cm_new[d] >= state->box[d][d])
4503 rvec_dec(cm_new, state->box[d]);
4504 for (k = k0; (k < k1); k++)
4506 rvec_dec(state->x[k], state->box[d]);
4509 while (cm_new[d] < 0)
4511 rvec_inc(cm_new, state->box[d]);
4512 for (k = k0; (k < k1); k++)
4514 rvec_inc(state->x[k], state->box[d]);
4520 copy_rvec(cm_new, cg_cm[cg]);
4522 /* Determine where this cg should go */
4525 for (d = 0; d < dd->ndim; d++)
4530 flag |= DD_FLAG_FW(d);
4536 else if (dev[dim] == -1)
4538 flag |= DD_FLAG_BW(d);
4541 if (dd->nc[dim] > 2)
4552 /* Temporarily store the flag in move */
4553 move[cg] = mc + flag;
4557 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4558 gmx_domdec_t *dd, ivec tric_dir,
4559 t_state *state, rvec **f,
4560 t_forcerec *fr, t_mdatoms *md,
4568 int ncg[DIM*2], nat[DIM*2];
4569 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4570 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4571 int sbuf[2], rbuf[2];
4572 int home_pos_cg, home_pos_at, buf_pos;
4574 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4577 real inv_ncg, pos_d;
4579 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4581 cginfo_mb_t *cginfo_mb;
4582 gmx_domdec_comm_t *comm;
4584 int nthread, thread;
4588 check_screw_box(state->box);
4592 if (fr->cutoff_scheme == ecutsGROUP)
4597 for (i = 0; i < estNR; i++)
4603 case estX: /* Always present */ break;
4604 case estV: bV = (state->flags & (1<<i)); break;
4605 case estSDX: bSDX = (state->flags & (1<<i)); break;
4606 case estCGP: bCGP = (state->flags & (1<<i)); break;
4609 case estDISRE_INITF:
4610 case estDISRE_RM3TAV:
4611 case estORIRE_INITF:
4613 /* No processing required */
4616 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4621 if (dd->ncg_tot > comm->nalloc_int)
4623 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4624 srenew(comm->buf_int, comm->nalloc_int);
4626 move = comm->buf_int;
4628 /* Clear the count */
4629 for (c = 0; c < dd->ndim*2; c++)
4635 npbcdim = dd->npbcdim;
4637 for (d = 0; (d < DIM); d++)
4639 limitd[d] = dd->comm->cellsize_min[d];
4640 if (d >= npbcdim && dd->ci[d] == 0)
4642 cell_x0[d] = -GMX_FLOAT_MAX;
4646 cell_x0[d] = comm->cell_x0[d];
4648 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4650 cell_x1[d] = GMX_FLOAT_MAX;
4654 cell_x1[d] = comm->cell_x1[d];
4658 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4659 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4663 /* We check after communication if a charge group moved
4664 * more than one cell. Set the pre-comm check limit to float_max.
4666 limit0[d] = -GMX_FLOAT_MAX;
4667 limit1[d] = GMX_FLOAT_MAX;
4671 make_tric_corr_matrix(npbcdim, state->box, tcm);
4673 cgindex = dd->cgindex;
4675 nthread = gmx_omp_nthreads_get(emntDomdec);
4677 /* Compute the center of geometry for all home charge groups
4678 * and put them in the box and determine where they should go.
4680 #pragma omp parallel for num_threads(nthread) schedule(static)
4681 for (thread = 0; thread < nthread; thread++)
4683 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4684 cell_x0, cell_x1, limitd, limit0, limit1,
4686 ( thread *dd->ncg_home)/nthread,
4687 ((thread+1)*dd->ncg_home)/nthread,
4688 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4692 for (cg = 0; cg < dd->ncg_home; cg++)
4697 flag = mc & ~DD_FLAG_NRCG;
4698 mc = mc & DD_FLAG_NRCG;
4701 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4703 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4704 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4706 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4707 /* We store the cg size in the lower 16 bits
4708 * and the place where the charge group should go
4709 * in the next 6 bits. This saves some communication volume.
4711 nrcg = cgindex[cg+1] - cgindex[cg];
4712 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4718 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4719 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4722 for (i = 0; i < dd->ndim*2; i++)
4724 *ncg_moved += ncg[i];
4741 /* Make sure the communication buffers are large enough */
4742 for (mc = 0; mc < dd->ndim*2; mc++)
4744 nvr = ncg[mc] + nat[mc]*nvec;
4745 if (nvr > comm->cgcm_state_nalloc[mc])
4747 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4748 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4752 switch (fr->cutoff_scheme)
4755 /* Recalculating cg_cm might be cheaper than communicating,
4756 * but that could give rise to rounding issues.
4759 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4760 nvec, cg_cm, comm, bCompact);
4763 /* Without charge groups we send the moved atom coordinates
4764 * over twice. This is so the code below can be used without
4765 * many conditionals for both for with and without charge groups.
4768 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4769 nvec, state->x, comm, FALSE);
4772 home_pos_cg -= *ncg_moved;
4776 gmx_incons("unimplemented");
4782 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4783 nvec, vec++, state->x, comm, bCompact);
4786 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4787 nvec, vec++, state->v, comm, bCompact);
4791 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4792 nvec, vec++, state->sd_X, comm, bCompact);
4796 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4797 nvec, vec++, state->cg_p, comm, bCompact);
4802 compact_ind(dd->ncg_home, move,
4803 dd->index_gl, dd->cgindex, dd->gatindex,
4804 dd->ga2la, comm->bLocalCG,
4809 if (fr->cutoff_scheme == ecutsVERLET)
4811 moved = get_moved(comm, dd->ncg_home);
4813 for (k = 0; k < dd->ncg_home; k++)
4820 moved = fr->ns.grid->cell_index;
4823 clear_and_mark_ind(dd->ncg_home, move,
4824 dd->index_gl, dd->cgindex, dd->gatindex,
4825 dd->ga2la, comm->bLocalCG,
4829 cginfo_mb = fr->cginfo_mb;
4831 *ncg_stay_home = home_pos_cg;
4832 for (d = 0; d < dd->ndim; d++)
4838 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4841 /* Communicate the cg and atom counts */
4846 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4847 d, dir, sbuf[0], sbuf[1]);
4849 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4851 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4853 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4854 srenew(comm->buf_int, comm->nalloc_int);
4857 /* Communicate the charge group indices, sizes and flags */
4858 dd_sendrecv_int(dd, d, dir,
4859 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4860 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4862 nvs = ncg[cdd] + nat[cdd]*nvec;
4863 i = rbuf[0] + rbuf[1] *nvec;
4864 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4866 /* Communicate cgcm and state */
4867 dd_sendrecv_rvec(dd, d, dir,
4868 comm->cgcm_state[cdd], nvs,
4869 comm->vbuf.v+nvr, i);
4870 ncg_recv += rbuf[0];
4871 nat_recv += rbuf[1];
4875 /* Process the received charge groups */
4877 for (cg = 0; cg < ncg_recv; cg++)
4879 flag = comm->buf_int[cg*DD_CGIBS+1];
4881 if (dim >= npbcdim && dd->nc[dim] > 2)
4883 /* No pbc in this dim and more than one domain boundary.
4884 * We do a separate check if a charge group didn't move too far.
4886 if (((flag & DD_FLAG_FW(d)) &&
4887 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4888 ((flag & DD_FLAG_BW(d)) &&
4889 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4891 cg_move_error(fplog, dd, step, cg, dim,
4892 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4894 comm->vbuf.v[buf_pos],
4895 comm->vbuf.v[buf_pos],
4896 comm->vbuf.v[buf_pos][dim]);
4903 /* Check which direction this cg should go */
4904 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4908 /* The cell boundaries for dimension d2 are not equal
4909 * for each cell row of the lower dimension(s),
4910 * therefore we might need to redetermine where
4911 * this cg should go.
4914 /* If this cg crosses the box boundary in dimension d2
4915 * we can use the communicated flag, so we do not
4916 * have to worry about pbc.
4918 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4919 (flag & DD_FLAG_FW(d2))) ||
4920 (dd->ci[dim2] == 0 &&
4921 (flag & DD_FLAG_BW(d2)))))
4923 /* Clear the two flags for this dimension */
4924 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4925 /* Determine the location of this cg
4926 * in lattice coordinates
4928 pos_d = comm->vbuf.v[buf_pos][dim2];
4931 for (d3 = dim2+1; d3 < DIM; d3++)
4934 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4937 /* Check of we are not at the box edge.
4938 * pbc is only handled in the first step above,
4939 * but this check could move over pbc while
4940 * the first step did not due to different rounding.
4942 if (pos_d >= cell_x1[dim2] &&
4943 dd->ci[dim2] != dd->nc[dim2]-1)
4945 flag |= DD_FLAG_FW(d2);
4947 else if (pos_d < cell_x0[dim2] &&
4950 flag |= DD_FLAG_BW(d2);
4952 comm->buf_int[cg*DD_CGIBS+1] = flag;
4955 /* Set to which neighboring cell this cg should go */
4956 if (flag & DD_FLAG_FW(d2))
4960 else if (flag & DD_FLAG_BW(d2))
4962 if (dd->nc[dd->dim[d2]] > 2)
4974 nrcg = flag & DD_FLAG_NRCG;
4977 if (home_pos_cg+1 > dd->cg_nalloc)
4979 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4980 srenew(dd->index_gl, dd->cg_nalloc);
4981 srenew(dd->cgindex, dd->cg_nalloc+1);
4983 /* Set the global charge group index and size */
4984 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4985 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4986 /* Copy the state from the buffer */
4987 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4988 if (fr->cutoff_scheme == ecutsGROUP)
4991 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
4995 /* Set the cginfo */
4996 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4997 dd->index_gl[home_pos_cg]);
5000 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5003 if (home_pos_at+nrcg > state->nalloc)
5005 dd_realloc_state(state, f, home_pos_at+nrcg);
5007 for (i = 0; i < nrcg; i++)
5009 copy_rvec(comm->vbuf.v[buf_pos++],
5010 state->x[home_pos_at+i]);
5014 for (i = 0; i < nrcg; i++)
5016 copy_rvec(comm->vbuf.v[buf_pos++],
5017 state->v[home_pos_at+i]);
5022 for (i = 0; i < nrcg; i++)
5024 copy_rvec(comm->vbuf.v[buf_pos++],
5025 state->sd_X[home_pos_at+i]);
5030 for (i = 0; i < nrcg; i++)
5032 copy_rvec(comm->vbuf.v[buf_pos++],
5033 state->cg_p[home_pos_at+i]);
5037 home_pos_at += nrcg;
5041 /* Reallocate the buffers if necessary */
5042 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5044 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5045 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5047 nvr = ncg[mc] + nat[mc]*nvec;
5048 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5050 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5051 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5053 /* Copy from the receive to the send buffers */
5054 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5055 comm->buf_int + cg*DD_CGIBS,
5056 DD_CGIBS*sizeof(int));
5057 memcpy(comm->cgcm_state[mc][nvr],
5058 comm->vbuf.v[buf_pos],
5059 (1+nrcg*nvec)*sizeof(rvec));
5060 buf_pos += 1 + nrcg*nvec;
5067 /* With sorting (!bCompact) the indices are now only partially up to date
5068 * and ncg_home and nat_home are not the real count, since there are
5069 * "holes" in the arrays for the charge groups that moved to neighbors.
5071 if (fr->cutoff_scheme == ecutsVERLET)
5073 moved = get_moved(comm, home_pos_cg);
5075 for (i = dd->ncg_home; i < home_pos_cg; i++)
5080 dd->ncg_home = home_pos_cg;
5081 dd->nat_home = home_pos_at;
5086 "Finished repartitioning: cgs moved out %d, new home %d\n",
5087 *ncg_moved, dd->ncg_home-*ncg_moved);
5092 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5094 dd->comm->cycl[ddCycl] += cycles;
5095 dd->comm->cycl_n[ddCycl]++;
5096 if (cycles > dd->comm->cycl_max[ddCycl])
5098 dd->comm->cycl_max[ddCycl] = cycles;
5102 static double force_flop_count(t_nrnb *nrnb)
5109 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5111 /* To get closer to the real timings, we half the count
5112 * for the normal loops and again half it for water loops.
5115 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5117 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5121 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5124 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5127 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5129 sum += nrnb->n[i]*cost_nrnb(i);
5132 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5134 sum += nrnb->n[i]*cost_nrnb(i);
5140 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5142 if (dd->comm->eFlop)
5144 dd->comm->flop -= force_flop_count(nrnb);
5147 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5149 if (dd->comm->eFlop)
5151 dd->comm->flop += force_flop_count(nrnb);
5156 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5160 for (i = 0; i < ddCyclNr; i++)
5162 dd->comm->cycl[i] = 0;
5163 dd->comm->cycl_n[i] = 0;
5164 dd->comm->cycl_max[i] = 0;
5167 dd->comm->flop_n = 0;
5170 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5172 gmx_domdec_comm_t *comm;
5173 gmx_domdec_load_t *load;
5174 gmx_domdec_root_t *root = NULL;
5175 int d, dim, cid, i, pos;
5176 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5181 fprintf(debug, "get_load_distribution start\n");
5184 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5188 bSepPME = (dd->pme_nodeid >= 0);
5190 for (d = dd->ndim-1; d >= 0; d--)
5193 /* Check if we participate in the communication in this dimension */
5194 if (d == dd->ndim-1 ||
5195 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5197 load = &comm->load[d];
5200 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5203 if (d == dd->ndim-1)
5205 sbuf[pos++] = dd_force_load(comm);
5206 sbuf[pos++] = sbuf[0];
5209 sbuf[pos++] = sbuf[0];
5210 sbuf[pos++] = cell_frac;
5213 sbuf[pos++] = comm->cell_f_max0[d];
5214 sbuf[pos++] = comm->cell_f_min1[d];
5219 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5220 sbuf[pos++] = comm->cycl[ddCyclPME];
5225 sbuf[pos++] = comm->load[d+1].sum;
5226 sbuf[pos++] = comm->load[d+1].max;
5229 sbuf[pos++] = comm->load[d+1].sum_m;
5230 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5231 sbuf[pos++] = comm->load[d+1].flags;
5234 sbuf[pos++] = comm->cell_f_max0[d];
5235 sbuf[pos++] = comm->cell_f_min1[d];
5240 sbuf[pos++] = comm->load[d+1].mdf;
5241 sbuf[pos++] = comm->load[d+1].pme;
5245 /* Communicate a row in DD direction d.
5246 * The communicators are setup such that the root always has rank 0.
5249 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5250 load->load, load->nload*sizeof(float), MPI_BYTE,
5251 0, comm->mpi_comm_load[d]);
5253 if (dd->ci[dim] == dd->master_ci[dim])
5255 /* We are the root, process this row */
5256 if (comm->bDynLoadBal)
5258 root = comm->root[d];
5268 for (i = 0; i < dd->nc[dim]; i++)
5270 load->sum += load->load[pos++];
5271 load->max = max(load->max, load->load[pos]);
5277 /* This direction could not be load balanced properly,
5278 * therefore we need to use the maximum iso the average load.
5280 load->sum_m = max(load->sum_m, load->load[pos]);
5284 load->sum_m += load->load[pos];
5287 load->cvol_min = min(load->cvol_min, load->load[pos]);
5291 load->flags = (int)(load->load[pos++] + 0.5);
5295 root->cell_f_max0[i] = load->load[pos++];
5296 root->cell_f_min1[i] = load->load[pos++];
5301 load->mdf = max(load->mdf, load->load[pos]);
5303 load->pme = max(load->pme, load->load[pos]);
5307 if (comm->bDynLoadBal && root->bLimited)
5309 load->sum_m *= dd->nc[dim];
5310 load->flags |= (1<<d);
5318 comm->nload += dd_load_count(comm);
5319 comm->load_step += comm->cycl[ddCyclStep];
5320 comm->load_sum += comm->load[0].sum;
5321 comm->load_max += comm->load[0].max;
5322 if (comm->bDynLoadBal)
5324 for (d = 0; d < dd->ndim; d++)
5326 if (comm->load[0].flags & (1<<d))
5328 comm->load_lim[d]++;
5334 comm->load_mdf += comm->load[0].mdf;
5335 comm->load_pme += comm->load[0].pme;
5339 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5343 fprintf(debug, "get_load_distribution finished\n");
5347 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5349 /* Return the relative performance loss on the total run time
5350 * due to the force calculation load imbalance.
5352 if (dd->comm->nload > 0)
5355 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5356 (dd->comm->load_step*dd->nnodes);
5364 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5367 int npp, npme, nnodes, d, limp;
5368 float imbal, pme_f_ratio, lossf, lossp = 0;
5370 gmx_domdec_comm_t *comm;
5373 if (DDMASTER(dd) && comm->nload > 0)
5376 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5377 nnodes = npp + npme;
5378 imbal = comm->load_max*npp/comm->load_sum - 1;
5379 lossf = dd_force_imb_perf_loss(dd);
5380 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5381 fprintf(fplog, "%s", buf);
5382 fprintf(stderr, "\n");
5383 fprintf(stderr, "%s", buf);
5384 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5385 fprintf(fplog, "%s", buf);
5386 fprintf(stderr, "%s", buf);
5388 if (comm->bDynLoadBal)
5390 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5391 for (d = 0; d < dd->ndim; d++)
5393 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5394 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5400 sprintf(buf+strlen(buf), "\n");
5401 fprintf(fplog, "%s", buf);
5402 fprintf(stderr, "%s", buf);
5406 pme_f_ratio = comm->load_pme/comm->load_mdf;
5407 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5410 lossp *= (float)npme/(float)nnodes;
5414 lossp *= (float)npp/(float)nnodes;
5416 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5417 fprintf(fplog, "%s", buf);
5418 fprintf(stderr, "%s", buf);
5419 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5420 fprintf(fplog, "%s", buf);
5421 fprintf(stderr, "%s", buf);
5423 fprintf(fplog, "\n");
5424 fprintf(stderr, "\n");
5426 if (lossf >= DD_PERF_LOSS)
5429 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5430 " in the domain decomposition.\n", lossf*100);
5431 if (!comm->bDynLoadBal)
5433 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5437 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5439 fprintf(fplog, "%s\n", buf);
5440 fprintf(stderr, "%s\n", buf);
5442 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5445 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5446 " had %s work to do than the PP nodes.\n"
5447 " You might want to %s the number of PME nodes\n"
5448 " or %s the cut-off and the grid spacing.\n",
5450 (lossp < 0) ? "less" : "more",
5451 (lossp < 0) ? "decrease" : "increase",
5452 (lossp < 0) ? "decrease" : "increase");
5453 fprintf(fplog, "%s\n", buf);
5454 fprintf(stderr, "%s\n", buf);
5459 static float dd_vol_min(gmx_domdec_t *dd)
5461 return dd->comm->load[0].cvol_min*dd->nnodes;
5464 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5466 return dd->comm->load[0].flags;
5469 static float dd_f_imbal(gmx_domdec_t *dd)
5471 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5474 float dd_pme_f_ratio(gmx_domdec_t *dd)
5476 if (dd->comm->cycl_n[ddCyclPME] > 0)
5478 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5486 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5491 flags = dd_load_flags(dd);
5495 "DD load balancing is limited by minimum cell size in dimension");
5496 for (d = 0; d < dd->ndim; d++)
5500 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5503 fprintf(fplog, "\n");
5505 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5506 if (dd->comm->bDynLoadBal)
5508 fprintf(fplog, " vol min/aver %5.3f%c",
5509 dd_vol_min(dd), flags ? '!' : ' ');
5511 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5512 if (dd->comm->cycl_n[ddCyclPME])
5514 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5516 fprintf(fplog, "\n\n");
5519 static void dd_print_load_verbose(gmx_domdec_t *dd)
5521 if (dd->comm->bDynLoadBal)
5523 fprintf(stderr, "vol %4.2f%c ",
5524 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5526 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5527 if (dd->comm->cycl_n[ddCyclPME])
5529 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5534 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5539 gmx_domdec_root_t *root;
5540 gmx_bool bPartOfGroup = FALSE;
5542 dim = dd->dim[dim_ind];
5543 copy_ivec(loc, loc_c);
5544 for (i = 0; i < dd->nc[dim]; i++)
5547 rank = dd_index(dd->nc, loc_c);
5548 if (rank == dd->rank)
5550 /* This process is part of the group */
5551 bPartOfGroup = TRUE;
5554 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5558 dd->comm->mpi_comm_load[dim_ind] = c_row;
5559 if (dd->comm->eDLB != edlbNO)
5561 if (dd->ci[dim] == dd->master_ci[dim])
5563 /* This is the root process of this row */
5564 snew(dd->comm->root[dim_ind], 1);
5565 root = dd->comm->root[dim_ind];
5566 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5567 snew(root->old_cell_f, dd->nc[dim]+1);
5568 snew(root->bCellMin, dd->nc[dim]);
5571 snew(root->cell_f_max0, dd->nc[dim]);
5572 snew(root->cell_f_min1, dd->nc[dim]);
5573 snew(root->bound_min, dd->nc[dim]);
5574 snew(root->bound_max, dd->nc[dim]);
5576 snew(root->buf_ncd, dd->nc[dim]);
5580 /* This is not a root process, we only need to receive cell_f */
5581 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5584 if (dd->ci[dim] == dd->master_ci[dim])
5586 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5592 static void make_load_communicators(gmx_domdec_t *dd)
5595 int dim0, dim1, i, j;
5600 fprintf(debug, "Making load communicators\n");
5603 snew(dd->comm->load, dd->ndim);
5604 snew(dd->comm->mpi_comm_load, dd->ndim);
5607 make_load_communicator(dd, 0, loc);
5611 for (i = 0; i < dd->nc[dim0]; i++)
5614 make_load_communicator(dd, 1, loc);
5620 for (i = 0; i < dd->nc[dim0]; i++)
5624 for (j = 0; j < dd->nc[dim1]; j++)
5627 make_load_communicator(dd, 2, loc);
5634 fprintf(debug, "Finished making load communicators\n");
5639 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5642 int d, dim, i, j, m;
5645 ivec dd_zp[DD_MAXIZONE];
5646 gmx_domdec_zones_t *zones;
5647 gmx_domdec_ns_ranges_t *izone;
5649 for (d = 0; d < dd->ndim; d++)
5652 copy_ivec(dd->ci, tmp);
5653 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5654 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5655 copy_ivec(dd->ci, tmp);
5656 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5657 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5660 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5663 dd->neighbor[d][1]);
5669 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5671 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5672 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5679 for (i = 0; i < nzonep; i++)
5681 copy_ivec(dd_zp3[i], dd_zp[i]);
5687 for (i = 0; i < nzonep; i++)
5689 copy_ivec(dd_zp2[i], dd_zp[i]);
5695 for (i = 0; i < nzonep; i++)
5697 copy_ivec(dd_zp1[i], dd_zp[i]);
5701 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5706 zones = &dd->comm->zones;
5708 for (i = 0; i < nzone; i++)
5711 clear_ivec(zones->shift[i]);
5712 for (d = 0; d < dd->ndim; d++)
5714 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5719 for (i = 0; i < nzone; i++)
5721 for (d = 0; d < DIM; d++)
5723 s[d] = dd->ci[d] - zones->shift[i][d];
5728 else if (s[d] >= dd->nc[d])
5734 zones->nizone = nzonep;
5735 for (i = 0; i < zones->nizone; i++)
5737 if (dd_zp[i][0] != i)
5739 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5741 izone = &zones->izone[i];
5742 izone->j0 = dd_zp[i][1];
5743 izone->j1 = dd_zp[i][2];
5744 for (dim = 0; dim < DIM; dim++)
5746 if (dd->nc[dim] == 1)
5748 /* All shifts should be allowed */
5749 izone->shift0[dim] = -1;
5750 izone->shift1[dim] = 1;
5755 izone->shift0[d] = 0;
5756 izone->shift1[d] = 0;
5757 for(j=izone->j0; j<izone->j1; j++) {
5758 if (dd->shift[j][d] > dd->shift[i][d])
5759 izone->shift0[d] = -1;
5760 if (dd->shift[j][d] < dd->shift[i][d])
5761 izone->shift1[d] = 1;
5767 /* Assume the shift are not more than 1 cell */
5768 izone->shift0[dim] = 1;
5769 izone->shift1[dim] = -1;
5770 for (j = izone->j0; j < izone->j1; j++)
5772 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5773 if (shift_diff < izone->shift0[dim])
5775 izone->shift0[dim] = shift_diff;
5777 if (shift_diff > izone->shift1[dim])
5779 izone->shift1[dim] = shift_diff;
5786 if (dd->comm->eDLB != edlbNO)
5788 snew(dd->comm->root, dd->ndim);
5791 if (dd->comm->bRecordLoad)
5793 make_load_communicators(dd);
5797 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5800 gmx_domdec_comm_t *comm;
5811 if (comm->bCartesianPP)
5813 /* Set up cartesian communication for the particle-particle part */
5816 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5817 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5820 for (i = 0; i < DIM; i++)
5824 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5826 /* We overwrite the old communicator with the new cartesian one */
5827 cr->mpi_comm_mygroup = comm_cart;
5830 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5831 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5833 if (comm->bCartesianPP_PME)
5835 /* Since we want to use the original cartesian setup for sim,
5836 * and not the one after split, we need to make an index.
5838 snew(comm->ddindex2ddnodeid, dd->nnodes);
5839 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5840 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5841 /* Get the rank of the DD master,
5842 * above we made sure that the master node is a PP node.
5852 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5854 else if (comm->bCartesianPP)
5856 if (cr->npmenodes == 0)
5858 /* The PP communicator is also
5859 * the communicator for this simulation
5861 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5863 cr->nodeid = dd->rank;
5865 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5867 /* We need to make an index to go from the coordinates
5868 * to the nodeid of this simulation.
5870 snew(comm->ddindex2simnodeid, dd->nnodes);
5871 snew(buf, dd->nnodes);
5872 if (cr->duty & DUTY_PP)
5874 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5876 /* Communicate the ddindex to simulation nodeid index */
5877 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5878 cr->mpi_comm_mysim);
5881 /* Determine the master coordinates and rank.
5882 * The DD master should be the same node as the master of this sim.
5884 for (i = 0; i < dd->nnodes; i++)
5886 if (comm->ddindex2simnodeid[i] == 0)
5888 ddindex2xyz(dd->nc, i, dd->master_ci);
5889 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5894 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5899 /* No Cartesian communicators */
5900 /* We use the rank in dd->comm->all as DD index */
5901 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5902 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5904 clear_ivec(dd->master_ci);
5911 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5912 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5917 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5918 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5922 static void receive_ddindex2simnodeid(t_commrec *cr)
5926 gmx_domdec_comm_t *comm;
5933 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5935 snew(comm->ddindex2simnodeid, dd->nnodes);
5936 snew(buf, dd->nnodes);
5937 if (cr->duty & DUTY_PP)
5939 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5942 /* Communicate the ddindex to simulation nodeid index */
5943 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5944 cr->mpi_comm_mysim);
5951 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5952 int ncg, int natoms)
5954 gmx_domdec_master_t *ma;
5959 snew(ma->ncg, dd->nnodes);
5960 snew(ma->index, dd->nnodes+1);
5962 snew(ma->nat, dd->nnodes);
5963 snew(ma->ibuf, dd->nnodes*2);
5964 snew(ma->cell_x, DIM);
5965 for (i = 0; i < DIM; i++)
5967 snew(ma->cell_x[i], dd->nc[i]+1);
5970 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5976 snew(ma->vbuf, natoms);
5982 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5986 gmx_domdec_comm_t *comm;
5997 if (comm->bCartesianPP)
5999 for (i = 1; i < DIM; i++)
6001 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6003 if (bDiv[YY] || bDiv[ZZ])
6005 comm->bCartesianPP_PME = TRUE;
6006 /* If we have 2D PME decomposition, which is always in x+y,
6007 * we stack the PME only nodes in z.
6008 * Otherwise we choose the direction that provides the thinnest slab
6009 * of PME only nodes as this will have the least effect
6010 * on the PP communication.
6011 * But for the PME communication the opposite might be better.
6013 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6015 dd->nc[YY] > dd->nc[ZZ]))
6017 comm->cartpmedim = ZZ;
6021 comm->cartpmedim = YY;
6023 comm->ntot[comm->cartpmedim]
6024 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6028 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]);
6030 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6035 if (comm->bCartesianPP_PME)
6039 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]);
6042 for (i = 0; i < DIM; i++)
6046 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6049 MPI_Comm_rank(comm_cart, &rank);
6050 if (MASTERNODE(cr) && rank != 0)
6052 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6055 /* With this assigment we loose the link to the original communicator
6056 * which will usually be MPI_COMM_WORLD, unless have multisim.
6058 cr->mpi_comm_mysim = comm_cart;
6059 cr->sim_nodeid = rank;
6061 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6065 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6066 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6069 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6073 if (cr->npmenodes == 0 ||
6074 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6076 cr->duty = DUTY_PME;
6079 /* Split the sim communicator into PP and PME only nodes */
6080 MPI_Comm_split(cr->mpi_comm_mysim,
6082 dd_index(comm->ntot, dd->ci),
6083 &cr->mpi_comm_mygroup);
6087 switch (dd_node_order)
6092 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6095 case ddnoINTERLEAVE:
6096 /* Interleave the PP-only and PME-only nodes,
6097 * as on clusters with dual-core machines this will double
6098 * the communication bandwidth of the PME processes
6099 * and thus speed up the PP <-> PME and inter PME communication.
6103 fprintf(fplog, "Interleaving PP and PME nodes\n");
6105 comm->pmenodes = dd_pmenodes(cr);
6110 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6113 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6115 cr->duty = DUTY_PME;
6122 /* Split the sim communicator into PP and PME only nodes */
6123 MPI_Comm_split(cr->mpi_comm_mysim,
6126 &cr->mpi_comm_mygroup);
6127 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6133 fprintf(fplog, "This is a %s only node\n\n",
6134 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6138 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6141 gmx_domdec_comm_t *comm;
6147 copy_ivec(dd->nc, comm->ntot);
6149 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6150 comm->bCartesianPP_PME = FALSE;
6152 /* Reorder the nodes by default. This might change the MPI ranks.
6153 * Real reordering is only supported on very few architectures,
6154 * Blue Gene is one of them.
6156 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6158 if (cr->npmenodes > 0)
6160 /* Split the communicator into a PP and PME part */
6161 split_communicator(fplog, cr, dd_node_order, CartReorder);
6162 if (comm->bCartesianPP_PME)
6164 /* We (possibly) reordered the nodes in split_communicator,
6165 * so it is no longer required in make_pp_communicator.
6167 CartReorder = FALSE;
6172 /* All nodes do PP and PME */
6174 /* We do not require separate communicators */
6175 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6179 if (cr->duty & DUTY_PP)
6181 /* Copy or make a new PP communicator */
6182 make_pp_communicator(fplog, cr, CartReorder);
6186 receive_ddindex2simnodeid(cr);
6189 if (!(cr->duty & DUTY_PME))
6191 /* Set up the commnuication to our PME node */
6192 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6193 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6196 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6197 dd->pme_nodeid, dd->pme_receive_vir_ener);
6202 dd->pme_nodeid = -1;
6207 dd->ma = init_gmx_domdec_master_t(dd,
6209 comm->cgs_gl.index[comm->cgs_gl.nr]);
6213 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6215 real *slb_frac, tot;
6220 if (nc > 1 && size_string != NULL)
6224 fprintf(fplog, "Using static load balancing for the %s direction\n",
6229 for (i = 0; i < nc; i++)
6232 sscanf(size_string, "%lf%n", &dbl, &n);
6235 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6244 fprintf(fplog, "Relative cell sizes:");
6246 for (i = 0; i < nc; i++)
6251 fprintf(fplog, " %5.3f", slb_frac[i]);
6256 fprintf(fplog, "\n");
6263 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6266 gmx_mtop_ilistloop_t iloop;
6270 iloop = gmx_mtop_ilistloop_init(mtop);
6271 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6273 for (ftype = 0; ftype < F_NRE; ftype++)
6275 if ((interaction_function[ftype].flags & IF_BOND) &&
6278 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6286 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6292 val = getenv(env_var);
6295 if (sscanf(val, "%d", &nst) <= 0)
6301 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6309 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6313 fprintf(stderr, "\n%s\n", warn_string);
6317 fprintf(fplog, "\n%s\n", warn_string);
6321 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6322 t_inputrec *ir, FILE *fplog)
6324 if (ir->ePBC == epbcSCREW &&
6325 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6327 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6330 if (ir->ns_type == ensSIMPLE)
6332 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6335 if (ir->nstlist == 0)
6337 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6340 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6342 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6346 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6351 r = ddbox->box_size[XX];
6352 for (di = 0; di < dd->ndim; di++)
6355 /* Check using the initial average cell size */
6356 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6362 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6363 const char *dlb_opt, gmx_bool bRecordLoad,
6364 unsigned long Flags, t_inputrec *ir)
6372 case 'a': eDLB = edlbAUTO; break;
6373 case 'n': eDLB = edlbNO; break;
6374 case 'y': eDLB = edlbYES; break;
6375 default: gmx_incons("Unknown dlb_opt");
6378 if (Flags & MD_RERUN)
6383 if (!EI_DYNAMICS(ir->eI))
6385 if (eDLB == edlbYES)
6387 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6388 dd_warning(cr, fplog, buf);
6396 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6401 if (Flags & MD_REPRODUCIBLE)
6408 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6412 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6415 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6423 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6428 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6430 /* Decomposition order z,y,x */
6433 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6435 for (dim = DIM-1; dim >= 0; dim--)
6437 if (dd->nc[dim] > 1)
6439 dd->dim[dd->ndim++] = dim;
6445 /* Decomposition order x,y,z */
6446 for (dim = 0; dim < DIM; dim++)
6448 if (dd->nc[dim] > 1)
6450 dd->dim[dd->ndim++] = dim;
6456 static gmx_domdec_comm_t *init_dd_comm()
6458 gmx_domdec_comm_t *comm;
6462 snew(comm->cggl_flag, DIM*2);
6463 snew(comm->cgcm_state, DIM*2);
6464 for (i = 0; i < DIM*2; i++)
6466 comm->cggl_flag_nalloc[i] = 0;
6467 comm->cgcm_state_nalloc[i] = 0;
6470 comm->nalloc_int = 0;
6471 comm->buf_int = NULL;
6473 vec_rvec_init(&comm->vbuf);
6475 comm->n_load_have = 0;
6476 comm->n_load_collect = 0;
6478 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6480 comm->sum_nat[i] = 0;
6484 comm->load_step = 0;
6487 clear_ivec(comm->load_lim);
6494 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6495 unsigned long Flags,
6497 real comm_distance_min, real rconstr,
6498 const char *dlb_opt, real dlb_scale,
6499 const char *sizex, const char *sizey, const char *sizez,
6500 gmx_mtop_t *mtop, t_inputrec *ir,
6501 matrix box, rvec *x,
6503 int *npme_x, int *npme_y)
6506 gmx_domdec_comm_t *comm;
6509 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6516 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6521 dd->comm = init_dd_comm();
6523 snew(comm->cggl_flag, DIM*2);
6524 snew(comm->cgcm_state, DIM*2);
6526 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6527 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6529 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6530 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6531 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6532 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6533 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6534 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6535 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6536 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6538 dd->pme_recv_f_alloc = 0;
6539 dd->pme_recv_f_buf = NULL;
6541 if (dd->bSendRecv2 && fplog)
6543 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");
6549 fprintf(fplog, "Will load balance based on FLOP count\n");
6551 if (comm->eFlop > 1)
6553 srand(1+cr->nodeid);
6555 comm->bRecordLoad = TRUE;
6559 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6563 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6565 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6568 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6570 dd->bGridJump = comm->bDynLoadBal;
6571 comm->bPMELoadBalDLBLimits = FALSE;
6573 if (comm->nstSortCG)
6577 if (comm->nstSortCG == 1)
6579 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6583 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6587 snew(comm->sort, 1);
6593 fprintf(fplog, "Will not sort the charge groups\n");
6597 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6599 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6600 if (comm->bInterCGBondeds)
6602 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6606 comm->bInterCGMultiBody = FALSE;
6609 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6610 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6612 if (ir->rlistlong == 0)
6614 /* Set the cut-off to some very large value,
6615 * so we don't need if statements everywhere in the code.
6616 * We use sqrt, since the cut-off is squared in some places.
6618 comm->cutoff = GMX_CUTOFF_INF;
6622 comm->cutoff = ir->rlistlong;
6624 comm->cutoff_mbody = 0;
6626 comm->cellsize_limit = 0;
6627 comm->bBondComm = FALSE;
6629 if (comm->bInterCGBondeds)
6631 if (comm_distance_min > 0)
6633 comm->cutoff_mbody = comm_distance_min;
6634 if (Flags & MD_DDBONDCOMM)
6636 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6640 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6642 r_bonded_limit = comm->cutoff_mbody;
6644 else if (ir->bPeriodicMols)
6646 /* Can not easily determine the required cut-off */
6647 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");
6648 comm->cutoff_mbody = comm->cutoff/2;
6649 r_bonded_limit = comm->cutoff_mbody;
6655 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6656 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6658 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6659 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6661 /* We use an initial margin of 10% for the minimum cell size,
6662 * except when we are just below the non-bonded cut-off.
6664 if (Flags & MD_DDBONDCOMM)
6666 if (max(r_2b, r_mb) > comm->cutoff)
6668 r_bonded = max(r_2b, r_mb);
6669 r_bonded_limit = 1.1*r_bonded;
6670 comm->bBondComm = TRUE;
6675 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6677 /* We determine cutoff_mbody later */
6681 /* No special bonded communication,
6682 * simply increase the DD cut-off.
6684 r_bonded_limit = 1.1*max(r_2b, r_mb);
6685 comm->cutoff_mbody = r_bonded_limit;
6686 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6689 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6693 "Minimum cell size due to bonded interactions: %.3f nm\n",
6694 comm->cellsize_limit);
6698 if (dd->bInterCGcons && rconstr <= 0)
6700 /* There is a cell size limit due to the constraints (P-LINCS) */
6701 rconstr = constr_r_max(fplog, mtop, ir);
6705 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6707 if (rconstr > comm->cellsize_limit)
6709 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6713 else if (rconstr > 0 && fplog)
6715 /* Here we do not check for dd->bInterCGcons,
6716 * because one can also set a cell size limit for virtual sites only
6717 * and at this point we don't know yet if there are intercg v-sites.
6720 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6723 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6725 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6729 copy_ivec(nc, dd->nc);
6730 set_dd_dim(fplog, dd);
6731 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6733 if (cr->npmenodes == -1)
6737 acs = average_cellsize_min(dd, ddbox);
6738 if (acs < comm->cellsize_limit)
6742 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6744 gmx_fatal_collective(FARGS, cr, NULL,
6745 "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",
6746 acs, comm->cellsize_limit);
6751 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6753 /* We need to choose the optimal DD grid and possibly PME nodes */
6754 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6755 comm->eDLB != edlbNO, dlb_scale,
6756 comm->cellsize_limit, comm->cutoff,
6757 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6759 if (dd->nc[XX] == 0)
6761 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6762 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6763 !bC ? "-rdd" : "-rcon",
6764 comm->eDLB != edlbNO ? " or -dds" : "",
6765 bC ? " or your LINCS settings" : "");
6767 gmx_fatal_collective(FARGS, cr, NULL,
6768 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6770 "Look in the log file for details on the domain decomposition",
6771 cr->nnodes-cr->npmenodes, limit, buf);
6773 set_dd_dim(fplog, dd);
6779 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6780 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6783 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6784 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6786 gmx_fatal_collective(FARGS, cr, NULL,
6787 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6788 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6790 if (cr->npmenodes > dd->nnodes)
6792 gmx_fatal_collective(FARGS, cr, NULL,
6793 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6795 if (cr->npmenodes > 0)
6797 comm->npmenodes = cr->npmenodes;
6801 comm->npmenodes = dd->nnodes;
6804 if (EEL_PME(ir->coulombtype))
6806 /* The following choices should match those
6807 * in comm_cost_est in domdec_setup.c.
6808 * Note that here the checks have to take into account
6809 * that the decomposition might occur in a different order than xyz
6810 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6811 * in which case they will not match those in comm_cost_est,
6812 * but since that is mainly for testing purposes that's fine.
6814 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6815 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6816 getenv("GMX_PMEONEDD") == NULL)
6818 comm->npmedecompdim = 2;
6819 comm->npmenodes_x = dd->nc[XX];
6820 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6824 /* In case nc is 1 in both x and y we could still choose to
6825 * decompose pme in y instead of x, but we use x for simplicity.
6827 comm->npmedecompdim = 1;
6828 if (dd->dim[0] == YY)
6830 comm->npmenodes_x = 1;
6831 comm->npmenodes_y = comm->npmenodes;
6835 comm->npmenodes_x = comm->npmenodes;
6836 comm->npmenodes_y = 1;
6841 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6842 comm->npmenodes_x, comm->npmenodes_y, 1);
6847 comm->npmedecompdim = 0;
6848 comm->npmenodes_x = 0;
6849 comm->npmenodes_y = 0;
6852 /* Technically we don't need both of these,
6853 * but it simplifies code not having to recalculate it.
6855 *npme_x = comm->npmenodes_x;
6856 *npme_y = comm->npmenodes_y;
6858 snew(comm->slb_frac, DIM);
6859 if (comm->eDLB == edlbNO)
6861 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6862 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6863 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6866 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6868 if (comm->bBondComm || comm->eDLB != edlbNO)
6870 /* Set the bonded communication distance to halfway
6871 * the minimum and the maximum,
6872 * since the extra communication cost is nearly zero.
6874 acs = average_cellsize_min(dd, ddbox);
6875 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6876 if (comm->eDLB != edlbNO)
6878 /* Check if this does not limit the scaling */
6879 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6881 if (!comm->bBondComm)
6883 /* Without bBondComm do not go beyond the n.b. cut-off */
6884 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6885 if (comm->cellsize_limit >= comm->cutoff)
6887 /* We don't loose a lot of efficieny
6888 * when increasing it to the n.b. cut-off.
6889 * It can even be slightly faster, because we need
6890 * less checks for the communication setup.
6892 comm->cutoff_mbody = comm->cutoff;
6895 /* Check if we did not end up below our original limit */
6896 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6898 if (comm->cutoff_mbody > comm->cellsize_limit)
6900 comm->cellsize_limit = comm->cutoff_mbody;
6903 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6908 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6909 "cellsize limit %f\n",
6910 comm->bBondComm, comm->cellsize_limit);
6915 check_dd_restrictions(cr, dd, ir, fplog);
6918 comm->partition_step = INT_MIN;
6921 clear_dd_cycle_counts(dd);
6926 static void set_dlb_limits(gmx_domdec_t *dd)
6931 for (d = 0; d < dd->ndim; d++)
6933 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6934 dd->comm->cellsize_min[dd->dim[d]] =
6935 dd->comm->cellsize_min_dlb[dd->dim[d]];
6940 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6943 gmx_domdec_comm_t *comm;
6953 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);
6956 cellsize_min = comm->cellsize_min[dd->dim[0]];
6957 for (d = 1; d < dd->ndim; d++)
6959 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6962 if (cellsize_min < comm->cellsize_limit*1.05)
6964 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");
6966 /* Change DLB from "auto" to "no". */
6967 comm->eDLB = edlbNO;
6972 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6973 comm->bDynLoadBal = TRUE;
6974 dd->bGridJump = TRUE;
6978 /* We can set the required cell size info here,
6979 * so we do not need to communicate this.
6980 * The grid is completely uniform.
6982 for (d = 0; d < dd->ndim; d++)
6986 comm->load[d].sum_m = comm->load[d].sum;
6988 nc = dd->nc[dd->dim[d]];
6989 for (i = 0; i < nc; i++)
6991 comm->root[d]->cell_f[i] = i/(real)nc;
6994 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6995 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6998 comm->root[d]->cell_f[nc] = 1.0;
7003 static char *init_bLocalCG(gmx_mtop_t *mtop)
7008 ncg = ncg_mtop(mtop);
7009 snew(bLocalCG, ncg);
7010 for (cg = 0; cg < ncg; cg++)
7012 bLocalCG[cg] = FALSE;
7018 void dd_init_bondeds(FILE *fplog,
7019 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7020 gmx_vsite_t *vsite, gmx_constr_t constr,
7021 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7023 gmx_domdec_comm_t *comm;
7027 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7031 if (comm->bBondComm)
7033 /* Communicate atoms beyond the cut-off for bonded interactions */
7036 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7038 comm->bLocalCG = init_bLocalCG(mtop);
7042 /* Only communicate atoms based on cut-off */
7043 comm->cglink = NULL;
7044 comm->bLocalCG = NULL;
7048 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7050 gmx_bool bDynLoadBal, real dlb_scale,
7053 gmx_domdec_comm_t *comm;
7068 fprintf(fplog, "The maximum number of communication pulses is:");
7069 for (d = 0; d < dd->ndim; d++)
7071 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7073 fprintf(fplog, "\n");
7074 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7075 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7076 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7077 for (d = 0; d < DIM; d++)
7081 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7088 comm->cellsize_min_dlb[d]/
7089 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7091 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7094 fprintf(fplog, "\n");
7098 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7099 fprintf(fplog, "The initial number of communication pulses is:");
7100 for (d = 0; d < dd->ndim; d++)
7102 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7104 fprintf(fplog, "\n");
7105 fprintf(fplog, "The initial domain decomposition cell size is:");
7106 for (d = 0; d < DIM; d++)
7110 fprintf(fplog, " %c %.2f nm",
7111 dim2char(d), dd->comm->cellsize_min[d]);
7114 fprintf(fplog, "\n\n");
7117 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7119 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7120 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7121 "non-bonded interactions", "", comm->cutoff);
7125 limit = dd->comm->cellsize_limit;
7129 if (dynamic_dd_box(ddbox, ir))
7131 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7133 limit = dd->comm->cellsize_min[XX];
7134 for (d = 1; d < DIM; d++)
7136 limit = min(limit, dd->comm->cellsize_min[d]);
7140 if (comm->bInterCGBondeds)
7142 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7143 "two-body bonded interactions", "(-rdd)",
7144 max(comm->cutoff, comm->cutoff_mbody));
7145 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7146 "multi-body bonded interactions", "(-rdd)",
7147 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7151 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7152 "virtual site constructions", "(-rcon)", limit);
7154 if (dd->constraint_comm)
7156 sprintf(buf, "atoms separated by up to %d constraints",
7158 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7159 buf, "(-rcon)", limit);
7161 fprintf(fplog, "\n");
7167 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7169 const t_inputrec *ir,
7170 const gmx_ddbox_t *ddbox)
7172 gmx_domdec_comm_t *comm;
7173 int d, dim, npulse, npulse_d_max, npulse_d;
7178 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7180 /* Determine the maximum number of comm. pulses in one dimension */
7182 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7184 /* Determine the maximum required number of grid pulses */
7185 if (comm->cellsize_limit >= comm->cutoff)
7187 /* Only a single pulse is required */
7190 else if (!bNoCutOff && comm->cellsize_limit > 0)
7192 /* We round down slightly here to avoid overhead due to the latency
7193 * of extra communication calls when the cut-off
7194 * would be only slightly longer than the cell size.
7195 * Later cellsize_limit is redetermined,
7196 * so we can not miss interactions due to this rounding.
7198 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7202 /* There is no cell size limit */
7203 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7206 if (!bNoCutOff && npulse > 1)
7208 /* See if we can do with less pulses, based on dlb_scale */
7210 for (d = 0; d < dd->ndim; d++)
7213 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7214 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7215 npulse_d_max = max(npulse_d_max, npulse_d);
7217 npulse = min(npulse, npulse_d_max);
7220 /* This env var can override npulse */
7221 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7228 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7229 for (d = 0; d < dd->ndim; d++)
7231 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7232 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7233 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7234 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7235 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7237 comm->bVacDLBNoLimit = FALSE;
7241 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7242 if (!comm->bVacDLBNoLimit)
7244 comm->cellsize_limit = max(comm->cellsize_limit,
7245 comm->cutoff/comm->maxpulse);
7247 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7248 /* Set the minimum cell size for each DD dimension */
7249 for (d = 0; d < dd->ndim; d++)
7251 if (comm->bVacDLBNoLimit ||
7252 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7254 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7258 comm->cellsize_min_dlb[dd->dim[d]] =
7259 comm->cutoff/comm->cd[d].np_dlb;
7262 if (comm->cutoff_mbody <= 0)
7264 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7266 if (comm->bDynLoadBal)
7272 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7274 /* If each molecule is a single charge group
7275 * or we use domain decomposition for each periodic dimension,
7276 * we do not need to take pbc into account for the bonded interactions.
7278 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7281 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7284 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7285 t_inputrec *ir, t_forcerec *fr,
7288 gmx_domdec_comm_t *comm;
7294 /* Initialize the thread data.
7295 * This can not be done in init_domain_decomposition,
7296 * as the numbers of threads is determined later.
7298 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7301 snew(comm->dth, comm->nth);
7304 if (EEL_PME(ir->coulombtype))
7306 init_ddpme(dd, &comm->ddpme[0], 0);
7307 if (comm->npmedecompdim >= 2)
7309 init_ddpme(dd, &comm->ddpme[1], 1);
7314 comm->npmenodes = 0;
7315 if (dd->pme_nodeid >= 0)
7317 gmx_fatal_collective(FARGS, NULL, dd,
7318 "Can not have separate PME nodes without PME electrostatics");
7324 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7326 if (comm->eDLB != edlbNO)
7328 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7331 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7332 if (comm->eDLB == edlbAUTO)
7336 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7338 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7341 if (ir->ePBC == epbcNONE)
7343 vol_frac = 1 - 1/(double)dd->nnodes;
7348 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7352 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7354 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7356 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7359 static gmx_bool test_dd_cutoff(t_commrec *cr,
7360 t_state *state, t_inputrec *ir,
7371 set_ddbox(dd, FALSE, cr, ir, state->box,
7372 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7376 for (d = 0; d < dd->ndim; d++)
7380 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7381 if (dynamic_dd_box(&ddbox, ir))
7383 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7386 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7388 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7389 dd->comm->cd[d].np_dlb > 0)
7391 if (np > dd->comm->cd[d].np_dlb)
7396 /* If a current local cell size is smaller than the requested
7397 * cut-off, we could still fix it, but this gets very complicated.
7398 * Without fixing here, we might actually need more checks.
7400 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7407 if (dd->comm->eDLB != edlbNO)
7409 /* If DLB is not active yet, we don't need to check the grid jumps.
7410 * Actually we shouldn't, because then the grid jump data is not set.
7412 if (dd->comm->bDynLoadBal &&
7413 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7418 gmx_sumi(1, &LocallyLimited, cr);
7420 if (LocallyLimited > 0)
7429 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7432 gmx_bool bCutoffAllowed;
7434 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7438 cr->dd->comm->cutoff = cutoff_req;
7441 return bCutoffAllowed;
7444 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7446 gmx_domdec_comm_t *comm;
7448 comm = cr->dd->comm;
7450 /* Turn on the DLB limiting (might have been on already) */
7451 comm->bPMELoadBalDLBLimits = TRUE;
7453 /* Change the cut-off limit */
7454 comm->PMELoadBal_max_cutoff = comm->cutoff;
7457 static void merge_cg_buffers(int ncell,
7458 gmx_domdec_comm_dim_t *cd, int pulse,
7460 int *index_gl, int *recv_i,
7461 rvec *cg_cm, rvec *recv_vr,
7463 cginfo_mb_t *cginfo_mb, int *cginfo)
7465 gmx_domdec_ind_t *ind, *ind_p;
7466 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7467 int shift, shift_at;
7469 ind = &cd->ind[pulse];
7471 /* First correct the already stored data */
7472 shift = ind->nrecv[ncell];
7473 for (cell = ncell-1; cell >= 0; cell--)
7475 shift -= ind->nrecv[cell];
7478 /* Move the cg's present from previous grid pulses */
7479 cg0 = ncg_cell[ncell+cell];
7480 cg1 = ncg_cell[ncell+cell+1];
7481 cgindex[cg1+shift] = cgindex[cg1];
7482 for (cg = cg1-1; cg >= cg0; cg--)
7484 index_gl[cg+shift] = index_gl[cg];
7485 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7486 cgindex[cg+shift] = cgindex[cg];
7487 cginfo[cg+shift] = cginfo[cg];
7489 /* Correct the already stored send indices for the shift */
7490 for (p = 1; p <= pulse; p++)
7492 ind_p = &cd->ind[p];
7494 for (c = 0; c < cell; c++)
7496 cg0 += ind_p->nsend[c];
7498 cg1 = cg0 + ind_p->nsend[cell];
7499 for (cg = cg0; cg < cg1; cg++)
7501 ind_p->index[cg] += shift;
7507 /* Merge in the communicated buffers */
7511 for (cell = 0; cell < ncell; cell++)
7513 cg1 = ncg_cell[ncell+cell+1] + shift;
7516 /* Correct the old cg indices */
7517 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7519 cgindex[cg+1] += shift_at;
7522 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7524 /* Copy this charge group from the buffer */
7525 index_gl[cg1] = recv_i[cg0];
7526 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7527 /* Add it to the cgindex */
7528 cg_gl = index_gl[cg1];
7529 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7530 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7531 cgindex[cg1+1] = cgindex[cg1] + nat;
7536 shift += ind->nrecv[cell];
7537 ncg_cell[ncell+cell+1] = cg1;
7541 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7542 int nzone, int cg0, const int *cgindex)
7546 /* Store the atom block boundaries for easy copying of communication buffers
7549 for (zone = 0; zone < nzone; zone++)
7551 for (p = 0; p < cd->np; p++)
7553 cd->ind[p].cell2at0[zone] = cgindex[cg];
7554 cg += cd->ind[p].nrecv[zone];
7555 cd->ind[p].cell2at1[zone] = cgindex[cg];
7560 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7566 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7568 if (!bLocalCG[link->a[i]])
7577 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7579 real c[DIM][4]; /* the corners for the non-bonded communication */
7580 real cr0; /* corner for rounding */
7581 real cr1[4]; /* corners for rounding */
7582 real bc[DIM]; /* corners for bounded communication */
7583 real bcr1; /* corner for rounding for bonded communication */
7586 /* Determine the corners of the domain(s) we are communicating with */
7588 set_dd_corners(const gmx_domdec_t *dd,
7589 int dim0, int dim1, int dim2,
7593 const gmx_domdec_comm_t *comm;
7594 const gmx_domdec_zones_t *zones;
7599 zones = &comm->zones;
7601 /* Keep the compiler happy */
7605 /* The first dimension is equal for all cells */
7606 c->c[0][0] = comm->cell_x0[dim0];
7609 c->bc[0] = c->c[0][0];
7614 /* This cell row is only seen from the first row */
7615 c->c[1][0] = comm->cell_x0[dim1];
7616 /* All rows can see this row */
7617 c->c[1][1] = comm->cell_x0[dim1];
7620 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7623 /* For the multi-body distance we need the maximum */
7624 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7627 /* Set the upper-right corner for rounding */
7628 c->cr0 = comm->cell_x1[dim0];
7633 for (j = 0; j < 4; j++)
7635 c->c[2][j] = comm->cell_x0[dim2];
7639 /* Use the maximum of the i-cells that see a j-cell */
7640 for (i = 0; i < zones->nizone; i++)
7642 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7648 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7654 /* For the multi-body distance we need the maximum */
7655 c->bc[2] = comm->cell_x0[dim2];
7656 for (i = 0; i < 2; i++)
7658 for (j = 0; j < 2; j++)
7660 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7666 /* Set the upper-right corner for rounding */
7667 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7668 * Only cell (0,0,0) can see cell 7 (1,1,1)
7670 c->cr1[0] = comm->cell_x1[dim1];
7671 c->cr1[3] = comm->cell_x1[dim1];
7674 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7677 /* For the multi-body distance we need the maximum */
7678 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7685 /* Determine which cg's we need to send in this pulse from this zone */
7687 get_zone_pulse_cgs(gmx_domdec_t *dd,
7688 int zonei, int zone,
7690 const int *index_gl,
7692 int dim, int dim_ind,
7693 int dim0, int dim1, int dim2,
7694 real r_comm2, real r_bcomm2,
7698 real skew_fac2_d, real skew_fac_01,
7699 rvec *v_d, rvec *v_0, rvec *v_1,
7700 const dd_corners_t *c,
7702 gmx_bool bDistBonded,
7708 gmx_domdec_ind_t *ind,
7709 int **ibuf, int *ibuf_nalloc,
7715 gmx_domdec_comm_t *comm;
7717 gmx_bool bDistMB_pulse;
7719 real r2, rb2, r, tric_sh;
7722 int nsend_z, nsend, nat;
7726 bScrew = (dd->bScrewPBC && dim == XX);
7728 bDistMB_pulse = (bDistMB && bDistBonded);
7734 for (cg = cg0; cg < cg1; cg++)
7738 if (tric_dist[dim_ind] == 0)
7740 /* Rectangular direction, easy */
7741 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7748 r = cg_cm[cg][dim] - c->bc[dim_ind];
7754 /* Rounding gives at most a 16% reduction
7755 * in communicated atoms
7757 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7759 r = cg_cm[cg][dim0] - c->cr0;
7760 /* This is the first dimension, so always r >= 0 */
7767 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7769 r = cg_cm[cg][dim1] - c->cr1[zone];
7776 r = cg_cm[cg][dim1] - c->bcr1;
7786 /* Triclinic direction, more complicated */
7789 /* Rounding, conservative as the skew_fac multiplication
7790 * will slightly underestimate the distance.
7792 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7794 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7795 for (i = dim0+1; i < DIM; i++)
7797 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7799 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7802 rb[dim0] = rn[dim0];
7805 /* Take care that the cell planes along dim0 might not
7806 * be orthogonal to those along dim1 and dim2.
7808 for (i = 1; i <= dim_ind; i++)
7811 if (normal[dim0][dimd] > 0)
7813 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7816 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7821 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7823 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7825 for (i = dim1+1; i < DIM; i++)
7827 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7829 rn[dim1] += tric_sh;
7832 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7833 /* Take care of coupling of the distances
7834 * to the planes along dim0 and dim1 through dim2.
7836 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7837 /* Take care that the cell planes along dim1
7838 * might not be orthogonal to that along dim2.
7840 if (normal[dim1][dim2] > 0)
7842 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7848 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7851 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7852 /* Take care of coupling of the distances
7853 * to the planes along dim0 and dim1 through dim2.
7855 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7856 /* Take care that the cell planes along dim1
7857 * might not be orthogonal to that along dim2.
7859 if (normal[dim1][dim2] > 0)
7861 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7866 /* The distance along the communication direction */
7867 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7869 for (i = dim+1; i < DIM; i++)
7871 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7876 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7877 /* Take care of coupling of the distances
7878 * to the planes along dim0 and dim1 through dim2.
7880 if (dim_ind == 1 && zonei == 1)
7882 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7888 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7891 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7892 /* Take care of coupling of the distances
7893 * to the planes along dim0 and dim1 through dim2.
7895 if (dim_ind == 1 && zonei == 1)
7897 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7905 ((bDistMB && rb2 < r_bcomm2) ||
7906 (bDist2B && r2 < r_bcomm2)) &&
7908 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7909 missing_link(comm->cglink, index_gl[cg],
7912 /* Make an index to the local charge groups */
7913 if (nsend+1 > ind->nalloc)
7915 ind->nalloc = over_alloc_large(nsend+1);
7916 srenew(ind->index, ind->nalloc);
7918 if (nsend+1 > *ibuf_nalloc)
7920 *ibuf_nalloc = over_alloc_large(nsend+1);
7921 srenew(*ibuf, *ibuf_nalloc);
7923 ind->index[nsend] = cg;
7924 (*ibuf)[nsend] = index_gl[cg];
7926 vec_rvec_check_alloc(vbuf, nsend+1);
7928 if (dd->ci[dim] == 0)
7930 /* Correct cg_cm for pbc */
7931 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7934 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7935 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7940 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7943 nat += cgindex[cg+1] - cgindex[cg];
7949 *nsend_z_ptr = nsend_z;
7952 static void setup_dd_communication(gmx_domdec_t *dd,
7953 matrix box, gmx_ddbox_t *ddbox,
7954 t_forcerec *fr, t_state *state, rvec **f)
7956 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7957 int nzone, nzone_send, zone, zonei, cg0, cg1;
7958 int c, i, j, cg, cg_gl, nrcg;
7959 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7960 gmx_domdec_comm_t *comm;
7961 gmx_domdec_zones_t *zones;
7962 gmx_domdec_comm_dim_t *cd;
7963 gmx_domdec_ind_t *ind;
7964 cginfo_mb_t *cginfo_mb;
7965 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7966 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7967 dd_corners_t corners;
7969 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7970 real skew_fac2_d, skew_fac_01;
7977 fprintf(debug, "Setting up DD communication\n");
7982 switch (fr->cutoff_scheme)
7991 gmx_incons("unimplemented");
7995 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
7997 dim = dd->dim[dim_ind];
7999 /* Check if we need to use triclinic distances */
8000 tric_dist[dim_ind] = 0;
8001 for (i = 0; i <= dim_ind; i++)
8003 if (ddbox->tric_dir[dd->dim[i]])
8005 tric_dist[dim_ind] = 1;
8010 bBondComm = comm->bBondComm;
8012 /* Do we need to determine extra distances for multi-body bondeds? */
8013 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8015 /* Do we need to determine extra distances for only two-body bondeds? */
8016 bDist2B = (bBondComm && !bDistMB);
8018 r_comm2 = sqr(comm->cutoff);
8019 r_bcomm2 = sqr(comm->cutoff_mbody);
8023 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8026 zones = &comm->zones;
8029 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8030 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8032 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8034 /* Triclinic stuff */
8035 normal = ddbox->normal;
8039 v_0 = ddbox->v[dim0];
8040 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8042 /* Determine the coupling coefficient for the distances
8043 * to the cell planes along dim0 and dim1 through dim2.
8044 * This is required for correct rounding.
8047 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8050 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8056 v_1 = ddbox->v[dim1];
8059 zone_cg_range = zones->cg_range;
8060 index_gl = dd->index_gl;
8061 cgindex = dd->cgindex;
8062 cginfo_mb = fr->cginfo_mb;
8064 zone_cg_range[0] = 0;
8065 zone_cg_range[1] = dd->ncg_home;
8066 comm->zone_ncg1[0] = dd->ncg_home;
8067 pos_cg = dd->ncg_home;
8069 nat_tot = dd->nat_home;
8071 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8073 dim = dd->dim[dim_ind];
8074 cd = &comm->cd[dim_ind];
8076 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8078 /* No pbc in this dimension, the first node should not comm. */
8086 v_d = ddbox->v[dim];
8087 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8089 cd->bInPlace = TRUE;
8090 for (p = 0; p < cd->np; p++)
8092 /* Only atoms communicated in the first pulse are used
8093 * for multi-body bonded interactions or for bBondComm.
8095 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8100 for (zone = 0; zone < nzone_send; zone++)
8102 if (tric_dist[dim_ind] && dim_ind > 0)
8104 /* Determine slightly more optimized skew_fac's
8106 * This reduces the number of communicated atoms
8107 * by about 10% for 3D DD of rhombic dodecahedra.
8109 for (dimd = 0; dimd < dim; dimd++)
8111 sf2_round[dimd] = 1;
8112 if (ddbox->tric_dir[dimd])
8114 for (i = dd->dim[dimd]+1; i < DIM; i++)
8116 /* If we are shifted in dimension i
8117 * and the cell plane is tilted forward
8118 * in dimension i, skip this coupling.
8120 if (!(zones->shift[nzone+zone][i] &&
8121 ddbox->v[dimd][i][dimd] >= 0))
8124 sqr(ddbox->v[dimd][i][dimd]);
8127 sf2_round[dimd] = 1/sf2_round[dimd];
8132 zonei = zone_perm[dim_ind][zone];
8135 /* Here we permutate the zones to obtain a convenient order
8136 * for neighbor searching
8138 cg0 = zone_cg_range[zonei];
8139 cg1 = zone_cg_range[zonei+1];
8143 /* Look only at the cg's received in the previous grid pulse
8145 cg1 = zone_cg_range[nzone+zone+1];
8146 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8149 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8150 for (th = 0; th < comm->nth; th++)
8152 gmx_domdec_ind_t *ind_p;
8153 int **ibuf_p, *ibuf_nalloc_p;
8155 int *nsend_p, *nat_p;
8161 /* Thread 0 writes in the comm buffers */
8163 ibuf_p = &comm->buf_int;
8164 ibuf_nalloc_p = &comm->nalloc_int;
8165 vbuf_p = &comm->vbuf;
8168 nsend_zone_p = &ind->nsend[zone];
8172 /* Other threads write into temp buffers */
8173 ind_p = &comm->dth[th].ind;
8174 ibuf_p = &comm->dth[th].ibuf;
8175 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8176 vbuf_p = &comm->dth[th].vbuf;
8177 nsend_p = &comm->dth[th].nsend;
8178 nat_p = &comm->dth[th].nat;
8179 nsend_zone_p = &comm->dth[th].nsend_zone;
8181 comm->dth[th].nsend = 0;
8182 comm->dth[th].nat = 0;
8183 comm->dth[th].nsend_zone = 0;
8193 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8194 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8197 /* Get the cg's for this pulse in this zone */
8198 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8200 dim, dim_ind, dim0, dim1, dim2,
8203 normal, skew_fac2_d, skew_fac_01,
8204 v_d, v_0, v_1, &corners, sf2_round,
8205 bDistBonded, bBondComm,
8209 ibuf_p, ibuf_nalloc_p,
8215 /* Append data of threads>=1 to the communication buffers */
8216 for (th = 1; th < comm->nth; th++)
8218 dd_comm_setup_work_t *dth;
8221 dth = &comm->dth[th];
8223 ns1 = nsend + dth->nsend_zone;
8224 if (ns1 > ind->nalloc)
8226 ind->nalloc = over_alloc_dd(ns1);
8227 srenew(ind->index, ind->nalloc);
8229 if (ns1 > comm->nalloc_int)
8231 comm->nalloc_int = over_alloc_dd(ns1);
8232 srenew(comm->buf_int, comm->nalloc_int);
8234 if (ns1 > comm->vbuf.nalloc)
8236 comm->vbuf.nalloc = over_alloc_dd(ns1);
8237 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8240 for (i = 0; i < dth->nsend_zone; i++)
8242 ind->index[nsend] = dth->ind.index[i];
8243 comm->buf_int[nsend] = dth->ibuf[i];
8244 copy_rvec(dth->vbuf.v[i],
8245 comm->vbuf.v[nsend]);
8249 ind->nsend[zone] += dth->nsend_zone;
8252 /* Clear the counts in case we do not have pbc */
8253 for (zone = nzone_send; zone < nzone; zone++)
8255 ind->nsend[zone] = 0;
8257 ind->nsend[nzone] = nsend;
8258 ind->nsend[nzone+1] = nat;
8259 /* Communicate the number of cg's and atoms to receive */
8260 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8261 ind->nsend, nzone+2,
8262 ind->nrecv, nzone+2);
8264 /* The rvec buffer is also required for atom buffers of size nsend
8265 * in dd_move_x and dd_move_f.
8267 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8271 /* We can receive in place if only the last zone is not empty */
8272 for (zone = 0; zone < nzone-1; zone++)
8274 if (ind->nrecv[zone] > 0)
8276 cd->bInPlace = FALSE;
8281 /* The int buffer is only required here for the cg indices */
8282 if (ind->nrecv[nzone] > comm->nalloc_int2)
8284 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8285 srenew(comm->buf_int2, comm->nalloc_int2);
8287 /* The rvec buffer is also required for atom buffers
8288 * of size nrecv in dd_move_x and dd_move_f.
8290 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8291 vec_rvec_check_alloc(&comm->vbuf2, i);
8295 /* Make space for the global cg indices */
8296 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8297 || dd->cg_nalloc == 0)
8299 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8300 srenew(index_gl, dd->cg_nalloc);
8301 srenew(cgindex, dd->cg_nalloc+1);
8303 /* Communicate the global cg indices */
8306 recv_i = index_gl + pos_cg;
8310 recv_i = comm->buf_int2;
8312 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8313 comm->buf_int, nsend,
8314 recv_i, ind->nrecv[nzone]);
8316 /* Make space for cg_cm */
8317 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8318 if (fr->cutoff_scheme == ecutsGROUP)
8326 /* Communicate cg_cm */
8329 recv_vr = cg_cm + pos_cg;
8333 recv_vr = comm->vbuf2.v;
8335 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8336 comm->vbuf.v, nsend,
8337 recv_vr, ind->nrecv[nzone]);
8339 /* Make the charge group index */
8342 zone = (p == 0 ? 0 : nzone - 1);
8343 while (zone < nzone)
8345 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8347 cg_gl = index_gl[pos_cg];
8348 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8349 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8350 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8353 /* Update the charge group presence,
8354 * so we can use it in the next pass of the loop.
8356 comm->bLocalCG[cg_gl] = TRUE;
8362 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8365 zone_cg_range[nzone+zone] = pos_cg;
8370 /* This part of the code is never executed with bBondComm. */
8371 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8372 index_gl, recv_i, cg_cm, recv_vr,
8373 cgindex, fr->cginfo_mb, fr->cginfo);
8374 pos_cg += ind->nrecv[nzone];
8376 nat_tot += ind->nrecv[nzone+1];
8380 /* Store the atom block for easy copying of communication buffers */
8381 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8385 dd->index_gl = index_gl;
8386 dd->cgindex = cgindex;
8388 dd->ncg_tot = zone_cg_range[zones->n];
8389 dd->nat_tot = nat_tot;
8390 comm->nat[ddnatHOME] = dd->nat_home;
8391 for (i = ddnatZONE; i < ddnatNR; i++)
8393 comm->nat[i] = dd->nat_tot;
8398 /* We don't need to update cginfo, since that was alrady done above.
8399 * So we pass NULL for the forcerec.
8401 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8402 NULL, comm->bLocalCG);
8407 fprintf(debug, "Finished setting up DD communication, zones:");
8408 for (c = 0; c < zones->n; c++)
8410 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8412 fprintf(debug, "\n");
8416 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8420 for (c = 0; c < zones->nizone; c++)
8422 zones->izone[c].cg1 = zones->cg_range[c+1];
8423 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8424 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8428 static void set_zones_size(gmx_domdec_t *dd,
8429 matrix box, const gmx_ddbox_t *ddbox,
8430 int zone_start, int zone_end)
8432 gmx_domdec_comm_t *comm;
8433 gmx_domdec_zones_t *zones;
8435 int z, zi, zj0, zj1, d, dim;
8438 real size_j, add_tric;
8443 zones = &comm->zones;
8445 /* Do we need to determine extra distances for multi-body bondeds? */
8446 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8448 for (z = zone_start; z < zone_end; z++)
8450 /* Copy cell limits to zone limits.
8451 * Valid for non-DD dims and non-shifted dims.
8453 copy_rvec(comm->cell_x0, zones->size[z].x0);
8454 copy_rvec(comm->cell_x1, zones->size[z].x1);
8457 for (d = 0; d < dd->ndim; d++)
8461 for (z = 0; z < zones->n; z++)
8463 /* With a staggered grid we have different sizes
8464 * for non-shifted dimensions.
8466 if (dd->bGridJump && zones->shift[z][dim] == 0)
8470 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8471 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8475 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8476 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8482 rcmbs = comm->cutoff_mbody;
8483 if (ddbox->tric_dir[dim])
8485 rcs /= ddbox->skew_fac[dim];
8486 rcmbs /= ddbox->skew_fac[dim];
8489 /* Set the lower limit for the shifted zone dimensions */
8490 for (z = zone_start; z < zone_end; z++)
8492 if (zones->shift[z][dim] > 0)
8495 if (!dd->bGridJump || d == 0)
8497 zones->size[z].x0[dim] = comm->cell_x1[dim];
8498 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8502 /* Here we take the lower limit of the zone from
8503 * the lowest domain of the zone below.
8507 zones->size[z].x0[dim] =
8508 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8514 zones->size[z].x0[dim] =
8515 zones->size[zone_perm[2][z-4]].x0[dim];
8519 zones->size[z].x0[dim] =
8520 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8523 /* A temporary limit, is updated below */
8524 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8528 for (zi = 0; zi < zones->nizone; zi++)
8530 if (zones->shift[zi][dim] == 0)
8532 /* This takes the whole zone into account.
8533 * With multiple pulses this will lead
8534 * to a larger zone then strictly necessary.
8536 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8537 zones->size[zi].x1[dim]+rcmbs);
8545 /* Loop over the i-zones to set the upper limit of each
8548 for (zi = 0; zi < zones->nizone; zi++)
8550 if (zones->shift[zi][dim] == 0)
8552 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8554 if (zones->shift[z][dim] > 0)
8556 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8557 zones->size[zi].x1[dim]+rcs);
8564 for (z = zone_start; z < zone_end; z++)
8566 /* Initialization only required to keep the compiler happy */
8567 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8570 /* To determine the bounding box for a zone we need to find
8571 * the extreme corners of 4, 2 or 1 corners.
8573 nc = 1 << (ddbox->npbcdim - 1);
8575 for (c = 0; c < nc; c++)
8577 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8581 corner[YY] = zones->size[z].x0[YY];
8585 corner[YY] = zones->size[z].x1[YY];
8589 corner[ZZ] = zones->size[z].x0[ZZ];
8593 corner[ZZ] = zones->size[z].x1[ZZ];
8595 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8597 /* With 1D domain decomposition the cg's are not in
8598 * the triclinic box, but triclinic x-y and rectangular y-z.
8599 * Shift y back, so it will later end up at 0.
8601 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8603 /* Apply the triclinic couplings */
8604 for (i = YY; i < ddbox->npbcdim; i++)
8606 for (j = XX; j < i; j++)
8608 corner[j] += corner[i]*box[i][j]/box[i][i];
8613 copy_rvec(corner, corner_min);
8614 copy_rvec(corner, corner_max);
8618 for (i = 0; i < DIM; i++)
8620 corner_min[i] = min(corner_min[i], corner[i]);
8621 corner_max[i] = max(corner_max[i], corner[i]);
8625 /* Copy the extreme cornes without offset along x */
8626 for (i = 0; i < DIM; i++)
8628 zones->size[z].bb_x0[i] = corner_min[i];
8629 zones->size[z].bb_x1[i] = corner_max[i];
8631 /* Add the offset along x */
8632 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8633 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8636 if (zone_start == 0)
8639 for (dim = 0; dim < DIM; dim++)
8641 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8643 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8648 for (z = zone_start; z < zone_end; z++)
8650 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8652 zones->size[z].x0[XX], zones->size[z].x1[XX],
8653 zones->size[z].x0[YY], zones->size[z].x1[YY],
8654 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8655 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8657 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8658 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8659 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8664 static int comp_cgsort(const void *a, const void *b)
8668 gmx_cgsort_t *cga, *cgb;
8669 cga = (gmx_cgsort_t *)a;
8670 cgb = (gmx_cgsort_t *)b;
8672 comp = cga->nsc - cgb->nsc;
8675 comp = cga->ind_gl - cgb->ind_gl;
8681 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8686 /* Order the data */
8687 for (i = 0; i < n; i++)
8689 buf[i] = a[sort[i].ind];
8692 /* Copy back to the original array */
8693 for (i = 0; i < n; i++)
8699 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8704 /* Order the data */
8705 for (i = 0; i < n; i++)
8707 copy_rvec(v[sort[i].ind], buf[i]);
8710 /* Copy back to the original array */
8711 for (i = 0; i < n; i++)
8713 copy_rvec(buf[i], v[i]);
8717 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8720 int a, atot, cg, cg0, cg1, i;
8722 if (cgindex == NULL)
8724 /* Avoid the useless loop of the atoms within a cg */
8725 order_vec_cg(ncg, sort, v, buf);
8730 /* Order the data */
8732 for (cg = 0; cg < ncg; cg++)
8734 cg0 = cgindex[sort[cg].ind];
8735 cg1 = cgindex[sort[cg].ind+1];
8736 for (i = cg0; i < cg1; i++)
8738 copy_rvec(v[i], buf[a]);
8744 /* Copy back to the original array */
8745 for (a = 0; a < atot; a++)
8747 copy_rvec(buf[a], v[a]);
8751 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8752 int nsort_new, gmx_cgsort_t *sort_new,
8753 gmx_cgsort_t *sort1)
8757 /* The new indices are not very ordered, so we qsort them */
8758 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8760 /* sort2 is already ordered, so now we can merge the two arrays */
8764 while (i2 < nsort2 || i_new < nsort_new)
8768 sort1[i1++] = sort_new[i_new++];
8770 else if (i_new == nsort_new)
8772 sort1[i1++] = sort2[i2++];
8774 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8775 (sort2[i2].nsc == sort_new[i_new].nsc &&
8776 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8778 sort1[i1++] = sort2[i2++];
8782 sort1[i1++] = sort_new[i_new++];
8787 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8789 gmx_domdec_sort_t *sort;
8790 gmx_cgsort_t *cgsort, *sort_i;
8791 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8792 int sort_last, sort_skip;
8794 sort = dd->comm->sort;
8796 a = fr->ns.grid->cell_index;
8798 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8800 if (ncg_home_old >= 0)
8802 /* The charge groups that remained in the same ns grid cell
8803 * are completely ordered. So we can sort efficiently by sorting
8804 * the charge groups that did move into the stationary list.
8809 for (i = 0; i < dd->ncg_home; i++)
8811 /* Check if this cg did not move to another node */
8814 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8816 /* This cg is new on this node or moved ns grid cell */
8817 if (nsort_new >= sort->sort_new_nalloc)
8819 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8820 srenew(sort->sort_new, sort->sort_new_nalloc);
8822 sort_i = &(sort->sort_new[nsort_new++]);
8826 /* This cg did not move */
8827 sort_i = &(sort->sort2[nsort2++]);
8829 /* Sort on the ns grid cell indices
8830 * and the global topology index.
8831 * index_gl is irrelevant with cell ns,
8832 * but we set it here anyhow to avoid a conditional.
8835 sort_i->ind_gl = dd->index_gl[i];
8842 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8845 /* Sort efficiently */
8846 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8851 cgsort = sort->sort;
8853 for (i = 0; i < dd->ncg_home; i++)
8855 /* Sort on the ns grid cell indices
8856 * and the global topology index
8858 cgsort[i].nsc = a[i];
8859 cgsort[i].ind_gl = dd->index_gl[i];
8861 if (cgsort[i].nsc < moved)
8868 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8870 /* Determine the order of the charge groups using qsort */
8871 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8877 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8880 int ncg_new, i, *a, na;
8882 sort = dd->comm->sort->sort;
8884 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8887 for (i = 0; i < na; i++)
8891 sort[ncg_new].ind = a[i];
8899 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8900 rvec *cgcm, t_forcerec *fr, t_state *state,
8903 gmx_domdec_sort_t *sort;
8904 gmx_cgsort_t *cgsort, *sort_i;
8906 int ncg_new, i, *ibuf, cgsize;
8909 sort = dd->comm->sort;
8911 if (dd->ncg_home > sort->sort_nalloc)
8913 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8914 srenew(sort->sort, sort->sort_nalloc);
8915 srenew(sort->sort2, sort->sort_nalloc);
8917 cgsort = sort->sort;
8919 switch (fr->cutoff_scheme)
8922 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8925 ncg_new = dd_sort_order_nbnxn(dd, fr);
8928 gmx_incons("unimplemented");
8932 /* We alloc with the old size, since cgindex is still old */
8933 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8934 vbuf = dd->comm->vbuf.v;
8938 cgindex = dd->cgindex;
8945 /* Remove the charge groups which are no longer at home here */
8946 dd->ncg_home = ncg_new;
8949 fprintf(debug, "Set the new home charge group count to %d\n",
8953 /* Reorder the state */
8954 for (i = 0; i < estNR; i++)
8956 if (EST_DISTR(i) && (state->flags & (1<<i)))
8961 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8964 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8967 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8970 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8974 case estDISRE_INITF:
8975 case estDISRE_RM3TAV:
8976 case estORIRE_INITF:
8978 /* No ordering required */
8981 gmx_incons("Unknown state entry encountered in dd_sort_state");
8986 if (fr->cutoff_scheme == ecutsGROUP)
8989 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
8992 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8994 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8995 srenew(sort->ibuf, sort->ibuf_nalloc);
8998 /* Reorder the global cg index */
8999 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9000 /* Reorder the cginfo */
9001 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9002 /* Rebuild the local cg index */
9006 for (i = 0; i < dd->ncg_home; i++)
9008 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9009 ibuf[i+1] = ibuf[i] + cgsize;
9011 for (i = 0; i < dd->ncg_home+1; i++)
9013 dd->cgindex[i] = ibuf[i];
9018 for (i = 0; i < dd->ncg_home+1; i++)
9023 /* Set the home atom number */
9024 dd->nat_home = dd->cgindex[dd->ncg_home];
9026 if (fr->cutoff_scheme == ecutsVERLET)
9028 /* The atoms are now exactly in grid order, update the grid order */
9029 nbnxn_set_atomorder(fr->nbv->nbs);
9033 /* Copy the sorted ns cell indices back to the ns grid struct */
9034 for (i = 0; i < dd->ncg_home; i++)
9036 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9038 fr->ns.grid->nr = dd->ncg_home;
9042 static void add_dd_statistics(gmx_domdec_t *dd)
9044 gmx_domdec_comm_t *comm;
9049 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9051 comm->sum_nat[ddnat-ddnatZONE] +=
9052 comm->nat[ddnat] - comm->nat[ddnat-1];
9057 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9059 gmx_domdec_comm_t *comm;
9064 /* Reset all the statistics and counters for total run counting */
9065 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9067 comm->sum_nat[ddnat-ddnatZONE] = 0;
9071 comm->load_step = 0;
9074 clear_ivec(comm->load_lim);
9079 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9081 gmx_domdec_comm_t *comm;
9085 comm = cr->dd->comm;
9087 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9094 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");
9096 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9098 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9103 " av. #atoms communicated per step for force: %d x %.1f\n",
9107 if (cr->dd->vsite_comm)
9110 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9111 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9116 if (cr->dd->constraint_comm)
9119 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9120 1 + ir->nLincsIter, av);
9124 gmx_incons(" Unknown type for DD statistics");
9127 fprintf(fplog, "\n");
9129 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9131 print_dd_load_av(fplog, cr->dd);
9135 void dd_partition_system(FILE *fplog,
9136 gmx_large_int_t step,
9138 gmx_bool bMasterState,
9140 t_state *state_global,
9141 gmx_mtop_t *top_global,
9143 t_state *state_local,
9146 gmx_localtop_t *top_local,
9149 gmx_shellfc_t shellfc,
9150 gmx_constr_t constr,
9152 gmx_wallcycle_t wcycle,
9156 gmx_domdec_comm_t *comm;
9157 gmx_ddbox_t ddbox = {0};
9159 gmx_large_int_t step_pcoupl;
9160 rvec cell_ns_x0, cell_ns_x1;
9161 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9162 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9163 gmx_bool bRedist, bSortCG, bResortAll;
9164 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9171 bBoxChanged = (bMasterState || DEFORM(*ir));
9172 if (ir->epc != epcNO)
9174 /* With nstpcouple > 1 pressure coupling happens.
9175 * one step after calculating the pressure.
9176 * Box scaling happens at the end of the MD step,
9177 * after the DD partitioning.
9178 * We therefore have to do DLB in the first partitioning
9179 * after an MD step where P-coupling occured.
9180 * We need to determine the last step in which p-coupling occurred.
9181 * MRS -- need to validate this for vv?
9186 step_pcoupl = step - 1;
9190 step_pcoupl = ((step - 1)/n)*n + 1;
9192 if (step_pcoupl >= comm->partition_step)
9198 bNStGlobalComm = (step % nstglobalcomm == 0);
9200 if (!comm->bDynLoadBal)
9206 /* Should we do dynamic load balacing this step?
9207 * Since it requires (possibly expensive) global communication,
9208 * we might want to do DLB less frequently.
9210 if (bBoxChanged || ir->epc != epcNO)
9212 bDoDLB = bBoxChanged;
9216 bDoDLB = bNStGlobalComm;
9220 /* Check if we have recorded loads on the nodes */
9221 if (comm->bRecordLoad && dd_load_count(comm))
9223 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9225 /* Check if we should use DLB at the second partitioning
9226 * and every 100 partitionings,
9227 * so the extra communication cost is negligible.
9229 n = max(100, nstglobalcomm);
9230 bCheckDLB = (comm->n_load_collect == 0 ||
9231 comm->n_load_have % n == n-1);
9238 /* Print load every nstlog, first and last step to the log file */
9239 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9240 comm->n_load_collect == 0 ||
9242 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9244 /* Avoid extra communication due to verbose screen output
9245 * when nstglobalcomm is set.
9247 if (bDoDLB || bLogLoad || bCheckDLB ||
9248 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9250 get_load_distribution(dd, wcycle);
9255 dd_print_load(fplog, dd, step-1);
9259 dd_print_load_verbose(dd);
9262 comm->n_load_collect++;
9266 /* Since the timings are node dependent, the master decides */
9270 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9273 fprintf(debug, "step %s, imb loss %f\n",
9274 gmx_step_str(step, sbuf),
9275 dd_force_imb_perf_loss(dd));
9278 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9281 turn_on_dlb(fplog, cr, step);
9286 comm->n_load_have++;
9289 cgs_gl = &comm->cgs_gl;
9294 /* Clear the old state */
9295 clear_dd_indices(dd, 0, 0);
9297 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9298 TRUE, cgs_gl, state_global->x, &ddbox);
9300 get_cg_distribution(fplog, step, dd, cgs_gl,
9301 state_global->box, &ddbox, state_global->x);
9303 dd_distribute_state(dd, cgs_gl,
9304 state_global, state_local, f);
9306 dd_make_local_cgs(dd, &top_local->cgs);
9308 /* Ensure that we have space for the new distribution */
9309 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9311 if (fr->cutoff_scheme == ecutsGROUP)
9313 calc_cgcm(fplog, 0, dd->ncg_home,
9314 &top_local->cgs, state_local->x, fr->cg_cm);
9317 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9319 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9323 else if (state_local->ddp_count != dd->ddp_count)
9325 if (state_local->ddp_count > dd->ddp_count)
9327 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9330 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9332 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);
9335 /* Clear the old state */
9336 clear_dd_indices(dd, 0, 0);
9338 /* Build the new indices */
9339 rebuild_cgindex(dd, cgs_gl->index, state_local);
9340 make_dd_indices(dd, cgs_gl->index, 0);
9342 if (fr->cutoff_scheme == ecutsGROUP)
9344 /* Redetermine the cg COMs */
9345 calc_cgcm(fplog, 0, dd->ncg_home,
9346 &top_local->cgs, state_local->x, fr->cg_cm);
9349 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9351 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9353 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9354 TRUE, &top_local->cgs, state_local->x, &ddbox);
9356 bRedist = comm->bDynLoadBal;
9360 /* We have the full state, only redistribute the cgs */
9362 /* Clear the non-home indices */
9363 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9365 /* Avoid global communication for dim's without pbc and -gcom */
9366 if (!bNStGlobalComm)
9368 copy_rvec(comm->box0, ddbox.box0 );
9369 copy_rvec(comm->box_size, ddbox.box_size);
9371 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9372 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9377 /* For dim's without pbc and -gcom */
9378 copy_rvec(ddbox.box0, comm->box0 );
9379 copy_rvec(ddbox.box_size, comm->box_size);
9381 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9384 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9386 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9389 /* Check if we should sort the charge groups */
9390 if (comm->nstSortCG > 0)
9392 bSortCG = (bMasterState ||
9393 (bRedist && (step % comm->nstSortCG == 0)));
9400 ncg_home_old = dd->ncg_home;
9405 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9407 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9408 state_local, f, fr, mdatoms,
9409 !bSortCG, nrnb, &cg0, &ncg_moved);
9411 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9414 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9416 &comm->cell_x0, &comm->cell_x1,
9417 dd->ncg_home, fr->cg_cm,
9418 cell_ns_x0, cell_ns_x1, &grid_density);
9422 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9425 switch (fr->cutoff_scheme)
9428 copy_ivec(fr->ns.grid->n, ncells_old);
9429 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9430 state_local->box, cell_ns_x0, cell_ns_x1,
9431 fr->rlistlong, grid_density);
9434 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9437 gmx_incons("unimplemented");
9439 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9440 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9444 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9446 /* Sort the state on charge group position.
9447 * This enables exact restarts from this step.
9448 * It also improves performance by about 15% with larger numbers
9449 * of atoms per node.
9452 /* Fill the ns grid with the home cell,
9453 * so we can sort with the indices.
9455 set_zones_ncg_home(dd);
9457 switch (fr->cutoff_scheme)
9460 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9462 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9464 comm->zones.size[0].bb_x0,
9465 comm->zones.size[0].bb_x1,
9467 comm->zones.dens_zone0,
9470 ncg_moved, bRedist ? comm->moved : NULL,
9471 fr->nbv->grp[eintLocal].kernel_type,
9472 fr->nbv->grp[eintLocal].nbat);
9474 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9477 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9478 0, dd->ncg_home, fr->cg_cm);
9480 copy_ivec(fr->ns.grid->n, ncells_new);
9483 gmx_incons("unimplemented");
9486 bResortAll = bMasterState;
9488 /* Check if we can user the old order and ns grid cell indices
9489 * of the charge groups to sort the charge groups efficiently.
9491 if (ncells_new[XX] != ncells_old[XX] ||
9492 ncells_new[YY] != ncells_old[YY] ||
9493 ncells_new[ZZ] != ncells_old[ZZ])
9500 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9501 gmx_step_str(step, sbuf), dd->ncg_home);
9503 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9504 bResortAll ? -1 : ncg_home_old);
9505 /* Rebuild all the indices */
9507 ga2la_clear(dd->ga2la);
9509 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9512 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9514 /* Setup up the communication and communicate the coordinates */
9515 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9517 /* Set the indices */
9518 make_dd_indices(dd, cgs_gl->index, cg0);
9520 /* Set the charge group boundaries for neighbor searching */
9521 set_cg_boundaries(&comm->zones);
9523 if (fr->cutoff_scheme == ecutsVERLET)
9525 set_zones_size(dd, state_local->box, &ddbox,
9526 bSortCG ? 1 : 0, comm->zones.n);
9529 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9532 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9533 -1,state_local->x,state_local->box);
9536 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9538 /* Extract a local topology from the global topology */
9539 for (i = 0; i < dd->ndim; i++)
9541 np[dd->dim[i]] = comm->cd[i].np;
9543 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9544 comm->cellsize_min, np,
9546 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9547 vsite, top_global, top_local);
9549 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9551 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9553 /* Set up the special atom communication */
9554 n = comm->nat[ddnatZONE];
9555 for (i = ddnatZONE+1; i < ddnatNR; i++)
9560 if (vsite && vsite->n_intercg_vsite)
9562 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9566 if (dd->bInterCGcons || dd->bInterCGsettles)
9568 /* Only for inter-cg constraints we need special code */
9569 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9570 constr, ir->nProjOrder,
9571 top_local->idef.il);
9575 gmx_incons("Unknown special atom type setup");
9580 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9582 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9584 /* Make space for the extra coordinates for virtual site
9585 * or constraint communication.
9587 state_local->natoms = comm->nat[ddnatNR-1];
9588 if (state_local->natoms > state_local->nalloc)
9590 dd_realloc_state(state_local, f, state_local->natoms);
9593 if (fr->bF_NoVirSum)
9595 if (vsite && vsite->n_intercg_vsite)
9597 nat_f_novirsum = comm->nat[ddnatVSITE];
9601 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9603 nat_f_novirsum = dd->nat_tot;
9607 nat_f_novirsum = dd->nat_home;
9616 /* Set the number of atoms required for the force calculation.
9617 * Forces need to be constrained when using a twin-range setup
9618 * or with energy minimization. For simple simulations we could
9619 * avoid some allocation, zeroing and copying, but this is
9620 * probably not worth the complications ande checking.
9622 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9623 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9625 /* We make the all mdatoms up to nat_tot_con.
9626 * We could save some work by only setting invmass
9627 * between nat_tot and nat_tot_con.
9629 /* This call also sets the new number of home particles to dd->nat_home */
9630 atoms2md(top_global, ir,
9631 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9633 /* Now we have the charges we can sort the FE interactions */
9634 dd_sort_local_top(dd, mdatoms, top_local);
9638 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9639 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9644 /* Make the local shell stuff, currently no communication is done */
9645 make_local_shells(cr, mdatoms, shellfc);
9648 if (ir->implicit_solvent)
9650 make_local_gb(cr, fr->born, ir->gb_algorithm);
9653 init_bonded_thread_force_reduction(fr, &top_local->idef);
9655 if (!(cr->duty & DUTY_PME))
9657 /* Send the charges to our PME only node */
9658 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9659 mdatoms->chargeA, mdatoms->chargeB,
9660 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9665 set_constraints(constr, top_local, ir, mdatoms, cr);
9668 if (ir->ePull != epullNO)
9670 /* Update the local pull groups */
9671 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9676 /* Update the local rotation groups */
9677 dd_make_local_rotation_groups(dd, ir->rot);
9681 add_dd_statistics(dd);
9683 /* Make sure we only count the cycles for this DD partitioning */
9684 clear_dd_cycle_counts(dd);
9686 /* Because the order of the atoms might have changed since
9687 * the last vsite construction, we need to communicate the constructing
9688 * atom coordinates again (for spreading the forces this MD step).
9690 dd_move_x_vsites(dd, state_local->box, state_local->x);
9692 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9694 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9696 dd_move_x(dd, state_local->box, state_local->x);
9697 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9698 -1, state_local->x, state_local->box);
9701 /* Store the partitioning step */
9702 comm->partition_step = step;
9704 /* Increase the DD partitioning counter */
9706 /* The state currently matches this DD partitioning count, store it */
9707 state_local->ddp_count = dd->ddp_count;
9710 /* The DD master node knows the complete cg distribution,
9711 * store the count so we can possibly skip the cg info communication.
9713 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9716 if (comm->DD_debug > 0)
9718 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9719 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9720 "after partitioning");