2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/utility/smalloc.h"
49 #include "gmx_fatal.h"
50 #include "gmx_fatal_collective.h"
53 #include "domdec_network.h"
56 #include "chargegroup.h"
65 #include "mtop_util.h"
66 #include "gmx_ga2la.h"
68 #include "nbnxn_search.h"
70 #include "gmx_omp_nthreads.h"
71 #include "gpu_utils.h"
73 #include "gromacs/fileio/futil.h"
74 #include "gromacs/fileio/gmxfio.h"
75 #include "gromacs/fileio/pdbio.h"
76 #include "gromacs/timing/wallcycle.h"
77 #include "gromacs/utility/gmxmpi.h"
78 #include "gromacs/swap/swapcoords.h"
79 #include "gromacs/utility/qsort_threadsafe.h"
80 #include "gromacs/pulling/pull.h"
81 #include "gromacs/pulling/pull_rotation.h"
82 #include "gromacs/imd/imd.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 /* Is eDLB=edlbAUTO locked such that we currently can't turn it on? */
273 gmx_bool bDLB_locked;
274 /* Are we actually using DLB? */
275 gmx_bool bDynLoadBal;
277 /* Cell sizes for static load balancing, first index cartesian */
280 /* The width of the communicated boundaries */
283 /* The minimum cell size (including triclinic correction) */
285 /* For dlb, for use with edlbAUTO */
286 rvec cellsize_min_dlb;
287 /* The lower limit for the DD cell size with DLB */
289 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
290 gmx_bool bVacDLBNoLimit;
292 /* With PME load balancing we set limits on DLB */
293 gmx_bool bPMELoadBalDLBLimits;
294 /* DLB needs to take into account that we want to allow this maximum
295 * cut-off (for PME load balancing), this could limit cell boundaries.
297 real PMELoadBal_max_cutoff;
299 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
301 /* box0 and box_size are required with dim's without pbc and -gcom */
305 /* The cell boundaries */
309 /* The old location of the cell boundaries, to check cg displacements */
313 /* The communication setup and charge group boundaries for the zones */
314 gmx_domdec_zones_t zones;
316 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
317 * cell boundaries of neighboring cells for dynamic load balancing.
319 gmx_ddzone_t zone_d1[2];
320 gmx_ddzone_t zone_d2[2][2];
322 /* The coordinate/force communication setup and indices */
323 gmx_domdec_comm_dim_t cd[DIM];
324 /* The maximum number of cells to communicate with in one dimension */
327 /* Which cg distribution is stored on the master node */
328 int master_cg_ddp_count;
330 /* The number of cg's received from the direct neighbors */
331 int zone_ncg1[DD_MAXZONE];
333 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
336 /* Array for signalling if atoms have moved to another domain */
340 /* Communication buffer for general use */
344 /* Communication buffer for general use */
347 /* Temporary storage for thread parallel communication setup */
349 dd_comm_setup_work_t *dth;
351 /* Communication buffers only used with multiple grid pulses */
356 /* Communication buffers for local redistribution */
358 int cggl_flag_nalloc[DIM*2];
360 int cgcm_state_nalloc[DIM*2];
362 /* Cell sizes for dynamic load balancing */
363 gmx_domdec_root_t **root;
367 real cell_f_max0[DIM];
368 real cell_f_min1[DIM];
370 /* Stuff for load communication */
371 gmx_bool bRecordLoad;
372 gmx_domdec_load_t *load;
373 int nrank_gpu_shared;
375 MPI_Comm *mpi_comm_load;
376 MPI_Comm mpi_comm_gpu_shared;
379 /* Maximum DLB scaling per load balancing step in percent */
383 float cycl[ddCyclNr];
384 int cycl_n[ddCyclNr];
385 float cycl_max[ddCyclNr];
386 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
390 /* How many times have did we have load measurements */
392 /* How many times have we collected the load measurements */
396 double sum_nat[ddnatNR-ddnatZONE];
406 /* The last partition step */
407 gmx_int64_t partition_step;
415 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
418 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
419 #define DD_FLAG_NRCG 65535
420 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
421 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
423 /* Zone permutation required to obtain consecutive charge groups
424 * for neighbor searching.
426 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
428 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
429 * components see only j zones with that component 0.
432 /* The DD zone order */
433 static const ivec dd_zo[DD_MAXZONE] =
434 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
439 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
444 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
449 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
451 /* Factors used to avoid problems due to rounding issues */
452 #define DD_CELL_MARGIN 1.0001
453 #define DD_CELL_MARGIN2 1.00005
454 /* Factor to account for pressure scaling during nstlist steps */
455 #define DD_PRES_SCALE_MARGIN 1.02
457 /* Turn on DLB when the load imbalance causes this amount of total loss.
458 * There is a bit of overhead with DLB and it's difficult to achieve
459 * a load imbalance of less than 2% with DLB.
461 #define DD_PERF_LOSS_DLB_ON 0.02
463 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
464 #define DD_PERF_LOSS_WARN 0.05
466 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
468 /* Use separate MPI send and receive commands
469 * when nnodes <= GMX_DD_NNODES_SENDRECV.
470 * This saves memory (and some copying for small nnodes).
471 * For high parallelization scatter and gather calls are used.
473 #define GMX_DD_NNODES_SENDRECV 4
477 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
479 static void index2xyz(ivec nc,int ind,ivec xyz)
481 xyz[XX] = ind % nc[XX];
482 xyz[YY] = (ind / nc[XX]) % nc[YY];
483 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
487 /* This order is required to minimize the coordinate communication in PME
488 * which uses decomposition in the x direction.
490 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
492 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
494 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
495 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
496 xyz[ZZ] = ind % nc[ZZ];
499 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
504 ddindex = dd_index(dd->nc, c);
505 if (dd->comm->bCartesianPP_PME)
507 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
509 else if (dd->comm->bCartesianPP)
512 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
523 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
525 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
528 int ddglatnr(gmx_domdec_t *dd, int i)
538 if (i >= dd->comm->nat[ddnatNR-1])
540 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
542 atnr = dd->gatindex[i] + 1;
548 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
550 return &dd->comm->cgs_gl;
553 static void vec_rvec_init(vec_rvec_t *v)
559 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
563 v->nalloc = over_alloc_dd(n);
564 srenew(v->v, v->nalloc);
568 void dd_store_state(gmx_domdec_t *dd, t_state *state)
572 if (state->ddp_count != dd->ddp_count)
574 gmx_incons("The state does not the domain decomposition state");
577 state->ncg_gl = dd->ncg_home;
578 if (state->ncg_gl > state->cg_gl_nalloc)
580 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
581 srenew(state->cg_gl, state->cg_gl_nalloc);
583 for (i = 0; i < state->ncg_gl; i++)
585 state->cg_gl[i] = dd->index_gl[i];
588 state->ddp_count_cg_gl = dd->ddp_count;
591 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
593 return &dd->comm->zones;
596 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
597 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
599 gmx_domdec_zones_t *zones;
602 zones = &dd->comm->zones;
605 while (icg >= zones->izone[izone].cg1)
614 else if (izone < zones->nizone)
616 *jcg0 = zones->izone[izone].jcg0;
620 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
621 icg, izone, zones->nizone);
624 *jcg1 = zones->izone[izone].jcg1;
626 for (d = 0; d < dd->ndim; d++)
629 shift0[dim] = zones->izone[izone].shift0[dim];
630 shift1[dim] = zones->izone[izone].shift1[dim];
631 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
633 /* A conservative approach, this can be optimized */
640 int dd_natoms_vsite(gmx_domdec_t *dd)
642 return dd->comm->nat[ddnatVSITE];
645 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
647 *at_start = dd->comm->nat[ddnatCON-1];
648 *at_end = dd->comm->nat[ddnatCON];
651 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
653 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
654 int *index, *cgindex;
655 gmx_domdec_comm_t *comm;
656 gmx_domdec_comm_dim_t *cd;
657 gmx_domdec_ind_t *ind;
658 rvec shift = {0, 0, 0}, *buf, *rbuf;
659 gmx_bool bPBC, bScrew;
663 cgindex = dd->cgindex;
668 nat_tot = dd->nat_home;
669 for (d = 0; d < dd->ndim; d++)
671 bPBC = (dd->ci[dd->dim[d]] == 0);
672 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
675 copy_rvec(box[dd->dim[d]], shift);
678 for (p = 0; p < cd->np; p++)
685 for (i = 0; i < ind->nsend[nzone]; i++)
687 at0 = cgindex[index[i]];
688 at1 = cgindex[index[i]+1];
689 for (j = at0; j < at1; j++)
691 copy_rvec(x[j], buf[n]);
698 for (i = 0; i < ind->nsend[nzone]; i++)
700 at0 = cgindex[index[i]];
701 at1 = cgindex[index[i]+1];
702 for (j = at0; j < at1; j++)
704 /* We need to shift the coordinates */
705 rvec_add(x[j], shift, buf[n]);
712 for (i = 0; i < ind->nsend[nzone]; i++)
714 at0 = cgindex[index[i]];
715 at1 = cgindex[index[i]+1];
716 for (j = at0; j < at1; j++)
719 buf[n][XX] = x[j][XX] + shift[XX];
721 * This operation requires a special shift force
722 * treatment, which is performed in calc_vir.
724 buf[n][YY] = box[YY][YY] - x[j][YY];
725 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
737 rbuf = comm->vbuf2.v;
739 /* Send and receive the coordinates */
740 dd_sendrecv_rvec(dd, d, dddirBackward,
741 buf, ind->nsend[nzone+1],
742 rbuf, ind->nrecv[nzone+1]);
746 for (zone = 0; zone < nzone; zone++)
748 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
750 copy_rvec(rbuf[j], x[i]);
755 nat_tot += ind->nrecv[nzone+1];
761 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
763 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
764 int *index, *cgindex;
765 gmx_domdec_comm_t *comm;
766 gmx_domdec_comm_dim_t *cd;
767 gmx_domdec_ind_t *ind;
771 gmx_bool bPBC, bScrew;
775 cgindex = dd->cgindex;
780 nzone = comm->zones.n/2;
781 nat_tot = dd->nat_tot;
782 for (d = dd->ndim-1; d >= 0; d--)
784 bPBC = (dd->ci[dd->dim[d]] == 0);
785 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
786 if (fshift == NULL && !bScrew)
790 /* Determine which shift vector we need */
796 for (p = cd->np-1; p >= 0; p--)
799 nat_tot -= ind->nrecv[nzone+1];
806 sbuf = comm->vbuf2.v;
808 for (zone = 0; zone < nzone; zone++)
810 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
812 copy_rvec(f[i], sbuf[j]);
817 /* Communicate the forces */
818 dd_sendrecv_rvec(dd, d, dddirForward,
819 sbuf, ind->nrecv[nzone+1],
820 buf, ind->nsend[nzone+1]);
822 /* Add the received forces */
826 for (i = 0; i < ind->nsend[nzone]; i++)
828 at0 = cgindex[index[i]];
829 at1 = cgindex[index[i]+1];
830 for (j = at0; j < at1; j++)
832 rvec_inc(f[j], buf[n]);
839 for (i = 0; i < ind->nsend[nzone]; i++)
841 at0 = cgindex[index[i]];
842 at1 = cgindex[index[i]+1];
843 for (j = at0; j < at1; j++)
845 rvec_inc(f[j], buf[n]);
846 /* Add this force to the shift force */
847 rvec_inc(fshift[is], buf[n]);
854 for (i = 0; i < ind->nsend[nzone]; i++)
856 at0 = cgindex[index[i]];
857 at1 = cgindex[index[i]+1];
858 for (j = at0; j < at1; j++)
860 /* Rotate the force */
861 f[j][XX] += buf[n][XX];
862 f[j][YY] -= buf[n][YY];
863 f[j][ZZ] -= buf[n][ZZ];
866 /* Add this force to the shift force */
867 rvec_inc(fshift[is], buf[n]);
878 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
880 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
881 int *index, *cgindex;
882 gmx_domdec_comm_t *comm;
883 gmx_domdec_comm_dim_t *cd;
884 gmx_domdec_ind_t *ind;
889 cgindex = dd->cgindex;
891 buf = &comm->vbuf.v[0][0];
894 nat_tot = dd->nat_home;
895 for (d = 0; d < dd->ndim; d++)
898 for (p = 0; p < cd->np; p++)
903 for (i = 0; i < ind->nsend[nzone]; i++)
905 at0 = cgindex[index[i]];
906 at1 = cgindex[index[i]+1];
907 for (j = at0; j < at1; j++)
920 rbuf = &comm->vbuf2.v[0][0];
922 /* Send and receive the coordinates */
923 dd_sendrecv_real(dd, d, dddirBackward,
924 buf, ind->nsend[nzone+1],
925 rbuf, ind->nrecv[nzone+1]);
929 for (zone = 0; zone < nzone; zone++)
931 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
938 nat_tot += ind->nrecv[nzone+1];
944 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
946 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
947 int *index, *cgindex;
948 gmx_domdec_comm_t *comm;
949 gmx_domdec_comm_dim_t *cd;
950 gmx_domdec_ind_t *ind;
955 cgindex = dd->cgindex;
957 buf = &comm->vbuf.v[0][0];
960 nzone = comm->zones.n/2;
961 nat_tot = dd->nat_tot;
962 for (d = dd->ndim-1; d >= 0; d--)
965 for (p = cd->np-1; p >= 0; p--)
968 nat_tot -= ind->nrecv[nzone+1];
975 sbuf = &comm->vbuf2.v[0][0];
977 for (zone = 0; zone < nzone; zone++)
979 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
986 /* Communicate the forces */
987 dd_sendrecv_real(dd, d, dddirForward,
988 sbuf, ind->nrecv[nzone+1],
989 buf, ind->nsend[nzone+1]);
991 /* Add the received forces */
993 for (i = 0; i < ind->nsend[nzone]; i++)
995 at0 = cgindex[index[i]];
996 at1 = cgindex[index[i]+1];
997 for (j = at0; j < at1; j++)
1008 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1010 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",
1012 zone->min0, zone->max1,
1013 zone->mch0, zone->mch0,
1014 zone->p1_0, zone->p1_1);
1018 #define DDZONECOMM_MAXZONE 5
1019 #define DDZONECOMM_BUFSIZE 3
1021 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1022 int ddimind, int direction,
1023 gmx_ddzone_t *buf_s, int n_s,
1024 gmx_ddzone_t *buf_r, int n_r)
1026 #define ZBS DDZONECOMM_BUFSIZE
1027 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1028 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1031 for (i = 0; i < n_s; i++)
1033 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1034 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1035 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1036 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1037 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1038 vbuf_s[i*ZBS+1][2] = 0;
1039 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1040 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1041 vbuf_s[i*ZBS+2][2] = 0;
1044 dd_sendrecv_rvec(dd, ddimind, direction,
1048 for (i = 0; i < n_r; i++)
1050 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1051 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1052 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1053 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1054 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1055 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1056 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1062 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1063 rvec cell_ns_x0, rvec cell_ns_x1)
1065 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1067 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1068 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1069 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1070 rvec extr_s[2], extr_r[2];
1072 real dist_d, c = 0, det;
1073 gmx_domdec_comm_t *comm;
1074 gmx_bool bPBC, bUse;
1078 for (d = 1; d < dd->ndim; d++)
1081 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1082 zp->min0 = cell_ns_x0[dim];
1083 zp->max1 = cell_ns_x1[dim];
1084 zp->min1 = cell_ns_x1[dim];
1085 zp->mch0 = cell_ns_x0[dim];
1086 zp->mch1 = cell_ns_x1[dim];
1087 zp->p1_0 = cell_ns_x0[dim];
1088 zp->p1_1 = cell_ns_x1[dim];
1091 for (d = dd->ndim-2; d >= 0; d--)
1094 bPBC = (dim < ddbox->npbcdim);
1096 /* Use an rvec to store two reals */
1097 extr_s[d][0] = comm->cell_f0[d+1];
1098 extr_s[d][1] = comm->cell_f1[d+1];
1099 extr_s[d][2] = comm->cell_f1[d+1];
1102 /* Store the extremes in the backward sending buffer,
1103 * so the get updated separately from the forward communication.
1105 for (d1 = d; d1 < dd->ndim-1; d1++)
1107 /* We invert the order to be able to use the same loop for buf_e */
1108 buf_s[pos].min0 = extr_s[d1][1];
1109 buf_s[pos].max1 = extr_s[d1][0];
1110 buf_s[pos].min1 = extr_s[d1][2];
1111 buf_s[pos].mch0 = 0;
1112 buf_s[pos].mch1 = 0;
1113 /* Store the cell corner of the dimension we communicate along */
1114 buf_s[pos].p1_0 = comm->cell_x0[dim];
1115 buf_s[pos].p1_1 = 0;
1119 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1122 if (dd->ndim == 3 && d == 0)
1124 buf_s[pos] = comm->zone_d2[0][1];
1126 buf_s[pos] = comm->zone_d1[0];
1130 /* We only need to communicate the extremes
1131 * in the forward direction
1133 npulse = comm->cd[d].np;
1136 /* Take the minimum to avoid double communication */
1137 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1141 /* Without PBC we should really not communicate over
1142 * the boundaries, but implementing that complicates
1143 * the communication setup and therefore we simply
1144 * do all communication, but ignore some data.
1146 npulse_min = npulse;
1148 for (p = 0; p < npulse_min; p++)
1150 /* Communicate the extremes forward */
1151 bUse = (bPBC || dd->ci[dim] > 0);
1153 dd_sendrecv_rvec(dd, d, dddirForward,
1154 extr_s+d, dd->ndim-d-1,
1155 extr_r+d, dd->ndim-d-1);
1159 for (d1 = d; d1 < dd->ndim-1; d1++)
1161 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1162 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1163 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1169 for (p = 0; p < npulse; p++)
1171 /* Communicate all the zone information backward */
1172 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1174 dd_sendrecv_ddzone(dd, d, dddirBackward,
1181 for (d1 = d+1; d1 < dd->ndim; d1++)
1183 /* Determine the decrease of maximum required
1184 * communication height along d1 due to the distance along d,
1185 * this avoids a lot of useless atom communication.
1187 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1189 if (ddbox->tric_dir[dim])
1191 /* c is the off-diagonal coupling between the cell planes
1192 * along directions d and d1.
1194 c = ddbox->v[dim][dd->dim[d1]][dim];
1200 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1203 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1207 /* A negative value signals out of range */
1213 /* Accumulate the extremes over all pulses */
1214 for (i = 0; i < buf_size; i++)
1218 buf_e[i] = buf_r[i];
1224 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1225 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1226 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1229 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1237 if (bUse && dh[d1] >= 0)
1239 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1240 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1243 /* Copy the received buffer to the send buffer,
1244 * to pass the data through with the next pulse.
1246 buf_s[i] = buf_r[i];
1248 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1249 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1251 /* Store the extremes */
1254 for (d1 = d; d1 < dd->ndim-1; d1++)
1256 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1257 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1258 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1262 if (d == 1 || (d == 0 && dd->ndim == 3))
1264 for (i = d; i < 2; i++)
1266 comm->zone_d2[1-d][i] = buf_e[pos];
1272 comm->zone_d1[1] = buf_e[pos];
1282 for (i = 0; i < 2; i++)
1286 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1288 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1289 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1295 for (i = 0; i < 2; i++)
1297 for (j = 0; j < 2; j++)
1301 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1303 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1304 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1308 for (d = 1; d < dd->ndim; d++)
1310 comm->cell_f_max0[d] = extr_s[d-1][0];
1311 comm->cell_f_min1[d] = extr_s[d-1][1];
1314 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1315 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1320 static void dd_collect_cg(gmx_domdec_t *dd,
1321 t_state *state_local)
1323 gmx_domdec_master_t *ma = NULL;
1324 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1326 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1328 /* The master has the correct distribution */
1332 if (state_local->ddp_count == dd->ddp_count)
1334 /* The local state and DD are in sync, use the DD indices */
1335 ncg_home = dd->ncg_home;
1337 nat_home = dd->nat_home;
1339 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1341 /* The DD is out of sync with the local state, but we have stored
1342 * the cg indices with the local state, so we can use those.
1346 cgs_gl = &dd->comm->cgs_gl;
1348 ncg_home = state_local->ncg_gl;
1349 cg = state_local->cg_gl;
1351 for (i = 0; i < ncg_home; i++)
1353 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1358 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1372 /* Collect the charge group and atom counts on the master */
1373 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1378 for (i = 0; i < dd->nnodes; i++)
1380 ma->ncg[i] = ma->ibuf[2*i];
1381 ma->nat[i] = ma->ibuf[2*i+1];
1382 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1385 /* Make byte counts and indices */
1386 for (i = 0; i < dd->nnodes; i++)
1388 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1389 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1393 fprintf(debug, "Initial charge group distribution: ");
1394 for (i = 0; i < dd->nnodes; i++)
1396 fprintf(debug, " %d", ma->ncg[i]);
1398 fprintf(debug, "\n");
1402 /* Collect the charge group indices on the master */
1404 ncg_home*sizeof(int), cg,
1405 DDMASTER(dd) ? ma->ibuf : NULL,
1406 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1407 DDMASTER(dd) ? ma->cg : NULL);
1409 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1412 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1415 gmx_domdec_master_t *ma;
1416 int n, i, c, a, nalloc = 0;
1425 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1426 dd->rank, dd->mpi_comm_all);
1431 /* Copy the master coordinates to the global array */
1432 cgs_gl = &dd->comm->cgs_gl;
1434 n = DDMASTERRANK(dd);
1436 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1438 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1440 copy_rvec(lv[a++], v[c]);
1444 for (n = 0; n < dd->nnodes; n++)
1448 if (ma->nat[n] > nalloc)
1450 nalloc = over_alloc_dd(ma->nat[n]);
1451 srenew(buf, nalloc);
1454 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1455 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1458 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1460 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1462 copy_rvec(buf[a++], v[c]);
1471 static void get_commbuffer_counts(gmx_domdec_t *dd,
1472 int **counts, int **disps)
1474 gmx_domdec_master_t *ma;
1479 /* Make the rvec count and displacment arrays */
1481 *disps = ma->ibuf + dd->nnodes;
1482 for (n = 0; n < dd->nnodes; n++)
1484 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1485 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1489 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1492 gmx_domdec_master_t *ma;
1493 int *rcounts = NULL, *disps = NULL;
1502 get_commbuffer_counts(dd, &rcounts, &disps);
1507 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1511 cgs_gl = &dd->comm->cgs_gl;
1514 for (n = 0; n < dd->nnodes; n++)
1516 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1518 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1520 copy_rvec(buf[a++], v[c]);
1527 void dd_collect_vec(gmx_domdec_t *dd,
1528 t_state *state_local, rvec *lv, rvec *v)
1530 gmx_domdec_master_t *ma;
1531 int n, i, c, a, nalloc = 0;
1534 dd_collect_cg(dd, state_local);
1536 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1538 dd_collect_vec_sendrecv(dd, lv, v);
1542 dd_collect_vec_gatherv(dd, lv, v);
1547 void dd_collect_state(gmx_domdec_t *dd,
1548 t_state *state_local, t_state *state)
1552 nh = state->nhchainlength;
1556 for (i = 0; i < efptNR; i++)
1558 state->lambda[i] = state_local->lambda[i];
1560 state->fep_state = state_local->fep_state;
1561 state->veta = state_local->veta;
1562 state->vol0 = state_local->vol0;
1563 copy_mat(state_local->box, state->box);
1564 copy_mat(state_local->boxv, state->boxv);
1565 copy_mat(state_local->svir_prev, state->svir_prev);
1566 copy_mat(state_local->fvir_prev, state->fvir_prev);
1567 copy_mat(state_local->pres_prev, state->pres_prev);
1569 for (i = 0; i < state_local->ngtc; i++)
1571 for (j = 0; j < nh; j++)
1573 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1574 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1576 state->therm_integral[i] = state_local->therm_integral[i];
1578 for (i = 0; i < state_local->nnhpres; i++)
1580 for (j = 0; j < nh; j++)
1582 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1583 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1587 for (est = 0; est < estNR; est++)
1589 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1594 dd_collect_vec(dd, state_local, state_local->x, state->x);
1597 dd_collect_vec(dd, state_local, state_local->v, state->v);
1600 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1603 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1605 case estDISRE_INITF:
1606 case estDISRE_RM3TAV:
1607 case estORIRE_INITF:
1611 gmx_incons("Unknown state entry encountered in dd_collect_state");
1617 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1623 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1626 state->nalloc = over_alloc_dd(nalloc);
1628 for (est = 0; est < estNR; est++)
1630 if (EST_DISTR(est) && (state->flags & (1<<est)))
1635 srenew(state->x, state->nalloc);
1638 srenew(state->v, state->nalloc);
1641 srenew(state->sd_X, state->nalloc);
1644 srenew(state->cg_p, state->nalloc);
1646 case estDISRE_INITF:
1647 case estDISRE_RM3TAV:
1648 case estORIRE_INITF:
1650 /* No reallocation required */
1653 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1660 srenew(*f, state->nalloc);
1664 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1667 if (nalloc > fr->cg_nalloc)
1671 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1673 fr->cg_nalloc = over_alloc_dd(nalloc);
1674 srenew(fr->cginfo, fr->cg_nalloc);
1675 if (fr->cutoff_scheme == ecutsGROUP)
1677 srenew(fr->cg_cm, fr->cg_nalloc);
1680 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1682 /* We don't use charge groups, we use x in state to set up
1683 * the atom communication.
1685 dd_realloc_state(state, f, nalloc);
1689 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1692 gmx_domdec_master_t *ma;
1693 int n, i, c, a, nalloc = 0;
1700 for (n = 0; n < dd->nnodes; n++)
1704 if (ma->nat[n] > nalloc)
1706 nalloc = over_alloc_dd(ma->nat[n]);
1707 srenew(buf, nalloc);
1709 /* Use lv as a temporary buffer */
1711 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1713 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1715 copy_rvec(v[c], buf[a++]);
1718 if (a != ma->nat[n])
1720 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1725 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1726 DDRANK(dd, n), n, dd->mpi_comm_all);
1731 n = DDMASTERRANK(dd);
1733 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1735 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1737 copy_rvec(v[c], lv[a++]);
1744 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1745 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1750 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1753 gmx_domdec_master_t *ma;
1754 int *scounts = NULL, *disps = NULL;
1755 int n, i, c, a, nalloc = 0;
1762 get_commbuffer_counts(dd, &scounts, &disps);
1766 for (n = 0; n < dd->nnodes; n++)
1768 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1770 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1772 copy_rvec(v[c], buf[a++]);
1778 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1781 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1783 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1785 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1789 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1793 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1796 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1797 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1798 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1800 if (dfhist->nlambda > 0)
1802 int nlam = dfhist->nlambda;
1803 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1804 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1805 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1806 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1807 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1808 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1810 for (i = 0; i < nlam; i++)
1812 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1813 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1814 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1815 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1816 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1817 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1822 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1823 t_state *state, t_state *state_local,
1828 nh = state->nhchainlength;
1832 for (i = 0; i < efptNR; i++)
1834 state_local->lambda[i] = state->lambda[i];
1836 state_local->fep_state = state->fep_state;
1837 state_local->veta = state->veta;
1838 state_local->vol0 = state->vol0;
1839 copy_mat(state->box, state_local->box);
1840 copy_mat(state->box_rel, state_local->box_rel);
1841 copy_mat(state->boxv, state_local->boxv);
1842 copy_mat(state->svir_prev, state_local->svir_prev);
1843 copy_mat(state->fvir_prev, state_local->fvir_prev);
1844 copy_df_history(&state_local->dfhist, &state->dfhist);
1845 for (i = 0; i < state_local->ngtc; i++)
1847 for (j = 0; j < nh; j++)
1849 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1850 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1852 state_local->therm_integral[i] = state->therm_integral[i];
1854 for (i = 0; i < state_local->nnhpres; i++)
1856 for (j = 0; j < nh; j++)
1858 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1859 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1863 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1864 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1865 dd_bcast(dd, sizeof(real), &state_local->veta);
1866 dd_bcast(dd, sizeof(real), &state_local->vol0);
1867 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1868 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1869 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1870 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1871 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1872 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1873 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1874 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1875 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1876 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1878 /* communicate df_history -- required for restarting from checkpoint */
1879 dd_distribute_dfhist(dd, &state_local->dfhist);
1881 if (dd->nat_home > state_local->nalloc)
1883 dd_realloc_state(state_local, f, dd->nat_home);
1885 for (i = 0; i < estNR; i++)
1887 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1892 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1895 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1898 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1901 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1903 case estDISRE_INITF:
1904 case estDISRE_RM3TAV:
1905 case estORIRE_INITF:
1907 /* Not implemented yet */
1910 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1916 static char dim2char(int dim)
1922 case XX: c = 'X'; break;
1923 case YY: c = 'Y'; break;
1924 case ZZ: c = 'Z'; break;
1925 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1931 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1932 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1934 rvec grid_s[2], *grid_r = NULL, cx, r;
1935 char fname[STRLEN], buf[22];
1937 int a, i, d, z, y, x;
1941 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1942 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1946 snew(grid_r, 2*dd->nnodes);
1949 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1953 for (d = 0; d < DIM; d++)
1955 for (i = 0; i < DIM; i++)
1963 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1965 tric[d][i] = box[i][d]/box[i][i];
1974 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1975 out = gmx_fio_fopen(fname, "w");
1976 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1978 for (i = 0; i < dd->nnodes; i++)
1980 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1981 for (d = 0; d < DIM; d++)
1983 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1985 for (z = 0; z < 2; z++)
1987 for (y = 0; y < 2; y++)
1989 for (x = 0; x < 2; x++)
1991 cx[XX] = grid_r[i*2+x][XX];
1992 cx[YY] = grid_r[i*2+y][YY];
1993 cx[ZZ] = grid_r[i*2+z][ZZ];
1995 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
1996 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2000 for (d = 0; d < DIM; d++)
2002 for (x = 0; x < 4; x++)
2006 case 0: y = 1 + i*8 + 2*x; break;
2007 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2008 case 2: y = 1 + i*8 + x; break;
2010 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2014 gmx_fio_fclose(out);
2019 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2020 gmx_mtop_t *mtop, t_commrec *cr,
2021 int natoms, rvec x[], matrix box)
2023 char fname[STRLEN], buf[22];
2025 int i, ii, resnr, c;
2026 char *atomname, *resname;
2033 natoms = dd->comm->nat[ddnatVSITE];
2036 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2038 out = gmx_fio_fopen(fname, "w");
2040 fprintf(out, "TITLE %s\n", title);
2041 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2042 for (i = 0; i < natoms; i++)
2044 ii = dd->gatindex[i];
2045 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2046 if (i < dd->comm->nat[ddnatZONE])
2049 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2055 else if (i < dd->comm->nat[ddnatVSITE])
2057 b = dd->comm->zones.n;
2061 b = dd->comm->zones.n + 1;
2063 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2064 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2066 fprintf(out, "TER\n");
2068 gmx_fio_fclose(out);
2071 real dd_cutoff_mbody(gmx_domdec_t *dd)
2073 gmx_domdec_comm_t *comm;
2080 if (comm->bInterCGBondeds)
2082 if (comm->cutoff_mbody > 0)
2084 r = comm->cutoff_mbody;
2088 /* cutoff_mbody=0 means we do not have DLB */
2089 r = comm->cellsize_min[dd->dim[0]];
2090 for (di = 1; di < dd->ndim; di++)
2092 r = min(r, comm->cellsize_min[dd->dim[di]]);
2094 if (comm->bBondComm)
2096 r = max(r, comm->cutoff_mbody);
2100 r = min(r, comm->cutoff);
2108 real dd_cutoff_twobody(gmx_domdec_t *dd)
2112 r_mb = dd_cutoff_mbody(dd);
2114 return max(dd->comm->cutoff, r_mb);
2118 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2122 nc = dd->nc[dd->comm->cartpmedim];
2123 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2124 copy_ivec(coord, coord_pme);
2125 coord_pme[dd->comm->cartpmedim] =
2126 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2129 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2131 /* Here we assign a PME node to communicate with this DD node
2132 * by assuming that the major index of both is x.
2133 * We add cr->npmenodes/2 to obtain an even distribution.
2135 return (ddindex*npme + npme/2)/ndd;
2138 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2140 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2143 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2145 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2148 static int *dd_pmenodes(t_commrec *cr)
2153 snew(pmenodes, cr->npmenodes);
2155 for (i = 0; i < cr->dd->nnodes; i++)
2157 p0 = cr_ddindex2pmeindex(cr, i);
2158 p1 = cr_ddindex2pmeindex(cr, i+1);
2159 if (i+1 == cr->dd->nnodes || p1 > p0)
2163 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2165 pmenodes[n] = i + 1 + n;
2173 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2176 ivec coords, coords_pme, nc;
2181 if (dd->comm->bCartesian) {
2182 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2183 dd_coords2pmecoords(dd,coords,coords_pme);
2184 copy_ivec(dd->ntot,nc);
2185 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2186 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2188 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2190 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2196 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2201 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2203 gmx_domdec_comm_t *comm;
2205 int ddindex, nodeid = -1;
2207 comm = cr->dd->comm;
2212 if (comm->bCartesianPP_PME)
2215 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2220 ddindex = dd_index(cr->dd->nc, coords);
2221 if (comm->bCartesianPP)
2223 nodeid = comm->ddindex2simnodeid[ddindex];
2229 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2241 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2244 gmx_domdec_comm_t *comm;
2245 ivec coord, coord_pme;
2252 /* This assumes a uniform x domain decomposition grid cell size */
2253 if (comm->bCartesianPP_PME)
2256 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2257 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2259 /* This is a PP node */
2260 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2261 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2265 else if (comm->bCartesianPP)
2267 if (sim_nodeid < dd->nnodes)
2269 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2274 /* This assumes DD cells with identical x coordinates
2275 * are numbered sequentially.
2277 if (dd->comm->pmenodes == NULL)
2279 if (sim_nodeid < dd->nnodes)
2281 /* The DD index equals the nodeid */
2282 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2288 while (sim_nodeid > dd->comm->pmenodes[i])
2292 if (sim_nodeid < dd->comm->pmenodes[i])
2294 pmenode = dd->comm->pmenodes[i];
2302 void get_pme_nnodes(const gmx_domdec_t *dd,
2303 int *npmenodes_x, int *npmenodes_y)
2307 *npmenodes_x = dd->comm->npmenodes_x;
2308 *npmenodes_y = dd->comm->npmenodes_y;
2317 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2319 gmx_bool bPMEOnlyNode;
2321 if (DOMAINDECOMP(cr))
2323 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2327 bPMEOnlyNode = FALSE;
2330 return bPMEOnlyNode;
2333 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2334 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2338 ivec coord, coord_pme;
2342 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2345 for (x = 0; x < dd->nc[XX]; x++)
2347 for (y = 0; y < dd->nc[YY]; y++)
2349 for (z = 0; z < dd->nc[ZZ]; z++)
2351 if (dd->comm->bCartesianPP_PME)
2356 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2357 if (dd->ci[XX] == coord_pme[XX] &&
2358 dd->ci[YY] == coord_pme[YY] &&
2359 dd->ci[ZZ] == coord_pme[ZZ])
2361 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2366 /* The slab corresponds to the nodeid in the PME group */
2367 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2369 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2376 /* The last PP-only node is the peer node */
2377 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2381 fprintf(debug, "Receive coordinates from PP ranks:");
2382 for (x = 0; x < *nmy_ddnodes; x++)
2384 fprintf(debug, " %d", (*my_ddnodes)[x]);
2386 fprintf(debug, "\n");
2390 static gmx_bool receive_vir_ener(t_commrec *cr)
2392 gmx_domdec_comm_t *comm;
2393 int pmenode, coords[DIM], rank;
2397 if (cr->npmenodes < cr->dd->nnodes)
2399 comm = cr->dd->comm;
2400 if (comm->bCartesianPP_PME)
2402 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2404 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2405 coords[comm->cartpmedim]++;
2406 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2408 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2409 if (dd_simnode2pmenode(cr, rank) == pmenode)
2411 /* This is not the last PP node for pmenode */
2419 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2420 if (cr->sim_nodeid+1 < cr->nnodes &&
2421 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2423 /* This is not the last PP node for pmenode */
2432 static void set_zones_ncg_home(gmx_domdec_t *dd)
2434 gmx_domdec_zones_t *zones;
2437 zones = &dd->comm->zones;
2439 zones->cg_range[0] = 0;
2440 for (i = 1; i < zones->n+1; i++)
2442 zones->cg_range[i] = dd->ncg_home;
2444 /* zone_ncg1[0] should always be equal to ncg_home */
2445 dd->comm->zone_ncg1[0] = dd->ncg_home;
2448 static void rebuild_cgindex(gmx_domdec_t *dd,
2449 const int *gcgs_index, t_state *state)
2451 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2454 dd_cg_gl = dd->index_gl;
2455 cgindex = dd->cgindex;
2458 for (i = 0; i < state->ncg_gl; i++)
2462 dd_cg_gl[i] = cg_gl;
2463 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2467 dd->ncg_home = state->ncg_gl;
2470 set_zones_ncg_home(dd);
2473 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2475 while (cg >= cginfo_mb->cg_end)
2480 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2483 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2484 t_forcerec *fr, char *bLocalCG)
2486 cginfo_mb_t *cginfo_mb;
2492 cginfo_mb = fr->cginfo_mb;
2493 cginfo = fr->cginfo;
2495 for (cg = cg0; cg < cg1; cg++)
2497 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2501 if (bLocalCG != NULL)
2503 for (cg = cg0; cg < cg1; cg++)
2505 bLocalCG[index_gl[cg]] = TRUE;
2510 static void make_dd_indices(gmx_domdec_t *dd,
2511 const int *gcgs_index, int cg_start)
2513 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2514 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2519 bLocalCG = dd->comm->bLocalCG;
2521 if (dd->nat_tot > dd->gatindex_nalloc)
2523 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2524 srenew(dd->gatindex, dd->gatindex_nalloc);
2527 nzone = dd->comm->zones.n;
2528 zone2cg = dd->comm->zones.cg_range;
2529 zone_ncg1 = dd->comm->zone_ncg1;
2530 index_gl = dd->index_gl;
2531 gatindex = dd->gatindex;
2532 bCGs = dd->comm->bCGs;
2534 if (zone2cg[1] != dd->ncg_home)
2536 gmx_incons("dd->ncg_zone is not up to date");
2539 /* Make the local to global and global to local atom index */
2540 a = dd->cgindex[cg_start];
2541 for (zone = 0; zone < nzone; zone++)
2549 cg0 = zone2cg[zone];
2551 cg1 = zone2cg[zone+1];
2552 cg1_p1 = cg0 + zone_ncg1[zone];
2554 for (cg = cg0; cg < cg1; cg++)
2559 /* Signal that this cg is from more than one pulse away */
2562 cg_gl = index_gl[cg];
2565 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2568 ga2la_set(dd->ga2la, a_gl, a, zone1);
2574 gatindex[a] = cg_gl;
2575 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2582 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2585 int ncg, i, ngl, nerr;
2588 if (bLocalCG == NULL)
2592 for (i = 0; i < dd->ncg_tot; i++)
2594 if (!bLocalCG[dd->index_gl[i]])
2597 "DD rank %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);
2602 for (i = 0; i < ncg_sys; i++)
2609 if (ngl != dd->ncg_tot)
2611 fprintf(stderr, "DD rank %d, %s: In bLocalCG %d cgs are marked as local, whereas there are %d\n", dd->rank, where, ngl, dd->ncg_tot);
2618 static void check_index_consistency(gmx_domdec_t *dd,
2619 int natoms_sys, int ncg_sys,
2622 int nerr, ngl, i, a, cell;
2627 if (dd->comm->DD_debug > 1)
2629 snew(have, natoms_sys);
2630 for (a = 0; a < dd->nat_tot; a++)
2632 if (have[dd->gatindex[a]] > 0)
2634 fprintf(stderr, "DD rank %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
2638 have[dd->gatindex[a]] = a + 1;
2644 snew(have, dd->nat_tot);
2647 for (i = 0; i < natoms_sys; i++)
2649 if (ga2la_get(dd->ga2la, i, &a, &cell))
2651 if (a >= dd->nat_tot)
2653 fprintf(stderr, "DD rank %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);
2659 if (dd->gatindex[a] != i)
2661 fprintf(stderr, "DD rank %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);
2668 if (ngl != dd->nat_tot)
2671 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2672 dd->rank, where, ngl, dd->nat_tot);
2674 for (a = 0; a < dd->nat_tot; a++)
2679 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2680 dd->rank, where, a+1, dd->gatindex[a]+1);
2685 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2689 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2690 dd->rank, where, nerr);
2694 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2701 /* Clear the whole list without searching */
2702 ga2la_clear(dd->ga2la);
2706 for (i = a_start; i < dd->nat_tot; i++)
2708 ga2la_del(dd->ga2la, dd->gatindex[i]);
2712 bLocalCG = dd->comm->bLocalCG;
2715 for (i = cg_start; i < dd->ncg_tot; i++)
2717 bLocalCG[dd->index_gl[i]] = FALSE;
2721 dd_clear_local_vsite_indices(dd);
2723 if (dd->constraints)
2725 dd_clear_local_constraint_indices(dd);
2729 /* This function should be used for moving the domain boudaries during DLB,
2730 * for obtaining the minimum cell size. It checks the initially set limit
2731 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2732 * and, possibly, a longer cut-off limit set for PME load balancing.
2734 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2738 cellsize_min = comm->cellsize_min[dim];
2740 if (!comm->bVacDLBNoLimit)
2742 /* The cut-off might have changed, e.g. by PME load balacning,
2743 * from the value used to set comm->cellsize_min, so check it.
2745 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2747 if (comm->bPMELoadBalDLBLimits)
2749 /* Check for the cut-off limit set by the PME load balancing */
2750 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2754 return cellsize_min;
2757 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2760 real grid_jump_limit;
2762 /* The distance between the boundaries of cells at distance
2763 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2764 * and by the fact that cells should not be shifted by more than
2765 * half their size, such that cg's only shift by one cell
2766 * at redecomposition.
2768 grid_jump_limit = comm->cellsize_limit;
2769 if (!comm->bVacDLBNoLimit)
2771 if (comm->bPMELoadBalDLBLimits)
2773 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2775 grid_jump_limit = max(grid_jump_limit,
2776 cutoff/comm->cd[dim_ind].np);
2779 return grid_jump_limit;
2782 static gmx_bool check_grid_jump(gmx_int64_t step,
2788 gmx_domdec_comm_t *comm;
2797 for (d = 1; d < dd->ndim; d++)
2800 limit = grid_jump_limit(comm, cutoff, d);
2801 bfac = ddbox->box_size[dim];
2802 if (ddbox->tric_dir[dim])
2804 bfac *= ddbox->skew_fac[dim];
2806 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2807 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2815 /* This error should never be triggered under normal
2816 * circumstances, but you never know ...
2818 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 fewer ranks might avoid this issue.",
2819 gmx_step_str(step, buf),
2820 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2828 static int dd_load_count(gmx_domdec_comm_t *comm)
2830 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2833 static float dd_force_load(gmx_domdec_comm_t *comm)
2840 if (comm->eFlop > 1)
2842 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2847 load = comm->cycl[ddCyclF];
2848 if (comm->cycl_n[ddCyclF] > 1)
2850 /* Subtract the maximum of the last n cycle counts
2851 * to get rid of possible high counts due to other sources,
2852 * for instance system activity, that would otherwise
2853 * affect the dynamic load balancing.
2855 load -= comm->cycl_max[ddCyclF];
2859 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2861 float gpu_wait, gpu_wait_sum;
2863 gpu_wait = comm->cycl[ddCyclWaitGPU];
2864 if (comm->cycl_n[ddCyclF] > 1)
2866 /* We should remove the WaitGPU time of the same MD step
2867 * as the one with the maximum F time, since the F time
2868 * and the wait time are not independent.
2869 * Furthermore, the step for the max F time should be chosen
2870 * the same on all ranks that share the same GPU.
2871 * But to keep the code simple, we remove the average instead.
2872 * The main reason for artificially long times at some steps
2873 * is spurious CPU activity or MPI time, so we don't expect
2874 * that changes in the GPU wait time matter a lot here.
2876 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2878 /* Sum the wait times over the ranks that share the same GPU */
2879 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2880 comm->mpi_comm_gpu_shared);
2881 /* Replace the wait time by the average over the ranks */
2882 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2890 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2892 gmx_domdec_comm_t *comm;
2897 snew(*dim_f, dd->nc[dim]+1);
2899 for (i = 1; i < dd->nc[dim]; i++)
2901 if (comm->slb_frac[dim])
2903 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2907 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2910 (*dim_f)[dd->nc[dim]] = 1;
2913 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2915 int pmeindex, slab, nso, i;
2918 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2924 ddpme->dim = dimind;
2926 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2928 ddpme->nslab = (ddpme->dim == 0 ?
2929 dd->comm->npmenodes_x :
2930 dd->comm->npmenodes_y);
2932 if (ddpme->nslab <= 1)
2937 nso = dd->comm->npmenodes/ddpme->nslab;
2938 /* Determine for each PME slab the PP location range for dimension dim */
2939 snew(ddpme->pp_min, ddpme->nslab);
2940 snew(ddpme->pp_max, ddpme->nslab);
2941 for (slab = 0; slab < ddpme->nslab; slab++)
2943 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2944 ddpme->pp_max[slab] = 0;
2946 for (i = 0; i < dd->nnodes; i++)
2948 ddindex2xyz(dd->nc, i, xyz);
2949 /* For y only use our y/z slab.
2950 * This assumes that the PME x grid size matches the DD grid size.
2952 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2954 pmeindex = ddindex2pmeindex(dd, i);
2957 slab = pmeindex/nso;
2961 slab = pmeindex % ddpme->nslab;
2963 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2964 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2968 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2971 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2973 if (dd->comm->ddpme[0].dim == XX)
2975 return dd->comm->ddpme[0].maxshift;
2983 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2985 if (dd->comm->ddpme[0].dim == YY)
2987 return dd->comm->ddpme[0].maxshift;
2989 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2991 return dd->comm->ddpme[1].maxshift;
2999 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3000 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3002 gmx_domdec_comm_t *comm;
3005 real range, pme_boundary;
3009 nc = dd->nc[ddpme->dim];
3012 if (!ddpme->dim_match)
3014 /* PP decomposition is not along dim: the worst situation */
3017 else if (ns <= 3 || (bUniform && ns == nc))
3019 /* The optimal situation */
3024 /* We need to check for all pme nodes which nodes they
3025 * could possibly need to communicate with.
3027 xmin = ddpme->pp_min;
3028 xmax = ddpme->pp_max;
3029 /* Allow for atoms to be maximally 2/3 times the cut-off
3030 * out of their DD cell. This is a reasonable balance between
3031 * between performance and support for most charge-group/cut-off
3034 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3035 /* Avoid extra communication when we are exactly at a boundary */
3039 for (s = 0; s < ns; s++)
3041 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3042 pme_boundary = (real)s/ns;
3045 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3047 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3051 pme_boundary = (real)(s+1)/ns;
3054 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3056 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3063 ddpme->maxshift = sh;
3067 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3068 ddpme->dim, ddpme->maxshift);
3072 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3076 for (d = 0; d < dd->ndim; d++)
3079 if (dim < ddbox->nboundeddim &&
3080 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3081 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3083 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",
3084 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3085 dd->nc[dim], dd->comm->cellsize_limit);
3091 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3094 /* Set the domain boundaries. Use for static (or no) load balancing,
3095 * and also for the starting state for dynamic load balancing.
3096 * setmode determine if and where the boundaries are stored, use enum above.
3097 * Returns the number communication pulses in npulse.
3099 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3100 int setmode, ivec npulse)
3102 gmx_domdec_comm_t *comm;
3105 real *cell_x, cell_dx, cellsize;
3109 for (d = 0; d < DIM; d++)
3111 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3113 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3116 cell_dx = ddbox->box_size[d]/dd->nc[d];
3119 case setcellsizeslbMASTER:
3120 for (j = 0; j < dd->nc[d]+1; j++)
3122 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3125 case setcellsizeslbLOCAL:
3126 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3127 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3132 cellsize = cell_dx*ddbox->skew_fac[d];
3133 while (cellsize*npulse[d] < comm->cutoff)
3137 cellsize_min[d] = cellsize;
3141 /* Statically load balanced grid */
3142 /* Also when we are not doing a master distribution we determine
3143 * all cell borders in a loop to obtain identical values
3144 * to the master distribution case and to determine npulse.
3146 if (setmode == setcellsizeslbMASTER)
3148 cell_x = dd->ma->cell_x[d];
3152 snew(cell_x, dd->nc[d]+1);
3154 cell_x[0] = ddbox->box0[d];
3155 for (j = 0; j < dd->nc[d]; j++)
3157 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3158 cell_x[j+1] = cell_x[j] + cell_dx;
3159 cellsize = cell_dx*ddbox->skew_fac[d];
3160 while (cellsize*npulse[d] < comm->cutoff &&
3161 npulse[d] < dd->nc[d]-1)
3165 cellsize_min[d] = min(cellsize_min[d], cellsize);
3167 if (setmode == setcellsizeslbLOCAL)
3169 comm->cell_x0[d] = cell_x[dd->ci[d]];
3170 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3172 if (setmode != setcellsizeslbMASTER)
3177 /* The following limitation is to avoid that a cell would receive
3178 * some of its own home charge groups back over the periodic boundary.
3179 * Double charge groups cause trouble with the global indices.
3181 if (d < ddbox->npbcdim &&
3182 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3184 char error_string[STRLEN];
3186 sprintf(error_string,
3187 "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",
3188 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3190 dd->nc[d], dd->nc[d],
3191 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3193 if (setmode == setcellsizeslbLOCAL)
3195 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3199 gmx_fatal(FARGS, error_string);
3204 if (!comm->bDynLoadBal)
3206 copy_rvec(cellsize_min, comm->cellsize_min);
3209 for (d = 0; d < comm->npmedecompdim; d++)
3211 set_pme_maxshift(dd, &comm->ddpme[d],
3212 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3213 comm->ddpme[d].slb_dim_f);
3218 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3219 int d, int dim, gmx_domdec_root_t *root,
3221 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3223 gmx_domdec_comm_t *comm;
3224 int ncd, i, j, nmin, nmin_old;
3225 gmx_bool bLimLo, bLimHi;
3227 real fac, halfway, cellsize_limit_f_i, region_size;
3228 gmx_bool bPBC, bLastHi = FALSE;
3229 int nrange[] = {range[0], range[1]};
3231 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3237 bPBC = (dim < ddbox->npbcdim);
3239 cell_size = root->buf_ncd;
3243 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3246 /* First we need to check if the scaling does not make cells
3247 * smaller than the smallest allowed size.
3248 * We need to do this iteratively, since if a cell is too small,
3249 * it needs to be enlarged, which makes all the other cells smaller,
3250 * which could in turn make another cell smaller than allowed.
3252 for (i = range[0]; i < range[1]; i++)
3254 root->bCellMin[i] = FALSE;
3260 /* We need the total for normalization */
3262 for (i = range[0]; i < range[1]; i++)
3264 if (root->bCellMin[i] == FALSE)
3266 fac += cell_size[i];
3269 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3270 /* Determine the cell boundaries */
3271 for (i = range[0]; i < range[1]; i++)
3273 if (root->bCellMin[i] == FALSE)
3275 cell_size[i] *= fac;
3276 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3278 cellsize_limit_f_i = 0;
3282 cellsize_limit_f_i = cellsize_limit_f;
3284 if (cell_size[i] < cellsize_limit_f_i)
3286 root->bCellMin[i] = TRUE;
3287 cell_size[i] = cellsize_limit_f_i;
3291 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3294 while (nmin > nmin_old);
3297 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3298 /* For this check we should not use DD_CELL_MARGIN,
3299 * but a slightly smaller factor,
3300 * since rounding could get use below the limit.
3302 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3305 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",
3306 gmx_step_str(step, buf),
3307 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3308 ncd, comm->cellsize_min[dim]);
3311 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3315 /* Check if the boundary did not displace more than halfway
3316 * each of the cells it bounds, as this could cause problems,
3317 * especially when the differences between cell sizes are large.
3318 * If changes are applied, they will not make cells smaller
3319 * than the cut-off, as we check all the boundaries which
3320 * might be affected by a change and if the old state was ok,
3321 * the cells will at most be shrunk back to their old size.
3323 for (i = range[0]+1; i < range[1]; i++)
3325 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3326 if (root->cell_f[i] < halfway)
3328 root->cell_f[i] = halfway;
3329 /* Check if the change also causes shifts of the next boundaries */
3330 for (j = i+1; j < range[1]; j++)
3332 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3334 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3338 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3339 if (root->cell_f[i] > halfway)
3341 root->cell_f[i] = halfway;
3342 /* Check if the change also causes shifts of the next boundaries */
3343 for (j = i-1; j >= range[0]+1; j--)
3345 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3347 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3354 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3355 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3356 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3357 * for a and b nrange is used */
3360 /* Take care of the staggering of the cell boundaries */
3363 for (i = range[0]; i < range[1]; i++)
3365 root->cell_f_max0[i] = root->cell_f[i];
3366 root->cell_f_min1[i] = root->cell_f[i+1];
3371 for (i = range[0]+1; i < range[1]; i++)
3373 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3374 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3375 if (bLimLo && bLimHi)
3377 /* Both limits violated, try the best we can */
3378 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3379 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3380 nrange[0] = range[0];
3382 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3385 nrange[1] = range[1];
3386 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3392 /* root->cell_f[i] = root->bound_min[i]; */
3393 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3396 else if (bLimHi && !bLastHi)
3399 if (nrange[1] < range[1]) /* found a LimLo before */
3401 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3402 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3403 nrange[0] = nrange[1];
3405 root->cell_f[i] = root->bound_max[i];
3407 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3409 nrange[1] = range[1];
3412 if (nrange[1] < range[1]) /* found last a LimLo */
3414 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3415 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3416 nrange[0] = nrange[1];
3417 nrange[1] = range[1];
3418 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3420 else if (nrange[0] > range[0]) /* found at least one LimHi */
3422 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3429 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3430 int d, int dim, gmx_domdec_root_t *root,
3431 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3432 gmx_bool bUniform, gmx_int64_t step)
3434 gmx_domdec_comm_t *comm;
3435 int ncd, d1, i, j, pos;
3437 real load_aver, load_i, imbalance, change, change_max, sc;
3438 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3442 int range[] = { 0, 0 };
3446 /* Convert the maximum change from the input percentage to a fraction */
3447 change_limit = comm->dlb_scale_lim*0.01;
3451 bPBC = (dim < ddbox->npbcdim);
3453 cell_size = root->buf_ncd;
3455 /* Store the original boundaries */
3456 for (i = 0; i < ncd+1; i++)
3458 root->old_cell_f[i] = root->cell_f[i];
3462 for (i = 0; i < ncd; i++)
3464 cell_size[i] = 1.0/ncd;
3467 else if (dd_load_count(comm) > 0)
3469 load_aver = comm->load[d].sum_m/ncd;
3471 for (i = 0; i < ncd; i++)
3473 /* Determine the relative imbalance of cell i */
3474 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3475 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3476 /* Determine the change of the cell size using underrelaxation */
3477 change = -relax*imbalance;
3478 change_max = max(change_max, max(change, -change));
3480 /* Limit the amount of scaling.
3481 * We need to use the same rescaling for all cells in one row,
3482 * otherwise the load balancing might not converge.
3485 if (change_max > change_limit)
3487 sc *= change_limit/change_max;
3489 for (i = 0; i < ncd; i++)
3491 /* Determine the relative imbalance of cell i */
3492 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3493 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3494 /* Determine the change of the cell size using underrelaxation */
3495 change = -sc*imbalance;
3496 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3500 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3501 cellsize_limit_f *= DD_CELL_MARGIN;
3502 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3503 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3504 if (ddbox->tric_dir[dim])
3506 cellsize_limit_f /= ddbox->skew_fac[dim];
3507 dist_min_f /= ddbox->skew_fac[dim];
3509 if (bDynamicBox && d > 0)
3511 dist_min_f *= DD_PRES_SCALE_MARGIN;
3513 if (d > 0 && !bUniform)
3515 /* Make sure that the grid is not shifted too much */
3516 for (i = 1; i < ncd; i++)
3518 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3520 gmx_incons("Inconsistent DD boundary staggering limits!");
3522 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3523 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3526 root->bound_min[i] += 0.5*space;
3528 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3529 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3532 root->bound_max[i] += 0.5*space;
3537 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3539 root->cell_f_max0[i-1] + dist_min_f,
3540 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3541 root->cell_f_min1[i] - dist_min_f);
3546 root->cell_f[0] = 0;
3547 root->cell_f[ncd] = 1;
3548 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3551 /* After the checks above, the cells should obey the cut-off
3552 * restrictions, but it does not hurt to check.
3554 for (i = 0; i < ncd; i++)
3558 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3559 dim, i, root->cell_f[i], root->cell_f[i+1]);
3562 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3563 root->cell_f[i+1] - root->cell_f[i] <
3564 cellsize_limit_f/DD_CELL_MARGIN)
3568 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3569 gmx_step_str(step, buf), dim2char(dim), i,
3570 (root->cell_f[i+1] - root->cell_f[i])
3571 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3576 /* Store the cell boundaries of the lower dimensions at the end */
3577 for (d1 = 0; d1 < d; d1++)
3579 root->cell_f[pos++] = comm->cell_f0[d1];
3580 root->cell_f[pos++] = comm->cell_f1[d1];
3583 if (d < comm->npmedecompdim)
3585 /* The master determines the maximum shift for
3586 * the coordinate communication between separate PME nodes.
3588 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3590 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3593 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3597 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3598 gmx_ddbox_t *ddbox, int dimind)
3600 gmx_domdec_comm_t *comm;
3605 /* Set the cell dimensions */
3606 dim = dd->dim[dimind];
3607 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3608 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3609 if (dim >= ddbox->nboundeddim)
3611 comm->cell_x0[dim] += ddbox->box0[dim];
3612 comm->cell_x1[dim] += ddbox->box0[dim];
3616 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3617 int d, int dim, real *cell_f_row,
3620 gmx_domdec_comm_t *comm;
3626 /* Each node would only need to know two fractions,
3627 * but it is probably cheaper to broadcast the whole array.
3629 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3630 0, comm->mpi_comm_load[d]);
3632 /* Copy the fractions for this dimension from the buffer */
3633 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3634 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3635 /* The whole array was communicated, so set the buffer position */
3636 pos = dd->nc[dim] + 1;
3637 for (d1 = 0; d1 <= d; d1++)
3641 /* Copy the cell fractions of the lower dimensions */
3642 comm->cell_f0[d1] = cell_f_row[pos++];
3643 comm->cell_f1[d1] = cell_f_row[pos++];
3645 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3647 /* Convert the communicated shift from float to int */
3648 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3651 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3655 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3656 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3657 gmx_bool bUniform, gmx_int64_t step)
3659 gmx_domdec_comm_t *comm;
3661 gmx_bool bRowMember, bRowRoot;
3666 for (d = 0; d < dd->ndim; d++)
3671 for (d1 = d; d1 < dd->ndim; d1++)
3673 if (dd->ci[dd->dim[d1]] > 0)
3686 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3687 ddbox, bDynamicBox, bUniform, step);
3688 cell_f_row = comm->root[d]->cell_f;
3692 cell_f_row = comm->cell_f_row;
3694 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3699 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3703 /* This function assumes the box is static and should therefore
3704 * not be called when the box has changed since the last
3705 * call to dd_partition_system.
3707 for (d = 0; d < dd->ndim; d++)
3709 relative_to_absolute_cell_bounds(dd, ddbox, d);
3715 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3716 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3717 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3718 gmx_wallcycle_t wcycle)
3720 gmx_domdec_comm_t *comm;
3727 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3728 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3729 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3731 else if (bDynamicBox)
3733 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3736 /* Set the dimensions for which no DD is used */
3737 for (dim = 0; dim < DIM; dim++)
3739 if (dd->nc[dim] == 1)
3741 comm->cell_x0[dim] = 0;
3742 comm->cell_x1[dim] = ddbox->box_size[dim];
3743 if (dim >= ddbox->nboundeddim)
3745 comm->cell_x0[dim] += ddbox->box0[dim];
3746 comm->cell_x1[dim] += ddbox->box0[dim];
3752 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3755 gmx_domdec_comm_dim_t *cd;
3757 for (d = 0; d < dd->ndim; d++)
3759 cd = &dd->comm->cd[d];
3760 np = npulse[dd->dim[d]];
3761 if (np > cd->np_nalloc)
3765 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3766 dim2char(dd->dim[d]), np);
3768 if (DDMASTER(dd) && cd->np_nalloc > 0)
3770 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3772 srenew(cd->ind, np);
3773 for (i = cd->np_nalloc; i < np; i++)
3775 cd->ind[i].index = NULL;
3776 cd->ind[i].nalloc = 0;
3785 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3786 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3787 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3788 gmx_wallcycle_t wcycle)
3790 gmx_domdec_comm_t *comm;
3796 /* Copy the old cell boundaries for the cg displacement check */
3797 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3798 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3800 if (comm->bDynLoadBal)
3804 check_box_size(dd, ddbox);
3806 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3810 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3811 realloc_comm_ind(dd, npulse);
3816 for (d = 0; d < DIM; d++)
3818 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3819 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3824 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3826 rvec cell_ns_x0, rvec cell_ns_x1,
3829 gmx_domdec_comm_t *comm;
3834 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3836 dim = dd->dim[dim_ind];
3838 /* Without PBC we don't have restrictions on the outer cells */
3839 if (!(dim >= ddbox->npbcdim &&
3840 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3841 comm->bDynLoadBal &&
3842 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3843 comm->cellsize_min[dim])
3846 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",
3847 gmx_step_str(step, buf), dim2char(dim),
3848 comm->cell_x1[dim] - comm->cell_x0[dim],
3849 ddbox->skew_fac[dim],
3850 dd->comm->cellsize_min[dim],
3851 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3855 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3857 /* Communicate the boundaries and update cell_ns_x0/1 */
3858 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3859 if (dd->bGridJump && dd->ndim > 1)
3861 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3866 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3870 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3878 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3879 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3888 static void check_screw_box(matrix box)
3890 /* Mathematical limitation */
3891 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3893 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3896 /* Limitation due to the asymmetry of the eighth shell method */
3897 if (box[ZZ][YY] != 0)
3899 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3903 static void distribute_cg(FILE *fplog, gmx_int64_t step,
3904 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3907 gmx_domdec_master_t *ma;
3908 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3909 int i, icg, j, k, k0, k1, d, npbcdim;
3911 rvec box_size, cg_cm;
3913 real nrcg, inv_ncg, pos_d;
3915 gmx_bool bUnbounded, bScrew;
3919 if (tmp_ind == NULL)
3921 snew(tmp_nalloc, dd->nnodes);
3922 snew(tmp_ind, dd->nnodes);
3923 for (i = 0; i < dd->nnodes; i++)
3925 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3926 snew(tmp_ind[i], tmp_nalloc[i]);
3930 /* Clear the count */
3931 for (i = 0; i < dd->nnodes; i++)
3937 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3939 cgindex = cgs->index;
3941 /* Compute the center of geometry for all charge groups */
3942 for (icg = 0; icg < cgs->nr; icg++)
3945 k1 = cgindex[icg+1];
3949 copy_rvec(pos[k0], cg_cm);
3956 for (k = k0; (k < k1); k++)
3958 rvec_inc(cg_cm, pos[k]);
3960 for (d = 0; (d < DIM); d++)
3962 cg_cm[d] *= inv_ncg;
3965 /* Put the charge group in the box and determine the cell index */
3966 for (d = DIM-1; d >= 0; d--)
3969 if (d < dd->npbcdim)
3971 bScrew = (dd->bScrewPBC && d == XX);
3972 if (tric_dir[d] && dd->nc[d] > 1)
3974 /* Use triclinic coordintates for this dimension */
3975 for (j = d+1; j < DIM; j++)
3977 pos_d += cg_cm[j]*tcm[j][d];
3980 while (pos_d >= box[d][d])
3983 rvec_dec(cg_cm, box[d]);
3986 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3987 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3989 for (k = k0; (k < k1); k++)
3991 rvec_dec(pos[k], box[d]);
3994 pos[k][YY] = box[YY][YY] - pos[k][YY];
3995 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4002 rvec_inc(cg_cm, box[d]);
4005 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4006 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4008 for (k = k0; (k < k1); k++)
4010 rvec_inc(pos[k], box[d]);
4013 pos[k][YY] = box[YY][YY] - pos[k][YY];
4014 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4019 /* This could be done more efficiently */
4021 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4026 i = dd_index(dd->nc, ind);
4027 if (ma->ncg[i] == tmp_nalloc[i])
4029 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4030 srenew(tmp_ind[i], tmp_nalloc[i]);
4032 tmp_ind[i][ma->ncg[i]] = icg;
4034 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4038 for (i = 0; i < dd->nnodes; i++)
4041 for (k = 0; k < ma->ncg[i]; k++)
4043 ma->cg[k1++] = tmp_ind[i][k];
4046 ma->index[dd->nnodes] = k1;
4048 for (i = 0; i < dd->nnodes; i++)
4058 fprintf(fplog, "Charge group distribution at step %s:",
4059 gmx_step_str(step, buf));
4060 for (i = 0; i < dd->nnodes; i++)
4062 fprintf(fplog, " %d", ma->ncg[i]);
4064 fprintf(fplog, "\n");
4068 static void get_cg_distribution(FILE *fplog, gmx_int64_t step, gmx_domdec_t *dd,
4069 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4072 gmx_domdec_master_t *ma = NULL;
4075 int *ibuf, buf2[2] = { 0, 0 };
4076 gmx_bool bMaster = DDMASTER(dd);
4084 check_screw_box(box);
4087 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4089 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4090 for (i = 0; i < dd->nnodes; i++)
4092 ma->ibuf[2*i] = ma->ncg[i];
4093 ma->ibuf[2*i+1] = ma->nat[i];
4101 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4103 dd->ncg_home = buf2[0];
4104 dd->nat_home = buf2[1];
4105 dd->ncg_tot = dd->ncg_home;
4106 dd->nat_tot = dd->nat_home;
4107 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4109 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4110 srenew(dd->index_gl, dd->cg_nalloc);
4111 srenew(dd->cgindex, dd->cg_nalloc+1);
4115 for (i = 0; i < dd->nnodes; i++)
4117 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4118 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4123 DDMASTER(dd) ? ma->ibuf : NULL,
4124 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4125 DDMASTER(dd) ? ma->cg : NULL,
4126 dd->ncg_home*sizeof(int), dd->index_gl);
4128 /* Determine the home charge group sizes */
4130 for (i = 0; i < dd->ncg_home; i++)
4132 cg_gl = dd->index_gl[i];
4134 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4139 fprintf(debug, "Home charge groups:\n");
4140 for (i = 0; i < dd->ncg_home; i++)
4142 fprintf(debug, " %d", dd->index_gl[i]);
4145 fprintf(debug, "\n");
4148 fprintf(debug, "\n");
4152 static int compact_and_copy_vec_at(int ncg, int *move,
4155 rvec *src, gmx_domdec_comm_t *comm,
4158 int m, icg, i, i0, i1, nrcg;
4164 for (m = 0; m < DIM*2; m++)
4170 for (icg = 0; icg < ncg; icg++)
4172 i1 = cgindex[icg+1];
4178 /* Compact the home array in place */
4179 for (i = i0; i < i1; i++)
4181 copy_rvec(src[i], src[home_pos++]);
4187 /* Copy to the communication buffer */
4189 pos_vec[m] += 1 + vec*nrcg;
4190 for (i = i0; i < i1; i++)
4192 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4194 pos_vec[m] += (nvec - vec - 1)*nrcg;
4198 home_pos += i1 - i0;
4206 static int compact_and_copy_vec_cg(int ncg, int *move,
4208 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4211 int m, icg, i0, i1, nrcg;
4217 for (m = 0; m < DIM*2; m++)
4223 for (icg = 0; icg < ncg; icg++)
4225 i1 = cgindex[icg+1];
4231 /* Compact the home array in place */
4232 copy_rvec(src[icg], src[home_pos++]);
4238 /* Copy to the communication buffer */
4239 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4240 pos_vec[m] += 1 + nrcg*nvec;
4252 static int compact_ind(int ncg, int *move,
4253 int *index_gl, int *cgindex,
4255 gmx_ga2la_t ga2la, char *bLocalCG,
4258 int cg, nat, a0, a1, a, a_gl;
4263 for (cg = 0; cg < ncg; cg++)
4269 /* Compact the home arrays in place.
4270 * Anything that can be done here avoids access to global arrays.
4272 cgindex[home_pos] = nat;
4273 for (a = a0; a < a1; a++)
4276 gatindex[nat] = a_gl;
4277 /* The cell number stays 0, so we don't need to set it */
4278 ga2la_change_la(ga2la, a_gl, nat);
4281 index_gl[home_pos] = index_gl[cg];
4282 cginfo[home_pos] = cginfo[cg];
4283 /* The charge group remains local, so bLocalCG does not change */
4288 /* Clear the global indices */
4289 for (a = a0; a < a1; a++)
4291 ga2la_del(ga2la, gatindex[a]);
4295 bLocalCG[index_gl[cg]] = FALSE;
4299 cgindex[home_pos] = nat;
4304 static void clear_and_mark_ind(int ncg, int *move,
4305 int *index_gl, int *cgindex, int *gatindex,
4306 gmx_ga2la_t ga2la, char *bLocalCG,
4311 for (cg = 0; cg < ncg; cg++)
4317 /* Clear the global indices */
4318 for (a = a0; a < a1; a++)
4320 ga2la_del(ga2la, gatindex[a]);
4324 bLocalCG[index_gl[cg]] = FALSE;
4326 /* Signal that this cg has moved using the ns cell index.
4327 * Here we set it to -1. fill_grid will change it
4328 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4330 cell_index[cg] = -1;
4335 static void print_cg_move(FILE *fplog,
4337 gmx_int64_t step, int cg, int dim, int dir,
4338 gmx_bool bHaveCgcmOld, real limitd,
4339 rvec cm_old, rvec cm_new, real pos_d)
4341 gmx_domdec_comm_t *comm;
4346 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4349 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4350 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4351 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4355 /* We don't have a limiting distance available: don't print it */
4356 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4357 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4358 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4360 fprintf(fplog, "distance out of cell %f\n",
4361 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4364 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4365 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4367 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4368 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4369 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4371 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4372 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4374 comm->cell_x0[dim], comm->cell_x1[dim]);
4377 static void cg_move_error(FILE *fplog,
4379 gmx_int64_t step, int cg, int dim, int dir,
4380 gmx_bool bHaveCgcmOld, real limitd,
4381 rvec cm_old, rvec cm_new, real pos_d)
4385 print_cg_move(fplog, dd, step, cg, dim, dir,
4386 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4388 print_cg_move(stderr, dd, step, cg, dim, dir,
4389 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4391 "%s moved too far between two domain decomposition steps\n"
4392 "This usually means that your system is not well equilibrated",
4393 dd->comm->bCGs ? "A charge group" : "An atom");
4396 static void rotate_state_atom(t_state *state, int a)
4400 for (est = 0; est < estNR; est++)
4402 if (EST_DISTR(est) && (state->flags & (1<<est)))
4407 /* Rotate the complete state; for a rectangular box only */
4408 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4409 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4412 state->v[a][YY] = -state->v[a][YY];
4413 state->v[a][ZZ] = -state->v[a][ZZ];
4416 state->sd_X[a][YY] = -state->sd_X[a][YY];
4417 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4420 state->cg_p[a][YY] = -state->cg_p[a][YY];
4421 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4423 case estDISRE_INITF:
4424 case estDISRE_RM3TAV:
4425 case estORIRE_INITF:
4427 /* These are distances, so not affected by rotation */
4430 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4436 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4438 if (natoms > comm->moved_nalloc)
4440 /* Contents should be preserved here */
4441 comm->moved_nalloc = over_alloc_dd(natoms);
4442 srenew(comm->moved, comm->moved_nalloc);
4448 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4451 ivec tric_dir, matrix tcm,
4452 rvec cell_x0, rvec cell_x1,
4453 rvec limitd, rvec limit0, rvec limit1,
4455 int cg_start, int cg_end,
4460 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4461 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4465 real inv_ncg, pos_d;
4468 npbcdim = dd->npbcdim;
4470 for (cg = cg_start; cg < cg_end; cg++)
4477 copy_rvec(state->x[k0], cm_new);
4484 for (k = k0; (k < k1); k++)
4486 rvec_inc(cm_new, state->x[k]);
4488 for (d = 0; (d < DIM); d++)
4490 cm_new[d] = inv_ncg*cm_new[d];
4495 /* Do pbc and check DD cell boundary crossings */
4496 for (d = DIM-1; d >= 0; d--)
4500 bScrew = (dd->bScrewPBC && d == XX);
4501 /* Determine the location of this cg in lattice coordinates */
4505 for (d2 = d+1; d2 < DIM; d2++)
4507 pos_d += cm_new[d2]*tcm[d2][d];
4510 /* Put the charge group in the triclinic unit-cell */
4511 if (pos_d >= cell_x1[d])
4513 if (pos_d >= limit1[d])
4515 cg_move_error(fplog, dd, step, cg, d, 1,
4516 cg_cm != state->x, limitd[d],
4517 cg_cm[cg], cm_new, pos_d);
4520 if (dd->ci[d] == dd->nc[d] - 1)
4522 rvec_dec(cm_new, state->box[d]);
4525 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4526 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4528 for (k = k0; (k < k1); k++)
4530 rvec_dec(state->x[k], state->box[d]);
4533 rotate_state_atom(state, k);
4538 else if (pos_d < cell_x0[d])
4540 if (pos_d < limit0[d])
4542 cg_move_error(fplog, dd, step, cg, d, -1,
4543 cg_cm != state->x, limitd[d],
4544 cg_cm[cg], cm_new, pos_d);
4549 rvec_inc(cm_new, state->box[d]);
4552 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4553 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4555 for (k = k0; (k < k1); k++)
4557 rvec_inc(state->x[k], state->box[d]);
4560 rotate_state_atom(state, k);
4566 else if (d < npbcdim)
4568 /* Put the charge group in the rectangular unit-cell */
4569 while (cm_new[d] >= state->box[d][d])
4571 rvec_dec(cm_new, state->box[d]);
4572 for (k = k0; (k < k1); k++)
4574 rvec_dec(state->x[k], state->box[d]);
4577 while (cm_new[d] < 0)
4579 rvec_inc(cm_new, state->box[d]);
4580 for (k = k0; (k < k1); k++)
4582 rvec_inc(state->x[k], state->box[d]);
4588 copy_rvec(cm_new, cg_cm[cg]);
4590 /* Determine where this cg should go */
4593 for (d = 0; d < dd->ndim; d++)
4598 flag |= DD_FLAG_FW(d);
4604 else if (dev[dim] == -1)
4606 flag |= DD_FLAG_BW(d);
4609 if (dd->nc[dim] > 2)
4620 /* Temporarily store the flag in move */
4621 move[cg] = mc + flag;
4625 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4626 gmx_domdec_t *dd, ivec tric_dir,
4627 t_state *state, rvec **f,
4636 int ncg[DIM*2], nat[DIM*2];
4637 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4638 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4639 int sbuf[2], rbuf[2];
4640 int home_pos_cg, home_pos_at, buf_pos;
4642 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4645 real inv_ncg, pos_d;
4647 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4649 cginfo_mb_t *cginfo_mb;
4650 gmx_domdec_comm_t *comm;
4652 int nthread, thread;
4656 check_screw_box(state->box);
4660 if (fr->cutoff_scheme == ecutsGROUP)
4665 for (i = 0; i < estNR; i++)
4671 case estX: /* Always present */ break;
4672 case estV: bV = (state->flags & (1<<i)); break;
4673 case estSDX: bSDX = (state->flags & (1<<i)); break;
4674 case estCGP: bCGP = (state->flags & (1<<i)); break;
4677 case estDISRE_INITF:
4678 case estDISRE_RM3TAV:
4679 case estORIRE_INITF:
4681 /* No processing required */
4684 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4689 if (dd->ncg_tot > comm->nalloc_int)
4691 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4692 srenew(comm->buf_int, comm->nalloc_int);
4694 move = comm->buf_int;
4696 /* Clear the count */
4697 for (c = 0; c < dd->ndim*2; c++)
4703 npbcdim = dd->npbcdim;
4705 for (d = 0; (d < DIM); d++)
4707 limitd[d] = dd->comm->cellsize_min[d];
4708 if (d >= npbcdim && dd->ci[d] == 0)
4710 cell_x0[d] = -GMX_FLOAT_MAX;
4714 cell_x0[d] = comm->cell_x0[d];
4716 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4718 cell_x1[d] = GMX_FLOAT_MAX;
4722 cell_x1[d] = comm->cell_x1[d];
4726 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4727 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4731 /* We check after communication if a charge group moved
4732 * more than one cell. Set the pre-comm check limit to float_max.
4734 limit0[d] = -GMX_FLOAT_MAX;
4735 limit1[d] = GMX_FLOAT_MAX;
4739 make_tric_corr_matrix(npbcdim, state->box, tcm);
4741 cgindex = dd->cgindex;
4743 nthread = gmx_omp_nthreads_get(emntDomdec);
4745 /* Compute the center of geometry for all home charge groups
4746 * and put them in the box and determine where they should go.
4748 #pragma omp parallel for num_threads(nthread) schedule(static)
4749 for (thread = 0; thread < nthread; thread++)
4751 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4752 cell_x0, cell_x1, limitd, limit0, limit1,
4754 ( thread *dd->ncg_home)/nthread,
4755 ((thread+1)*dd->ncg_home)/nthread,
4756 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4760 for (cg = 0; cg < dd->ncg_home; cg++)
4765 flag = mc & ~DD_FLAG_NRCG;
4766 mc = mc & DD_FLAG_NRCG;
4769 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4771 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4772 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4774 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4775 /* We store the cg size in the lower 16 bits
4776 * and the place where the charge group should go
4777 * in the next 6 bits. This saves some communication volume.
4779 nrcg = cgindex[cg+1] - cgindex[cg];
4780 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4786 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4787 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4790 for (i = 0; i < dd->ndim*2; i++)
4792 *ncg_moved += ncg[i];
4809 /* Make sure the communication buffers are large enough */
4810 for (mc = 0; mc < dd->ndim*2; mc++)
4812 nvr = ncg[mc] + nat[mc]*nvec;
4813 if (nvr > comm->cgcm_state_nalloc[mc])
4815 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4816 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4820 switch (fr->cutoff_scheme)
4823 /* Recalculating cg_cm might be cheaper than communicating,
4824 * but that could give rise to rounding issues.
4827 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4828 nvec, cg_cm, comm, bCompact);
4831 /* Without charge groups we send the moved atom coordinates
4832 * over twice. This is so the code below can be used without
4833 * many conditionals for both for with and without charge groups.
4836 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4837 nvec, state->x, comm, FALSE);
4840 home_pos_cg -= *ncg_moved;
4844 gmx_incons("unimplemented");
4850 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4851 nvec, vec++, state->x, comm, bCompact);
4854 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4855 nvec, vec++, state->v, comm, bCompact);
4859 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4860 nvec, vec++, state->sd_X, comm, bCompact);
4864 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4865 nvec, vec++, state->cg_p, comm, bCompact);
4870 compact_ind(dd->ncg_home, move,
4871 dd->index_gl, dd->cgindex, dd->gatindex,
4872 dd->ga2la, comm->bLocalCG,
4877 if (fr->cutoff_scheme == ecutsVERLET)
4879 moved = get_moved(comm, dd->ncg_home);
4881 for (k = 0; k < dd->ncg_home; k++)
4888 moved = fr->ns.grid->cell_index;
4891 clear_and_mark_ind(dd->ncg_home, move,
4892 dd->index_gl, dd->cgindex, dd->gatindex,
4893 dd->ga2la, comm->bLocalCG,
4897 cginfo_mb = fr->cginfo_mb;
4899 *ncg_stay_home = home_pos_cg;
4900 for (d = 0; d < dd->ndim; d++)
4906 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4909 /* Communicate the cg and atom counts */
4914 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4915 d, dir, sbuf[0], sbuf[1]);
4917 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4919 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4921 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4922 srenew(comm->buf_int, comm->nalloc_int);
4925 /* Communicate the charge group indices, sizes and flags */
4926 dd_sendrecv_int(dd, d, dir,
4927 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4928 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4930 nvs = ncg[cdd] + nat[cdd]*nvec;
4931 i = rbuf[0] + rbuf[1] *nvec;
4932 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4934 /* Communicate cgcm and state */
4935 dd_sendrecv_rvec(dd, d, dir,
4936 comm->cgcm_state[cdd], nvs,
4937 comm->vbuf.v+nvr, i);
4938 ncg_recv += rbuf[0];
4939 nat_recv += rbuf[1];
4943 /* Process the received charge groups */
4945 for (cg = 0; cg < ncg_recv; cg++)
4947 flag = comm->buf_int[cg*DD_CGIBS+1];
4949 if (dim >= npbcdim && dd->nc[dim] > 2)
4951 /* No pbc in this dim and more than one domain boundary.
4952 * We do a separate check if a charge group didn't move too far.
4954 if (((flag & DD_FLAG_FW(d)) &&
4955 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4956 ((flag & DD_FLAG_BW(d)) &&
4957 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4959 cg_move_error(fplog, dd, step, cg, dim,
4960 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4961 fr->cutoff_scheme == ecutsGROUP, 0,
4962 comm->vbuf.v[buf_pos],
4963 comm->vbuf.v[buf_pos],
4964 comm->vbuf.v[buf_pos][dim]);
4971 /* Check which direction this cg should go */
4972 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4976 /* The cell boundaries for dimension d2 are not equal
4977 * for each cell row of the lower dimension(s),
4978 * therefore we might need to redetermine where
4979 * this cg should go.
4982 /* If this cg crosses the box boundary in dimension d2
4983 * we can use the communicated flag, so we do not
4984 * have to worry about pbc.
4986 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4987 (flag & DD_FLAG_FW(d2))) ||
4988 (dd->ci[dim2] == 0 &&
4989 (flag & DD_FLAG_BW(d2)))))
4991 /* Clear the two flags for this dimension */
4992 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4993 /* Determine the location of this cg
4994 * in lattice coordinates
4996 pos_d = comm->vbuf.v[buf_pos][dim2];
4999 for (d3 = dim2+1; d3 < DIM; d3++)
5002 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5005 /* Check of we are not at the box edge.
5006 * pbc is only handled in the first step above,
5007 * but this check could move over pbc while
5008 * the first step did not due to different rounding.
5010 if (pos_d >= cell_x1[dim2] &&
5011 dd->ci[dim2] != dd->nc[dim2]-1)
5013 flag |= DD_FLAG_FW(d2);
5015 else if (pos_d < cell_x0[dim2] &&
5018 flag |= DD_FLAG_BW(d2);
5020 comm->buf_int[cg*DD_CGIBS+1] = flag;
5023 /* Set to which neighboring cell this cg should go */
5024 if (flag & DD_FLAG_FW(d2))
5028 else if (flag & DD_FLAG_BW(d2))
5030 if (dd->nc[dd->dim[d2]] > 2)
5042 nrcg = flag & DD_FLAG_NRCG;
5045 if (home_pos_cg+1 > dd->cg_nalloc)
5047 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5048 srenew(dd->index_gl, dd->cg_nalloc);
5049 srenew(dd->cgindex, dd->cg_nalloc+1);
5051 /* Set the global charge group index and size */
5052 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5053 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5054 /* Copy the state from the buffer */
5055 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5056 if (fr->cutoff_scheme == ecutsGROUP)
5059 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5063 /* Set the cginfo */
5064 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5065 dd->index_gl[home_pos_cg]);
5068 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5071 if (home_pos_at+nrcg > state->nalloc)
5073 dd_realloc_state(state, f, home_pos_at+nrcg);
5075 for (i = 0; i < nrcg; i++)
5077 copy_rvec(comm->vbuf.v[buf_pos++],
5078 state->x[home_pos_at+i]);
5082 for (i = 0; i < nrcg; i++)
5084 copy_rvec(comm->vbuf.v[buf_pos++],
5085 state->v[home_pos_at+i]);
5090 for (i = 0; i < nrcg; i++)
5092 copy_rvec(comm->vbuf.v[buf_pos++],
5093 state->sd_X[home_pos_at+i]);
5098 for (i = 0; i < nrcg; i++)
5100 copy_rvec(comm->vbuf.v[buf_pos++],
5101 state->cg_p[home_pos_at+i]);
5105 home_pos_at += nrcg;
5109 /* Reallocate the buffers if necessary */
5110 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5112 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5113 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5115 nvr = ncg[mc] + nat[mc]*nvec;
5116 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5118 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5119 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5121 /* Copy from the receive to the send buffers */
5122 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5123 comm->buf_int + cg*DD_CGIBS,
5124 DD_CGIBS*sizeof(int));
5125 memcpy(comm->cgcm_state[mc][nvr],
5126 comm->vbuf.v[buf_pos],
5127 (1+nrcg*nvec)*sizeof(rvec));
5128 buf_pos += 1 + nrcg*nvec;
5135 /* With sorting (!bCompact) the indices are now only partially up to date
5136 * and ncg_home and nat_home are not the real count, since there are
5137 * "holes" in the arrays for the charge groups that moved to neighbors.
5139 if (fr->cutoff_scheme == ecutsVERLET)
5141 moved = get_moved(comm, home_pos_cg);
5143 for (i = dd->ncg_home; i < home_pos_cg; i++)
5148 dd->ncg_home = home_pos_cg;
5149 dd->nat_home = home_pos_at;
5154 "Finished repartitioning: cgs moved out %d, new home %d\n",
5155 *ncg_moved, dd->ncg_home-*ncg_moved);
5160 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5162 dd->comm->cycl[ddCycl] += cycles;
5163 dd->comm->cycl_n[ddCycl]++;
5164 if (cycles > dd->comm->cycl_max[ddCycl])
5166 dd->comm->cycl_max[ddCycl] = cycles;
5170 static double force_flop_count(t_nrnb *nrnb)
5177 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5179 /* To get closer to the real timings, we half the count
5180 * for the normal loops and again half it for water loops.
5183 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5185 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5189 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5192 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5195 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5197 sum += nrnb->n[i]*cost_nrnb(i);
5200 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5202 sum += nrnb->n[i]*cost_nrnb(i);
5208 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5210 if (dd->comm->eFlop)
5212 dd->comm->flop -= force_flop_count(nrnb);
5215 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5217 if (dd->comm->eFlop)
5219 dd->comm->flop += force_flop_count(nrnb);
5224 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5228 for (i = 0; i < ddCyclNr; i++)
5230 dd->comm->cycl[i] = 0;
5231 dd->comm->cycl_n[i] = 0;
5232 dd->comm->cycl_max[i] = 0;
5235 dd->comm->flop_n = 0;
5238 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5240 gmx_domdec_comm_t *comm;
5241 gmx_domdec_load_t *load;
5242 gmx_domdec_root_t *root = NULL;
5243 int d, dim, cid, i, pos;
5244 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5249 fprintf(debug, "get_load_distribution start\n");
5252 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5256 bSepPME = (dd->pme_nodeid >= 0);
5258 for (d = dd->ndim-1; d >= 0; d--)
5261 /* Check if we participate in the communication in this dimension */
5262 if (d == dd->ndim-1 ||
5263 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5265 load = &comm->load[d];
5268 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5271 if (d == dd->ndim-1)
5273 sbuf[pos++] = dd_force_load(comm);
5274 sbuf[pos++] = sbuf[0];
5277 sbuf[pos++] = sbuf[0];
5278 sbuf[pos++] = cell_frac;
5281 sbuf[pos++] = comm->cell_f_max0[d];
5282 sbuf[pos++] = comm->cell_f_min1[d];
5287 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5288 sbuf[pos++] = comm->cycl[ddCyclPME];
5293 sbuf[pos++] = comm->load[d+1].sum;
5294 sbuf[pos++] = comm->load[d+1].max;
5297 sbuf[pos++] = comm->load[d+1].sum_m;
5298 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5299 sbuf[pos++] = comm->load[d+1].flags;
5302 sbuf[pos++] = comm->cell_f_max0[d];
5303 sbuf[pos++] = comm->cell_f_min1[d];
5308 sbuf[pos++] = comm->load[d+1].mdf;
5309 sbuf[pos++] = comm->load[d+1].pme;
5313 /* Communicate a row in DD direction d.
5314 * The communicators are setup such that the root always has rank 0.
5317 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5318 load->load, load->nload*sizeof(float), MPI_BYTE,
5319 0, comm->mpi_comm_load[d]);
5321 if (dd->ci[dim] == dd->master_ci[dim])
5323 /* We are the root, process this row */
5324 if (comm->bDynLoadBal)
5326 root = comm->root[d];
5336 for (i = 0; i < dd->nc[dim]; i++)
5338 load->sum += load->load[pos++];
5339 load->max = max(load->max, load->load[pos]);
5345 /* This direction could not be load balanced properly,
5346 * therefore we need to use the maximum iso the average load.
5348 load->sum_m = max(load->sum_m, load->load[pos]);
5352 load->sum_m += load->load[pos];
5355 load->cvol_min = min(load->cvol_min, load->load[pos]);
5359 load->flags = (int)(load->load[pos++] + 0.5);
5363 root->cell_f_max0[i] = load->load[pos++];
5364 root->cell_f_min1[i] = load->load[pos++];
5369 load->mdf = max(load->mdf, load->load[pos]);
5371 load->pme = max(load->pme, load->load[pos]);
5375 if (comm->bDynLoadBal && root->bLimited)
5377 load->sum_m *= dd->nc[dim];
5378 load->flags |= (1<<d);
5386 comm->nload += dd_load_count(comm);
5387 comm->load_step += comm->cycl[ddCyclStep];
5388 comm->load_sum += comm->load[0].sum;
5389 comm->load_max += comm->load[0].max;
5390 if (comm->bDynLoadBal)
5392 for (d = 0; d < dd->ndim; d++)
5394 if (comm->load[0].flags & (1<<d))
5396 comm->load_lim[d]++;
5402 comm->load_mdf += comm->load[0].mdf;
5403 comm->load_pme += comm->load[0].pme;
5407 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5411 fprintf(debug, "get_load_distribution finished\n");
5415 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5417 /* Return the relative performance loss on the total run time
5418 * due to the force calculation load imbalance.
5420 if (dd->comm->nload > 0)
5423 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5424 (dd->comm->load_step*dd->nnodes);
5432 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5435 int npp, npme, nnodes, d, limp;
5436 float imbal, pme_f_ratio, lossf, lossp = 0;
5438 gmx_domdec_comm_t *comm;
5441 if (DDMASTER(dd) && comm->nload > 0)
5444 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5445 nnodes = npp + npme;
5446 imbal = comm->load_max*npp/comm->load_sum - 1;
5447 lossf = dd_force_imb_perf_loss(dd);
5448 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5449 fprintf(fplog, "%s", buf);
5450 fprintf(stderr, "\n");
5451 fprintf(stderr, "%s", buf);
5452 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5453 fprintf(fplog, "%s", buf);
5454 fprintf(stderr, "%s", buf);
5456 if (comm->bDynLoadBal)
5458 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5459 for (d = 0; d < dd->ndim; d++)
5461 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5462 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5468 sprintf(buf+strlen(buf), "\n");
5469 fprintf(fplog, "%s", buf);
5470 fprintf(stderr, "%s", buf);
5474 pme_f_ratio = comm->load_pme/comm->load_mdf;
5475 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5478 lossp *= (float)npme/(float)nnodes;
5482 lossp *= (float)npp/(float)nnodes;
5484 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5485 fprintf(fplog, "%s", buf);
5486 fprintf(stderr, "%s", buf);
5487 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5488 fprintf(fplog, "%s", buf);
5489 fprintf(stderr, "%s", buf);
5491 fprintf(fplog, "\n");
5492 fprintf(stderr, "\n");
5494 if (lossf >= DD_PERF_LOSS_WARN)
5497 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5498 " in the domain decomposition.\n", lossf*100);
5499 if (!comm->bDynLoadBal)
5501 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5505 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5507 fprintf(fplog, "%s\n", buf);
5508 fprintf(stderr, "%s\n", buf);
5510 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5513 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5514 " had %s work to do than the PP ranks.\n"
5515 " You might want to %s the number of PME ranks\n"
5516 " or %s the cut-off and the grid spacing.\n",
5518 (lossp < 0) ? "less" : "more",
5519 (lossp < 0) ? "decrease" : "increase",
5520 (lossp < 0) ? "decrease" : "increase");
5521 fprintf(fplog, "%s\n", buf);
5522 fprintf(stderr, "%s\n", buf);
5527 static float dd_vol_min(gmx_domdec_t *dd)
5529 return dd->comm->load[0].cvol_min*dd->nnodes;
5532 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5534 return dd->comm->load[0].flags;
5537 static float dd_f_imbal(gmx_domdec_t *dd)
5539 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5542 float dd_pme_f_ratio(gmx_domdec_t *dd)
5544 if (dd->comm->cycl_n[ddCyclPME] > 0)
5546 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5554 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5559 flags = dd_load_flags(dd);
5563 "DD load balancing is limited by minimum cell size in dimension");
5564 for (d = 0; d < dd->ndim; d++)
5568 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5571 fprintf(fplog, "\n");
5573 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5574 if (dd->comm->bDynLoadBal)
5576 fprintf(fplog, " vol min/aver %5.3f%c",
5577 dd_vol_min(dd), flags ? '!' : ' ');
5579 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5580 if (dd->comm->cycl_n[ddCyclPME])
5582 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5584 fprintf(fplog, "\n\n");
5587 static void dd_print_load_verbose(gmx_domdec_t *dd)
5589 if (dd->comm->bDynLoadBal)
5591 fprintf(stderr, "vol %4.2f%c ",
5592 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5594 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5595 if (dd->comm->cycl_n[ddCyclPME])
5597 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5602 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5607 gmx_domdec_root_t *root;
5608 gmx_bool bPartOfGroup = FALSE;
5610 dim = dd->dim[dim_ind];
5611 copy_ivec(loc, loc_c);
5612 for (i = 0; i < dd->nc[dim]; i++)
5615 rank = dd_index(dd->nc, loc_c);
5616 if (rank == dd->rank)
5618 /* This process is part of the group */
5619 bPartOfGroup = TRUE;
5622 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5626 dd->comm->mpi_comm_load[dim_ind] = c_row;
5627 if (dd->comm->eDLB != edlbNO)
5629 if (dd->ci[dim] == dd->master_ci[dim])
5631 /* This is the root process of this row */
5632 snew(dd->comm->root[dim_ind], 1);
5633 root = dd->comm->root[dim_ind];
5634 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5635 snew(root->old_cell_f, dd->nc[dim]+1);
5636 snew(root->bCellMin, dd->nc[dim]);
5639 snew(root->cell_f_max0, dd->nc[dim]);
5640 snew(root->cell_f_min1, dd->nc[dim]);
5641 snew(root->bound_min, dd->nc[dim]);
5642 snew(root->bound_max, dd->nc[dim]);
5644 snew(root->buf_ncd, dd->nc[dim]);
5648 /* This is not a root process, we only need to receive cell_f */
5649 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5652 if (dd->ci[dim] == dd->master_ci[dim])
5654 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5660 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5661 const gmx_hw_info_t gmx_unused *hwinfo,
5662 const gmx_hw_opt_t gmx_unused *hw_opt)
5665 int physicalnode_id_hash;
5668 MPI_Comm mpi_comm_pp_physicalnode;
5670 if (!(cr->duty & DUTY_PP) ||
5671 hw_opt->gpu_opt.ncuda_dev_use == 0)
5673 /* Only PP nodes (currently) use GPUs.
5674 * If we don't have GPUs, there are no resources to share.
5679 physicalnode_id_hash = gmx_physicalnode_id_hash();
5681 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5687 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5688 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5689 dd->rank, physicalnode_id_hash, gpu_id);
5691 /* Split the PP communicator over the physical nodes */
5692 /* TODO: See if we should store this (before), as it's also used for
5693 * for the nodecomm summution.
5695 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5696 &mpi_comm_pp_physicalnode);
5697 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5698 &dd->comm->mpi_comm_gpu_shared);
5699 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5700 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5704 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5707 /* Note that some ranks could share a GPU, while others don't */
5709 if (dd->comm->nrank_gpu_shared == 1)
5711 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5716 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5719 int dim0, dim1, i, j;
5724 fprintf(debug, "Making load communicators\n");
5727 snew(dd->comm->load, dd->ndim);
5728 snew(dd->comm->mpi_comm_load, dd->ndim);
5731 make_load_communicator(dd, 0, loc);
5735 for (i = 0; i < dd->nc[dim0]; i++)
5738 make_load_communicator(dd, 1, loc);
5744 for (i = 0; i < dd->nc[dim0]; i++)
5748 for (j = 0; j < dd->nc[dim1]; j++)
5751 make_load_communicator(dd, 2, loc);
5758 fprintf(debug, "Finished making load communicators\n");
5763 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5766 int d, dim, i, j, m;
5769 ivec dd_zp[DD_MAXIZONE];
5770 gmx_domdec_zones_t *zones;
5771 gmx_domdec_ns_ranges_t *izone;
5773 for (d = 0; d < dd->ndim; d++)
5776 copy_ivec(dd->ci, tmp);
5777 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5778 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5779 copy_ivec(dd->ci, tmp);
5780 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5781 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5784 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5787 dd->neighbor[d][1]);
5793 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5795 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5796 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5803 for (i = 0; i < nzonep; i++)
5805 copy_ivec(dd_zp3[i], dd_zp[i]);
5811 for (i = 0; i < nzonep; i++)
5813 copy_ivec(dd_zp2[i], dd_zp[i]);
5819 for (i = 0; i < nzonep; i++)
5821 copy_ivec(dd_zp1[i], dd_zp[i]);
5825 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5830 zones = &dd->comm->zones;
5832 for (i = 0; i < nzone; i++)
5835 clear_ivec(zones->shift[i]);
5836 for (d = 0; d < dd->ndim; d++)
5838 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5843 for (i = 0; i < nzone; i++)
5845 for (d = 0; d < DIM; d++)
5847 s[d] = dd->ci[d] - zones->shift[i][d];
5852 else if (s[d] >= dd->nc[d])
5858 zones->nizone = nzonep;
5859 for (i = 0; i < zones->nizone; i++)
5861 if (dd_zp[i][0] != i)
5863 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5865 izone = &zones->izone[i];
5866 izone->j0 = dd_zp[i][1];
5867 izone->j1 = dd_zp[i][2];
5868 for (dim = 0; dim < DIM; dim++)
5870 if (dd->nc[dim] == 1)
5872 /* All shifts should be allowed */
5873 izone->shift0[dim] = -1;
5874 izone->shift1[dim] = 1;
5879 izone->shift0[d] = 0;
5880 izone->shift1[d] = 0;
5881 for(j=izone->j0; j<izone->j1; j++) {
5882 if (dd->shift[j][d] > dd->shift[i][d])
5883 izone->shift0[d] = -1;
5884 if (dd->shift[j][d] < dd->shift[i][d])
5885 izone->shift1[d] = 1;
5891 /* Assume the shift are not more than 1 cell */
5892 izone->shift0[dim] = 1;
5893 izone->shift1[dim] = -1;
5894 for (j = izone->j0; j < izone->j1; j++)
5896 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5897 if (shift_diff < izone->shift0[dim])
5899 izone->shift0[dim] = shift_diff;
5901 if (shift_diff > izone->shift1[dim])
5903 izone->shift1[dim] = shift_diff;
5910 if (dd->comm->eDLB != edlbNO)
5912 snew(dd->comm->root, dd->ndim);
5915 if (dd->comm->bRecordLoad)
5917 make_load_communicators(dd);
5921 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5924 gmx_domdec_comm_t *comm;
5935 if (comm->bCartesianPP)
5937 /* Set up cartesian communication for the particle-particle part */
5940 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5941 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5944 for (i = 0; i < DIM; i++)
5948 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5950 /* We overwrite the old communicator with the new cartesian one */
5951 cr->mpi_comm_mygroup = comm_cart;
5954 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5955 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5957 if (comm->bCartesianPP_PME)
5959 /* Since we want to use the original cartesian setup for sim,
5960 * and not the one after split, we need to make an index.
5962 snew(comm->ddindex2ddnodeid, dd->nnodes);
5963 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5964 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5965 /* Get the rank of the DD master,
5966 * above we made sure that the master node is a PP node.
5976 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5978 else if (comm->bCartesianPP)
5980 if (cr->npmenodes == 0)
5982 /* The PP communicator is also
5983 * the communicator for this simulation
5985 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5987 cr->nodeid = dd->rank;
5989 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5991 /* We need to make an index to go from the coordinates
5992 * to the nodeid of this simulation.
5994 snew(comm->ddindex2simnodeid, dd->nnodes);
5995 snew(buf, dd->nnodes);
5996 if (cr->duty & DUTY_PP)
5998 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6000 /* Communicate the ddindex to simulation nodeid index */
6001 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6002 cr->mpi_comm_mysim);
6005 /* Determine the master coordinates and rank.
6006 * The DD master should be the same node as the master of this sim.
6008 for (i = 0; i < dd->nnodes; i++)
6010 if (comm->ddindex2simnodeid[i] == 0)
6012 ddindex2xyz(dd->nc, i, dd->master_ci);
6013 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6018 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6023 /* No Cartesian communicators */
6024 /* We use the rank in dd->comm->all as DD index */
6025 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6026 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6028 clear_ivec(dd->master_ci);
6035 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6036 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6041 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6042 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6046 static void receive_ddindex2simnodeid(t_commrec *cr)
6050 gmx_domdec_comm_t *comm;
6057 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6059 snew(comm->ddindex2simnodeid, dd->nnodes);
6060 snew(buf, dd->nnodes);
6061 if (cr->duty & DUTY_PP)
6063 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6066 /* Communicate the ddindex to simulation nodeid index */
6067 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6068 cr->mpi_comm_mysim);
6075 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6076 int ncg, int natoms)
6078 gmx_domdec_master_t *ma;
6083 snew(ma->ncg, dd->nnodes);
6084 snew(ma->index, dd->nnodes+1);
6086 snew(ma->nat, dd->nnodes);
6087 snew(ma->ibuf, dd->nnodes*2);
6088 snew(ma->cell_x, DIM);
6089 for (i = 0; i < DIM; i++)
6091 snew(ma->cell_x[i], dd->nc[i]+1);
6094 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6100 snew(ma->vbuf, natoms);
6106 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6107 int gmx_unused reorder)
6110 gmx_domdec_comm_t *comm;
6121 if (comm->bCartesianPP)
6123 for (i = 1; i < DIM; i++)
6125 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6127 if (bDiv[YY] || bDiv[ZZ])
6129 comm->bCartesianPP_PME = TRUE;
6130 /* If we have 2D PME decomposition, which is always in x+y,
6131 * we stack the PME only nodes in z.
6132 * Otherwise we choose the direction that provides the thinnest slab
6133 * of PME only nodes as this will have the least effect
6134 * on the PP communication.
6135 * But for the PME communication the opposite might be better.
6137 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6139 dd->nc[YY] > dd->nc[ZZ]))
6141 comm->cartpmedim = ZZ;
6145 comm->cartpmedim = YY;
6147 comm->ntot[comm->cartpmedim]
6148 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6152 fprintf(fplog, "Number of PME-only ranks (%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]);
6154 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6159 if (comm->bCartesianPP_PME)
6163 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]);
6166 for (i = 0; i < DIM; i++)
6170 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6173 MPI_Comm_rank(comm_cart, &rank);
6174 if (MASTERNODE(cr) && rank != 0)
6176 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6179 /* With this assigment we loose the link to the original communicator
6180 * which will usually be MPI_COMM_WORLD, unless have multisim.
6182 cr->mpi_comm_mysim = comm_cart;
6183 cr->sim_nodeid = rank;
6185 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6189 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6190 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6193 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6197 if (cr->npmenodes == 0 ||
6198 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6200 cr->duty = DUTY_PME;
6203 /* Split the sim communicator into PP and PME only nodes */
6204 MPI_Comm_split(cr->mpi_comm_mysim,
6206 dd_index(comm->ntot, dd->ci),
6207 &cr->mpi_comm_mygroup);
6211 switch (dd_node_order)
6216 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6219 case ddnoINTERLEAVE:
6220 /* Interleave the PP-only and PME-only nodes,
6221 * as on clusters with dual-core machines this will double
6222 * the communication bandwidth of the PME processes
6223 * and thus speed up the PP <-> PME and inter PME communication.
6227 fprintf(fplog, "Interleaving PP and PME ranks\n");
6229 comm->pmenodes = dd_pmenodes(cr);
6234 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6237 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6239 cr->duty = DUTY_PME;
6246 /* Split the sim communicator into PP and PME only nodes */
6247 MPI_Comm_split(cr->mpi_comm_mysim,
6250 &cr->mpi_comm_mygroup);
6251 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6257 fprintf(fplog, "This rank does only %s work.\n\n",
6258 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6262 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6265 gmx_domdec_comm_t *comm;
6271 copy_ivec(dd->nc, comm->ntot);
6273 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6274 comm->bCartesianPP_PME = FALSE;
6276 /* Reorder the nodes by default. This might change the MPI ranks.
6277 * Real reordering is only supported on very few architectures,
6278 * Blue Gene is one of them.
6280 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6282 if (cr->npmenodes > 0)
6284 /* Split the communicator into a PP and PME part */
6285 split_communicator(fplog, cr, dd_node_order, CartReorder);
6286 if (comm->bCartesianPP_PME)
6288 /* We (possibly) reordered the nodes in split_communicator,
6289 * so it is no longer required in make_pp_communicator.
6291 CartReorder = FALSE;
6296 /* All nodes do PP and PME */
6298 /* We do not require separate communicators */
6299 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6303 if (cr->duty & DUTY_PP)
6305 /* Copy or make a new PP communicator */
6306 make_pp_communicator(fplog, cr, CartReorder);
6310 receive_ddindex2simnodeid(cr);
6313 if (!(cr->duty & DUTY_PME))
6315 /* Set up the commnuication to our PME node */
6316 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6317 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6320 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6321 dd->pme_nodeid, dd->pme_receive_vir_ener);
6326 dd->pme_nodeid = -1;
6331 dd->ma = init_gmx_domdec_master_t(dd,
6333 comm->cgs_gl.index[comm->cgs_gl.nr]);
6337 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6339 real *slb_frac, tot;
6344 if (nc > 1 && size_string != NULL)
6348 fprintf(fplog, "Using static load balancing for the %s direction\n",
6353 for (i = 0; i < nc; i++)
6356 sscanf(size_string, "%lf%n", &dbl, &n);
6359 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6368 fprintf(fplog, "Relative cell sizes:");
6370 for (i = 0; i < nc; i++)
6375 fprintf(fplog, " %5.3f", slb_frac[i]);
6380 fprintf(fplog, "\n");
6387 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6390 gmx_mtop_ilistloop_t iloop;
6394 iloop = gmx_mtop_ilistloop_init(mtop);
6395 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6397 for (ftype = 0; ftype < F_NRE; ftype++)
6399 if ((interaction_function[ftype].flags & IF_BOND) &&
6402 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6410 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6416 val = getenv(env_var);
6419 if (sscanf(val, "%d", &nst) <= 0)
6425 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6433 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6437 fprintf(stderr, "\n%s\n", warn_string);
6441 fprintf(fplog, "\n%s\n", warn_string);
6445 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6446 t_inputrec *ir, FILE *fplog)
6448 if (ir->ePBC == epbcSCREW &&
6449 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6451 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6454 if (ir->ns_type == ensSIMPLE)
6456 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6459 if (ir->nstlist == 0)
6461 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6464 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6466 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6470 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6475 r = ddbox->box_size[XX];
6476 for (di = 0; di < dd->ndim; di++)
6479 /* Check using the initial average cell size */
6480 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6486 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6487 const char *dlb_opt, gmx_bool bRecordLoad,
6488 unsigned long Flags, t_inputrec *ir)
6496 case 'a': eDLB = edlbAUTO; break;
6497 case 'n': eDLB = edlbNO; break;
6498 case 'y': eDLB = edlbYES; break;
6499 default: gmx_incons("Unknown dlb_opt");
6502 if (Flags & MD_RERUN)
6507 if (!EI_DYNAMICS(ir->eI))
6509 if (eDLB == edlbYES)
6511 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6512 dd_warning(cr, fplog, buf);
6520 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6525 if (Flags & MD_REPRODUCIBLE)
6532 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6536 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6539 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6547 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6552 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6554 /* Decomposition order z,y,x */
6557 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6559 for (dim = DIM-1; dim >= 0; dim--)
6561 if (dd->nc[dim] > 1)
6563 dd->dim[dd->ndim++] = dim;
6569 /* Decomposition order x,y,z */
6570 for (dim = 0; dim < DIM; dim++)
6572 if (dd->nc[dim] > 1)
6574 dd->dim[dd->ndim++] = dim;
6580 static gmx_domdec_comm_t *init_dd_comm()
6582 gmx_domdec_comm_t *comm;
6586 snew(comm->cggl_flag, DIM*2);
6587 snew(comm->cgcm_state, DIM*2);
6588 for (i = 0; i < DIM*2; i++)
6590 comm->cggl_flag_nalloc[i] = 0;
6591 comm->cgcm_state_nalloc[i] = 0;
6594 comm->nalloc_int = 0;
6595 comm->buf_int = NULL;
6597 vec_rvec_init(&comm->vbuf);
6599 comm->n_load_have = 0;
6600 comm->n_load_collect = 0;
6602 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6604 comm->sum_nat[i] = 0;
6608 comm->load_step = 0;
6611 clear_ivec(comm->load_lim);
6618 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6619 unsigned long Flags,
6621 real comm_distance_min, real rconstr,
6622 const char *dlb_opt, real dlb_scale,
6623 const char *sizex, const char *sizey, const char *sizez,
6624 gmx_mtop_t *mtop, t_inputrec *ir,
6625 matrix box, rvec *x,
6627 int *npme_x, int *npme_y)
6630 gmx_domdec_comm_t *comm;
6633 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6640 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6645 dd->comm = init_dd_comm();
6647 snew(comm->cggl_flag, DIM*2);
6648 snew(comm->cgcm_state, DIM*2);
6650 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6651 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6653 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6654 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6655 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6656 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6657 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6658 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6659 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6660 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6662 dd->pme_recv_f_alloc = 0;
6663 dd->pme_recv_f_buf = NULL;
6665 if (dd->bSendRecv2 && fplog)
6667 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");
6673 fprintf(fplog, "Will load balance based on FLOP count\n");
6675 if (comm->eFlop > 1)
6677 srand(1+cr->nodeid);
6679 comm->bRecordLoad = TRUE;
6683 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6687 /* Initialize to GPU share count to 0, might change later */
6688 comm->nrank_gpu_shared = 0;
6690 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6691 comm->bDLB_locked = FALSE;
6693 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6696 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6698 dd->bGridJump = comm->bDynLoadBal;
6699 comm->bPMELoadBalDLBLimits = FALSE;
6701 if (comm->nstSortCG)
6705 if (comm->nstSortCG == 1)
6707 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6711 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6715 snew(comm->sort, 1);
6721 fprintf(fplog, "Will not sort the charge groups\n");
6725 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6727 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6728 if (comm->bInterCGBondeds)
6730 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6734 comm->bInterCGMultiBody = FALSE;
6737 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6738 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6740 if (ir->rlistlong == 0)
6742 /* Set the cut-off to some very large value,
6743 * so we don't need if statements everywhere in the code.
6744 * We use sqrt, since the cut-off is squared in some places.
6746 comm->cutoff = GMX_CUTOFF_INF;
6750 comm->cutoff = ir->rlistlong;
6752 comm->cutoff_mbody = 0;
6754 comm->cellsize_limit = 0;
6755 comm->bBondComm = FALSE;
6757 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6758 * within nstlist steps. Since boundaries are allowed to displace by half
6759 * a cell size, DD cells should be at least the size of the list buffer.
6761 comm->cellsize_limit = max(comm->cellsize_limit,
6762 ir->rlistlong - max(ir->rvdw, ir->rcoulomb));
6764 if (comm->bInterCGBondeds)
6766 if (comm_distance_min > 0)
6768 comm->cutoff_mbody = comm_distance_min;
6769 if (Flags & MD_DDBONDCOMM)
6771 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6775 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6777 r_bonded_limit = comm->cutoff_mbody;
6779 else if (ir->bPeriodicMols)
6781 /* Can not easily determine the required cut-off */
6782 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");
6783 comm->cutoff_mbody = comm->cutoff/2;
6784 r_bonded_limit = comm->cutoff_mbody;
6790 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6791 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6793 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6794 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6796 /* We use an initial margin of 10% for the minimum cell size,
6797 * except when we are just below the non-bonded cut-off.
6799 if (Flags & MD_DDBONDCOMM)
6801 if (max(r_2b, r_mb) > comm->cutoff)
6803 r_bonded = max(r_2b, r_mb);
6804 r_bonded_limit = 1.1*r_bonded;
6805 comm->bBondComm = TRUE;
6810 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6812 /* We determine cutoff_mbody later */
6816 /* No special bonded communication,
6817 * simply increase the DD cut-off.
6819 r_bonded_limit = 1.1*max(r_2b, r_mb);
6820 comm->cutoff_mbody = r_bonded_limit;
6821 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6824 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6828 "Minimum cell size due to bonded interactions: %.3f nm\n",
6829 comm->cellsize_limit);
6833 if (dd->bInterCGcons && rconstr <= 0)
6835 /* There is a cell size limit due to the constraints (P-LINCS) */
6836 rconstr = constr_r_max(fplog, mtop, ir);
6840 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6842 if (rconstr > comm->cellsize_limit)
6844 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6848 else if (rconstr > 0 && fplog)
6850 /* Here we do not check for dd->bInterCGcons,
6851 * because one can also set a cell size limit for virtual sites only
6852 * and at this point we don't know yet if there are intercg v-sites.
6855 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6858 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6860 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6864 copy_ivec(nc, dd->nc);
6865 set_dd_dim(fplog, dd);
6866 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6868 if (cr->npmenodes == -1)
6872 acs = average_cellsize_min(dd, ddbox);
6873 if (acs < comm->cellsize_limit)
6877 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6879 gmx_fatal_collective(FARGS, cr, NULL,
6880 "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",
6881 acs, comm->cellsize_limit);
6886 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6888 /* We need to choose the optimal DD grid and possibly PME nodes */
6889 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6890 comm->eDLB != edlbNO, dlb_scale,
6891 comm->cellsize_limit, comm->cutoff,
6892 comm->bInterCGBondeds);
6894 if (dd->nc[XX] == 0)
6896 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6897 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6898 !bC ? "-rdd" : "-rcon",
6899 comm->eDLB != edlbNO ? " or -dds" : "",
6900 bC ? " or your LINCS settings" : "");
6902 gmx_fatal_collective(FARGS, cr, NULL,
6903 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6905 "Look in the log file for details on the domain decomposition",
6906 cr->nnodes-cr->npmenodes, limit, buf);
6908 set_dd_dim(fplog, dd);
6914 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6915 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6918 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6919 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6921 gmx_fatal_collective(FARGS, cr, NULL,
6922 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6923 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6925 if (cr->npmenodes > dd->nnodes)
6927 gmx_fatal_collective(FARGS, cr, NULL,
6928 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6930 if (cr->npmenodes > 0)
6932 comm->npmenodes = cr->npmenodes;
6936 comm->npmenodes = dd->nnodes;
6939 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6941 /* The following choices should match those
6942 * in comm_cost_est in domdec_setup.c.
6943 * Note that here the checks have to take into account
6944 * that the decomposition might occur in a different order than xyz
6945 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6946 * in which case they will not match those in comm_cost_est,
6947 * but since that is mainly for testing purposes that's fine.
6949 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6950 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6951 getenv("GMX_PMEONEDD") == NULL)
6953 comm->npmedecompdim = 2;
6954 comm->npmenodes_x = dd->nc[XX];
6955 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6959 /* In case nc is 1 in both x and y we could still choose to
6960 * decompose pme in y instead of x, but we use x for simplicity.
6962 comm->npmedecompdim = 1;
6963 if (dd->dim[0] == YY)
6965 comm->npmenodes_x = 1;
6966 comm->npmenodes_y = comm->npmenodes;
6970 comm->npmenodes_x = comm->npmenodes;
6971 comm->npmenodes_y = 1;
6976 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6977 comm->npmenodes_x, comm->npmenodes_y, 1);
6982 comm->npmedecompdim = 0;
6983 comm->npmenodes_x = 0;
6984 comm->npmenodes_y = 0;
6987 /* Technically we don't need both of these,
6988 * but it simplifies code not having to recalculate it.
6990 *npme_x = comm->npmenodes_x;
6991 *npme_y = comm->npmenodes_y;
6993 snew(comm->slb_frac, DIM);
6994 if (comm->eDLB == edlbNO)
6996 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6997 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6998 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7001 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7003 if (comm->bBondComm || comm->eDLB != edlbNO)
7005 /* Set the bonded communication distance to halfway
7006 * the minimum and the maximum,
7007 * since the extra communication cost is nearly zero.
7009 acs = average_cellsize_min(dd, ddbox);
7010 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7011 if (comm->eDLB != edlbNO)
7013 /* Check if this does not limit the scaling */
7014 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
7016 if (!comm->bBondComm)
7018 /* Without bBondComm do not go beyond the n.b. cut-off */
7019 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
7020 if (comm->cellsize_limit >= comm->cutoff)
7022 /* We don't loose a lot of efficieny
7023 * when increasing it to the n.b. cut-off.
7024 * It can even be slightly faster, because we need
7025 * less checks for the communication setup.
7027 comm->cutoff_mbody = comm->cutoff;
7030 /* Check if we did not end up below our original limit */
7031 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
7033 if (comm->cutoff_mbody > comm->cellsize_limit)
7035 comm->cellsize_limit = comm->cutoff_mbody;
7038 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7043 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7044 "cellsize limit %f\n",
7045 comm->bBondComm, comm->cellsize_limit);
7050 check_dd_restrictions(cr, dd, ir, fplog);
7053 comm->partition_step = INT_MIN;
7056 clear_dd_cycle_counts(dd);
7061 static void set_dlb_limits(gmx_domdec_t *dd)
7066 for (d = 0; d < dd->ndim; d++)
7068 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7069 dd->comm->cellsize_min[dd->dim[d]] =
7070 dd->comm->cellsize_min_dlb[dd->dim[d]];
7075 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7078 gmx_domdec_comm_t *comm;
7088 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);
7091 cellsize_min = comm->cellsize_min[dd->dim[0]];
7092 for (d = 1; d < dd->ndim; d++)
7094 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7097 if (cellsize_min < comm->cellsize_limit*1.05)
7099 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");
7101 /* Change DLB from "auto" to "no". */
7102 comm->eDLB = edlbNO;
7107 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7108 comm->bDynLoadBal = TRUE;
7109 dd->bGridJump = TRUE;
7113 /* We can set the required cell size info here,
7114 * so we do not need to communicate this.
7115 * The grid is completely uniform.
7117 for (d = 0; d < dd->ndim; d++)
7121 comm->load[d].sum_m = comm->load[d].sum;
7123 nc = dd->nc[dd->dim[d]];
7124 for (i = 0; i < nc; i++)
7126 comm->root[d]->cell_f[i] = i/(real)nc;
7129 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7130 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7133 comm->root[d]->cell_f[nc] = 1.0;
7138 static char *init_bLocalCG(gmx_mtop_t *mtop)
7143 ncg = ncg_mtop(mtop);
7144 snew(bLocalCG, ncg);
7145 for (cg = 0; cg < ncg; cg++)
7147 bLocalCG[cg] = FALSE;
7153 void dd_init_bondeds(FILE *fplog,
7154 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7156 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7158 gmx_domdec_comm_t *comm;
7162 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7166 if (comm->bBondComm)
7168 /* Communicate atoms beyond the cut-off for bonded interactions */
7171 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7173 comm->bLocalCG = init_bLocalCG(mtop);
7177 /* Only communicate atoms based on cut-off */
7178 comm->cglink = NULL;
7179 comm->bLocalCG = NULL;
7183 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7185 gmx_bool bDynLoadBal, real dlb_scale,
7188 gmx_domdec_comm_t *comm;
7203 fprintf(fplog, "The maximum number of communication pulses is:");
7204 for (d = 0; d < dd->ndim; d++)
7206 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7208 fprintf(fplog, "\n");
7209 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7210 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7211 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7212 for (d = 0; d < DIM; d++)
7216 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7223 comm->cellsize_min_dlb[d]/
7224 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7226 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7229 fprintf(fplog, "\n");
7233 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7234 fprintf(fplog, "The initial number of communication pulses is:");
7235 for (d = 0; d < dd->ndim; d++)
7237 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7239 fprintf(fplog, "\n");
7240 fprintf(fplog, "The initial domain decomposition cell size is:");
7241 for (d = 0; d < DIM; d++)
7245 fprintf(fplog, " %c %.2f nm",
7246 dim2char(d), dd->comm->cellsize_min[d]);
7249 fprintf(fplog, "\n\n");
7252 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7254 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7255 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7256 "non-bonded interactions", "", comm->cutoff);
7260 limit = dd->comm->cellsize_limit;
7264 if (dynamic_dd_box(ddbox, ir))
7266 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7268 limit = dd->comm->cellsize_min[XX];
7269 for (d = 1; d < DIM; d++)
7271 limit = min(limit, dd->comm->cellsize_min[d]);
7275 if (comm->bInterCGBondeds)
7277 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7278 "two-body bonded interactions", "(-rdd)",
7279 max(comm->cutoff, comm->cutoff_mbody));
7280 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7281 "multi-body bonded interactions", "(-rdd)",
7282 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7286 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7287 "virtual site constructions", "(-rcon)", limit);
7289 if (dd->constraint_comm)
7291 sprintf(buf, "atoms separated by up to %d constraints",
7293 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7294 buf, "(-rcon)", limit);
7296 fprintf(fplog, "\n");
7302 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7304 const t_inputrec *ir,
7305 const gmx_ddbox_t *ddbox)
7307 gmx_domdec_comm_t *comm;
7308 int d, dim, npulse, npulse_d_max, npulse_d;
7313 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7315 /* Determine the maximum number of comm. pulses in one dimension */
7317 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7319 /* Determine the maximum required number of grid pulses */
7320 if (comm->cellsize_limit >= comm->cutoff)
7322 /* Only a single pulse is required */
7325 else if (!bNoCutOff && comm->cellsize_limit > 0)
7327 /* We round down slightly here to avoid overhead due to the latency
7328 * of extra communication calls when the cut-off
7329 * would be only slightly longer than the cell size.
7330 * Later cellsize_limit is redetermined,
7331 * so we can not miss interactions due to this rounding.
7333 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7337 /* There is no cell size limit */
7338 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7341 if (!bNoCutOff && npulse > 1)
7343 /* See if we can do with less pulses, based on dlb_scale */
7345 for (d = 0; d < dd->ndim; d++)
7348 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7349 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7350 npulse_d_max = max(npulse_d_max, npulse_d);
7352 npulse = min(npulse, npulse_d_max);
7355 /* This env var can override npulse */
7356 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7363 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7364 for (d = 0; d < dd->ndim; d++)
7366 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7367 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7368 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7369 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7370 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7372 comm->bVacDLBNoLimit = FALSE;
7376 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7377 if (!comm->bVacDLBNoLimit)
7379 comm->cellsize_limit = max(comm->cellsize_limit,
7380 comm->cutoff/comm->maxpulse);
7382 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7383 /* Set the minimum cell size for each DD dimension */
7384 for (d = 0; d < dd->ndim; d++)
7386 if (comm->bVacDLBNoLimit ||
7387 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7389 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7393 comm->cellsize_min_dlb[dd->dim[d]] =
7394 comm->cutoff/comm->cd[d].np_dlb;
7397 if (comm->cutoff_mbody <= 0)
7399 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7401 if (comm->bDynLoadBal)
7407 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7409 /* If each molecule is a single charge group
7410 * or we use domain decomposition for each periodic dimension,
7411 * we do not need to take pbc into account for the bonded interactions.
7413 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7416 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7419 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7420 t_inputrec *ir, gmx_ddbox_t *ddbox)
7422 gmx_domdec_comm_t *comm;
7428 /* Initialize the thread data.
7429 * This can not be done in init_domain_decomposition,
7430 * as the numbers of threads is determined later.
7432 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7435 snew(comm->dth, comm->nth);
7438 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7440 init_ddpme(dd, &comm->ddpme[0], 0);
7441 if (comm->npmedecompdim >= 2)
7443 init_ddpme(dd, &comm->ddpme[1], 1);
7448 comm->npmenodes = 0;
7449 if (dd->pme_nodeid >= 0)
7451 gmx_fatal_collective(FARGS, NULL, dd,
7452 "Can not have separate PME ranks without PME electrostatics");
7458 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7460 if (comm->eDLB != edlbNO)
7462 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7465 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7466 if (comm->eDLB == edlbAUTO)
7470 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7472 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7475 if (ir->ePBC == epbcNONE)
7477 vol_frac = 1 - 1/(double)dd->nnodes;
7482 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7486 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7488 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7490 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7493 static gmx_bool test_dd_cutoff(t_commrec *cr,
7494 t_state *state, t_inputrec *ir,
7505 set_ddbox(dd, FALSE, cr, ir, state->box,
7506 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7510 for (d = 0; d < dd->ndim; d++)
7514 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7515 if (dynamic_dd_box(&ddbox, ir))
7517 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7520 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7522 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7523 dd->comm->cd[d].np_dlb > 0)
7525 if (np > dd->comm->cd[d].np_dlb)
7530 /* If a current local cell size is smaller than the requested
7531 * cut-off, we could still fix it, but this gets very complicated.
7532 * Without fixing here, we might actually need more checks.
7534 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7541 if (dd->comm->eDLB != edlbNO)
7543 /* If DLB is not active yet, we don't need to check the grid jumps.
7544 * Actually we shouldn't, because then the grid jump data is not set.
7546 if (dd->comm->bDynLoadBal &&
7547 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7552 gmx_sumi(1, &LocallyLimited, cr);
7554 if (LocallyLimited > 0)
7563 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7566 gmx_bool bCutoffAllowed;
7568 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7572 cr->dd->comm->cutoff = cutoff_req;
7575 return bCutoffAllowed;
7578 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7580 gmx_domdec_comm_t *comm;
7582 comm = cr->dd->comm;
7584 /* Turn on the DLB limiting (might have been on already) */
7585 comm->bPMELoadBalDLBLimits = TRUE;
7587 /* Change the cut-off limit */
7588 comm->PMELoadBal_max_cutoff = comm->cutoff;
7591 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7593 return dd->comm->bDLB_locked;
7596 void dd_dlb_set_lock(gmx_domdec_t *dd, gmx_bool bValue)
7598 /* We can only lock the DLB when it is set to auto, otherwise don't lock */
7599 if (dd->comm->eDLB == edlbAUTO)
7601 dd->comm->bDLB_locked = bValue;
7605 static void merge_cg_buffers(int ncell,
7606 gmx_domdec_comm_dim_t *cd, int pulse,
7608 int *index_gl, int *recv_i,
7609 rvec *cg_cm, rvec *recv_vr,
7611 cginfo_mb_t *cginfo_mb, int *cginfo)
7613 gmx_domdec_ind_t *ind, *ind_p;
7614 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7615 int shift, shift_at;
7617 ind = &cd->ind[pulse];
7619 /* First correct the already stored data */
7620 shift = ind->nrecv[ncell];
7621 for (cell = ncell-1; cell >= 0; cell--)
7623 shift -= ind->nrecv[cell];
7626 /* Move the cg's present from previous grid pulses */
7627 cg0 = ncg_cell[ncell+cell];
7628 cg1 = ncg_cell[ncell+cell+1];
7629 cgindex[cg1+shift] = cgindex[cg1];
7630 for (cg = cg1-1; cg >= cg0; cg--)
7632 index_gl[cg+shift] = index_gl[cg];
7633 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7634 cgindex[cg+shift] = cgindex[cg];
7635 cginfo[cg+shift] = cginfo[cg];
7637 /* Correct the already stored send indices for the shift */
7638 for (p = 1; p <= pulse; p++)
7640 ind_p = &cd->ind[p];
7642 for (c = 0; c < cell; c++)
7644 cg0 += ind_p->nsend[c];
7646 cg1 = cg0 + ind_p->nsend[cell];
7647 for (cg = cg0; cg < cg1; cg++)
7649 ind_p->index[cg] += shift;
7655 /* Merge in the communicated buffers */
7659 for (cell = 0; cell < ncell; cell++)
7661 cg1 = ncg_cell[ncell+cell+1] + shift;
7664 /* Correct the old cg indices */
7665 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7667 cgindex[cg+1] += shift_at;
7670 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7672 /* Copy this charge group from the buffer */
7673 index_gl[cg1] = recv_i[cg0];
7674 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7675 /* Add it to the cgindex */
7676 cg_gl = index_gl[cg1];
7677 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7678 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7679 cgindex[cg1+1] = cgindex[cg1] + nat;
7684 shift += ind->nrecv[cell];
7685 ncg_cell[ncell+cell+1] = cg1;
7689 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7690 int nzone, int cg0, const int *cgindex)
7694 /* Store the atom block boundaries for easy copying of communication buffers
7697 for (zone = 0; zone < nzone; zone++)
7699 for (p = 0; p < cd->np; p++)
7701 cd->ind[p].cell2at0[zone] = cgindex[cg];
7702 cg += cd->ind[p].nrecv[zone];
7703 cd->ind[p].cell2at1[zone] = cgindex[cg];
7708 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7714 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7716 if (!bLocalCG[link->a[i]])
7725 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7727 real c[DIM][4]; /* the corners for the non-bonded communication */
7728 real cr0; /* corner for rounding */
7729 real cr1[4]; /* corners for rounding */
7730 real bc[DIM]; /* corners for bounded communication */
7731 real bcr1; /* corner for rounding for bonded communication */
7734 /* Determine the corners of the domain(s) we are communicating with */
7736 set_dd_corners(const gmx_domdec_t *dd,
7737 int dim0, int dim1, int dim2,
7741 const gmx_domdec_comm_t *comm;
7742 const gmx_domdec_zones_t *zones;
7747 zones = &comm->zones;
7749 /* Keep the compiler happy */
7753 /* The first dimension is equal for all cells */
7754 c->c[0][0] = comm->cell_x0[dim0];
7757 c->bc[0] = c->c[0][0];
7762 /* This cell row is only seen from the first row */
7763 c->c[1][0] = comm->cell_x0[dim1];
7764 /* All rows can see this row */
7765 c->c[1][1] = comm->cell_x0[dim1];
7768 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7771 /* For the multi-body distance we need the maximum */
7772 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7775 /* Set the upper-right corner for rounding */
7776 c->cr0 = comm->cell_x1[dim0];
7781 for (j = 0; j < 4; j++)
7783 c->c[2][j] = comm->cell_x0[dim2];
7787 /* Use the maximum of the i-cells that see a j-cell */
7788 for (i = 0; i < zones->nizone; i++)
7790 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7796 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7802 /* For the multi-body distance we need the maximum */
7803 c->bc[2] = comm->cell_x0[dim2];
7804 for (i = 0; i < 2; i++)
7806 for (j = 0; j < 2; j++)
7808 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7814 /* Set the upper-right corner for rounding */
7815 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7816 * Only cell (0,0,0) can see cell 7 (1,1,1)
7818 c->cr1[0] = comm->cell_x1[dim1];
7819 c->cr1[3] = comm->cell_x1[dim1];
7822 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7825 /* For the multi-body distance we need the maximum */
7826 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7833 /* Determine which cg's we need to send in this pulse from this zone */
7835 get_zone_pulse_cgs(gmx_domdec_t *dd,
7836 int zonei, int zone,
7838 const int *index_gl,
7840 int dim, int dim_ind,
7841 int dim0, int dim1, int dim2,
7842 real r_comm2, real r_bcomm2,
7846 real skew_fac2_d, real skew_fac_01,
7847 rvec *v_d, rvec *v_0, rvec *v_1,
7848 const dd_corners_t *c,
7850 gmx_bool bDistBonded,
7856 gmx_domdec_ind_t *ind,
7857 int **ibuf, int *ibuf_nalloc,
7863 gmx_domdec_comm_t *comm;
7865 gmx_bool bDistMB_pulse;
7867 real r2, rb2, r, tric_sh;
7870 int nsend_z, nsend, nat;
7874 bScrew = (dd->bScrewPBC && dim == XX);
7876 bDistMB_pulse = (bDistMB && bDistBonded);
7882 for (cg = cg0; cg < cg1; cg++)
7886 if (tric_dist[dim_ind] == 0)
7888 /* Rectangular direction, easy */
7889 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7896 r = cg_cm[cg][dim] - c->bc[dim_ind];
7902 /* Rounding gives at most a 16% reduction
7903 * in communicated atoms
7905 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7907 r = cg_cm[cg][dim0] - c->cr0;
7908 /* This is the first dimension, so always r >= 0 */
7915 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7917 r = cg_cm[cg][dim1] - c->cr1[zone];
7924 r = cg_cm[cg][dim1] - c->bcr1;
7934 /* Triclinic direction, more complicated */
7937 /* Rounding, conservative as the skew_fac multiplication
7938 * will slightly underestimate the distance.
7940 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7942 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7943 for (i = dim0+1; i < DIM; i++)
7945 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7947 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7950 rb[dim0] = rn[dim0];
7953 /* Take care that the cell planes along dim0 might not
7954 * be orthogonal to those along dim1 and dim2.
7956 for (i = 1; i <= dim_ind; i++)
7959 if (normal[dim0][dimd] > 0)
7961 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7964 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7969 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7971 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7973 for (i = dim1+1; i < DIM; i++)
7975 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7977 rn[dim1] += tric_sh;
7980 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7981 /* Take care of coupling of the distances
7982 * to the planes along dim0 and dim1 through dim2.
7984 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7985 /* Take care that the cell planes along dim1
7986 * might not be orthogonal to that along dim2.
7988 if (normal[dim1][dim2] > 0)
7990 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7996 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7999 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8000 /* Take care of coupling of the distances
8001 * to the planes along dim0 and dim1 through dim2.
8003 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8004 /* Take care that the cell planes along dim1
8005 * might not be orthogonal to that along dim2.
8007 if (normal[dim1][dim2] > 0)
8009 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8014 /* The distance along the communication direction */
8015 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8017 for (i = dim+1; i < DIM; i++)
8019 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8024 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8025 /* Take care of coupling of the distances
8026 * to the planes along dim0 and dim1 through dim2.
8028 if (dim_ind == 1 && zonei == 1)
8030 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8036 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8039 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8040 /* Take care of coupling of the distances
8041 * to the planes along dim0 and dim1 through dim2.
8043 if (dim_ind == 1 && zonei == 1)
8045 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8053 ((bDistMB && rb2 < r_bcomm2) ||
8054 (bDist2B && r2 < r_bcomm2)) &&
8056 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8057 missing_link(comm->cglink, index_gl[cg],
8060 /* Make an index to the local charge groups */
8061 if (nsend+1 > ind->nalloc)
8063 ind->nalloc = over_alloc_large(nsend+1);
8064 srenew(ind->index, ind->nalloc);
8066 if (nsend+1 > *ibuf_nalloc)
8068 *ibuf_nalloc = over_alloc_large(nsend+1);
8069 srenew(*ibuf, *ibuf_nalloc);
8071 ind->index[nsend] = cg;
8072 (*ibuf)[nsend] = index_gl[cg];
8074 vec_rvec_check_alloc(vbuf, nsend+1);
8076 if (dd->ci[dim] == 0)
8078 /* Correct cg_cm for pbc */
8079 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8082 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8083 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8088 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8091 nat += cgindex[cg+1] - cgindex[cg];
8097 *nsend_z_ptr = nsend_z;
8100 static void setup_dd_communication(gmx_domdec_t *dd,
8101 matrix box, gmx_ddbox_t *ddbox,
8102 t_forcerec *fr, t_state *state, rvec **f)
8104 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8105 int nzone, nzone_send, zone, zonei, cg0, cg1;
8106 int c, i, j, cg, cg_gl, nrcg;
8107 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8108 gmx_domdec_comm_t *comm;
8109 gmx_domdec_zones_t *zones;
8110 gmx_domdec_comm_dim_t *cd;
8111 gmx_domdec_ind_t *ind;
8112 cginfo_mb_t *cginfo_mb;
8113 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8114 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
8115 dd_corners_t corners;
8117 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8118 real skew_fac2_d, skew_fac_01;
8125 fprintf(debug, "Setting up DD communication\n");
8130 switch (fr->cutoff_scheme)
8139 gmx_incons("unimplemented");
8143 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8145 dim = dd->dim[dim_ind];
8147 /* Check if we need to use triclinic distances */
8148 tric_dist[dim_ind] = 0;
8149 for (i = 0; i <= dim_ind; i++)
8151 if (ddbox->tric_dir[dd->dim[i]])
8153 tric_dist[dim_ind] = 1;
8158 bBondComm = comm->bBondComm;
8160 /* Do we need to determine extra distances for multi-body bondeds? */
8161 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8163 /* Do we need to determine extra distances for only two-body bondeds? */
8164 bDist2B = (bBondComm && !bDistMB);
8166 r_comm2 = sqr(comm->cutoff);
8167 r_bcomm2 = sqr(comm->cutoff_mbody);
8171 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8174 zones = &comm->zones;
8177 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8178 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8180 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8182 /* Triclinic stuff */
8183 normal = ddbox->normal;
8187 v_0 = ddbox->v[dim0];
8188 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8190 /* Determine the coupling coefficient for the distances
8191 * to the cell planes along dim0 and dim1 through dim2.
8192 * This is required for correct rounding.
8195 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8198 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8204 v_1 = ddbox->v[dim1];
8207 zone_cg_range = zones->cg_range;
8208 index_gl = dd->index_gl;
8209 cgindex = dd->cgindex;
8210 cginfo_mb = fr->cginfo_mb;
8212 zone_cg_range[0] = 0;
8213 zone_cg_range[1] = dd->ncg_home;
8214 comm->zone_ncg1[0] = dd->ncg_home;
8215 pos_cg = dd->ncg_home;
8217 nat_tot = dd->nat_home;
8219 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8221 dim = dd->dim[dim_ind];
8222 cd = &comm->cd[dim_ind];
8224 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8226 /* No pbc in this dimension, the first node should not comm. */
8234 v_d = ddbox->v[dim];
8235 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8237 cd->bInPlace = TRUE;
8238 for (p = 0; p < cd->np; p++)
8240 /* Only atoms communicated in the first pulse are used
8241 * for multi-body bonded interactions or for bBondComm.
8243 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8248 for (zone = 0; zone < nzone_send; zone++)
8250 if (tric_dist[dim_ind] && dim_ind > 0)
8252 /* Determine slightly more optimized skew_fac's
8254 * This reduces the number of communicated atoms
8255 * by about 10% for 3D DD of rhombic dodecahedra.
8257 for (dimd = 0; dimd < dim; dimd++)
8259 sf2_round[dimd] = 1;
8260 if (ddbox->tric_dir[dimd])
8262 for (i = dd->dim[dimd]+1; i < DIM; i++)
8264 /* If we are shifted in dimension i
8265 * and the cell plane is tilted forward
8266 * in dimension i, skip this coupling.
8268 if (!(zones->shift[nzone+zone][i] &&
8269 ddbox->v[dimd][i][dimd] >= 0))
8272 sqr(ddbox->v[dimd][i][dimd]);
8275 sf2_round[dimd] = 1/sf2_round[dimd];
8280 zonei = zone_perm[dim_ind][zone];
8283 /* Here we permutate the zones to obtain a convenient order
8284 * for neighbor searching
8286 cg0 = zone_cg_range[zonei];
8287 cg1 = zone_cg_range[zonei+1];
8291 /* Look only at the cg's received in the previous grid pulse
8293 cg1 = zone_cg_range[nzone+zone+1];
8294 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8297 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8298 for (th = 0; th < comm->nth; th++)
8300 gmx_domdec_ind_t *ind_p;
8301 int **ibuf_p, *ibuf_nalloc_p;
8303 int *nsend_p, *nat_p;
8309 /* Thread 0 writes in the comm buffers */
8311 ibuf_p = &comm->buf_int;
8312 ibuf_nalloc_p = &comm->nalloc_int;
8313 vbuf_p = &comm->vbuf;
8316 nsend_zone_p = &ind->nsend[zone];
8320 /* Other threads write into temp buffers */
8321 ind_p = &comm->dth[th].ind;
8322 ibuf_p = &comm->dth[th].ibuf;
8323 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8324 vbuf_p = &comm->dth[th].vbuf;
8325 nsend_p = &comm->dth[th].nsend;
8326 nat_p = &comm->dth[th].nat;
8327 nsend_zone_p = &comm->dth[th].nsend_zone;
8329 comm->dth[th].nsend = 0;
8330 comm->dth[th].nat = 0;
8331 comm->dth[th].nsend_zone = 0;
8341 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8342 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8345 /* Get the cg's for this pulse in this zone */
8346 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8348 dim, dim_ind, dim0, dim1, dim2,
8351 normal, skew_fac2_d, skew_fac_01,
8352 v_d, v_0, v_1, &corners, sf2_round,
8353 bDistBonded, bBondComm,
8357 ibuf_p, ibuf_nalloc_p,
8363 /* Append data of threads>=1 to the communication buffers */
8364 for (th = 1; th < comm->nth; th++)
8366 dd_comm_setup_work_t *dth;
8369 dth = &comm->dth[th];
8371 ns1 = nsend + dth->nsend_zone;
8372 if (ns1 > ind->nalloc)
8374 ind->nalloc = over_alloc_dd(ns1);
8375 srenew(ind->index, ind->nalloc);
8377 if (ns1 > comm->nalloc_int)
8379 comm->nalloc_int = over_alloc_dd(ns1);
8380 srenew(comm->buf_int, comm->nalloc_int);
8382 if (ns1 > comm->vbuf.nalloc)
8384 comm->vbuf.nalloc = over_alloc_dd(ns1);
8385 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8388 for (i = 0; i < dth->nsend_zone; i++)
8390 ind->index[nsend] = dth->ind.index[i];
8391 comm->buf_int[nsend] = dth->ibuf[i];
8392 copy_rvec(dth->vbuf.v[i],
8393 comm->vbuf.v[nsend]);
8397 ind->nsend[zone] += dth->nsend_zone;
8400 /* Clear the counts in case we do not have pbc */
8401 for (zone = nzone_send; zone < nzone; zone++)
8403 ind->nsend[zone] = 0;
8405 ind->nsend[nzone] = nsend;
8406 ind->nsend[nzone+1] = nat;
8407 /* Communicate the number of cg's and atoms to receive */
8408 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8409 ind->nsend, nzone+2,
8410 ind->nrecv, nzone+2);
8412 /* The rvec buffer is also required for atom buffers of size nsend
8413 * in dd_move_x and dd_move_f.
8415 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8419 /* We can receive in place if only the last zone is not empty */
8420 for (zone = 0; zone < nzone-1; zone++)
8422 if (ind->nrecv[zone] > 0)
8424 cd->bInPlace = FALSE;
8429 /* The int buffer is only required here for the cg indices */
8430 if (ind->nrecv[nzone] > comm->nalloc_int2)
8432 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8433 srenew(comm->buf_int2, comm->nalloc_int2);
8435 /* The rvec buffer is also required for atom buffers
8436 * of size nrecv in dd_move_x and dd_move_f.
8438 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8439 vec_rvec_check_alloc(&comm->vbuf2, i);
8443 /* Make space for the global cg indices */
8444 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8445 || dd->cg_nalloc == 0)
8447 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8448 srenew(index_gl, dd->cg_nalloc);
8449 srenew(cgindex, dd->cg_nalloc+1);
8451 /* Communicate the global cg indices */
8454 recv_i = index_gl + pos_cg;
8458 recv_i = comm->buf_int2;
8460 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8461 comm->buf_int, nsend,
8462 recv_i, ind->nrecv[nzone]);
8464 /* Make space for cg_cm */
8465 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8466 if (fr->cutoff_scheme == ecutsGROUP)
8474 /* Communicate cg_cm */
8477 recv_vr = cg_cm + pos_cg;
8481 recv_vr = comm->vbuf2.v;
8483 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8484 comm->vbuf.v, nsend,
8485 recv_vr, ind->nrecv[nzone]);
8487 /* Make the charge group index */
8490 zone = (p == 0 ? 0 : nzone - 1);
8491 while (zone < nzone)
8493 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8495 cg_gl = index_gl[pos_cg];
8496 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8497 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8498 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8501 /* Update the charge group presence,
8502 * so we can use it in the next pass of the loop.
8504 comm->bLocalCG[cg_gl] = TRUE;
8510 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8513 zone_cg_range[nzone+zone] = pos_cg;
8518 /* This part of the code is never executed with bBondComm. */
8519 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8520 index_gl, recv_i, cg_cm, recv_vr,
8521 cgindex, fr->cginfo_mb, fr->cginfo);
8522 pos_cg += ind->nrecv[nzone];
8524 nat_tot += ind->nrecv[nzone+1];
8528 /* Store the atom block for easy copying of communication buffers */
8529 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8533 dd->index_gl = index_gl;
8534 dd->cgindex = cgindex;
8536 dd->ncg_tot = zone_cg_range[zones->n];
8537 dd->nat_tot = nat_tot;
8538 comm->nat[ddnatHOME] = dd->nat_home;
8539 for (i = ddnatZONE; i < ddnatNR; i++)
8541 comm->nat[i] = dd->nat_tot;
8546 /* We don't need to update cginfo, since that was alrady done above.
8547 * So we pass NULL for the forcerec.
8549 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8550 NULL, comm->bLocalCG);
8555 fprintf(debug, "Finished setting up DD communication, zones:");
8556 for (c = 0; c < zones->n; c++)
8558 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8560 fprintf(debug, "\n");
8564 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8568 for (c = 0; c < zones->nizone; c++)
8570 zones->izone[c].cg1 = zones->cg_range[c+1];
8571 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8572 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8576 static void set_zones_size(gmx_domdec_t *dd,
8577 matrix box, const gmx_ddbox_t *ddbox,
8578 int zone_start, int zone_end)
8580 gmx_domdec_comm_t *comm;
8581 gmx_domdec_zones_t *zones;
8583 int z, zi, zj0, zj1, d, dim;
8586 real size_j, add_tric;
8591 zones = &comm->zones;
8593 /* Do we need to determine extra distances for multi-body bondeds? */
8594 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8596 for (z = zone_start; z < zone_end; z++)
8598 /* Copy cell limits to zone limits.
8599 * Valid for non-DD dims and non-shifted dims.
8601 copy_rvec(comm->cell_x0, zones->size[z].x0);
8602 copy_rvec(comm->cell_x1, zones->size[z].x1);
8605 for (d = 0; d < dd->ndim; d++)
8609 for (z = 0; z < zones->n; z++)
8611 /* With a staggered grid we have different sizes
8612 * for non-shifted dimensions.
8614 if (dd->bGridJump && zones->shift[z][dim] == 0)
8618 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8619 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8623 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8624 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8630 rcmbs = comm->cutoff_mbody;
8631 if (ddbox->tric_dir[dim])
8633 rcs /= ddbox->skew_fac[dim];
8634 rcmbs /= ddbox->skew_fac[dim];
8637 /* Set the lower limit for the shifted zone dimensions */
8638 for (z = zone_start; z < zone_end; z++)
8640 if (zones->shift[z][dim] > 0)
8643 if (!dd->bGridJump || d == 0)
8645 zones->size[z].x0[dim] = comm->cell_x1[dim];
8646 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8650 /* Here we take the lower limit of the zone from
8651 * the lowest domain of the zone below.
8655 zones->size[z].x0[dim] =
8656 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8662 zones->size[z].x0[dim] =
8663 zones->size[zone_perm[2][z-4]].x0[dim];
8667 zones->size[z].x0[dim] =
8668 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8671 /* A temporary limit, is updated below */
8672 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8676 for (zi = 0; zi < zones->nizone; zi++)
8678 if (zones->shift[zi][dim] == 0)
8680 /* This takes the whole zone into account.
8681 * With multiple pulses this will lead
8682 * to a larger zone then strictly necessary.
8684 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8685 zones->size[zi].x1[dim]+rcmbs);
8693 /* Loop over the i-zones to set the upper limit of each
8696 for (zi = 0; zi < zones->nizone; zi++)
8698 if (zones->shift[zi][dim] == 0)
8700 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8702 if (zones->shift[z][dim] > 0)
8704 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8705 zones->size[zi].x1[dim]+rcs);
8712 for (z = zone_start; z < zone_end; z++)
8714 /* Initialization only required to keep the compiler happy */
8715 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8718 /* To determine the bounding box for a zone we need to find
8719 * the extreme corners of 4, 2 or 1 corners.
8721 nc = 1 << (ddbox->npbcdim - 1);
8723 for (c = 0; c < nc; c++)
8725 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8729 corner[YY] = zones->size[z].x0[YY];
8733 corner[YY] = zones->size[z].x1[YY];
8737 corner[ZZ] = zones->size[z].x0[ZZ];
8741 corner[ZZ] = zones->size[z].x1[ZZ];
8743 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8745 /* With 1D domain decomposition the cg's are not in
8746 * the triclinic box, but triclinic x-y and rectangular y-z.
8747 * Shift y back, so it will later end up at 0.
8749 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8751 /* Apply the triclinic couplings */
8752 assert(ddbox->npbcdim <= DIM);
8753 for (i = YY; i < ddbox->npbcdim; i++)
8755 for (j = XX; j < i; j++)
8757 corner[j] += corner[i]*box[i][j]/box[i][i];
8762 copy_rvec(corner, corner_min);
8763 copy_rvec(corner, corner_max);
8767 for (i = 0; i < DIM; i++)
8769 corner_min[i] = min(corner_min[i], corner[i]);
8770 corner_max[i] = max(corner_max[i], corner[i]);
8774 /* Copy the extreme cornes without offset along x */
8775 for (i = 0; i < DIM; i++)
8777 zones->size[z].bb_x0[i] = corner_min[i];
8778 zones->size[z].bb_x1[i] = corner_max[i];
8780 /* Add the offset along x */
8781 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8782 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8785 if (zone_start == 0)
8788 for (dim = 0; dim < DIM; dim++)
8790 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8792 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8797 for (z = zone_start; z < zone_end; z++)
8799 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8801 zones->size[z].x0[XX], zones->size[z].x1[XX],
8802 zones->size[z].x0[YY], zones->size[z].x1[YY],
8803 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8804 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8806 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8807 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8808 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8813 static int comp_cgsort(const void *a, const void *b)
8817 gmx_cgsort_t *cga, *cgb;
8818 cga = (gmx_cgsort_t *)a;
8819 cgb = (gmx_cgsort_t *)b;
8821 comp = cga->nsc - cgb->nsc;
8824 comp = cga->ind_gl - cgb->ind_gl;
8830 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8835 /* Order the data */
8836 for (i = 0; i < n; i++)
8838 buf[i] = a[sort[i].ind];
8841 /* Copy back to the original array */
8842 for (i = 0; i < n; i++)
8848 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8853 /* Order the data */
8854 for (i = 0; i < n; i++)
8856 copy_rvec(v[sort[i].ind], buf[i]);
8859 /* Copy back to the original array */
8860 for (i = 0; i < n; i++)
8862 copy_rvec(buf[i], v[i]);
8866 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8869 int a, atot, cg, cg0, cg1, i;
8871 if (cgindex == NULL)
8873 /* Avoid the useless loop of the atoms within a cg */
8874 order_vec_cg(ncg, sort, v, buf);
8879 /* Order the data */
8881 for (cg = 0; cg < ncg; cg++)
8883 cg0 = cgindex[sort[cg].ind];
8884 cg1 = cgindex[sort[cg].ind+1];
8885 for (i = cg0; i < cg1; i++)
8887 copy_rvec(v[i], buf[a]);
8893 /* Copy back to the original array */
8894 for (a = 0; a < atot; a++)
8896 copy_rvec(buf[a], v[a]);
8900 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8901 int nsort_new, gmx_cgsort_t *sort_new,
8902 gmx_cgsort_t *sort1)
8906 /* The new indices are not very ordered, so we qsort them */
8907 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8909 /* sort2 is already ordered, so now we can merge the two arrays */
8913 while (i2 < nsort2 || i_new < nsort_new)
8917 sort1[i1++] = sort_new[i_new++];
8919 else if (i_new == nsort_new)
8921 sort1[i1++] = sort2[i2++];
8923 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8924 (sort2[i2].nsc == sort_new[i_new].nsc &&
8925 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8927 sort1[i1++] = sort2[i2++];
8931 sort1[i1++] = sort_new[i_new++];
8936 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8938 gmx_domdec_sort_t *sort;
8939 gmx_cgsort_t *cgsort, *sort_i;
8940 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8941 int sort_last, sort_skip;
8943 sort = dd->comm->sort;
8945 a = fr->ns.grid->cell_index;
8947 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8949 if (ncg_home_old >= 0)
8951 /* The charge groups that remained in the same ns grid cell
8952 * are completely ordered. So we can sort efficiently by sorting
8953 * the charge groups that did move into the stationary list.
8958 for (i = 0; i < dd->ncg_home; i++)
8960 /* Check if this cg did not move to another node */
8963 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8965 /* This cg is new on this node or moved ns grid cell */
8966 if (nsort_new >= sort->sort_new_nalloc)
8968 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8969 srenew(sort->sort_new, sort->sort_new_nalloc);
8971 sort_i = &(sort->sort_new[nsort_new++]);
8975 /* This cg did not move */
8976 sort_i = &(sort->sort2[nsort2++]);
8978 /* Sort on the ns grid cell indices
8979 * and the global topology index.
8980 * index_gl is irrelevant with cell ns,
8981 * but we set it here anyhow to avoid a conditional.
8984 sort_i->ind_gl = dd->index_gl[i];
8991 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8994 /* Sort efficiently */
8995 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9000 cgsort = sort->sort;
9002 for (i = 0; i < dd->ncg_home; i++)
9004 /* Sort on the ns grid cell indices
9005 * and the global topology index
9007 cgsort[i].nsc = a[i];
9008 cgsort[i].ind_gl = dd->index_gl[i];
9010 if (cgsort[i].nsc < moved)
9017 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9019 /* Determine the order of the charge groups using qsort */
9020 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9026 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9029 int ncg_new, i, *a, na;
9031 sort = dd->comm->sort->sort;
9033 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9036 for (i = 0; i < na; i++)
9040 sort[ncg_new].ind = a[i];
9048 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9051 gmx_domdec_sort_t *sort;
9052 gmx_cgsort_t *cgsort, *sort_i;
9054 int ncg_new, i, *ibuf, cgsize;
9057 sort = dd->comm->sort;
9059 if (dd->ncg_home > sort->sort_nalloc)
9061 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9062 srenew(sort->sort, sort->sort_nalloc);
9063 srenew(sort->sort2, sort->sort_nalloc);
9065 cgsort = sort->sort;
9067 switch (fr->cutoff_scheme)
9070 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9073 ncg_new = dd_sort_order_nbnxn(dd, fr);
9076 gmx_incons("unimplemented");
9080 /* We alloc with the old size, since cgindex is still old */
9081 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9082 vbuf = dd->comm->vbuf.v;
9086 cgindex = dd->cgindex;
9093 /* Remove the charge groups which are no longer at home here */
9094 dd->ncg_home = ncg_new;
9097 fprintf(debug, "Set the new home charge group count to %d\n",
9101 /* Reorder the state */
9102 for (i = 0; i < estNR; i++)
9104 if (EST_DISTR(i) && (state->flags & (1<<i)))
9109 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9112 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9115 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9118 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9122 case estDISRE_INITF:
9123 case estDISRE_RM3TAV:
9124 case estORIRE_INITF:
9126 /* No ordering required */
9129 gmx_incons("Unknown state entry encountered in dd_sort_state");
9134 if (fr->cutoff_scheme == ecutsGROUP)
9137 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9140 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9142 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9143 srenew(sort->ibuf, sort->ibuf_nalloc);
9146 /* Reorder the global cg index */
9147 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9148 /* Reorder the cginfo */
9149 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9150 /* Rebuild the local cg index */
9154 for (i = 0; i < dd->ncg_home; i++)
9156 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9157 ibuf[i+1] = ibuf[i] + cgsize;
9159 for (i = 0; i < dd->ncg_home+1; i++)
9161 dd->cgindex[i] = ibuf[i];
9166 for (i = 0; i < dd->ncg_home+1; i++)
9171 /* Set the home atom number */
9172 dd->nat_home = dd->cgindex[dd->ncg_home];
9174 if (fr->cutoff_scheme == ecutsVERLET)
9176 /* The atoms are now exactly in grid order, update the grid order */
9177 nbnxn_set_atomorder(fr->nbv->nbs);
9181 /* Copy the sorted ns cell indices back to the ns grid struct */
9182 for (i = 0; i < dd->ncg_home; i++)
9184 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9186 fr->ns.grid->nr = dd->ncg_home;
9190 static void add_dd_statistics(gmx_domdec_t *dd)
9192 gmx_domdec_comm_t *comm;
9197 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9199 comm->sum_nat[ddnat-ddnatZONE] +=
9200 comm->nat[ddnat] - comm->nat[ddnat-1];
9205 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9207 gmx_domdec_comm_t *comm;
9212 /* Reset all the statistics and counters for total run counting */
9213 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9215 comm->sum_nat[ddnat-ddnatZONE] = 0;
9219 comm->load_step = 0;
9222 clear_ivec(comm->load_lim);
9227 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9229 gmx_domdec_comm_t *comm;
9233 comm = cr->dd->comm;
9235 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9242 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");
9244 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9246 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9251 " av. #atoms communicated per step for force: %d x %.1f\n",
9255 if (cr->dd->vsite_comm)
9258 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9259 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9264 if (cr->dd->constraint_comm)
9267 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9268 1 + ir->nLincsIter, av);
9272 gmx_incons(" Unknown type for DD statistics");
9275 fprintf(fplog, "\n");
9277 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9279 print_dd_load_av(fplog, cr->dd);
9283 void dd_partition_system(FILE *fplog,
9286 gmx_bool bMasterState,
9288 t_state *state_global,
9289 gmx_mtop_t *top_global,
9291 t_state *state_local,
9294 gmx_localtop_t *top_local,
9297 gmx_shellfc_t shellfc,
9298 gmx_constr_t constr,
9300 gmx_wallcycle_t wcycle,
9304 gmx_domdec_comm_t *comm;
9305 gmx_ddbox_t ddbox = {0};
9307 gmx_int64_t step_pcoupl;
9308 rvec cell_ns_x0, cell_ns_x1;
9309 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9310 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9311 gmx_bool bRedist, bSortCG, bResortAll;
9312 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9319 bBoxChanged = (bMasterState || DEFORM(*ir));
9320 if (ir->epc != epcNO)
9322 /* With nstpcouple > 1 pressure coupling happens.
9323 * one step after calculating the pressure.
9324 * Box scaling happens at the end of the MD step,
9325 * after the DD partitioning.
9326 * We therefore have to do DLB in the first partitioning
9327 * after an MD step where P-coupling occured.
9328 * We need to determine the last step in which p-coupling occurred.
9329 * MRS -- need to validate this for vv?
9334 step_pcoupl = step - 1;
9338 step_pcoupl = ((step - 1)/n)*n + 1;
9340 if (step_pcoupl >= comm->partition_step)
9346 bNStGlobalComm = (step % nstglobalcomm == 0);
9348 if (!comm->bDynLoadBal)
9354 /* Should we do dynamic load balacing this step?
9355 * Since it requires (possibly expensive) global communication,
9356 * we might want to do DLB less frequently.
9358 if (bBoxChanged || ir->epc != epcNO)
9360 bDoDLB = bBoxChanged;
9364 bDoDLB = bNStGlobalComm;
9368 /* Check if we have recorded loads on the nodes */
9369 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9371 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal && !dd_dlb_is_locked(dd))
9373 /* Check if we should use DLB at the second partitioning
9374 * and every 100 partitionings,
9375 * so the extra communication cost is negligible.
9377 const int nddp_chk_dlb = 100;
9379 bCheckDLB = (comm->n_load_collect == 0 ||
9380 comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
9387 /* Print load every nstlog, first and last step to the log file */
9388 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9389 comm->n_load_collect == 0 ||
9391 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9393 /* Avoid extra communication due to verbose screen output
9394 * when nstglobalcomm is set.
9396 if (bDoDLB || bLogLoad || bCheckDLB ||
9397 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9399 get_load_distribution(dd, wcycle);
9404 dd_print_load(fplog, dd, step-1);
9408 dd_print_load_verbose(dd);
9411 comm->n_load_collect++;
9415 /* Since the timings are node dependent, the master decides */
9418 /* Here we check if the max PME rank load is more than 0.98
9419 * the max PP force load. If so, PP DLB will not help,
9420 * since we are (almost) limited by PME. Furthermore,
9421 * DLB will cause a significant extra x/f redistribution
9422 * cost on the PME ranks, which will then surely result
9423 * in lower total performance.
9424 * This check might be fragile, since one measurement
9425 * below 0.98 (although only done once every 100 DD part.)
9426 * could turn on DLB for the rest of the run.
9428 if (cr->npmenodes > 0 &&
9429 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9436 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9440 fprintf(debug, "step %s, imb loss %f\n",
9441 gmx_step_str(step, sbuf),
9442 dd_force_imb_perf_loss(dd));
9445 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9448 turn_on_dlb(fplog, cr, step);
9453 comm->n_load_have++;
9456 cgs_gl = &comm->cgs_gl;
9461 /* Clear the old state */
9462 clear_dd_indices(dd, 0, 0);
9465 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9466 TRUE, cgs_gl, state_global->x, &ddbox);
9468 get_cg_distribution(fplog, step, dd, cgs_gl,
9469 state_global->box, &ddbox, state_global->x);
9471 dd_distribute_state(dd, cgs_gl,
9472 state_global, state_local, f);
9474 dd_make_local_cgs(dd, &top_local->cgs);
9476 /* Ensure that we have space for the new distribution */
9477 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9479 if (fr->cutoff_scheme == ecutsGROUP)
9481 calc_cgcm(fplog, 0, dd->ncg_home,
9482 &top_local->cgs, state_local->x, fr->cg_cm);
9485 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9487 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9489 else if (state_local->ddp_count != dd->ddp_count)
9491 if (state_local->ddp_count > dd->ddp_count)
9493 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9496 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9498 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);
9501 /* Clear the old state */
9502 clear_dd_indices(dd, 0, 0);
9504 /* Build the new indices */
9505 rebuild_cgindex(dd, cgs_gl->index, state_local);
9506 make_dd_indices(dd, cgs_gl->index, 0);
9507 ncgindex_set = dd->ncg_home;
9509 if (fr->cutoff_scheme == ecutsGROUP)
9511 /* Redetermine the cg COMs */
9512 calc_cgcm(fplog, 0, dd->ncg_home,
9513 &top_local->cgs, state_local->x, fr->cg_cm);
9516 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9518 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9520 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9521 TRUE, &top_local->cgs, state_local->x, &ddbox);
9523 bRedist = comm->bDynLoadBal;
9527 /* We have the full state, only redistribute the cgs */
9529 /* Clear the non-home indices */
9530 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9533 /* Avoid global communication for dim's without pbc and -gcom */
9534 if (!bNStGlobalComm)
9536 copy_rvec(comm->box0, ddbox.box0 );
9537 copy_rvec(comm->box_size, ddbox.box_size);
9539 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9540 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9545 /* For dim's without pbc and -gcom */
9546 copy_rvec(ddbox.box0, comm->box0 );
9547 copy_rvec(ddbox.box_size, comm->box_size);
9549 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9552 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9554 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9557 /* Check if we should sort the charge groups */
9558 if (comm->nstSortCG > 0)
9560 bSortCG = (bMasterState ||
9561 (bRedist && (step % comm->nstSortCG == 0)));
9568 ncg_home_old = dd->ncg_home;
9573 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9575 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9577 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9579 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9582 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9584 &comm->cell_x0, &comm->cell_x1,
9585 dd->ncg_home, fr->cg_cm,
9586 cell_ns_x0, cell_ns_x1, &grid_density);
9590 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9593 switch (fr->cutoff_scheme)
9596 copy_ivec(fr->ns.grid->n, ncells_old);
9597 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9598 state_local->box, cell_ns_x0, cell_ns_x1,
9599 fr->rlistlong, grid_density);
9602 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9605 gmx_incons("unimplemented");
9607 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9608 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9612 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9614 /* Sort the state on charge group position.
9615 * This enables exact restarts from this step.
9616 * It also improves performance by about 15% with larger numbers
9617 * of atoms per node.
9620 /* Fill the ns grid with the home cell,
9621 * so we can sort with the indices.
9623 set_zones_ncg_home(dd);
9625 switch (fr->cutoff_scheme)
9628 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9630 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9632 comm->zones.size[0].bb_x0,
9633 comm->zones.size[0].bb_x1,
9635 comm->zones.dens_zone0,
9638 ncg_moved, bRedist ? comm->moved : NULL,
9639 fr->nbv->grp[eintLocal].kernel_type,
9640 fr->nbv->grp[eintLocal].nbat);
9642 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9645 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9646 0, dd->ncg_home, fr->cg_cm);
9648 copy_ivec(fr->ns.grid->n, ncells_new);
9651 gmx_incons("unimplemented");
9654 bResortAll = bMasterState;
9656 /* Check if we can user the old order and ns grid cell indices
9657 * of the charge groups to sort the charge groups efficiently.
9659 if (ncells_new[XX] != ncells_old[XX] ||
9660 ncells_new[YY] != ncells_old[YY] ||
9661 ncells_new[ZZ] != ncells_old[ZZ])
9668 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9669 gmx_step_str(step, sbuf), dd->ncg_home);
9671 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9672 bResortAll ? -1 : ncg_home_old);
9673 /* Rebuild all the indices */
9674 ga2la_clear(dd->ga2la);
9677 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9680 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9682 /* Setup up the communication and communicate the coordinates */
9683 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9685 /* Set the indices */
9686 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9688 /* Set the charge group boundaries for neighbor searching */
9689 set_cg_boundaries(&comm->zones);
9691 if (fr->cutoff_scheme == ecutsVERLET)
9693 set_zones_size(dd, state_local->box, &ddbox,
9694 bSortCG ? 1 : 0, comm->zones.n);
9697 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9700 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9701 -1,state_local->x,state_local->box);
9704 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9706 /* Extract a local topology from the global topology */
9707 for (i = 0; i < dd->ndim; i++)
9709 np[dd->dim[i]] = comm->cd[i].np;
9711 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9712 comm->cellsize_min, np,
9714 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9715 vsite, top_global, top_local);
9717 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9719 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9721 /* Set up the special atom communication */
9722 n = comm->nat[ddnatZONE];
9723 for (i = ddnatZONE+1; i < ddnatNR; i++)
9728 if (vsite && vsite->n_intercg_vsite)
9730 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9734 if (dd->bInterCGcons || dd->bInterCGsettles)
9736 /* Only for inter-cg constraints we need special code */
9737 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9738 constr, ir->nProjOrder,
9739 top_local->idef.il);
9743 gmx_incons("Unknown special atom type setup");
9748 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9750 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9752 /* Make space for the extra coordinates for virtual site
9753 * or constraint communication.
9755 state_local->natoms = comm->nat[ddnatNR-1];
9756 if (state_local->natoms > state_local->nalloc)
9758 dd_realloc_state(state_local, f, state_local->natoms);
9761 if (fr->bF_NoVirSum)
9763 if (vsite && vsite->n_intercg_vsite)
9765 nat_f_novirsum = comm->nat[ddnatVSITE];
9769 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9771 nat_f_novirsum = dd->nat_tot;
9775 nat_f_novirsum = dd->nat_home;
9784 /* Set the number of atoms required for the force calculation.
9785 * Forces need to be constrained when using a twin-range setup
9786 * or with energy minimization. For simple simulations we could
9787 * avoid some allocation, zeroing and copying, but this is
9788 * probably not worth the complications ande checking.
9790 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9791 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9793 /* We make the all mdatoms up to nat_tot_con.
9794 * We could save some work by only setting invmass
9795 * between nat_tot and nat_tot_con.
9797 /* This call also sets the new number of home particles to dd->nat_home */
9798 atoms2md(top_global, ir,
9799 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9801 /* Now we have the charges we can sort the FE interactions */
9802 dd_sort_local_top(dd, mdatoms, top_local);
9806 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9807 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9808 mdatoms, FALSE, vsite);
9813 /* Make the local shell stuff, currently no communication is done */
9814 make_local_shells(cr, mdatoms, shellfc);
9817 if (ir->implicit_solvent)
9819 make_local_gb(cr, fr->born, ir->gb_algorithm);
9822 setup_bonded_threading(fr, &top_local->idef);
9824 if (!(cr->duty & DUTY_PME))
9826 /* Send the charges and/or c6/sigmas to our PME only node */
9827 gmx_pme_send_parameters(cr,
9829 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9830 mdatoms->chargeA, mdatoms->chargeB,
9831 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9832 mdatoms->sigmaA, mdatoms->sigmaB,
9833 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9838 set_constraints(constr, top_local, ir, mdatoms, cr);
9841 if (ir->ePull != epullNO)
9843 /* Update the local pull groups */
9844 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9849 /* Update the local rotation groups */
9850 dd_make_local_rotation_groups(dd, ir->rot);
9853 if (ir->eSwapCoords != eswapNO)
9855 /* Update the local groups needed for ion swapping */
9856 dd_make_local_swap_groups(dd, ir->swap);
9859 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9860 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9862 add_dd_statistics(dd);
9864 /* Make sure we only count the cycles for this DD partitioning */
9865 clear_dd_cycle_counts(dd);
9867 /* Because the order of the atoms might have changed since
9868 * the last vsite construction, we need to communicate the constructing
9869 * atom coordinates again (for spreading the forces this MD step).
9871 dd_move_x_vsites(dd, state_local->box, state_local->x);
9873 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9875 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9877 dd_move_x(dd, state_local->box, state_local->x);
9878 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9879 -1, state_local->x, state_local->box);
9882 /* Store the partitioning step */
9883 comm->partition_step = step;
9885 /* Increase the DD partitioning counter */
9887 /* The state currently matches this DD partitioning count, store it */
9888 state_local->ddp_count = dd->ddp_count;
9891 /* The DD master node knows the complete cg distribution,
9892 * store the count so we can possibly skip the cg info communication.
9894 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9897 if (comm->DD_debug > 0)
9899 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9900 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9901 "after partitioning");