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 void get_pme_nnodes(const gmx_domdec_t *dd,
2321 int *npmenodes_x, int *npmenodes_y)
2325 *npmenodes_x = dd->comm->npmenodes_x;
2326 *npmenodes_y = dd->comm->npmenodes_y;
2335 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2337 gmx_bool bPMEOnlyNode;
2339 if (DOMAINDECOMP(cr))
2341 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2345 bPMEOnlyNode = FALSE;
2348 return bPMEOnlyNode;
2351 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2352 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2356 ivec coord, coord_pme;
2360 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2363 for (x = 0; x < dd->nc[XX]; x++)
2365 for (y = 0; y < dd->nc[YY]; y++)
2367 for (z = 0; z < dd->nc[ZZ]; z++)
2369 if (dd->comm->bCartesianPP_PME)
2374 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2375 if (dd->ci[XX] == coord_pme[XX] &&
2376 dd->ci[YY] == coord_pme[YY] &&
2377 dd->ci[ZZ] == coord_pme[ZZ])
2379 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2384 /* The slab corresponds to the nodeid in the PME group */
2385 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2387 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2394 /* The last PP-only node is the peer node */
2395 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2399 fprintf(debug, "Receive coordinates from PP nodes:");
2400 for (x = 0; x < *nmy_ddnodes; x++)
2402 fprintf(debug, " %d", (*my_ddnodes)[x]);
2404 fprintf(debug, "\n");
2408 static gmx_bool receive_vir_ener(t_commrec *cr)
2410 gmx_domdec_comm_t *comm;
2411 int pmenode, coords[DIM], rank;
2415 if (cr->npmenodes < cr->dd->nnodes)
2417 comm = cr->dd->comm;
2418 if (comm->bCartesianPP_PME)
2420 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2422 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2423 coords[comm->cartpmedim]++;
2424 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2426 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2427 if (dd_simnode2pmenode(cr, rank) == pmenode)
2429 /* This is not the last PP node for pmenode */
2437 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2438 if (cr->sim_nodeid+1 < cr->nnodes &&
2439 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2441 /* This is not the last PP node for pmenode */
2450 static void set_zones_ncg_home(gmx_domdec_t *dd)
2452 gmx_domdec_zones_t *zones;
2455 zones = &dd->comm->zones;
2457 zones->cg_range[0] = 0;
2458 for (i = 1; i < zones->n+1; i++)
2460 zones->cg_range[i] = dd->ncg_home;
2464 static void rebuild_cgindex(gmx_domdec_t *dd,
2465 const int *gcgs_index, t_state *state)
2467 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2470 dd_cg_gl = dd->index_gl;
2471 cgindex = dd->cgindex;
2474 for (i = 0; i < state->ncg_gl; i++)
2478 dd_cg_gl[i] = cg_gl;
2479 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2483 dd->ncg_home = state->ncg_gl;
2486 set_zones_ncg_home(dd);
2489 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2491 while (cg >= cginfo_mb->cg_end)
2496 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2499 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2500 t_forcerec *fr, char *bLocalCG)
2502 cginfo_mb_t *cginfo_mb;
2508 cginfo_mb = fr->cginfo_mb;
2509 cginfo = fr->cginfo;
2511 for (cg = cg0; cg < cg1; cg++)
2513 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2517 if (bLocalCG != NULL)
2519 for (cg = cg0; cg < cg1; cg++)
2521 bLocalCG[index_gl[cg]] = TRUE;
2526 static void make_dd_indices(gmx_domdec_t *dd,
2527 const int *gcgs_index, int cg_start)
2529 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2530 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2535 bLocalCG = dd->comm->bLocalCG;
2537 if (dd->nat_tot > dd->gatindex_nalloc)
2539 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2540 srenew(dd->gatindex, dd->gatindex_nalloc);
2543 nzone = dd->comm->zones.n;
2544 zone2cg = dd->comm->zones.cg_range;
2545 zone_ncg1 = dd->comm->zone_ncg1;
2546 index_gl = dd->index_gl;
2547 gatindex = dd->gatindex;
2548 bCGs = dd->comm->bCGs;
2550 if (zone2cg[1] != dd->ncg_home)
2552 gmx_incons("dd->ncg_zone is not up to date");
2555 /* Make the local to global and global to local atom index */
2556 a = dd->cgindex[cg_start];
2557 for (zone = 0; zone < nzone; zone++)
2565 cg0 = zone2cg[zone];
2567 cg1 = zone2cg[zone+1];
2568 cg1_p1 = cg0 + zone_ncg1[zone];
2570 for (cg = cg0; cg < cg1; cg++)
2575 /* Signal that this cg is from more than one pulse away */
2578 cg_gl = index_gl[cg];
2581 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2584 ga2la_set(dd->ga2la, a_gl, a, zone1);
2590 gatindex[a] = cg_gl;
2591 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2598 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2601 int ncg, i, ngl, nerr;
2604 if (bLocalCG == NULL)
2608 for (i = 0; i < dd->ncg_tot; i++)
2610 if (!bLocalCG[dd->index_gl[i]])
2613 "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);
2618 for (i = 0; i < ncg_sys; i++)
2625 if (ngl != dd->ncg_tot)
2627 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);
2634 static void check_index_consistency(gmx_domdec_t *dd,
2635 int natoms_sys, int ncg_sys,
2638 int nerr, ngl, i, a, cell;
2643 if (dd->comm->DD_debug > 1)
2645 snew(have, natoms_sys);
2646 for (a = 0; a < dd->nat_tot; a++)
2648 if (have[dd->gatindex[a]] > 0)
2650 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);
2654 have[dd->gatindex[a]] = a + 1;
2660 snew(have, dd->nat_tot);
2663 for (i = 0; i < natoms_sys; i++)
2665 if (ga2la_get(dd->ga2la, i, &a, &cell))
2667 if (a >= dd->nat_tot)
2669 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);
2675 if (dd->gatindex[a] != i)
2677 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);
2684 if (ngl != dd->nat_tot)
2687 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2688 dd->rank, where, ngl, dd->nat_tot);
2690 for (a = 0; a < dd->nat_tot; a++)
2695 "DD node %d, %s: local atom %d, global %d has no global index\n",
2696 dd->rank, where, a+1, dd->gatindex[a]+1);
2701 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2705 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2706 dd->rank, where, nerr);
2710 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2717 /* Clear the whole list without searching */
2718 ga2la_clear(dd->ga2la);
2722 for (i = a_start; i < dd->nat_tot; i++)
2724 ga2la_del(dd->ga2la, dd->gatindex[i]);
2728 bLocalCG = dd->comm->bLocalCG;
2731 for (i = cg_start; i < dd->ncg_tot; i++)
2733 bLocalCG[dd->index_gl[i]] = FALSE;
2737 dd_clear_local_vsite_indices(dd);
2739 if (dd->constraints)
2741 dd_clear_local_constraint_indices(dd);
2745 /* This function should be used for moving the domain boudaries during DLB,
2746 * for obtaining the minimum cell size. It checks the initially set limit
2747 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2748 * and, possibly, a longer cut-off limit set for PME load balancing.
2750 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2754 cellsize_min = comm->cellsize_min[dim];
2756 if (!comm->bVacDLBNoLimit)
2758 /* The cut-off might have changed, e.g. by PME load balacning,
2759 * from the value used to set comm->cellsize_min, so check it.
2761 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2763 if (comm->bPMELoadBalDLBLimits)
2765 /* Check for the cut-off limit set by the PME load balancing */
2766 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2770 return cellsize_min;
2773 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2776 real grid_jump_limit;
2778 /* The distance between the boundaries of cells at distance
2779 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2780 * and by the fact that cells should not be shifted by more than
2781 * half their size, such that cg's only shift by one cell
2782 * at redecomposition.
2784 grid_jump_limit = comm->cellsize_limit;
2785 if (!comm->bVacDLBNoLimit)
2787 if (comm->bPMELoadBalDLBLimits)
2789 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2791 grid_jump_limit = max(grid_jump_limit,
2792 cutoff/comm->cd[dim_ind].np);
2795 return grid_jump_limit;
2798 static gmx_bool check_grid_jump(gmx_large_int_t step,
2804 gmx_domdec_comm_t *comm;
2813 for (d = 1; d < dd->ndim; d++)
2816 limit = grid_jump_limit(comm, cutoff, d);
2817 bfac = ddbox->box_size[dim];
2818 if (ddbox->tric_dir[dim])
2820 bfac *= ddbox->skew_fac[dim];
2822 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2823 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2831 /* This error should never be triggered under normal
2832 * circumstances, but you never know ...
2834 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.",
2835 gmx_step_str(step, buf),
2836 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2844 static int dd_load_count(gmx_domdec_comm_t *comm)
2846 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2849 static float dd_force_load(gmx_domdec_comm_t *comm)
2856 if (comm->eFlop > 1)
2858 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2863 load = comm->cycl[ddCyclF];
2864 if (comm->cycl_n[ddCyclF] > 1)
2866 /* Subtract the maximum of the last n cycle counts
2867 * to get rid of possible high counts due to other soures,
2868 * for instance system activity, that would otherwise
2869 * affect the dynamic load balancing.
2871 load -= comm->cycl_max[ddCyclF];
2878 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2880 gmx_domdec_comm_t *comm;
2885 snew(*dim_f, dd->nc[dim]+1);
2887 for (i = 1; i < dd->nc[dim]; i++)
2889 if (comm->slb_frac[dim])
2891 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2895 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2898 (*dim_f)[dd->nc[dim]] = 1;
2901 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2903 int pmeindex, slab, nso, i;
2906 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2912 ddpme->dim = dimind;
2914 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2916 ddpme->nslab = (ddpme->dim == 0 ?
2917 dd->comm->npmenodes_x :
2918 dd->comm->npmenodes_y);
2920 if (ddpme->nslab <= 1)
2925 nso = dd->comm->npmenodes/ddpme->nslab;
2926 /* Determine for each PME slab the PP location range for dimension dim */
2927 snew(ddpme->pp_min, ddpme->nslab);
2928 snew(ddpme->pp_max, ddpme->nslab);
2929 for (slab = 0; slab < ddpme->nslab; slab++)
2931 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2932 ddpme->pp_max[slab] = 0;
2934 for (i = 0; i < dd->nnodes; i++)
2936 ddindex2xyz(dd->nc, i, xyz);
2937 /* For y only use our y/z slab.
2938 * This assumes that the PME x grid size matches the DD grid size.
2940 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2942 pmeindex = ddindex2pmeindex(dd, i);
2945 slab = pmeindex/nso;
2949 slab = pmeindex % ddpme->nslab;
2951 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2952 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2956 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2959 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2961 if (dd->comm->ddpme[0].dim == XX)
2963 return dd->comm->ddpme[0].maxshift;
2971 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2973 if (dd->comm->ddpme[0].dim == YY)
2975 return dd->comm->ddpme[0].maxshift;
2977 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2979 return dd->comm->ddpme[1].maxshift;
2987 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2988 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2990 gmx_domdec_comm_t *comm;
2993 real range, pme_boundary;
2997 nc = dd->nc[ddpme->dim];
3000 if (!ddpme->dim_match)
3002 /* PP decomposition is not along dim: the worst situation */
3005 else if (ns <= 3 || (bUniform && ns == nc))
3007 /* The optimal situation */
3012 /* We need to check for all pme nodes which nodes they
3013 * could possibly need to communicate with.
3015 xmin = ddpme->pp_min;
3016 xmax = ddpme->pp_max;
3017 /* Allow for atoms to be maximally 2/3 times the cut-off
3018 * out of their DD cell. This is a reasonable balance between
3019 * between performance and support for most charge-group/cut-off
3022 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3023 /* Avoid extra communication when we are exactly at a boundary */
3027 for (s = 0; s < ns; s++)
3029 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3030 pme_boundary = (real)s/ns;
3033 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3035 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3039 pme_boundary = (real)(s+1)/ns;
3042 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3044 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3051 ddpme->maxshift = sh;
3055 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3056 ddpme->dim, ddpme->maxshift);
3060 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3064 for (d = 0; d < dd->ndim; d++)
3067 if (dim < ddbox->nboundeddim &&
3068 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3069 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3071 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",
3072 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3073 dd->nc[dim], dd->comm->cellsize_limit);
3078 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3079 gmx_bool bMaster, ivec npulse)
3081 gmx_domdec_comm_t *comm;
3084 real *cell_x, cell_dx, cellsize;
3088 for (d = 0; d < DIM; d++)
3090 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3092 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3095 cell_dx = ddbox->box_size[d]/dd->nc[d];
3098 for (j = 0; j < dd->nc[d]+1; j++)
3100 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3105 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3106 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3108 cellsize = cell_dx*ddbox->skew_fac[d];
3109 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3113 cellsize_min[d] = cellsize;
3117 /* Statically load balanced grid */
3118 /* Also when we are not doing a master distribution we determine
3119 * all cell borders in a loop to obtain identical values
3120 * to the master distribution case and to determine npulse.
3124 cell_x = dd->ma->cell_x[d];
3128 snew(cell_x, dd->nc[d]+1);
3130 cell_x[0] = ddbox->box0[d];
3131 for (j = 0; j < dd->nc[d]; j++)
3133 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3134 cell_x[j+1] = cell_x[j] + cell_dx;
3135 cellsize = cell_dx*ddbox->skew_fac[d];
3136 while (cellsize*npulse[d] < comm->cutoff &&
3137 npulse[d] < dd->nc[d]-1)
3141 cellsize_min[d] = min(cellsize_min[d], cellsize);
3145 comm->cell_x0[d] = cell_x[dd->ci[d]];
3146 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3150 /* The following limitation is to avoid that a cell would receive
3151 * some of its own home charge groups back over the periodic boundary.
3152 * Double charge groups cause trouble with the global indices.
3154 if (d < ddbox->npbcdim &&
3155 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3157 gmx_fatal_collective(FARGS, NULL, dd,
3158 "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",
3159 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3161 dd->nc[d], dd->nc[d],
3162 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3166 if (!comm->bDynLoadBal)
3168 copy_rvec(cellsize_min, comm->cellsize_min);
3171 for (d = 0; d < comm->npmedecompdim; d++)
3173 set_pme_maxshift(dd, &comm->ddpme[d],
3174 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3175 comm->ddpme[d].slb_dim_f);
3180 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3181 int d, int dim, gmx_domdec_root_t *root,
3183 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3185 gmx_domdec_comm_t *comm;
3186 int ncd, i, j, nmin, nmin_old;
3187 gmx_bool bLimLo, bLimHi;
3189 real fac, halfway, cellsize_limit_f_i, region_size;
3190 gmx_bool bPBC, bLastHi = FALSE;
3191 int nrange[] = {range[0], range[1]};
3193 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3199 bPBC = (dim < ddbox->npbcdim);
3201 cell_size = root->buf_ncd;
3205 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3208 /* First we need to check if the scaling does not make cells
3209 * smaller than the smallest allowed size.
3210 * We need to do this iteratively, since if a cell is too small,
3211 * it needs to be enlarged, which makes all the other cells smaller,
3212 * which could in turn make another cell smaller than allowed.
3214 for (i = range[0]; i < range[1]; i++)
3216 root->bCellMin[i] = FALSE;
3222 /* We need the total for normalization */
3224 for (i = range[0]; i < range[1]; i++)
3226 if (root->bCellMin[i] == FALSE)
3228 fac += cell_size[i];
3231 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3232 /* Determine the cell boundaries */
3233 for (i = range[0]; i < range[1]; i++)
3235 if (root->bCellMin[i] == FALSE)
3237 cell_size[i] *= fac;
3238 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3240 cellsize_limit_f_i = 0;
3244 cellsize_limit_f_i = cellsize_limit_f;
3246 if (cell_size[i] < cellsize_limit_f_i)
3248 root->bCellMin[i] = TRUE;
3249 cell_size[i] = cellsize_limit_f_i;
3253 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3256 while (nmin > nmin_old);
3259 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3260 /* For this check we should not use DD_CELL_MARGIN,
3261 * but a slightly smaller factor,
3262 * since rounding could get use below the limit.
3264 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3267 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",
3268 gmx_step_str(step, buf),
3269 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3270 ncd, comm->cellsize_min[dim]);
3273 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3277 /* Check if the boundary did not displace more than halfway
3278 * each of the cells it bounds, as this could cause problems,
3279 * especially when the differences between cell sizes are large.
3280 * If changes are applied, they will not make cells smaller
3281 * than the cut-off, as we check all the boundaries which
3282 * might be affected by a change and if the old state was ok,
3283 * the cells will at most be shrunk back to their old size.
3285 for (i = range[0]+1; i < range[1]; i++)
3287 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3288 if (root->cell_f[i] < halfway)
3290 root->cell_f[i] = halfway;
3291 /* Check if the change also causes shifts of the next boundaries */
3292 for (j = i+1; j < range[1]; j++)
3294 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3296 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3300 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3301 if (root->cell_f[i] > halfway)
3303 root->cell_f[i] = halfway;
3304 /* Check if the change also causes shifts of the next boundaries */
3305 for (j = i-1; j >= range[0]+1; j--)
3307 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3309 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3316 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3317 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3318 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3319 * for a and b nrange is used */
3322 /* Take care of the staggering of the cell boundaries */
3325 for (i = range[0]; i < range[1]; i++)
3327 root->cell_f_max0[i] = root->cell_f[i];
3328 root->cell_f_min1[i] = root->cell_f[i+1];
3333 for (i = range[0]+1; i < range[1]; i++)
3335 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3336 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3337 if (bLimLo && bLimHi)
3339 /* Both limits violated, try the best we can */
3340 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3341 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3342 nrange[0] = range[0];
3344 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3347 nrange[1] = range[1];
3348 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3354 /* root->cell_f[i] = root->bound_min[i]; */
3355 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3358 else if (bLimHi && !bLastHi)
3361 if (nrange[1] < range[1]) /* found a LimLo before */
3363 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3364 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3365 nrange[0] = nrange[1];
3367 root->cell_f[i] = root->bound_max[i];
3369 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3371 nrange[1] = range[1];
3374 if (nrange[1] < range[1]) /* found last a LimLo */
3376 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3377 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3378 nrange[0] = nrange[1];
3379 nrange[1] = range[1];
3380 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3382 else if (nrange[0] > range[0]) /* found at least one LimHi */
3384 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3391 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3392 int d, int dim, gmx_domdec_root_t *root,
3393 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3394 gmx_bool bUniform, gmx_large_int_t step)
3396 gmx_domdec_comm_t *comm;
3397 int ncd, d1, i, j, pos;
3399 real load_aver, load_i, imbalance, change, change_max, sc;
3400 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3404 int range[] = { 0, 0 };
3408 /* Convert the maximum change from the input percentage to a fraction */
3409 change_limit = comm->dlb_scale_lim*0.01;
3413 bPBC = (dim < ddbox->npbcdim);
3415 cell_size = root->buf_ncd;
3417 /* Store the original boundaries */
3418 for (i = 0; i < ncd+1; i++)
3420 root->old_cell_f[i] = root->cell_f[i];
3424 for (i = 0; i < ncd; i++)
3426 cell_size[i] = 1.0/ncd;
3429 else if (dd_load_count(comm))
3431 load_aver = comm->load[d].sum_m/ncd;
3433 for (i = 0; i < ncd; i++)
3435 /* Determine the relative imbalance of cell i */
3436 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3437 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3438 /* Determine the change of the cell size using underrelaxation */
3439 change = -relax*imbalance;
3440 change_max = max(change_max, max(change, -change));
3442 /* Limit the amount of scaling.
3443 * We need to use the same rescaling for all cells in one row,
3444 * otherwise the load balancing might not converge.
3447 if (change_max > change_limit)
3449 sc *= change_limit/change_max;
3451 for (i = 0; i < ncd; i++)
3453 /* Determine the relative imbalance of cell i */
3454 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3455 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3456 /* Determine the change of the cell size using underrelaxation */
3457 change = -sc*imbalance;
3458 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3462 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3463 cellsize_limit_f *= DD_CELL_MARGIN;
3464 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3465 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3466 if (ddbox->tric_dir[dim])
3468 cellsize_limit_f /= ddbox->skew_fac[dim];
3469 dist_min_f /= ddbox->skew_fac[dim];
3471 if (bDynamicBox && d > 0)
3473 dist_min_f *= DD_PRES_SCALE_MARGIN;
3475 if (d > 0 && !bUniform)
3477 /* Make sure that the grid is not shifted too much */
3478 for (i = 1; i < ncd; i++)
3480 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3482 gmx_incons("Inconsistent DD boundary staggering limits!");
3484 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3485 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3488 root->bound_min[i] += 0.5*space;
3490 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3491 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3494 root->bound_max[i] += 0.5*space;
3499 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3501 root->cell_f_max0[i-1] + dist_min_f,
3502 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3503 root->cell_f_min1[i] - dist_min_f);
3508 root->cell_f[0] = 0;
3509 root->cell_f[ncd] = 1;
3510 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3513 /* After the checks above, the cells should obey the cut-off
3514 * restrictions, but it does not hurt to check.
3516 for (i = 0; i < ncd; i++)
3520 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3521 dim, i, root->cell_f[i], root->cell_f[i+1]);
3524 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3525 root->cell_f[i+1] - root->cell_f[i] <
3526 cellsize_limit_f/DD_CELL_MARGIN)
3530 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3531 gmx_step_str(step, buf), dim2char(dim), i,
3532 (root->cell_f[i+1] - root->cell_f[i])
3533 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3538 /* Store the cell boundaries of the lower dimensions at the end */
3539 for (d1 = 0; d1 < d; d1++)
3541 root->cell_f[pos++] = comm->cell_f0[d1];
3542 root->cell_f[pos++] = comm->cell_f1[d1];
3545 if (d < comm->npmedecompdim)
3547 /* The master determines the maximum shift for
3548 * the coordinate communication between separate PME nodes.
3550 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3552 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3555 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3559 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3560 gmx_ddbox_t *ddbox, int dimind)
3562 gmx_domdec_comm_t *comm;
3567 /* Set the cell dimensions */
3568 dim = dd->dim[dimind];
3569 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3570 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3571 if (dim >= ddbox->nboundeddim)
3573 comm->cell_x0[dim] += ddbox->box0[dim];
3574 comm->cell_x1[dim] += ddbox->box0[dim];
3578 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3579 int d, int dim, real *cell_f_row,
3582 gmx_domdec_comm_t *comm;
3588 /* Each node would only need to know two fractions,
3589 * but it is probably cheaper to broadcast the whole array.
3591 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3592 0, comm->mpi_comm_load[d]);
3594 /* Copy the fractions for this dimension from the buffer */
3595 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3596 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3597 /* The whole array was communicated, so set the buffer position */
3598 pos = dd->nc[dim] + 1;
3599 for (d1 = 0; d1 <= d; d1++)
3603 /* Copy the cell fractions of the lower dimensions */
3604 comm->cell_f0[d1] = cell_f_row[pos++];
3605 comm->cell_f1[d1] = cell_f_row[pos++];
3607 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3609 /* Convert the communicated shift from float to int */
3610 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3613 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3617 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3618 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3619 gmx_bool bUniform, gmx_large_int_t step)
3621 gmx_domdec_comm_t *comm;
3623 gmx_bool bRowMember, bRowRoot;
3628 for (d = 0; d < dd->ndim; d++)
3633 for (d1 = d; d1 < dd->ndim; d1++)
3635 if (dd->ci[dd->dim[d1]] > 0)
3648 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3649 ddbox, bDynamicBox, bUniform, step);
3650 cell_f_row = comm->root[d]->cell_f;
3654 cell_f_row = comm->cell_f_row;
3656 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3661 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3665 /* This function assumes the box is static and should therefore
3666 * not be called when the box has changed since the last
3667 * call to dd_partition_system.
3669 for (d = 0; d < dd->ndim; d++)
3671 relative_to_absolute_cell_bounds(dd, ddbox, d);
3677 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3678 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3679 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3680 gmx_wallcycle_t wcycle)
3682 gmx_domdec_comm_t *comm;
3689 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3690 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3691 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3693 else if (bDynamicBox)
3695 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3698 /* Set the dimensions for which no DD is used */
3699 for (dim = 0; dim < DIM; dim++)
3701 if (dd->nc[dim] == 1)
3703 comm->cell_x0[dim] = 0;
3704 comm->cell_x1[dim] = ddbox->box_size[dim];
3705 if (dim >= ddbox->nboundeddim)
3707 comm->cell_x0[dim] += ddbox->box0[dim];
3708 comm->cell_x1[dim] += ddbox->box0[dim];
3714 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3717 gmx_domdec_comm_dim_t *cd;
3719 for (d = 0; d < dd->ndim; d++)
3721 cd = &dd->comm->cd[d];
3722 np = npulse[dd->dim[d]];
3723 if (np > cd->np_nalloc)
3727 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3728 dim2char(dd->dim[d]), np);
3730 if (DDMASTER(dd) && cd->np_nalloc > 0)
3732 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3734 srenew(cd->ind, np);
3735 for (i = cd->np_nalloc; i < np; i++)
3737 cd->ind[i].index = NULL;
3738 cd->ind[i].nalloc = 0;
3747 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3748 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3749 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3750 gmx_wallcycle_t wcycle)
3752 gmx_domdec_comm_t *comm;
3758 /* Copy the old cell boundaries for the cg displacement check */
3759 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3760 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3762 if (comm->bDynLoadBal)
3766 check_box_size(dd, ddbox);
3768 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3772 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3773 realloc_comm_ind(dd, npulse);
3778 for (d = 0; d < DIM; d++)
3780 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3781 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3786 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3788 rvec cell_ns_x0, rvec cell_ns_x1,
3789 gmx_large_int_t step)
3791 gmx_domdec_comm_t *comm;
3796 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3798 dim = dd->dim[dim_ind];
3800 /* Without PBC we don't have restrictions on the outer cells */
3801 if (!(dim >= ddbox->npbcdim &&
3802 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3803 comm->bDynLoadBal &&
3804 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3805 comm->cellsize_min[dim])
3808 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",
3809 gmx_step_str(step, buf), dim2char(dim),
3810 comm->cell_x1[dim] - comm->cell_x0[dim],
3811 ddbox->skew_fac[dim],
3812 dd->comm->cellsize_min[dim],
3813 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3817 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3819 /* Communicate the boundaries and update cell_ns_x0/1 */
3820 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3821 if (dd->bGridJump && dd->ndim > 1)
3823 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3828 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3832 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3840 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3841 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3850 static void check_screw_box(matrix box)
3852 /* Mathematical limitation */
3853 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3855 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3858 /* Limitation due to the asymmetry of the eighth shell method */
3859 if (box[ZZ][YY] != 0)
3861 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3865 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3866 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3869 gmx_domdec_master_t *ma;
3870 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3871 int i, icg, j, k, k0, k1, d, npbcdim;
3873 rvec box_size, cg_cm;
3875 real nrcg, inv_ncg, pos_d;
3877 gmx_bool bUnbounded, bScrew;
3881 if (tmp_ind == NULL)
3883 snew(tmp_nalloc, dd->nnodes);
3884 snew(tmp_ind, dd->nnodes);
3885 for (i = 0; i < dd->nnodes; i++)
3887 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3888 snew(tmp_ind[i], tmp_nalloc[i]);
3892 /* Clear the count */
3893 for (i = 0; i < dd->nnodes; i++)
3899 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3901 cgindex = cgs->index;
3903 /* Compute the center of geometry for all charge groups */
3904 for (icg = 0; icg < cgs->nr; icg++)
3907 k1 = cgindex[icg+1];
3911 copy_rvec(pos[k0], cg_cm);
3918 for (k = k0; (k < k1); k++)
3920 rvec_inc(cg_cm, pos[k]);
3922 for (d = 0; (d < DIM); d++)
3924 cg_cm[d] *= inv_ncg;
3927 /* Put the charge group in the box and determine the cell index */
3928 for (d = DIM-1; d >= 0; d--)
3931 if (d < dd->npbcdim)
3933 bScrew = (dd->bScrewPBC && d == XX);
3934 if (tric_dir[d] && dd->nc[d] > 1)
3936 /* Use triclinic coordintates for this dimension */
3937 for (j = d+1; j < DIM; j++)
3939 pos_d += cg_cm[j]*tcm[j][d];
3942 while (pos_d >= box[d][d])
3945 rvec_dec(cg_cm, box[d]);
3948 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3949 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3951 for (k = k0; (k < k1); k++)
3953 rvec_dec(pos[k], box[d]);
3956 pos[k][YY] = box[YY][YY] - pos[k][YY];
3957 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3964 rvec_inc(cg_cm, box[d]);
3967 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3968 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3970 for (k = k0; (k < k1); k++)
3972 rvec_inc(pos[k], box[d]);
3975 pos[k][YY] = box[YY][YY] - pos[k][YY];
3976 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3981 /* This could be done more efficiently */
3983 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3988 i = dd_index(dd->nc, ind);
3989 if (ma->ncg[i] == tmp_nalloc[i])
3991 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3992 srenew(tmp_ind[i], tmp_nalloc[i]);
3994 tmp_ind[i][ma->ncg[i]] = icg;
3996 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4000 for (i = 0; i < dd->nnodes; i++)
4003 for (k = 0; k < ma->ncg[i]; k++)
4005 ma->cg[k1++] = tmp_ind[i][k];
4008 ma->index[dd->nnodes] = k1;
4010 for (i = 0; i < dd->nnodes; i++)
4020 fprintf(fplog, "Charge group distribution at step %s:",
4021 gmx_step_str(step, buf));
4022 for (i = 0; i < dd->nnodes; i++)
4024 fprintf(fplog, " %d", ma->ncg[i]);
4026 fprintf(fplog, "\n");
4030 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4031 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4034 gmx_domdec_master_t *ma = NULL;
4037 int *ibuf, buf2[2] = { 0, 0 };
4038 gmx_bool bMaster = DDMASTER(dd);
4045 check_screw_box(box);
4048 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4050 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4051 for (i = 0; i < dd->nnodes; i++)
4053 ma->ibuf[2*i] = ma->ncg[i];
4054 ma->ibuf[2*i+1] = ma->nat[i];
4062 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4064 dd->ncg_home = buf2[0];
4065 dd->nat_home = buf2[1];
4066 dd->ncg_tot = dd->ncg_home;
4067 dd->nat_tot = dd->nat_home;
4068 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4070 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4071 srenew(dd->index_gl, dd->cg_nalloc);
4072 srenew(dd->cgindex, dd->cg_nalloc+1);
4076 for (i = 0; i < dd->nnodes; i++)
4078 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4079 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4084 DDMASTER(dd) ? ma->ibuf : NULL,
4085 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4086 DDMASTER(dd) ? ma->cg : NULL,
4087 dd->ncg_home*sizeof(int), dd->index_gl);
4089 /* Determine the home charge group sizes */
4091 for (i = 0; i < dd->ncg_home; i++)
4093 cg_gl = dd->index_gl[i];
4095 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4100 fprintf(debug, "Home charge groups:\n");
4101 for (i = 0; i < dd->ncg_home; i++)
4103 fprintf(debug, " %d", dd->index_gl[i]);
4106 fprintf(debug, "\n");
4109 fprintf(debug, "\n");
4113 static int compact_and_copy_vec_at(int ncg, int *move,
4116 rvec *src, gmx_domdec_comm_t *comm,
4119 int m, icg, i, i0, i1, nrcg;
4125 for (m = 0; m < DIM*2; m++)
4131 for (icg = 0; icg < ncg; icg++)
4133 i1 = cgindex[icg+1];
4139 /* Compact the home array in place */
4140 for (i = i0; i < i1; i++)
4142 copy_rvec(src[i], src[home_pos++]);
4148 /* Copy to the communication buffer */
4150 pos_vec[m] += 1 + vec*nrcg;
4151 for (i = i0; i < i1; i++)
4153 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4155 pos_vec[m] += (nvec - vec - 1)*nrcg;
4159 home_pos += i1 - i0;
4167 static int compact_and_copy_vec_cg(int ncg, int *move,
4169 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4172 int m, icg, i0, i1, nrcg;
4178 for (m = 0; m < DIM*2; m++)
4184 for (icg = 0; icg < ncg; icg++)
4186 i1 = cgindex[icg+1];
4192 /* Compact the home array in place */
4193 copy_rvec(src[icg], src[home_pos++]);
4199 /* Copy to the communication buffer */
4200 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4201 pos_vec[m] += 1 + nrcg*nvec;
4213 static int compact_ind(int ncg, int *move,
4214 int *index_gl, int *cgindex,
4216 gmx_ga2la_t ga2la, char *bLocalCG,
4219 int cg, nat, a0, a1, a, a_gl;
4224 for (cg = 0; cg < ncg; cg++)
4230 /* Compact the home arrays in place.
4231 * Anything that can be done here avoids access to global arrays.
4233 cgindex[home_pos] = nat;
4234 for (a = a0; a < a1; a++)
4237 gatindex[nat] = a_gl;
4238 /* The cell number stays 0, so we don't need to set it */
4239 ga2la_change_la(ga2la, a_gl, nat);
4242 index_gl[home_pos] = index_gl[cg];
4243 cginfo[home_pos] = cginfo[cg];
4244 /* The charge group remains local, so bLocalCG does not change */
4249 /* Clear the global indices */
4250 for (a = a0; a < a1; a++)
4252 ga2la_del(ga2la, gatindex[a]);
4256 bLocalCG[index_gl[cg]] = FALSE;
4260 cgindex[home_pos] = nat;
4265 static void clear_and_mark_ind(int ncg, int *move,
4266 int *index_gl, int *cgindex, int *gatindex,
4267 gmx_ga2la_t ga2la, char *bLocalCG,
4272 for (cg = 0; cg < ncg; cg++)
4278 /* Clear the global indices */
4279 for (a = a0; a < a1; a++)
4281 ga2la_del(ga2la, gatindex[a]);
4285 bLocalCG[index_gl[cg]] = FALSE;
4287 /* Signal that this cg has moved using the ns cell index.
4288 * Here we set it to -1. fill_grid will change it
4289 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4291 cell_index[cg] = -1;
4296 static void print_cg_move(FILE *fplog,
4298 gmx_large_int_t step, int cg, int dim, int dir,
4299 gmx_bool bHaveLimitdAndCMOld, real limitd,
4300 rvec cm_old, rvec cm_new, real pos_d)
4302 gmx_domdec_comm_t *comm;
4307 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4308 if (bHaveLimitdAndCMOld)
4310 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4311 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4315 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4316 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4318 fprintf(fplog, "distance out of cell %f\n",
4319 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4320 if (bHaveLimitdAndCMOld)
4322 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4323 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4325 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4326 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4327 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4329 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4330 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4332 comm->cell_x0[dim], comm->cell_x1[dim]);
4335 static void cg_move_error(FILE *fplog,
4337 gmx_large_int_t step, int cg, int dim, int dir,
4338 gmx_bool bHaveLimitdAndCMOld, real limitd,
4339 rvec cm_old, rvec cm_new, real pos_d)
4343 print_cg_move(fplog, dd, step, cg, dim, dir,
4344 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4346 print_cg_move(stderr, dd, step, cg, dim, dir,
4347 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4349 "A charge group moved too far between two domain decomposition steps\n"
4350 "This usually means that your system is not well equilibrated");
4353 static void rotate_state_atom(t_state *state, int a)
4357 for (est = 0; est < estNR; est++)
4359 if (EST_DISTR(est) && (state->flags & (1<<est)))
4364 /* Rotate the complete state; for a rectangular box only */
4365 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4366 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4369 state->v[a][YY] = -state->v[a][YY];
4370 state->v[a][ZZ] = -state->v[a][ZZ];
4373 state->sd_X[a][YY] = -state->sd_X[a][YY];
4374 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4377 state->cg_p[a][YY] = -state->cg_p[a][YY];
4378 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4380 case estDISRE_INITF:
4381 case estDISRE_RM3TAV:
4382 case estORIRE_INITF:
4384 /* These are distances, so not affected by rotation */
4387 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4393 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4395 if (natoms > comm->moved_nalloc)
4397 /* Contents should be preserved here */
4398 comm->moved_nalloc = over_alloc_dd(natoms);
4399 srenew(comm->moved, comm->moved_nalloc);
4405 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4408 ivec tric_dir, matrix tcm,
4409 rvec cell_x0, rvec cell_x1,
4410 rvec limitd, rvec limit0, rvec limit1,
4412 int cg_start, int cg_end,
4417 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4418 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4422 real inv_ncg, pos_d;
4425 npbcdim = dd->npbcdim;
4427 for (cg = cg_start; cg < cg_end; cg++)
4434 copy_rvec(state->x[k0], cm_new);
4441 for (k = k0; (k < k1); k++)
4443 rvec_inc(cm_new, state->x[k]);
4445 for (d = 0; (d < DIM); d++)
4447 cm_new[d] = inv_ncg*cm_new[d];
4452 /* Do pbc and check DD cell boundary crossings */
4453 for (d = DIM-1; d >= 0; d--)
4457 bScrew = (dd->bScrewPBC && d == XX);
4458 /* Determine the location of this cg in lattice coordinates */
4462 for (d2 = d+1; d2 < DIM; d2++)
4464 pos_d += cm_new[d2]*tcm[d2][d];
4467 /* Put the charge group in the triclinic unit-cell */
4468 if (pos_d >= cell_x1[d])
4470 if (pos_d >= limit1[d])
4472 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4473 cg_cm[cg], cm_new, pos_d);
4476 if (dd->ci[d] == dd->nc[d] - 1)
4478 rvec_dec(cm_new, state->box[d]);
4481 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4482 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4484 for (k = k0; (k < k1); k++)
4486 rvec_dec(state->x[k], state->box[d]);
4489 rotate_state_atom(state, k);
4494 else if (pos_d < cell_x0[d])
4496 if (pos_d < limit0[d])
4498 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4499 cg_cm[cg], cm_new, pos_d);
4504 rvec_inc(cm_new, state->box[d]);
4507 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4508 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4510 for (k = k0; (k < k1); k++)
4512 rvec_inc(state->x[k], state->box[d]);
4515 rotate_state_atom(state, k);
4521 else if (d < npbcdim)
4523 /* Put the charge group in the rectangular unit-cell */
4524 while (cm_new[d] >= state->box[d][d])
4526 rvec_dec(cm_new, state->box[d]);
4527 for (k = k0; (k < k1); k++)
4529 rvec_dec(state->x[k], state->box[d]);
4532 while (cm_new[d] < 0)
4534 rvec_inc(cm_new, state->box[d]);
4535 for (k = k0; (k < k1); k++)
4537 rvec_inc(state->x[k], state->box[d]);
4543 copy_rvec(cm_new, cg_cm[cg]);
4545 /* Determine where this cg should go */
4548 for (d = 0; d < dd->ndim; d++)
4553 flag |= DD_FLAG_FW(d);
4559 else if (dev[dim] == -1)
4561 flag |= DD_FLAG_BW(d);
4564 if (dd->nc[dim] > 2)
4575 /* Temporarily store the flag in move */
4576 move[cg] = mc + flag;
4580 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4581 gmx_domdec_t *dd, ivec tric_dir,
4582 t_state *state, rvec **f,
4583 t_forcerec *fr, t_mdatoms *md,
4591 int ncg[DIM*2], nat[DIM*2];
4592 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4593 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4594 int sbuf[2], rbuf[2];
4595 int home_pos_cg, home_pos_at, buf_pos;
4597 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4600 real inv_ncg, pos_d;
4602 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4604 cginfo_mb_t *cginfo_mb;
4605 gmx_domdec_comm_t *comm;
4607 int nthread, thread;
4611 check_screw_box(state->box);
4615 if (fr->cutoff_scheme == ecutsGROUP)
4620 for (i = 0; i < estNR; i++)
4626 case estX: /* Always present */ break;
4627 case estV: bV = (state->flags & (1<<i)); break;
4628 case estSDX: bSDX = (state->flags & (1<<i)); break;
4629 case estCGP: bCGP = (state->flags & (1<<i)); break;
4632 case estDISRE_INITF:
4633 case estDISRE_RM3TAV:
4634 case estORIRE_INITF:
4636 /* No processing required */
4639 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4644 if (dd->ncg_tot > comm->nalloc_int)
4646 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4647 srenew(comm->buf_int, comm->nalloc_int);
4649 move = comm->buf_int;
4651 /* Clear the count */
4652 for (c = 0; c < dd->ndim*2; c++)
4658 npbcdim = dd->npbcdim;
4660 for (d = 0; (d < DIM); d++)
4662 limitd[d] = dd->comm->cellsize_min[d];
4663 if (d >= npbcdim && dd->ci[d] == 0)
4665 cell_x0[d] = -GMX_FLOAT_MAX;
4669 cell_x0[d] = comm->cell_x0[d];
4671 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4673 cell_x1[d] = GMX_FLOAT_MAX;
4677 cell_x1[d] = comm->cell_x1[d];
4681 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4682 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4686 /* We check after communication if a charge group moved
4687 * more than one cell. Set the pre-comm check limit to float_max.
4689 limit0[d] = -GMX_FLOAT_MAX;
4690 limit1[d] = GMX_FLOAT_MAX;
4694 make_tric_corr_matrix(npbcdim, state->box, tcm);
4696 cgindex = dd->cgindex;
4698 nthread = gmx_omp_nthreads_get(emntDomdec);
4700 /* Compute the center of geometry for all home charge groups
4701 * and put them in the box and determine where they should go.
4703 #pragma omp parallel for num_threads(nthread) schedule(static)
4704 for (thread = 0; thread < nthread; thread++)
4706 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4707 cell_x0, cell_x1, limitd, limit0, limit1,
4709 ( thread *dd->ncg_home)/nthread,
4710 ((thread+1)*dd->ncg_home)/nthread,
4711 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4715 for (cg = 0; cg < dd->ncg_home; cg++)
4720 flag = mc & ~DD_FLAG_NRCG;
4721 mc = mc & DD_FLAG_NRCG;
4724 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4726 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4727 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4729 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4730 /* We store the cg size in the lower 16 bits
4731 * and the place where the charge group should go
4732 * in the next 6 bits. This saves some communication volume.
4734 nrcg = cgindex[cg+1] - cgindex[cg];
4735 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4741 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4742 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4745 for (i = 0; i < dd->ndim*2; i++)
4747 *ncg_moved += ncg[i];
4764 /* Make sure the communication buffers are large enough */
4765 for (mc = 0; mc < dd->ndim*2; mc++)
4767 nvr = ncg[mc] + nat[mc]*nvec;
4768 if (nvr > comm->cgcm_state_nalloc[mc])
4770 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4771 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4775 switch (fr->cutoff_scheme)
4778 /* Recalculating cg_cm might be cheaper than communicating,
4779 * but that could give rise to rounding issues.
4782 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4783 nvec, cg_cm, comm, bCompact);
4786 /* Without charge groups we send the moved atom coordinates
4787 * over twice. This is so the code below can be used without
4788 * many conditionals for both for with and without charge groups.
4791 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4792 nvec, state->x, comm, FALSE);
4795 home_pos_cg -= *ncg_moved;
4799 gmx_incons("unimplemented");
4805 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4806 nvec, vec++, state->x, comm, bCompact);
4809 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4810 nvec, vec++, state->v, comm, bCompact);
4814 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4815 nvec, vec++, state->sd_X, comm, bCompact);
4819 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4820 nvec, vec++, state->cg_p, comm, bCompact);
4825 compact_ind(dd->ncg_home, move,
4826 dd->index_gl, dd->cgindex, dd->gatindex,
4827 dd->ga2la, comm->bLocalCG,
4832 if (fr->cutoff_scheme == ecutsVERLET)
4834 moved = get_moved(comm, dd->ncg_home);
4836 for (k = 0; k < dd->ncg_home; k++)
4843 moved = fr->ns.grid->cell_index;
4846 clear_and_mark_ind(dd->ncg_home, move,
4847 dd->index_gl, dd->cgindex, dd->gatindex,
4848 dd->ga2la, comm->bLocalCG,
4852 cginfo_mb = fr->cginfo_mb;
4854 *ncg_stay_home = home_pos_cg;
4855 for (d = 0; d < dd->ndim; d++)
4861 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4864 /* Communicate the cg and atom counts */
4869 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4870 d, dir, sbuf[0], sbuf[1]);
4872 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4874 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4876 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4877 srenew(comm->buf_int, comm->nalloc_int);
4880 /* Communicate the charge group indices, sizes and flags */
4881 dd_sendrecv_int(dd, d, dir,
4882 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4883 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4885 nvs = ncg[cdd] + nat[cdd]*nvec;
4886 i = rbuf[0] + rbuf[1] *nvec;
4887 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4889 /* Communicate cgcm and state */
4890 dd_sendrecv_rvec(dd, d, dir,
4891 comm->cgcm_state[cdd], nvs,
4892 comm->vbuf.v+nvr, i);
4893 ncg_recv += rbuf[0];
4894 nat_recv += rbuf[1];
4898 /* Process the received charge groups */
4900 for (cg = 0; cg < ncg_recv; cg++)
4902 flag = comm->buf_int[cg*DD_CGIBS+1];
4904 if (dim >= npbcdim && dd->nc[dim] > 2)
4906 /* No pbc in this dim and more than one domain boundary.
4907 * We do a separate check if a charge group didn't move too far.
4909 if (((flag & DD_FLAG_FW(d)) &&
4910 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4911 ((flag & DD_FLAG_BW(d)) &&
4912 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4914 cg_move_error(fplog, dd, step, cg, dim,
4915 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4917 comm->vbuf.v[buf_pos],
4918 comm->vbuf.v[buf_pos],
4919 comm->vbuf.v[buf_pos][dim]);
4926 /* Check which direction this cg should go */
4927 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4931 /* The cell boundaries for dimension d2 are not equal
4932 * for each cell row of the lower dimension(s),
4933 * therefore we might need to redetermine where
4934 * this cg should go.
4937 /* If this cg crosses the box boundary in dimension d2
4938 * we can use the communicated flag, so we do not
4939 * have to worry about pbc.
4941 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4942 (flag & DD_FLAG_FW(d2))) ||
4943 (dd->ci[dim2] == 0 &&
4944 (flag & DD_FLAG_BW(d2)))))
4946 /* Clear the two flags for this dimension */
4947 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4948 /* Determine the location of this cg
4949 * in lattice coordinates
4951 pos_d = comm->vbuf.v[buf_pos][dim2];
4954 for (d3 = dim2+1; d3 < DIM; d3++)
4957 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4960 /* Check of we are not at the box edge.
4961 * pbc is only handled in the first step above,
4962 * but this check could move over pbc while
4963 * the first step did not due to different rounding.
4965 if (pos_d >= cell_x1[dim2] &&
4966 dd->ci[dim2] != dd->nc[dim2]-1)
4968 flag |= DD_FLAG_FW(d2);
4970 else if (pos_d < cell_x0[dim2] &&
4973 flag |= DD_FLAG_BW(d2);
4975 comm->buf_int[cg*DD_CGIBS+1] = flag;
4978 /* Set to which neighboring cell this cg should go */
4979 if (flag & DD_FLAG_FW(d2))
4983 else if (flag & DD_FLAG_BW(d2))
4985 if (dd->nc[dd->dim[d2]] > 2)
4997 nrcg = flag & DD_FLAG_NRCG;
5000 if (home_pos_cg+1 > dd->cg_nalloc)
5002 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5003 srenew(dd->index_gl, dd->cg_nalloc);
5004 srenew(dd->cgindex, dd->cg_nalloc+1);
5006 /* Set the global charge group index and size */
5007 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5008 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5009 /* Copy the state from the buffer */
5010 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5011 if (fr->cutoff_scheme == ecutsGROUP)
5014 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5018 /* Set the cginfo */
5019 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5020 dd->index_gl[home_pos_cg]);
5023 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5026 if (home_pos_at+nrcg > state->nalloc)
5028 dd_realloc_state(state, f, home_pos_at+nrcg);
5030 for (i = 0; i < nrcg; i++)
5032 copy_rvec(comm->vbuf.v[buf_pos++],
5033 state->x[home_pos_at+i]);
5037 for (i = 0; i < nrcg; i++)
5039 copy_rvec(comm->vbuf.v[buf_pos++],
5040 state->v[home_pos_at+i]);
5045 for (i = 0; i < nrcg; i++)
5047 copy_rvec(comm->vbuf.v[buf_pos++],
5048 state->sd_X[home_pos_at+i]);
5053 for (i = 0; i < nrcg; i++)
5055 copy_rvec(comm->vbuf.v[buf_pos++],
5056 state->cg_p[home_pos_at+i]);
5060 home_pos_at += nrcg;
5064 /* Reallocate the buffers if necessary */
5065 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5067 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5068 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5070 nvr = ncg[mc] + nat[mc]*nvec;
5071 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5073 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5074 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5076 /* Copy from the receive to the send buffers */
5077 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5078 comm->buf_int + cg*DD_CGIBS,
5079 DD_CGIBS*sizeof(int));
5080 memcpy(comm->cgcm_state[mc][nvr],
5081 comm->vbuf.v[buf_pos],
5082 (1+nrcg*nvec)*sizeof(rvec));
5083 buf_pos += 1 + nrcg*nvec;
5090 /* With sorting (!bCompact) the indices are now only partially up to date
5091 * and ncg_home and nat_home are not the real count, since there are
5092 * "holes" in the arrays for the charge groups that moved to neighbors.
5094 if (fr->cutoff_scheme == ecutsVERLET)
5096 moved = get_moved(comm, home_pos_cg);
5098 for (i = dd->ncg_home; i < home_pos_cg; i++)
5103 dd->ncg_home = home_pos_cg;
5104 dd->nat_home = home_pos_at;
5109 "Finished repartitioning: cgs moved out %d, new home %d\n",
5110 *ncg_moved, dd->ncg_home-*ncg_moved);
5115 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5117 dd->comm->cycl[ddCycl] += cycles;
5118 dd->comm->cycl_n[ddCycl]++;
5119 if (cycles > dd->comm->cycl_max[ddCycl])
5121 dd->comm->cycl_max[ddCycl] = cycles;
5125 static double force_flop_count(t_nrnb *nrnb)
5132 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5134 /* To get closer to the real timings, we half the count
5135 * for the normal loops and again half it for water loops.
5138 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5140 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5144 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5147 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5150 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5152 sum += nrnb->n[i]*cost_nrnb(i);
5155 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5157 sum += nrnb->n[i]*cost_nrnb(i);
5163 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5165 if (dd->comm->eFlop)
5167 dd->comm->flop -= force_flop_count(nrnb);
5170 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5172 if (dd->comm->eFlop)
5174 dd->comm->flop += force_flop_count(nrnb);
5179 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5183 for (i = 0; i < ddCyclNr; i++)
5185 dd->comm->cycl[i] = 0;
5186 dd->comm->cycl_n[i] = 0;
5187 dd->comm->cycl_max[i] = 0;
5190 dd->comm->flop_n = 0;
5193 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5195 gmx_domdec_comm_t *comm;
5196 gmx_domdec_load_t *load;
5197 gmx_domdec_root_t *root = NULL;
5198 int d, dim, cid, i, pos;
5199 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5204 fprintf(debug, "get_load_distribution start\n");
5207 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5211 bSepPME = (dd->pme_nodeid >= 0);
5213 for (d = dd->ndim-1; d >= 0; d--)
5216 /* Check if we participate in the communication in this dimension */
5217 if (d == dd->ndim-1 ||
5218 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5220 load = &comm->load[d];
5223 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5226 if (d == dd->ndim-1)
5228 sbuf[pos++] = dd_force_load(comm);
5229 sbuf[pos++] = sbuf[0];
5232 sbuf[pos++] = sbuf[0];
5233 sbuf[pos++] = cell_frac;
5236 sbuf[pos++] = comm->cell_f_max0[d];
5237 sbuf[pos++] = comm->cell_f_min1[d];
5242 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5243 sbuf[pos++] = comm->cycl[ddCyclPME];
5248 sbuf[pos++] = comm->load[d+1].sum;
5249 sbuf[pos++] = comm->load[d+1].max;
5252 sbuf[pos++] = comm->load[d+1].sum_m;
5253 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5254 sbuf[pos++] = comm->load[d+1].flags;
5257 sbuf[pos++] = comm->cell_f_max0[d];
5258 sbuf[pos++] = comm->cell_f_min1[d];
5263 sbuf[pos++] = comm->load[d+1].mdf;
5264 sbuf[pos++] = comm->load[d+1].pme;
5268 /* Communicate a row in DD direction d.
5269 * The communicators are setup such that the root always has rank 0.
5272 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5273 load->load, load->nload*sizeof(float), MPI_BYTE,
5274 0, comm->mpi_comm_load[d]);
5276 if (dd->ci[dim] == dd->master_ci[dim])
5278 /* We are the root, process this row */
5279 if (comm->bDynLoadBal)
5281 root = comm->root[d];
5291 for (i = 0; i < dd->nc[dim]; i++)
5293 load->sum += load->load[pos++];
5294 load->max = max(load->max, load->load[pos]);
5300 /* This direction could not be load balanced properly,
5301 * therefore we need to use the maximum iso the average load.
5303 load->sum_m = max(load->sum_m, load->load[pos]);
5307 load->sum_m += load->load[pos];
5310 load->cvol_min = min(load->cvol_min, load->load[pos]);
5314 load->flags = (int)(load->load[pos++] + 0.5);
5318 root->cell_f_max0[i] = load->load[pos++];
5319 root->cell_f_min1[i] = load->load[pos++];
5324 load->mdf = max(load->mdf, load->load[pos]);
5326 load->pme = max(load->pme, load->load[pos]);
5330 if (comm->bDynLoadBal && root->bLimited)
5332 load->sum_m *= dd->nc[dim];
5333 load->flags |= (1<<d);
5341 comm->nload += dd_load_count(comm);
5342 comm->load_step += comm->cycl[ddCyclStep];
5343 comm->load_sum += comm->load[0].sum;
5344 comm->load_max += comm->load[0].max;
5345 if (comm->bDynLoadBal)
5347 for (d = 0; d < dd->ndim; d++)
5349 if (comm->load[0].flags & (1<<d))
5351 comm->load_lim[d]++;
5357 comm->load_mdf += comm->load[0].mdf;
5358 comm->load_pme += comm->load[0].pme;
5362 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5366 fprintf(debug, "get_load_distribution finished\n");
5370 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5372 /* Return the relative performance loss on the total run time
5373 * due to the force calculation load imbalance.
5375 if (dd->comm->nload > 0)
5378 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5379 (dd->comm->load_step*dd->nnodes);
5387 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5390 int npp, npme, nnodes, d, limp;
5391 float imbal, pme_f_ratio, lossf, lossp = 0;
5393 gmx_domdec_comm_t *comm;
5396 if (DDMASTER(dd) && comm->nload > 0)
5399 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5400 nnodes = npp + npme;
5401 imbal = comm->load_max*npp/comm->load_sum - 1;
5402 lossf = dd_force_imb_perf_loss(dd);
5403 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5404 fprintf(fplog, "%s", buf);
5405 fprintf(stderr, "\n");
5406 fprintf(stderr, "%s", buf);
5407 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5408 fprintf(fplog, "%s", buf);
5409 fprintf(stderr, "%s", buf);
5411 if (comm->bDynLoadBal)
5413 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5414 for (d = 0; d < dd->ndim; d++)
5416 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5417 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5423 sprintf(buf+strlen(buf), "\n");
5424 fprintf(fplog, "%s", buf);
5425 fprintf(stderr, "%s", buf);
5429 pme_f_ratio = comm->load_pme/comm->load_mdf;
5430 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5433 lossp *= (float)npme/(float)nnodes;
5437 lossp *= (float)npp/(float)nnodes;
5439 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5440 fprintf(fplog, "%s", buf);
5441 fprintf(stderr, "%s", buf);
5442 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5443 fprintf(fplog, "%s", buf);
5444 fprintf(stderr, "%s", buf);
5446 fprintf(fplog, "\n");
5447 fprintf(stderr, "\n");
5449 if (lossf >= DD_PERF_LOSS)
5452 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5453 " in the domain decomposition.\n", lossf*100);
5454 if (!comm->bDynLoadBal)
5456 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5460 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5462 fprintf(fplog, "%s\n", buf);
5463 fprintf(stderr, "%s\n", buf);
5465 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5468 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5469 " had %s work to do than the PP nodes.\n"
5470 " You might want to %s the number of PME nodes\n"
5471 " or %s the cut-off and the grid spacing.\n",
5473 (lossp < 0) ? "less" : "more",
5474 (lossp < 0) ? "decrease" : "increase",
5475 (lossp < 0) ? "decrease" : "increase");
5476 fprintf(fplog, "%s\n", buf);
5477 fprintf(stderr, "%s\n", buf);
5482 static float dd_vol_min(gmx_domdec_t *dd)
5484 return dd->comm->load[0].cvol_min*dd->nnodes;
5487 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5489 return dd->comm->load[0].flags;
5492 static float dd_f_imbal(gmx_domdec_t *dd)
5494 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5497 float dd_pme_f_ratio(gmx_domdec_t *dd)
5499 if (dd->comm->cycl_n[ddCyclPME] > 0)
5501 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5509 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5514 flags = dd_load_flags(dd);
5518 "DD load balancing is limited by minimum cell size in dimension");
5519 for (d = 0; d < dd->ndim; d++)
5523 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5526 fprintf(fplog, "\n");
5528 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5529 if (dd->comm->bDynLoadBal)
5531 fprintf(fplog, " vol min/aver %5.3f%c",
5532 dd_vol_min(dd), flags ? '!' : ' ');
5534 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5535 if (dd->comm->cycl_n[ddCyclPME])
5537 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5539 fprintf(fplog, "\n\n");
5542 static void dd_print_load_verbose(gmx_domdec_t *dd)
5544 if (dd->comm->bDynLoadBal)
5546 fprintf(stderr, "vol %4.2f%c ",
5547 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5549 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5550 if (dd->comm->cycl_n[ddCyclPME])
5552 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5557 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5562 gmx_domdec_root_t *root;
5563 gmx_bool bPartOfGroup = FALSE;
5565 dim = dd->dim[dim_ind];
5566 copy_ivec(loc, loc_c);
5567 for (i = 0; i < dd->nc[dim]; i++)
5570 rank = dd_index(dd->nc, loc_c);
5571 if (rank == dd->rank)
5573 /* This process is part of the group */
5574 bPartOfGroup = TRUE;
5577 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5581 dd->comm->mpi_comm_load[dim_ind] = c_row;
5582 if (dd->comm->eDLB != edlbNO)
5584 if (dd->ci[dim] == dd->master_ci[dim])
5586 /* This is the root process of this row */
5587 snew(dd->comm->root[dim_ind], 1);
5588 root = dd->comm->root[dim_ind];
5589 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5590 snew(root->old_cell_f, dd->nc[dim]+1);
5591 snew(root->bCellMin, dd->nc[dim]);
5594 snew(root->cell_f_max0, dd->nc[dim]);
5595 snew(root->cell_f_min1, dd->nc[dim]);
5596 snew(root->bound_min, dd->nc[dim]);
5597 snew(root->bound_max, dd->nc[dim]);
5599 snew(root->buf_ncd, dd->nc[dim]);
5603 /* This is not a root process, we only need to receive cell_f */
5604 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5607 if (dd->ci[dim] == dd->master_ci[dim])
5609 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5615 static void make_load_communicators(gmx_domdec_t *dd)
5618 int dim0, dim1, i, j;
5623 fprintf(debug, "Making load communicators\n");
5626 snew(dd->comm->load, dd->ndim);
5627 snew(dd->comm->mpi_comm_load, dd->ndim);
5630 make_load_communicator(dd, 0, loc);
5634 for (i = 0; i < dd->nc[dim0]; i++)
5637 make_load_communicator(dd, 1, loc);
5643 for (i = 0; i < dd->nc[dim0]; i++)
5647 for (j = 0; j < dd->nc[dim1]; j++)
5650 make_load_communicator(dd, 2, loc);
5657 fprintf(debug, "Finished making load communicators\n");
5662 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5665 int d, dim, i, j, m;
5668 ivec dd_zp[DD_MAXIZONE];
5669 gmx_domdec_zones_t *zones;
5670 gmx_domdec_ns_ranges_t *izone;
5672 for (d = 0; d < dd->ndim; d++)
5675 copy_ivec(dd->ci, tmp);
5676 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5677 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5678 copy_ivec(dd->ci, tmp);
5679 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5680 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5683 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5686 dd->neighbor[d][1]);
5692 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5694 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5695 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5702 for (i = 0; i < nzonep; i++)
5704 copy_ivec(dd_zp3[i], dd_zp[i]);
5710 for (i = 0; i < nzonep; i++)
5712 copy_ivec(dd_zp2[i], dd_zp[i]);
5718 for (i = 0; i < nzonep; i++)
5720 copy_ivec(dd_zp1[i], dd_zp[i]);
5724 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5729 zones = &dd->comm->zones;
5731 for (i = 0; i < nzone; i++)
5734 clear_ivec(zones->shift[i]);
5735 for (d = 0; d < dd->ndim; d++)
5737 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5742 for (i = 0; i < nzone; i++)
5744 for (d = 0; d < DIM; d++)
5746 s[d] = dd->ci[d] - zones->shift[i][d];
5751 else if (s[d] >= dd->nc[d])
5757 zones->nizone = nzonep;
5758 for (i = 0; i < zones->nizone; i++)
5760 if (dd_zp[i][0] != i)
5762 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5764 izone = &zones->izone[i];
5765 izone->j0 = dd_zp[i][1];
5766 izone->j1 = dd_zp[i][2];
5767 for (dim = 0; dim < DIM; dim++)
5769 if (dd->nc[dim] == 1)
5771 /* All shifts should be allowed */
5772 izone->shift0[dim] = -1;
5773 izone->shift1[dim] = 1;
5778 izone->shift0[d] = 0;
5779 izone->shift1[d] = 0;
5780 for(j=izone->j0; j<izone->j1; j++) {
5781 if (dd->shift[j][d] > dd->shift[i][d])
5782 izone->shift0[d] = -1;
5783 if (dd->shift[j][d] < dd->shift[i][d])
5784 izone->shift1[d] = 1;
5790 /* Assume the shift are not more than 1 cell */
5791 izone->shift0[dim] = 1;
5792 izone->shift1[dim] = -1;
5793 for (j = izone->j0; j < izone->j1; j++)
5795 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5796 if (shift_diff < izone->shift0[dim])
5798 izone->shift0[dim] = shift_diff;
5800 if (shift_diff > izone->shift1[dim])
5802 izone->shift1[dim] = shift_diff;
5809 if (dd->comm->eDLB != edlbNO)
5811 snew(dd->comm->root, dd->ndim);
5814 if (dd->comm->bRecordLoad)
5816 make_load_communicators(dd);
5820 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5823 gmx_domdec_comm_t *comm;
5834 if (comm->bCartesianPP)
5836 /* Set up cartesian communication for the particle-particle part */
5839 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5840 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5843 for (i = 0; i < DIM; i++)
5847 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5849 /* We overwrite the old communicator with the new cartesian one */
5850 cr->mpi_comm_mygroup = comm_cart;
5853 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5854 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5856 if (comm->bCartesianPP_PME)
5858 /* Since we want to use the original cartesian setup for sim,
5859 * and not the one after split, we need to make an index.
5861 snew(comm->ddindex2ddnodeid, dd->nnodes);
5862 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5863 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5864 /* Get the rank of the DD master,
5865 * above we made sure that the master node is a PP node.
5875 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5877 else if (comm->bCartesianPP)
5879 if (cr->npmenodes == 0)
5881 /* The PP communicator is also
5882 * the communicator for this simulation
5884 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5886 cr->nodeid = dd->rank;
5888 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5890 /* We need to make an index to go from the coordinates
5891 * to the nodeid of this simulation.
5893 snew(comm->ddindex2simnodeid, dd->nnodes);
5894 snew(buf, dd->nnodes);
5895 if (cr->duty & DUTY_PP)
5897 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5899 /* Communicate the ddindex to simulation nodeid index */
5900 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5901 cr->mpi_comm_mysim);
5904 /* Determine the master coordinates and rank.
5905 * The DD master should be the same node as the master of this sim.
5907 for (i = 0; i < dd->nnodes; i++)
5909 if (comm->ddindex2simnodeid[i] == 0)
5911 ddindex2xyz(dd->nc, i, dd->master_ci);
5912 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5917 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5922 /* No Cartesian communicators */
5923 /* We use the rank in dd->comm->all as DD index */
5924 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5925 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5927 clear_ivec(dd->master_ci);
5934 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5935 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5940 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5941 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5945 static void receive_ddindex2simnodeid(t_commrec *cr)
5949 gmx_domdec_comm_t *comm;
5956 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5958 snew(comm->ddindex2simnodeid, dd->nnodes);
5959 snew(buf, dd->nnodes);
5960 if (cr->duty & DUTY_PP)
5962 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5965 /* Communicate the ddindex to simulation nodeid index */
5966 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5967 cr->mpi_comm_mysim);
5974 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5975 int ncg, int natoms)
5977 gmx_domdec_master_t *ma;
5982 snew(ma->ncg, dd->nnodes);
5983 snew(ma->index, dd->nnodes+1);
5985 snew(ma->nat, dd->nnodes);
5986 snew(ma->ibuf, dd->nnodes*2);
5987 snew(ma->cell_x, DIM);
5988 for (i = 0; i < DIM; i++)
5990 snew(ma->cell_x[i], dd->nc[i]+1);
5993 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5999 snew(ma->vbuf, natoms);
6005 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
6009 gmx_domdec_comm_t *comm;
6020 if (comm->bCartesianPP)
6022 for (i = 1; i < DIM; i++)
6024 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6026 if (bDiv[YY] || bDiv[ZZ])
6028 comm->bCartesianPP_PME = TRUE;
6029 /* If we have 2D PME decomposition, which is always in x+y,
6030 * we stack the PME only nodes in z.
6031 * Otherwise we choose the direction that provides the thinnest slab
6032 * of PME only nodes as this will have the least effect
6033 * on the PP communication.
6034 * But for the PME communication the opposite might be better.
6036 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6038 dd->nc[YY] > dd->nc[ZZ]))
6040 comm->cartpmedim = ZZ;
6044 comm->cartpmedim = YY;
6046 comm->ntot[comm->cartpmedim]
6047 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6051 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]);
6053 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6058 if (comm->bCartesianPP_PME)
6062 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]);
6065 for (i = 0; i < DIM; i++)
6069 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6072 MPI_Comm_rank(comm_cart, &rank);
6073 if (MASTERNODE(cr) && rank != 0)
6075 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6078 /* With this assigment we loose the link to the original communicator
6079 * which will usually be MPI_COMM_WORLD, unless have multisim.
6081 cr->mpi_comm_mysim = comm_cart;
6082 cr->sim_nodeid = rank;
6084 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6088 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6089 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6092 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6096 if (cr->npmenodes == 0 ||
6097 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6099 cr->duty = DUTY_PME;
6102 /* Split the sim communicator into PP and PME only nodes */
6103 MPI_Comm_split(cr->mpi_comm_mysim,
6105 dd_index(comm->ntot, dd->ci),
6106 &cr->mpi_comm_mygroup);
6110 switch (dd_node_order)
6115 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6118 case ddnoINTERLEAVE:
6119 /* Interleave the PP-only and PME-only nodes,
6120 * as on clusters with dual-core machines this will double
6121 * the communication bandwidth of the PME processes
6122 * and thus speed up the PP <-> PME and inter PME communication.
6126 fprintf(fplog, "Interleaving PP and PME nodes\n");
6128 comm->pmenodes = dd_pmenodes(cr);
6133 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6136 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6138 cr->duty = DUTY_PME;
6145 /* Split the sim communicator into PP and PME only nodes */
6146 MPI_Comm_split(cr->mpi_comm_mysim,
6149 &cr->mpi_comm_mygroup);
6150 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6156 fprintf(fplog, "This is a %s only node\n\n",
6157 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6161 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6164 gmx_domdec_comm_t *comm;
6170 copy_ivec(dd->nc, comm->ntot);
6172 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6173 comm->bCartesianPP_PME = FALSE;
6175 /* Reorder the nodes by default. This might change the MPI ranks.
6176 * Real reordering is only supported on very few architectures,
6177 * Blue Gene is one of them.
6179 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6181 if (cr->npmenodes > 0)
6183 /* Split the communicator into a PP and PME part */
6184 split_communicator(fplog, cr, dd_node_order, CartReorder);
6185 if (comm->bCartesianPP_PME)
6187 /* We (possibly) reordered the nodes in split_communicator,
6188 * so it is no longer required in make_pp_communicator.
6190 CartReorder = FALSE;
6195 /* All nodes do PP and PME */
6197 /* We do not require separate communicators */
6198 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6202 if (cr->duty & DUTY_PP)
6204 /* Copy or make a new PP communicator */
6205 make_pp_communicator(fplog, cr, CartReorder);
6209 receive_ddindex2simnodeid(cr);
6212 if (!(cr->duty & DUTY_PME))
6214 /* Set up the commnuication to our PME node */
6215 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6216 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6219 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6220 dd->pme_nodeid, dd->pme_receive_vir_ener);
6225 dd->pme_nodeid = -1;
6230 dd->ma = init_gmx_domdec_master_t(dd,
6232 comm->cgs_gl.index[comm->cgs_gl.nr]);
6236 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6238 real *slb_frac, tot;
6243 if (nc > 1 && size_string != NULL)
6247 fprintf(fplog, "Using static load balancing for the %s direction\n",
6252 for (i = 0; i < nc; i++)
6255 sscanf(size_string, "%lf%n", &dbl, &n);
6258 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6267 fprintf(fplog, "Relative cell sizes:");
6269 for (i = 0; i < nc; i++)
6274 fprintf(fplog, " %5.3f", slb_frac[i]);
6279 fprintf(fplog, "\n");
6286 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6289 gmx_mtop_ilistloop_t iloop;
6293 iloop = gmx_mtop_ilistloop_init(mtop);
6294 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6296 for (ftype = 0; ftype < F_NRE; ftype++)
6298 if ((interaction_function[ftype].flags & IF_BOND) &&
6301 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6309 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6315 val = getenv(env_var);
6318 if (sscanf(val, "%d", &nst) <= 0)
6324 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6332 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6336 fprintf(stderr, "\n%s\n", warn_string);
6340 fprintf(fplog, "\n%s\n", warn_string);
6344 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6345 t_inputrec *ir, FILE *fplog)
6347 if (ir->ePBC == epbcSCREW &&
6348 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6350 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6353 if (ir->ns_type == ensSIMPLE)
6355 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6358 if (ir->nstlist == 0)
6360 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6363 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6365 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6369 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6374 r = ddbox->box_size[XX];
6375 for (di = 0; di < dd->ndim; di++)
6378 /* Check using the initial average cell size */
6379 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6385 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6386 const char *dlb_opt, gmx_bool bRecordLoad,
6387 unsigned long Flags, t_inputrec *ir)
6395 case 'a': eDLB = edlbAUTO; break;
6396 case 'n': eDLB = edlbNO; break;
6397 case 'y': eDLB = edlbYES; break;
6398 default: gmx_incons("Unknown dlb_opt");
6401 if (Flags & MD_RERUN)
6406 if (!EI_DYNAMICS(ir->eI))
6408 if (eDLB == edlbYES)
6410 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6411 dd_warning(cr, fplog, buf);
6419 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6424 if (Flags & MD_REPRODUCIBLE)
6431 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6435 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6438 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6446 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6451 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6453 /* Decomposition order z,y,x */
6456 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6458 for (dim = DIM-1; dim >= 0; dim--)
6460 if (dd->nc[dim] > 1)
6462 dd->dim[dd->ndim++] = dim;
6468 /* Decomposition order x,y,z */
6469 for (dim = 0; dim < DIM; dim++)
6471 if (dd->nc[dim] > 1)
6473 dd->dim[dd->ndim++] = dim;
6479 static gmx_domdec_comm_t *init_dd_comm()
6481 gmx_domdec_comm_t *comm;
6485 snew(comm->cggl_flag, DIM*2);
6486 snew(comm->cgcm_state, DIM*2);
6487 for (i = 0; i < DIM*2; i++)
6489 comm->cggl_flag_nalloc[i] = 0;
6490 comm->cgcm_state_nalloc[i] = 0;
6493 comm->nalloc_int = 0;
6494 comm->buf_int = NULL;
6496 vec_rvec_init(&comm->vbuf);
6498 comm->n_load_have = 0;
6499 comm->n_load_collect = 0;
6501 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6503 comm->sum_nat[i] = 0;
6507 comm->load_step = 0;
6510 clear_ivec(comm->load_lim);
6517 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6518 unsigned long Flags,
6520 real comm_distance_min, real rconstr,
6521 const char *dlb_opt, real dlb_scale,
6522 const char *sizex, const char *sizey, const char *sizez,
6523 gmx_mtop_t *mtop, t_inputrec *ir,
6524 matrix box, rvec *x,
6526 int *npme_x, int *npme_y)
6529 gmx_domdec_comm_t *comm;
6532 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6539 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6544 dd->comm = init_dd_comm();
6546 snew(comm->cggl_flag, DIM*2);
6547 snew(comm->cgcm_state, DIM*2);
6549 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6550 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6552 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6553 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6554 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6555 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6556 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6557 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6558 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6559 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6561 dd->pme_recv_f_alloc = 0;
6562 dd->pme_recv_f_buf = NULL;
6564 if (dd->bSendRecv2 && fplog)
6566 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");
6572 fprintf(fplog, "Will load balance based on FLOP count\n");
6574 if (comm->eFlop > 1)
6576 srand(1+cr->nodeid);
6578 comm->bRecordLoad = TRUE;
6582 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6586 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6588 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6591 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6593 dd->bGridJump = comm->bDynLoadBal;
6594 comm->bPMELoadBalDLBLimits = FALSE;
6596 if (comm->nstSortCG)
6600 if (comm->nstSortCG == 1)
6602 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6606 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6610 snew(comm->sort, 1);
6616 fprintf(fplog, "Will not sort the charge groups\n");
6620 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6622 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6623 if (comm->bInterCGBondeds)
6625 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6629 comm->bInterCGMultiBody = FALSE;
6632 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6633 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6635 if (ir->rlistlong == 0)
6637 /* Set the cut-off to some very large value,
6638 * so we don't need if statements everywhere in the code.
6639 * We use sqrt, since the cut-off is squared in some places.
6641 comm->cutoff = GMX_CUTOFF_INF;
6645 comm->cutoff = ir->rlistlong;
6647 comm->cutoff_mbody = 0;
6649 comm->cellsize_limit = 0;
6650 comm->bBondComm = FALSE;
6652 if (comm->bInterCGBondeds)
6654 if (comm_distance_min > 0)
6656 comm->cutoff_mbody = comm_distance_min;
6657 if (Flags & MD_DDBONDCOMM)
6659 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6663 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6665 r_bonded_limit = comm->cutoff_mbody;
6667 else if (ir->bPeriodicMols)
6669 /* Can not easily determine the required cut-off */
6670 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");
6671 comm->cutoff_mbody = comm->cutoff/2;
6672 r_bonded_limit = comm->cutoff_mbody;
6678 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6679 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6681 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6682 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6684 /* We use an initial margin of 10% for the minimum cell size,
6685 * except when we are just below the non-bonded cut-off.
6687 if (Flags & MD_DDBONDCOMM)
6689 if (max(r_2b, r_mb) > comm->cutoff)
6691 r_bonded = max(r_2b, r_mb);
6692 r_bonded_limit = 1.1*r_bonded;
6693 comm->bBondComm = TRUE;
6698 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6700 /* We determine cutoff_mbody later */
6704 /* No special bonded communication,
6705 * simply increase the DD cut-off.
6707 r_bonded_limit = 1.1*max(r_2b, r_mb);
6708 comm->cutoff_mbody = r_bonded_limit;
6709 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6712 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6716 "Minimum cell size due to bonded interactions: %.3f nm\n",
6717 comm->cellsize_limit);
6721 if (dd->bInterCGcons && rconstr <= 0)
6723 /* There is a cell size limit due to the constraints (P-LINCS) */
6724 rconstr = constr_r_max(fplog, mtop, ir);
6728 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6730 if (rconstr > comm->cellsize_limit)
6732 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6736 else if (rconstr > 0 && fplog)
6738 /* Here we do not check for dd->bInterCGcons,
6739 * because one can also set a cell size limit for virtual sites only
6740 * and at this point we don't know yet if there are intercg v-sites.
6743 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6746 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6748 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6752 copy_ivec(nc, dd->nc);
6753 set_dd_dim(fplog, dd);
6754 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6756 if (cr->npmenodes == -1)
6760 acs = average_cellsize_min(dd, ddbox);
6761 if (acs < comm->cellsize_limit)
6765 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6767 gmx_fatal_collective(FARGS, cr, NULL,
6768 "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",
6769 acs, comm->cellsize_limit);
6774 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6776 /* We need to choose the optimal DD grid and possibly PME nodes */
6777 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6778 comm->eDLB != edlbNO, dlb_scale,
6779 comm->cellsize_limit, comm->cutoff,
6780 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6782 if (dd->nc[XX] == 0)
6784 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6785 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6786 !bC ? "-rdd" : "-rcon",
6787 comm->eDLB != edlbNO ? " or -dds" : "",
6788 bC ? " or your LINCS settings" : "");
6790 gmx_fatal_collective(FARGS, cr, NULL,
6791 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6793 "Look in the log file for details on the domain decomposition",
6794 cr->nnodes-cr->npmenodes, limit, buf);
6796 set_dd_dim(fplog, dd);
6802 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6803 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6806 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6807 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6809 gmx_fatal_collective(FARGS, cr, NULL,
6810 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6811 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6813 if (cr->npmenodes > dd->nnodes)
6815 gmx_fatal_collective(FARGS, cr, NULL,
6816 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6818 if (cr->npmenodes > 0)
6820 comm->npmenodes = cr->npmenodes;
6824 comm->npmenodes = dd->nnodes;
6827 if (EEL_PME(ir->coulombtype))
6829 /* The following choices should match those
6830 * in comm_cost_est in domdec_setup.c.
6831 * Note that here the checks have to take into account
6832 * that the decomposition might occur in a different order than xyz
6833 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6834 * in which case they will not match those in comm_cost_est,
6835 * but since that is mainly for testing purposes that's fine.
6837 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6838 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6839 getenv("GMX_PMEONEDD") == NULL)
6841 comm->npmedecompdim = 2;
6842 comm->npmenodes_x = dd->nc[XX];
6843 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6847 /* In case nc is 1 in both x and y we could still choose to
6848 * decompose pme in y instead of x, but we use x for simplicity.
6850 comm->npmedecompdim = 1;
6851 if (dd->dim[0] == YY)
6853 comm->npmenodes_x = 1;
6854 comm->npmenodes_y = comm->npmenodes;
6858 comm->npmenodes_x = comm->npmenodes;
6859 comm->npmenodes_y = 1;
6864 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6865 comm->npmenodes_x, comm->npmenodes_y, 1);
6870 comm->npmedecompdim = 0;
6871 comm->npmenodes_x = 0;
6872 comm->npmenodes_y = 0;
6875 /* Technically we don't need both of these,
6876 * but it simplifies code not having to recalculate it.
6878 *npme_x = comm->npmenodes_x;
6879 *npme_y = comm->npmenodes_y;
6881 snew(comm->slb_frac, DIM);
6882 if (comm->eDLB == edlbNO)
6884 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6885 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6886 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6889 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6891 if (comm->bBondComm || comm->eDLB != edlbNO)
6893 /* Set the bonded communication distance to halfway
6894 * the minimum and the maximum,
6895 * since the extra communication cost is nearly zero.
6897 acs = average_cellsize_min(dd, ddbox);
6898 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6899 if (comm->eDLB != edlbNO)
6901 /* Check if this does not limit the scaling */
6902 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6904 if (!comm->bBondComm)
6906 /* Without bBondComm do not go beyond the n.b. cut-off */
6907 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6908 if (comm->cellsize_limit >= comm->cutoff)
6910 /* We don't loose a lot of efficieny
6911 * when increasing it to the n.b. cut-off.
6912 * It can even be slightly faster, because we need
6913 * less checks for the communication setup.
6915 comm->cutoff_mbody = comm->cutoff;
6918 /* Check if we did not end up below our original limit */
6919 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6921 if (comm->cutoff_mbody > comm->cellsize_limit)
6923 comm->cellsize_limit = comm->cutoff_mbody;
6926 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6931 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6932 "cellsize limit %f\n",
6933 comm->bBondComm, comm->cellsize_limit);
6938 check_dd_restrictions(cr, dd, ir, fplog);
6941 comm->partition_step = INT_MIN;
6944 clear_dd_cycle_counts(dd);
6949 static void set_dlb_limits(gmx_domdec_t *dd)
6954 for (d = 0; d < dd->ndim; d++)
6956 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6957 dd->comm->cellsize_min[dd->dim[d]] =
6958 dd->comm->cellsize_min_dlb[dd->dim[d]];
6963 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6966 gmx_domdec_comm_t *comm;
6976 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);
6979 cellsize_min = comm->cellsize_min[dd->dim[0]];
6980 for (d = 1; d < dd->ndim; d++)
6982 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6985 if (cellsize_min < comm->cellsize_limit*1.05)
6987 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");
6989 /* Change DLB from "auto" to "no". */
6990 comm->eDLB = edlbNO;
6995 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6996 comm->bDynLoadBal = TRUE;
6997 dd->bGridJump = TRUE;
7001 /* We can set the required cell size info here,
7002 * so we do not need to communicate this.
7003 * The grid is completely uniform.
7005 for (d = 0; d < dd->ndim; d++)
7009 comm->load[d].sum_m = comm->load[d].sum;
7011 nc = dd->nc[dd->dim[d]];
7012 for (i = 0; i < nc; i++)
7014 comm->root[d]->cell_f[i] = i/(real)nc;
7017 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7018 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7021 comm->root[d]->cell_f[nc] = 1.0;
7026 static char *init_bLocalCG(gmx_mtop_t *mtop)
7031 ncg = ncg_mtop(mtop);
7032 snew(bLocalCG, ncg);
7033 for (cg = 0; cg < ncg; cg++)
7035 bLocalCG[cg] = FALSE;
7041 void dd_init_bondeds(FILE *fplog,
7042 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7043 gmx_vsite_t *vsite, gmx_constr_t constr,
7044 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7046 gmx_domdec_comm_t *comm;
7050 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7054 if (comm->bBondComm)
7056 /* Communicate atoms beyond the cut-off for bonded interactions */
7059 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7061 comm->bLocalCG = init_bLocalCG(mtop);
7065 /* Only communicate atoms based on cut-off */
7066 comm->cglink = NULL;
7067 comm->bLocalCG = NULL;
7071 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7073 gmx_bool bDynLoadBal, real dlb_scale,
7076 gmx_domdec_comm_t *comm;
7091 fprintf(fplog, "The maximum number of communication pulses is:");
7092 for (d = 0; d < dd->ndim; d++)
7094 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7096 fprintf(fplog, "\n");
7097 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7098 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7099 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7100 for (d = 0; d < DIM; d++)
7104 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7111 comm->cellsize_min_dlb[d]/
7112 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7114 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7117 fprintf(fplog, "\n");
7121 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7122 fprintf(fplog, "The initial number of communication pulses is:");
7123 for (d = 0; d < dd->ndim; d++)
7125 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7127 fprintf(fplog, "\n");
7128 fprintf(fplog, "The initial domain decomposition cell size is:");
7129 for (d = 0; d < DIM; d++)
7133 fprintf(fplog, " %c %.2f nm",
7134 dim2char(d), dd->comm->cellsize_min[d]);
7137 fprintf(fplog, "\n\n");
7140 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7142 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7143 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7144 "non-bonded interactions", "", comm->cutoff);
7148 limit = dd->comm->cellsize_limit;
7152 if (dynamic_dd_box(ddbox, ir))
7154 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7156 limit = dd->comm->cellsize_min[XX];
7157 for (d = 1; d < DIM; d++)
7159 limit = min(limit, dd->comm->cellsize_min[d]);
7163 if (comm->bInterCGBondeds)
7165 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7166 "two-body bonded interactions", "(-rdd)",
7167 max(comm->cutoff, comm->cutoff_mbody));
7168 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7169 "multi-body bonded interactions", "(-rdd)",
7170 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7174 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7175 "virtual site constructions", "(-rcon)", limit);
7177 if (dd->constraint_comm)
7179 sprintf(buf, "atoms separated by up to %d constraints",
7181 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7182 buf, "(-rcon)", limit);
7184 fprintf(fplog, "\n");
7190 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7192 const t_inputrec *ir,
7193 const gmx_ddbox_t *ddbox)
7195 gmx_domdec_comm_t *comm;
7196 int d, dim, npulse, npulse_d_max, npulse_d;
7201 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7203 /* Determine the maximum number of comm. pulses in one dimension */
7205 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7207 /* Determine the maximum required number of grid pulses */
7208 if (comm->cellsize_limit >= comm->cutoff)
7210 /* Only a single pulse is required */
7213 else if (!bNoCutOff && comm->cellsize_limit > 0)
7215 /* We round down slightly here to avoid overhead due to the latency
7216 * of extra communication calls when the cut-off
7217 * would be only slightly longer than the cell size.
7218 * Later cellsize_limit is redetermined,
7219 * so we can not miss interactions due to this rounding.
7221 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7225 /* There is no cell size limit */
7226 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7229 if (!bNoCutOff && npulse > 1)
7231 /* See if we can do with less pulses, based on dlb_scale */
7233 for (d = 0; d < dd->ndim; d++)
7236 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7237 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7238 npulse_d_max = max(npulse_d_max, npulse_d);
7240 npulse = min(npulse, npulse_d_max);
7243 /* This env var can override npulse */
7244 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7251 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7252 for (d = 0; d < dd->ndim; d++)
7254 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7255 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7256 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7257 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7258 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7260 comm->bVacDLBNoLimit = FALSE;
7264 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7265 if (!comm->bVacDLBNoLimit)
7267 comm->cellsize_limit = max(comm->cellsize_limit,
7268 comm->cutoff/comm->maxpulse);
7270 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7271 /* Set the minimum cell size for each DD dimension */
7272 for (d = 0; d < dd->ndim; d++)
7274 if (comm->bVacDLBNoLimit ||
7275 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7277 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7281 comm->cellsize_min_dlb[dd->dim[d]] =
7282 comm->cutoff/comm->cd[d].np_dlb;
7285 if (comm->cutoff_mbody <= 0)
7287 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7289 if (comm->bDynLoadBal)
7295 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7297 /* If each molecule is a single charge group
7298 * or we use domain decomposition for each periodic dimension,
7299 * we do not need to take pbc into account for the bonded interactions.
7301 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7304 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7307 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7308 t_inputrec *ir, t_forcerec *fr,
7311 gmx_domdec_comm_t *comm;
7317 /* Initialize the thread data.
7318 * This can not be done in init_domain_decomposition,
7319 * as the numbers of threads is determined later.
7321 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7324 snew(comm->dth, comm->nth);
7327 if (EEL_PME(ir->coulombtype))
7329 init_ddpme(dd, &comm->ddpme[0], 0);
7330 if (comm->npmedecompdim >= 2)
7332 init_ddpme(dd, &comm->ddpme[1], 1);
7337 comm->npmenodes = 0;
7338 if (dd->pme_nodeid >= 0)
7340 gmx_fatal_collective(FARGS, NULL, dd,
7341 "Can not have separate PME nodes without PME electrostatics");
7347 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7349 if (comm->eDLB != edlbNO)
7351 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7354 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7355 if (comm->eDLB == edlbAUTO)
7359 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7361 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7364 if (ir->ePBC == epbcNONE)
7366 vol_frac = 1 - 1/(double)dd->nnodes;
7371 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7375 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7377 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7379 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7382 static gmx_bool test_dd_cutoff(t_commrec *cr,
7383 t_state *state, t_inputrec *ir,
7394 set_ddbox(dd, FALSE, cr, ir, state->box,
7395 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7399 for (d = 0; d < dd->ndim; d++)
7403 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7404 if (dynamic_dd_box(&ddbox, ir))
7406 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7409 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7411 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7412 dd->comm->cd[d].np_dlb > 0)
7414 if (np > dd->comm->cd[d].np_dlb)
7419 /* If a current local cell size is smaller than the requested
7420 * cut-off, we could still fix it, but this gets very complicated.
7421 * Without fixing here, we might actually need more checks.
7423 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7430 if (dd->comm->eDLB != edlbNO)
7432 /* If DLB is not active yet, we don't need to check the grid jumps.
7433 * Actually we shouldn't, because then the grid jump data is not set.
7435 if (dd->comm->bDynLoadBal &&
7436 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7441 gmx_sumi(1, &LocallyLimited, cr);
7443 if (LocallyLimited > 0)
7452 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7455 gmx_bool bCutoffAllowed;
7457 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7461 cr->dd->comm->cutoff = cutoff_req;
7464 return bCutoffAllowed;
7467 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7469 gmx_domdec_comm_t *comm;
7471 comm = cr->dd->comm;
7473 /* Turn on the DLB limiting (might have been on already) */
7474 comm->bPMELoadBalDLBLimits = TRUE;
7476 /* Change the cut-off limit */
7477 comm->PMELoadBal_max_cutoff = comm->cutoff;
7480 static void merge_cg_buffers(int ncell,
7481 gmx_domdec_comm_dim_t *cd, int pulse,
7483 int *index_gl, int *recv_i,
7484 rvec *cg_cm, rvec *recv_vr,
7486 cginfo_mb_t *cginfo_mb, int *cginfo)
7488 gmx_domdec_ind_t *ind, *ind_p;
7489 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7490 int shift, shift_at;
7492 ind = &cd->ind[pulse];
7494 /* First correct the already stored data */
7495 shift = ind->nrecv[ncell];
7496 for (cell = ncell-1; cell >= 0; cell--)
7498 shift -= ind->nrecv[cell];
7501 /* Move the cg's present from previous grid pulses */
7502 cg0 = ncg_cell[ncell+cell];
7503 cg1 = ncg_cell[ncell+cell+1];
7504 cgindex[cg1+shift] = cgindex[cg1];
7505 for (cg = cg1-1; cg >= cg0; cg--)
7507 index_gl[cg+shift] = index_gl[cg];
7508 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7509 cgindex[cg+shift] = cgindex[cg];
7510 cginfo[cg+shift] = cginfo[cg];
7512 /* Correct the already stored send indices for the shift */
7513 for (p = 1; p <= pulse; p++)
7515 ind_p = &cd->ind[p];
7517 for (c = 0; c < cell; c++)
7519 cg0 += ind_p->nsend[c];
7521 cg1 = cg0 + ind_p->nsend[cell];
7522 for (cg = cg0; cg < cg1; cg++)
7524 ind_p->index[cg] += shift;
7530 /* Merge in the communicated buffers */
7534 for (cell = 0; cell < ncell; cell++)
7536 cg1 = ncg_cell[ncell+cell+1] + shift;
7539 /* Correct the old cg indices */
7540 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7542 cgindex[cg+1] += shift_at;
7545 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7547 /* Copy this charge group from the buffer */
7548 index_gl[cg1] = recv_i[cg0];
7549 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7550 /* Add it to the cgindex */
7551 cg_gl = index_gl[cg1];
7552 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7553 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7554 cgindex[cg1+1] = cgindex[cg1] + nat;
7559 shift += ind->nrecv[cell];
7560 ncg_cell[ncell+cell+1] = cg1;
7564 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7565 int nzone, int cg0, const int *cgindex)
7569 /* Store the atom block boundaries for easy copying of communication buffers
7572 for (zone = 0; zone < nzone; zone++)
7574 for (p = 0; p < cd->np; p++)
7576 cd->ind[p].cell2at0[zone] = cgindex[cg];
7577 cg += cd->ind[p].nrecv[zone];
7578 cd->ind[p].cell2at1[zone] = cgindex[cg];
7583 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7589 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7591 if (!bLocalCG[link->a[i]])
7600 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7602 real c[DIM][4]; /* the corners for the non-bonded communication */
7603 real cr0; /* corner for rounding */
7604 real cr1[4]; /* corners for rounding */
7605 real bc[DIM]; /* corners for bounded communication */
7606 real bcr1; /* corner for rounding for bonded communication */
7609 /* Determine the corners of the domain(s) we are communicating with */
7611 set_dd_corners(const gmx_domdec_t *dd,
7612 int dim0, int dim1, int dim2,
7616 const gmx_domdec_comm_t *comm;
7617 const gmx_domdec_zones_t *zones;
7622 zones = &comm->zones;
7624 /* Keep the compiler happy */
7628 /* The first dimension is equal for all cells */
7629 c->c[0][0] = comm->cell_x0[dim0];
7632 c->bc[0] = c->c[0][0];
7637 /* This cell row is only seen from the first row */
7638 c->c[1][0] = comm->cell_x0[dim1];
7639 /* All rows can see this row */
7640 c->c[1][1] = comm->cell_x0[dim1];
7643 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7646 /* For the multi-body distance we need the maximum */
7647 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7650 /* Set the upper-right corner for rounding */
7651 c->cr0 = comm->cell_x1[dim0];
7656 for (j = 0; j < 4; j++)
7658 c->c[2][j] = comm->cell_x0[dim2];
7662 /* Use the maximum of the i-cells that see a j-cell */
7663 for (i = 0; i < zones->nizone; i++)
7665 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7671 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7677 /* For the multi-body distance we need the maximum */
7678 c->bc[2] = comm->cell_x0[dim2];
7679 for (i = 0; i < 2; i++)
7681 for (j = 0; j < 2; j++)
7683 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7689 /* Set the upper-right corner for rounding */
7690 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7691 * Only cell (0,0,0) can see cell 7 (1,1,1)
7693 c->cr1[0] = comm->cell_x1[dim1];
7694 c->cr1[3] = comm->cell_x1[dim1];
7697 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7700 /* For the multi-body distance we need the maximum */
7701 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7708 /* Determine which cg's we need to send in this pulse from this zone */
7710 get_zone_pulse_cgs(gmx_domdec_t *dd,
7711 int zonei, int zone,
7713 const int *index_gl,
7715 int dim, int dim_ind,
7716 int dim0, int dim1, int dim2,
7717 real r_comm2, real r_bcomm2,
7721 real skew_fac2_d, real skew_fac_01,
7722 rvec *v_d, rvec *v_0, rvec *v_1,
7723 const dd_corners_t *c,
7725 gmx_bool bDistBonded,
7731 gmx_domdec_ind_t *ind,
7732 int **ibuf, int *ibuf_nalloc,
7738 gmx_domdec_comm_t *comm;
7740 gmx_bool bDistMB_pulse;
7742 real r2, rb2, r, tric_sh;
7745 int nsend_z, nsend, nat;
7749 bScrew = (dd->bScrewPBC && dim == XX);
7751 bDistMB_pulse = (bDistMB && bDistBonded);
7757 for (cg = cg0; cg < cg1; cg++)
7761 if (tric_dist[dim_ind] == 0)
7763 /* Rectangular direction, easy */
7764 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7771 r = cg_cm[cg][dim] - c->bc[dim_ind];
7777 /* Rounding gives at most a 16% reduction
7778 * in communicated atoms
7780 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7782 r = cg_cm[cg][dim0] - c->cr0;
7783 /* This is the first dimension, so always r >= 0 */
7790 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7792 r = cg_cm[cg][dim1] - c->cr1[zone];
7799 r = cg_cm[cg][dim1] - c->bcr1;
7809 /* Triclinic direction, more complicated */
7812 /* Rounding, conservative as the skew_fac multiplication
7813 * will slightly underestimate the distance.
7815 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7817 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7818 for (i = dim0+1; i < DIM; i++)
7820 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7822 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7825 rb[dim0] = rn[dim0];
7828 /* Take care that the cell planes along dim0 might not
7829 * be orthogonal to those along dim1 and dim2.
7831 for (i = 1; i <= dim_ind; i++)
7834 if (normal[dim0][dimd] > 0)
7836 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7839 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7844 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7846 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7848 for (i = dim1+1; i < DIM; i++)
7850 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7852 rn[dim1] += tric_sh;
7855 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7856 /* Take care of coupling of the distances
7857 * to the planes along dim0 and dim1 through dim2.
7859 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7860 /* Take care that the cell planes along dim1
7861 * might not be orthogonal to that along dim2.
7863 if (normal[dim1][dim2] > 0)
7865 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7871 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7874 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7875 /* Take care of coupling of the distances
7876 * to the planes along dim0 and dim1 through dim2.
7878 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7879 /* Take care that the cell planes along dim1
7880 * might not be orthogonal to that along dim2.
7882 if (normal[dim1][dim2] > 0)
7884 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7889 /* The distance along the communication direction */
7890 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7892 for (i = dim+1; i < DIM; i++)
7894 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7899 r2 += rn[dim]*rn[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 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7911 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7914 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7915 /* Take care of coupling of the distances
7916 * to the planes along dim0 and dim1 through dim2.
7918 if (dim_ind == 1 && zonei == 1)
7920 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7928 ((bDistMB && rb2 < r_bcomm2) ||
7929 (bDist2B && r2 < r_bcomm2)) &&
7931 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7932 missing_link(comm->cglink, index_gl[cg],
7935 /* Make an index to the local charge groups */
7936 if (nsend+1 > ind->nalloc)
7938 ind->nalloc = over_alloc_large(nsend+1);
7939 srenew(ind->index, ind->nalloc);
7941 if (nsend+1 > *ibuf_nalloc)
7943 *ibuf_nalloc = over_alloc_large(nsend+1);
7944 srenew(*ibuf, *ibuf_nalloc);
7946 ind->index[nsend] = cg;
7947 (*ibuf)[nsend] = index_gl[cg];
7949 vec_rvec_check_alloc(vbuf, nsend+1);
7951 if (dd->ci[dim] == 0)
7953 /* Correct cg_cm for pbc */
7954 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7957 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7958 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7963 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7966 nat += cgindex[cg+1] - cgindex[cg];
7972 *nsend_z_ptr = nsend_z;
7975 static void setup_dd_communication(gmx_domdec_t *dd,
7976 matrix box, gmx_ddbox_t *ddbox,
7977 t_forcerec *fr, t_state *state, rvec **f)
7979 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7980 int nzone, nzone_send, zone, zonei, cg0, cg1;
7981 int c, i, j, cg, cg_gl, nrcg;
7982 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7983 gmx_domdec_comm_t *comm;
7984 gmx_domdec_zones_t *zones;
7985 gmx_domdec_comm_dim_t *cd;
7986 gmx_domdec_ind_t *ind;
7987 cginfo_mb_t *cginfo_mb;
7988 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7989 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7990 dd_corners_t corners;
7992 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7993 real skew_fac2_d, skew_fac_01;
8000 fprintf(debug, "Setting up DD communication\n");
8005 switch (fr->cutoff_scheme)
8014 gmx_incons("unimplemented");
8018 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8020 dim = dd->dim[dim_ind];
8022 /* Check if we need to use triclinic distances */
8023 tric_dist[dim_ind] = 0;
8024 for (i = 0; i <= dim_ind; i++)
8026 if (ddbox->tric_dir[dd->dim[i]])
8028 tric_dist[dim_ind] = 1;
8033 bBondComm = comm->bBondComm;
8035 /* Do we need to determine extra distances for multi-body bondeds? */
8036 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8038 /* Do we need to determine extra distances for only two-body bondeds? */
8039 bDist2B = (bBondComm && !bDistMB);
8041 r_comm2 = sqr(comm->cutoff);
8042 r_bcomm2 = sqr(comm->cutoff_mbody);
8046 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8049 zones = &comm->zones;
8052 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8053 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8055 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8057 /* Triclinic stuff */
8058 normal = ddbox->normal;
8062 v_0 = ddbox->v[dim0];
8063 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8065 /* Determine the coupling coefficient for the distances
8066 * to the cell planes along dim0 and dim1 through dim2.
8067 * This is required for correct rounding.
8070 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8073 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8079 v_1 = ddbox->v[dim1];
8082 zone_cg_range = zones->cg_range;
8083 index_gl = dd->index_gl;
8084 cgindex = dd->cgindex;
8085 cginfo_mb = fr->cginfo_mb;
8087 zone_cg_range[0] = 0;
8088 zone_cg_range[1] = dd->ncg_home;
8089 comm->zone_ncg1[0] = dd->ncg_home;
8090 pos_cg = dd->ncg_home;
8092 nat_tot = dd->nat_home;
8094 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8096 dim = dd->dim[dim_ind];
8097 cd = &comm->cd[dim_ind];
8099 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8101 /* No pbc in this dimension, the first node should not comm. */
8109 v_d = ddbox->v[dim];
8110 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8112 cd->bInPlace = TRUE;
8113 for (p = 0; p < cd->np; p++)
8115 /* Only atoms communicated in the first pulse are used
8116 * for multi-body bonded interactions or for bBondComm.
8118 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8123 for (zone = 0; zone < nzone_send; zone++)
8125 if (tric_dist[dim_ind] && dim_ind > 0)
8127 /* Determine slightly more optimized skew_fac's
8129 * This reduces the number of communicated atoms
8130 * by about 10% for 3D DD of rhombic dodecahedra.
8132 for (dimd = 0; dimd < dim; dimd++)
8134 sf2_round[dimd] = 1;
8135 if (ddbox->tric_dir[dimd])
8137 for (i = dd->dim[dimd]+1; i < DIM; i++)
8139 /* If we are shifted in dimension i
8140 * and the cell plane is tilted forward
8141 * in dimension i, skip this coupling.
8143 if (!(zones->shift[nzone+zone][i] &&
8144 ddbox->v[dimd][i][dimd] >= 0))
8147 sqr(ddbox->v[dimd][i][dimd]);
8150 sf2_round[dimd] = 1/sf2_round[dimd];
8155 zonei = zone_perm[dim_ind][zone];
8158 /* Here we permutate the zones to obtain a convenient order
8159 * for neighbor searching
8161 cg0 = zone_cg_range[zonei];
8162 cg1 = zone_cg_range[zonei+1];
8166 /* Look only at the cg's received in the previous grid pulse
8168 cg1 = zone_cg_range[nzone+zone+1];
8169 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8172 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8173 for (th = 0; th < comm->nth; th++)
8175 gmx_domdec_ind_t *ind_p;
8176 int **ibuf_p, *ibuf_nalloc_p;
8178 int *nsend_p, *nat_p;
8184 /* Thread 0 writes in the comm buffers */
8186 ibuf_p = &comm->buf_int;
8187 ibuf_nalloc_p = &comm->nalloc_int;
8188 vbuf_p = &comm->vbuf;
8191 nsend_zone_p = &ind->nsend[zone];
8195 /* Other threads write into temp buffers */
8196 ind_p = &comm->dth[th].ind;
8197 ibuf_p = &comm->dth[th].ibuf;
8198 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8199 vbuf_p = &comm->dth[th].vbuf;
8200 nsend_p = &comm->dth[th].nsend;
8201 nat_p = &comm->dth[th].nat;
8202 nsend_zone_p = &comm->dth[th].nsend_zone;
8204 comm->dth[th].nsend = 0;
8205 comm->dth[th].nat = 0;
8206 comm->dth[th].nsend_zone = 0;
8216 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8217 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8220 /* Get the cg's for this pulse in this zone */
8221 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8223 dim, dim_ind, dim0, dim1, dim2,
8226 normal, skew_fac2_d, skew_fac_01,
8227 v_d, v_0, v_1, &corners, sf2_round,
8228 bDistBonded, bBondComm,
8232 ibuf_p, ibuf_nalloc_p,
8238 /* Append data of threads>=1 to the communication buffers */
8239 for (th = 1; th < comm->nth; th++)
8241 dd_comm_setup_work_t *dth;
8244 dth = &comm->dth[th];
8246 ns1 = nsend + dth->nsend_zone;
8247 if (ns1 > ind->nalloc)
8249 ind->nalloc = over_alloc_dd(ns1);
8250 srenew(ind->index, ind->nalloc);
8252 if (ns1 > comm->nalloc_int)
8254 comm->nalloc_int = over_alloc_dd(ns1);
8255 srenew(comm->buf_int, comm->nalloc_int);
8257 if (ns1 > comm->vbuf.nalloc)
8259 comm->vbuf.nalloc = over_alloc_dd(ns1);
8260 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8263 for (i = 0; i < dth->nsend_zone; i++)
8265 ind->index[nsend] = dth->ind.index[i];
8266 comm->buf_int[nsend] = dth->ibuf[i];
8267 copy_rvec(dth->vbuf.v[i],
8268 comm->vbuf.v[nsend]);
8272 ind->nsend[zone] += dth->nsend_zone;
8275 /* Clear the counts in case we do not have pbc */
8276 for (zone = nzone_send; zone < nzone; zone++)
8278 ind->nsend[zone] = 0;
8280 ind->nsend[nzone] = nsend;
8281 ind->nsend[nzone+1] = nat;
8282 /* Communicate the number of cg's and atoms to receive */
8283 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8284 ind->nsend, nzone+2,
8285 ind->nrecv, nzone+2);
8287 /* The rvec buffer is also required for atom buffers of size nsend
8288 * in dd_move_x and dd_move_f.
8290 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8294 /* We can receive in place if only the last zone is not empty */
8295 for (zone = 0; zone < nzone-1; zone++)
8297 if (ind->nrecv[zone] > 0)
8299 cd->bInPlace = FALSE;
8304 /* The int buffer is only required here for the cg indices */
8305 if (ind->nrecv[nzone] > comm->nalloc_int2)
8307 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8308 srenew(comm->buf_int2, comm->nalloc_int2);
8310 /* The rvec buffer is also required for atom buffers
8311 * of size nrecv in dd_move_x and dd_move_f.
8313 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8314 vec_rvec_check_alloc(&comm->vbuf2, i);
8318 /* Make space for the global cg indices */
8319 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8320 || dd->cg_nalloc == 0)
8322 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8323 srenew(index_gl, dd->cg_nalloc);
8324 srenew(cgindex, dd->cg_nalloc+1);
8326 /* Communicate the global cg indices */
8329 recv_i = index_gl + pos_cg;
8333 recv_i = comm->buf_int2;
8335 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8336 comm->buf_int, nsend,
8337 recv_i, ind->nrecv[nzone]);
8339 /* Make space for cg_cm */
8340 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8341 if (fr->cutoff_scheme == ecutsGROUP)
8349 /* Communicate cg_cm */
8352 recv_vr = cg_cm + pos_cg;
8356 recv_vr = comm->vbuf2.v;
8358 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8359 comm->vbuf.v, nsend,
8360 recv_vr, ind->nrecv[nzone]);
8362 /* Make the charge group index */
8365 zone = (p == 0 ? 0 : nzone - 1);
8366 while (zone < nzone)
8368 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8370 cg_gl = index_gl[pos_cg];
8371 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8372 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8373 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8376 /* Update the charge group presence,
8377 * so we can use it in the next pass of the loop.
8379 comm->bLocalCG[cg_gl] = TRUE;
8385 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8388 zone_cg_range[nzone+zone] = pos_cg;
8393 /* This part of the code is never executed with bBondComm. */
8394 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8395 index_gl, recv_i, cg_cm, recv_vr,
8396 cgindex, fr->cginfo_mb, fr->cginfo);
8397 pos_cg += ind->nrecv[nzone];
8399 nat_tot += ind->nrecv[nzone+1];
8403 /* Store the atom block for easy copying of communication buffers */
8404 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8408 dd->index_gl = index_gl;
8409 dd->cgindex = cgindex;
8411 dd->ncg_tot = zone_cg_range[zones->n];
8412 dd->nat_tot = nat_tot;
8413 comm->nat[ddnatHOME] = dd->nat_home;
8414 for (i = ddnatZONE; i < ddnatNR; i++)
8416 comm->nat[i] = dd->nat_tot;
8421 /* We don't need to update cginfo, since that was alrady done above.
8422 * So we pass NULL for the forcerec.
8424 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8425 NULL, comm->bLocalCG);
8430 fprintf(debug, "Finished setting up DD communication, zones:");
8431 for (c = 0; c < zones->n; c++)
8433 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8435 fprintf(debug, "\n");
8439 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8443 for (c = 0; c < zones->nizone; c++)
8445 zones->izone[c].cg1 = zones->cg_range[c+1];
8446 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8447 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8451 static void set_zones_size(gmx_domdec_t *dd,
8452 matrix box, const gmx_ddbox_t *ddbox,
8453 int zone_start, int zone_end)
8455 gmx_domdec_comm_t *comm;
8456 gmx_domdec_zones_t *zones;
8458 int z, zi, zj0, zj1, d, dim;
8461 real size_j, add_tric;
8466 zones = &comm->zones;
8468 /* Do we need to determine extra distances for multi-body bondeds? */
8469 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8471 for (z = zone_start; z < zone_end; z++)
8473 /* Copy cell limits to zone limits.
8474 * Valid for non-DD dims and non-shifted dims.
8476 copy_rvec(comm->cell_x0, zones->size[z].x0);
8477 copy_rvec(comm->cell_x1, zones->size[z].x1);
8480 for (d = 0; d < dd->ndim; d++)
8484 for (z = 0; z < zones->n; z++)
8486 /* With a staggered grid we have different sizes
8487 * for non-shifted dimensions.
8489 if (dd->bGridJump && zones->shift[z][dim] == 0)
8493 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8494 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8498 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8499 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8505 rcmbs = comm->cutoff_mbody;
8506 if (ddbox->tric_dir[dim])
8508 rcs /= ddbox->skew_fac[dim];
8509 rcmbs /= ddbox->skew_fac[dim];
8512 /* Set the lower limit for the shifted zone dimensions */
8513 for (z = zone_start; z < zone_end; z++)
8515 if (zones->shift[z][dim] > 0)
8518 if (!dd->bGridJump || d == 0)
8520 zones->size[z].x0[dim] = comm->cell_x1[dim];
8521 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8525 /* Here we take the lower limit of the zone from
8526 * the lowest domain of the zone below.
8530 zones->size[z].x0[dim] =
8531 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8537 zones->size[z].x0[dim] =
8538 zones->size[zone_perm[2][z-4]].x0[dim];
8542 zones->size[z].x0[dim] =
8543 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8546 /* A temporary limit, is updated below */
8547 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8551 for (zi = 0; zi < zones->nizone; zi++)
8553 if (zones->shift[zi][dim] == 0)
8555 /* This takes the whole zone into account.
8556 * With multiple pulses this will lead
8557 * to a larger zone then strictly necessary.
8559 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8560 zones->size[zi].x1[dim]+rcmbs);
8568 /* Loop over the i-zones to set the upper limit of each
8571 for (zi = 0; zi < zones->nizone; zi++)
8573 if (zones->shift[zi][dim] == 0)
8575 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8577 if (zones->shift[z][dim] > 0)
8579 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8580 zones->size[zi].x1[dim]+rcs);
8587 for (z = zone_start; z < zone_end; z++)
8589 /* Initialization only required to keep the compiler happy */
8590 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8593 /* To determine the bounding box for a zone we need to find
8594 * the extreme corners of 4, 2 or 1 corners.
8596 nc = 1 << (ddbox->npbcdim - 1);
8598 for (c = 0; c < nc; c++)
8600 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8604 corner[YY] = zones->size[z].x0[YY];
8608 corner[YY] = zones->size[z].x1[YY];
8612 corner[ZZ] = zones->size[z].x0[ZZ];
8616 corner[ZZ] = zones->size[z].x1[ZZ];
8618 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8620 /* With 1D domain decomposition the cg's are not in
8621 * the triclinic box, but triclinic x-y and rectangular y-z.
8622 * Shift y back, so it will later end up at 0.
8624 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8626 /* Apply the triclinic couplings */
8627 for (i = YY; i < ddbox->npbcdim; i++)
8629 for (j = XX; j < i; j++)
8631 corner[j] += corner[i]*box[i][j]/box[i][i];
8636 copy_rvec(corner, corner_min);
8637 copy_rvec(corner, corner_max);
8641 for (i = 0; i < DIM; i++)
8643 corner_min[i] = min(corner_min[i], corner[i]);
8644 corner_max[i] = max(corner_max[i], corner[i]);
8648 /* Copy the extreme cornes without offset along x */
8649 for (i = 0; i < DIM; i++)
8651 zones->size[z].bb_x0[i] = corner_min[i];
8652 zones->size[z].bb_x1[i] = corner_max[i];
8654 /* Add the offset along x */
8655 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8656 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8659 if (zone_start == 0)
8662 for (dim = 0; dim < DIM; dim++)
8664 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8666 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8671 for (z = zone_start; z < zone_end; z++)
8673 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8675 zones->size[z].x0[XX], zones->size[z].x1[XX],
8676 zones->size[z].x0[YY], zones->size[z].x1[YY],
8677 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8678 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8680 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8681 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8682 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8687 static int comp_cgsort(const void *a, const void *b)
8691 gmx_cgsort_t *cga, *cgb;
8692 cga = (gmx_cgsort_t *)a;
8693 cgb = (gmx_cgsort_t *)b;
8695 comp = cga->nsc - cgb->nsc;
8698 comp = cga->ind_gl - cgb->ind_gl;
8704 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8709 /* Order the data */
8710 for (i = 0; i < n; i++)
8712 buf[i] = a[sort[i].ind];
8715 /* Copy back to the original array */
8716 for (i = 0; i < n; i++)
8722 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8727 /* Order the data */
8728 for (i = 0; i < n; i++)
8730 copy_rvec(v[sort[i].ind], buf[i]);
8733 /* Copy back to the original array */
8734 for (i = 0; i < n; i++)
8736 copy_rvec(buf[i], v[i]);
8740 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8743 int a, atot, cg, cg0, cg1, i;
8745 if (cgindex == NULL)
8747 /* Avoid the useless loop of the atoms within a cg */
8748 order_vec_cg(ncg, sort, v, buf);
8753 /* Order the data */
8755 for (cg = 0; cg < ncg; cg++)
8757 cg0 = cgindex[sort[cg].ind];
8758 cg1 = cgindex[sort[cg].ind+1];
8759 for (i = cg0; i < cg1; i++)
8761 copy_rvec(v[i], buf[a]);
8767 /* Copy back to the original array */
8768 for (a = 0; a < atot; a++)
8770 copy_rvec(buf[a], v[a]);
8774 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8775 int nsort_new, gmx_cgsort_t *sort_new,
8776 gmx_cgsort_t *sort1)
8780 /* The new indices are not very ordered, so we qsort them */
8781 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8783 /* sort2 is already ordered, so now we can merge the two arrays */
8787 while (i2 < nsort2 || i_new < nsort_new)
8791 sort1[i1++] = sort_new[i_new++];
8793 else if (i_new == nsort_new)
8795 sort1[i1++] = sort2[i2++];
8797 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8798 (sort2[i2].nsc == sort_new[i_new].nsc &&
8799 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8801 sort1[i1++] = sort2[i2++];
8805 sort1[i1++] = sort_new[i_new++];
8810 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8812 gmx_domdec_sort_t *sort;
8813 gmx_cgsort_t *cgsort, *sort_i;
8814 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8815 int sort_last, sort_skip;
8817 sort = dd->comm->sort;
8819 a = fr->ns.grid->cell_index;
8821 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8823 if (ncg_home_old >= 0)
8825 /* The charge groups that remained in the same ns grid cell
8826 * are completely ordered. So we can sort efficiently by sorting
8827 * the charge groups that did move into the stationary list.
8832 for (i = 0; i < dd->ncg_home; i++)
8834 /* Check if this cg did not move to another node */
8837 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8839 /* This cg is new on this node or moved ns grid cell */
8840 if (nsort_new >= sort->sort_new_nalloc)
8842 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8843 srenew(sort->sort_new, sort->sort_new_nalloc);
8845 sort_i = &(sort->sort_new[nsort_new++]);
8849 /* This cg did not move */
8850 sort_i = &(sort->sort2[nsort2++]);
8852 /* Sort on the ns grid cell indices
8853 * and the global topology index.
8854 * index_gl is irrelevant with cell ns,
8855 * but we set it here anyhow to avoid a conditional.
8858 sort_i->ind_gl = dd->index_gl[i];
8865 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8868 /* Sort efficiently */
8869 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8874 cgsort = sort->sort;
8876 for (i = 0; i < dd->ncg_home; i++)
8878 /* Sort on the ns grid cell indices
8879 * and the global topology index
8881 cgsort[i].nsc = a[i];
8882 cgsort[i].ind_gl = dd->index_gl[i];
8884 if (cgsort[i].nsc < moved)
8891 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8893 /* Determine the order of the charge groups using qsort */
8894 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8900 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8903 int ncg_new, i, *a, na;
8905 sort = dd->comm->sort->sort;
8907 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8910 for (i = 0; i < na; i++)
8914 sort[ncg_new].ind = a[i];
8922 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8923 rvec *cgcm, t_forcerec *fr, t_state *state,
8926 gmx_domdec_sort_t *sort;
8927 gmx_cgsort_t *cgsort, *sort_i;
8929 int ncg_new, i, *ibuf, cgsize;
8932 sort = dd->comm->sort;
8934 if (dd->ncg_home > sort->sort_nalloc)
8936 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8937 srenew(sort->sort, sort->sort_nalloc);
8938 srenew(sort->sort2, sort->sort_nalloc);
8940 cgsort = sort->sort;
8942 switch (fr->cutoff_scheme)
8945 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8948 ncg_new = dd_sort_order_nbnxn(dd, fr);
8951 gmx_incons("unimplemented");
8955 /* We alloc with the old size, since cgindex is still old */
8956 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8957 vbuf = dd->comm->vbuf.v;
8961 cgindex = dd->cgindex;
8968 /* Remove the charge groups which are no longer at home here */
8969 dd->ncg_home = ncg_new;
8972 fprintf(debug, "Set the new home charge group count to %d\n",
8976 /* Reorder the state */
8977 for (i = 0; i < estNR; i++)
8979 if (EST_DISTR(i) && (state->flags & (1<<i)))
8984 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8987 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8990 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8993 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8997 case estDISRE_INITF:
8998 case estDISRE_RM3TAV:
8999 case estORIRE_INITF:
9001 /* No ordering required */
9004 gmx_incons("Unknown state entry encountered in dd_sort_state");
9009 if (fr->cutoff_scheme == ecutsGROUP)
9012 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9015 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9017 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9018 srenew(sort->ibuf, sort->ibuf_nalloc);
9021 /* Reorder the global cg index */
9022 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9023 /* Reorder the cginfo */
9024 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9025 /* Rebuild the local cg index */
9029 for (i = 0; i < dd->ncg_home; i++)
9031 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9032 ibuf[i+1] = ibuf[i] + cgsize;
9034 for (i = 0; i < dd->ncg_home+1; i++)
9036 dd->cgindex[i] = ibuf[i];
9041 for (i = 0; i < dd->ncg_home+1; i++)
9046 /* Set the home atom number */
9047 dd->nat_home = dd->cgindex[dd->ncg_home];
9049 if (fr->cutoff_scheme == ecutsVERLET)
9051 /* The atoms are now exactly in grid order, update the grid order */
9052 nbnxn_set_atomorder(fr->nbv->nbs);
9056 /* Copy the sorted ns cell indices back to the ns grid struct */
9057 for (i = 0; i < dd->ncg_home; i++)
9059 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9061 fr->ns.grid->nr = dd->ncg_home;
9065 static void add_dd_statistics(gmx_domdec_t *dd)
9067 gmx_domdec_comm_t *comm;
9072 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9074 comm->sum_nat[ddnat-ddnatZONE] +=
9075 comm->nat[ddnat] - comm->nat[ddnat-1];
9080 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9082 gmx_domdec_comm_t *comm;
9087 /* Reset all the statistics and counters for total run counting */
9088 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9090 comm->sum_nat[ddnat-ddnatZONE] = 0;
9094 comm->load_step = 0;
9097 clear_ivec(comm->load_lim);
9102 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9104 gmx_domdec_comm_t *comm;
9108 comm = cr->dd->comm;
9110 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9117 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");
9119 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9121 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9126 " av. #atoms communicated per step for force: %d x %.1f\n",
9130 if (cr->dd->vsite_comm)
9133 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9134 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9139 if (cr->dd->constraint_comm)
9142 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9143 1 + ir->nLincsIter, av);
9147 gmx_incons(" Unknown type for DD statistics");
9150 fprintf(fplog, "\n");
9152 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9154 print_dd_load_av(fplog, cr->dd);
9158 void dd_partition_system(FILE *fplog,
9159 gmx_large_int_t step,
9161 gmx_bool bMasterState,
9163 t_state *state_global,
9164 gmx_mtop_t *top_global,
9166 t_state *state_local,
9169 gmx_localtop_t *top_local,
9172 gmx_shellfc_t shellfc,
9173 gmx_constr_t constr,
9175 gmx_wallcycle_t wcycle,
9179 gmx_domdec_comm_t *comm;
9180 gmx_ddbox_t ddbox = {0};
9182 gmx_large_int_t step_pcoupl;
9183 rvec cell_ns_x0, cell_ns_x1;
9184 int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9185 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9186 gmx_bool bRedist, bSortCG, bResortAll;
9187 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9194 bBoxChanged = (bMasterState || DEFORM(*ir));
9195 if (ir->epc != epcNO)
9197 /* With nstpcouple > 1 pressure coupling happens.
9198 * one step after calculating the pressure.
9199 * Box scaling happens at the end of the MD step,
9200 * after the DD partitioning.
9201 * We therefore have to do DLB in the first partitioning
9202 * after an MD step where P-coupling occured.
9203 * We need to determine the last step in which p-coupling occurred.
9204 * MRS -- need to validate this for vv?
9209 step_pcoupl = step - 1;
9213 step_pcoupl = ((step - 1)/n)*n + 1;
9215 if (step_pcoupl >= comm->partition_step)
9221 bNStGlobalComm = (step % nstglobalcomm == 0);
9223 if (!comm->bDynLoadBal)
9229 /* Should we do dynamic load balacing this step?
9230 * Since it requires (possibly expensive) global communication,
9231 * we might want to do DLB less frequently.
9233 if (bBoxChanged || ir->epc != epcNO)
9235 bDoDLB = bBoxChanged;
9239 bDoDLB = bNStGlobalComm;
9243 /* Check if we have recorded loads on the nodes */
9244 if (comm->bRecordLoad && dd_load_count(comm))
9246 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9248 /* Check if we should use DLB at the second partitioning
9249 * and every 100 partitionings,
9250 * so the extra communication cost is negligible.
9252 n = max(100, nstglobalcomm);
9253 bCheckDLB = (comm->n_load_collect == 0 ||
9254 comm->n_load_have % n == n-1);
9261 /* Print load every nstlog, first and last step to the log file */
9262 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9263 comm->n_load_collect == 0 ||
9265 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9267 /* Avoid extra communication due to verbose screen output
9268 * when nstglobalcomm is set.
9270 if (bDoDLB || bLogLoad || bCheckDLB ||
9271 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9273 get_load_distribution(dd, wcycle);
9278 dd_print_load(fplog, dd, step-1);
9282 dd_print_load_verbose(dd);
9285 comm->n_load_collect++;
9289 /* Since the timings are node dependent, the master decides */
9293 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9296 fprintf(debug, "step %s, imb loss %f\n",
9297 gmx_step_str(step, sbuf),
9298 dd_force_imb_perf_loss(dd));
9301 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9304 turn_on_dlb(fplog, cr, step);
9309 comm->n_load_have++;
9312 cgs_gl = &comm->cgs_gl;
9317 /* Clear the old state */
9318 clear_dd_indices(dd, 0, 0);
9320 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9321 TRUE, cgs_gl, state_global->x, &ddbox);
9323 get_cg_distribution(fplog, step, dd, cgs_gl,
9324 state_global->box, &ddbox, state_global->x);
9326 dd_distribute_state(dd, cgs_gl,
9327 state_global, state_local, f);
9329 dd_make_local_cgs(dd, &top_local->cgs);
9331 /* Ensure that we have space for the new distribution */
9332 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9334 if (fr->cutoff_scheme == ecutsGROUP)
9336 calc_cgcm(fplog, 0, dd->ncg_home,
9337 &top_local->cgs, state_local->x, fr->cg_cm);
9340 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9342 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9346 else if (state_local->ddp_count != dd->ddp_count)
9348 if (state_local->ddp_count > dd->ddp_count)
9350 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9353 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9355 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);
9358 /* Clear the old state */
9359 clear_dd_indices(dd, 0, 0);
9361 /* Build the new indices */
9362 rebuild_cgindex(dd, cgs_gl->index, state_local);
9363 make_dd_indices(dd, cgs_gl->index, 0);
9365 if (fr->cutoff_scheme == ecutsGROUP)
9367 /* Redetermine the cg COMs */
9368 calc_cgcm(fplog, 0, dd->ncg_home,
9369 &top_local->cgs, state_local->x, fr->cg_cm);
9372 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9374 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9376 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9377 TRUE, &top_local->cgs, state_local->x, &ddbox);
9379 bRedist = comm->bDynLoadBal;
9383 /* We have the full state, only redistribute the cgs */
9385 /* Clear the non-home indices */
9386 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9388 /* Avoid global communication for dim's without pbc and -gcom */
9389 if (!bNStGlobalComm)
9391 copy_rvec(comm->box0, ddbox.box0 );
9392 copy_rvec(comm->box_size, ddbox.box_size);
9394 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9395 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9400 /* For dim's without pbc and -gcom */
9401 copy_rvec(ddbox.box0, comm->box0 );
9402 copy_rvec(ddbox.box_size, comm->box_size);
9404 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9407 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9409 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9412 /* Check if we should sort the charge groups */
9413 if (comm->nstSortCG > 0)
9415 bSortCG = (bMasterState ||
9416 (bRedist && (step % comm->nstSortCG == 0)));
9423 ncg_home_old = dd->ncg_home;
9428 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9430 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9431 state_local, f, fr, mdatoms,
9432 !bSortCG, nrnb, &cg0, &ncg_moved);
9434 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9437 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9439 &comm->cell_x0, &comm->cell_x1,
9440 dd->ncg_home, fr->cg_cm,
9441 cell_ns_x0, cell_ns_x1, &grid_density);
9445 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9448 switch (fr->cutoff_scheme)
9451 copy_ivec(fr->ns.grid->n, ncells_old);
9452 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9453 state_local->box, cell_ns_x0, cell_ns_x1,
9454 fr->rlistlong, grid_density);
9457 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9460 gmx_incons("unimplemented");
9462 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9463 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9467 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9469 /* Sort the state on charge group position.
9470 * This enables exact restarts from this step.
9471 * It also improves performance by about 15% with larger numbers
9472 * of atoms per node.
9475 /* Fill the ns grid with the home cell,
9476 * so we can sort with the indices.
9478 set_zones_ncg_home(dd);
9480 switch (fr->cutoff_scheme)
9483 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9485 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9487 comm->zones.size[0].bb_x0,
9488 comm->zones.size[0].bb_x1,
9490 comm->zones.dens_zone0,
9493 ncg_moved, bRedist ? comm->moved : NULL,
9494 fr->nbv->grp[eintLocal].kernel_type,
9495 fr->nbv->grp[eintLocal].nbat);
9497 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9500 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9501 0, dd->ncg_home, fr->cg_cm);
9503 copy_ivec(fr->ns.grid->n, ncells_new);
9506 gmx_incons("unimplemented");
9509 bResortAll = bMasterState;
9511 /* Check if we can user the old order and ns grid cell indices
9512 * of the charge groups to sort the charge groups efficiently.
9514 if (ncells_new[XX] != ncells_old[XX] ||
9515 ncells_new[YY] != ncells_old[YY] ||
9516 ncells_new[ZZ] != ncells_old[ZZ])
9523 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9524 gmx_step_str(step, sbuf), dd->ncg_home);
9526 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9527 bResortAll ? -1 : ncg_home_old);
9528 /* Rebuild all the indices */
9530 ga2la_clear(dd->ga2la);
9532 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9535 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9537 /* Setup up the communication and communicate the coordinates */
9538 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9540 /* Set the indices */
9541 make_dd_indices(dd, cgs_gl->index, cg0);
9543 /* Set the charge group boundaries for neighbor searching */
9544 set_cg_boundaries(&comm->zones);
9546 if (fr->cutoff_scheme == ecutsVERLET)
9548 set_zones_size(dd, state_local->box, &ddbox,
9549 bSortCG ? 1 : 0, comm->zones.n);
9552 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9555 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9556 -1,state_local->x,state_local->box);
9559 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9561 /* Extract a local topology from the global topology */
9562 for (i = 0; i < dd->ndim; i++)
9564 np[dd->dim[i]] = comm->cd[i].np;
9566 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9567 comm->cellsize_min, np,
9569 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9570 vsite, top_global, top_local);
9572 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9574 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9576 /* Set up the special atom communication */
9577 n = comm->nat[ddnatZONE];
9578 for (i = ddnatZONE+1; i < ddnatNR; i++)
9583 if (vsite && vsite->n_intercg_vsite)
9585 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9589 if (dd->bInterCGcons || dd->bInterCGsettles)
9591 /* Only for inter-cg constraints we need special code */
9592 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9593 constr, ir->nProjOrder,
9594 top_local->idef.il);
9598 gmx_incons("Unknown special atom type setup");
9603 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9605 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9607 /* Make space for the extra coordinates for virtual site
9608 * or constraint communication.
9610 state_local->natoms = comm->nat[ddnatNR-1];
9611 if (state_local->natoms > state_local->nalloc)
9613 dd_realloc_state(state_local, f, state_local->natoms);
9616 if (fr->bF_NoVirSum)
9618 if (vsite && vsite->n_intercg_vsite)
9620 nat_f_novirsum = comm->nat[ddnatVSITE];
9624 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9626 nat_f_novirsum = dd->nat_tot;
9630 nat_f_novirsum = dd->nat_home;
9639 /* Set the number of atoms required for the force calculation.
9640 * Forces need to be constrained when using a twin-range setup
9641 * or with energy minimization. For simple simulations we could
9642 * avoid some allocation, zeroing and copying, but this is
9643 * probably not worth the complications ande checking.
9645 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9646 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9648 /* We make the all mdatoms up to nat_tot_con.
9649 * We could save some work by only setting invmass
9650 * between nat_tot and nat_tot_con.
9652 /* This call also sets the new number of home particles to dd->nat_home */
9653 atoms2md(top_global, ir,
9654 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9656 /* Now we have the charges we can sort the FE interactions */
9657 dd_sort_local_top(dd, mdatoms, top_local);
9661 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9662 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9667 /* Make the local shell stuff, currently no communication is done */
9668 make_local_shells(cr, mdatoms, shellfc);
9671 if (ir->implicit_solvent)
9673 make_local_gb(cr, fr->born, ir->gb_algorithm);
9676 init_bonded_thread_force_reduction(fr, &top_local->idef);
9678 if (!(cr->duty & DUTY_PME))
9680 /* Send the charges to our PME only node */
9681 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9682 mdatoms->chargeA, mdatoms->chargeB,
9683 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9688 set_constraints(constr, top_local, ir, mdatoms, cr);
9691 if (ir->ePull != epullNO)
9693 /* Update the local pull groups */
9694 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9699 /* Update the local rotation groups */
9700 dd_make_local_rotation_groups(dd, ir->rot);
9704 add_dd_statistics(dd);
9706 /* Make sure we only count the cycles for this DD partitioning */
9707 clear_dd_cycle_counts(dd);
9709 /* Because the order of the atoms might have changed since
9710 * the last vsite construction, we need to communicate the constructing
9711 * atom coordinates again (for spreading the forces this MD step).
9713 dd_move_x_vsites(dd, state_local->box, state_local->x);
9715 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9717 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9719 dd_move_x(dd, state_local->box, state_local->x);
9720 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9721 -1, state_local->x, state_local->box);
9724 /* Store the partitioning step */
9725 comm->partition_step = step;
9727 /* Increase the DD partitioning counter */
9729 /* The state currently matches this DD partitioning count, store it */
9730 state_local->ddp_count = dd->ddp_count;
9733 /* The DD master node knows the complete cg distribution,
9734 * store the count so we can possibly skip the cg info communication.
9736 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9739 if (comm->DD_debug > 0)
9741 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9742 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9743 "after partitioning");