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)
2743 /* The cut-off might have changed, e.g. by PME load balacning,
2744 * from the value used to set comm->cellsize_min, so check it.
2746 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2748 if (comm->bPMELoadBalDLBLimits)
2750 /* Check for the cut-off limit set by the PME load balancing */
2751 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2755 return cellsize_min;
2758 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2761 real grid_jump_limit;
2763 /* The distance between the boundaries of cells at distance
2764 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2765 * and by the fact that cells should not be shifted by more than
2766 * half their size, such that cg's only shift by one cell
2767 * at redecomposition.
2769 grid_jump_limit = comm->cellsize_limit;
2770 if (!comm->bVacDLBNoLimit)
2772 if (comm->bPMELoadBalDLBLimits)
2774 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2776 grid_jump_limit = max(grid_jump_limit,
2777 cutoff/comm->cd[dim_ind].np);
2780 return grid_jump_limit;
2783 static gmx_bool check_grid_jump(gmx_large_int_t step,
2789 gmx_domdec_comm_t *comm;
2798 for (d = 1; d < dd->ndim; d++)
2801 limit = grid_jump_limit(comm, cutoff, d);
2802 bfac = ddbox->box_size[dim];
2803 if (ddbox->tric_dir[dim])
2805 bfac *= ddbox->skew_fac[dim];
2807 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2808 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2816 /* This error should never be triggered under normal
2817 * circumstances, but you never know ...
2819 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.",
2820 gmx_step_str(step, buf),
2821 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2829 static int dd_load_count(gmx_domdec_comm_t *comm)
2831 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2834 static float dd_force_load(gmx_domdec_comm_t *comm)
2841 if (comm->eFlop > 1)
2843 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2848 load = comm->cycl[ddCyclF];
2849 if (comm->cycl_n[ddCyclF] > 1)
2851 /* Subtract the maximum of the last n cycle counts
2852 * to get rid of possible high counts due to other soures,
2853 * for instance system activity, that would otherwise
2854 * affect the dynamic load balancing.
2856 load -= comm->cycl_max[ddCyclF];
2863 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2865 gmx_domdec_comm_t *comm;
2870 snew(*dim_f, dd->nc[dim]+1);
2872 for (i = 1; i < dd->nc[dim]; i++)
2874 if (comm->slb_frac[dim])
2876 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2880 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2883 (*dim_f)[dd->nc[dim]] = 1;
2886 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2888 int pmeindex, slab, nso, i;
2891 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2897 ddpme->dim = dimind;
2899 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2901 ddpme->nslab = (ddpme->dim == 0 ?
2902 dd->comm->npmenodes_x :
2903 dd->comm->npmenodes_y);
2905 if (ddpme->nslab <= 1)
2910 nso = dd->comm->npmenodes/ddpme->nslab;
2911 /* Determine for each PME slab the PP location range for dimension dim */
2912 snew(ddpme->pp_min, ddpme->nslab);
2913 snew(ddpme->pp_max, ddpme->nslab);
2914 for (slab = 0; slab < ddpme->nslab; slab++)
2916 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2917 ddpme->pp_max[slab] = 0;
2919 for (i = 0; i < dd->nnodes; i++)
2921 ddindex2xyz(dd->nc, i, xyz);
2922 /* For y only use our y/z slab.
2923 * This assumes that the PME x grid size matches the DD grid size.
2925 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2927 pmeindex = ddindex2pmeindex(dd, i);
2930 slab = pmeindex/nso;
2934 slab = pmeindex % ddpme->nslab;
2936 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2937 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2941 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2944 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2946 if (dd->comm->ddpme[0].dim == XX)
2948 return dd->comm->ddpme[0].maxshift;
2956 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2958 if (dd->comm->ddpme[0].dim == YY)
2960 return dd->comm->ddpme[0].maxshift;
2962 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2964 return dd->comm->ddpme[1].maxshift;
2972 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2973 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2975 gmx_domdec_comm_t *comm;
2978 real range, pme_boundary;
2982 nc = dd->nc[ddpme->dim];
2985 if (!ddpme->dim_match)
2987 /* PP decomposition is not along dim: the worst situation */
2990 else if (ns <= 3 || (bUniform && ns == nc))
2992 /* The optimal situation */
2997 /* We need to check for all pme nodes which nodes they
2998 * could possibly need to communicate with.
3000 xmin = ddpme->pp_min;
3001 xmax = ddpme->pp_max;
3002 /* Allow for atoms to be maximally 2/3 times the cut-off
3003 * out of their DD cell. This is a reasonable balance between
3004 * between performance and support for most charge-group/cut-off
3007 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3008 /* Avoid extra communication when we are exactly at a boundary */
3012 for (s = 0; s < ns; s++)
3014 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3015 pme_boundary = (real)s/ns;
3018 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3020 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3024 pme_boundary = (real)(s+1)/ns;
3027 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3029 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3036 ddpme->maxshift = sh;
3040 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3041 ddpme->dim, ddpme->maxshift);
3045 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3049 for (d = 0; d < dd->ndim; d++)
3052 if (dim < ddbox->nboundeddim &&
3053 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3054 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3056 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",
3057 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3058 dd->nc[dim], dd->comm->cellsize_limit);
3063 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3064 gmx_bool bMaster, ivec npulse)
3066 gmx_domdec_comm_t *comm;
3069 real *cell_x, cell_dx, cellsize;
3073 for (d = 0; d < DIM; d++)
3075 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3077 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3080 cell_dx = ddbox->box_size[d]/dd->nc[d];
3083 for (j = 0; j < dd->nc[d]+1; j++)
3085 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3090 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3091 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3093 cellsize = cell_dx*ddbox->skew_fac[d];
3094 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3098 cellsize_min[d] = cellsize;
3102 /* Statically load balanced grid */
3103 /* Also when we are not doing a master distribution we determine
3104 * all cell borders in a loop to obtain identical values
3105 * to the master distribution case and to determine npulse.
3109 cell_x = dd->ma->cell_x[d];
3113 snew(cell_x, dd->nc[d]+1);
3115 cell_x[0] = ddbox->box0[d];
3116 for (j = 0; j < dd->nc[d]; j++)
3118 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3119 cell_x[j+1] = cell_x[j] + cell_dx;
3120 cellsize = cell_dx*ddbox->skew_fac[d];
3121 while (cellsize*npulse[d] < comm->cutoff &&
3122 npulse[d] < dd->nc[d]-1)
3126 cellsize_min[d] = min(cellsize_min[d], cellsize);
3130 comm->cell_x0[d] = cell_x[dd->ci[d]];
3131 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3135 /* The following limitation is to avoid that a cell would receive
3136 * some of its own home charge groups back over the periodic boundary.
3137 * Double charge groups cause trouble with the global indices.
3139 if (d < ddbox->npbcdim &&
3140 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3142 gmx_fatal_collective(FARGS, NULL, dd,
3143 "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",
3144 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3146 dd->nc[d], dd->nc[d],
3147 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3151 if (!comm->bDynLoadBal)
3153 copy_rvec(cellsize_min, comm->cellsize_min);
3156 for (d = 0; d < comm->npmedecompdim; d++)
3158 set_pme_maxshift(dd, &comm->ddpme[d],
3159 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3160 comm->ddpme[d].slb_dim_f);
3165 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3166 int d, int dim, gmx_domdec_root_t *root,
3168 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3170 gmx_domdec_comm_t *comm;
3171 int ncd, i, j, nmin, nmin_old;
3172 gmx_bool bLimLo, bLimHi;
3174 real fac, halfway, cellsize_limit_f_i, region_size;
3175 gmx_bool bPBC, bLastHi = FALSE;
3176 int nrange[] = {range[0], range[1]};
3178 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3184 bPBC = (dim < ddbox->npbcdim);
3186 cell_size = root->buf_ncd;
3190 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3193 /* First we need to check if the scaling does not make cells
3194 * smaller than the smallest allowed size.
3195 * We need to do this iteratively, since if a cell is too small,
3196 * it needs to be enlarged, which makes all the other cells smaller,
3197 * which could in turn make another cell smaller than allowed.
3199 for (i = range[0]; i < range[1]; i++)
3201 root->bCellMin[i] = FALSE;
3207 /* We need the total for normalization */
3209 for (i = range[0]; i < range[1]; i++)
3211 if (root->bCellMin[i] == FALSE)
3213 fac += cell_size[i];
3216 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3217 /* Determine the cell boundaries */
3218 for (i = range[0]; i < range[1]; i++)
3220 if (root->bCellMin[i] == FALSE)
3222 cell_size[i] *= fac;
3223 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3225 cellsize_limit_f_i = 0;
3229 cellsize_limit_f_i = cellsize_limit_f;
3231 if (cell_size[i] < cellsize_limit_f_i)
3233 root->bCellMin[i] = TRUE;
3234 cell_size[i] = cellsize_limit_f_i;
3238 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3241 while (nmin > nmin_old);
3244 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3245 /* For this check we should not use DD_CELL_MARGIN,
3246 * but a slightly smaller factor,
3247 * since rounding could get use below the limit.
3249 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3252 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",
3253 gmx_step_str(step, buf),
3254 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3255 ncd, comm->cellsize_min[dim]);
3258 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3262 /* Check if the boundary did not displace more than halfway
3263 * each of the cells it bounds, as this could cause problems,
3264 * especially when the differences between cell sizes are large.
3265 * If changes are applied, they will not make cells smaller
3266 * than the cut-off, as we check all the boundaries which
3267 * might be affected by a change and if the old state was ok,
3268 * the cells will at most be shrunk back to their old size.
3270 for (i = range[0]+1; i < range[1]; i++)
3272 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3273 if (root->cell_f[i] < halfway)
3275 root->cell_f[i] = halfway;
3276 /* Check if the change also causes shifts of the next boundaries */
3277 for (j = i+1; j < range[1]; j++)
3279 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3281 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3285 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3286 if (root->cell_f[i] > halfway)
3288 root->cell_f[i] = halfway;
3289 /* Check if the change also causes shifts of the next boundaries */
3290 for (j = i-1; j >= range[0]+1; j--)
3292 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3294 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3301 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3302 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3303 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3304 * for a and b nrange is used */
3307 /* Take care of the staggering of the cell boundaries */
3310 for (i = range[0]; i < range[1]; i++)
3312 root->cell_f_max0[i] = root->cell_f[i];
3313 root->cell_f_min1[i] = root->cell_f[i+1];
3318 for (i = range[0]+1; i < range[1]; i++)
3320 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3321 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3322 if (bLimLo && bLimHi)
3324 /* Both limits violated, try the best we can */
3325 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3326 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3327 nrange[0] = range[0];
3329 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3332 nrange[1] = range[1];
3333 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3339 /* root->cell_f[i] = root->bound_min[i]; */
3340 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3343 else if (bLimHi && !bLastHi)
3346 if (nrange[1] < range[1]) /* found a LimLo before */
3348 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3349 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3350 nrange[0] = nrange[1];
3352 root->cell_f[i] = root->bound_max[i];
3354 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3356 nrange[1] = range[1];
3359 if (nrange[1] < range[1]) /* found last a LimLo */
3361 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3362 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3363 nrange[0] = nrange[1];
3364 nrange[1] = range[1];
3365 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3367 else if (nrange[0] > range[0]) /* found at least one LimHi */
3369 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3376 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3377 int d, int dim, gmx_domdec_root_t *root,
3378 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3379 gmx_bool bUniform, gmx_large_int_t step)
3381 gmx_domdec_comm_t *comm;
3382 int ncd, d1, i, j, pos;
3384 real load_aver, load_i, imbalance, change, change_max, sc;
3385 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3389 int range[] = { 0, 0 };
3393 /* Convert the maximum change from the input percentage to a fraction */
3394 change_limit = comm->dlb_scale_lim*0.01;
3398 bPBC = (dim < ddbox->npbcdim);
3400 cell_size = root->buf_ncd;
3402 /* Store the original boundaries */
3403 for (i = 0; i < ncd+1; i++)
3405 root->old_cell_f[i] = root->cell_f[i];
3409 for (i = 0; i < ncd; i++)
3411 cell_size[i] = 1.0/ncd;
3414 else if (dd_load_count(comm))
3416 load_aver = comm->load[d].sum_m/ncd;
3418 for (i = 0; i < ncd; i++)
3420 /* Determine the relative imbalance of cell i */
3421 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3422 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3423 /* Determine the change of the cell size using underrelaxation */
3424 change = -relax*imbalance;
3425 change_max = max(change_max, max(change, -change));
3427 /* Limit the amount of scaling.
3428 * We need to use the same rescaling for all cells in one row,
3429 * otherwise the load balancing might not converge.
3432 if (change_max > change_limit)
3434 sc *= change_limit/change_max;
3436 for (i = 0; i < ncd; i++)
3438 /* Determine the relative imbalance of cell i */
3439 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3440 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3441 /* Determine the change of the cell size using underrelaxation */
3442 change = -sc*imbalance;
3443 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3447 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3448 cellsize_limit_f *= DD_CELL_MARGIN;
3449 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3450 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3451 if (ddbox->tric_dir[dim])
3453 cellsize_limit_f /= ddbox->skew_fac[dim];
3454 dist_min_f /= ddbox->skew_fac[dim];
3456 if (bDynamicBox && d > 0)
3458 dist_min_f *= DD_PRES_SCALE_MARGIN;
3460 if (d > 0 && !bUniform)
3462 /* Make sure that the grid is not shifted too much */
3463 for (i = 1; i < ncd; i++)
3465 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3467 gmx_incons("Inconsistent DD boundary staggering limits!");
3469 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3470 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3473 root->bound_min[i] += 0.5*space;
3475 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3476 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3479 root->bound_max[i] += 0.5*space;
3484 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3486 root->cell_f_max0[i-1] + dist_min_f,
3487 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3488 root->cell_f_min1[i] - dist_min_f);
3493 root->cell_f[0] = 0;
3494 root->cell_f[ncd] = 1;
3495 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3498 /* After the checks above, the cells should obey the cut-off
3499 * restrictions, but it does not hurt to check.
3501 for (i = 0; i < ncd; i++)
3505 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3506 dim, i, root->cell_f[i], root->cell_f[i+1]);
3509 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3510 root->cell_f[i+1] - root->cell_f[i] <
3511 cellsize_limit_f/DD_CELL_MARGIN)
3515 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3516 gmx_step_str(step, buf), dim2char(dim), i,
3517 (root->cell_f[i+1] - root->cell_f[i])
3518 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3523 /* Store the cell boundaries of the lower dimensions at the end */
3524 for (d1 = 0; d1 < d; d1++)
3526 root->cell_f[pos++] = comm->cell_f0[d1];
3527 root->cell_f[pos++] = comm->cell_f1[d1];
3530 if (d < comm->npmedecompdim)
3532 /* The master determines the maximum shift for
3533 * the coordinate communication between separate PME nodes.
3535 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3537 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3540 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3544 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3545 gmx_ddbox_t *ddbox, int dimind)
3547 gmx_domdec_comm_t *comm;
3552 /* Set the cell dimensions */
3553 dim = dd->dim[dimind];
3554 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3555 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3556 if (dim >= ddbox->nboundeddim)
3558 comm->cell_x0[dim] += ddbox->box0[dim];
3559 comm->cell_x1[dim] += ddbox->box0[dim];
3563 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3564 int d, int dim, real *cell_f_row,
3567 gmx_domdec_comm_t *comm;
3573 /* Each node would only need to know two fractions,
3574 * but it is probably cheaper to broadcast the whole array.
3576 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3577 0, comm->mpi_comm_load[d]);
3579 /* Copy the fractions for this dimension from the buffer */
3580 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3581 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3582 /* The whole array was communicated, so set the buffer position */
3583 pos = dd->nc[dim] + 1;
3584 for (d1 = 0; d1 <= d; d1++)
3588 /* Copy the cell fractions of the lower dimensions */
3589 comm->cell_f0[d1] = cell_f_row[pos++];
3590 comm->cell_f1[d1] = cell_f_row[pos++];
3592 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3594 /* Convert the communicated shift from float to int */
3595 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3598 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3602 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3603 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3604 gmx_bool bUniform, gmx_large_int_t step)
3606 gmx_domdec_comm_t *comm;
3608 gmx_bool bRowMember, bRowRoot;
3613 for (d = 0; d < dd->ndim; d++)
3618 for (d1 = d; d1 < dd->ndim; d1++)
3620 if (dd->ci[dd->dim[d1]] > 0)
3633 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3634 ddbox, bDynamicBox, bUniform, step);
3635 cell_f_row = comm->root[d]->cell_f;
3639 cell_f_row = comm->cell_f_row;
3641 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3646 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3650 /* This function assumes the box is static and should therefore
3651 * not be called when the box has changed since the last
3652 * call to dd_partition_system.
3654 for (d = 0; d < dd->ndim; d++)
3656 relative_to_absolute_cell_bounds(dd, ddbox, d);
3662 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3663 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3664 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3665 gmx_wallcycle_t wcycle)
3667 gmx_domdec_comm_t *comm;
3674 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3675 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3676 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3678 else if (bDynamicBox)
3680 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3683 /* Set the dimensions for which no DD is used */
3684 for (dim = 0; dim < DIM; dim++)
3686 if (dd->nc[dim] == 1)
3688 comm->cell_x0[dim] = 0;
3689 comm->cell_x1[dim] = ddbox->box_size[dim];
3690 if (dim >= ddbox->nboundeddim)
3692 comm->cell_x0[dim] += ddbox->box0[dim];
3693 comm->cell_x1[dim] += ddbox->box0[dim];
3699 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3702 gmx_domdec_comm_dim_t *cd;
3704 for (d = 0; d < dd->ndim; d++)
3706 cd = &dd->comm->cd[d];
3707 np = npulse[dd->dim[d]];
3708 if (np > cd->np_nalloc)
3712 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3713 dim2char(dd->dim[d]), np);
3715 if (DDMASTER(dd) && cd->np_nalloc > 0)
3717 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3719 srenew(cd->ind, np);
3720 for (i = cd->np_nalloc; i < np; i++)
3722 cd->ind[i].index = NULL;
3723 cd->ind[i].nalloc = 0;
3732 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3733 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3734 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3735 gmx_wallcycle_t wcycle)
3737 gmx_domdec_comm_t *comm;
3743 /* Copy the old cell boundaries for the cg displacement check */
3744 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3745 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3747 if (comm->bDynLoadBal)
3751 check_box_size(dd, ddbox);
3753 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3757 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3758 realloc_comm_ind(dd, npulse);
3763 for (d = 0; d < DIM; d++)
3765 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3766 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3771 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3773 rvec cell_ns_x0, rvec cell_ns_x1,
3774 gmx_large_int_t step)
3776 gmx_domdec_comm_t *comm;
3781 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3783 dim = dd->dim[dim_ind];
3785 /* Without PBC we don't have restrictions on the outer cells */
3786 if (!(dim >= ddbox->npbcdim &&
3787 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3788 comm->bDynLoadBal &&
3789 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3790 comm->cellsize_min[dim])
3793 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",
3794 gmx_step_str(step, buf), dim2char(dim),
3795 comm->cell_x1[dim] - comm->cell_x0[dim],
3796 ddbox->skew_fac[dim],
3797 dd->comm->cellsize_min[dim],
3798 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3802 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3804 /* Communicate the boundaries and update cell_ns_x0/1 */
3805 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3806 if (dd->bGridJump && dd->ndim > 1)
3808 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3813 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3817 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3825 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3826 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3835 static void check_screw_box(matrix box)
3837 /* Mathematical limitation */
3838 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3840 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3843 /* Limitation due to the asymmetry of the eighth shell method */
3844 if (box[ZZ][YY] != 0)
3846 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3850 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3851 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3854 gmx_domdec_master_t *ma;
3855 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3856 int i, icg, j, k, k0, k1, d, npbcdim;
3858 rvec box_size, cg_cm;
3860 real nrcg, inv_ncg, pos_d;
3862 gmx_bool bUnbounded, bScrew;
3866 if (tmp_ind == NULL)
3868 snew(tmp_nalloc, dd->nnodes);
3869 snew(tmp_ind, dd->nnodes);
3870 for (i = 0; i < dd->nnodes; i++)
3872 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3873 snew(tmp_ind[i], tmp_nalloc[i]);
3877 /* Clear the count */
3878 for (i = 0; i < dd->nnodes; i++)
3884 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3886 cgindex = cgs->index;
3888 /* Compute the center of geometry for all charge groups */
3889 for (icg = 0; icg < cgs->nr; icg++)
3892 k1 = cgindex[icg+1];
3896 copy_rvec(pos[k0], cg_cm);
3903 for (k = k0; (k < k1); k++)
3905 rvec_inc(cg_cm, pos[k]);
3907 for (d = 0; (d < DIM); d++)
3909 cg_cm[d] *= inv_ncg;
3912 /* Put the charge group in the box and determine the cell index */
3913 for (d = DIM-1; d >= 0; d--)
3916 if (d < dd->npbcdim)
3918 bScrew = (dd->bScrewPBC && d == XX);
3919 if (tric_dir[d] && dd->nc[d] > 1)
3921 /* Use triclinic coordintates for this dimension */
3922 for (j = d+1; j < DIM; j++)
3924 pos_d += cg_cm[j]*tcm[j][d];
3927 while (pos_d >= box[d][d])
3930 rvec_dec(cg_cm, box[d]);
3933 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3934 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3936 for (k = k0; (k < k1); k++)
3938 rvec_dec(pos[k], box[d]);
3941 pos[k][YY] = box[YY][YY] - pos[k][YY];
3942 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3949 rvec_inc(cg_cm, box[d]);
3952 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3953 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3955 for (k = k0; (k < k1); k++)
3957 rvec_inc(pos[k], box[d]);
3960 pos[k][YY] = box[YY][YY] - pos[k][YY];
3961 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3966 /* This could be done more efficiently */
3968 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3973 i = dd_index(dd->nc, ind);
3974 if (ma->ncg[i] == tmp_nalloc[i])
3976 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3977 srenew(tmp_ind[i], tmp_nalloc[i]);
3979 tmp_ind[i][ma->ncg[i]] = icg;
3981 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3985 for (i = 0; i < dd->nnodes; i++)
3988 for (k = 0; k < ma->ncg[i]; k++)
3990 ma->cg[k1++] = tmp_ind[i][k];
3993 ma->index[dd->nnodes] = k1;
3995 for (i = 0; i < dd->nnodes; i++)
4005 fprintf(fplog, "Charge group distribution at step %s:",
4006 gmx_step_str(step, buf));
4007 for (i = 0; i < dd->nnodes; i++)
4009 fprintf(fplog, " %d", ma->ncg[i]);
4011 fprintf(fplog, "\n");
4015 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4016 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4019 gmx_domdec_master_t *ma = NULL;
4022 int *ibuf, buf2[2] = { 0, 0 };
4023 gmx_bool bMaster = DDMASTER(dd);
4030 check_screw_box(box);
4033 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4035 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4036 for (i = 0; i < dd->nnodes; i++)
4038 ma->ibuf[2*i] = ma->ncg[i];
4039 ma->ibuf[2*i+1] = ma->nat[i];
4047 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4049 dd->ncg_home = buf2[0];
4050 dd->nat_home = buf2[1];
4051 dd->ncg_tot = dd->ncg_home;
4052 dd->nat_tot = dd->nat_home;
4053 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4055 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4056 srenew(dd->index_gl, dd->cg_nalloc);
4057 srenew(dd->cgindex, dd->cg_nalloc+1);
4061 for (i = 0; i < dd->nnodes; i++)
4063 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4064 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4069 DDMASTER(dd) ? ma->ibuf : NULL,
4070 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4071 DDMASTER(dd) ? ma->cg : NULL,
4072 dd->ncg_home*sizeof(int), dd->index_gl);
4074 /* Determine the home charge group sizes */
4076 for (i = 0; i < dd->ncg_home; i++)
4078 cg_gl = dd->index_gl[i];
4080 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4085 fprintf(debug, "Home charge groups:\n");
4086 for (i = 0; i < dd->ncg_home; i++)
4088 fprintf(debug, " %d", dd->index_gl[i]);
4091 fprintf(debug, "\n");
4094 fprintf(debug, "\n");
4098 static int compact_and_copy_vec_at(int ncg, int *move,
4101 rvec *src, gmx_domdec_comm_t *comm,
4104 int m, icg, i, i0, i1, nrcg;
4110 for (m = 0; m < DIM*2; m++)
4116 for (icg = 0; icg < ncg; icg++)
4118 i1 = cgindex[icg+1];
4124 /* Compact the home array in place */
4125 for (i = i0; i < i1; i++)
4127 copy_rvec(src[i], src[home_pos++]);
4133 /* Copy to the communication buffer */
4135 pos_vec[m] += 1 + vec*nrcg;
4136 for (i = i0; i < i1; i++)
4138 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4140 pos_vec[m] += (nvec - vec - 1)*nrcg;
4144 home_pos += i1 - i0;
4152 static int compact_and_copy_vec_cg(int ncg, int *move,
4154 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4157 int m, icg, i0, i1, nrcg;
4163 for (m = 0; m < DIM*2; m++)
4169 for (icg = 0; icg < ncg; icg++)
4171 i1 = cgindex[icg+1];
4177 /* Compact the home array in place */
4178 copy_rvec(src[icg], src[home_pos++]);
4184 /* Copy to the communication buffer */
4185 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4186 pos_vec[m] += 1 + nrcg*nvec;
4198 static int compact_ind(int ncg, int *move,
4199 int *index_gl, int *cgindex,
4201 gmx_ga2la_t ga2la, char *bLocalCG,
4204 int cg, nat, a0, a1, a, a_gl;
4209 for (cg = 0; cg < ncg; cg++)
4215 /* Compact the home arrays in place.
4216 * Anything that can be done here avoids access to global arrays.
4218 cgindex[home_pos] = nat;
4219 for (a = a0; a < a1; a++)
4222 gatindex[nat] = a_gl;
4223 /* The cell number stays 0, so we don't need to set it */
4224 ga2la_change_la(ga2la, a_gl, nat);
4227 index_gl[home_pos] = index_gl[cg];
4228 cginfo[home_pos] = cginfo[cg];
4229 /* The charge group remains local, so bLocalCG does not change */
4234 /* Clear the global indices */
4235 for (a = a0; a < a1; a++)
4237 ga2la_del(ga2la, gatindex[a]);
4241 bLocalCG[index_gl[cg]] = FALSE;
4245 cgindex[home_pos] = nat;
4250 static void clear_and_mark_ind(int ncg, int *move,
4251 int *index_gl, int *cgindex, int *gatindex,
4252 gmx_ga2la_t ga2la, char *bLocalCG,
4257 for (cg = 0; cg < ncg; cg++)
4263 /* Clear the global indices */
4264 for (a = a0; a < a1; a++)
4266 ga2la_del(ga2la, gatindex[a]);
4270 bLocalCG[index_gl[cg]] = FALSE;
4272 /* Signal that this cg has moved using the ns cell index.
4273 * Here we set it to -1. fill_grid will change it
4274 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4276 cell_index[cg] = -1;
4281 static void print_cg_move(FILE *fplog,
4283 gmx_large_int_t step, int cg, int dim, int dir,
4284 gmx_bool bHaveLimitdAndCMOld, real limitd,
4285 rvec cm_old, rvec cm_new, real pos_d)
4287 gmx_domdec_comm_t *comm;
4292 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4293 if (bHaveLimitdAndCMOld)
4295 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4296 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4300 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4301 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4303 fprintf(fplog, "distance out of cell %f\n",
4304 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4305 if (bHaveLimitdAndCMOld)
4307 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4308 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4310 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4311 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4312 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4314 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4315 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4317 comm->cell_x0[dim], comm->cell_x1[dim]);
4320 static void cg_move_error(FILE *fplog,
4322 gmx_large_int_t step, int cg, int dim, int dir,
4323 gmx_bool bHaveLimitdAndCMOld, real limitd,
4324 rvec cm_old, rvec cm_new, real pos_d)
4328 print_cg_move(fplog, dd, step, cg, dim, dir,
4329 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4331 print_cg_move(stderr, dd, step, cg, dim, dir,
4332 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4334 "A charge group moved too far between two domain decomposition steps\n"
4335 "This usually means that your system is not well equilibrated");
4338 static void rotate_state_atom(t_state *state, int a)
4342 for (est = 0; est < estNR; est++)
4344 if (EST_DISTR(est) && (state->flags & (1<<est)))
4349 /* Rotate the complete state; for a rectangular box only */
4350 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4351 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4354 state->v[a][YY] = -state->v[a][YY];
4355 state->v[a][ZZ] = -state->v[a][ZZ];
4358 state->sd_X[a][YY] = -state->sd_X[a][YY];
4359 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4362 state->cg_p[a][YY] = -state->cg_p[a][YY];
4363 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4365 case estDISRE_INITF:
4366 case estDISRE_RM3TAV:
4367 case estORIRE_INITF:
4369 /* These are distances, so not affected by rotation */
4372 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4378 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4380 if (natoms > comm->moved_nalloc)
4382 /* Contents should be preserved here */
4383 comm->moved_nalloc = over_alloc_dd(natoms);
4384 srenew(comm->moved, comm->moved_nalloc);
4390 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4393 ivec tric_dir, matrix tcm,
4394 rvec cell_x0, rvec cell_x1,
4395 rvec limitd, rvec limit0, rvec limit1,
4397 int cg_start, int cg_end,
4402 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4403 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4407 real inv_ncg, pos_d;
4410 npbcdim = dd->npbcdim;
4412 for (cg = cg_start; cg < cg_end; cg++)
4419 copy_rvec(state->x[k0], cm_new);
4426 for (k = k0; (k < k1); k++)
4428 rvec_inc(cm_new, state->x[k]);
4430 for (d = 0; (d < DIM); d++)
4432 cm_new[d] = inv_ncg*cm_new[d];
4437 /* Do pbc and check DD cell boundary crossings */
4438 for (d = DIM-1; d >= 0; d--)
4442 bScrew = (dd->bScrewPBC && d == XX);
4443 /* Determine the location of this cg in lattice coordinates */
4447 for (d2 = d+1; d2 < DIM; d2++)
4449 pos_d += cm_new[d2]*tcm[d2][d];
4452 /* Put the charge group in the triclinic unit-cell */
4453 if (pos_d >= cell_x1[d])
4455 if (pos_d >= limit1[d])
4457 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4458 cg_cm[cg], cm_new, pos_d);
4461 if (dd->ci[d] == dd->nc[d] - 1)
4463 rvec_dec(cm_new, state->box[d]);
4466 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4467 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4469 for (k = k0; (k < k1); k++)
4471 rvec_dec(state->x[k], state->box[d]);
4474 rotate_state_atom(state, k);
4479 else if (pos_d < cell_x0[d])
4481 if (pos_d < limit0[d])
4483 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4484 cg_cm[cg], cm_new, pos_d);
4489 rvec_inc(cm_new, state->box[d]);
4492 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4493 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4495 for (k = k0; (k < k1); k++)
4497 rvec_inc(state->x[k], state->box[d]);
4500 rotate_state_atom(state, k);
4506 else if (d < npbcdim)
4508 /* Put the charge group in the rectangular unit-cell */
4509 while (cm_new[d] >= state->box[d][d])
4511 rvec_dec(cm_new, state->box[d]);
4512 for (k = k0; (k < k1); k++)
4514 rvec_dec(state->x[k], state->box[d]);
4517 while (cm_new[d] < 0)
4519 rvec_inc(cm_new, state->box[d]);
4520 for (k = k0; (k < k1); k++)
4522 rvec_inc(state->x[k], state->box[d]);
4528 copy_rvec(cm_new, cg_cm[cg]);
4530 /* Determine where this cg should go */
4533 for (d = 0; d < dd->ndim; d++)
4538 flag |= DD_FLAG_FW(d);
4544 else if (dev[dim] == -1)
4546 flag |= DD_FLAG_BW(d);
4549 if (dd->nc[dim] > 2)
4560 /* Temporarily store the flag in move */
4561 move[cg] = mc + flag;
4565 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4566 gmx_domdec_t *dd, ivec tric_dir,
4567 t_state *state, rvec **f,
4568 t_forcerec *fr, t_mdatoms *md,
4576 int ncg[DIM*2], nat[DIM*2];
4577 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4578 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4579 int sbuf[2], rbuf[2];
4580 int home_pos_cg, home_pos_at, buf_pos;
4582 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4585 real inv_ncg, pos_d;
4587 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4589 cginfo_mb_t *cginfo_mb;
4590 gmx_domdec_comm_t *comm;
4592 int nthread, thread;
4596 check_screw_box(state->box);
4600 if (fr->cutoff_scheme == ecutsGROUP)
4605 for (i = 0; i < estNR; i++)
4611 case estX: /* Always present */ break;
4612 case estV: bV = (state->flags & (1<<i)); break;
4613 case estSDX: bSDX = (state->flags & (1<<i)); break;
4614 case estCGP: bCGP = (state->flags & (1<<i)); break;
4617 case estDISRE_INITF:
4618 case estDISRE_RM3TAV:
4619 case estORIRE_INITF:
4621 /* No processing required */
4624 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4629 if (dd->ncg_tot > comm->nalloc_int)
4631 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4632 srenew(comm->buf_int, comm->nalloc_int);
4634 move = comm->buf_int;
4636 /* Clear the count */
4637 for (c = 0; c < dd->ndim*2; c++)
4643 npbcdim = dd->npbcdim;
4645 for (d = 0; (d < DIM); d++)
4647 limitd[d] = dd->comm->cellsize_min[d];
4648 if (d >= npbcdim && dd->ci[d] == 0)
4650 cell_x0[d] = -GMX_FLOAT_MAX;
4654 cell_x0[d] = comm->cell_x0[d];
4656 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4658 cell_x1[d] = GMX_FLOAT_MAX;
4662 cell_x1[d] = comm->cell_x1[d];
4666 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4667 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4671 /* We check after communication if a charge group moved
4672 * more than one cell. Set the pre-comm check limit to float_max.
4674 limit0[d] = -GMX_FLOAT_MAX;
4675 limit1[d] = GMX_FLOAT_MAX;
4679 make_tric_corr_matrix(npbcdim, state->box, tcm);
4681 cgindex = dd->cgindex;
4683 nthread = gmx_omp_nthreads_get(emntDomdec);
4685 /* Compute the center of geometry for all home charge groups
4686 * and put them in the box and determine where they should go.
4688 #pragma omp parallel for num_threads(nthread) schedule(static)
4689 for (thread = 0; thread < nthread; thread++)
4691 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4692 cell_x0, cell_x1, limitd, limit0, limit1,
4694 ( thread *dd->ncg_home)/nthread,
4695 ((thread+1)*dd->ncg_home)/nthread,
4696 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4700 for (cg = 0; cg < dd->ncg_home; cg++)
4705 flag = mc & ~DD_FLAG_NRCG;
4706 mc = mc & DD_FLAG_NRCG;
4709 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4711 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4712 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4714 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4715 /* We store the cg size in the lower 16 bits
4716 * and the place where the charge group should go
4717 * in the next 6 bits. This saves some communication volume.
4719 nrcg = cgindex[cg+1] - cgindex[cg];
4720 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4726 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4727 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4730 for (i = 0; i < dd->ndim*2; i++)
4732 *ncg_moved += ncg[i];
4749 /* Make sure the communication buffers are large enough */
4750 for (mc = 0; mc < dd->ndim*2; mc++)
4752 nvr = ncg[mc] + nat[mc]*nvec;
4753 if (nvr > comm->cgcm_state_nalloc[mc])
4755 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4756 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4760 switch (fr->cutoff_scheme)
4763 /* Recalculating cg_cm might be cheaper than communicating,
4764 * but that could give rise to rounding issues.
4767 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4768 nvec, cg_cm, comm, bCompact);
4771 /* Without charge groups we send the moved atom coordinates
4772 * over twice. This is so the code below can be used without
4773 * many conditionals for both for with and without charge groups.
4776 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4777 nvec, state->x, comm, FALSE);
4780 home_pos_cg -= *ncg_moved;
4784 gmx_incons("unimplemented");
4790 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4791 nvec, vec++, state->x, comm, bCompact);
4794 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4795 nvec, vec++, state->v, comm, bCompact);
4799 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4800 nvec, vec++, state->sd_X, comm, bCompact);
4804 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4805 nvec, vec++, state->cg_p, comm, bCompact);
4810 compact_ind(dd->ncg_home, move,
4811 dd->index_gl, dd->cgindex, dd->gatindex,
4812 dd->ga2la, comm->bLocalCG,
4817 if (fr->cutoff_scheme == ecutsVERLET)
4819 moved = get_moved(comm, dd->ncg_home);
4821 for (k = 0; k < dd->ncg_home; k++)
4828 moved = fr->ns.grid->cell_index;
4831 clear_and_mark_ind(dd->ncg_home, move,
4832 dd->index_gl, dd->cgindex, dd->gatindex,
4833 dd->ga2la, comm->bLocalCG,
4837 cginfo_mb = fr->cginfo_mb;
4839 *ncg_stay_home = home_pos_cg;
4840 for (d = 0; d < dd->ndim; d++)
4846 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4849 /* Communicate the cg and atom counts */
4854 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4855 d, dir, sbuf[0], sbuf[1]);
4857 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4859 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4861 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4862 srenew(comm->buf_int, comm->nalloc_int);
4865 /* Communicate the charge group indices, sizes and flags */
4866 dd_sendrecv_int(dd, d, dir,
4867 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4868 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4870 nvs = ncg[cdd] + nat[cdd]*nvec;
4871 i = rbuf[0] + rbuf[1] *nvec;
4872 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4874 /* Communicate cgcm and state */
4875 dd_sendrecv_rvec(dd, d, dir,
4876 comm->cgcm_state[cdd], nvs,
4877 comm->vbuf.v+nvr, i);
4878 ncg_recv += rbuf[0];
4879 nat_recv += rbuf[1];
4883 /* Process the received charge groups */
4885 for (cg = 0; cg < ncg_recv; cg++)
4887 flag = comm->buf_int[cg*DD_CGIBS+1];
4889 if (dim >= npbcdim && dd->nc[dim] > 2)
4891 /* No pbc in this dim and more than one domain boundary.
4892 * We do a separate check if a charge group didn't move too far.
4894 if (((flag & DD_FLAG_FW(d)) &&
4895 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4896 ((flag & DD_FLAG_BW(d)) &&
4897 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4899 cg_move_error(fplog, dd, step, cg, dim,
4900 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4902 comm->vbuf.v[buf_pos],
4903 comm->vbuf.v[buf_pos],
4904 comm->vbuf.v[buf_pos][dim]);
4911 /* Check which direction this cg should go */
4912 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4916 /* The cell boundaries for dimension d2 are not equal
4917 * for each cell row of the lower dimension(s),
4918 * therefore we might need to redetermine where
4919 * this cg should go.
4922 /* If this cg crosses the box boundary in dimension d2
4923 * we can use the communicated flag, so we do not
4924 * have to worry about pbc.
4926 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4927 (flag & DD_FLAG_FW(d2))) ||
4928 (dd->ci[dim2] == 0 &&
4929 (flag & DD_FLAG_BW(d2)))))
4931 /* Clear the two flags for this dimension */
4932 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4933 /* Determine the location of this cg
4934 * in lattice coordinates
4936 pos_d = comm->vbuf.v[buf_pos][dim2];
4939 for (d3 = dim2+1; d3 < DIM; d3++)
4942 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4945 /* Check of we are not at the box edge.
4946 * pbc is only handled in the first step above,
4947 * but this check could move over pbc while
4948 * the first step did not due to different rounding.
4950 if (pos_d >= cell_x1[dim2] &&
4951 dd->ci[dim2] != dd->nc[dim2]-1)
4953 flag |= DD_FLAG_FW(d2);
4955 else if (pos_d < cell_x0[dim2] &&
4958 flag |= DD_FLAG_BW(d2);
4960 comm->buf_int[cg*DD_CGIBS+1] = flag;
4963 /* Set to which neighboring cell this cg should go */
4964 if (flag & DD_FLAG_FW(d2))
4968 else if (flag & DD_FLAG_BW(d2))
4970 if (dd->nc[dd->dim[d2]] > 2)
4982 nrcg = flag & DD_FLAG_NRCG;
4985 if (home_pos_cg+1 > dd->cg_nalloc)
4987 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4988 srenew(dd->index_gl, dd->cg_nalloc);
4989 srenew(dd->cgindex, dd->cg_nalloc+1);
4991 /* Set the global charge group index and size */
4992 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4993 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4994 /* Copy the state from the buffer */
4995 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
4996 if (fr->cutoff_scheme == ecutsGROUP)
4999 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5003 /* Set the cginfo */
5004 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5005 dd->index_gl[home_pos_cg]);
5008 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5011 if (home_pos_at+nrcg > state->nalloc)
5013 dd_realloc_state(state, f, home_pos_at+nrcg);
5015 for (i = 0; i < nrcg; i++)
5017 copy_rvec(comm->vbuf.v[buf_pos++],
5018 state->x[home_pos_at+i]);
5022 for (i = 0; i < nrcg; i++)
5024 copy_rvec(comm->vbuf.v[buf_pos++],
5025 state->v[home_pos_at+i]);
5030 for (i = 0; i < nrcg; i++)
5032 copy_rvec(comm->vbuf.v[buf_pos++],
5033 state->sd_X[home_pos_at+i]);
5038 for (i = 0; i < nrcg; i++)
5040 copy_rvec(comm->vbuf.v[buf_pos++],
5041 state->cg_p[home_pos_at+i]);
5045 home_pos_at += nrcg;
5049 /* Reallocate the buffers if necessary */
5050 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5052 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5053 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5055 nvr = ncg[mc] + nat[mc]*nvec;
5056 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5058 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5059 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5061 /* Copy from the receive to the send buffers */
5062 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5063 comm->buf_int + cg*DD_CGIBS,
5064 DD_CGIBS*sizeof(int));
5065 memcpy(comm->cgcm_state[mc][nvr],
5066 comm->vbuf.v[buf_pos],
5067 (1+nrcg*nvec)*sizeof(rvec));
5068 buf_pos += 1 + nrcg*nvec;
5075 /* With sorting (!bCompact) the indices are now only partially up to date
5076 * and ncg_home and nat_home are not the real count, since there are
5077 * "holes" in the arrays for the charge groups that moved to neighbors.
5079 if (fr->cutoff_scheme == ecutsVERLET)
5081 moved = get_moved(comm, home_pos_cg);
5083 for (i = dd->ncg_home; i < home_pos_cg; i++)
5088 dd->ncg_home = home_pos_cg;
5089 dd->nat_home = home_pos_at;
5094 "Finished repartitioning: cgs moved out %d, new home %d\n",
5095 *ncg_moved, dd->ncg_home-*ncg_moved);
5100 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5102 dd->comm->cycl[ddCycl] += cycles;
5103 dd->comm->cycl_n[ddCycl]++;
5104 if (cycles > dd->comm->cycl_max[ddCycl])
5106 dd->comm->cycl_max[ddCycl] = cycles;
5110 static double force_flop_count(t_nrnb *nrnb)
5117 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5119 /* To get closer to the real timings, we half the count
5120 * for the normal loops and again half it for water loops.
5123 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5125 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5129 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5132 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5135 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5137 sum += nrnb->n[i]*cost_nrnb(i);
5140 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5142 sum += nrnb->n[i]*cost_nrnb(i);
5148 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5150 if (dd->comm->eFlop)
5152 dd->comm->flop -= force_flop_count(nrnb);
5155 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5157 if (dd->comm->eFlop)
5159 dd->comm->flop += force_flop_count(nrnb);
5164 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5168 for (i = 0; i < ddCyclNr; i++)
5170 dd->comm->cycl[i] = 0;
5171 dd->comm->cycl_n[i] = 0;
5172 dd->comm->cycl_max[i] = 0;
5175 dd->comm->flop_n = 0;
5178 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5180 gmx_domdec_comm_t *comm;
5181 gmx_domdec_load_t *load;
5182 gmx_domdec_root_t *root = NULL;
5183 int d, dim, cid, i, pos;
5184 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5189 fprintf(debug, "get_load_distribution start\n");
5192 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5196 bSepPME = (dd->pme_nodeid >= 0);
5198 for (d = dd->ndim-1; d >= 0; d--)
5201 /* Check if we participate in the communication in this dimension */
5202 if (d == dd->ndim-1 ||
5203 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5205 load = &comm->load[d];
5208 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5211 if (d == dd->ndim-1)
5213 sbuf[pos++] = dd_force_load(comm);
5214 sbuf[pos++] = sbuf[0];
5217 sbuf[pos++] = sbuf[0];
5218 sbuf[pos++] = cell_frac;
5221 sbuf[pos++] = comm->cell_f_max0[d];
5222 sbuf[pos++] = comm->cell_f_min1[d];
5227 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5228 sbuf[pos++] = comm->cycl[ddCyclPME];
5233 sbuf[pos++] = comm->load[d+1].sum;
5234 sbuf[pos++] = comm->load[d+1].max;
5237 sbuf[pos++] = comm->load[d+1].sum_m;
5238 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5239 sbuf[pos++] = comm->load[d+1].flags;
5242 sbuf[pos++] = comm->cell_f_max0[d];
5243 sbuf[pos++] = comm->cell_f_min1[d];
5248 sbuf[pos++] = comm->load[d+1].mdf;
5249 sbuf[pos++] = comm->load[d+1].pme;
5253 /* Communicate a row in DD direction d.
5254 * The communicators are setup such that the root always has rank 0.
5257 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5258 load->load, load->nload*sizeof(float), MPI_BYTE,
5259 0, comm->mpi_comm_load[d]);
5261 if (dd->ci[dim] == dd->master_ci[dim])
5263 /* We are the root, process this row */
5264 if (comm->bDynLoadBal)
5266 root = comm->root[d];
5276 for (i = 0; i < dd->nc[dim]; i++)
5278 load->sum += load->load[pos++];
5279 load->max = max(load->max, load->load[pos]);
5285 /* This direction could not be load balanced properly,
5286 * therefore we need to use the maximum iso the average load.
5288 load->sum_m = max(load->sum_m, load->load[pos]);
5292 load->sum_m += load->load[pos];
5295 load->cvol_min = min(load->cvol_min, load->load[pos]);
5299 load->flags = (int)(load->load[pos++] + 0.5);
5303 root->cell_f_max0[i] = load->load[pos++];
5304 root->cell_f_min1[i] = load->load[pos++];
5309 load->mdf = max(load->mdf, load->load[pos]);
5311 load->pme = max(load->pme, load->load[pos]);
5315 if (comm->bDynLoadBal && root->bLimited)
5317 load->sum_m *= dd->nc[dim];
5318 load->flags |= (1<<d);
5326 comm->nload += dd_load_count(comm);
5327 comm->load_step += comm->cycl[ddCyclStep];
5328 comm->load_sum += comm->load[0].sum;
5329 comm->load_max += comm->load[0].max;
5330 if (comm->bDynLoadBal)
5332 for (d = 0; d < dd->ndim; d++)
5334 if (comm->load[0].flags & (1<<d))
5336 comm->load_lim[d]++;
5342 comm->load_mdf += comm->load[0].mdf;
5343 comm->load_pme += comm->load[0].pme;
5347 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5351 fprintf(debug, "get_load_distribution finished\n");
5355 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5357 /* Return the relative performance loss on the total run time
5358 * due to the force calculation load imbalance.
5360 if (dd->comm->nload > 0)
5363 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5364 (dd->comm->load_step*dd->nnodes);
5372 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5375 int npp, npme, nnodes, d, limp;
5376 float imbal, pme_f_ratio, lossf, lossp = 0;
5378 gmx_domdec_comm_t *comm;
5381 if (DDMASTER(dd) && comm->nload > 0)
5384 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5385 nnodes = npp + npme;
5386 imbal = comm->load_max*npp/comm->load_sum - 1;
5387 lossf = dd_force_imb_perf_loss(dd);
5388 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5389 fprintf(fplog, "%s", buf);
5390 fprintf(stderr, "\n");
5391 fprintf(stderr, "%s", buf);
5392 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5393 fprintf(fplog, "%s", buf);
5394 fprintf(stderr, "%s", buf);
5396 if (comm->bDynLoadBal)
5398 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5399 for (d = 0; d < dd->ndim; d++)
5401 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5402 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5408 sprintf(buf+strlen(buf), "\n");
5409 fprintf(fplog, "%s", buf);
5410 fprintf(stderr, "%s", buf);
5414 pme_f_ratio = comm->load_pme/comm->load_mdf;
5415 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5418 lossp *= (float)npme/(float)nnodes;
5422 lossp *= (float)npp/(float)nnodes;
5424 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5425 fprintf(fplog, "%s", buf);
5426 fprintf(stderr, "%s", buf);
5427 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5428 fprintf(fplog, "%s", buf);
5429 fprintf(stderr, "%s", buf);
5431 fprintf(fplog, "\n");
5432 fprintf(stderr, "\n");
5434 if (lossf >= DD_PERF_LOSS)
5437 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5438 " in the domain decomposition.\n", lossf*100);
5439 if (!comm->bDynLoadBal)
5441 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5445 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5447 fprintf(fplog, "%s\n", buf);
5448 fprintf(stderr, "%s\n", buf);
5450 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5453 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5454 " had %s work to do than the PP nodes.\n"
5455 " You might want to %s the number of PME nodes\n"
5456 " or %s the cut-off and the grid spacing.\n",
5458 (lossp < 0) ? "less" : "more",
5459 (lossp < 0) ? "decrease" : "increase",
5460 (lossp < 0) ? "decrease" : "increase");
5461 fprintf(fplog, "%s\n", buf);
5462 fprintf(stderr, "%s\n", buf);
5467 static float dd_vol_min(gmx_domdec_t *dd)
5469 return dd->comm->load[0].cvol_min*dd->nnodes;
5472 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5474 return dd->comm->load[0].flags;
5477 static float dd_f_imbal(gmx_domdec_t *dd)
5479 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5482 float dd_pme_f_ratio(gmx_domdec_t *dd)
5484 if (dd->comm->cycl_n[ddCyclPME] > 0)
5486 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5494 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5499 flags = dd_load_flags(dd);
5503 "DD load balancing is limited by minimum cell size in dimension");
5504 for (d = 0; d < dd->ndim; d++)
5508 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5511 fprintf(fplog, "\n");
5513 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5514 if (dd->comm->bDynLoadBal)
5516 fprintf(fplog, " vol min/aver %5.3f%c",
5517 dd_vol_min(dd), flags ? '!' : ' ');
5519 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5520 if (dd->comm->cycl_n[ddCyclPME])
5522 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5524 fprintf(fplog, "\n\n");
5527 static void dd_print_load_verbose(gmx_domdec_t *dd)
5529 if (dd->comm->bDynLoadBal)
5531 fprintf(stderr, "vol %4.2f%c ",
5532 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5534 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5535 if (dd->comm->cycl_n[ddCyclPME])
5537 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5542 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5547 gmx_domdec_root_t *root;
5548 gmx_bool bPartOfGroup = FALSE;
5550 dim = dd->dim[dim_ind];
5551 copy_ivec(loc, loc_c);
5552 for (i = 0; i < dd->nc[dim]; i++)
5555 rank = dd_index(dd->nc, loc_c);
5556 if (rank == dd->rank)
5558 /* This process is part of the group */
5559 bPartOfGroup = TRUE;
5562 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5566 dd->comm->mpi_comm_load[dim_ind] = c_row;
5567 if (dd->comm->eDLB != edlbNO)
5569 if (dd->ci[dim] == dd->master_ci[dim])
5571 /* This is the root process of this row */
5572 snew(dd->comm->root[dim_ind], 1);
5573 root = dd->comm->root[dim_ind];
5574 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5575 snew(root->old_cell_f, dd->nc[dim]+1);
5576 snew(root->bCellMin, dd->nc[dim]);
5579 snew(root->cell_f_max0, dd->nc[dim]);
5580 snew(root->cell_f_min1, dd->nc[dim]);
5581 snew(root->bound_min, dd->nc[dim]);
5582 snew(root->bound_max, dd->nc[dim]);
5584 snew(root->buf_ncd, dd->nc[dim]);
5588 /* This is not a root process, we only need to receive cell_f */
5589 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5592 if (dd->ci[dim] == dd->master_ci[dim])
5594 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5600 static void make_load_communicators(gmx_domdec_t *dd)
5603 int dim0, dim1, i, j;
5608 fprintf(debug, "Making load communicators\n");
5611 snew(dd->comm->load, dd->ndim);
5612 snew(dd->comm->mpi_comm_load, dd->ndim);
5615 make_load_communicator(dd, 0, loc);
5619 for (i = 0; i < dd->nc[dim0]; i++)
5622 make_load_communicator(dd, 1, loc);
5628 for (i = 0; i < dd->nc[dim0]; i++)
5632 for (j = 0; j < dd->nc[dim1]; j++)
5635 make_load_communicator(dd, 2, loc);
5642 fprintf(debug, "Finished making load communicators\n");
5647 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5650 int d, dim, i, j, m;
5653 ivec dd_zp[DD_MAXIZONE];
5654 gmx_domdec_zones_t *zones;
5655 gmx_domdec_ns_ranges_t *izone;
5657 for (d = 0; d < dd->ndim; d++)
5660 copy_ivec(dd->ci, tmp);
5661 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5662 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5663 copy_ivec(dd->ci, tmp);
5664 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5665 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5668 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5671 dd->neighbor[d][1]);
5677 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5679 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5680 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5687 for (i = 0; i < nzonep; i++)
5689 copy_ivec(dd_zp3[i], dd_zp[i]);
5695 for (i = 0; i < nzonep; i++)
5697 copy_ivec(dd_zp2[i], dd_zp[i]);
5703 for (i = 0; i < nzonep; i++)
5705 copy_ivec(dd_zp1[i], dd_zp[i]);
5709 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5714 zones = &dd->comm->zones;
5716 for (i = 0; i < nzone; i++)
5719 clear_ivec(zones->shift[i]);
5720 for (d = 0; d < dd->ndim; d++)
5722 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5727 for (i = 0; i < nzone; i++)
5729 for (d = 0; d < DIM; d++)
5731 s[d] = dd->ci[d] - zones->shift[i][d];
5736 else if (s[d] >= dd->nc[d])
5742 zones->nizone = nzonep;
5743 for (i = 0; i < zones->nizone; i++)
5745 if (dd_zp[i][0] != i)
5747 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5749 izone = &zones->izone[i];
5750 izone->j0 = dd_zp[i][1];
5751 izone->j1 = dd_zp[i][2];
5752 for (dim = 0; dim < DIM; dim++)
5754 if (dd->nc[dim] == 1)
5756 /* All shifts should be allowed */
5757 izone->shift0[dim] = -1;
5758 izone->shift1[dim] = 1;
5763 izone->shift0[d] = 0;
5764 izone->shift1[d] = 0;
5765 for(j=izone->j0; j<izone->j1; j++) {
5766 if (dd->shift[j][d] > dd->shift[i][d])
5767 izone->shift0[d] = -1;
5768 if (dd->shift[j][d] < dd->shift[i][d])
5769 izone->shift1[d] = 1;
5775 /* Assume the shift are not more than 1 cell */
5776 izone->shift0[dim] = 1;
5777 izone->shift1[dim] = -1;
5778 for (j = izone->j0; j < izone->j1; j++)
5780 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5781 if (shift_diff < izone->shift0[dim])
5783 izone->shift0[dim] = shift_diff;
5785 if (shift_diff > izone->shift1[dim])
5787 izone->shift1[dim] = shift_diff;
5794 if (dd->comm->eDLB != edlbNO)
5796 snew(dd->comm->root, dd->ndim);
5799 if (dd->comm->bRecordLoad)
5801 make_load_communicators(dd);
5805 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5808 gmx_domdec_comm_t *comm;
5819 if (comm->bCartesianPP)
5821 /* Set up cartesian communication for the particle-particle part */
5824 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5825 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5828 for (i = 0; i < DIM; i++)
5832 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5834 /* We overwrite the old communicator with the new cartesian one */
5835 cr->mpi_comm_mygroup = comm_cart;
5838 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5839 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5841 if (comm->bCartesianPP_PME)
5843 /* Since we want to use the original cartesian setup for sim,
5844 * and not the one after split, we need to make an index.
5846 snew(comm->ddindex2ddnodeid, dd->nnodes);
5847 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5848 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5849 /* Get the rank of the DD master,
5850 * above we made sure that the master node is a PP node.
5860 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5862 else if (comm->bCartesianPP)
5864 if (cr->npmenodes == 0)
5866 /* The PP communicator is also
5867 * the communicator for this simulation
5869 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5871 cr->nodeid = dd->rank;
5873 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5875 /* We need to make an index to go from the coordinates
5876 * to the nodeid of this simulation.
5878 snew(comm->ddindex2simnodeid, dd->nnodes);
5879 snew(buf, dd->nnodes);
5880 if (cr->duty & DUTY_PP)
5882 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5884 /* Communicate the ddindex to simulation nodeid index */
5885 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5886 cr->mpi_comm_mysim);
5889 /* Determine the master coordinates and rank.
5890 * The DD master should be the same node as the master of this sim.
5892 for (i = 0; i < dd->nnodes; i++)
5894 if (comm->ddindex2simnodeid[i] == 0)
5896 ddindex2xyz(dd->nc, i, dd->master_ci);
5897 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5902 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5907 /* No Cartesian communicators */
5908 /* We use the rank in dd->comm->all as DD index */
5909 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5910 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5912 clear_ivec(dd->master_ci);
5919 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5920 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5925 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5926 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5930 static void receive_ddindex2simnodeid(t_commrec *cr)
5934 gmx_domdec_comm_t *comm;
5941 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5943 snew(comm->ddindex2simnodeid, dd->nnodes);
5944 snew(buf, dd->nnodes);
5945 if (cr->duty & DUTY_PP)
5947 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5950 /* Communicate the ddindex to simulation nodeid index */
5951 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5952 cr->mpi_comm_mysim);
5959 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5960 int ncg, int natoms)
5962 gmx_domdec_master_t *ma;
5967 snew(ma->ncg, dd->nnodes);
5968 snew(ma->index, dd->nnodes+1);
5970 snew(ma->nat, dd->nnodes);
5971 snew(ma->ibuf, dd->nnodes*2);
5972 snew(ma->cell_x, DIM);
5973 for (i = 0; i < DIM; i++)
5975 snew(ma->cell_x[i], dd->nc[i]+1);
5978 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5984 snew(ma->vbuf, natoms);
5990 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
5994 gmx_domdec_comm_t *comm;
6005 if (comm->bCartesianPP)
6007 for (i = 1; i < DIM; i++)
6009 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6011 if (bDiv[YY] || bDiv[ZZ])
6013 comm->bCartesianPP_PME = TRUE;
6014 /* If we have 2D PME decomposition, which is always in x+y,
6015 * we stack the PME only nodes in z.
6016 * Otherwise we choose the direction that provides the thinnest slab
6017 * of PME only nodes as this will have the least effect
6018 * on the PP communication.
6019 * But for the PME communication the opposite might be better.
6021 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6023 dd->nc[YY] > dd->nc[ZZ]))
6025 comm->cartpmedim = ZZ;
6029 comm->cartpmedim = YY;
6031 comm->ntot[comm->cartpmedim]
6032 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6036 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]);
6038 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6043 if (comm->bCartesianPP_PME)
6047 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]);
6050 for (i = 0; i < DIM; i++)
6054 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6057 MPI_Comm_rank(comm_cart, &rank);
6058 if (MASTERNODE(cr) && rank != 0)
6060 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6063 /* With this assigment we loose the link to the original communicator
6064 * which will usually be MPI_COMM_WORLD, unless have multisim.
6066 cr->mpi_comm_mysim = comm_cart;
6067 cr->sim_nodeid = rank;
6069 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6073 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6074 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6077 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6081 if (cr->npmenodes == 0 ||
6082 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6084 cr->duty = DUTY_PME;
6087 /* Split the sim communicator into PP and PME only nodes */
6088 MPI_Comm_split(cr->mpi_comm_mysim,
6090 dd_index(comm->ntot, dd->ci),
6091 &cr->mpi_comm_mygroup);
6095 switch (dd_node_order)
6100 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6103 case ddnoINTERLEAVE:
6104 /* Interleave the PP-only and PME-only nodes,
6105 * as on clusters with dual-core machines this will double
6106 * the communication bandwidth of the PME processes
6107 * and thus speed up the PP <-> PME and inter PME communication.
6111 fprintf(fplog, "Interleaving PP and PME nodes\n");
6113 comm->pmenodes = dd_pmenodes(cr);
6118 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6121 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6123 cr->duty = DUTY_PME;
6130 /* Split the sim communicator into PP and PME only nodes */
6131 MPI_Comm_split(cr->mpi_comm_mysim,
6134 &cr->mpi_comm_mygroup);
6135 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6141 fprintf(fplog, "This is a %s only node\n\n",
6142 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6146 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6149 gmx_domdec_comm_t *comm;
6155 copy_ivec(dd->nc, comm->ntot);
6157 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6158 comm->bCartesianPP_PME = FALSE;
6160 /* Reorder the nodes by default. This might change the MPI ranks.
6161 * Real reordering is only supported on very few architectures,
6162 * Blue Gene is one of them.
6164 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6166 if (cr->npmenodes > 0)
6168 /* Split the communicator into a PP and PME part */
6169 split_communicator(fplog, cr, dd_node_order, CartReorder);
6170 if (comm->bCartesianPP_PME)
6172 /* We (possibly) reordered the nodes in split_communicator,
6173 * so it is no longer required in make_pp_communicator.
6175 CartReorder = FALSE;
6180 /* All nodes do PP and PME */
6182 /* We do not require separate communicators */
6183 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6187 if (cr->duty & DUTY_PP)
6189 /* Copy or make a new PP communicator */
6190 make_pp_communicator(fplog, cr, CartReorder);
6194 receive_ddindex2simnodeid(cr);
6197 if (!(cr->duty & DUTY_PME))
6199 /* Set up the commnuication to our PME node */
6200 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6201 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6204 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6205 dd->pme_nodeid, dd->pme_receive_vir_ener);
6210 dd->pme_nodeid = -1;
6215 dd->ma = init_gmx_domdec_master_t(dd,
6217 comm->cgs_gl.index[comm->cgs_gl.nr]);
6221 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6223 real *slb_frac, tot;
6228 if (nc > 1 && size_string != NULL)
6232 fprintf(fplog, "Using static load balancing for the %s direction\n",
6237 for (i = 0; i < nc; i++)
6240 sscanf(size_string, "%lf%n", &dbl, &n);
6243 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6252 fprintf(fplog, "Relative cell sizes:");
6254 for (i = 0; i < nc; i++)
6259 fprintf(fplog, " %5.3f", slb_frac[i]);
6264 fprintf(fplog, "\n");
6271 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6274 gmx_mtop_ilistloop_t iloop;
6278 iloop = gmx_mtop_ilistloop_init(mtop);
6279 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6281 for (ftype = 0; ftype < F_NRE; ftype++)
6283 if ((interaction_function[ftype].flags & IF_BOND) &&
6286 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6294 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6300 val = getenv(env_var);
6303 if (sscanf(val, "%d", &nst) <= 0)
6309 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6317 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6321 fprintf(stderr, "\n%s\n", warn_string);
6325 fprintf(fplog, "\n%s\n", warn_string);
6329 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6330 t_inputrec *ir, FILE *fplog)
6332 if (ir->ePBC == epbcSCREW &&
6333 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6335 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6338 if (ir->ns_type == ensSIMPLE)
6340 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6343 if (ir->nstlist == 0)
6345 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6348 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6350 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6354 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6359 r = ddbox->box_size[XX];
6360 for (di = 0; di < dd->ndim; di++)
6363 /* Check using the initial average cell size */
6364 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6370 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6371 const char *dlb_opt, gmx_bool bRecordLoad,
6372 unsigned long Flags, t_inputrec *ir)
6380 case 'a': eDLB = edlbAUTO; break;
6381 case 'n': eDLB = edlbNO; break;
6382 case 'y': eDLB = edlbYES; break;
6383 default: gmx_incons("Unknown dlb_opt");
6386 if (Flags & MD_RERUN)
6391 if (!EI_DYNAMICS(ir->eI))
6393 if (eDLB == edlbYES)
6395 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6396 dd_warning(cr, fplog, buf);
6404 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6409 if (Flags & MD_REPRODUCIBLE)
6416 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6420 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6423 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6431 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6436 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6438 /* Decomposition order z,y,x */
6441 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6443 for (dim = DIM-1; dim >= 0; dim--)
6445 if (dd->nc[dim] > 1)
6447 dd->dim[dd->ndim++] = dim;
6453 /* Decomposition order x,y,z */
6454 for (dim = 0; dim < DIM; dim++)
6456 if (dd->nc[dim] > 1)
6458 dd->dim[dd->ndim++] = dim;
6464 static gmx_domdec_comm_t *init_dd_comm()
6466 gmx_domdec_comm_t *comm;
6470 snew(comm->cggl_flag, DIM*2);
6471 snew(comm->cgcm_state, DIM*2);
6472 for (i = 0; i < DIM*2; i++)
6474 comm->cggl_flag_nalloc[i] = 0;
6475 comm->cgcm_state_nalloc[i] = 0;
6478 comm->nalloc_int = 0;
6479 comm->buf_int = NULL;
6481 vec_rvec_init(&comm->vbuf);
6483 comm->n_load_have = 0;
6484 comm->n_load_collect = 0;
6486 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6488 comm->sum_nat[i] = 0;
6492 comm->load_step = 0;
6495 clear_ivec(comm->load_lim);
6502 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6503 unsigned long Flags,
6505 real comm_distance_min, real rconstr,
6506 const char *dlb_opt, real dlb_scale,
6507 const char *sizex, const char *sizey, const char *sizez,
6508 gmx_mtop_t *mtop, t_inputrec *ir,
6509 matrix box, rvec *x,
6511 int *npme_x, int *npme_y)
6514 gmx_domdec_comm_t *comm;
6517 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6524 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6529 dd->comm = init_dd_comm();
6531 snew(comm->cggl_flag, DIM*2);
6532 snew(comm->cgcm_state, DIM*2);
6534 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6535 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6537 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6538 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6539 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6540 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6541 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6542 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6543 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6544 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6546 dd->pme_recv_f_alloc = 0;
6547 dd->pme_recv_f_buf = NULL;
6549 if (dd->bSendRecv2 && fplog)
6551 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");
6557 fprintf(fplog, "Will load balance based on FLOP count\n");
6559 if (comm->eFlop > 1)
6561 srand(1+cr->nodeid);
6563 comm->bRecordLoad = TRUE;
6567 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6571 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6573 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6576 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6578 dd->bGridJump = comm->bDynLoadBal;
6579 comm->bPMELoadBalDLBLimits = FALSE;
6581 if (comm->nstSortCG)
6585 if (comm->nstSortCG == 1)
6587 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6591 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6595 snew(comm->sort, 1);
6601 fprintf(fplog, "Will not sort the charge groups\n");
6605 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6607 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6608 if (comm->bInterCGBondeds)
6610 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6614 comm->bInterCGMultiBody = FALSE;
6617 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6618 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6620 if (ir->rlistlong == 0)
6622 /* Set the cut-off to some very large value,
6623 * so we don't need if statements everywhere in the code.
6624 * We use sqrt, since the cut-off is squared in some places.
6626 comm->cutoff = GMX_CUTOFF_INF;
6630 comm->cutoff = ir->rlistlong;
6632 comm->cutoff_mbody = 0;
6634 comm->cellsize_limit = 0;
6635 comm->bBondComm = FALSE;
6637 if (comm->bInterCGBondeds)
6639 if (comm_distance_min > 0)
6641 comm->cutoff_mbody = comm_distance_min;
6642 if (Flags & MD_DDBONDCOMM)
6644 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6648 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6650 r_bonded_limit = comm->cutoff_mbody;
6652 else if (ir->bPeriodicMols)
6654 /* Can not easily determine the required cut-off */
6655 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");
6656 comm->cutoff_mbody = comm->cutoff/2;
6657 r_bonded_limit = comm->cutoff_mbody;
6663 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6664 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6666 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6667 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6669 /* We use an initial margin of 10% for the minimum cell size,
6670 * except when we are just below the non-bonded cut-off.
6672 if (Flags & MD_DDBONDCOMM)
6674 if (max(r_2b, r_mb) > comm->cutoff)
6676 r_bonded = max(r_2b, r_mb);
6677 r_bonded_limit = 1.1*r_bonded;
6678 comm->bBondComm = TRUE;
6683 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6685 /* We determine cutoff_mbody later */
6689 /* No special bonded communication,
6690 * simply increase the DD cut-off.
6692 r_bonded_limit = 1.1*max(r_2b, r_mb);
6693 comm->cutoff_mbody = r_bonded_limit;
6694 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6697 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6701 "Minimum cell size due to bonded interactions: %.3f nm\n",
6702 comm->cellsize_limit);
6706 if (dd->bInterCGcons && rconstr <= 0)
6708 /* There is a cell size limit due to the constraints (P-LINCS) */
6709 rconstr = constr_r_max(fplog, mtop, ir);
6713 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6715 if (rconstr > comm->cellsize_limit)
6717 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6721 else if (rconstr > 0 && fplog)
6723 /* Here we do not check for dd->bInterCGcons,
6724 * because one can also set a cell size limit for virtual sites only
6725 * and at this point we don't know yet if there are intercg v-sites.
6728 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6731 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6733 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6737 copy_ivec(nc, dd->nc);
6738 set_dd_dim(fplog, dd);
6739 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6741 if (cr->npmenodes == -1)
6745 acs = average_cellsize_min(dd, ddbox);
6746 if (acs < comm->cellsize_limit)
6750 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6752 gmx_fatal_collective(FARGS, cr, NULL,
6753 "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",
6754 acs, comm->cellsize_limit);
6759 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6761 /* We need to choose the optimal DD grid and possibly PME nodes */
6762 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6763 comm->eDLB != edlbNO, dlb_scale,
6764 comm->cellsize_limit, comm->cutoff,
6765 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6767 if (dd->nc[XX] == 0)
6769 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6770 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6771 !bC ? "-rdd" : "-rcon",
6772 comm->eDLB != edlbNO ? " or -dds" : "",
6773 bC ? " or your LINCS settings" : "");
6775 gmx_fatal_collective(FARGS, cr, NULL,
6776 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6778 "Look in the log file for details on the domain decomposition",
6779 cr->nnodes-cr->npmenodes, limit, buf);
6781 set_dd_dim(fplog, dd);
6787 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6788 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6791 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6792 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6794 gmx_fatal_collective(FARGS, cr, NULL,
6795 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6796 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6798 if (cr->npmenodes > dd->nnodes)
6800 gmx_fatal_collective(FARGS, cr, NULL,
6801 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6803 if (cr->npmenodes > 0)
6805 comm->npmenodes = cr->npmenodes;
6809 comm->npmenodes = dd->nnodes;
6812 if (EEL_PME(ir->coulombtype))
6814 /* The following choices should match those
6815 * in comm_cost_est in domdec_setup.c.
6816 * Note that here the checks have to take into account
6817 * that the decomposition might occur in a different order than xyz
6818 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6819 * in which case they will not match those in comm_cost_est,
6820 * but since that is mainly for testing purposes that's fine.
6822 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6823 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6824 getenv("GMX_PMEONEDD") == NULL)
6826 comm->npmedecompdim = 2;
6827 comm->npmenodes_x = dd->nc[XX];
6828 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6832 /* In case nc is 1 in both x and y we could still choose to
6833 * decompose pme in y instead of x, but we use x for simplicity.
6835 comm->npmedecompdim = 1;
6836 if (dd->dim[0] == YY)
6838 comm->npmenodes_x = 1;
6839 comm->npmenodes_y = comm->npmenodes;
6843 comm->npmenodes_x = comm->npmenodes;
6844 comm->npmenodes_y = 1;
6849 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6850 comm->npmenodes_x, comm->npmenodes_y, 1);
6855 comm->npmedecompdim = 0;
6856 comm->npmenodes_x = 0;
6857 comm->npmenodes_y = 0;
6860 /* Technically we don't need both of these,
6861 * but it simplifies code not having to recalculate it.
6863 *npme_x = comm->npmenodes_x;
6864 *npme_y = comm->npmenodes_y;
6866 snew(comm->slb_frac, DIM);
6867 if (comm->eDLB == edlbNO)
6869 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6870 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6871 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6874 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6876 if (comm->bBondComm || comm->eDLB != edlbNO)
6878 /* Set the bonded communication distance to halfway
6879 * the minimum and the maximum,
6880 * since the extra communication cost is nearly zero.
6882 acs = average_cellsize_min(dd, ddbox);
6883 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6884 if (comm->eDLB != edlbNO)
6886 /* Check if this does not limit the scaling */
6887 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6889 if (!comm->bBondComm)
6891 /* Without bBondComm do not go beyond the n.b. cut-off */
6892 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6893 if (comm->cellsize_limit >= comm->cutoff)
6895 /* We don't loose a lot of efficieny
6896 * when increasing it to the n.b. cut-off.
6897 * It can even be slightly faster, because we need
6898 * less checks for the communication setup.
6900 comm->cutoff_mbody = comm->cutoff;
6903 /* Check if we did not end up below our original limit */
6904 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6906 if (comm->cutoff_mbody > comm->cellsize_limit)
6908 comm->cellsize_limit = comm->cutoff_mbody;
6911 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6916 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6917 "cellsize limit %f\n",
6918 comm->bBondComm, comm->cellsize_limit);
6923 check_dd_restrictions(cr, dd, ir, fplog);
6926 comm->partition_step = INT_MIN;
6929 clear_dd_cycle_counts(dd);
6934 static void set_dlb_limits(gmx_domdec_t *dd)
6939 for (d = 0; d < dd->ndim; d++)
6941 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6942 dd->comm->cellsize_min[dd->dim[d]] =
6943 dd->comm->cellsize_min_dlb[dd->dim[d]];
6948 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6951 gmx_domdec_comm_t *comm;
6961 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);
6964 cellsize_min = comm->cellsize_min[dd->dim[0]];
6965 for (d = 1; d < dd->ndim; d++)
6967 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6970 if (cellsize_min < comm->cellsize_limit*1.05)
6972 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");
6974 /* Change DLB from "auto" to "no". */
6975 comm->eDLB = edlbNO;
6980 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6981 comm->bDynLoadBal = TRUE;
6982 dd->bGridJump = TRUE;
6986 /* We can set the required cell size info here,
6987 * so we do not need to communicate this.
6988 * The grid is completely uniform.
6990 for (d = 0; d < dd->ndim; d++)
6994 comm->load[d].sum_m = comm->load[d].sum;
6996 nc = dd->nc[dd->dim[d]];
6997 for (i = 0; i < nc; i++)
6999 comm->root[d]->cell_f[i] = i/(real)nc;
7002 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7003 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7006 comm->root[d]->cell_f[nc] = 1.0;
7011 static char *init_bLocalCG(gmx_mtop_t *mtop)
7016 ncg = ncg_mtop(mtop);
7017 snew(bLocalCG, ncg);
7018 for (cg = 0; cg < ncg; cg++)
7020 bLocalCG[cg] = FALSE;
7026 void dd_init_bondeds(FILE *fplog,
7027 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7028 gmx_vsite_t *vsite, gmx_constr_t constr,
7029 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7031 gmx_domdec_comm_t *comm;
7035 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7039 if (comm->bBondComm)
7041 /* Communicate atoms beyond the cut-off for bonded interactions */
7044 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7046 comm->bLocalCG = init_bLocalCG(mtop);
7050 /* Only communicate atoms based on cut-off */
7051 comm->cglink = NULL;
7052 comm->bLocalCG = NULL;
7056 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7058 gmx_bool bDynLoadBal, real dlb_scale,
7061 gmx_domdec_comm_t *comm;
7076 fprintf(fplog, "The maximum number of communication pulses is:");
7077 for (d = 0; d < dd->ndim; d++)
7079 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7081 fprintf(fplog, "\n");
7082 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7083 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7084 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7085 for (d = 0; d < DIM; d++)
7089 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7096 comm->cellsize_min_dlb[d]/
7097 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7099 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7102 fprintf(fplog, "\n");
7106 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7107 fprintf(fplog, "The initial number of communication pulses is:");
7108 for (d = 0; d < dd->ndim; d++)
7110 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7112 fprintf(fplog, "\n");
7113 fprintf(fplog, "The initial domain decomposition cell size is:");
7114 for (d = 0; d < DIM; d++)
7118 fprintf(fplog, " %c %.2f nm",
7119 dim2char(d), dd->comm->cellsize_min[d]);
7122 fprintf(fplog, "\n\n");
7125 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7127 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7128 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7129 "non-bonded interactions", "", comm->cutoff);
7133 limit = dd->comm->cellsize_limit;
7137 if (dynamic_dd_box(ddbox, ir))
7139 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7141 limit = dd->comm->cellsize_min[XX];
7142 for (d = 1; d < DIM; d++)
7144 limit = min(limit, dd->comm->cellsize_min[d]);
7148 if (comm->bInterCGBondeds)
7150 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7151 "two-body bonded interactions", "(-rdd)",
7152 max(comm->cutoff, comm->cutoff_mbody));
7153 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7154 "multi-body bonded interactions", "(-rdd)",
7155 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7159 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7160 "virtual site constructions", "(-rcon)", limit);
7162 if (dd->constraint_comm)
7164 sprintf(buf, "atoms separated by up to %d constraints",
7166 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7167 buf, "(-rcon)", limit);
7169 fprintf(fplog, "\n");
7175 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7177 const t_inputrec *ir,
7178 const gmx_ddbox_t *ddbox)
7180 gmx_domdec_comm_t *comm;
7181 int d, dim, npulse, npulse_d_max, npulse_d;
7186 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7188 /* Determine the maximum number of comm. pulses in one dimension */
7190 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7192 /* Determine the maximum required number of grid pulses */
7193 if (comm->cellsize_limit >= comm->cutoff)
7195 /* Only a single pulse is required */
7198 else if (!bNoCutOff && comm->cellsize_limit > 0)
7200 /* We round down slightly here to avoid overhead due to the latency
7201 * of extra communication calls when the cut-off
7202 * would be only slightly longer than the cell size.
7203 * Later cellsize_limit is redetermined,
7204 * so we can not miss interactions due to this rounding.
7206 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7210 /* There is no cell size limit */
7211 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7214 if (!bNoCutOff && npulse > 1)
7216 /* See if we can do with less pulses, based on dlb_scale */
7218 for (d = 0; d < dd->ndim; d++)
7221 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7222 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7223 npulse_d_max = max(npulse_d_max, npulse_d);
7225 npulse = min(npulse, npulse_d_max);
7228 /* This env var can override npulse */
7229 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7236 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7237 for (d = 0; d < dd->ndim; d++)
7239 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7240 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7241 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7242 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7243 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7245 comm->bVacDLBNoLimit = FALSE;
7249 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7250 if (!comm->bVacDLBNoLimit)
7252 comm->cellsize_limit = max(comm->cellsize_limit,
7253 comm->cutoff/comm->maxpulse);
7255 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7256 /* Set the minimum cell size for each DD dimension */
7257 for (d = 0; d < dd->ndim; d++)
7259 if (comm->bVacDLBNoLimit ||
7260 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7262 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7266 comm->cellsize_min_dlb[dd->dim[d]] =
7267 comm->cutoff/comm->cd[d].np_dlb;
7270 if (comm->cutoff_mbody <= 0)
7272 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7274 if (comm->bDynLoadBal)
7280 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7282 /* If each molecule is a single charge group
7283 * or we use domain decomposition for each periodic dimension,
7284 * we do not need to take pbc into account for the bonded interactions.
7286 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7289 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7292 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7293 t_inputrec *ir, t_forcerec *fr,
7296 gmx_domdec_comm_t *comm;
7302 /* Initialize the thread data.
7303 * This can not be done in init_domain_decomposition,
7304 * as the numbers of threads is determined later.
7306 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7309 snew(comm->dth, comm->nth);
7312 if (EEL_PME(ir->coulombtype))
7314 init_ddpme(dd, &comm->ddpme[0], 0);
7315 if (comm->npmedecompdim >= 2)
7317 init_ddpme(dd, &comm->ddpme[1], 1);
7322 comm->npmenodes = 0;
7323 if (dd->pme_nodeid >= 0)
7325 gmx_fatal_collective(FARGS, NULL, dd,
7326 "Can not have separate PME nodes without PME electrostatics");
7332 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7334 if (comm->eDLB != edlbNO)
7336 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7339 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7340 if (comm->eDLB == edlbAUTO)
7344 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7346 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7349 if (ir->ePBC == epbcNONE)
7351 vol_frac = 1 - 1/(double)dd->nnodes;
7356 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7360 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7362 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7364 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7367 static gmx_bool test_dd_cutoff(t_commrec *cr,
7368 t_state *state, t_inputrec *ir,
7379 set_ddbox(dd, FALSE, cr, ir, state->box,
7380 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7384 for (d = 0; d < dd->ndim; d++)
7388 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7389 if (dynamic_dd_box(&ddbox, ir))
7391 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7394 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7396 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7397 dd->comm->cd[d].np_dlb > 0)
7399 if (np > dd->comm->cd[d].np_dlb)
7404 /* If a current local cell size is smaller than the requested
7405 * cut-off, we could still fix it, but this gets very complicated.
7406 * Without fixing here, we might actually need more checks.
7408 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7415 if (dd->comm->eDLB != edlbNO)
7417 /* If DLB is not active yet, we don't need to check the grid jumps.
7418 * Actually we shouldn't, because then the grid jump data is not set.
7420 if (dd->comm->bDynLoadBal &&
7421 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7426 gmx_sumi(1, &LocallyLimited, cr);
7428 if (LocallyLimited > 0)
7437 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7440 gmx_bool bCutoffAllowed;
7442 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7446 cr->dd->comm->cutoff = cutoff_req;
7449 return bCutoffAllowed;
7452 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7454 gmx_domdec_comm_t *comm;
7456 comm = cr->dd->comm;
7458 /* Turn on the DLB limiting (might have been on already) */
7459 comm->bPMELoadBalDLBLimits = TRUE;
7461 /* Change the cut-off limit */
7462 comm->PMELoadBal_max_cutoff = comm->cutoff;
7465 static void merge_cg_buffers(int ncell,
7466 gmx_domdec_comm_dim_t *cd, int pulse,
7468 int *index_gl, int *recv_i,
7469 rvec *cg_cm, rvec *recv_vr,
7471 cginfo_mb_t *cginfo_mb, int *cginfo)
7473 gmx_domdec_ind_t *ind, *ind_p;
7474 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7475 int shift, shift_at;
7477 ind = &cd->ind[pulse];
7479 /* First correct the already stored data */
7480 shift = ind->nrecv[ncell];
7481 for (cell = ncell-1; cell >= 0; cell--)
7483 shift -= ind->nrecv[cell];
7486 /* Move the cg's present from previous grid pulses */
7487 cg0 = ncg_cell[ncell+cell];
7488 cg1 = ncg_cell[ncell+cell+1];
7489 cgindex[cg1+shift] = cgindex[cg1];
7490 for (cg = cg1-1; cg >= cg0; cg--)
7492 index_gl[cg+shift] = index_gl[cg];
7493 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7494 cgindex[cg+shift] = cgindex[cg];
7495 cginfo[cg+shift] = cginfo[cg];
7497 /* Correct the already stored send indices for the shift */
7498 for (p = 1; p <= pulse; p++)
7500 ind_p = &cd->ind[p];
7502 for (c = 0; c < cell; c++)
7504 cg0 += ind_p->nsend[c];
7506 cg1 = cg0 + ind_p->nsend[cell];
7507 for (cg = cg0; cg < cg1; cg++)
7509 ind_p->index[cg] += shift;
7515 /* Merge in the communicated buffers */
7519 for (cell = 0; cell < ncell; cell++)
7521 cg1 = ncg_cell[ncell+cell+1] + shift;
7524 /* Correct the old cg indices */
7525 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7527 cgindex[cg+1] += shift_at;
7530 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7532 /* Copy this charge group from the buffer */
7533 index_gl[cg1] = recv_i[cg0];
7534 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7535 /* Add it to the cgindex */
7536 cg_gl = index_gl[cg1];
7537 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7538 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7539 cgindex[cg1+1] = cgindex[cg1] + nat;
7544 shift += ind->nrecv[cell];
7545 ncg_cell[ncell+cell+1] = cg1;
7549 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7550 int nzone, int cg0, const int *cgindex)
7554 /* Store the atom block boundaries for easy copying of communication buffers
7557 for (zone = 0; zone < nzone; zone++)
7559 for (p = 0; p < cd->np; p++)
7561 cd->ind[p].cell2at0[zone] = cgindex[cg];
7562 cg += cd->ind[p].nrecv[zone];
7563 cd->ind[p].cell2at1[zone] = cgindex[cg];
7568 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7574 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7576 if (!bLocalCG[link->a[i]])
7585 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7587 real c[DIM][4]; /* the corners for the non-bonded communication */
7588 real cr0; /* corner for rounding */
7589 real cr1[4]; /* corners for rounding */
7590 real bc[DIM]; /* corners for bounded communication */
7591 real bcr1; /* corner for rounding for bonded communication */
7594 /* Determine the corners of the domain(s) we are communicating with */
7596 set_dd_corners(const gmx_domdec_t *dd,
7597 int dim0, int dim1, int dim2,
7601 const gmx_domdec_comm_t *comm;
7602 const gmx_domdec_zones_t *zones;
7607 zones = &comm->zones;
7609 /* Keep the compiler happy */
7613 /* The first dimension is equal for all cells */
7614 c->c[0][0] = comm->cell_x0[dim0];
7617 c->bc[0] = c->c[0][0];
7622 /* This cell row is only seen from the first row */
7623 c->c[1][0] = comm->cell_x0[dim1];
7624 /* All rows can see this row */
7625 c->c[1][1] = comm->cell_x0[dim1];
7628 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7631 /* For the multi-body distance we need the maximum */
7632 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7635 /* Set the upper-right corner for rounding */
7636 c->cr0 = comm->cell_x1[dim0];
7641 for (j = 0; j < 4; j++)
7643 c->c[2][j] = comm->cell_x0[dim2];
7647 /* Use the maximum of the i-cells that see a j-cell */
7648 for (i = 0; i < zones->nizone; i++)
7650 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7656 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7662 /* For the multi-body distance we need the maximum */
7663 c->bc[2] = comm->cell_x0[dim2];
7664 for (i = 0; i < 2; i++)
7666 for (j = 0; j < 2; j++)
7668 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7674 /* Set the upper-right corner for rounding */
7675 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7676 * Only cell (0,0,0) can see cell 7 (1,1,1)
7678 c->cr1[0] = comm->cell_x1[dim1];
7679 c->cr1[3] = comm->cell_x1[dim1];
7682 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7685 /* For the multi-body distance we need the maximum */
7686 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7693 /* Determine which cg's we need to send in this pulse from this zone */
7695 get_zone_pulse_cgs(gmx_domdec_t *dd,
7696 int zonei, int zone,
7698 const int *index_gl,
7700 int dim, int dim_ind,
7701 int dim0, int dim1, int dim2,
7702 real r_comm2, real r_bcomm2,
7706 real skew_fac2_d, real skew_fac_01,
7707 rvec *v_d, rvec *v_0, rvec *v_1,
7708 const dd_corners_t *c,
7710 gmx_bool bDistBonded,
7716 gmx_domdec_ind_t *ind,
7717 int **ibuf, int *ibuf_nalloc,
7723 gmx_domdec_comm_t *comm;
7725 gmx_bool bDistMB_pulse;
7727 real r2, rb2, r, tric_sh;
7730 int nsend_z, nsend, nat;
7734 bScrew = (dd->bScrewPBC && dim == XX);
7736 bDistMB_pulse = (bDistMB && bDistBonded);
7742 for (cg = cg0; cg < cg1; cg++)
7746 if (tric_dist[dim_ind] == 0)
7748 /* Rectangular direction, easy */
7749 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7756 r = cg_cm[cg][dim] - c->bc[dim_ind];
7762 /* Rounding gives at most a 16% reduction
7763 * in communicated atoms
7765 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7767 r = cg_cm[cg][dim0] - c->cr0;
7768 /* This is the first dimension, so always r >= 0 */
7775 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7777 r = cg_cm[cg][dim1] - c->cr1[zone];
7784 r = cg_cm[cg][dim1] - c->bcr1;
7794 /* Triclinic direction, more complicated */
7797 /* Rounding, conservative as the skew_fac multiplication
7798 * will slightly underestimate the distance.
7800 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7802 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7803 for (i = dim0+1; i < DIM; i++)
7805 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7807 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7810 rb[dim0] = rn[dim0];
7813 /* Take care that the cell planes along dim0 might not
7814 * be orthogonal to those along dim1 and dim2.
7816 for (i = 1; i <= dim_ind; i++)
7819 if (normal[dim0][dimd] > 0)
7821 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7824 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7829 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7831 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7833 for (i = dim1+1; i < DIM; i++)
7835 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7837 rn[dim1] += tric_sh;
7840 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7841 /* Take care of coupling of the distances
7842 * to the planes along dim0 and dim1 through dim2.
7844 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7845 /* Take care that the cell planes along dim1
7846 * might not be orthogonal to that along dim2.
7848 if (normal[dim1][dim2] > 0)
7850 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7856 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7859 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7860 /* Take care of coupling of the distances
7861 * to the planes along dim0 and dim1 through dim2.
7863 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7864 /* Take care that the cell planes along dim1
7865 * might not be orthogonal to that along dim2.
7867 if (normal[dim1][dim2] > 0)
7869 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7874 /* The distance along the communication direction */
7875 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7877 for (i = dim+1; i < DIM; i++)
7879 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7884 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7885 /* Take care of coupling of the distances
7886 * to the planes along dim0 and dim1 through dim2.
7888 if (dim_ind == 1 && zonei == 1)
7890 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7896 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7899 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7900 /* Take care of coupling of the distances
7901 * to the planes along dim0 and dim1 through dim2.
7903 if (dim_ind == 1 && zonei == 1)
7905 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7913 ((bDistMB && rb2 < r_bcomm2) ||
7914 (bDist2B && r2 < r_bcomm2)) &&
7916 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7917 missing_link(comm->cglink, index_gl[cg],
7920 /* Make an index to the local charge groups */
7921 if (nsend+1 > ind->nalloc)
7923 ind->nalloc = over_alloc_large(nsend+1);
7924 srenew(ind->index, ind->nalloc);
7926 if (nsend+1 > *ibuf_nalloc)
7928 *ibuf_nalloc = over_alloc_large(nsend+1);
7929 srenew(*ibuf, *ibuf_nalloc);
7931 ind->index[nsend] = cg;
7932 (*ibuf)[nsend] = index_gl[cg];
7934 vec_rvec_check_alloc(vbuf, nsend+1);
7936 if (dd->ci[dim] == 0)
7938 /* Correct cg_cm for pbc */
7939 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7942 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7943 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7948 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7951 nat += cgindex[cg+1] - cgindex[cg];
7957 *nsend_z_ptr = nsend_z;
7960 static void setup_dd_communication(gmx_domdec_t *dd,
7961 matrix box, gmx_ddbox_t *ddbox,
7962 t_forcerec *fr, t_state *state, rvec **f)
7964 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7965 int nzone, nzone_send, zone, zonei, cg0, cg1;
7966 int c, i, j, cg, cg_gl, nrcg;
7967 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7968 gmx_domdec_comm_t *comm;
7969 gmx_domdec_zones_t *zones;
7970 gmx_domdec_comm_dim_t *cd;
7971 gmx_domdec_ind_t *ind;
7972 cginfo_mb_t *cginfo_mb;
7973 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7974 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7975 dd_corners_t corners;
7977 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7978 real skew_fac2_d, skew_fac_01;
7985 fprintf(debug, "Setting up DD communication\n");
7990 switch (fr->cutoff_scheme)
7999 gmx_incons("unimplemented");
8003 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8005 dim = dd->dim[dim_ind];
8007 /* Check if we need to use triclinic distances */
8008 tric_dist[dim_ind] = 0;
8009 for (i = 0; i <= dim_ind; i++)
8011 if (ddbox->tric_dir[dd->dim[i]])
8013 tric_dist[dim_ind] = 1;
8018 bBondComm = comm->bBondComm;
8020 /* Do we need to determine extra distances for multi-body bondeds? */
8021 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8023 /* Do we need to determine extra distances for only two-body bondeds? */
8024 bDist2B = (bBondComm && !bDistMB);
8026 r_comm2 = sqr(comm->cutoff);
8027 r_bcomm2 = sqr(comm->cutoff_mbody);
8031 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8034 zones = &comm->zones;
8037 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8038 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8040 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8042 /* Triclinic stuff */
8043 normal = ddbox->normal;
8047 v_0 = ddbox->v[dim0];
8048 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8050 /* Determine the coupling coefficient for the distances
8051 * to the cell planes along dim0 and dim1 through dim2.
8052 * This is required for correct rounding.
8055 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8058 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8064 v_1 = ddbox->v[dim1];
8067 zone_cg_range = zones->cg_range;
8068 index_gl = dd->index_gl;
8069 cgindex = dd->cgindex;
8070 cginfo_mb = fr->cginfo_mb;
8072 zone_cg_range[0] = 0;
8073 zone_cg_range[1] = dd->ncg_home;
8074 comm->zone_ncg1[0] = dd->ncg_home;
8075 pos_cg = dd->ncg_home;
8077 nat_tot = dd->nat_home;
8079 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8081 dim = dd->dim[dim_ind];
8082 cd = &comm->cd[dim_ind];
8084 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8086 /* No pbc in this dimension, the first node should not comm. */
8094 v_d = ddbox->v[dim];
8095 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8097 cd->bInPlace = TRUE;
8098 for (p = 0; p < cd->np; p++)
8100 /* Only atoms communicated in the first pulse are used
8101 * for multi-body bonded interactions or for bBondComm.
8103 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8108 for (zone = 0; zone < nzone_send; zone++)
8110 if (tric_dist[dim_ind] && dim_ind > 0)
8112 /* Determine slightly more optimized skew_fac's
8114 * This reduces the number of communicated atoms
8115 * by about 10% for 3D DD of rhombic dodecahedra.
8117 for (dimd = 0; dimd < dim; dimd++)
8119 sf2_round[dimd] = 1;
8120 if (ddbox->tric_dir[dimd])
8122 for (i = dd->dim[dimd]+1; i < DIM; i++)
8124 /* If we are shifted in dimension i
8125 * and the cell plane is tilted forward
8126 * in dimension i, skip this coupling.
8128 if (!(zones->shift[nzone+zone][i] &&
8129 ddbox->v[dimd][i][dimd] >= 0))
8132 sqr(ddbox->v[dimd][i][dimd]);
8135 sf2_round[dimd] = 1/sf2_round[dimd];
8140 zonei = zone_perm[dim_ind][zone];
8143 /* Here we permutate the zones to obtain a convenient order
8144 * for neighbor searching
8146 cg0 = zone_cg_range[zonei];
8147 cg1 = zone_cg_range[zonei+1];
8151 /* Look only at the cg's received in the previous grid pulse
8153 cg1 = zone_cg_range[nzone+zone+1];
8154 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8157 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8158 for (th = 0; th < comm->nth; th++)
8160 gmx_domdec_ind_t *ind_p;
8161 int **ibuf_p, *ibuf_nalloc_p;
8163 int *nsend_p, *nat_p;
8169 /* Thread 0 writes in the comm buffers */
8171 ibuf_p = &comm->buf_int;
8172 ibuf_nalloc_p = &comm->nalloc_int;
8173 vbuf_p = &comm->vbuf;
8176 nsend_zone_p = &ind->nsend[zone];
8180 /* Other threads write into temp buffers */
8181 ind_p = &comm->dth[th].ind;
8182 ibuf_p = &comm->dth[th].ibuf;
8183 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8184 vbuf_p = &comm->dth[th].vbuf;
8185 nsend_p = &comm->dth[th].nsend;
8186 nat_p = &comm->dth[th].nat;
8187 nsend_zone_p = &comm->dth[th].nsend_zone;
8189 comm->dth[th].nsend = 0;
8190 comm->dth[th].nat = 0;
8191 comm->dth[th].nsend_zone = 0;
8201 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8202 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8205 /* Get the cg's for this pulse in this zone */
8206 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8208 dim, dim_ind, dim0, dim1, dim2,
8211 normal, skew_fac2_d, skew_fac_01,
8212 v_d, v_0, v_1, &corners, sf2_round,
8213 bDistBonded, bBondComm,
8217 ibuf_p, ibuf_nalloc_p,
8223 /* Append data of threads>=1 to the communication buffers */
8224 for (th = 1; th < comm->nth; th++)
8226 dd_comm_setup_work_t *dth;
8229 dth = &comm->dth[th];
8231 ns1 = nsend + dth->nsend_zone;
8232 if (ns1 > ind->nalloc)
8234 ind->nalloc = over_alloc_dd(ns1);
8235 srenew(ind->index, ind->nalloc);
8237 if (ns1 > comm->nalloc_int)
8239 comm->nalloc_int = over_alloc_dd(ns1);
8240 srenew(comm->buf_int, comm->nalloc_int);
8242 if (ns1 > comm->vbuf.nalloc)
8244 comm->vbuf.nalloc = over_alloc_dd(ns1);
8245 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8248 for (i = 0; i < dth->nsend_zone; i++)
8250 ind->index[nsend] = dth->ind.index[i];
8251 comm->buf_int[nsend] = dth->ibuf[i];
8252 copy_rvec(dth->vbuf.v[i],
8253 comm->vbuf.v[nsend]);
8257 ind->nsend[zone] += dth->nsend_zone;
8260 /* Clear the counts in case we do not have pbc */
8261 for (zone = nzone_send; zone < nzone; zone++)
8263 ind->nsend[zone] = 0;
8265 ind->nsend[nzone] = nsend;
8266 ind->nsend[nzone+1] = nat;
8267 /* Communicate the number of cg's and atoms to receive */
8268 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8269 ind->nsend, nzone+2,
8270 ind->nrecv, nzone+2);
8272 /* The rvec buffer is also required for atom buffers of size nsend
8273 * in dd_move_x and dd_move_f.
8275 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8279 /* We can receive in place if only the last zone is not empty */
8280 for (zone = 0; zone < nzone-1; zone++)
8282 if (ind->nrecv[zone] > 0)
8284 cd->bInPlace = FALSE;
8289 /* The int buffer is only required here for the cg indices */
8290 if (ind->nrecv[nzone] > comm->nalloc_int2)
8292 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8293 srenew(comm->buf_int2, comm->nalloc_int2);
8295 /* The rvec buffer is also required for atom buffers
8296 * of size nrecv in dd_move_x and dd_move_f.
8298 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8299 vec_rvec_check_alloc(&comm->vbuf2, i);
8303 /* Make space for the global cg indices */
8304 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8305 || dd->cg_nalloc == 0)
8307 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8308 srenew(index_gl, dd->cg_nalloc);
8309 srenew(cgindex, dd->cg_nalloc+1);
8311 /* Communicate the global cg indices */
8314 recv_i = index_gl + pos_cg;
8318 recv_i = comm->buf_int2;
8320 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8321 comm->buf_int, nsend,
8322 recv_i, ind->nrecv[nzone]);
8324 /* Make space for cg_cm */
8325 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8326 if (fr->cutoff_scheme == ecutsGROUP)
8334 /* Communicate cg_cm */
8337 recv_vr = cg_cm + pos_cg;
8341 recv_vr = comm->vbuf2.v;
8343 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8344 comm->vbuf.v, nsend,
8345 recv_vr, ind->nrecv[nzone]);
8347 /* Make the charge group index */
8350 zone = (p == 0 ? 0 : nzone - 1);
8351 while (zone < nzone)
8353 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8355 cg_gl = index_gl[pos_cg];
8356 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8357 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8358 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8361 /* Update the charge group presence,
8362 * so we can use it in the next pass of the loop.
8364 comm->bLocalCG[cg_gl] = TRUE;
8370 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8373 zone_cg_range[nzone+zone] = pos_cg;
8378 /* This part of the code is never executed with bBondComm. */
8379 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8380 index_gl, recv_i, cg_cm, recv_vr,
8381 cgindex, fr->cginfo_mb, fr->cginfo);
8382 pos_cg += ind->nrecv[nzone];
8384 nat_tot += ind->nrecv[nzone+1];
8388 /* Store the atom block for easy copying of communication buffers */
8389 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8393 dd->index_gl = index_gl;
8394 dd->cgindex = cgindex;
8396 dd->ncg_tot = zone_cg_range[zones->n];
8397 dd->nat_tot = nat_tot;
8398 comm->nat[ddnatHOME] = dd->nat_home;
8399 for (i = ddnatZONE; i < ddnatNR; i++)
8401 comm->nat[i] = dd->nat_tot;
8406 /* We don't need to update cginfo, since that was alrady done above.
8407 * So we pass NULL for the forcerec.
8409 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8410 NULL, comm->bLocalCG);
8415 fprintf(debug, "Finished setting up DD communication, zones:");
8416 for (c = 0; c < zones->n; c++)
8418 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8420 fprintf(debug, "\n");
8424 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8428 for (c = 0; c < zones->nizone; c++)
8430 zones->izone[c].cg1 = zones->cg_range[c+1];
8431 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8432 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8436 static void set_zones_size(gmx_domdec_t *dd,
8437 matrix box, const gmx_ddbox_t *ddbox,
8438 int zone_start, int zone_end)
8440 gmx_domdec_comm_t *comm;
8441 gmx_domdec_zones_t *zones;
8443 int z, zi, zj0, zj1, d, dim;
8446 real size_j, add_tric;
8451 zones = &comm->zones;
8453 /* Do we need to determine extra distances for multi-body bondeds? */
8454 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8456 for (z = zone_start; z < zone_end; z++)
8458 /* Copy cell limits to zone limits.
8459 * Valid for non-DD dims and non-shifted dims.
8461 copy_rvec(comm->cell_x0, zones->size[z].x0);
8462 copy_rvec(comm->cell_x1, zones->size[z].x1);
8465 for (d = 0; d < dd->ndim; d++)
8469 for (z = 0; z < zones->n; z++)
8471 /* With a staggered grid we have different sizes
8472 * for non-shifted dimensions.
8474 if (dd->bGridJump && zones->shift[z][dim] == 0)
8478 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8479 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8483 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8484 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8490 rcmbs = comm->cutoff_mbody;
8491 if (ddbox->tric_dir[dim])
8493 rcs /= ddbox->skew_fac[dim];
8494 rcmbs /= ddbox->skew_fac[dim];
8497 /* Set the lower limit for the shifted zone dimensions */
8498 for (z = zone_start; z < zone_end; z++)
8500 if (zones->shift[z][dim] > 0)
8503 if (!dd->bGridJump || d == 0)
8505 zones->size[z].x0[dim] = comm->cell_x1[dim];
8506 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8510 /* Here we take the lower limit of the zone from
8511 * the lowest domain of the zone below.
8515 zones->size[z].x0[dim] =
8516 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8522 zones->size[z].x0[dim] =
8523 zones->size[zone_perm[2][z-4]].x0[dim];
8527 zones->size[z].x0[dim] =
8528 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8531 /* A temporary limit, is updated below */
8532 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8536 for (zi = 0; zi < zones->nizone; zi++)
8538 if (zones->shift[zi][dim] == 0)
8540 /* This takes the whole zone into account.
8541 * With multiple pulses this will lead
8542 * to a larger zone then strictly necessary.
8544 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8545 zones->size[zi].x1[dim]+rcmbs);
8553 /* Loop over the i-zones to set the upper limit of each
8556 for (zi = 0; zi < zones->nizone; zi++)
8558 if (zones->shift[zi][dim] == 0)
8560 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8562 if (zones->shift[z][dim] > 0)
8564 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8565 zones->size[zi].x1[dim]+rcs);
8572 for (z = zone_start; z < zone_end; z++)
8574 /* Initialization only required to keep the compiler happy */
8575 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8578 /* To determine the bounding box for a zone we need to find
8579 * the extreme corners of 4, 2 or 1 corners.
8581 nc = 1 << (ddbox->npbcdim - 1);
8583 for (c = 0; c < nc; c++)
8585 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8589 corner[YY] = zones->size[z].x0[YY];
8593 corner[YY] = zones->size[z].x1[YY];
8597 corner[ZZ] = zones->size[z].x0[ZZ];
8601 corner[ZZ] = zones->size[z].x1[ZZ];
8603 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8605 /* With 1D domain decomposition the cg's are not in
8606 * the triclinic box, but triclinic x-y and rectangular y-z.
8607 * Shift y back, so it will later end up at 0.
8609 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8611 /* Apply the triclinic couplings */
8612 for (i = YY; i < ddbox->npbcdim; i++)
8614 for (j = XX; j < i; j++)
8616 corner[j] += corner[i]*box[i][j]/box[i][i];
8621 copy_rvec(corner, corner_min);
8622 copy_rvec(corner, corner_max);
8626 for (i = 0; i < DIM; i++)
8628 corner_min[i] = min(corner_min[i], corner[i]);
8629 corner_max[i] = max(corner_max[i], corner[i]);
8633 /* Copy the extreme cornes without offset along x */
8634 for (i = 0; i < DIM; i++)
8636 zones->size[z].bb_x0[i] = corner_min[i];
8637 zones->size[z].bb_x1[i] = corner_max[i];
8639 /* Add the offset along x */
8640 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8641 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8644 if (zone_start == 0)
8647 for (dim = 0; dim < DIM; dim++)
8649 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8651 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8656 for (z = zone_start; z < zone_end; z++)
8658 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8660 zones->size[z].x0[XX], zones->size[z].x1[XX],
8661 zones->size[z].x0[YY], zones->size[z].x1[YY],
8662 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8663 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8665 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8666 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8667 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8672 static int comp_cgsort(const void *a, const void *b)
8676 gmx_cgsort_t *cga, *cgb;
8677 cga = (gmx_cgsort_t *)a;
8678 cgb = (gmx_cgsort_t *)b;
8680 comp = cga->nsc - cgb->nsc;
8683 comp = cga->ind_gl - cgb->ind_gl;
8689 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8694 /* Order the data */
8695 for (i = 0; i < n; i++)
8697 buf[i] = a[sort[i].ind];
8700 /* Copy back to the original array */
8701 for (i = 0; i < n; i++)
8707 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8712 /* Order the data */
8713 for (i = 0; i < n; i++)
8715 copy_rvec(v[sort[i].ind], buf[i]);
8718 /* Copy back to the original array */
8719 for (i = 0; i < n; i++)
8721 copy_rvec(buf[i], v[i]);
8725 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8728 int a, atot, cg, cg0, cg1, i;
8730 if (cgindex == NULL)
8732 /* Avoid the useless loop of the atoms within a cg */
8733 order_vec_cg(ncg, sort, v, buf);
8738 /* Order the data */
8740 for (cg = 0; cg < ncg; cg++)
8742 cg0 = cgindex[sort[cg].ind];
8743 cg1 = cgindex[sort[cg].ind+1];
8744 for (i = cg0; i < cg1; i++)
8746 copy_rvec(v[i], buf[a]);
8752 /* Copy back to the original array */
8753 for (a = 0; a < atot; a++)
8755 copy_rvec(buf[a], v[a]);
8759 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8760 int nsort_new, gmx_cgsort_t *sort_new,
8761 gmx_cgsort_t *sort1)
8765 /* The new indices are not very ordered, so we qsort them */
8766 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8768 /* sort2 is already ordered, so now we can merge the two arrays */
8772 while (i2 < nsort2 || i_new < nsort_new)
8776 sort1[i1++] = sort_new[i_new++];
8778 else if (i_new == nsort_new)
8780 sort1[i1++] = sort2[i2++];
8782 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8783 (sort2[i2].nsc == sort_new[i_new].nsc &&
8784 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8786 sort1[i1++] = sort2[i2++];
8790 sort1[i1++] = sort_new[i_new++];
8795 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8797 gmx_domdec_sort_t *sort;
8798 gmx_cgsort_t *cgsort, *sort_i;
8799 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8800 int sort_last, sort_skip;
8802 sort = dd->comm->sort;
8804 a = fr->ns.grid->cell_index;
8806 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8808 if (ncg_home_old >= 0)
8810 /* The charge groups that remained in the same ns grid cell
8811 * are completely ordered. So we can sort efficiently by sorting
8812 * the charge groups that did move into the stationary list.
8817 for (i = 0; i < dd->ncg_home; i++)
8819 /* Check if this cg did not move to another node */
8822 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8824 /* This cg is new on this node or moved ns grid cell */
8825 if (nsort_new >= sort->sort_new_nalloc)
8827 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8828 srenew(sort->sort_new, sort->sort_new_nalloc);
8830 sort_i = &(sort->sort_new[nsort_new++]);
8834 /* This cg did not move */
8835 sort_i = &(sort->sort2[nsort2++]);
8837 /* Sort on the ns grid cell indices
8838 * and the global topology index.
8839 * index_gl is irrelevant with cell ns,
8840 * but we set it here anyhow to avoid a conditional.
8843 sort_i->ind_gl = dd->index_gl[i];
8850 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8853 /* Sort efficiently */
8854 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8859 cgsort = sort->sort;
8861 for (i = 0; i < dd->ncg_home; i++)
8863 /* Sort on the ns grid cell indices
8864 * and the global topology index
8866 cgsort[i].nsc = a[i];
8867 cgsort[i].ind_gl = dd->index_gl[i];
8869 if (cgsort[i].nsc < moved)
8876 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8878 /* Determine the order of the charge groups using qsort */
8879 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8885 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8888 int ncg_new, i, *a, na;
8890 sort = dd->comm->sort->sort;
8892 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8895 for (i = 0; i < na; i++)
8899 sort[ncg_new].ind = a[i];
8907 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8908 rvec *cgcm, t_forcerec *fr, t_state *state,
8911 gmx_domdec_sort_t *sort;
8912 gmx_cgsort_t *cgsort, *sort_i;
8914 int ncg_new, i, *ibuf, cgsize;
8917 sort = dd->comm->sort;
8919 if (dd->ncg_home > sort->sort_nalloc)
8921 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8922 srenew(sort->sort, sort->sort_nalloc);
8923 srenew(sort->sort2, sort->sort_nalloc);
8925 cgsort = sort->sort;
8927 switch (fr->cutoff_scheme)
8930 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8933 ncg_new = dd_sort_order_nbnxn(dd, fr);
8936 gmx_incons("unimplemented");
8940 /* We alloc with the old size, since cgindex is still old */
8941 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8942 vbuf = dd->comm->vbuf.v;
8946 cgindex = dd->cgindex;
8953 /* Remove the charge groups which are no longer at home here */
8954 dd->ncg_home = ncg_new;
8957 fprintf(debug, "Set the new home charge group count to %d\n",
8961 /* Reorder the state */
8962 for (i = 0; i < estNR; i++)
8964 if (EST_DISTR(i) && (state->flags & (1<<i)))
8969 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8972 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8975 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8978 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8982 case estDISRE_INITF:
8983 case estDISRE_RM3TAV:
8984 case estORIRE_INITF:
8986 /* No ordering required */
8989 gmx_incons("Unknown state entry encountered in dd_sort_state");
8994 if (fr->cutoff_scheme == ecutsGROUP)
8997 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9000 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9002 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9003 srenew(sort->ibuf, sort->ibuf_nalloc);
9006 /* Reorder the global cg index */
9007 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9008 /* Reorder the cginfo */
9009 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9010 /* Rebuild the local cg index */
9014 for (i = 0; i < dd->ncg_home; i++)
9016 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9017 ibuf[i+1] = ibuf[i] + cgsize;
9019 for (i = 0; i < dd->ncg_home+1; i++)
9021 dd->cgindex[i] = ibuf[i];
9026 for (i = 0; i < dd->ncg_home+1; i++)
9031 /* Set the home atom number */
9032 dd->nat_home = dd->cgindex[dd->ncg_home];
9034 if (fr->cutoff_scheme == ecutsVERLET)
9036 /* The atoms are now exactly in grid order, update the grid order */
9037 nbnxn_set_atomorder(fr->nbv->nbs);
9041 /* Copy the sorted ns cell indices back to the ns grid struct */
9042 for (i = 0; i < dd->ncg_home; i++)
9044 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9046 fr->ns.grid->nr = dd->ncg_home;
9050 static void add_dd_statistics(gmx_domdec_t *dd)
9052 gmx_domdec_comm_t *comm;
9057 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9059 comm->sum_nat[ddnat-ddnatZONE] +=
9060 comm->nat[ddnat] - comm->nat[ddnat-1];
9065 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9067 gmx_domdec_comm_t *comm;
9072 /* Reset all the statistics and counters for total run counting */
9073 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9075 comm->sum_nat[ddnat-ddnatZONE] = 0;
9079 comm->load_step = 0;
9082 clear_ivec(comm->load_lim);
9087 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9089 gmx_domdec_comm_t *comm;
9093 comm = cr->dd->comm;
9095 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9102 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");
9104 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9106 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9111 " av. #atoms communicated per step for force: %d x %.1f\n",
9115 if (cr->dd->vsite_comm)
9118 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9119 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9124 if (cr->dd->constraint_comm)
9127 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9128 1 + ir->nLincsIter, av);
9132 gmx_incons(" Unknown type for DD statistics");
9135 fprintf(fplog, "\n");
9137 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9139 print_dd_load_av(fplog, cr->dd);
9143 void dd_partition_system(FILE *fplog,
9144 gmx_large_int_t step,
9146 gmx_bool bMasterState,
9148 t_state *state_global,
9149 gmx_mtop_t *top_global,
9151 t_state *state_local,
9154 gmx_localtop_t *top_local,
9157 gmx_shellfc_t shellfc,
9158 gmx_constr_t constr,
9160 gmx_wallcycle_t wcycle,
9164 gmx_domdec_comm_t *comm;
9165 gmx_ddbox_t ddbox = {0};
9167 gmx_large_int_t step_pcoupl;
9168 rvec cell_ns_x0, cell_ns_x1;
9169 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9170 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9171 gmx_bool bRedist, bSortCG, bResortAll;
9172 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9179 bBoxChanged = (bMasterState || DEFORM(*ir));
9180 if (ir->epc != epcNO)
9182 /* With nstpcouple > 1 pressure coupling happens.
9183 * one step after calculating the pressure.
9184 * Box scaling happens at the end of the MD step,
9185 * after the DD partitioning.
9186 * We therefore have to do DLB in the first partitioning
9187 * after an MD step where P-coupling occured.
9188 * We need to determine the last step in which p-coupling occurred.
9189 * MRS -- need to validate this for vv?
9194 step_pcoupl = step - 1;
9198 step_pcoupl = ((step - 1)/n)*n + 1;
9200 if (step_pcoupl >= comm->partition_step)
9206 bNStGlobalComm = (step % nstglobalcomm == 0);
9208 if (!comm->bDynLoadBal)
9214 /* Should we do dynamic load balacing this step?
9215 * Since it requires (possibly expensive) global communication,
9216 * we might want to do DLB less frequently.
9218 if (bBoxChanged || ir->epc != epcNO)
9220 bDoDLB = bBoxChanged;
9224 bDoDLB = bNStGlobalComm;
9228 /* Check if we have recorded loads on the nodes */
9229 if (comm->bRecordLoad && dd_load_count(comm))
9231 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9233 /* Check if we should use DLB at the second partitioning
9234 * and every 100 partitionings,
9235 * so the extra communication cost is negligible.
9237 n = max(100, nstglobalcomm);
9238 bCheckDLB = (comm->n_load_collect == 0 ||
9239 comm->n_load_have % n == n-1);
9246 /* Print load every nstlog, first and last step to the log file */
9247 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9248 comm->n_load_collect == 0 ||
9250 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9252 /* Avoid extra communication due to verbose screen output
9253 * when nstglobalcomm is set.
9255 if (bDoDLB || bLogLoad || bCheckDLB ||
9256 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9258 get_load_distribution(dd, wcycle);
9263 dd_print_load(fplog, dd, step-1);
9267 dd_print_load_verbose(dd);
9270 comm->n_load_collect++;
9274 /* Since the timings are node dependent, the master decides */
9278 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9281 fprintf(debug, "step %s, imb loss %f\n",
9282 gmx_step_str(step, sbuf),
9283 dd_force_imb_perf_loss(dd));
9286 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9289 turn_on_dlb(fplog, cr, step);
9294 comm->n_load_have++;
9297 cgs_gl = &comm->cgs_gl;
9302 /* Clear the old state */
9303 clear_dd_indices(dd, 0, 0);
9305 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9306 TRUE, cgs_gl, state_global->x, &ddbox);
9308 get_cg_distribution(fplog, step, dd, cgs_gl,
9309 state_global->box, &ddbox, state_global->x);
9311 dd_distribute_state(dd, cgs_gl,
9312 state_global, state_local, f);
9314 dd_make_local_cgs(dd, &top_local->cgs);
9316 /* Ensure that we have space for the new distribution */
9317 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9319 if (fr->cutoff_scheme == ecutsGROUP)
9321 calc_cgcm(fplog, 0, dd->ncg_home,
9322 &top_local->cgs, state_local->x, fr->cg_cm);
9325 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9327 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9331 else if (state_local->ddp_count != dd->ddp_count)
9333 if (state_local->ddp_count > dd->ddp_count)
9335 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9338 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9340 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);
9343 /* Clear the old state */
9344 clear_dd_indices(dd, 0, 0);
9346 /* Build the new indices */
9347 rebuild_cgindex(dd, cgs_gl->index, state_local);
9348 make_dd_indices(dd, cgs_gl->index, 0);
9350 if (fr->cutoff_scheme == ecutsGROUP)
9352 /* Redetermine the cg COMs */
9353 calc_cgcm(fplog, 0, dd->ncg_home,
9354 &top_local->cgs, state_local->x, fr->cg_cm);
9357 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9359 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9361 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9362 TRUE, &top_local->cgs, state_local->x, &ddbox);
9364 bRedist = comm->bDynLoadBal;
9368 /* We have the full state, only redistribute the cgs */
9370 /* Clear the non-home indices */
9371 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9373 /* Avoid global communication for dim's without pbc and -gcom */
9374 if (!bNStGlobalComm)
9376 copy_rvec(comm->box0, ddbox.box0 );
9377 copy_rvec(comm->box_size, ddbox.box_size);
9379 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9380 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9385 /* For dim's without pbc and -gcom */
9386 copy_rvec(ddbox.box0, comm->box0 );
9387 copy_rvec(ddbox.box_size, comm->box_size);
9389 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9392 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9394 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9397 /* Check if we should sort the charge groups */
9398 if (comm->nstSortCG > 0)
9400 bSortCG = (bMasterState ||
9401 (bRedist && (step % comm->nstSortCG == 0)));
9408 ncg_home_old = dd->ncg_home;
9413 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9415 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9416 state_local, f, fr, mdatoms,
9417 !bSortCG, nrnb, &cg0, &ncg_moved);
9419 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9422 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9424 &comm->cell_x0, &comm->cell_x1,
9425 dd->ncg_home, fr->cg_cm,
9426 cell_ns_x0, cell_ns_x1, &grid_density);
9430 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9433 switch (fr->cutoff_scheme)
9436 copy_ivec(fr->ns.grid->n, ncells_old);
9437 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9438 state_local->box, cell_ns_x0, cell_ns_x1,
9439 fr->rlistlong, grid_density);
9442 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9445 gmx_incons("unimplemented");
9447 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9448 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9452 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9454 /* Sort the state on charge group position.
9455 * This enables exact restarts from this step.
9456 * It also improves performance by about 15% with larger numbers
9457 * of atoms per node.
9460 /* Fill the ns grid with the home cell,
9461 * so we can sort with the indices.
9463 set_zones_ncg_home(dd);
9465 switch (fr->cutoff_scheme)
9468 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9470 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9472 comm->zones.size[0].bb_x0,
9473 comm->zones.size[0].bb_x1,
9475 comm->zones.dens_zone0,
9478 ncg_moved, bRedist ? comm->moved : NULL,
9479 fr->nbv->grp[eintLocal].kernel_type,
9480 fr->nbv->grp[eintLocal].nbat);
9482 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9485 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9486 0, dd->ncg_home, fr->cg_cm);
9488 copy_ivec(fr->ns.grid->n, ncells_new);
9491 gmx_incons("unimplemented");
9494 bResortAll = bMasterState;
9496 /* Check if we can user the old order and ns grid cell indices
9497 * of the charge groups to sort the charge groups efficiently.
9499 if (ncells_new[XX] != ncells_old[XX] ||
9500 ncells_new[YY] != ncells_old[YY] ||
9501 ncells_new[ZZ] != ncells_old[ZZ])
9508 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9509 gmx_step_str(step, sbuf), dd->ncg_home);
9511 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9512 bResortAll ? -1 : ncg_home_old);
9513 /* Rebuild all the indices */
9515 ga2la_clear(dd->ga2la);
9517 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9520 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9522 /* Setup up the communication and communicate the coordinates */
9523 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9525 /* Set the indices */
9526 make_dd_indices(dd, cgs_gl->index, cg0);
9528 /* Set the charge group boundaries for neighbor searching */
9529 set_cg_boundaries(&comm->zones);
9531 if (fr->cutoff_scheme == ecutsVERLET)
9533 set_zones_size(dd, state_local->box, &ddbox,
9534 bSortCG ? 1 : 0, comm->zones.n);
9537 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9540 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9541 -1,state_local->x,state_local->box);
9544 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9546 /* Extract a local topology from the global topology */
9547 for (i = 0; i < dd->ndim; i++)
9549 np[dd->dim[i]] = comm->cd[i].np;
9551 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9552 comm->cellsize_min, np,
9554 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9555 vsite, top_global, top_local);
9557 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9559 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9561 /* Set up the special atom communication */
9562 n = comm->nat[ddnatZONE];
9563 for (i = ddnatZONE+1; i < ddnatNR; i++)
9568 if (vsite && vsite->n_intercg_vsite)
9570 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9574 if (dd->bInterCGcons || dd->bInterCGsettles)
9576 /* Only for inter-cg constraints we need special code */
9577 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9578 constr, ir->nProjOrder,
9579 top_local->idef.il);
9583 gmx_incons("Unknown special atom type setup");
9588 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9590 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9592 /* Make space for the extra coordinates for virtual site
9593 * or constraint communication.
9595 state_local->natoms = comm->nat[ddnatNR-1];
9596 if (state_local->natoms > state_local->nalloc)
9598 dd_realloc_state(state_local, f, state_local->natoms);
9601 if (fr->bF_NoVirSum)
9603 if (vsite && vsite->n_intercg_vsite)
9605 nat_f_novirsum = comm->nat[ddnatVSITE];
9609 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9611 nat_f_novirsum = dd->nat_tot;
9615 nat_f_novirsum = dd->nat_home;
9624 /* Set the number of atoms required for the force calculation.
9625 * Forces need to be constrained when using a twin-range setup
9626 * or with energy minimization. For simple simulations we could
9627 * avoid some allocation, zeroing and copying, but this is
9628 * probably not worth the complications ande checking.
9630 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9631 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9633 /* We make the all mdatoms up to nat_tot_con.
9634 * We could save some work by only setting invmass
9635 * between nat_tot and nat_tot_con.
9637 /* This call also sets the new number of home particles to dd->nat_home */
9638 atoms2md(top_global, ir,
9639 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9641 /* Now we have the charges we can sort the FE interactions */
9642 dd_sort_local_top(dd, mdatoms, top_local);
9646 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9647 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9652 /* Make the local shell stuff, currently no communication is done */
9653 make_local_shells(cr, mdatoms, shellfc);
9656 if (ir->implicit_solvent)
9658 make_local_gb(cr, fr->born, ir->gb_algorithm);
9661 init_bonded_thread_force_reduction(fr, &top_local->idef);
9663 if (!(cr->duty & DUTY_PME))
9665 /* Send the charges to our PME only node */
9666 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9667 mdatoms->chargeA, mdatoms->chargeB,
9668 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9673 set_constraints(constr, top_local, ir, mdatoms, cr);
9676 if (ir->ePull != epullNO)
9678 /* Update the local pull groups */
9679 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9684 /* Update the local rotation groups */
9685 dd_make_local_rotation_groups(dd, ir->rot);
9689 add_dd_statistics(dd);
9691 /* Make sure we only count the cycles for this DD partitioning */
9692 clear_dd_cycle_counts(dd);
9694 /* Because the order of the atoms might have changed since
9695 * the last vsite construction, we need to communicate the constructing
9696 * atom coordinates again (for spreading the forces this MD step).
9698 dd_move_x_vsites(dd, state_local->box, state_local->x);
9700 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9702 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9704 dd_move_x(dd, state_local->box, state_local->x);
9705 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9706 -1, state_local->x, state_local->box);
9709 /* Store the partitioning step */
9710 comm->partition_step = step;
9712 /* Increase the DD partitioning counter */
9714 /* The state currently matches this DD partitioning count, store it */
9715 state_local->ddp_count = dd->ddp_count;
9718 /* The DD master node knows the complete cg distribution,
9719 * store the count so we can possibly skip the cg info communication.
9721 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9724 if (comm->DD_debug > 0)
9726 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9727 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9728 "after partitioning");