2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015, 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.
51 #include "gromacs/domdec/domdec_network.h"
52 #include "gromacs/ewald/pme.h"
53 #include "gromacs/fileio/gmxfio.h"
54 #include "gromacs/fileio/pdbio.h"
55 #include "gromacs/gmxlib/gpu_utils/gpu_utils.h"
56 #include "gromacs/imd/imd.h"
57 #include "gromacs/legacyheaders/chargegroup.h"
58 #include "gromacs/legacyheaders/constr.h"
59 #include "gromacs/legacyheaders/force.h"
60 #include "gromacs/legacyheaders/genborn.h"
61 #include "gromacs/legacyheaders/gmx_ga2la.h"
62 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
63 #include "gromacs/legacyheaders/mdatoms.h"
64 #include "gromacs/legacyheaders/mdrun.h"
65 #include "gromacs/legacyheaders/names.h"
66 #include "gromacs/legacyheaders/network.h"
67 #include "gromacs/legacyheaders/nrnb.h"
68 #include "gromacs/legacyheaders/nsgrid.h"
69 #include "gromacs/legacyheaders/shellfc.h"
70 #include "gromacs/legacyheaders/typedefs.h"
71 #include "gromacs/legacyheaders/vsite.h"
72 #include "gromacs/legacyheaders/types/commrec.h"
73 #include "gromacs/legacyheaders/types/constr.h"
74 #include "gromacs/legacyheaders/types/enums.h"
75 #include "gromacs/legacyheaders/types/forcerec.h"
76 #include "gromacs/legacyheaders/types/hw_info.h"
77 #include "gromacs/legacyheaders/types/ifunc.h"
78 #include "gromacs/legacyheaders/types/inputrec.h"
79 #include "gromacs/legacyheaders/types/mdatom.h"
80 #include "gromacs/legacyheaders/types/nrnb.h"
81 #include "gromacs/legacyheaders/types/ns.h"
82 #include "gromacs/legacyheaders/types/nsgrid.h"
83 #include "gromacs/legacyheaders/types/shellfc.h"
84 #include "gromacs/legacyheaders/types/simple.h"
85 #include "gromacs/legacyheaders/types/state.h"
86 #include "gromacs/listed-forces/manage-threading.h"
87 #include "gromacs/math/vec.h"
88 #include "gromacs/math/vectypes.h"
89 #include "gromacs/mdlib/nb_verlet.h"
90 #include "gromacs/mdlib/nbnxn_search.h"
91 #include "gromacs/pbcutil/ishift.h"
92 #include "gromacs/pbcutil/pbc.h"
93 #include "gromacs/pulling/pull.h"
94 #include "gromacs/pulling/pull_rotation.h"
95 #include "gromacs/swap/swapcoords.h"
96 #include "gromacs/timing/wallcycle.h"
97 #include "gromacs/topology/block.h"
98 #include "gromacs/topology/idef.h"
99 #include "gromacs/topology/mtop_util.h"
100 #include "gromacs/topology/topology.h"
101 #include "gromacs/utility/basedefinitions.h"
102 #include "gromacs/utility/basenetwork.h"
103 #include "gromacs/utility/cstringutil.h"
104 #include "gromacs/utility/fatalerror.h"
105 #include "gromacs/utility/gmxmpi.h"
106 #include "gromacs/utility/qsort_threadsafe.h"
107 #include "gromacs/utility/real.h"
108 #include "gromacs/utility/smalloc.h"
110 #include "domdec_constraints.h"
111 #include "domdec_internal.h"
112 #include "domdec_vsite.h"
114 #define DDRANK(dd, rank) (rank)
115 #define DDMASTERRANK(dd) (dd->masterrank)
117 typedef struct gmx_domdec_master
119 /* The cell boundaries */
121 /* The global charge group division */
122 int *ncg; /* Number of home charge groups for each node */
123 int *index; /* Index of nnodes+1 into cg */
124 int *cg; /* Global charge group index */
125 int *nat; /* Number of home atoms for each node. */
126 int *ibuf; /* Buffer for communication */
127 rvec *vbuf; /* Buffer for state scattering and gathering */
128 } gmx_domdec_master_t;
132 /* The numbers of charge groups to send and receive for each cell
133 * that requires communication, the last entry contains the total
134 * number of atoms that needs to be communicated.
136 int nsend[DD_MAXIZONE+2];
137 int nrecv[DD_MAXIZONE+2];
138 /* The charge groups to send */
141 /* The atom range for non-in-place communication */
142 int cell2at0[DD_MAXIZONE];
143 int cell2at1[DD_MAXIZONE];
148 int np; /* Number of grid pulses in this dimension */
149 int np_dlb; /* For dlb, for use with edlbAUTO */
150 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
152 gmx_bool bInPlace; /* Can we communicate in place? */
153 } gmx_domdec_comm_dim_t;
157 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
158 real *cell_f; /* State var.: cell boundaries, box relative */
159 real *old_cell_f; /* Temp. var.: old cell size */
160 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
161 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
162 real *bound_min; /* Temp. var.: lower limit for cell boundary */
163 real *bound_max; /* Temp. var.: upper limit for cell boundary */
164 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
165 real *buf_ncd; /* Temp. var. */
168 #define DD_NLOAD_MAX 9
170 /* Here floats are accurate enough, since these variables
171 * only influence the load balancing, not the actual MD results.
198 gmx_cgsort_t *sort_new;
210 /* This enum determines the order of the coordinates.
211 * ddnatHOME and ddnatZONE should be first and second,
212 * the others can be ordered as wanted.
215 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
219 edlbAUTO, edlbNO, edlbYES, edlbNR
221 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
225 int dim; /* The dimension */
226 gmx_bool dim_match; /* Tells if DD and PME dims match */
227 int nslab; /* The number of PME slabs in this dimension */
228 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
229 int *pp_min; /* The minimum pp node location, size nslab */
230 int *pp_max; /* The maximum pp node location,size nslab */
231 int maxshift; /* The maximum shift for coordinate redistribution in PME */
236 real min0; /* The minimum bottom of this zone */
237 real max1; /* The maximum top of this zone */
238 real min1; /* The minimum top of this zone */
239 real mch0; /* The maximum bottom communicaton height for this zone */
240 real mch1; /* The maximum top communicaton height for this zone */
241 real p1_0; /* The bottom value of the first cell in this zone */
242 real p1_1; /* The top value of the first cell in this zone */
247 gmx_domdec_ind_t ind;
254 } dd_comm_setup_work_t;
256 typedef struct gmx_domdec_comm
258 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
259 * unless stated otherwise.
262 /* The number of decomposition dimensions for PME, 0: no PME */
264 /* The number of nodes doing PME (PP/PME or only PME) */
268 /* The communication setup including the PME only nodes */
269 gmx_bool bCartesianPP_PME;
272 int *pmenodes; /* size npmenodes */
273 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
274 * but with bCartesianPP_PME */
275 gmx_ddpme_t ddpme[2];
277 /* The DD particle-particle nodes only */
278 gmx_bool bCartesianPP;
279 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
281 /* The global charge groups */
284 /* Should we sort the cgs */
286 gmx_domdec_sort_t *sort;
288 /* Are there charge groups? */
291 /* Are there bonded and multi-body interactions between charge groups? */
292 gmx_bool bInterCGBondeds;
293 gmx_bool bInterCGMultiBody;
295 /* Data for the optional bonded interaction atom communication range */
302 /* Is eDLB=edlbAUTO locked such that we currently can't turn it on? */
303 gmx_bool bDLB_locked;
304 /* Are we actually using DLB? */
305 gmx_bool bDynLoadBal;
307 /* Cell sizes for static load balancing, first index cartesian */
310 /* The width of the communicated boundaries */
313 /* The minimum cell size (including triclinic correction) */
315 /* For dlb, for use with edlbAUTO */
316 rvec cellsize_min_dlb;
317 /* The lower limit for the DD cell size with DLB */
319 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
320 gmx_bool bVacDLBNoLimit;
322 /* With PME load balancing we set limits on DLB */
323 gmx_bool bPMELoadBalDLBLimits;
324 /* DLB needs to take into account that we want to allow this maximum
325 * cut-off (for PME load balancing), this could limit cell boundaries.
327 real PMELoadBal_max_cutoff;
329 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
331 /* box0 and box_size are required with dim's without pbc and -gcom */
335 /* The cell boundaries */
339 /* The old location of the cell boundaries, to check cg displacements */
343 /* The communication setup and charge group boundaries for the zones */
344 gmx_domdec_zones_t zones;
346 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
347 * cell boundaries of neighboring cells for dynamic load balancing.
349 gmx_ddzone_t zone_d1[2];
350 gmx_ddzone_t zone_d2[2][2];
352 /* The coordinate/force communication setup and indices */
353 gmx_domdec_comm_dim_t cd[DIM];
354 /* The maximum number of cells to communicate with in one dimension */
357 /* Which cg distribution is stored on the master node */
358 int master_cg_ddp_count;
360 /* The number of cg's received from the direct neighbors */
361 int zone_ncg1[DD_MAXZONE];
363 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
366 /* Array for signalling if atoms have moved to another domain */
370 /* Communication buffer for general use */
374 /* Communication buffer for general use */
377 /* Temporary storage for thread parallel communication setup */
379 dd_comm_setup_work_t *dth;
381 /* Communication buffers only used with multiple grid pulses */
386 /* Communication buffers for local redistribution */
388 int cggl_flag_nalloc[DIM*2];
390 int cgcm_state_nalloc[DIM*2];
392 /* Cell sizes for dynamic load balancing */
393 gmx_domdec_root_t **root;
397 real cell_f_max0[DIM];
398 real cell_f_min1[DIM];
400 /* Stuff for load communication */
401 gmx_bool bRecordLoad;
402 gmx_domdec_load_t *load;
403 int nrank_gpu_shared;
405 MPI_Comm *mpi_comm_load;
406 MPI_Comm mpi_comm_gpu_shared;
409 /* Maximum DLB scaling per load balancing step in percent */
413 float cycl[ddCyclNr];
414 int cycl_n[ddCyclNr];
415 float cycl_max[ddCyclNr];
416 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
420 /* How many times have did we have load measurements */
422 /* How many times have we collected the load measurements */
426 double sum_nat[ddnatNR-ddnatZONE];
436 /* The last partition step */
437 gmx_int64_t partition_step;
445 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
448 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
449 #define DD_FLAG_NRCG 65535
450 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
451 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
453 /* Zone permutation required to obtain consecutive charge groups
454 * for neighbor searching.
456 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
458 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
459 * components see only j zones with that component 0.
462 /* The DD zone order */
463 static const ivec dd_zo[DD_MAXZONE] =
464 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
469 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
474 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
479 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
481 /* Factors used to avoid problems due to rounding issues */
482 #define DD_CELL_MARGIN 1.0001
483 #define DD_CELL_MARGIN2 1.00005
484 /* Factor to account for pressure scaling during nstlist steps */
485 #define DD_PRES_SCALE_MARGIN 1.02
487 /* Turn on DLB when the load imbalance causes this amount of total loss.
488 * There is a bit of overhead with DLB and it's difficult to achieve
489 * a load imbalance of less than 2% with DLB.
491 #define DD_PERF_LOSS_DLB_ON 0.02
493 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
494 #define DD_PERF_LOSS_WARN 0.05
496 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
498 /* Use separate MPI send and receive commands
499 * when nnodes <= GMX_DD_NNODES_SENDRECV.
500 * This saves memory (and some copying for small nnodes).
501 * For high parallelization scatter and gather calls are used.
503 #define GMX_DD_NNODES_SENDRECV 4
507 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
509 static void index2xyz(ivec nc,int ind,ivec xyz)
511 xyz[XX] = ind % nc[XX];
512 xyz[YY] = (ind / nc[XX]) % nc[YY];
513 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
517 /* This order is required to minimize the coordinate communication in PME
518 * which uses decomposition in the x direction.
520 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
522 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
524 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
525 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
526 xyz[ZZ] = ind % nc[ZZ];
529 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
534 ddindex = dd_index(dd->nc, c);
535 if (dd->comm->bCartesianPP_PME)
537 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
539 else if (dd->comm->bCartesianPP)
542 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
553 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
555 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
558 int ddglatnr(gmx_domdec_t *dd, int i)
568 if (i >= dd->comm->nat[ddnatNR-1])
570 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
572 atnr = dd->gatindex[i] + 1;
578 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
580 return &dd->comm->cgs_gl;
583 static void vec_rvec_init(vec_rvec_t *v)
589 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
593 v->nalloc = over_alloc_dd(n);
594 srenew(v->v, v->nalloc);
598 void dd_store_state(gmx_domdec_t *dd, t_state *state)
602 if (state->ddp_count != dd->ddp_count)
604 gmx_incons("The state does not the domain decomposition state");
607 state->ncg_gl = dd->ncg_home;
608 if (state->ncg_gl > state->cg_gl_nalloc)
610 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
611 srenew(state->cg_gl, state->cg_gl_nalloc);
613 for (i = 0; i < state->ncg_gl; i++)
615 state->cg_gl[i] = dd->index_gl[i];
618 state->ddp_count_cg_gl = dd->ddp_count;
621 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
623 return &dd->comm->zones;
626 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
627 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
629 gmx_domdec_zones_t *zones;
632 zones = &dd->comm->zones;
635 while (icg >= zones->izone[izone].cg1)
644 else if (izone < zones->nizone)
646 *jcg0 = zones->izone[izone].jcg0;
650 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
651 icg, izone, zones->nizone);
654 *jcg1 = zones->izone[izone].jcg1;
656 for (d = 0; d < dd->ndim; d++)
659 shift0[dim] = zones->izone[izone].shift0[dim];
660 shift1[dim] = zones->izone[izone].shift1[dim];
661 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
663 /* A conservative approach, this can be optimized */
670 int dd_natoms_vsite(gmx_domdec_t *dd)
672 return dd->comm->nat[ddnatVSITE];
675 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
677 *at_start = dd->comm->nat[ddnatCON-1];
678 *at_end = dd->comm->nat[ddnatCON];
681 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
683 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
684 int *index, *cgindex;
685 gmx_domdec_comm_t *comm;
686 gmx_domdec_comm_dim_t *cd;
687 gmx_domdec_ind_t *ind;
688 rvec shift = {0, 0, 0}, *buf, *rbuf;
689 gmx_bool bPBC, bScrew;
693 cgindex = dd->cgindex;
698 nat_tot = dd->nat_home;
699 for (d = 0; d < dd->ndim; d++)
701 bPBC = (dd->ci[dd->dim[d]] == 0);
702 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
705 copy_rvec(box[dd->dim[d]], shift);
708 for (p = 0; p < cd->np; p++)
715 for (i = 0; i < ind->nsend[nzone]; i++)
717 at0 = cgindex[index[i]];
718 at1 = cgindex[index[i]+1];
719 for (j = at0; j < at1; j++)
721 copy_rvec(x[j], buf[n]);
728 for (i = 0; i < ind->nsend[nzone]; i++)
730 at0 = cgindex[index[i]];
731 at1 = cgindex[index[i]+1];
732 for (j = at0; j < at1; j++)
734 /* We need to shift the coordinates */
735 rvec_add(x[j], shift, buf[n]);
742 for (i = 0; i < ind->nsend[nzone]; i++)
744 at0 = cgindex[index[i]];
745 at1 = cgindex[index[i]+1];
746 for (j = at0; j < at1; j++)
749 buf[n][XX] = x[j][XX] + shift[XX];
751 * This operation requires a special shift force
752 * treatment, which is performed in calc_vir.
754 buf[n][YY] = box[YY][YY] - x[j][YY];
755 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
767 rbuf = comm->vbuf2.v;
769 /* Send and receive the coordinates */
770 dd_sendrecv_rvec(dd, d, dddirBackward,
771 buf, ind->nsend[nzone+1],
772 rbuf, ind->nrecv[nzone+1]);
776 for (zone = 0; zone < nzone; zone++)
778 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
780 copy_rvec(rbuf[j], x[i]);
785 nat_tot += ind->nrecv[nzone+1];
791 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
793 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
794 int *index, *cgindex;
795 gmx_domdec_comm_t *comm;
796 gmx_domdec_comm_dim_t *cd;
797 gmx_domdec_ind_t *ind;
801 gmx_bool bShiftForcesNeedPbc, bScrew;
805 cgindex = dd->cgindex;
809 nzone = comm->zones.n/2;
810 nat_tot = dd->nat_tot;
811 for (d = dd->ndim-1; d >= 0; d--)
813 /* Only forces in domains near the PBC boundaries need to
814 consider PBC in the treatment of fshift */
815 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
816 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
817 if (fshift == NULL && !bScrew)
819 bShiftForcesNeedPbc = FALSE;
821 /* Determine which shift vector we need */
827 for (p = cd->np-1; p >= 0; p--)
830 nat_tot -= ind->nrecv[nzone+1];
837 sbuf = comm->vbuf2.v;
839 for (zone = 0; zone < nzone; zone++)
841 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
843 copy_rvec(f[i], sbuf[j]);
848 /* Communicate the forces */
849 dd_sendrecv_rvec(dd, d, dddirForward,
850 sbuf, ind->nrecv[nzone+1],
851 buf, ind->nsend[nzone+1]);
853 /* Add the received forces */
855 if (!bShiftForcesNeedPbc)
857 for (i = 0; i < ind->nsend[nzone]; i++)
859 at0 = cgindex[index[i]];
860 at1 = cgindex[index[i]+1];
861 for (j = at0; j < at1; j++)
863 rvec_inc(f[j], buf[n]);
870 /* fshift should always be defined if this function is
871 * called when bShiftForcesNeedPbc is true */
872 assert(NULL != fshift);
873 for (i = 0; i < ind->nsend[nzone]; i++)
875 at0 = cgindex[index[i]];
876 at1 = cgindex[index[i]+1];
877 for (j = at0; j < at1; j++)
879 rvec_inc(f[j], buf[n]);
880 /* Add this force to the shift force */
881 rvec_inc(fshift[is], buf[n]);
888 for (i = 0; i < ind->nsend[nzone]; i++)
890 at0 = cgindex[index[i]];
891 at1 = cgindex[index[i]+1];
892 for (j = at0; j < at1; j++)
894 /* Rotate the force */
895 f[j][XX] += buf[n][XX];
896 f[j][YY] -= buf[n][YY];
897 f[j][ZZ] -= buf[n][ZZ];
900 /* Add this force to the shift force */
901 rvec_inc(fshift[is], buf[n]);
912 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
914 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
915 int *index, *cgindex;
916 gmx_domdec_comm_t *comm;
917 gmx_domdec_comm_dim_t *cd;
918 gmx_domdec_ind_t *ind;
923 cgindex = dd->cgindex;
925 buf = &comm->vbuf.v[0][0];
928 nat_tot = dd->nat_home;
929 for (d = 0; d < dd->ndim; d++)
932 for (p = 0; p < cd->np; p++)
937 for (i = 0; i < ind->nsend[nzone]; i++)
939 at0 = cgindex[index[i]];
940 at1 = cgindex[index[i]+1];
941 for (j = at0; j < at1; j++)
954 rbuf = &comm->vbuf2.v[0][0];
956 /* Send and receive the coordinates */
957 dd_sendrecv_real(dd, d, dddirBackward,
958 buf, ind->nsend[nzone+1],
959 rbuf, ind->nrecv[nzone+1]);
963 for (zone = 0; zone < nzone; zone++)
965 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
972 nat_tot += ind->nrecv[nzone+1];
978 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
980 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
981 int *index, *cgindex;
982 gmx_domdec_comm_t *comm;
983 gmx_domdec_comm_dim_t *cd;
984 gmx_domdec_ind_t *ind;
989 cgindex = dd->cgindex;
991 buf = &comm->vbuf.v[0][0];
993 nzone = comm->zones.n/2;
994 nat_tot = dd->nat_tot;
995 for (d = dd->ndim-1; d >= 0; d--)
998 for (p = cd->np-1; p >= 0; p--)
1001 nat_tot -= ind->nrecv[nzone+1];
1008 sbuf = &comm->vbuf2.v[0][0];
1010 for (zone = 0; zone < nzone; zone++)
1012 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
1019 /* Communicate the forces */
1020 dd_sendrecv_real(dd, d, dddirForward,
1021 sbuf, ind->nrecv[nzone+1],
1022 buf, ind->nsend[nzone+1]);
1024 /* Add the received forces */
1026 for (i = 0; i < ind->nsend[nzone]; i++)
1028 at0 = cgindex[index[i]];
1029 at1 = cgindex[index[i]+1];
1030 for (j = at0; j < at1; j++)
1041 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1043 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",
1045 zone->min0, zone->max1,
1046 zone->mch0, zone->mch0,
1047 zone->p1_0, zone->p1_1);
1051 #define DDZONECOMM_MAXZONE 5
1052 #define DDZONECOMM_BUFSIZE 3
1054 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1055 int ddimind, int direction,
1056 gmx_ddzone_t *buf_s, int n_s,
1057 gmx_ddzone_t *buf_r, int n_r)
1059 #define ZBS DDZONECOMM_BUFSIZE
1060 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1061 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1064 for (i = 0; i < n_s; i++)
1066 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1067 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1068 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1069 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1070 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1071 vbuf_s[i*ZBS+1][2] = 0;
1072 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1073 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1074 vbuf_s[i*ZBS+2][2] = 0;
1077 dd_sendrecv_rvec(dd, ddimind, direction,
1081 for (i = 0; i < n_r; i++)
1083 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1084 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1085 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1086 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1087 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1088 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1089 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1095 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1096 rvec cell_ns_x0, rvec cell_ns_x1)
1098 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1100 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1101 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1102 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1103 rvec extr_s[2], extr_r[2];
1105 real dist_d, c = 0, det;
1106 gmx_domdec_comm_t *comm;
1107 gmx_bool bPBC, bUse;
1111 for (d = 1; d < dd->ndim; d++)
1114 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1115 zp->min0 = cell_ns_x0[dim];
1116 zp->max1 = cell_ns_x1[dim];
1117 zp->min1 = cell_ns_x1[dim];
1118 zp->mch0 = cell_ns_x0[dim];
1119 zp->mch1 = cell_ns_x1[dim];
1120 zp->p1_0 = cell_ns_x0[dim];
1121 zp->p1_1 = cell_ns_x1[dim];
1124 for (d = dd->ndim-2; d >= 0; d--)
1127 bPBC = (dim < ddbox->npbcdim);
1129 /* Use an rvec to store two reals */
1130 extr_s[d][0] = comm->cell_f0[d+1];
1131 extr_s[d][1] = comm->cell_f1[d+1];
1132 extr_s[d][2] = comm->cell_f1[d+1];
1135 /* Store the extremes in the backward sending buffer,
1136 * so the get updated separately from the forward communication.
1138 for (d1 = d; d1 < dd->ndim-1; d1++)
1140 /* We invert the order to be able to use the same loop for buf_e */
1141 buf_s[pos].min0 = extr_s[d1][1];
1142 buf_s[pos].max1 = extr_s[d1][0];
1143 buf_s[pos].min1 = extr_s[d1][2];
1144 buf_s[pos].mch0 = 0;
1145 buf_s[pos].mch1 = 0;
1146 /* Store the cell corner of the dimension we communicate along */
1147 buf_s[pos].p1_0 = comm->cell_x0[dim];
1148 buf_s[pos].p1_1 = 0;
1152 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1155 if (dd->ndim == 3 && d == 0)
1157 buf_s[pos] = comm->zone_d2[0][1];
1159 buf_s[pos] = comm->zone_d1[0];
1163 /* We only need to communicate the extremes
1164 * in the forward direction
1166 npulse = comm->cd[d].np;
1169 /* Take the minimum to avoid double communication */
1170 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1174 /* Without PBC we should really not communicate over
1175 * the boundaries, but implementing that complicates
1176 * the communication setup and therefore we simply
1177 * do all communication, but ignore some data.
1179 npulse_min = npulse;
1181 for (p = 0; p < npulse_min; p++)
1183 /* Communicate the extremes forward */
1184 bUse = (bPBC || dd->ci[dim] > 0);
1186 dd_sendrecv_rvec(dd, d, dddirForward,
1187 extr_s+d, dd->ndim-d-1,
1188 extr_r+d, dd->ndim-d-1);
1192 for (d1 = d; d1 < dd->ndim-1; d1++)
1194 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1195 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1196 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1202 for (p = 0; p < npulse; p++)
1204 /* Communicate all the zone information backward */
1205 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1207 dd_sendrecv_ddzone(dd, d, dddirBackward,
1214 for (d1 = d+1; d1 < dd->ndim; d1++)
1216 /* Determine the decrease of maximum required
1217 * communication height along d1 due to the distance along d,
1218 * this avoids a lot of useless atom communication.
1220 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1222 if (ddbox->tric_dir[dim])
1224 /* c is the off-diagonal coupling between the cell planes
1225 * along directions d and d1.
1227 c = ddbox->v[dim][dd->dim[d1]][dim];
1233 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1236 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1240 /* A negative value signals out of range */
1246 /* Accumulate the extremes over all pulses */
1247 for (i = 0; i < buf_size; i++)
1251 buf_e[i] = buf_r[i];
1257 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1258 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1259 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1262 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1270 if (bUse && dh[d1] >= 0)
1272 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1273 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1276 /* Copy the received buffer to the send buffer,
1277 * to pass the data through with the next pulse.
1279 buf_s[i] = buf_r[i];
1281 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1282 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1284 /* Store the extremes */
1287 for (d1 = d; d1 < dd->ndim-1; d1++)
1289 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1290 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1291 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1295 if (d == 1 || (d == 0 && dd->ndim == 3))
1297 for (i = d; i < 2; i++)
1299 comm->zone_d2[1-d][i] = buf_e[pos];
1305 comm->zone_d1[1] = buf_e[pos];
1315 for (i = 0; i < 2; i++)
1319 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1321 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1322 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1328 for (i = 0; i < 2; i++)
1330 for (j = 0; j < 2; j++)
1334 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1336 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1337 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1341 for (d = 1; d < dd->ndim; d++)
1343 comm->cell_f_max0[d] = extr_s[d-1][0];
1344 comm->cell_f_min1[d] = extr_s[d-1][1];
1347 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1348 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1353 static void dd_collect_cg(gmx_domdec_t *dd,
1354 t_state *state_local)
1356 gmx_domdec_master_t *ma = NULL;
1357 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1359 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1361 /* The master has the correct distribution */
1365 if (state_local->ddp_count == dd->ddp_count)
1367 /* The local state and DD are in sync, use the DD indices */
1368 ncg_home = dd->ncg_home;
1370 nat_home = dd->nat_home;
1372 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1374 /* The DD is out of sync with the local state, but we have stored
1375 * the cg indices with the local state, so we can use those.
1379 cgs_gl = &dd->comm->cgs_gl;
1381 ncg_home = state_local->ncg_gl;
1382 cg = state_local->cg_gl;
1384 for (i = 0; i < ncg_home; i++)
1386 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1391 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1405 /* Collect the charge group and atom counts on the master */
1406 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1411 for (i = 0; i < dd->nnodes; i++)
1413 ma->ncg[i] = ma->ibuf[2*i];
1414 ma->nat[i] = ma->ibuf[2*i+1];
1415 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1418 /* Make byte counts and indices */
1419 for (i = 0; i < dd->nnodes; i++)
1421 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1422 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1426 fprintf(debug, "Initial charge group distribution: ");
1427 for (i = 0; i < dd->nnodes; i++)
1429 fprintf(debug, " %d", ma->ncg[i]);
1431 fprintf(debug, "\n");
1435 /* Collect the charge group indices on the master */
1437 ncg_home*sizeof(int), cg,
1438 DDMASTER(dd) ? ma->ibuf : NULL,
1439 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1440 DDMASTER(dd) ? ma->cg : NULL);
1442 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1445 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1448 gmx_domdec_master_t *ma;
1449 int n, i, c, a, nalloc = 0;
1458 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1459 dd->rank, dd->mpi_comm_all);
1464 /* Copy the master coordinates to the global array */
1465 cgs_gl = &dd->comm->cgs_gl;
1467 n = DDMASTERRANK(dd);
1469 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1471 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1473 copy_rvec(lv[a++], v[c]);
1477 for (n = 0; n < dd->nnodes; n++)
1481 if (ma->nat[n] > nalloc)
1483 nalloc = over_alloc_dd(ma->nat[n]);
1484 srenew(buf, nalloc);
1487 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1488 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1491 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1493 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1495 copy_rvec(buf[a++], v[c]);
1504 static void get_commbuffer_counts(gmx_domdec_t *dd,
1505 int **counts, int **disps)
1507 gmx_domdec_master_t *ma;
1512 /* Make the rvec count and displacment arrays */
1514 *disps = ma->ibuf + dd->nnodes;
1515 for (n = 0; n < dd->nnodes; n++)
1517 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1518 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1522 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1525 gmx_domdec_master_t *ma;
1526 int *rcounts = NULL, *disps = NULL;
1535 get_commbuffer_counts(dd, &rcounts, &disps);
1540 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1544 cgs_gl = &dd->comm->cgs_gl;
1547 for (n = 0; n < dd->nnodes; n++)
1549 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1551 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1553 copy_rvec(buf[a++], v[c]);
1560 void dd_collect_vec(gmx_domdec_t *dd,
1561 t_state *state_local, rvec *lv, rvec *v)
1563 dd_collect_cg(dd, state_local);
1565 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1567 dd_collect_vec_sendrecv(dd, lv, v);
1571 dd_collect_vec_gatherv(dd, lv, v);
1576 void dd_collect_state(gmx_domdec_t *dd,
1577 t_state *state_local, t_state *state)
1581 nh = state->nhchainlength;
1585 for (i = 0; i < efptNR; i++)
1587 state->lambda[i] = state_local->lambda[i];
1589 state->fep_state = state_local->fep_state;
1590 state->veta = state_local->veta;
1591 state->vol0 = state_local->vol0;
1592 copy_mat(state_local->box, state->box);
1593 copy_mat(state_local->boxv, state->boxv);
1594 copy_mat(state_local->svir_prev, state->svir_prev);
1595 copy_mat(state_local->fvir_prev, state->fvir_prev);
1596 copy_mat(state_local->pres_prev, state->pres_prev);
1598 for (i = 0; i < state_local->ngtc; i++)
1600 for (j = 0; j < nh; j++)
1602 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1603 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1605 state->therm_integral[i] = state_local->therm_integral[i];
1607 for (i = 0; i < state_local->nnhpres; i++)
1609 for (j = 0; j < nh; j++)
1611 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1612 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1616 for (est = 0; est < estNR; est++)
1618 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1623 dd_collect_vec(dd, state_local, state_local->x, state->x);
1626 dd_collect_vec(dd, state_local, state_local->v, state->v);
1629 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1632 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1634 case estDISRE_INITF:
1635 case estDISRE_RM3TAV:
1636 case estORIRE_INITF:
1640 gmx_incons("Unknown state entry encountered in dd_collect_state");
1646 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1652 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1655 state->nalloc = over_alloc_dd(nalloc);
1657 for (est = 0; est < estNR; est++)
1659 if (EST_DISTR(est) && (state->flags & (1<<est)))
1664 srenew(state->x, state->nalloc);
1667 srenew(state->v, state->nalloc);
1670 srenew(state->sd_X, state->nalloc);
1673 srenew(state->cg_p, state->nalloc);
1675 case estDISRE_INITF:
1676 case estDISRE_RM3TAV:
1677 case estORIRE_INITF:
1679 /* No reallocation required */
1682 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1689 srenew(*f, state->nalloc);
1693 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1696 if (nalloc > fr->cg_nalloc)
1700 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1702 fr->cg_nalloc = over_alloc_dd(nalloc);
1703 srenew(fr->cginfo, fr->cg_nalloc);
1704 if (fr->cutoff_scheme == ecutsGROUP)
1706 srenew(fr->cg_cm, fr->cg_nalloc);
1709 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1711 /* We don't use charge groups, we use x in state to set up
1712 * the atom communication.
1714 dd_realloc_state(state, f, nalloc);
1718 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1721 gmx_domdec_master_t *ma;
1722 int n, i, c, a, nalloc = 0;
1729 for (n = 0; n < dd->nnodes; n++)
1733 if (ma->nat[n] > nalloc)
1735 nalloc = over_alloc_dd(ma->nat[n]);
1736 srenew(buf, nalloc);
1738 /* Use lv as a temporary buffer */
1740 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1742 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1744 copy_rvec(v[c], buf[a++]);
1747 if (a != ma->nat[n])
1749 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1754 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1755 DDRANK(dd, n), n, dd->mpi_comm_all);
1760 n = DDMASTERRANK(dd);
1762 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1764 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1766 copy_rvec(v[c], lv[a++]);
1773 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1774 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1779 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1782 gmx_domdec_master_t *ma;
1783 int *scounts = NULL, *disps = NULL;
1791 get_commbuffer_counts(dd, &scounts, &disps);
1795 for (n = 0; n < dd->nnodes; n++)
1797 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1799 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1801 copy_rvec(v[c], buf[a++]);
1807 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1810 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1812 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1814 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1818 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1822 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1825 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1826 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1827 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1829 if (dfhist->nlambda > 0)
1831 int nlam = dfhist->nlambda;
1832 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1833 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1834 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1835 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1836 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1837 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1839 for (i = 0; i < nlam; i++)
1841 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1842 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1843 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1844 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1845 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1846 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1851 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1852 t_state *state, t_state *state_local,
1857 nh = state->nhchainlength;
1861 for (i = 0; i < efptNR; i++)
1863 state_local->lambda[i] = state->lambda[i];
1865 state_local->fep_state = state->fep_state;
1866 state_local->veta = state->veta;
1867 state_local->vol0 = state->vol0;
1868 copy_mat(state->box, state_local->box);
1869 copy_mat(state->box_rel, state_local->box_rel);
1870 copy_mat(state->boxv, state_local->boxv);
1871 copy_mat(state->svir_prev, state_local->svir_prev);
1872 copy_mat(state->fvir_prev, state_local->fvir_prev);
1873 copy_df_history(&state_local->dfhist, &state->dfhist);
1874 for (i = 0; i < state_local->ngtc; i++)
1876 for (j = 0; j < nh; j++)
1878 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1879 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1881 state_local->therm_integral[i] = state->therm_integral[i];
1883 for (i = 0; i < state_local->nnhpres; i++)
1885 for (j = 0; j < nh; j++)
1887 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1888 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1892 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1893 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1894 dd_bcast(dd, sizeof(real), &state_local->veta);
1895 dd_bcast(dd, sizeof(real), &state_local->vol0);
1896 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1897 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1898 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1899 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1900 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1901 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1902 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1903 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1904 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1905 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1907 /* communicate df_history -- required for restarting from checkpoint */
1908 dd_distribute_dfhist(dd, &state_local->dfhist);
1910 if (dd->nat_home > state_local->nalloc)
1912 dd_realloc_state(state_local, f, dd->nat_home);
1914 for (i = 0; i < estNR; i++)
1916 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1921 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1924 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1927 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1930 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1932 case estDISRE_INITF:
1933 case estDISRE_RM3TAV:
1934 case estORIRE_INITF:
1936 /* Not implemented yet */
1939 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1945 static char dim2char(int dim)
1951 case XX: c = 'X'; break;
1952 case YY: c = 'Y'; break;
1953 case ZZ: c = 'Z'; break;
1954 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1960 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1961 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1963 rvec grid_s[2], *grid_r = NULL, cx, r;
1964 char fname[STRLEN], buf[22];
1966 int a, i, d, z, y, x;
1970 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1971 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1975 snew(grid_r, 2*dd->nnodes);
1978 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1982 for (d = 0; d < DIM; d++)
1984 for (i = 0; i < DIM; i++)
1992 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1994 tric[d][i] = box[i][d]/box[i][i];
2003 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2004 out = gmx_fio_fopen(fname, "w");
2005 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2007 for (i = 0; i < dd->nnodes; i++)
2009 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2010 for (d = 0; d < DIM; d++)
2012 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2014 for (z = 0; z < 2; z++)
2016 for (y = 0; y < 2; y++)
2018 for (x = 0; x < 2; x++)
2020 cx[XX] = grid_r[i*2+x][XX];
2021 cx[YY] = grid_r[i*2+y][YY];
2022 cx[ZZ] = grid_r[i*2+z][ZZ];
2024 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2025 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2029 for (d = 0; d < DIM; d++)
2031 for (x = 0; x < 4; x++)
2035 case 0: y = 1 + i*8 + 2*x; break;
2036 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2037 case 2: y = 1 + i*8 + x; break;
2039 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2043 gmx_fio_fclose(out);
2048 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2049 gmx_mtop_t *mtop, t_commrec *cr,
2050 int natoms, rvec x[], matrix box)
2052 char fname[STRLEN], buf[22];
2054 int i, ii, resnr, c;
2055 char *atomname, *resname;
2062 natoms = dd->comm->nat[ddnatVSITE];
2065 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2067 out = gmx_fio_fopen(fname, "w");
2069 fprintf(out, "TITLE %s\n", title);
2070 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2071 for (i = 0; i < natoms; i++)
2073 ii = dd->gatindex[i];
2074 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2075 if (i < dd->comm->nat[ddnatZONE])
2078 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2084 else if (i < dd->comm->nat[ddnatVSITE])
2086 b = dd->comm->zones.n;
2090 b = dd->comm->zones.n + 1;
2092 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2093 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2095 fprintf(out, "TER\n");
2097 gmx_fio_fclose(out);
2100 real dd_cutoff_multibody(const gmx_domdec_t *dd)
2102 gmx_domdec_comm_t *comm;
2109 if (comm->bInterCGBondeds)
2111 if (comm->cutoff_mbody > 0)
2113 r = comm->cutoff_mbody;
2117 /* cutoff_mbody=0 means we do not have DLB */
2118 r = comm->cellsize_min[dd->dim[0]];
2119 for (di = 1; di < dd->ndim; di++)
2121 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2123 if (comm->bBondComm)
2125 r = std::max(r, comm->cutoff_mbody);
2129 r = std::min(r, comm->cutoff);
2137 real dd_cutoff_twobody(const gmx_domdec_t *dd)
2141 r_mb = dd_cutoff_multibody(dd);
2143 return std::max(dd->comm->cutoff, r_mb);
2147 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2151 nc = dd->nc[dd->comm->cartpmedim];
2152 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2153 copy_ivec(coord, coord_pme);
2154 coord_pme[dd->comm->cartpmedim] =
2155 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2158 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2160 /* Here we assign a PME node to communicate with this DD node
2161 * by assuming that the major index of both is x.
2162 * We add cr->npmenodes/2 to obtain an even distribution.
2164 return (ddindex*npme + npme/2)/ndd;
2167 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2169 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2172 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2174 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2177 static int *dd_pmenodes(t_commrec *cr)
2182 snew(pmenodes, cr->npmenodes);
2184 for (i = 0; i < cr->dd->nnodes; i++)
2186 p0 = cr_ddindex2pmeindex(cr, i);
2187 p1 = cr_ddindex2pmeindex(cr, i+1);
2188 if (i+1 == cr->dd->nnodes || p1 > p0)
2192 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2194 pmenodes[n] = i + 1 + n;
2202 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2210 if (dd->comm->bCartesian) {
2211 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2212 dd_coords2pmecoords(dd,coords,coords_pme);
2213 copy_ivec(dd->ntot,nc);
2214 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2215 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2217 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2219 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2225 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2230 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2232 gmx_domdec_comm_t *comm;
2234 int ddindex, nodeid = -1;
2236 comm = cr->dd->comm;
2241 if (comm->bCartesianPP_PME)
2244 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2249 ddindex = dd_index(cr->dd->nc, coords);
2250 if (comm->bCartesianPP)
2252 nodeid = comm->ddindex2simnodeid[ddindex];
2258 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2270 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2273 gmx_domdec_comm_t *comm;
2280 /* This assumes a uniform x domain decomposition grid cell size */
2281 if (comm->bCartesianPP_PME)
2284 ivec coord, coord_pme;
2285 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2286 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2288 /* This is a PP node */
2289 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2290 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2294 else if (comm->bCartesianPP)
2296 if (sim_nodeid < dd->nnodes)
2298 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2303 /* This assumes DD cells with identical x coordinates
2304 * are numbered sequentially.
2306 if (dd->comm->pmenodes == NULL)
2308 if (sim_nodeid < dd->nnodes)
2310 /* The DD index equals the nodeid */
2311 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2317 while (sim_nodeid > dd->comm->pmenodes[i])
2321 if (sim_nodeid < dd->comm->pmenodes[i])
2323 pmenode = dd->comm->pmenodes[i];
2331 void get_pme_nnodes(const gmx_domdec_t *dd,
2332 int *npmenodes_x, int *npmenodes_y)
2336 *npmenodes_x = dd->comm->npmenodes_x;
2337 *npmenodes_y = dd->comm->npmenodes_y;
2346 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2347 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2351 ivec coord, coord_pme;
2355 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2358 for (x = 0; x < dd->nc[XX]; x++)
2360 for (y = 0; y < dd->nc[YY]; y++)
2362 for (z = 0; z < dd->nc[ZZ]; z++)
2364 if (dd->comm->bCartesianPP_PME)
2369 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2370 if (dd->ci[XX] == coord_pme[XX] &&
2371 dd->ci[YY] == coord_pme[YY] &&
2372 dd->ci[ZZ] == coord_pme[ZZ])
2374 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2379 /* The slab corresponds to the nodeid in the PME group */
2380 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2382 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2389 /* The last PP-only node is the peer node */
2390 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2394 fprintf(debug, "Receive coordinates from PP ranks:");
2395 for (x = 0; x < *nmy_ddnodes; x++)
2397 fprintf(debug, " %d", (*my_ddnodes)[x]);
2399 fprintf(debug, "\n");
2403 static gmx_bool receive_vir_ener(t_commrec *cr)
2405 gmx_domdec_comm_t *comm;
2410 if (cr->npmenodes < cr->dd->nnodes)
2412 comm = cr->dd->comm;
2413 if (comm->bCartesianPP_PME)
2415 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2418 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2419 coords[comm->cartpmedim]++;
2420 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2423 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2424 if (dd_simnode2pmenode(cr, rank) == pmenode)
2426 /* This is not the last PP node for pmenode */
2434 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2435 if (cr->sim_nodeid+1 < cr->nnodes &&
2436 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2438 /* This is not the last PP node for pmenode */
2447 static void set_zones_ncg_home(gmx_domdec_t *dd)
2449 gmx_domdec_zones_t *zones;
2452 zones = &dd->comm->zones;
2454 zones->cg_range[0] = 0;
2455 for (i = 1; i < zones->n+1; i++)
2457 zones->cg_range[i] = dd->ncg_home;
2459 /* zone_ncg1[0] should always be equal to ncg_home */
2460 dd->comm->zone_ncg1[0] = dd->ncg_home;
2463 static void rebuild_cgindex(gmx_domdec_t *dd,
2464 const int *gcgs_index, t_state *state)
2466 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2469 dd_cg_gl = dd->index_gl;
2470 cgindex = dd->cgindex;
2473 for (i = 0; i < state->ncg_gl; i++)
2477 dd_cg_gl[i] = cg_gl;
2478 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2482 dd->ncg_home = state->ncg_gl;
2485 set_zones_ncg_home(dd);
2488 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2490 while (cg >= cginfo_mb->cg_end)
2495 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2498 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2499 t_forcerec *fr, char *bLocalCG)
2501 cginfo_mb_t *cginfo_mb;
2507 cginfo_mb = fr->cginfo_mb;
2508 cginfo = fr->cginfo;
2510 for (cg = cg0; cg < cg1; cg++)
2512 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2516 if (bLocalCG != NULL)
2518 for (cg = cg0; cg < cg1; cg++)
2520 bLocalCG[index_gl[cg]] = TRUE;
2525 static void make_dd_indices(gmx_domdec_t *dd,
2526 const int *gcgs_index, int cg_start)
2528 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2529 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2532 if (dd->nat_tot > dd->gatindex_nalloc)
2534 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2535 srenew(dd->gatindex, dd->gatindex_nalloc);
2538 nzone = dd->comm->zones.n;
2539 zone2cg = dd->comm->zones.cg_range;
2540 zone_ncg1 = dd->comm->zone_ncg1;
2541 index_gl = dd->index_gl;
2542 gatindex = dd->gatindex;
2543 bCGs = dd->comm->bCGs;
2545 if (zone2cg[1] != dd->ncg_home)
2547 gmx_incons("dd->ncg_zone is not up to date");
2550 /* Make the local to global and global to local atom index */
2551 a = dd->cgindex[cg_start];
2552 for (zone = 0; zone < nzone; zone++)
2560 cg0 = zone2cg[zone];
2562 cg1 = zone2cg[zone+1];
2563 cg1_p1 = cg0 + zone_ncg1[zone];
2565 for (cg = cg0; cg < cg1; cg++)
2570 /* Signal that this cg is from more than one pulse away */
2573 cg_gl = index_gl[cg];
2576 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2579 ga2la_set(dd->ga2la, a_gl, a, zone1);
2585 gatindex[a] = cg_gl;
2586 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2593 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2599 if (bLocalCG == NULL)
2603 for (i = 0; i < dd->ncg_tot; i++)
2605 if (!bLocalCG[dd->index_gl[i]])
2608 "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);
2613 for (i = 0; i < ncg_sys; i++)
2620 if (ngl != dd->ncg_tot)
2622 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);
2629 static void check_index_consistency(gmx_domdec_t *dd,
2630 int natoms_sys, int ncg_sys,
2633 int nerr, ngl, i, a, cell;
2638 if (dd->comm->DD_debug > 1)
2640 snew(have, natoms_sys);
2641 for (a = 0; a < dd->nat_tot; a++)
2643 if (have[dd->gatindex[a]] > 0)
2645 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);
2649 have[dd->gatindex[a]] = a + 1;
2655 snew(have, dd->nat_tot);
2658 for (i = 0; i < natoms_sys; i++)
2660 if (ga2la_get(dd->ga2la, i, &a, &cell))
2662 if (a >= dd->nat_tot)
2664 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);
2670 if (dd->gatindex[a] != i)
2672 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);
2679 if (ngl != dd->nat_tot)
2682 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2683 dd->rank, where, ngl, dd->nat_tot);
2685 for (a = 0; a < dd->nat_tot; a++)
2690 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2691 dd->rank, where, a+1, dd->gatindex[a]+1);
2696 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2700 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2701 dd->rank, where, nerr);
2705 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2712 /* Clear the whole list without searching */
2713 ga2la_clear(dd->ga2la);
2717 for (i = a_start; i < dd->nat_tot; i++)
2719 ga2la_del(dd->ga2la, dd->gatindex[i]);
2723 bLocalCG = dd->comm->bLocalCG;
2726 for (i = cg_start; i < dd->ncg_tot; i++)
2728 bLocalCG[dd->index_gl[i]] = FALSE;
2732 dd_clear_local_vsite_indices(dd);
2734 if (dd->constraints)
2736 dd_clear_local_constraint_indices(dd);
2740 /* This function should be used for moving the domain boudaries during DLB,
2741 * for obtaining the minimum cell size. It checks the initially set limit
2742 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2743 * and, possibly, a longer cut-off limit set for PME load balancing.
2745 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2749 cellsize_min = comm->cellsize_min[dim];
2751 if (!comm->bVacDLBNoLimit)
2753 /* The cut-off might have changed, e.g. by PME load balacning,
2754 * from the value used to set comm->cellsize_min, so check it.
2756 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2758 if (comm->bPMELoadBalDLBLimits)
2760 /* Check for the cut-off limit set by the PME load balancing */
2761 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2765 return cellsize_min;
2768 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2771 real grid_jump_limit;
2773 /* The distance between the boundaries of cells at distance
2774 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2775 * and by the fact that cells should not be shifted by more than
2776 * half their size, such that cg's only shift by one cell
2777 * at redecomposition.
2779 grid_jump_limit = comm->cellsize_limit;
2780 if (!comm->bVacDLBNoLimit)
2782 if (comm->bPMELoadBalDLBLimits)
2784 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2786 grid_jump_limit = std::max(grid_jump_limit,
2787 cutoff/comm->cd[dim_ind].np);
2790 return grid_jump_limit;
2793 static gmx_bool check_grid_jump(gmx_int64_t step,
2799 gmx_domdec_comm_t *comm;
2808 for (d = 1; d < dd->ndim; d++)
2811 limit = grid_jump_limit(comm, cutoff, d);
2812 bfac = ddbox->box_size[dim];
2813 if (ddbox->tric_dir[dim])
2815 bfac *= ddbox->skew_fac[dim];
2817 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2818 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2826 /* This error should never be triggered under normal
2827 * circumstances, but you never know ...
2829 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.",
2830 gmx_step_str(step, buf),
2831 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2839 static int dd_load_count(gmx_domdec_comm_t *comm)
2841 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2844 static float dd_force_load(gmx_domdec_comm_t *comm)
2851 if (comm->eFlop > 1)
2853 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2858 load = comm->cycl[ddCyclF];
2859 if (comm->cycl_n[ddCyclF] > 1)
2861 /* Subtract the maximum of the last n cycle counts
2862 * to get rid of possible high counts due to other sources,
2863 * for instance system activity, that would otherwise
2864 * affect the dynamic load balancing.
2866 load -= comm->cycl_max[ddCyclF];
2870 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2872 float gpu_wait, gpu_wait_sum;
2874 gpu_wait = comm->cycl[ddCyclWaitGPU];
2875 if (comm->cycl_n[ddCyclF] > 1)
2877 /* We should remove the WaitGPU time of the same MD step
2878 * as the one with the maximum F time, since the F time
2879 * and the wait time are not independent.
2880 * Furthermore, the step for the max F time should be chosen
2881 * the same on all ranks that share the same GPU.
2882 * But to keep the code simple, we remove the average instead.
2883 * The main reason for artificially long times at some steps
2884 * is spurious CPU activity or MPI time, so we don't expect
2885 * that changes in the GPU wait time matter a lot here.
2887 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2889 /* Sum the wait times over the ranks that share the same GPU */
2890 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2891 comm->mpi_comm_gpu_shared);
2892 /* Replace the wait time by the average over the ranks */
2893 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2901 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2903 gmx_domdec_comm_t *comm;
2908 snew(*dim_f, dd->nc[dim]+1);
2910 for (i = 1; i < dd->nc[dim]; i++)
2912 if (comm->slb_frac[dim])
2914 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2918 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2921 (*dim_f)[dd->nc[dim]] = 1;
2924 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2926 int pmeindex, slab, nso, i;
2929 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2935 ddpme->dim = dimind;
2937 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2939 ddpme->nslab = (ddpme->dim == 0 ?
2940 dd->comm->npmenodes_x :
2941 dd->comm->npmenodes_y);
2943 if (ddpme->nslab <= 1)
2948 nso = dd->comm->npmenodes/ddpme->nslab;
2949 /* Determine for each PME slab the PP location range for dimension dim */
2950 snew(ddpme->pp_min, ddpme->nslab);
2951 snew(ddpme->pp_max, ddpme->nslab);
2952 for (slab = 0; slab < ddpme->nslab; slab++)
2954 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2955 ddpme->pp_max[slab] = 0;
2957 for (i = 0; i < dd->nnodes; i++)
2959 ddindex2xyz(dd->nc, i, xyz);
2960 /* For y only use our y/z slab.
2961 * This assumes that the PME x grid size matches the DD grid size.
2963 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2965 pmeindex = ddindex2pmeindex(dd, i);
2968 slab = pmeindex/nso;
2972 slab = pmeindex % ddpme->nslab;
2974 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2975 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2979 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2982 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2984 if (dd->comm->ddpme[0].dim == XX)
2986 return dd->comm->ddpme[0].maxshift;
2994 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2996 if (dd->comm->ddpme[0].dim == YY)
2998 return dd->comm->ddpme[0].maxshift;
3000 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3002 return dd->comm->ddpme[1].maxshift;
3010 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3011 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3013 gmx_domdec_comm_t *comm;
3016 real range, pme_boundary;
3020 nc = dd->nc[ddpme->dim];
3023 if (!ddpme->dim_match)
3025 /* PP decomposition is not along dim: the worst situation */
3028 else if (ns <= 3 || (bUniform && ns == nc))
3030 /* The optimal situation */
3035 /* We need to check for all pme nodes which nodes they
3036 * could possibly need to communicate with.
3038 xmin = ddpme->pp_min;
3039 xmax = ddpme->pp_max;
3040 /* Allow for atoms to be maximally 2/3 times the cut-off
3041 * out of their DD cell. This is a reasonable balance between
3042 * between performance and support for most charge-group/cut-off
3045 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3046 /* Avoid extra communication when we are exactly at a boundary */
3050 for (s = 0; s < ns; s++)
3052 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3053 pme_boundary = (real)s/ns;
3056 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3058 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3062 pme_boundary = (real)(s+1)/ns;
3065 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3067 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3074 ddpme->maxshift = sh;
3078 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3079 ddpme->dim, ddpme->maxshift);
3083 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3087 for (d = 0; d < dd->ndim; d++)
3090 if (dim < ddbox->nboundeddim &&
3091 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3092 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3094 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",
3095 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3096 dd->nc[dim], dd->comm->cellsize_limit);
3102 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3105 /* Set the domain boundaries. Use for static (or no) load balancing,
3106 * and also for the starting state for dynamic load balancing.
3107 * setmode determine if and where the boundaries are stored, use enum above.
3108 * Returns the number communication pulses in npulse.
3110 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3111 int setmode, ivec npulse)
3113 gmx_domdec_comm_t *comm;
3116 real *cell_x, cell_dx, cellsize;
3120 for (d = 0; d < DIM; d++)
3122 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3124 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3127 cell_dx = ddbox->box_size[d]/dd->nc[d];
3130 case setcellsizeslbMASTER:
3131 for (j = 0; j < dd->nc[d]+1; j++)
3133 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3136 case setcellsizeslbLOCAL:
3137 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3138 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3143 cellsize = cell_dx*ddbox->skew_fac[d];
3144 while (cellsize*npulse[d] < comm->cutoff)
3148 cellsize_min[d] = cellsize;
3152 /* Statically load balanced grid */
3153 /* Also when we are not doing a master distribution we determine
3154 * all cell borders in a loop to obtain identical values
3155 * to the master distribution case and to determine npulse.
3157 if (setmode == setcellsizeslbMASTER)
3159 cell_x = dd->ma->cell_x[d];
3163 snew(cell_x, dd->nc[d]+1);
3165 cell_x[0] = ddbox->box0[d];
3166 for (j = 0; j < dd->nc[d]; j++)
3168 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3169 cell_x[j+1] = cell_x[j] + cell_dx;
3170 cellsize = cell_dx*ddbox->skew_fac[d];
3171 while (cellsize*npulse[d] < comm->cutoff &&
3172 npulse[d] < dd->nc[d]-1)
3176 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3178 if (setmode == setcellsizeslbLOCAL)
3180 comm->cell_x0[d] = cell_x[dd->ci[d]];
3181 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3183 if (setmode != setcellsizeslbMASTER)
3188 /* The following limitation is to avoid that a cell would receive
3189 * some of its own home charge groups back over the periodic boundary.
3190 * Double charge groups cause trouble with the global indices.
3192 if (d < ddbox->npbcdim &&
3193 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3195 char error_string[STRLEN];
3197 sprintf(error_string,
3198 "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",
3199 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3201 dd->nc[d], dd->nc[d],
3202 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3204 if (setmode == setcellsizeslbLOCAL)
3206 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3210 gmx_fatal(FARGS, error_string);
3215 if (!comm->bDynLoadBal)
3217 copy_rvec(cellsize_min, comm->cellsize_min);
3220 for (d = 0; d < comm->npmedecompdim; d++)
3222 set_pme_maxshift(dd, &comm->ddpme[d],
3223 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3224 comm->ddpme[d].slb_dim_f);
3229 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3230 int d, int dim, gmx_domdec_root_t *root,
3232 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3234 gmx_domdec_comm_t *comm;
3235 int ncd, i, j, nmin, nmin_old;
3236 gmx_bool bLimLo, bLimHi;
3238 real fac, halfway, cellsize_limit_f_i, region_size;
3239 gmx_bool bPBC, bLastHi = FALSE;
3240 int nrange[] = {range[0], range[1]};
3242 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3248 bPBC = (dim < ddbox->npbcdim);
3250 cell_size = root->buf_ncd;
3254 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3257 /* First we need to check if the scaling does not make cells
3258 * smaller than the smallest allowed size.
3259 * We need to do this iteratively, since if a cell is too small,
3260 * it needs to be enlarged, which makes all the other cells smaller,
3261 * which could in turn make another cell smaller than allowed.
3263 for (i = range[0]; i < range[1]; i++)
3265 root->bCellMin[i] = FALSE;
3271 /* We need the total for normalization */
3273 for (i = range[0]; i < range[1]; i++)
3275 if (root->bCellMin[i] == FALSE)
3277 fac += cell_size[i];
3280 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3281 /* Determine the cell boundaries */
3282 for (i = range[0]; i < range[1]; i++)
3284 if (root->bCellMin[i] == FALSE)
3286 cell_size[i] *= fac;
3287 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3289 cellsize_limit_f_i = 0;
3293 cellsize_limit_f_i = cellsize_limit_f;
3295 if (cell_size[i] < cellsize_limit_f_i)
3297 root->bCellMin[i] = TRUE;
3298 cell_size[i] = cellsize_limit_f_i;
3302 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3305 while (nmin > nmin_old);
3308 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3309 /* For this check we should not use DD_CELL_MARGIN,
3310 * but a slightly smaller factor,
3311 * since rounding could get use below the limit.
3313 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3316 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",
3317 gmx_step_str(step, buf),
3318 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3319 ncd, comm->cellsize_min[dim]);
3322 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3326 /* Check if the boundary did not displace more than halfway
3327 * each of the cells it bounds, as this could cause problems,
3328 * especially when the differences between cell sizes are large.
3329 * If changes are applied, they will not make cells smaller
3330 * than the cut-off, as we check all the boundaries which
3331 * might be affected by a change and if the old state was ok,
3332 * the cells will at most be shrunk back to their old size.
3334 for (i = range[0]+1; i < range[1]; i++)
3336 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3337 if (root->cell_f[i] < halfway)
3339 root->cell_f[i] = halfway;
3340 /* Check if the change also causes shifts of the next boundaries */
3341 for (j = i+1; j < range[1]; j++)
3343 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3345 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3349 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3350 if (root->cell_f[i] > halfway)
3352 root->cell_f[i] = halfway;
3353 /* Check if the change also causes shifts of the next boundaries */
3354 for (j = i-1; j >= range[0]+1; j--)
3356 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3358 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3365 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3366 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3367 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3368 * for a and b nrange is used */
3371 /* Take care of the staggering of the cell boundaries */
3374 for (i = range[0]; i < range[1]; i++)
3376 root->cell_f_max0[i] = root->cell_f[i];
3377 root->cell_f_min1[i] = root->cell_f[i+1];
3382 for (i = range[0]+1; i < range[1]; i++)
3384 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3385 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3386 if (bLimLo && bLimHi)
3388 /* Both limits violated, try the best we can */
3389 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3390 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3391 nrange[0] = range[0];
3393 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3396 nrange[1] = range[1];
3397 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3403 /* root->cell_f[i] = root->bound_min[i]; */
3404 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3407 else if (bLimHi && !bLastHi)
3410 if (nrange[1] < range[1]) /* found a LimLo before */
3412 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3413 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3414 nrange[0] = nrange[1];
3416 root->cell_f[i] = root->bound_max[i];
3418 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3420 nrange[1] = range[1];
3423 if (nrange[1] < range[1]) /* found last a LimLo */
3425 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3426 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3427 nrange[0] = nrange[1];
3428 nrange[1] = range[1];
3429 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3431 else if (nrange[0] > range[0]) /* found at least one LimHi */
3433 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3440 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3441 int d, int dim, gmx_domdec_root_t *root,
3442 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3443 gmx_bool bUniform, gmx_int64_t step)
3445 gmx_domdec_comm_t *comm;
3446 int ncd, d1, i, pos;
3448 real load_aver, load_i, imbalance, change, change_max, sc;
3449 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3453 int range[] = { 0, 0 };
3457 /* Convert the maximum change from the input percentage to a fraction */
3458 change_limit = comm->dlb_scale_lim*0.01;
3462 bPBC = (dim < ddbox->npbcdim);
3464 cell_size = root->buf_ncd;
3466 /* Store the original boundaries */
3467 for (i = 0; i < ncd+1; i++)
3469 root->old_cell_f[i] = root->cell_f[i];
3473 for (i = 0; i < ncd; i++)
3475 cell_size[i] = 1.0/ncd;
3478 else if (dd_load_count(comm) > 0)
3480 load_aver = comm->load[d].sum_m/ncd;
3482 for (i = 0; i < ncd; i++)
3484 /* Determine the relative imbalance of cell i */
3485 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3486 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3487 /* Determine the change of the cell size using underrelaxation */
3488 change = -relax*imbalance;
3489 change_max = std::max(change_max, std::max(change, -change));
3491 /* Limit the amount of scaling.
3492 * We need to use the same rescaling for all cells in one row,
3493 * otherwise the load balancing might not converge.
3496 if (change_max > change_limit)
3498 sc *= change_limit/change_max;
3500 for (i = 0; i < ncd; i++)
3502 /* Determine the relative imbalance of cell i */
3503 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3504 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3505 /* Determine the change of the cell size using underrelaxation */
3506 change = -sc*imbalance;
3507 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3511 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3512 cellsize_limit_f *= DD_CELL_MARGIN;
3513 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3514 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3515 if (ddbox->tric_dir[dim])
3517 cellsize_limit_f /= ddbox->skew_fac[dim];
3518 dist_min_f /= ddbox->skew_fac[dim];
3520 if (bDynamicBox && d > 0)
3522 dist_min_f *= DD_PRES_SCALE_MARGIN;
3524 if (d > 0 && !bUniform)
3526 /* Make sure that the grid is not shifted too much */
3527 for (i = 1; i < ncd; i++)
3529 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3531 gmx_incons("Inconsistent DD boundary staggering limits!");
3533 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3534 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3537 root->bound_min[i] += 0.5*space;
3539 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3540 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3543 root->bound_max[i] += 0.5*space;
3548 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3550 root->cell_f_max0[i-1] + dist_min_f,
3551 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3552 root->cell_f_min1[i] - dist_min_f);
3557 root->cell_f[0] = 0;
3558 root->cell_f[ncd] = 1;
3559 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3562 /* After the checks above, the cells should obey the cut-off
3563 * restrictions, but it does not hurt to check.
3565 for (i = 0; i < ncd; i++)
3569 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3570 dim, i, root->cell_f[i], root->cell_f[i+1]);
3573 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3574 root->cell_f[i+1] - root->cell_f[i] <
3575 cellsize_limit_f/DD_CELL_MARGIN)
3579 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3580 gmx_step_str(step, buf), dim2char(dim), i,
3581 (root->cell_f[i+1] - root->cell_f[i])
3582 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3587 /* Store the cell boundaries of the lower dimensions at the end */
3588 for (d1 = 0; d1 < d; d1++)
3590 root->cell_f[pos++] = comm->cell_f0[d1];
3591 root->cell_f[pos++] = comm->cell_f1[d1];
3594 if (d < comm->npmedecompdim)
3596 /* The master determines the maximum shift for
3597 * the coordinate communication between separate PME nodes.
3599 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3601 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3604 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3608 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3609 gmx_ddbox_t *ddbox, int dimind)
3611 gmx_domdec_comm_t *comm;
3616 /* Set the cell dimensions */
3617 dim = dd->dim[dimind];
3618 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3619 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3620 if (dim >= ddbox->nboundeddim)
3622 comm->cell_x0[dim] += ddbox->box0[dim];
3623 comm->cell_x1[dim] += ddbox->box0[dim];
3627 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3628 int d, int dim, real *cell_f_row,
3631 gmx_domdec_comm_t *comm;
3637 /* Each node would only need to know two fractions,
3638 * but it is probably cheaper to broadcast the whole array.
3640 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3641 0, comm->mpi_comm_load[d]);
3643 /* Copy the fractions for this dimension from the buffer */
3644 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3645 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3646 /* The whole array was communicated, so set the buffer position */
3647 pos = dd->nc[dim] + 1;
3648 for (d1 = 0; d1 <= d; d1++)
3652 /* Copy the cell fractions of the lower dimensions */
3653 comm->cell_f0[d1] = cell_f_row[pos++];
3654 comm->cell_f1[d1] = cell_f_row[pos++];
3656 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3658 /* Convert the communicated shift from float to int */
3659 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3662 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3666 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3667 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3668 gmx_bool bUniform, gmx_int64_t step)
3670 gmx_domdec_comm_t *comm;
3672 gmx_bool bRowMember, bRowRoot;
3677 for (d = 0; d < dd->ndim; d++)
3682 for (d1 = d; d1 < dd->ndim; d1++)
3684 if (dd->ci[dd->dim[d1]] > 0)
3697 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3698 ddbox, bDynamicBox, bUniform, step);
3699 cell_f_row = comm->root[d]->cell_f;
3703 cell_f_row = comm->cell_f_row;
3705 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3710 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3714 /* This function assumes the box is static and should therefore
3715 * not be called when the box has changed since the last
3716 * call to dd_partition_system.
3718 for (d = 0; d < dd->ndim; d++)
3720 relative_to_absolute_cell_bounds(dd, ddbox, d);
3726 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3727 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3728 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3729 gmx_wallcycle_t wcycle)
3731 gmx_domdec_comm_t *comm;
3738 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3739 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3740 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3742 else if (bDynamicBox)
3744 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3747 /* Set the dimensions for which no DD is used */
3748 for (dim = 0; dim < DIM; dim++)
3750 if (dd->nc[dim] == 1)
3752 comm->cell_x0[dim] = 0;
3753 comm->cell_x1[dim] = ddbox->box_size[dim];
3754 if (dim >= ddbox->nboundeddim)
3756 comm->cell_x0[dim] += ddbox->box0[dim];
3757 comm->cell_x1[dim] += ddbox->box0[dim];
3763 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3766 gmx_domdec_comm_dim_t *cd;
3768 for (d = 0; d < dd->ndim; d++)
3770 cd = &dd->comm->cd[d];
3771 np = npulse[dd->dim[d]];
3772 if (np > cd->np_nalloc)
3776 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3777 dim2char(dd->dim[d]), np);
3779 if (DDMASTER(dd) && cd->np_nalloc > 0)
3781 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3783 srenew(cd->ind, np);
3784 for (i = cd->np_nalloc; i < np; i++)
3786 cd->ind[i].index = NULL;
3787 cd->ind[i].nalloc = 0;
3796 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3797 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3798 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3799 gmx_wallcycle_t wcycle)
3801 gmx_domdec_comm_t *comm;
3807 /* Copy the old cell boundaries for the cg displacement check */
3808 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3809 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3811 if (comm->bDynLoadBal)
3815 check_box_size(dd, ddbox);
3817 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3821 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3822 realloc_comm_ind(dd, npulse);
3827 for (d = 0; d < DIM; d++)
3829 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3830 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3835 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3837 rvec cell_ns_x0, rvec cell_ns_x1,
3840 gmx_domdec_comm_t *comm;
3845 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3847 dim = dd->dim[dim_ind];
3849 /* Without PBC we don't have restrictions on the outer cells */
3850 if (!(dim >= ddbox->npbcdim &&
3851 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3852 comm->bDynLoadBal &&
3853 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3854 comm->cellsize_min[dim])
3857 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",
3858 gmx_step_str(step, buf), dim2char(dim),
3859 comm->cell_x1[dim] - comm->cell_x0[dim],
3860 ddbox->skew_fac[dim],
3861 dd->comm->cellsize_min[dim],
3862 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3866 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3868 /* Communicate the boundaries and update cell_ns_x0/1 */
3869 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3870 if (dd->bGridJump && dd->ndim > 1)
3872 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3877 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3881 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3889 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3890 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3899 static void check_screw_box(matrix box)
3901 /* Mathematical limitation */
3902 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3904 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3907 /* Limitation due to the asymmetry of the eighth shell method */
3908 if (box[ZZ][YY] != 0)
3910 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3914 static void distribute_cg(FILE *fplog,
3915 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3918 gmx_domdec_master_t *ma;
3919 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3920 int i, icg, j, k, k0, k1, d;
3924 real nrcg, inv_ncg, pos_d;
3930 if (tmp_ind == NULL)
3932 snew(tmp_nalloc, dd->nnodes);
3933 snew(tmp_ind, dd->nnodes);
3934 for (i = 0; i < dd->nnodes; i++)
3936 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3937 snew(tmp_ind[i], tmp_nalloc[i]);
3941 /* Clear the count */
3942 for (i = 0; i < dd->nnodes; i++)
3948 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3950 cgindex = cgs->index;
3952 /* Compute the center of geometry for all charge groups */
3953 for (icg = 0; icg < cgs->nr; icg++)
3956 k1 = cgindex[icg+1];
3960 copy_rvec(pos[k0], cg_cm);
3967 for (k = k0; (k < k1); k++)
3969 rvec_inc(cg_cm, pos[k]);
3971 for (d = 0; (d < DIM); d++)
3973 cg_cm[d] *= inv_ncg;
3976 /* Put the charge group in the box and determine the cell index */
3977 for (d = DIM-1; d >= 0; d--)
3980 if (d < dd->npbcdim)
3982 bScrew = (dd->bScrewPBC && d == XX);
3983 if (tric_dir[d] && dd->nc[d] > 1)
3985 /* Use triclinic coordintates for this dimension */
3986 for (j = d+1; j < DIM; j++)
3988 pos_d += cg_cm[j]*tcm[j][d];
3991 while (pos_d >= box[d][d])
3994 rvec_dec(cg_cm, box[d]);
3997 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3998 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4000 for (k = k0; (k < k1); k++)
4002 rvec_dec(pos[k], box[d]);
4005 pos[k][YY] = box[YY][YY] - pos[k][YY];
4006 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4013 rvec_inc(cg_cm, box[d]);
4016 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4017 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4019 for (k = k0; (k < k1); k++)
4021 rvec_inc(pos[k], box[d]);
4024 pos[k][YY] = box[YY][YY] - pos[k][YY];
4025 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4030 /* This could be done more efficiently */
4032 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4037 i = dd_index(dd->nc, ind);
4038 if (ma->ncg[i] == tmp_nalloc[i])
4040 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4041 srenew(tmp_ind[i], tmp_nalloc[i]);
4043 tmp_ind[i][ma->ncg[i]] = icg;
4045 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4049 for (i = 0; i < dd->nnodes; i++)
4052 for (k = 0; k < ma->ncg[i]; k++)
4054 ma->cg[k1++] = tmp_ind[i][k];
4057 ma->index[dd->nnodes] = k1;
4059 for (i = 0; i < dd->nnodes; i++)
4068 /* Here we avoid int overflows due to #atoms^2: use double, dsqr */
4069 int nat_sum, nat_min, nat_max;
4074 nat_min = ma->nat[0];
4075 nat_max = ma->nat[0];
4076 for (i = 0; i < dd->nnodes; i++)
4078 nat_sum += ma->nat[i];
4079 nat2_sum += dsqr(ma->nat[i]);
4080 nat_min = std::min(nat_min, ma->nat[i]);
4081 nat_max = std::max(nat_max, ma->nat[i]);
4083 nat_sum /= dd->nnodes;
4084 nat2_sum /= dd->nnodes;
4086 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
4089 static_cast<int>(sqrt(nat2_sum - dsqr(nat_sum) + 0.5)),
4094 static void get_cg_distribution(FILE *fplog, gmx_domdec_t *dd,
4095 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4098 gmx_domdec_master_t *ma = NULL;
4101 int *ibuf, buf2[2] = { 0, 0 };
4102 gmx_bool bMaster = DDMASTER(dd);
4110 check_screw_box(box);
4113 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4115 distribute_cg(fplog, box, ddbox->tric_dir, cgs, pos, dd);
4116 for (i = 0; i < dd->nnodes; i++)
4118 ma->ibuf[2*i] = ma->ncg[i];
4119 ma->ibuf[2*i+1] = ma->nat[i];
4127 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4129 dd->ncg_home = buf2[0];
4130 dd->nat_home = buf2[1];
4131 dd->ncg_tot = dd->ncg_home;
4132 dd->nat_tot = dd->nat_home;
4133 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4135 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4136 srenew(dd->index_gl, dd->cg_nalloc);
4137 srenew(dd->cgindex, dd->cg_nalloc+1);
4141 for (i = 0; i < dd->nnodes; i++)
4143 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4144 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4149 bMaster ? ma->ibuf : NULL,
4150 bMaster ? ma->ibuf+dd->nnodes : NULL,
4151 bMaster ? ma->cg : NULL,
4152 dd->ncg_home*sizeof(int), dd->index_gl);
4154 /* Determine the home charge group sizes */
4156 for (i = 0; i < dd->ncg_home; i++)
4158 cg_gl = dd->index_gl[i];
4160 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4165 fprintf(debug, "Home charge groups:\n");
4166 for (i = 0; i < dd->ncg_home; i++)
4168 fprintf(debug, " %d", dd->index_gl[i]);
4171 fprintf(debug, "\n");
4174 fprintf(debug, "\n");
4178 static int compact_and_copy_vec_at(int ncg, int *move,
4181 rvec *src, gmx_domdec_comm_t *comm,
4184 int m, icg, i, i0, i1, nrcg;
4190 for (m = 0; m < DIM*2; m++)
4196 for (icg = 0; icg < ncg; icg++)
4198 i1 = cgindex[icg+1];
4204 /* Compact the home array in place */
4205 for (i = i0; i < i1; i++)
4207 copy_rvec(src[i], src[home_pos++]);
4213 /* Copy to the communication buffer */
4215 pos_vec[m] += 1 + vec*nrcg;
4216 for (i = i0; i < i1; i++)
4218 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4220 pos_vec[m] += (nvec - vec - 1)*nrcg;
4224 home_pos += i1 - i0;
4232 static int compact_and_copy_vec_cg(int ncg, int *move,
4234 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4237 int m, icg, i0, i1, nrcg;
4243 for (m = 0; m < DIM*2; m++)
4249 for (icg = 0; icg < ncg; icg++)
4251 i1 = cgindex[icg+1];
4257 /* Compact the home array in place */
4258 copy_rvec(src[icg], src[home_pos++]);
4264 /* Copy to the communication buffer */
4265 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4266 pos_vec[m] += 1 + nrcg*nvec;
4278 static int compact_ind(int ncg, int *move,
4279 int *index_gl, int *cgindex,
4281 gmx_ga2la_t ga2la, char *bLocalCG,
4284 int cg, nat, a0, a1, a, a_gl;
4289 for (cg = 0; cg < ncg; cg++)
4295 /* Compact the home arrays in place.
4296 * Anything that can be done here avoids access to global arrays.
4298 cgindex[home_pos] = nat;
4299 for (a = a0; a < a1; a++)
4302 gatindex[nat] = a_gl;
4303 /* The cell number stays 0, so we don't need to set it */
4304 ga2la_change_la(ga2la, a_gl, nat);
4307 index_gl[home_pos] = index_gl[cg];
4308 cginfo[home_pos] = cginfo[cg];
4309 /* The charge group remains local, so bLocalCG does not change */
4314 /* Clear the global indices */
4315 for (a = a0; a < a1; a++)
4317 ga2la_del(ga2la, gatindex[a]);
4321 bLocalCG[index_gl[cg]] = FALSE;
4325 cgindex[home_pos] = nat;
4330 static void clear_and_mark_ind(int ncg, int *move,
4331 int *index_gl, int *cgindex, int *gatindex,
4332 gmx_ga2la_t ga2la, char *bLocalCG,
4337 for (cg = 0; cg < ncg; cg++)
4343 /* Clear the global indices */
4344 for (a = a0; a < a1; a++)
4346 ga2la_del(ga2la, gatindex[a]);
4350 bLocalCG[index_gl[cg]] = FALSE;
4352 /* Signal that this cg has moved using the ns cell index.
4353 * Here we set it to -1. fill_grid will change it
4354 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4356 cell_index[cg] = -1;
4361 static void print_cg_move(FILE *fplog,
4363 gmx_int64_t step, int cg, int dim, int dir,
4364 gmx_bool bHaveCgcmOld, real limitd,
4365 rvec cm_old, rvec cm_new, real pos_d)
4367 gmx_domdec_comm_t *comm;
4372 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4375 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4376 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4377 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4381 /* We don't have a limiting distance available: don't print it */
4382 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4383 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4384 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4386 fprintf(fplog, "distance out of cell %f\n",
4387 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4390 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4391 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4393 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4394 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4395 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4397 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4398 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4400 comm->cell_x0[dim], comm->cell_x1[dim]);
4403 static void cg_move_error(FILE *fplog,
4405 gmx_int64_t step, int cg, int dim, int dir,
4406 gmx_bool bHaveCgcmOld, real limitd,
4407 rvec cm_old, rvec cm_new, real pos_d)
4411 print_cg_move(fplog, dd, step, cg, dim, dir,
4412 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4414 print_cg_move(stderr, dd, step, cg, dim, dir,
4415 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4417 "%s moved too far between two domain decomposition steps\n"
4418 "This usually means that your system is not well equilibrated",
4419 dd->comm->bCGs ? "A charge group" : "An atom");
4422 static void rotate_state_atom(t_state *state, int a)
4426 for (est = 0; est < estNR; est++)
4428 if (EST_DISTR(est) && (state->flags & (1<<est)))
4433 /* Rotate the complete state; for a rectangular box only */
4434 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4435 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4438 state->v[a][YY] = -state->v[a][YY];
4439 state->v[a][ZZ] = -state->v[a][ZZ];
4442 state->sd_X[a][YY] = -state->sd_X[a][YY];
4443 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4446 state->cg_p[a][YY] = -state->cg_p[a][YY];
4447 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4449 case estDISRE_INITF:
4450 case estDISRE_RM3TAV:
4451 case estORIRE_INITF:
4453 /* These are distances, so not affected by rotation */
4456 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4462 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4464 if (natoms > comm->moved_nalloc)
4466 /* Contents should be preserved here */
4467 comm->moved_nalloc = over_alloc_dd(natoms);
4468 srenew(comm->moved, comm->moved_nalloc);
4474 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4477 ivec tric_dir, matrix tcm,
4478 rvec cell_x0, rvec cell_x1,
4479 rvec limitd, rvec limit0, rvec limit1,
4481 int cg_start, int cg_end,
4486 int cg, k, k0, k1, d, dim, d2;
4491 real inv_ncg, pos_d;
4494 npbcdim = dd->npbcdim;
4496 for (cg = cg_start; cg < cg_end; cg++)
4503 copy_rvec(state->x[k0], cm_new);
4510 for (k = k0; (k < k1); k++)
4512 rvec_inc(cm_new, state->x[k]);
4514 for (d = 0; (d < DIM); d++)
4516 cm_new[d] = inv_ncg*cm_new[d];
4521 /* Do pbc and check DD cell boundary crossings */
4522 for (d = DIM-1; d >= 0; d--)
4526 bScrew = (dd->bScrewPBC && d == XX);
4527 /* Determine the location of this cg in lattice coordinates */
4531 for (d2 = d+1; d2 < DIM; d2++)
4533 pos_d += cm_new[d2]*tcm[d2][d];
4536 /* Put the charge group in the triclinic unit-cell */
4537 if (pos_d >= cell_x1[d])
4539 if (pos_d >= limit1[d])
4541 cg_move_error(fplog, dd, step, cg, d, 1,
4542 cg_cm != state->x, limitd[d],
4543 cg_cm[cg], cm_new, pos_d);
4546 if (dd->ci[d] == dd->nc[d] - 1)
4548 rvec_dec(cm_new, state->box[d]);
4551 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4552 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4554 for (k = k0; (k < k1); k++)
4556 rvec_dec(state->x[k], state->box[d]);
4559 rotate_state_atom(state, k);
4564 else if (pos_d < cell_x0[d])
4566 if (pos_d < limit0[d])
4568 cg_move_error(fplog, dd, step, cg, d, -1,
4569 cg_cm != state->x, limitd[d],
4570 cg_cm[cg], cm_new, pos_d);
4575 rvec_inc(cm_new, state->box[d]);
4578 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4579 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4581 for (k = k0; (k < k1); k++)
4583 rvec_inc(state->x[k], state->box[d]);
4586 rotate_state_atom(state, k);
4592 else if (d < npbcdim)
4594 /* Put the charge group in the rectangular unit-cell */
4595 while (cm_new[d] >= state->box[d][d])
4597 rvec_dec(cm_new, state->box[d]);
4598 for (k = k0; (k < k1); k++)
4600 rvec_dec(state->x[k], state->box[d]);
4603 while (cm_new[d] < 0)
4605 rvec_inc(cm_new, state->box[d]);
4606 for (k = k0; (k < k1); k++)
4608 rvec_inc(state->x[k], state->box[d]);
4614 copy_rvec(cm_new, cg_cm[cg]);
4616 /* Determine where this cg should go */
4619 for (d = 0; d < dd->ndim; d++)
4624 flag |= DD_FLAG_FW(d);
4630 else if (dev[dim] == -1)
4632 flag |= DD_FLAG_BW(d);
4635 if (dd->nc[dim] > 2)
4646 /* Temporarily store the flag in move */
4647 move[cg] = mc + flag;
4651 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4652 gmx_domdec_t *dd, ivec tric_dir,
4653 t_state *state, rvec **f,
4662 int ncg[DIM*2], nat[DIM*2];
4663 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4664 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4665 int sbuf[2], rbuf[2];
4666 int home_pos_cg, home_pos_at, buf_pos;
4668 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4671 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4673 cginfo_mb_t *cginfo_mb;
4674 gmx_domdec_comm_t *comm;
4676 int nthread, thread;
4680 check_screw_box(state->box);
4684 if (fr->cutoff_scheme == ecutsGROUP)
4689 for (i = 0; i < estNR; i++)
4695 case estX: /* Always present */ break;
4696 case estV: bV = (state->flags & (1<<i)); break;
4697 case estSDX: bSDX = (state->flags & (1<<i)); break;
4698 case estCGP: bCGP = (state->flags & (1<<i)); break;
4701 case estDISRE_INITF:
4702 case estDISRE_RM3TAV:
4703 case estORIRE_INITF:
4705 /* No processing required */
4708 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4713 if (dd->ncg_tot > comm->nalloc_int)
4715 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4716 srenew(comm->buf_int, comm->nalloc_int);
4718 move = comm->buf_int;
4720 /* Clear the count */
4721 for (c = 0; c < dd->ndim*2; c++)
4727 npbcdim = dd->npbcdim;
4729 for (d = 0; (d < DIM); d++)
4731 limitd[d] = dd->comm->cellsize_min[d];
4732 if (d >= npbcdim && dd->ci[d] == 0)
4734 cell_x0[d] = -GMX_FLOAT_MAX;
4738 cell_x0[d] = comm->cell_x0[d];
4740 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4742 cell_x1[d] = GMX_FLOAT_MAX;
4746 cell_x1[d] = comm->cell_x1[d];
4750 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4751 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4755 /* We check after communication if a charge group moved
4756 * more than one cell. Set the pre-comm check limit to float_max.
4758 limit0[d] = -GMX_FLOAT_MAX;
4759 limit1[d] = GMX_FLOAT_MAX;
4763 make_tric_corr_matrix(npbcdim, state->box, tcm);
4765 cgindex = dd->cgindex;
4767 nthread = gmx_omp_nthreads_get(emntDomdec);
4769 /* Compute the center of geometry for all home charge groups
4770 * and put them in the box and determine where they should go.
4772 #pragma omp parallel for num_threads(nthread) schedule(static)
4773 for (thread = 0; thread < nthread; thread++)
4775 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4776 cell_x0, cell_x1, limitd, limit0, limit1,
4778 ( thread *dd->ncg_home)/nthread,
4779 ((thread+1)*dd->ncg_home)/nthread,
4780 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4784 for (cg = 0; cg < dd->ncg_home; cg++)
4789 flag = mc & ~DD_FLAG_NRCG;
4790 mc = mc & DD_FLAG_NRCG;
4793 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4795 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4796 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4798 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4799 /* We store the cg size in the lower 16 bits
4800 * and the place where the charge group should go
4801 * in the next 6 bits. This saves some communication volume.
4803 nrcg = cgindex[cg+1] - cgindex[cg];
4804 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4810 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4811 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4814 for (i = 0; i < dd->ndim*2; i++)
4816 *ncg_moved += ncg[i];
4833 /* Make sure the communication buffers are large enough */
4834 for (mc = 0; mc < dd->ndim*2; mc++)
4836 nvr = ncg[mc] + nat[mc]*nvec;
4837 if (nvr > comm->cgcm_state_nalloc[mc])
4839 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4840 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4844 switch (fr->cutoff_scheme)
4847 /* Recalculating cg_cm might be cheaper than communicating,
4848 * but that could give rise to rounding issues.
4851 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4852 nvec, cg_cm, comm, bCompact);
4855 /* Without charge groups we send the moved atom coordinates
4856 * over twice. This is so the code below can be used without
4857 * many conditionals for both for with and without charge groups.
4860 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4861 nvec, state->x, comm, FALSE);
4864 home_pos_cg -= *ncg_moved;
4868 gmx_incons("unimplemented");
4874 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4875 nvec, vec++, state->x, comm, bCompact);
4878 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4879 nvec, vec++, state->v, comm, bCompact);
4883 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4884 nvec, vec++, state->sd_X, comm, bCompact);
4888 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4889 nvec, vec++, state->cg_p, comm, bCompact);
4894 compact_ind(dd->ncg_home, move,
4895 dd->index_gl, dd->cgindex, dd->gatindex,
4896 dd->ga2la, comm->bLocalCG,
4901 if (fr->cutoff_scheme == ecutsVERLET)
4903 moved = get_moved(comm, dd->ncg_home);
4905 for (k = 0; k < dd->ncg_home; k++)
4912 moved = fr->ns.grid->cell_index;
4915 clear_and_mark_ind(dd->ncg_home, move,
4916 dd->index_gl, dd->cgindex, dd->gatindex,
4917 dd->ga2la, comm->bLocalCG,
4921 cginfo_mb = fr->cginfo_mb;
4923 *ncg_stay_home = home_pos_cg;
4924 for (d = 0; d < dd->ndim; d++)
4929 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4932 /* Communicate the cg and atom counts */
4937 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4938 d, dir, sbuf[0], sbuf[1]);
4940 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4942 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4944 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4945 srenew(comm->buf_int, comm->nalloc_int);
4948 /* Communicate the charge group indices, sizes and flags */
4949 dd_sendrecv_int(dd, d, dir,
4950 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4951 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4953 nvs = ncg[cdd] + nat[cdd]*nvec;
4954 i = rbuf[0] + rbuf[1] *nvec;
4955 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4957 /* Communicate cgcm and state */
4958 dd_sendrecv_rvec(dd, d, dir,
4959 comm->cgcm_state[cdd], nvs,
4960 comm->vbuf.v+nvr, i);
4961 ncg_recv += rbuf[0];
4965 /* Process the received charge groups */
4967 for (cg = 0; cg < ncg_recv; cg++)
4969 flag = comm->buf_int[cg*DD_CGIBS+1];
4971 if (dim >= npbcdim && dd->nc[dim] > 2)
4973 /* No pbc in this dim and more than one domain boundary.
4974 * We do a separate check if a charge group didn't move too far.
4976 if (((flag & DD_FLAG_FW(d)) &&
4977 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4978 ((flag & DD_FLAG_BW(d)) &&
4979 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4981 cg_move_error(fplog, dd, step, cg, dim,
4982 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4983 fr->cutoff_scheme == ecutsGROUP, 0,
4984 comm->vbuf.v[buf_pos],
4985 comm->vbuf.v[buf_pos],
4986 comm->vbuf.v[buf_pos][dim]);
4993 /* Check which direction this cg should go */
4994 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4998 /* The cell boundaries for dimension d2 are not equal
4999 * for each cell row of the lower dimension(s),
5000 * therefore we might need to redetermine where
5001 * this cg should go.
5004 /* If this cg crosses the box boundary in dimension d2
5005 * we can use the communicated flag, so we do not
5006 * have to worry about pbc.
5008 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5009 (flag & DD_FLAG_FW(d2))) ||
5010 (dd->ci[dim2] == 0 &&
5011 (flag & DD_FLAG_BW(d2)))))
5013 /* Clear the two flags for this dimension */
5014 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5015 /* Determine the location of this cg
5016 * in lattice coordinates
5018 pos_d = comm->vbuf.v[buf_pos][dim2];
5021 for (d3 = dim2+1; d3 < DIM; d3++)
5024 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5027 /* Check of we are not at the box edge.
5028 * pbc is only handled in the first step above,
5029 * but this check could move over pbc while
5030 * the first step did not due to different rounding.
5032 if (pos_d >= cell_x1[dim2] &&
5033 dd->ci[dim2] != dd->nc[dim2]-1)
5035 flag |= DD_FLAG_FW(d2);
5037 else if (pos_d < cell_x0[dim2] &&
5040 flag |= DD_FLAG_BW(d2);
5042 comm->buf_int[cg*DD_CGIBS+1] = flag;
5045 /* Set to which neighboring cell this cg should go */
5046 if (flag & DD_FLAG_FW(d2))
5050 else if (flag & DD_FLAG_BW(d2))
5052 if (dd->nc[dd->dim[d2]] > 2)
5064 nrcg = flag & DD_FLAG_NRCG;
5067 if (home_pos_cg+1 > dd->cg_nalloc)
5069 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5070 srenew(dd->index_gl, dd->cg_nalloc);
5071 srenew(dd->cgindex, dd->cg_nalloc+1);
5073 /* Set the global charge group index and size */
5074 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5075 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5076 /* Copy the state from the buffer */
5077 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5078 if (fr->cutoff_scheme == ecutsGROUP)
5081 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5085 /* Set the cginfo */
5086 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5087 dd->index_gl[home_pos_cg]);
5090 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5093 if (home_pos_at+nrcg > state->nalloc)
5095 dd_realloc_state(state, f, home_pos_at+nrcg);
5097 for (i = 0; i < nrcg; i++)
5099 copy_rvec(comm->vbuf.v[buf_pos++],
5100 state->x[home_pos_at+i]);
5104 for (i = 0; i < nrcg; i++)
5106 copy_rvec(comm->vbuf.v[buf_pos++],
5107 state->v[home_pos_at+i]);
5112 for (i = 0; i < nrcg; i++)
5114 copy_rvec(comm->vbuf.v[buf_pos++],
5115 state->sd_X[home_pos_at+i]);
5120 for (i = 0; i < nrcg; i++)
5122 copy_rvec(comm->vbuf.v[buf_pos++],
5123 state->cg_p[home_pos_at+i]);
5127 home_pos_at += nrcg;
5131 /* Reallocate the buffers if necessary */
5132 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5134 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5135 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5137 nvr = ncg[mc] + nat[mc]*nvec;
5138 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5140 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5141 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5143 /* Copy from the receive to the send buffers */
5144 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5145 comm->buf_int + cg*DD_CGIBS,
5146 DD_CGIBS*sizeof(int));
5147 memcpy(comm->cgcm_state[mc][nvr],
5148 comm->vbuf.v[buf_pos],
5149 (1+nrcg*nvec)*sizeof(rvec));
5150 buf_pos += 1 + nrcg*nvec;
5157 /* With sorting (!bCompact) the indices are now only partially up to date
5158 * and ncg_home and nat_home are not the real count, since there are
5159 * "holes" in the arrays for the charge groups that moved to neighbors.
5161 if (fr->cutoff_scheme == ecutsVERLET)
5163 moved = get_moved(comm, home_pos_cg);
5165 for (i = dd->ncg_home; i < home_pos_cg; i++)
5170 dd->ncg_home = home_pos_cg;
5171 dd->nat_home = home_pos_at;
5176 "Finished repartitioning: cgs moved out %d, new home %d\n",
5177 *ncg_moved, dd->ncg_home-*ncg_moved);
5182 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5184 dd->comm->cycl[ddCycl] += cycles;
5185 dd->comm->cycl_n[ddCycl]++;
5186 if (cycles > dd->comm->cycl_max[ddCycl])
5188 dd->comm->cycl_max[ddCycl] = cycles;
5192 static double force_flop_count(t_nrnb *nrnb)
5199 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5201 /* To get closer to the real timings, we half the count
5202 * for the normal loops and again half it for water loops.
5205 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5207 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5211 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5214 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5217 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5219 sum += nrnb->n[i]*cost_nrnb(i);
5222 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5224 sum += nrnb->n[i]*cost_nrnb(i);
5230 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5232 if (dd->comm->eFlop)
5234 dd->comm->flop -= force_flop_count(nrnb);
5237 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5239 if (dd->comm->eFlop)
5241 dd->comm->flop += force_flop_count(nrnb);
5246 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5250 for (i = 0; i < ddCyclNr; i++)
5252 dd->comm->cycl[i] = 0;
5253 dd->comm->cycl_n[i] = 0;
5254 dd->comm->cycl_max[i] = 0;
5257 dd->comm->flop_n = 0;
5260 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5262 gmx_domdec_comm_t *comm;
5263 gmx_domdec_load_t *load;
5264 gmx_domdec_root_t *root = NULL;
5266 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5271 fprintf(debug, "get_load_distribution start\n");
5274 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5278 bSepPME = (dd->pme_nodeid >= 0);
5280 for (d = dd->ndim-1; d >= 0; d--)
5283 /* Check if we participate in the communication in this dimension */
5284 if (d == dd->ndim-1 ||
5285 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5287 load = &comm->load[d];
5290 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5293 if (d == dd->ndim-1)
5295 sbuf[pos++] = dd_force_load(comm);
5296 sbuf[pos++] = sbuf[0];
5299 sbuf[pos++] = sbuf[0];
5300 sbuf[pos++] = cell_frac;
5303 sbuf[pos++] = comm->cell_f_max0[d];
5304 sbuf[pos++] = comm->cell_f_min1[d];
5309 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5310 sbuf[pos++] = comm->cycl[ddCyclPME];
5315 sbuf[pos++] = comm->load[d+1].sum;
5316 sbuf[pos++] = comm->load[d+1].max;
5319 sbuf[pos++] = comm->load[d+1].sum_m;
5320 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5321 sbuf[pos++] = comm->load[d+1].flags;
5324 sbuf[pos++] = comm->cell_f_max0[d];
5325 sbuf[pos++] = comm->cell_f_min1[d];
5330 sbuf[pos++] = comm->load[d+1].mdf;
5331 sbuf[pos++] = comm->load[d+1].pme;
5335 /* Communicate a row in DD direction d.
5336 * The communicators are setup such that the root always has rank 0.
5339 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5340 load->load, load->nload*sizeof(float), MPI_BYTE,
5341 0, comm->mpi_comm_load[d]);
5343 if (dd->ci[dim] == dd->master_ci[dim])
5345 /* We are the root, process this row */
5346 if (comm->bDynLoadBal)
5348 root = comm->root[d];
5358 for (i = 0; i < dd->nc[dim]; i++)
5360 load->sum += load->load[pos++];
5361 load->max = std::max(load->max, load->load[pos]);
5367 /* This direction could not be load balanced properly,
5368 * therefore we need to use the maximum iso the average load.
5370 load->sum_m = std::max(load->sum_m, load->load[pos]);
5374 load->sum_m += load->load[pos];
5377 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5381 load->flags = (int)(load->load[pos++] + 0.5);
5385 root->cell_f_max0[i] = load->load[pos++];
5386 root->cell_f_min1[i] = load->load[pos++];
5391 load->mdf = std::max(load->mdf, load->load[pos]);
5393 load->pme = std::max(load->pme, load->load[pos]);
5397 if (comm->bDynLoadBal && root->bLimited)
5399 load->sum_m *= dd->nc[dim];
5400 load->flags |= (1<<d);
5408 comm->nload += dd_load_count(comm);
5409 comm->load_step += comm->cycl[ddCyclStep];
5410 comm->load_sum += comm->load[0].sum;
5411 comm->load_max += comm->load[0].max;
5412 if (comm->bDynLoadBal)
5414 for (d = 0; d < dd->ndim; d++)
5416 if (comm->load[0].flags & (1<<d))
5418 comm->load_lim[d]++;
5424 comm->load_mdf += comm->load[0].mdf;
5425 comm->load_pme += comm->load[0].pme;
5429 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5433 fprintf(debug, "get_load_distribution finished\n");
5437 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5439 /* Return the relative performance loss on the total run time
5440 * due to the force calculation load imbalance.
5442 if (dd->comm->nload > 0)
5445 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5446 (dd->comm->load_step*dd->nnodes);
5454 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5457 int npp, npme, nnodes, d, limp;
5458 float imbal, pme_f_ratio, lossf, lossp = 0;
5460 gmx_domdec_comm_t *comm;
5463 if (DDMASTER(dd) && comm->nload > 0)
5466 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5467 nnodes = npp + npme;
5468 imbal = comm->load_max*npp/comm->load_sum - 1;
5469 lossf = dd_force_imb_perf_loss(dd);
5470 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5471 fprintf(fplog, "%s", buf);
5472 fprintf(stderr, "\n");
5473 fprintf(stderr, "%s", buf);
5474 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5475 fprintf(fplog, "%s", buf);
5476 fprintf(stderr, "%s", buf);
5478 if (comm->bDynLoadBal)
5480 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5481 for (d = 0; d < dd->ndim; d++)
5483 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5484 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5490 sprintf(buf+strlen(buf), "\n");
5491 fprintf(fplog, "%s", buf);
5492 fprintf(stderr, "%s", buf);
5496 pme_f_ratio = comm->load_pme/comm->load_mdf;
5497 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5500 lossp *= (float)npme/(float)nnodes;
5504 lossp *= (float)npp/(float)nnodes;
5506 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5507 fprintf(fplog, "%s", buf);
5508 fprintf(stderr, "%s", buf);
5509 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5510 fprintf(fplog, "%s", buf);
5511 fprintf(stderr, "%s", buf);
5513 fprintf(fplog, "\n");
5514 fprintf(stderr, "\n");
5516 if (lossf >= DD_PERF_LOSS_WARN)
5519 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5520 " in the domain decomposition.\n", lossf*100);
5521 if (!comm->bDynLoadBal)
5523 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5527 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5529 fprintf(fplog, "%s\n", buf);
5530 fprintf(stderr, "%s\n", buf);
5532 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5535 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5536 " had %s work to do than the PP ranks.\n"
5537 " You might want to %s the number of PME ranks\n"
5538 " or %s the cut-off and the grid spacing.\n",
5540 (lossp < 0) ? "less" : "more",
5541 (lossp < 0) ? "decrease" : "increase",
5542 (lossp < 0) ? "decrease" : "increase");
5543 fprintf(fplog, "%s\n", buf);
5544 fprintf(stderr, "%s\n", buf);
5549 static float dd_vol_min(gmx_domdec_t *dd)
5551 return dd->comm->load[0].cvol_min*dd->nnodes;
5554 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5556 return dd->comm->load[0].flags;
5559 static float dd_f_imbal(gmx_domdec_t *dd)
5561 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5564 float dd_pme_f_ratio(gmx_domdec_t *dd)
5566 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5568 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5576 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5581 flags = dd_load_flags(dd);
5585 "DD load balancing is limited by minimum cell size in dimension");
5586 for (d = 0; d < dd->ndim; d++)
5590 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5593 fprintf(fplog, "\n");
5595 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5596 if (dd->comm->bDynLoadBal)
5598 fprintf(fplog, " vol min/aver %5.3f%c",
5599 dd_vol_min(dd), flags ? '!' : ' ');
5601 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5602 if (dd->comm->cycl_n[ddCyclPME])
5604 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5606 fprintf(fplog, "\n\n");
5609 static void dd_print_load_verbose(gmx_domdec_t *dd)
5611 if (dd->comm->bDynLoadBal)
5613 fprintf(stderr, "vol %4.2f%c ",
5614 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5616 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5617 if (dd->comm->cycl_n[ddCyclPME])
5619 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5624 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5629 gmx_domdec_root_t *root;
5630 gmx_bool bPartOfGroup = FALSE;
5632 dim = dd->dim[dim_ind];
5633 copy_ivec(loc, loc_c);
5634 for (i = 0; i < dd->nc[dim]; i++)
5637 rank = dd_index(dd->nc, loc_c);
5638 if (rank == dd->rank)
5640 /* This process is part of the group */
5641 bPartOfGroup = TRUE;
5644 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5648 dd->comm->mpi_comm_load[dim_ind] = c_row;
5649 if (dd->comm->eDLB != edlbNO)
5651 if (dd->ci[dim] == dd->master_ci[dim])
5653 /* This is the root process of this row */
5654 snew(dd->comm->root[dim_ind], 1);
5655 root = dd->comm->root[dim_ind];
5656 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5657 snew(root->old_cell_f, dd->nc[dim]+1);
5658 snew(root->bCellMin, dd->nc[dim]);
5661 snew(root->cell_f_max0, dd->nc[dim]);
5662 snew(root->cell_f_min1, dd->nc[dim]);
5663 snew(root->bound_min, dd->nc[dim]);
5664 snew(root->bound_max, dd->nc[dim]);
5666 snew(root->buf_ncd, dd->nc[dim]);
5670 /* This is not a root process, we only need to receive cell_f */
5671 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5674 if (dd->ci[dim] == dd->master_ci[dim])
5676 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5682 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5683 const gmx_hw_info_t gmx_unused *hwinfo,
5684 const gmx_hw_opt_t gmx_unused *hw_opt)
5687 int physicalnode_id_hash;
5690 MPI_Comm mpi_comm_pp_physicalnode;
5692 if (!(cr->duty & DUTY_PP) ||
5693 hw_opt->gpu_opt.ncuda_dev_use == 0)
5695 /* Only PP nodes (currently) use GPUs.
5696 * If we don't have GPUs, there are no resources to share.
5701 physicalnode_id_hash = gmx_physicalnode_id_hash();
5703 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5709 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5710 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5711 dd->rank, physicalnode_id_hash, gpu_id);
5713 /* Split the PP communicator over the physical nodes */
5714 /* TODO: See if we should store this (before), as it's also used for
5715 * for the nodecomm summution.
5717 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5718 &mpi_comm_pp_physicalnode);
5719 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5720 &dd->comm->mpi_comm_gpu_shared);
5721 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5722 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5726 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5729 /* Note that some ranks could share a GPU, while others don't */
5731 if (dd->comm->nrank_gpu_shared == 1)
5733 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5738 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5741 int dim0, dim1, i, j;
5746 fprintf(debug, "Making load communicators\n");
5749 snew(dd->comm->load, dd->ndim);
5750 snew(dd->comm->mpi_comm_load, dd->ndim);
5753 make_load_communicator(dd, 0, loc);
5757 for (i = 0; i < dd->nc[dim0]; i++)
5760 make_load_communicator(dd, 1, loc);
5766 for (i = 0; i < dd->nc[dim0]; i++)
5770 for (j = 0; j < dd->nc[dim1]; j++)
5773 make_load_communicator(dd, 2, loc);
5780 fprintf(debug, "Finished making load communicators\n");
5785 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5787 int d, dim, i, j, m;
5790 ivec dd_zp[DD_MAXIZONE];
5791 gmx_domdec_zones_t *zones;
5792 gmx_domdec_ns_ranges_t *izone;
5794 for (d = 0; d < dd->ndim; d++)
5797 copy_ivec(dd->ci, tmp);
5798 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5799 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5800 copy_ivec(dd->ci, tmp);
5801 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5802 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5805 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5808 dd->neighbor[d][1]);
5814 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5816 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5817 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5824 for (i = 0; i < nzonep; i++)
5826 copy_ivec(dd_zp3[i], dd_zp[i]);
5832 for (i = 0; i < nzonep; i++)
5834 copy_ivec(dd_zp2[i], dd_zp[i]);
5840 for (i = 0; i < nzonep; i++)
5842 copy_ivec(dd_zp1[i], dd_zp[i]);
5846 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5851 zones = &dd->comm->zones;
5853 for (i = 0; i < nzone; i++)
5856 clear_ivec(zones->shift[i]);
5857 for (d = 0; d < dd->ndim; d++)
5859 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5864 for (i = 0; i < nzone; i++)
5866 for (d = 0; d < DIM; d++)
5868 s[d] = dd->ci[d] - zones->shift[i][d];
5873 else if (s[d] >= dd->nc[d])
5879 zones->nizone = nzonep;
5880 for (i = 0; i < zones->nizone; i++)
5882 if (dd_zp[i][0] != i)
5884 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5886 izone = &zones->izone[i];
5887 izone->j0 = dd_zp[i][1];
5888 izone->j1 = dd_zp[i][2];
5889 for (dim = 0; dim < DIM; dim++)
5891 if (dd->nc[dim] == 1)
5893 /* All shifts should be allowed */
5894 izone->shift0[dim] = -1;
5895 izone->shift1[dim] = 1;
5900 izone->shift0[d] = 0;
5901 izone->shift1[d] = 0;
5902 for(j=izone->j0; j<izone->j1; j++) {
5903 if (dd->shift[j][d] > dd->shift[i][d])
5904 izone->shift0[d] = -1;
5905 if (dd->shift[j][d] < dd->shift[i][d])
5906 izone->shift1[d] = 1;
5912 /* Assume the shift are not more than 1 cell */
5913 izone->shift0[dim] = 1;
5914 izone->shift1[dim] = -1;
5915 for (j = izone->j0; j < izone->j1; j++)
5917 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5918 if (shift_diff < izone->shift0[dim])
5920 izone->shift0[dim] = shift_diff;
5922 if (shift_diff > izone->shift1[dim])
5924 izone->shift1[dim] = shift_diff;
5931 if (dd->comm->eDLB != edlbNO)
5933 snew(dd->comm->root, dd->ndim);
5936 if (dd->comm->bRecordLoad)
5938 make_load_communicators(dd);
5942 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5948 gmx_domdec_comm_t *comm;
5955 if (comm->bCartesianPP)
5957 /* Set up cartesian communication for the particle-particle part */
5960 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5961 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5964 for (int i = 0; i < DIM; i++)
5968 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5970 /* We overwrite the old communicator with the new cartesian one */
5971 cr->mpi_comm_mygroup = comm_cart;
5974 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5975 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5977 if (comm->bCartesianPP_PME)
5979 /* Since we want to use the original cartesian setup for sim,
5980 * and not the one after split, we need to make an index.
5982 snew(comm->ddindex2ddnodeid, dd->nnodes);
5983 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5984 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5985 /* Get the rank of the DD master,
5986 * above we made sure that the master node is a PP node.
5996 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5998 else if (comm->bCartesianPP)
6000 if (cr->npmenodes == 0)
6002 /* The PP communicator is also
6003 * the communicator for this simulation
6005 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6007 cr->nodeid = dd->rank;
6009 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6011 /* We need to make an index to go from the coordinates
6012 * to the nodeid of this simulation.
6014 snew(comm->ddindex2simnodeid, dd->nnodes);
6015 snew(buf, dd->nnodes);
6016 if (cr->duty & DUTY_PP)
6018 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6020 /* Communicate the ddindex to simulation nodeid index */
6021 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6022 cr->mpi_comm_mysim);
6025 /* Determine the master coordinates and rank.
6026 * The DD master should be the same node as the master of this sim.
6028 for (int i = 0; i < dd->nnodes; i++)
6030 if (comm->ddindex2simnodeid[i] == 0)
6032 ddindex2xyz(dd->nc, i, dd->master_ci);
6033 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6038 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6043 /* No Cartesian communicators */
6044 /* We use the rank in dd->comm->all as DD index */
6045 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6046 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6048 clear_ivec(dd->master_ci);
6055 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6056 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6061 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6062 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6066 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6070 gmx_domdec_comm_t *comm;
6075 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6078 snew(comm->ddindex2simnodeid, dd->nnodes);
6079 snew(buf, dd->nnodes);
6080 if (cr->duty & DUTY_PP)
6082 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6084 /* Communicate the ddindex to simulation nodeid index */
6085 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6086 cr->mpi_comm_mysim);
6092 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6093 int ncg, int natoms)
6095 gmx_domdec_master_t *ma;
6100 snew(ma->ncg, dd->nnodes);
6101 snew(ma->index, dd->nnodes+1);
6103 snew(ma->nat, dd->nnodes);
6104 snew(ma->ibuf, dd->nnodes*2);
6105 snew(ma->cell_x, DIM);
6106 for (i = 0; i < DIM; i++)
6108 snew(ma->cell_x[i], dd->nc[i]+1);
6111 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6117 snew(ma->vbuf, natoms);
6123 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6124 int gmx_unused reorder)
6127 gmx_domdec_comm_t *comm;
6137 if (comm->bCartesianPP)
6139 for (i = 1; i < DIM; i++)
6141 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6143 if (bDiv[YY] || bDiv[ZZ])
6145 comm->bCartesianPP_PME = TRUE;
6146 /* If we have 2D PME decomposition, which is always in x+y,
6147 * we stack the PME only nodes in z.
6148 * Otherwise we choose the direction that provides the thinnest slab
6149 * of PME only nodes as this will have the least effect
6150 * on the PP communication.
6151 * But for the PME communication the opposite might be better.
6153 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6155 dd->nc[YY] > dd->nc[ZZ]))
6157 comm->cartpmedim = ZZ;
6161 comm->cartpmedim = YY;
6163 comm->ntot[comm->cartpmedim]
6164 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6168 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]);
6170 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6175 if (comm->bCartesianPP_PME)
6182 fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
6185 for (i = 0; i < DIM; i++)
6189 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6191 MPI_Comm_rank(comm_cart, &rank);
6192 if (MASTERNODE(cr) && rank != 0)
6194 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6197 /* With this assigment we loose the link to the original communicator
6198 * which will usually be MPI_COMM_WORLD, unless have multisim.
6200 cr->mpi_comm_mysim = comm_cart;
6201 cr->sim_nodeid = rank;
6203 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6207 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6208 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6211 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6215 if (cr->npmenodes == 0 ||
6216 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6218 cr->duty = DUTY_PME;
6221 /* Split the sim communicator into PP and PME only nodes */
6222 MPI_Comm_split(cr->mpi_comm_mysim,
6224 dd_index(comm->ntot, dd->ci),
6225 &cr->mpi_comm_mygroup);
6229 switch (dd_node_order)
6234 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6237 case ddnoINTERLEAVE:
6238 /* Interleave the PP-only and PME-only nodes,
6239 * as on clusters with dual-core machines this will double
6240 * the communication bandwidth of the PME processes
6241 * and thus speed up the PP <-> PME and inter PME communication.
6245 fprintf(fplog, "Interleaving PP and PME ranks\n");
6247 comm->pmenodes = dd_pmenodes(cr);
6252 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6255 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6257 cr->duty = DUTY_PME;
6264 /* Split the sim communicator into PP and PME only nodes */
6265 MPI_Comm_split(cr->mpi_comm_mysim,
6268 &cr->mpi_comm_mygroup);
6269 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6275 fprintf(fplog, "This rank does only %s work.\n\n",
6276 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6280 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6283 gmx_domdec_comm_t *comm;
6289 copy_ivec(dd->nc, comm->ntot);
6291 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6292 comm->bCartesianPP_PME = FALSE;
6294 /* Reorder the nodes by default. This might change the MPI ranks.
6295 * Real reordering is only supported on very few architectures,
6296 * Blue Gene is one of them.
6298 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6300 if (cr->npmenodes > 0)
6302 /* Split the communicator into a PP and PME part */
6303 split_communicator(fplog, cr, dd_node_order, CartReorder);
6304 if (comm->bCartesianPP_PME)
6306 /* We (possibly) reordered the nodes in split_communicator,
6307 * so it is no longer required in make_pp_communicator.
6309 CartReorder = FALSE;
6314 /* All nodes do PP and PME */
6316 /* We do not require separate communicators */
6317 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6321 if (cr->duty & DUTY_PP)
6323 /* Copy or make a new PP communicator */
6324 make_pp_communicator(fplog, cr, CartReorder);
6328 receive_ddindex2simnodeid(cr);
6331 if (!(cr->duty & DUTY_PME))
6333 /* Set up the commnuication to our PME node */
6334 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6335 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6338 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6339 dd->pme_nodeid, dd->pme_receive_vir_ener);
6344 dd->pme_nodeid = -1;
6349 dd->ma = init_gmx_domdec_master_t(dd,
6351 comm->cgs_gl.index[comm->cgs_gl.nr]);
6355 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6357 real *slb_frac, tot;
6362 if (nc > 1 && size_string != NULL)
6366 fprintf(fplog, "Using static load balancing for the %s direction\n",
6371 for (i = 0; i < nc; i++)
6374 sscanf(size_string, "%20lf%n", &dbl, &n);
6377 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6386 fprintf(fplog, "Relative cell sizes:");
6388 for (i = 0; i < nc; i++)
6393 fprintf(fplog, " %5.3f", slb_frac[i]);
6398 fprintf(fplog, "\n");
6405 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6408 gmx_mtop_ilistloop_t iloop;
6412 iloop = gmx_mtop_ilistloop_init(mtop);
6413 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6415 for (ftype = 0; ftype < F_NRE; ftype++)
6417 if ((interaction_function[ftype].flags & IF_BOND) &&
6420 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6428 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6434 val = getenv(env_var);
6437 if (sscanf(val, "%20d", &nst) <= 0)
6443 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6451 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6455 fprintf(stderr, "\n%s\n", warn_string);
6459 fprintf(fplog, "\n%s\n", warn_string);
6463 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6464 t_inputrec *ir, FILE *fplog)
6466 if (ir->ePBC == epbcSCREW &&
6467 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6469 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6472 if (ir->ns_type == ensSIMPLE)
6474 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6477 if (ir->nstlist == 0)
6479 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6482 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6484 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6488 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6493 r = ddbox->box_size[XX];
6494 for (di = 0; di < dd->ndim; di++)
6497 /* Check using the initial average cell size */
6498 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6504 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6505 const char *dlb_opt, gmx_bool bRecordLoad,
6506 unsigned long Flags, t_inputrec *ir)
6513 case 'a': eDLB = edlbAUTO; break;
6514 case 'n': eDLB = edlbNO; break;
6515 case 'y': eDLB = edlbYES; break;
6516 default: gmx_incons("Unknown dlb_opt");
6519 if (Flags & MD_RERUN)
6524 if (!EI_DYNAMICS(ir->eI))
6526 if (eDLB == edlbYES)
6528 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6529 dd_warning(cr, fplog, buf);
6537 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6542 if (Flags & MD_REPRODUCIBLE)
6549 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6553 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6556 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6564 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6569 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6571 /* Decomposition order z,y,x */
6574 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6576 for (dim = DIM-1; dim >= 0; dim--)
6578 if (dd->nc[dim] > 1)
6580 dd->dim[dd->ndim++] = dim;
6586 /* Decomposition order x,y,z */
6587 for (dim = 0; dim < DIM; dim++)
6589 if (dd->nc[dim] > 1)
6591 dd->dim[dd->ndim++] = dim;
6597 static gmx_domdec_comm_t *init_dd_comm()
6599 gmx_domdec_comm_t *comm;
6603 snew(comm->cggl_flag, DIM*2);
6604 snew(comm->cgcm_state, DIM*2);
6605 for (i = 0; i < DIM*2; i++)
6607 comm->cggl_flag_nalloc[i] = 0;
6608 comm->cgcm_state_nalloc[i] = 0;
6611 comm->nalloc_int = 0;
6612 comm->buf_int = NULL;
6614 vec_rvec_init(&comm->vbuf);
6616 comm->n_load_have = 0;
6617 comm->n_load_collect = 0;
6619 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6621 comm->sum_nat[i] = 0;
6625 comm->load_step = 0;
6628 clear_ivec(comm->load_lim);
6635 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6636 unsigned long Flags,
6638 real comm_distance_min, real rconstr,
6639 const char *dlb_opt, real dlb_scale,
6640 const char *sizex, const char *sizey, const char *sizez,
6641 gmx_mtop_t *mtop, t_inputrec *ir,
6642 matrix box, rvec *x,
6644 int *npme_x, int *npme_y)
6647 gmx_domdec_comm_t *comm;
6649 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6652 const real tenPercentMargin = 1.1;
6657 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6662 dd->comm = init_dd_comm();
6664 snew(comm->cggl_flag, DIM*2);
6665 snew(comm->cgcm_state, DIM*2);
6667 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6668 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6670 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6671 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6672 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6673 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6674 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6675 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6676 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6677 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6679 dd->pme_recv_f_alloc = 0;
6680 dd->pme_recv_f_buf = NULL;
6682 if (dd->bSendRecv2 && fplog)
6684 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");
6690 fprintf(fplog, "Will load balance based on FLOP count\n");
6692 if (comm->eFlop > 1)
6694 srand(1+cr->nodeid);
6696 comm->bRecordLoad = TRUE;
6700 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6704 /* Initialize to GPU share count to 0, might change later */
6705 comm->nrank_gpu_shared = 0;
6707 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6708 comm->bDLB_locked = FALSE;
6710 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6713 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6715 dd->bGridJump = comm->bDynLoadBal;
6716 comm->bPMELoadBalDLBLimits = FALSE;
6718 if (comm->nstSortCG)
6722 if (comm->nstSortCG == 1)
6724 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6728 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6732 snew(comm->sort, 1);
6738 fprintf(fplog, "Will not sort the charge groups\n");
6742 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6744 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6745 if (comm->bInterCGBondeds)
6747 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6751 comm->bInterCGMultiBody = FALSE;
6754 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6755 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6757 if (ir->rlistlong == 0)
6759 /* Set the cut-off to some very large value,
6760 * so we don't need if statements everywhere in the code.
6761 * We use sqrt, since the cut-off is squared in some places.
6763 comm->cutoff = GMX_CUTOFF_INF;
6767 comm->cutoff = ir->rlistlong;
6769 comm->cutoff_mbody = 0;
6771 comm->cellsize_limit = 0;
6772 comm->bBondComm = FALSE;
6774 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6775 * within nstlist steps. Since boundaries are allowed to displace by half
6776 * a cell size, DD cells should be at least the size of the list buffer.
6778 comm->cellsize_limit = std::max(comm->cellsize_limit,
6779 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6781 if (comm->bInterCGBondeds)
6783 if (comm_distance_min > 0)
6785 comm->cutoff_mbody = comm_distance_min;
6786 if (Flags & MD_DDBONDCOMM)
6788 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6792 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6794 r_bonded_limit = comm->cutoff_mbody;
6796 else if (ir->bPeriodicMols)
6798 /* Can not easily determine the required cut-off */
6799 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");
6800 comm->cutoff_mbody = comm->cutoff/2;
6801 r_bonded_limit = comm->cutoff_mbody;
6807 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6808 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6810 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6811 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6813 /* We use an initial margin of 10% for the minimum cell size,
6814 * except when we are just below the non-bonded cut-off.
6816 if (Flags & MD_DDBONDCOMM)
6818 if (std::max(r_2b, r_mb) > comm->cutoff)
6820 r_bonded = std::max(r_2b, r_mb);
6821 r_bonded_limit = tenPercentMargin*r_bonded;
6822 comm->bBondComm = TRUE;
6827 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6829 /* We determine cutoff_mbody later */
6833 /* No special bonded communication,
6834 * simply increase the DD cut-off.
6836 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6837 comm->cutoff_mbody = r_bonded_limit;
6838 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6844 "Minimum cell size due to bonded interactions: %.3f nm\n",
6847 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6850 if (dd->bInterCGcons && rconstr <= 0)
6852 /* There is a cell size limit due to the constraints (P-LINCS) */
6853 rconstr = constr_r_max(fplog, mtop, ir);
6857 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6859 if (rconstr > comm->cellsize_limit)
6861 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6865 else if (rconstr > 0 && fplog)
6867 /* Here we do not check for dd->bInterCGcons,
6868 * because one can also set a cell size limit for virtual sites only
6869 * and at this point we don't know yet if there are intercg v-sites.
6872 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6875 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6877 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6881 copy_ivec(nc, dd->nc);
6882 set_dd_dim(fplog, dd);
6883 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6885 if (cr->npmenodes == -1)
6889 acs = average_cellsize_min(dd, ddbox);
6890 if (acs < comm->cellsize_limit)
6894 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6896 gmx_fatal_collective(FARGS, cr, NULL,
6897 "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",
6898 acs, comm->cellsize_limit);
6903 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6905 /* We need to choose the optimal DD grid and possibly PME nodes */
6906 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6907 comm->eDLB != edlbNO, dlb_scale,
6908 comm->cellsize_limit, comm->cutoff,
6909 comm->bInterCGBondeds);
6911 if (dd->nc[XX] == 0)
6913 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6914 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6915 !bC ? "-rdd" : "-rcon",
6916 comm->eDLB != edlbNO ? " or -dds" : "",
6917 bC ? " or your LINCS settings" : "");
6919 gmx_fatal_collective(FARGS, cr, NULL,
6920 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6922 "Look in the log file for details on the domain decomposition",
6923 cr->nnodes-cr->npmenodes, limit, buf);
6925 set_dd_dim(fplog, dd);
6931 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6932 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6935 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6936 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6938 gmx_fatal_collective(FARGS, cr, NULL,
6939 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6940 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6942 if (cr->npmenodes > dd->nnodes)
6944 gmx_fatal_collective(FARGS, cr, NULL,
6945 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6947 if (cr->npmenodes > 0)
6949 comm->npmenodes = cr->npmenodes;
6953 comm->npmenodes = dd->nnodes;
6956 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6958 /* The following choices should match those
6959 * in comm_cost_est in domdec_setup.c.
6960 * Note that here the checks have to take into account
6961 * that the decomposition might occur in a different order than xyz
6962 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6963 * in which case they will not match those in comm_cost_est,
6964 * but since that is mainly for testing purposes that's fine.
6966 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6967 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6968 getenv("GMX_PMEONEDD") == NULL)
6970 comm->npmedecompdim = 2;
6971 comm->npmenodes_x = dd->nc[XX];
6972 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6976 /* In case nc is 1 in both x and y we could still choose to
6977 * decompose pme in y instead of x, but we use x for simplicity.
6979 comm->npmedecompdim = 1;
6980 if (dd->dim[0] == YY)
6982 comm->npmenodes_x = 1;
6983 comm->npmenodes_y = comm->npmenodes;
6987 comm->npmenodes_x = comm->npmenodes;
6988 comm->npmenodes_y = 1;
6993 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6994 comm->npmenodes_x, comm->npmenodes_y, 1);
6999 comm->npmedecompdim = 0;
7000 comm->npmenodes_x = 0;
7001 comm->npmenodes_y = 0;
7004 /* Technically we don't need both of these,
7005 * but it simplifies code not having to recalculate it.
7007 *npme_x = comm->npmenodes_x;
7008 *npme_y = comm->npmenodes_y;
7010 snew(comm->slb_frac, DIM);
7011 if (comm->eDLB == edlbNO)
7013 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7014 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7015 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7018 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7020 if (comm->bBondComm || comm->eDLB != edlbNO)
7022 /* Set the bonded communication distance to halfway
7023 * the minimum and the maximum,
7024 * since the extra communication cost is nearly zero.
7026 acs = average_cellsize_min(dd, ddbox);
7027 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7028 if (comm->eDLB != edlbNO)
7030 /* Check if this does not limit the scaling */
7031 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7033 if (!comm->bBondComm)
7035 /* Without bBondComm do not go beyond the n.b. cut-off */
7036 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7037 if (comm->cellsize_limit >= comm->cutoff)
7039 /* We don't loose a lot of efficieny
7040 * when increasing it to the n.b. cut-off.
7041 * It can even be slightly faster, because we need
7042 * less checks for the communication setup.
7044 comm->cutoff_mbody = comm->cutoff;
7047 /* Check if we did not end up below our original limit */
7048 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7050 if (comm->cutoff_mbody > comm->cellsize_limit)
7052 comm->cellsize_limit = comm->cutoff_mbody;
7055 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7060 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7061 "cellsize limit %f\n",
7062 comm->bBondComm, comm->cellsize_limit);
7067 check_dd_restrictions(cr, dd, ir, fplog);
7070 comm->partition_step = INT_MIN;
7073 clear_dd_cycle_counts(dd);
7078 static void set_dlb_limits(gmx_domdec_t *dd)
7083 for (d = 0; d < dd->ndim; d++)
7085 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7086 dd->comm->cellsize_min[dd->dim[d]] =
7087 dd->comm->cellsize_min_dlb[dd->dim[d]];
7092 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7095 gmx_domdec_comm_t *comm;
7105 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);
7108 cellsize_min = comm->cellsize_min[dd->dim[0]];
7109 for (d = 1; d < dd->ndim; d++)
7111 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7114 if (cellsize_min < comm->cellsize_limit*1.05)
7116 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");
7118 /* Change DLB from "auto" to "no". */
7119 comm->eDLB = edlbNO;
7124 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7125 comm->bDynLoadBal = TRUE;
7126 dd->bGridJump = TRUE;
7130 /* We can set the required cell size info here,
7131 * so we do not need to communicate this.
7132 * The grid is completely uniform.
7134 for (d = 0; d < dd->ndim; d++)
7138 comm->load[d].sum_m = comm->load[d].sum;
7140 nc = dd->nc[dd->dim[d]];
7141 for (i = 0; i < nc; i++)
7143 comm->root[d]->cell_f[i] = i/(real)nc;
7146 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7147 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7150 comm->root[d]->cell_f[nc] = 1.0;
7155 static char *init_bLocalCG(gmx_mtop_t *mtop)
7160 ncg = ncg_mtop(mtop);
7161 snew(bLocalCG, ncg);
7162 for (cg = 0; cg < ncg; cg++)
7164 bLocalCG[cg] = FALSE;
7170 void dd_init_bondeds(FILE *fplog,
7171 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7173 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7175 gmx_domdec_comm_t *comm;
7177 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7181 if (comm->bBondComm)
7183 /* Communicate atoms beyond the cut-off for bonded interactions */
7186 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7188 comm->bLocalCG = init_bLocalCG(mtop);
7192 /* Only communicate atoms based on cut-off */
7193 comm->cglink = NULL;
7194 comm->bLocalCG = NULL;
7198 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7200 gmx_bool bDynLoadBal, real dlb_scale,
7203 gmx_domdec_comm_t *comm;
7218 fprintf(fplog, "The maximum number of communication pulses is:");
7219 for (d = 0; d < dd->ndim; d++)
7221 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7223 fprintf(fplog, "\n");
7224 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7225 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7226 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7227 for (d = 0; d < DIM; d++)
7231 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7238 comm->cellsize_min_dlb[d]/
7239 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7241 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7244 fprintf(fplog, "\n");
7248 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7249 fprintf(fplog, "The initial number of communication pulses is:");
7250 for (d = 0; d < dd->ndim; d++)
7252 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7254 fprintf(fplog, "\n");
7255 fprintf(fplog, "The initial domain decomposition cell size is:");
7256 for (d = 0; d < DIM; d++)
7260 fprintf(fplog, " %c %.2f nm",
7261 dim2char(d), dd->comm->cellsize_min[d]);
7264 fprintf(fplog, "\n\n");
7267 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7269 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7270 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7271 "non-bonded interactions", "", comm->cutoff);
7275 limit = dd->comm->cellsize_limit;
7279 if (dynamic_dd_box(ddbox, ir))
7281 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7283 limit = dd->comm->cellsize_min[XX];
7284 for (d = 1; d < DIM; d++)
7286 limit = std::min(limit, dd->comm->cellsize_min[d]);
7290 if (comm->bInterCGBondeds)
7292 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7293 "two-body bonded interactions", "(-rdd)",
7294 std::max(comm->cutoff, comm->cutoff_mbody));
7295 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7296 "multi-body bonded interactions", "(-rdd)",
7297 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7301 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7302 "virtual site constructions", "(-rcon)", limit);
7304 if (dd->constraint_comm)
7306 sprintf(buf, "atoms separated by up to %d constraints",
7308 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7309 buf, "(-rcon)", limit);
7311 fprintf(fplog, "\n");
7317 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7319 const t_inputrec *ir,
7320 const gmx_ddbox_t *ddbox)
7322 gmx_domdec_comm_t *comm;
7323 int d, dim, npulse, npulse_d_max, npulse_d;
7328 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7330 /* Determine the maximum number of comm. pulses in one dimension */
7332 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7334 /* Determine the maximum required number of grid pulses */
7335 if (comm->cellsize_limit >= comm->cutoff)
7337 /* Only a single pulse is required */
7340 else if (!bNoCutOff && comm->cellsize_limit > 0)
7342 /* We round down slightly here to avoid overhead due to the latency
7343 * of extra communication calls when the cut-off
7344 * would be only slightly longer than the cell size.
7345 * Later cellsize_limit is redetermined,
7346 * so we can not miss interactions due to this rounding.
7348 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7352 /* There is no cell size limit */
7353 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7356 if (!bNoCutOff && npulse > 1)
7358 /* See if we can do with less pulses, based on dlb_scale */
7360 for (d = 0; d < dd->ndim; d++)
7363 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7364 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7365 npulse_d_max = std::max(npulse_d_max, npulse_d);
7367 npulse = std::min(npulse, npulse_d_max);
7370 /* This env var can override npulse */
7371 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7378 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7379 for (d = 0; d < dd->ndim; d++)
7381 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7382 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7383 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7384 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7385 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7387 comm->bVacDLBNoLimit = FALSE;
7391 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7392 if (!comm->bVacDLBNoLimit)
7394 comm->cellsize_limit = std::max(comm->cellsize_limit,
7395 comm->cutoff/comm->maxpulse);
7397 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7398 /* Set the minimum cell size for each DD dimension */
7399 for (d = 0; d < dd->ndim; d++)
7401 if (comm->bVacDLBNoLimit ||
7402 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7404 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7408 comm->cellsize_min_dlb[dd->dim[d]] =
7409 comm->cutoff/comm->cd[d].np_dlb;
7412 if (comm->cutoff_mbody <= 0)
7414 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7416 if (comm->bDynLoadBal)
7422 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7424 /* If each molecule is a single charge group
7425 * or we use domain decomposition for each periodic dimension,
7426 * we do not need to take pbc into account for the bonded interactions.
7428 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7431 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7434 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7435 t_inputrec *ir, gmx_ddbox_t *ddbox)
7437 gmx_domdec_comm_t *comm;
7443 /* Initialize the thread data.
7444 * This can not be done in init_domain_decomposition,
7445 * as the numbers of threads is determined later.
7447 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7450 snew(comm->dth, comm->nth);
7453 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7455 init_ddpme(dd, &comm->ddpme[0], 0);
7456 if (comm->npmedecompdim >= 2)
7458 init_ddpme(dd, &comm->ddpme[1], 1);
7463 comm->npmenodes = 0;
7464 if (dd->pme_nodeid >= 0)
7466 gmx_fatal_collective(FARGS, NULL, dd,
7467 "Can not have separate PME ranks without PME electrostatics");
7473 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7475 if (comm->eDLB != edlbNO)
7477 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7480 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7481 if (comm->eDLB == edlbAUTO)
7485 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7487 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7490 if (ir->ePBC == epbcNONE)
7492 vol_frac = 1 - 1/(double)dd->nnodes;
7497 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7501 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7503 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7505 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7508 static gmx_bool test_dd_cutoff(t_commrec *cr,
7509 t_state *state, t_inputrec *ir,
7520 set_ddbox(dd, FALSE, cr, ir, state->box,
7521 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7525 for (d = 0; d < dd->ndim; d++)
7529 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7530 if (dynamic_dd_box(&ddbox, ir))
7532 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7535 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7537 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7538 dd->comm->cd[d].np_dlb > 0)
7540 if (np > dd->comm->cd[d].np_dlb)
7545 /* If a current local cell size is smaller than the requested
7546 * cut-off, we could still fix it, but this gets very complicated.
7547 * Without fixing here, we might actually need more checks.
7549 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7556 if (dd->comm->eDLB != edlbNO)
7558 /* If DLB is not active yet, we don't need to check the grid jumps.
7559 * Actually we shouldn't, because then the grid jump data is not set.
7561 if (dd->comm->bDynLoadBal &&
7562 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7567 gmx_sumi(1, &LocallyLimited, cr);
7569 if (LocallyLimited > 0)
7578 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7581 gmx_bool bCutoffAllowed;
7583 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7587 cr->dd->comm->cutoff = cutoff_req;
7590 return bCutoffAllowed;
7593 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7595 gmx_domdec_comm_t *comm;
7597 comm = cr->dd->comm;
7599 /* Turn on the DLB limiting (might have been on already) */
7600 comm->bPMELoadBalDLBLimits = TRUE;
7602 /* Change the cut-off limit */
7603 comm->PMELoadBal_max_cutoff = comm->cutoff;
7606 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7608 return dd->comm->bDLB_locked;
7611 void dd_dlb_set_lock(gmx_domdec_t *dd, gmx_bool bValue)
7613 /* We can only lock the DLB when it is set to auto, otherwise don't lock */
7614 if (dd->comm->eDLB == edlbAUTO)
7616 dd->comm->bDLB_locked = bValue;
7620 static void merge_cg_buffers(int ncell,
7621 gmx_domdec_comm_dim_t *cd, int pulse,
7623 int *index_gl, int *recv_i,
7624 rvec *cg_cm, rvec *recv_vr,
7626 cginfo_mb_t *cginfo_mb, int *cginfo)
7628 gmx_domdec_ind_t *ind, *ind_p;
7629 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7630 int shift, shift_at;
7632 ind = &cd->ind[pulse];
7634 /* First correct the already stored data */
7635 shift = ind->nrecv[ncell];
7636 for (cell = ncell-1; cell >= 0; cell--)
7638 shift -= ind->nrecv[cell];
7641 /* Move the cg's present from previous grid pulses */
7642 cg0 = ncg_cell[ncell+cell];
7643 cg1 = ncg_cell[ncell+cell+1];
7644 cgindex[cg1+shift] = cgindex[cg1];
7645 for (cg = cg1-1; cg >= cg0; cg--)
7647 index_gl[cg+shift] = index_gl[cg];
7648 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7649 cgindex[cg+shift] = cgindex[cg];
7650 cginfo[cg+shift] = cginfo[cg];
7652 /* Correct the already stored send indices for the shift */
7653 for (p = 1; p <= pulse; p++)
7655 ind_p = &cd->ind[p];
7657 for (c = 0; c < cell; c++)
7659 cg0 += ind_p->nsend[c];
7661 cg1 = cg0 + ind_p->nsend[cell];
7662 for (cg = cg0; cg < cg1; cg++)
7664 ind_p->index[cg] += shift;
7670 /* Merge in the communicated buffers */
7674 for (cell = 0; cell < ncell; cell++)
7676 cg1 = ncg_cell[ncell+cell+1] + shift;
7679 /* Correct the old cg indices */
7680 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7682 cgindex[cg+1] += shift_at;
7685 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7687 /* Copy this charge group from the buffer */
7688 index_gl[cg1] = recv_i[cg0];
7689 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7690 /* Add it to the cgindex */
7691 cg_gl = index_gl[cg1];
7692 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7693 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7694 cgindex[cg1+1] = cgindex[cg1] + nat;
7699 shift += ind->nrecv[cell];
7700 ncg_cell[ncell+cell+1] = cg1;
7704 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7705 int nzone, int cg0, const int *cgindex)
7709 /* Store the atom block boundaries for easy copying of communication buffers
7712 for (zone = 0; zone < nzone; zone++)
7714 for (p = 0; p < cd->np; p++)
7716 cd->ind[p].cell2at0[zone] = cgindex[cg];
7717 cg += cd->ind[p].nrecv[zone];
7718 cd->ind[p].cell2at1[zone] = cgindex[cg];
7723 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7729 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7731 if (!bLocalCG[link->a[i]])
7740 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7742 real c[DIM][4]; /* the corners for the non-bonded communication */
7743 real cr0; /* corner for rounding */
7744 real cr1[4]; /* corners for rounding */
7745 real bc[DIM]; /* corners for bounded communication */
7746 real bcr1; /* corner for rounding for bonded communication */
7749 /* Determine the corners of the domain(s) we are communicating with */
7751 set_dd_corners(const gmx_domdec_t *dd,
7752 int dim0, int dim1, int dim2,
7756 const gmx_domdec_comm_t *comm;
7757 const gmx_domdec_zones_t *zones;
7762 zones = &comm->zones;
7764 /* Keep the compiler happy */
7768 /* The first dimension is equal for all cells */
7769 c->c[0][0] = comm->cell_x0[dim0];
7772 c->bc[0] = c->c[0][0];
7777 /* This cell row is only seen from the first row */
7778 c->c[1][0] = comm->cell_x0[dim1];
7779 /* All rows can see this row */
7780 c->c[1][1] = comm->cell_x0[dim1];
7783 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7786 /* For the multi-body distance we need the maximum */
7787 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7790 /* Set the upper-right corner for rounding */
7791 c->cr0 = comm->cell_x1[dim0];
7796 for (j = 0; j < 4; j++)
7798 c->c[2][j] = comm->cell_x0[dim2];
7802 /* Use the maximum of the i-cells that see a j-cell */
7803 for (i = 0; i < zones->nizone; i++)
7805 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7810 std::max(c->c[2][j-4],
7811 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7817 /* For the multi-body distance we need the maximum */
7818 c->bc[2] = comm->cell_x0[dim2];
7819 for (i = 0; i < 2; i++)
7821 for (j = 0; j < 2; j++)
7823 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7829 /* Set the upper-right corner for rounding */
7830 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7831 * Only cell (0,0,0) can see cell 7 (1,1,1)
7833 c->cr1[0] = comm->cell_x1[dim1];
7834 c->cr1[3] = comm->cell_x1[dim1];
7837 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7840 /* For the multi-body distance we need the maximum */
7841 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7848 /* Determine which cg's we need to send in this pulse from this zone */
7850 get_zone_pulse_cgs(gmx_domdec_t *dd,
7851 int zonei, int zone,
7853 const int *index_gl,
7855 int dim, int dim_ind,
7856 int dim0, int dim1, int dim2,
7857 real r_comm2, real r_bcomm2,
7861 real skew_fac2_d, real skew_fac_01,
7862 rvec *v_d, rvec *v_0, rvec *v_1,
7863 const dd_corners_t *c,
7865 gmx_bool bDistBonded,
7871 gmx_domdec_ind_t *ind,
7872 int **ibuf, int *ibuf_nalloc,
7878 gmx_domdec_comm_t *comm;
7880 gmx_bool bDistMB_pulse;
7882 real r2, rb2, r, tric_sh;
7885 int nsend_z, nsend, nat;
7889 bScrew = (dd->bScrewPBC && dim == XX);
7891 bDistMB_pulse = (bDistMB && bDistBonded);
7897 for (cg = cg0; cg < cg1; cg++)
7901 if (tric_dist[dim_ind] == 0)
7903 /* Rectangular direction, easy */
7904 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7911 r = cg_cm[cg][dim] - c->bc[dim_ind];
7917 /* Rounding gives at most a 16% reduction
7918 * in communicated atoms
7920 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7922 r = cg_cm[cg][dim0] - c->cr0;
7923 /* This is the first dimension, so always r >= 0 */
7930 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7932 r = cg_cm[cg][dim1] - c->cr1[zone];
7939 r = cg_cm[cg][dim1] - c->bcr1;
7949 /* Triclinic direction, more complicated */
7952 /* Rounding, conservative as the skew_fac multiplication
7953 * will slightly underestimate the distance.
7955 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7957 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7958 for (i = dim0+1; i < DIM; i++)
7960 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7962 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7965 rb[dim0] = rn[dim0];
7968 /* Take care that the cell planes along dim0 might not
7969 * be orthogonal to those along dim1 and dim2.
7971 for (i = 1; i <= dim_ind; i++)
7974 if (normal[dim0][dimd] > 0)
7976 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7979 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7984 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7986 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7988 for (i = dim1+1; i < DIM; i++)
7990 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7992 rn[dim1] += tric_sh;
7995 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7996 /* Take care of coupling of the distances
7997 * to the planes along dim0 and dim1 through dim2.
7999 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
8000 /* Take care that the cell planes along dim1
8001 * might not be orthogonal to that along dim2.
8003 if (normal[dim1][dim2] > 0)
8005 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8011 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8014 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8015 /* Take care of coupling of the distances
8016 * to the planes along dim0 and dim1 through dim2.
8018 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8019 /* Take care that the cell planes along dim1
8020 * might not be orthogonal to that along dim2.
8022 if (normal[dim1][dim2] > 0)
8024 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8029 /* The distance along the communication direction */
8030 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8032 for (i = dim+1; i < DIM; i++)
8034 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8039 r2 += rn[dim]*rn[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 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8051 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8054 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8055 /* Take care of coupling of the distances
8056 * to the planes along dim0 and dim1 through dim2.
8058 if (dim_ind == 1 && zonei == 1)
8060 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8068 ((bDistMB && rb2 < r_bcomm2) ||
8069 (bDist2B && r2 < r_bcomm2)) &&
8071 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8072 missing_link(comm->cglink, index_gl[cg],
8075 /* Make an index to the local charge groups */
8076 if (nsend+1 > ind->nalloc)
8078 ind->nalloc = over_alloc_large(nsend+1);
8079 srenew(ind->index, ind->nalloc);
8081 if (nsend+1 > *ibuf_nalloc)
8083 *ibuf_nalloc = over_alloc_large(nsend+1);
8084 srenew(*ibuf, *ibuf_nalloc);
8086 ind->index[nsend] = cg;
8087 (*ibuf)[nsend] = index_gl[cg];
8089 vec_rvec_check_alloc(vbuf, nsend+1);
8091 if (dd->ci[dim] == 0)
8093 /* Correct cg_cm for pbc */
8094 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8097 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8098 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8103 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8106 nat += cgindex[cg+1] - cgindex[cg];
8112 *nsend_z_ptr = nsend_z;
8115 static void setup_dd_communication(gmx_domdec_t *dd,
8116 matrix box, gmx_ddbox_t *ddbox,
8117 t_forcerec *fr, t_state *state, rvec **f)
8119 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8120 int nzone, nzone_send, zone, zonei, cg0, cg1;
8121 int c, i, cg, cg_gl, nrcg;
8122 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8123 gmx_domdec_comm_t *comm;
8124 gmx_domdec_zones_t *zones;
8125 gmx_domdec_comm_dim_t *cd;
8126 gmx_domdec_ind_t *ind;
8127 cginfo_mb_t *cginfo_mb;
8128 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8129 real r_comm2, r_bcomm2;
8130 dd_corners_t corners;
8132 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8133 real skew_fac2_d, skew_fac_01;
8140 fprintf(debug, "Setting up DD communication\n");
8145 switch (fr->cutoff_scheme)
8154 gmx_incons("unimplemented");
8158 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8160 /* Check if we need to use triclinic distances */
8161 tric_dist[dim_ind] = 0;
8162 for (i = 0; i <= dim_ind; i++)
8164 if (ddbox->tric_dir[dd->dim[i]])
8166 tric_dist[dim_ind] = 1;
8171 bBondComm = comm->bBondComm;
8173 /* Do we need to determine extra distances for multi-body bondeds? */
8174 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8176 /* Do we need to determine extra distances for only two-body bondeds? */
8177 bDist2B = (bBondComm && !bDistMB);
8179 r_comm2 = sqr(comm->cutoff);
8180 r_bcomm2 = sqr(comm->cutoff_mbody);
8184 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8187 zones = &comm->zones;
8190 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8191 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8193 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8195 /* Triclinic stuff */
8196 normal = ddbox->normal;
8200 v_0 = ddbox->v[dim0];
8201 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8203 /* Determine the coupling coefficient for the distances
8204 * to the cell planes along dim0 and dim1 through dim2.
8205 * This is required for correct rounding.
8208 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8211 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8217 v_1 = ddbox->v[dim1];
8220 zone_cg_range = zones->cg_range;
8221 index_gl = dd->index_gl;
8222 cgindex = dd->cgindex;
8223 cginfo_mb = fr->cginfo_mb;
8225 zone_cg_range[0] = 0;
8226 zone_cg_range[1] = dd->ncg_home;
8227 comm->zone_ncg1[0] = dd->ncg_home;
8228 pos_cg = dd->ncg_home;
8230 nat_tot = dd->nat_home;
8232 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8234 dim = dd->dim[dim_ind];
8235 cd = &comm->cd[dim_ind];
8237 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8239 /* No pbc in this dimension, the first node should not comm. */
8247 v_d = ddbox->v[dim];
8248 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8250 cd->bInPlace = TRUE;
8251 for (p = 0; p < cd->np; p++)
8253 /* Only atoms communicated in the first pulse are used
8254 * for multi-body bonded interactions or for bBondComm.
8256 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8261 for (zone = 0; zone < nzone_send; zone++)
8263 if (tric_dist[dim_ind] && dim_ind > 0)
8265 /* Determine slightly more optimized skew_fac's
8267 * This reduces the number of communicated atoms
8268 * by about 10% for 3D DD of rhombic dodecahedra.
8270 for (dimd = 0; dimd < dim; dimd++)
8272 sf2_round[dimd] = 1;
8273 if (ddbox->tric_dir[dimd])
8275 for (i = dd->dim[dimd]+1; i < DIM; i++)
8277 /* If we are shifted in dimension i
8278 * and the cell plane is tilted forward
8279 * in dimension i, skip this coupling.
8281 if (!(zones->shift[nzone+zone][i] &&
8282 ddbox->v[dimd][i][dimd] >= 0))
8285 sqr(ddbox->v[dimd][i][dimd]);
8288 sf2_round[dimd] = 1/sf2_round[dimd];
8293 zonei = zone_perm[dim_ind][zone];
8296 /* Here we permutate the zones to obtain a convenient order
8297 * for neighbor searching
8299 cg0 = zone_cg_range[zonei];
8300 cg1 = zone_cg_range[zonei+1];
8304 /* Look only at the cg's received in the previous grid pulse
8306 cg1 = zone_cg_range[nzone+zone+1];
8307 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8310 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8311 for (th = 0; th < comm->nth; th++)
8313 gmx_domdec_ind_t *ind_p;
8314 int **ibuf_p, *ibuf_nalloc_p;
8316 int *nsend_p, *nat_p;
8322 /* Thread 0 writes in the comm buffers */
8324 ibuf_p = &comm->buf_int;
8325 ibuf_nalloc_p = &comm->nalloc_int;
8326 vbuf_p = &comm->vbuf;
8329 nsend_zone_p = &ind->nsend[zone];
8333 /* Other threads write into temp buffers */
8334 ind_p = &comm->dth[th].ind;
8335 ibuf_p = &comm->dth[th].ibuf;
8336 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8337 vbuf_p = &comm->dth[th].vbuf;
8338 nsend_p = &comm->dth[th].nsend;
8339 nat_p = &comm->dth[th].nat;
8340 nsend_zone_p = &comm->dth[th].nsend_zone;
8342 comm->dth[th].nsend = 0;
8343 comm->dth[th].nat = 0;
8344 comm->dth[th].nsend_zone = 0;
8354 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8355 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8358 /* Get the cg's for this pulse in this zone */
8359 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8361 dim, dim_ind, dim0, dim1, dim2,
8364 normal, skew_fac2_d, skew_fac_01,
8365 v_d, v_0, v_1, &corners, sf2_round,
8366 bDistBonded, bBondComm,
8370 ibuf_p, ibuf_nalloc_p,
8376 /* Append data of threads>=1 to the communication buffers */
8377 for (th = 1; th < comm->nth; th++)
8379 dd_comm_setup_work_t *dth;
8382 dth = &comm->dth[th];
8384 ns1 = nsend + dth->nsend_zone;
8385 if (ns1 > ind->nalloc)
8387 ind->nalloc = over_alloc_dd(ns1);
8388 srenew(ind->index, ind->nalloc);
8390 if (ns1 > comm->nalloc_int)
8392 comm->nalloc_int = over_alloc_dd(ns1);
8393 srenew(comm->buf_int, comm->nalloc_int);
8395 if (ns1 > comm->vbuf.nalloc)
8397 comm->vbuf.nalloc = over_alloc_dd(ns1);
8398 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8401 for (i = 0; i < dth->nsend_zone; i++)
8403 ind->index[nsend] = dth->ind.index[i];
8404 comm->buf_int[nsend] = dth->ibuf[i];
8405 copy_rvec(dth->vbuf.v[i],
8406 comm->vbuf.v[nsend]);
8410 ind->nsend[zone] += dth->nsend_zone;
8413 /* Clear the counts in case we do not have pbc */
8414 for (zone = nzone_send; zone < nzone; zone++)
8416 ind->nsend[zone] = 0;
8418 ind->nsend[nzone] = nsend;
8419 ind->nsend[nzone+1] = nat;
8420 /* Communicate the number of cg's and atoms to receive */
8421 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8422 ind->nsend, nzone+2,
8423 ind->nrecv, nzone+2);
8425 /* The rvec buffer is also required for atom buffers of size nsend
8426 * in dd_move_x and dd_move_f.
8428 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8432 /* We can receive in place if only the last zone is not empty */
8433 for (zone = 0; zone < nzone-1; zone++)
8435 if (ind->nrecv[zone] > 0)
8437 cd->bInPlace = FALSE;
8442 /* The int buffer is only required here for the cg indices */
8443 if (ind->nrecv[nzone] > comm->nalloc_int2)
8445 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8446 srenew(comm->buf_int2, comm->nalloc_int2);
8448 /* The rvec buffer is also required for atom buffers
8449 * of size nrecv in dd_move_x and dd_move_f.
8451 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8452 vec_rvec_check_alloc(&comm->vbuf2, i);
8456 /* Make space for the global cg indices */
8457 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8458 || dd->cg_nalloc == 0)
8460 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8461 srenew(index_gl, dd->cg_nalloc);
8462 srenew(cgindex, dd->cg_nalloc+1);
8464 /* Communicate the global cg indices */
8467 recv_i = index_gl + pos_cg;
8471 recv_i = comm->buf_int2;
8473 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8474 comm->buf_int, nsend,
8475 recv_i, ind->nrecv[nzone]);
8477 /* Make space for cg_cm */
8478 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8479 if (fr->cutoff_scheme == ecutsGROUP)
8487 /* Communicate cg_cm */
8490 recv_vr = cg_cm + pos_cg;
8494 recv_vr = comm->vbuf2.v;
8496 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8497 comm->vbuf.v, nsend,
8498 recv_vr, ind->nrecv[nzone]);
8500 /* Make the charge group index */
8503 zone = (p == 0 ? 0 : nzone - 1);
8504 while (zone < nzone)
8506 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8508 cg_gl = index_gl[pos_cg];
8509 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8510 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8511 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8514 /* Update the charge group presence,
8515 * so we can use it in the next pass of the loop.
8517 comm->bLocalCG[cg_gl] = TRUE;
8523 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8526 zone_cg_range[nzone+zone] = pos_cg;
8531 /* This part of the code is never executed with bBondComm. */
8532 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8533 index_gl, recv_i, cg_cm, recv_vr,
8534 cgindex, fr->cginfo_mb, fr->cginfo);
8535 pos_cg += ind->nrecv[nzone];
8537 nat_tot += ind->nrecv[nzone+1];
8541 /* Store the atom block for easy copying of communication buffers */
8542 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8546 dd->index_gl = index_gl;
8547 dd->cgindex = cgindex;
8549 dd->ncg_tot = zone_cg_range[zones->n];
8550 dd->nat_tot = nat_tot;
8551 comm->nat[ddnatHOME] = dd->nat_home;
8552 for (i = ddnatZONE; i < ddnatNR; i++)
8554 comm->nat[i] = dd->nat_tot;
8559 /* We don't need to update cginfo, since that was alrady done above.
8560 * So we pass NULL for the forcerec.
8562 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8563 NULL, comm->bLocalCG);
8568 fprintf(debug, "Finished setting up DD communication, zones:");
8569 for (c = 0; c < zones->n; c++)
8571 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8573 fprintf(debug, "\n");
8577 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8581 for (c = 0; c < zones->nizone; c++)
8583 zones->izone[c].cg1 = zones->cg_range[c+1];
8584 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8585 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8589 static void set_zones_size(gmx_domdec_t *dd,
8590 matrix box, const gmx_ddbox_t *ddbox,
8591 int zone_start, int zone_end)
8593 gmx_domdec_comm_t *comm;
8594 gmx_domdec_zones_t *zones;
8603 zones = &comm->zones;
8605 /* Do we need to determine extra distances for multi-body bondeds? */
8606 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8608 for (z = zone_start; z < zone_end; z++)
8610 /* Copy cell limits to zone limits.
8611 * Valid for non-DD dims and non-shifted dims.
8613 copy_rvec(comm->cell_x0, zones->size[z].x0);
8614 copy_rvec(comm->cell_x1, zones->size[z].x1);
8617 for (d = 0; d < dd->ndim; d++)
8621 for (z = 0; z < zones->n; z++)
8623 /* With a staggered grid we have different sizes
8624 * for non-shifted dimensions.
8626 if (dd->bGridJump && zones->shift[z][dim] == 0)
8630 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8631 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8635 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8636 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8642 rcmbs = comm->cutoff_mbody;
8643 if (ddbox->tric_dir[dim])
8645 rcs /= ddbox->skew_fac[dim];
8646 rcmbs /= ddbox->skew_fac[dim];
8649 /* Set the lower limit for the shifted zone dimensions */
8650 for (z = zone_start; z < zone_end; z++)
8652 if (zones->shift[z][dim] > 0)
8655 if (!dd->bGridJump || d == 0)
8657 zones->size[z].x0[dim] = comm->cell_x1[dim];
8658 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8662 /* Here we take the lower limit of the zone from
8663 * the lowest domain of the zone below.
8667 zones->size[z].x0[dim] =
8668 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8674 zones->size[z].x0[dim] =
8675 zones->size[zone_perm[2][z-4]].x0[dim];
8679 zones->size[z].x0[dim] =
8680 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8683 /* A temporary limit, is updated below */
8684 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8688 for (zi = 0; zi < zones->nizone; zi++)
8690 if (zones->shift[zi][dim] == 0)
8692 /* This takes the whole zone into account.
8693 * With multiple pulses this will lead
8694 * to a larger zone then strictly necessary.
8696 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8697 zones->size[zi].x1[dim]+rcmbs);
8705 /* Loop over the i-zones to set the upper limit of each
8708 for (zi = 0; zi < zones->nizone; zi++)
8710 if (zones->shift[zi][dim] == 0)
8712 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8714 if (zones->shift[z][dim] > 0)
8716 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8717 zones->size[zi].x1[dim]+rcs);
8724 for (z = zone_start; z < zone_end; z++)
8726 /* Initialization only required to keep the compiler happy */
8727 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8730 /* To determine the bounding box for a zone we need to find
8731 * the extreme corners of 4, 2 or 1 corners.
8733 nc = 1 << (ddbox->nboundeddim - 1);
8735 for (c = 0; c < nc; c++)
8737 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8741 corner[YY] = zones->size[z].x0[YY];
8745 corner[YY] = zones->size[z].x1[YY];
8749 corner[ZZ] = zones->size[z].x0[ZZ];
8753 corner[ZZ] = zones->size[z].x1[ZZ];
8755 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8756 box[ZZ][1 - dd->dim[0]] != 0)
8758 /* With 1D domain decomposition the cg's are not in
8759 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8760 * Shift the corner of the z-vector back to along the box
8761 * vector of dimension d, so it will later end up at 0 along d.
8762 * This can affect the location of this corner along dd->dim[0]
8763 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8765 int d = 1 - dd->dim[0];
8767 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8769 /* Apply the triclinic couplings */
8770 assert(ddbox->npbcdim <= DIM);
8771 for (i = YY; i < ddbox->npbcdim; i++)
8773 for (j = XX; j < i; j++)
8775 corner[j] += corner[i]*box[i][j]/box[i][i];
8780 copy_rvec(corner, corner_min);
8781 copy_rvec(corner, corner_max);
8785 for (i = 0; i < DIM; i++)
8787 corner_min[i] = std::min(corner_min[i], corner[i]);
8788 corner_max[i] = std::max(corner_max[i], corner[i]);
8792 /* Copy the extreme cornes without offset along x */
8793 for (i = 0; i < DIM; i++)
8795 zones->size[z].bb_x0[i] = corner_min[i];
8796 zones->size[z].bb_x1[i] = corner_max[i];
8798 /* Add the offset along x */
8799 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8800 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8803 if (zone_start == 0)
8806 for (dim = 0; dim < DIM; dim++)
8808 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8810 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8815 for (z = zone_start; z < zone_end; z++)
8817 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8819 zones->size[z].x0[XX], zones->size[z].x1[XX],
8820 zones->size[z].x0[YY], zones->size[z].x1[YY],
8821 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8822 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8824 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8825 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8826 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8831 static int comp_cgsort(const void *a, const void *b)
8835 gmx_cgsort_t *cga, *cgb;
8836 cga = (gmx_cgsort_t *)a;
8837 cgb = (gmx_cgsort_t *)b;
8839 comp = cga->nsc - cgb->nsc;
8842 comp = cga->ind_gl - cgb->ind_gl;
8848 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8853 /* Order the data */
8854 for (i = 0; i < n; i++)
8856 buf[i] = a[sort[i].ind];
8859 /* Copy back to the original array */
8860 for (i = 0; i < n; i++)
8866 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8871 /* Order the data */
8872 for (i = 0; i < n; i++)
8874 copy_rvec(v[sort[i].ind], buf[i]);
8877 /* Copy back to the original array */
8878 for (i = 0; i < n; i++)
8880 copy_rvec(buf[i], v[i]);
8884 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8887 int a, atot, cg, cg0, cg1, i;
8889 if (cgindex == NULL)
8891 /* Avoid the useless loop of the atoms within a cg */
8892 order_vec_cg(ncg, sort, v, buf);
8897 /* Order the data */
8899 for (cg = 0; cg < ncg; cg++)
8901 cg0 = cgindex[sort[cg].ind];
8902 cg1 = cgindex[sort[cg].ind+1];
8903 for (i = cg0; i < cg1; i++)
8905 copy_rvec(v[i], buf[a]);
8911 /* Copy back to the original array */
8912 for (a = 0; a < atot; a++)
8914 copy_rvec(buf[a], v[a]);
8918 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8919 int nsort_new, gmx_cgsort_t *sort_new,
8920 gmx_cgsort_t *sort1)
8924 /* The new indices are not very ordered, so we qsort them */
8925 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8927 /* sort2 is already ordered, so now we can merge the two arrays */
8931 while (i2 < nsort2 || i_new < nsort_new)
8935 sort1[i1++] = sort_new[i_new++];
8937 else if (i_new == nsort_new)
8939 sort1[i1++] = sort2[i2++];
8941 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8942 (sort2[i2].nsc == sort_new[i_new].nsc &&
8943 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8945 sort1[i1++] = sort2[i2++];
8949 sort1[i1++] = sort_new[i_new++];
8954 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8956 gmx_domdec_sort_t *sort;
8957 gmx_cgsort_t *cgsort, *sort_i;
8958 int ncg_new, nsort2, nsort_new, i, *a, moved;
8960 sort = dd->comm->sort;
8962 a = fr->ns.grid->cell_index;
8964 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8966 if (ncg_home_old >= 0)
8968 /* The charge groups that remained in the same ns grid cell
8969 * are completely ordered. So we can sort efficiently by sorting
8970 * the charge groups that did move into the stationary list.
8975 for (i = 0; i < dd->ncg_home; i++)
8977 /* Check if this cg did not move to another node */
8980 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8982 /* This cg is new on this node or moved ns grid cell */
8983 if (nsort_new >= sort->sort_new_nalloc)
8985 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8986 srenew(sort->sort_new, sort->sort_new_nalloc);
8988 sort_i = &(sort->sort_new[nsort_new++]);
8992 /* This cg did not move */
8993 sort_i = &(sort->sort2[nsort2++]);
8995 /* Sort on the ns grid cell indices
8996 * and the global topology index.
8997 * index_gl is irrelevant with cell ns,
8998 * but we set it here anyhow to avoid a conditional.
9001 sort_i->ind_gl = dd->index_gl[i];
9008 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9011 /* Sort efficiently */
9012 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9017 cgsort = sort->sort;
9019 for (i = 0; i < dd->ncg_home; i++)
9021 /* Sort on the ns grid cell indices
9022 * and the global topology index
9024 cgsort[i].nsc = a[i];
9025 cgsort[i].ind_gl = dd->index_gl[i];
9027 if (cgsort[i].nsc < moved)
9034 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9036 /* Determine the order of the charge groups using qsort */
9037 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9043 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9046 int ncg_new, i, *a, na;
9048 sort = dd->comm->sort->sort;
9050 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9053 for (i = 0; i < na; i++)
9057 sort[ncg_new].ind = a[i];
9065 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9068 gmx_domdec_sort_t *sort;
9069 gmx_cgsort_t *cgsort;
9071 int ncg_new, i, *ibuf, cgsize;
9074 sort = dd->comm->sort;
9076 if (dd->ncg_home > sort->sort_nalloc)
9078 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9079 srenew(sort->sort, sort->sort_nalloc);
9080 srenew(sort->sort2, sort->sort_nalloc);
9082 cgsort = sort->sort;
9084 switch (fr->cutoff_scheme)
9087 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9090 ncg_new = dd_sort_order_nbnxn(dd, fr);
9093 gmx_incons("unimplemented");
9097 /* We alloc with the old size, since cgindex is still old */
9098 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9099 vbuf = dd->comm->vbuf.v;
9103 cgindex = dd->cgindex;
9110 /* Remove the charge groups which are no longer at home here */
9111 dd->ncg_home = ncg_new;
9114 fprintf(debug, "Set the new home charge group count to %d\n",
9118 /* Reorder the state */
9119 for (i = 0; i < estNR; i++)
9121 if (EST_DISTR(i) && (state->flags & (1<<i)))
9126 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9129 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9132 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9135 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9139 case estDISRE_INITF:
9140 case estDISRE_RM3TAV:
9141 case estORIRE_INITF:
9143 /* No ordering required */
9146 gmx_incons("Unknown state entry encountered in dd_sort_state");
9151 if (fr->cutoff_scheme == ecutsGROUP)
9154 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9157 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9159 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9160 srenew(sort->ibuf, sort->ibuf_nalloc);
9163 /* Reorder the global cg index */
9164 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9165 /* Reorder the cginfo */
9166 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9167 /* Rebuild the local cg index */
9171 for (i = 0; i < dd->ncg_home; i++)
9173 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9174 ibuf[i+1] = ibuf[i] + cgsize;
9176 for (i = 0; i < dd->ncg_home+1; i++)
9178 dd->cgindex[i] = ibuf[i];
9183 for (i = 0; i < dd->ncg_home+1; i++)
9188 /* Set the home atom number */
9189 dd->nat_home = dd->cgindex[dd->ncg_home];
9191 if (fr->cutoff_scheme == ecutsVERLET)
9193 /* The atoms are now exactly in grid order, update the grid order */
9194 nbnxn_set_atomorder(fr->nbv->nbs);
9198 /* Copy the sorted ns cell indices back to the ns grid struct */
9199 for (i = 0; i < dd->ncg_home; i++)
9201 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9203 fr->ns.grid->nr = dd->ncg_home;
9207 static void add_dd_statistics(gmx_domdec_t *dd)
9209 gmx_domdec_comm_t *comm;
9214 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9216 comm->sum_nat[ddnat-ddnatZONE] +=
9217 comm->nat[ddnat] - comm->nat[ddnat-1];
9222 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9224 gmx_domdec_comm_t *comm;
9229 /* Reset all the statistics and counters for total run counting */
9230 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9232 comm->sum_nat[ddnat-ddnatZONE] = 0;
9236 comm->load_step = 0;
9239 clear_ivec(comm->load_lim);
9244 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9246 gmx_domdec_comm_t *comm;
9250 comm = cr->dd->comm;
9252 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9259 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");
9261 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9263 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9268 " av. #atoms communicated per step for force: %d x %.1f\n",
9272 if (cr->dd->vsite_comm)
9275 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9276 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9281 if (cr->dd->constraint_comm)
9284 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9285 1 + ir->nLincsIter, av);
9289 gmx_incons(" Unknown type for DD statistics");
9292 fprintf(fplog, "\n");
9294 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9296 print_dd_load_av(fplog, cr->dd);
9300 void dd_partition_system(FILE *fplog,
9303 gmx_bool bMasterState,
9305 t_state *state_global,
9306 gmx_mtop_t *top_global,
9308 t_state *state_local,
9311 gmx_localtop_t *top_local,
9314 gmx_shellfc_t shellfc,
9315 gmx_constr_t constr,
9317 gmx_wallcycle_t wcycle,
9321 gmx_domdec_comm_t *comm;
9322 gmx_ddbox_t ddbox = {0};
9324 gmx_int64_t step_pcoupl;
9325 rvec cell_ns_x0, cell_ns_x1;
9326 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9327 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9328 gmx_bool bRedist, bSortCG, bResortAll;
9329 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9336 bBoxChanged = (bMasterState || DEFORM(*ir));
9337 if (ir->epc != epcNO)
9339 /* With nstpcouple > 1 pressure coupling happens.
9340 * one step after calculating the pressure.
9341 * Box scaling happens at the end of the MD step,
9342 * after the DD partitioning.
9343 * We therefore have to do DLB in the first partitioning
9344 * after an MD step where P-coupling occured.
9345 * We need to determine the last step in which p-coupling occurred.
9346 * MRS -- need to validate this for vv?
9351 step_pcoupl = step - 1;
9355 step_pcoupl = ((step - 1)/n)*n + 1;
9357 if (step_pcoupl >= comm->partition_step)
9363 bNStGlobalComm = (step % nstglobalcomm == 0);
9365 if (!comm->bDynLoadBal)
9371 /* Should we do dynamic load balacing this step?
9372 * Since it requires (possibly expensive) global communication,
9373 * we might want to do DLB less frequently.
9375 if (bBoxChanged || ir->epc != epcNO)
9377 bDoDLB = bBoxChanged;
9381 bDoDLB = bNStGlobalComm;
9385 /* Check if we have recorded loads on the nodes */
9386 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9388 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal && !dd_dlb_is_locked(dd))
9390 /* Check if we should use DLB at the second partitioning
9391 * and every 100 partitionings,
9392 * so the extra communication cost is negligible.
9394 const int nddp_chk_dlb = 100;
9395 bCheckDLB = (comm->n_load_collect == 0 ||
9396 comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
9403 /* Print load every nstlog, first and last step to the log file */
9404 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9405 comm->n_load_collect == 0 ||
9407 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9409 /* Avoid extra communication due to verbose screen output
9410 * when nstglobalcomm is set.
9412 if (bDoDLB || bLogLoad || bCheckDLB ||
9413 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9415 get_load_distribution(dd, wcycle);
9420 dd_print_load(fplog, dd, step-1);
9424 dd_print_load_verbose(dd);
9427 comm->n_load_collect++;
9431 /* Since the timings are node dependent, the master decides */
9434 /* Here we check if the max PME rank load is more than 0.98
9435 * the max PP force load. If so, PP DLB will not help,
9436 * since we are (almost) limited by PME. Furthermore,
9437 * DLB will cause a significant extra x/f redistribution
9438 * cost on the PME ranks, which will then surely result
9439 * in lower total performance.
9440 * This check might be fragile, since one measurement
9441 * below 0.98 (although only done once every 100 DD part.)
9442 * could turn on DLB for the rest of the run.
9444 if (cr->npmenodes > 0 &&
9445 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9452 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9456 fprintf(debug, "step %s, imb loss %f\n",
9457 gmx_step_str(step, sbuf),
9458 dd_force_imb_perf_loss(dd));
9461 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9464 turn_on_dlb(fplog, cr, step);
9469 comm->n_load_have++;
9472 cgs_gl = &comm->cgs_gl;
9477 /* Clear the old state */
9478 clear_dd_indices(dd, 0, 0);
9481 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9482 TRUE, cgs_gl, state_global->x, &ddbox);
9484 get_cg_distribution(fplog, dd, cgs_gl,
9485 state_global->box, &ddbox, state_global->x);
9487 dd_distribute_state(dd, cgs_gl,
9488 state_global, state_local, f);
9490 dd_make_local_cgs(dd, &top_local->cgs);
9492 /* Ensure that we have space for the new distribution */
9493 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9495 if (fr->cutoff_scheme == ecutsGROUP)
9497 calc_cgcm(fplog, 0, dd->ncg_home,
9498 &top_local->cgs, state_local->x, fr->cg_cm);
9501 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9503 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9505 else if (state_local->ddp_count != dd->ddp_count)
9507 if (state_local->ddp_count > dd->ddp_count)
9509 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9512 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9514 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);
9517 /* Clear the old state */
9518 clear_dd_indices(dd, 0, 0);
9520 /* Build the new indices */
9521 rebuild_cgindex(dd, cgs_gl->index, state_local);
9522 make_dd_indices(dd, cgs_gl->index, 0);
9523 ncgindex_set = dd->ncg_home;
9525 if (fr->cutoff_scheme == ecutsGROUP)
9527 /* Redetermine the cg COMs */
9528 calc_cgcm(fplog, 0, dd->ncg_home,
9529 &top_local->cgs, state_local->x, fr->cg_cm);
9532 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9534 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9536 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9537 TRUE, &top_local->cgs, state_local->x, &ddbox);
9539 bRedist = comm->bDynLoadBal;
9543 /* We have the full state, only redistribute the cgs */
9545 /* Clear the non-home indices */
9546 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9549 /* Avoid global communication for dim's without pbc and -gcom */
9550 if (!bNStGlobalComm)
9552 copy_rvec(comm->box0, ddbox.box0 );
9553 copy_rvec(comm->box_size, ddbox.box_size);
9555 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9556 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9561 /* For dim's without pbc and -gcom */
9562 copy_rvec(ddbox.box0, comm->box0 );
9563 copy_rvec(ddbox.box_size, comm->box_size);
9565 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9568 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9570 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9573 /* Check if we should sort the charge groups */
9574 if (comm->nstSortCG > 0)
9576 bSortCG = (bMasterState ||
9577 (bRedist && (step % comm->nstSortCG == 0)));
9584 ncg_home_old = dd->ncg_home;
9589 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9591 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9593 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9595 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9598 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9600 &comm->cell_x0, &comm->cell_x1,
9601 dd->ncg_home, fr->cg_cm,
9602 cell_ns_x0, cell_ns_x1, &grid_density);
9606 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9609 switch (fr->cutoff_scheme)
9612 copy_ivec(fr->ns.grid->n, ncells_old);
9613 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9614 state_local->box, cell_ns_x0, cell_ns_x1,
9615 fr->rlistlong, grid_density);
9618 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9621 gmx_incons("unimplemented");
9623 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9624 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9628 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9630 /* Sort the state on charge group position.
9631 * This enables exact restarts from this step.
9632 * It also improves performance by about 15% with larger numbers
9633 * of atoms per node.
9636 /* Fill the ns grid with the home cell,
9637 * so we can sort with the indices.
9639 set_zones_ncg_home(dd);
9641 switch (fr->cutoff_scheme)
9644 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9646 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9648 comm->zones.size[0].bb_x0,
9649 comm->zones.size[0].bb_x1,
9651 comm->zones.dens_zone0,
9654 ncg_moved, bRedist ? comm->moved : NULL,
9655 fr->nbv->grp[eintLocal].kernel_type,
9656 fr->nbv->grp[eintLocal].nbat);
9658 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9661 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9662 0, dd->ncg_home, fr->cg_cm);
9664 copy_ivec(fr->ns.grid->n, ncells_new);
9667 gmx_incons("unimplemented");
9670 bResortAll = bMasterState;
9672 /* Check if we can user the old order and ns grid cell indices
9673 * of the charge groups to sort the charge groups efficiently.
9675 if (ncells_new[XX] != ncells_old[XX] ||
9676 ncells_new[YY] != ncells_old[YY] ||
9677 ncells_new[ZZ] != ncells_old[ZZ])
9684 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9685 gmx_step_str(step, sbuf), dd->ncg_home);
9687 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9688 bResortAll ? -1 : ncg_home_old);
9689 /* Rebuild all the indices */
9690 ga2la_clear(dd->ga2la);
9693 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9696 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9698 /* Setup up the communication and communicate the coordinates */
9699 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9701 /* Set the indices */
9702 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9704 /* Set the charge group boundaries for neighbor searching */
9705 set_cg_boundaries(&comm->zones);
9707 if (fr->cutoff_scheme == ecutsVERLET)
9709 set_zones_size(dd, state_local->box, &ddbox,
9710 bSortCG ? 1 : 0, comm->zones.n);
9713 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9716 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9717 -1,state_local->x,state_local->box);
9720 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9722 /* Extract a local topology from the global topology */
9723 for (i = 0; i < dd->ndim; i++)
9725 np[dd->dim[i]] = comm->cd[i].np;
9727 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9728 comm->cellsize_min, np,
9730 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9731 vsite, top_global, top_local);
9733 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9735 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9737 /* Set up the special atom communication */
9738 n = comm->nat[ddnatZONE];
9739 for (i = ddnatZONE+1; i < ddnatNR; i++)
9744 if (vsite && vsite->n_intercg_vsite)
9746 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9750 if (dd->bInterCGcons || dd->bInterCGsettles)
9752 /* Only for inter-cg constraints we need special code */
9753 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9754 constr, ir->nProjOrder,
9755 top_local->idef.il);
9759 gmx_incons("Unknown special atom type setup");
9764 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9766 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9768 /* Make space for the extra coordinates for virtual site
9769 * or constraint communication.
9771 state_local->natoms = comm->nat[ddnatNR-1];
9772 if (state_local->natoms > state_local->nalloc)
9774 dd_realloc_state(state_local, f, state_local->natoms);
9777 if (fr->bF_NoVirSum)
9779 if (vsite && vsite->n_intercg_vsite)
9781 nat_f_novirsum = comm->nat[ddnatVSITE];
9785 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9787 nat_f_novirsum = dd->nat_tot;
9791 nat_f_novirsum = dd->nat_home;
9800 /* Set the number of atoms required for the force calculation.
9801 * Forces need to be constrained when using a twin-range setup
9802 * or with energy minimization. For simple simulations we could
9803 * avoid some allocation, zeroing and copying, but this is
9804 * probably not worth the complications ande checking.
9806 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9807 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9809 /* We make the all mdatoms up to nat_tot_con.
9810 * We could save some work by only setting invmass
9811 * between nat_tot and nat_tot_con.
9813 /* This call also sets the new number of home particles to dd->nat_home */
9814 atoms2md(top_global, ir,
9815 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9817 /* Now we have the charges we can sort the FE interactions */
9818 dd_sort_local_top(dd, mdatoms, top_local);
9822 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9823 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9824 mdatoms, FALSE, vsite);
9829 /* Make the local shell stuff, currently no communication is done */
9830 make_local_shells(cr, mdatoms, shellfc);
9833 if (ir->implicit_solvent)
9835 make_local_gb(cr, fr->born, ir->gb_algorithm);
9838 setup_bonded_threading(fr, &top_local->idef);
9840 if (!(cr->duty & DUTY_PME))
9842 /* Send the charges and/or c6/sigmas to our PME only node */
9843 gmx_pme_send_parameters(cr,
9845 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9846 mdatoms->chargeA, mdatoms->chargeB,
9847 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9848 mdatoms->sigmaA, mdatoms->sigmaB,
9849 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9854 set_constraints(constr, top_local, ir, mdatoms, cr);
9857 if (ir->ePull != epullNO)
9859 /* Update the local pull groups */
9860 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9865 /* Update the local rotation groups */
9866 dd_make_local_rotation_groups(dd, ir->rot);
9869 if (ir->eSwapCoords != eswapNO)
9871 /* Update the local groups needed for ion swapping */
9872 dd_make_local_swap_groups(dd, ir->swap);
9875 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9876 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9878 add_dd_statistics(dd);
9880 /* Make sure we only count the cycles for this DD partitioning */
9881 clear_dd_cycle_counts(dd);
9883 /* Because the order of the atoms might have changed since
9884 * the last vsite construction, we need to communicate the constructing
9885 * atom coordinates again (for spreading the forces this MD step).
9887 dd_move_x_vsites(dd, state_local->box, state_local->x);
9889 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9891 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9893 dd_move_x(dd, state_local->box, state_local->x);
9894 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9895 -1, state_local->x, state_local->box);
9898 /* Store the partitioning step */
9899 comm->partition_step = step;
9901 /* Increase the DD partitioning counter */
9903 /* The state currently matches this DD partitioning count, store it */
9904 state_local->ddp_count = dd->ddp_count;
9907 /* The DD master node knows the complete cg distribution,
9908 * store the count so we can possibly skip the cg info communication.
9910 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9913 if (comm->DD_debug > 0)
9915 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9916 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9917 "after partitioning");