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 /* With eDLB=edlbAUTO, should we check if to DLB on at the next DD? */
305 gmx_bool bCheckWhetherToTurnDlbOn;
306 /* Are we actually using DLB? */
307 gmx_bool bDynLoadBal;
309 /* Cell sizes for static load balancing, first index cartesian */
312 /* The width of the communicated boundaries */
315 /* The minimum cell size (including triclinic correction) */
317 /* For dlb, for use with edlbAUTO */
318 rvec cellsize_min_dlb;
319 /* The lower limit for the DD cell size with DLB */
321 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
322 gmx_bool bVacDLBNoLimit;
324 /* With PME load balancing we set limits on DLB */
325 gmx_bool bPMELoadBalDLBLimits;
326 /* DLB needs to take into account that we want to allow this maximum
327 * cut-off (for PME load balancing), this could limit cell boundaries.
329 real PMELoadBal_max_cutoff;
331 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
333 /* box0 and box_size are required with dim's without pbc and -gcom */
337 /* The cell boundaries */
341 /* The old location of the cell boundaries, to check cg displacements */
345 /* The communication setup and charge group boundaries for the zones */
346 gmx_domdec_zones_t zones;
348 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
349 * cell boundaries of neighboring cells for dynamic load balancing.
351 gmx_ddzone_t zone_d1[2];
352 gmx_ddzone_t zone_d2[2][2];
354 /* The coordinate/force communication setup and indices */
355 gmx_domdec_comm_dim_t cd[DIM];
356 /* The maximum number of cells to communicate with in one dimension */
359 /* Which cg distribution is stored on the master node */
360 int master_cg_ddp_count;
362 /* The number of cg's received from the direct neighbors */
363 int zone_ncg1[DD_MAXZONE];
365 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
368 /* Array for signalling if atoms have moved to another domain */
372 /* Communication buffer for general use */
376 /* Communication buffer for general use */
379 /* Temporary storage for thread parallel communication setup */
381 dd_comm_setup_work_t *dth;
383 /* Communication buffers only used with multiple grid pulses */
388 /* Communication buffers for local redistribution */
390 int cggl_flag_nalloc[DIM*2];
392 int cgcm_state_nalloc[DIM*2];
394 /* Cell sizes for dynamic load balancing */
395 gmx_domdec_root_t **root;
399 real cell_f_max0[DIM];
400 real cell_f_min1[DIM];
402 /* Stuff for load communication */
403 gmx_bool bRecordLoad;
404 gmx_domdec_load_t *load;
405 int nrank_gpu_shared;
407 MPI_Comm *mpi_comm_load;
408 MPI_Comm mpi_comm_gpu_shared;
411 /* Maximum DLB scaling per load balancing step in percent */
415 float cycl[ddCyclNr];
416 int cycl_n[ddCyclNr];
417 float cycl_max[ddCyclNr];
418 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
422 /* How many times have did we have load measurements */
424 /* How many times have we collected the load measurements */
428 double sum_nat[ddnatNR-ddnatZONE];
438 /* The last partition step */
439 gmx_int64_t partition_step;
447 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
450 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
451 #define DD_FLAG_NRCG 65535
452 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
453 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
455 /* Zone permutation required to obtain consecutive charge groups
456 * for neighbor searching.
458 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
460 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
461 * components see only j zones with that component 0.
464 /* The DD zone order */
465 static const ivec dd_zo[DD_MAXZONE] =
466 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
471 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
476 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
481 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
486 static const ivec dd_zp0[dd_zp0n] = {{0, 0, 1}};
488 /* Factors used to avoid problems due to rounding issues */
489 #define DD_CELL_MARGIN 1.0001
490 #define DD_CELL_MARGIN2 1.00005
491 /* Factor to account for pressure scaling during nstlist steps */
492 #define DD_PRES_SCALE_MARGIN 1.02
494 /* Turn on DLB when the load imbalance causes this amount of total loss.
495 * There is a bit of overhead with DLB and it's difficult to achieve
496 * a load imbalance of less than 2% with DLB.
498 #define DD_PERF_LOSS_DLB_ON 0.02
500 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
501 #define DD_PERF_LOSS_WARN 0.05
503 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
505 /* Use separate MPI send and receive commands
506 * when nnodes <= GMX_DD_NNODES_SENDRECV.
507 * This saves memory (and some copying for small nnodes).
508 * For high parallelization scatter and gather calls are used.
510 #define GMX_DD_NNODES_SENDRECV 4
514 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
516 static void index2xyz(ivec nc,int ind,ivec xyz)
518 xyz[XX] = ind % nc[XX];
519 xyz[YY] = (ind / nc[XX]) % nc[YY];
520 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
524 /* This order is required to minimize the coordinate communication in PME
525 * which uses decomposition in the x direction.
527 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
529 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
531 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
532 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
533 xyz[ZZ] = ind % nc[ZZ];
536 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
541 ddindex = dd_index(dd->nc, c);
542 if (dd->comm->bCartesianPP_PME)
544 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
546 else if (dd->comm->bCartesianPP)
549 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
560 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
562 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
565 int ddglatnr(gmx_domdec_t *dd, int i)
575 if (i >= dd->comm->nat[ddnatNR-1])
577 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
579 atnr = dd->gatindex[i] + 1;
585 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
587 return &dd->comm->cgs_gl;
590 static void vec_rvec_init(vec_rvec_t *v)
596 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
600 v->nalloc = over_alloc_dd(n);
601 srenew(v->v, v->nalloc);
605 void dd_store_state(gmx_domdec_t *dd, t_state *state)
609 if (state->ddp_count != dd->ddp_count)
611 gmx_incons("The state does not the domain decomposition state");
614 state->ncg_gl = dd->ncg_home;
615 if (state->ncg_gl > state->cg_gl_nalloc)
617 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
618 srenew(state->cg_gl, state->cg_gl_nalloc);
620 for (i = 0; i < state->ncg_gl; i++)
622 state->cg_gl[i] = dd->index_gl[i];
625 state->ddp_count_cg_gl = dd->ddp_count;
628 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
630 return &dd->comm->zones;
633 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
634 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
636 gmx_domdec_zones_t *zones;
639 zones = &dd->comm->zones;
642 while (icg >= zones->izone[izone].cg1)
651 else if (izone < zones->nizone)
653 *jcg0 = zones->izone[izone].jcg0;
657 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
658 icg, izone, zones->nizone);
661 *jcg1 = zones->izone[izone].jcg1;
663 for (d = 0; d < dd->ndim; d++)
666 shift0[dim] = zones->izone[izone].shift0[dim];
667 shift1[dim] = zones->izone[izone].shift1[dim];
668 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
670 /* A conservative approach, this can be optimized */
677 int dd_natoms_vsite(gmx_domdec_t *dd)
679 return dd->comm->nat[ddnatVSITE];
682 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
684 *at_start = dd->comm->nat[ddnatCON-1];
685 *at_end = dd->comm->nat[ddnatCON];
688 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
690 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
691 int *index, *cgindex;
692 gmx_domdec_comm_t *comm;
693 gmx_domdec_comm_dim_t *cd;
694 gmx_domdec_ind_t *ind;
695 rvec shift = {0, 0, 0}, *buf, *rbuf;
696 gmx_bool bPBC, bScrew;
700 cgindex = dd->cgindex;
705 nat_tot = dd->nat_home;
706 for (d = 0; d < dd->ndim; d++)
708 bPBC = (dd->ci[dd->dim[d]] == 0);
709 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
712 copy_rvec(box[dd->dim[d]], shift);
715 for (p = 0; p < cd->np; p++)
722 for (i = 0; i < ind->nsend[nzone]; i++)
724 at0 = cgindex[index[i]];
725 at1 = cgindex[index[i]+1];
726 for (j = at0; j < at1; j++)
728 copy_rvec(x[j], buf[n]);
735 for (i = 0; i < ind->nsend[nzone]; i++)
737 at0 = cgindex[index[i]];
738 at1 = cgindex[index[i]+1];
739 for (j = at0; j < at1; j++)
741 /* We need to shift the coordinates */
742 rvec_add(x[j], shift, buf[n]);
749 for (i = 0; i < ind->nsend[nzone]; i++)
751 at0 = cgindex[index[i]];
752 at1 = cgindex[index[i]+1];
753 for (j = at0; j < at1; j++)
756 buf[n][XX] = x[j][XX] + shift[XX];
758 * This operation requires a special shift force
759 * treatment, which is performed in calc_vir.
761 buf[n][YY] = box[YY][YY] - x[j][YY];
762 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
774 rbuf = comm->vbuf2.v;
776 /* Send and receive the coordinates */
777 dd_sendrecv_rvec(dd, d, dddirBackward,
778 buf, ind->nsend[nzone+1],
779 rbuf, ind->nrecv[nzone+1]);
783 for (zone = 0; zone < nzone; zone++)
785 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
787 copy_rvec(rbuf[j], x[i]);
792 nat_tot += ind->nrecv[nzone+1];
798 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
800 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
801 int *index, *cgindex;
802 gmx_domdec_comm_t *comm;
803 gmx_domdec_comm_dim_t *cd;
804 gmx_domdec_ind_t *ind;
808 gmx_bool bShiftForcesNeedPbc, bScrew;
812 cgindex = dd->cgindex;
816 nzone = comm->zones.n/2;
817 nat_tot = dd->nat_tot;
818 for (d = dd->ndim-1; d >= 0; d--)
820 /* Only forces in domains near the PBC boundaries need to
821 consider PBC in the treatment of fshift */
822 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
823 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
824 if (fshift == NULL && !bScrew)
826 bShiftForcesNeedPbc = FALSE;
828 /* Determine which shift vector we need */
834 for (p = cd->np-1; p >= 0; p--)
837 nat_tot -= ind->nrecv[nzone+1];
844 sbuf = comm->vbuf2.v;
846 for (zone = 0; zone < nzone; zone++)
848 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
850 copy_rvec(f[i], sbuf[j]);
855 /* Communicate the forces */
856 dd_sendrecv_rvec(dd, d, dddirForward,
857 sbuf, ind->nrecv[nzone+1],
858 buf, ind->nsend[nzone+1]);
860 /* Add the received forces */
862 if (!bShiftForcesNeedPbc)
864 for (i = 0; i < ind->nsend[nzone]; i++)
866 at0 = cgindex[index[i]];
867 at1 = cgindex[index[i]+1];
868 for (j = at0; j < at1; j++)
870 rvec_inc(f[j], buf[n]);
877 /* fshift should always be defined if this function is
878 * called when bShiftForcesNeedPbc is true */
879 assert(NULL != fshift);
880 for (i = 0; i < ind->nsend[nzone]; i++)
882 at0 = cgindex[index[i]];
883 at1 = cgindex[index[i]+1];
884 for (j = at0; j < at1; j++)
886 rvec_inc(f[j], buf[n]);
887 /* Add this force to the shift force */
888 rvec_inc(fshift[is], buf[n]);
895 for (i = 0; i < ind->nsend[nzone]; i++)
897 at0 = cgindex[index[i]];
898 at1 = cgindex[index[i]+1];
899 for (j = at0; j < at1; j++)
901 /* Rotate the force */
902 f[j][XX] += buf[n][XX];
903 f[j][YY] -= buf[n][YY];
904 f[j][ZZ] -= buf[n][ZZ];
907 /* Add this force to the shift force */
908 rvec_inc(fshift[is], buf[n]);
919 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
921 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
922 int *index, *cgindex;
923 gmx_domdec_comm_t *comm;
924 gmx_domdec_comm_dim_t *cd;
925 gmx_domdec_ind_t *ind;
930 cgindex = dd->cgindex;
932 buf = &comm->vbuf.v[0][0];
935 nat_tot = dd->nat_home;
936 for (d = 0; d < dd->ndim; d++)
939 for (p = 0; p < cd->np; p++)
944 for (i = 0; i < ind->nsend[nzone]; i++)
946 at0 = cgindex[index[i]];
947 at1 = cgindex[index[i]+1];
948 for (j = at0; j < at1; j++)
961 rbuf = &comm->vbuf2.v[0][0];
963 /* Send and receive the coordinates */
964 dd_sendrecv_real(dd, d, dddirBackward,
965 buf, ind->nsend[nzone+1],
966 rbuf, ind->nrecv[nzone+1]);
970 for (zone = 0; zone < nzone; zone++)
972 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
979 nat_tot += ind->nrecv[nzone+1];
985 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
987 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
988 int *index, *cgindex;
989 gmx_domdec_comm_t *comm;
990 gmx_domdec_comm_dim_t *cd;
991 gmx_domdec_ind_t *ind;
996 cgindex = dd->cgindex;
998 buf = &comm->vbuf.v[0][0];
1000 nzone = comm->zones.n/2;
1001 nat_tot = dd->nat_tot;
1002 for (d = dd->ndim-1; d >= 0; d--)
1005 for (p = cd->np-1; p >= 0; p--)
1008 nat_tot -= ind->nrecv[nzone+1];
1015 sbuf = &comm->vbuf2.v[0][0];
1017 for (zone = 0; zone < nzone; zone++)
1019 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
1026 /* Communicate the forces */
1027 dd_sendrecv_real(dd, d, dddirForward,
1028 sbuf, ind->nrecv[nzone+1],
1029 buf, ind->nsend[nzone+1]);
1031 /* Add the received forces */
1033 for (i = 0; i < ind->nsend[nzone]; i++)
1035 at0 = cgindex[index[i]];
1036 at1 = cgindex[index[i]+1];
1037 for (j = at0; j < at1; j++)
1048 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1050 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",
1052 zone->min0, zone->max1,
1053 zone->mch0, zone->mch0,
1054 zone->p1_0, zone->p1_1);
1058 #define DDZONECOMM_MAXZONE 5
1059 #define DDZONECOMM_BUFSIZE 3
1061 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1062 int ddimind, int direction,
1063 gmx_ddzone_t *buf_s, int n_s,
1064 gmx_ddzone_t *buf_r, int n_r)
1066 #define ZBS DDZONECOMM_BUFSIZE
1067 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1068 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1071 for (i = 0; i < n_s; i++)
1073 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1074 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1075 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1076 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1077 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1078 vbuf_s[i*ZBS+1][2] = 0;
1079 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1080 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1081 vbuf_s[i*ZBS+2][2] = 0;
1084 dd_sendrecv_rvec(dd, ddimind, direction,
1088 for (i = 0; i < n_r; i++)
1090 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1091 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1092 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1093 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1094 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1095 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1096 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1102 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1103 rvec cell_ns_x0, rvec cell_ns_x1)
1105 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1107 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1108 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1109 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1110 rvec extr_s[2], extr_r[2];
1112 real dist_d, c = 0, det;
1113 gmx_domdec_comm_t *comm;
1114 gmx_bool bPBC, bUse;
1118 for (d = 1; d < dd->ndim; d++)
1121 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1122 zp->min0 = cell_ns_x0[dim];
1123 zp->max1 = cell_ns_x1[dim];
1124 zp->min1 = cell_ns_x1[dim];
1125 zp->mch0 = cell_ns_x0[dim];
1126 zp->mch1 = cell_ns_x1[dim];
1127 zp->p1_0 = cell_ns_x0[dim];
1128 zp->p1_1 = cell_ns_x1[dim];
1131 for (d = dd->ndim-2; d >= 0; d--)
1134 bPBC = (dim < ddbox->npbcdim);
1136 /* Use an rvec to store two reals */
1137 extr_s[d][0] = comm->cell_f0[d+1];
1138 extr_s[d][1] = comm->cell_f1[d+1];
1139 extr_s[d][2] = comm->cell_f1[d+1];
1142 /* Store the extremes in the backward sending buffer,
1143 * so the get updated separately from the forward communication.
1145 for (d1 = d; d1 < dd->ndim-1; d1++)
1147 /* We invert the order to be able to use the same loop for buf_e */
1148 buf_s[pos].min0 = extr_s[d1][1];
1149 buf_s[pos].max1 = extr_s[d1][0];
1150 buf_s[pos].min1 = extr_s[d1][2];
1151 buf_s[pos].mch0 = 0;
1152 buf_s[pos].mch1 = 0;
1153 /* Store the cell corner of the dimension we communicate along */
1154 buf_s[pos].p1_0 = comm->cell_x0[dim];
1155 buf_s[pos].p1_1 = 0;
1159 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1162 if (dd->ndim == 3 && d == 0)
1164 buf_s[pos] = comm->zone_d2[0][1];
1166 buf_s[pos] = comm->zone_d1[0];
1170 /* We only need to communicate the extremes
1171 * in the forward direction
1173 npulse = comm->cd[d].np;
1176 /* Take the minimum to avoid double communication */
1177 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1181 /* Without PBC we should really not communicate over
1182 * the boundaries, but implementing that complicates
1183 * the communication setup and therefore we simply
1184 * do all communication, but ignore some data.
1186 npulse_min = npulse;
1188 for (p = 0; p < npulse_min; p++)
1190 /* Communicate the extremes forward */
1191 bUse = (bPBC || dd->ci[dim] > 0);
1193 dd_sendrecv_rvec(dd, d, dddirForward,
1194 extr_s+d, dd->ndim-d-1,
1195 extr_r+d, dd->ndim-d-1);
1199 for (d1 = d; d1 < dd->ndim-1; d1++)
1201 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1202 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1203 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1209 for (p = 0; p < npulse; p++)
1211 /* Communicate all the zone information backward */
1212 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1214 dd_sendrecv_ddzone(dd, d, dddirBackward,
1221 for (d1 = d+1; d1 < dd->ndim; d1++)
1223 /* Determine the decrease of maximum required
1224 * communication height along d1 due to the distance along d,
1225 * this avoids a lot of useless atom communication.
1227 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1229 if (ddbox->tric_dir[dim])
1231 /* c is the off-diagonal coupling between the cell planes
1232 * along directions d and d1.
1234 c = ddbox->v[dim][dd->dim[d1]][dim];
1240 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1243 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1247 /* A negative value signals out of range */
1253 /* Accumulate the extremes over all pulses */
1254 for (i = 0; i < buf_size; i++)
1258 buf_e[i] = buf_r[i];
1264 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1265 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1266 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1269 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1277 if (bUse && dh[d1] >= 0)
1279 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1280 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1283 /* Copy the received buffer to the send buffer,
1284 * to pass the data through with the next pulse.
1286 buf_s[i] = buf_r[i];
1288 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1289 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1291 /* Store the extremes */
1294 for (d1 = d; d1 < dd->ndim-1; d1++)
1296 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1297 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1298 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1302 if (d == 1 || (d == 0 && dd->ndim == 3))
1304 for (i = d; i < 2; i++)
1306 comm->zone_d2[1-d][i] = buf_e[pos];
1312 comm->zone_d1[1] = buf_e[pos];
1322 for (i = 0; i < 2; i++)
1326 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1328 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1329 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1335 for (i = 0; i < 2; i++)
1337 for (j = 0; j < 2; j++)
1341 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1343 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1344 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1348 for (d = 1; d < dd->ndim; d++)
1350 comm->cell_f_max0[d] = extr_s[d-1][0];
1351 comm->cell_f_min1[d] = extr_s[d-1][1];
1354 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1355 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1360 static void dd_collect_cg(gmx_domdec_t *dd,
1361 t_state *state_local)
1363 gmx_domdec_master_t *ma = NULL;
1364 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1366 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1368 /* The master has the correct distribution */
1372 if (state_local->ddp_count == dd->ddp_count)
1374 /* The local state and DD are in sync, use the DD indices */
1375 ncg_home = dd->ncg_home;
1377 nat_home = dd->nat_home;
1379 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1381 /* The DD is out of sync with the local state, but we have stored
1382 * the cg indices with the local state, so we can use those.
1386 cgs_gl = &dd->comm->cgs_gl;
1388 ncg_home = state_local->ncg_gl;
1389 cg = state_local->cg_gl;
1391 for (i = 0; i < ncg_home; i++)
1393 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1398 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1412 /* Collect the charge group and atom counts on the master */
1413 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1418 for (i = 0; i < dd->nnodes; i++)
1420 ma->ncg[i] = ma->ibuf[2*i];
1421 ma->nat[i] = ma->ibuf[2*i+1];
1422 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1425 /* Make byte counts and indices */
1426 for (i = 0; i < dd->nnodes; i++)
1428 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1429 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1433 fprintf(debug, "Initial charge group distribution: ");
1434 for (i = 0; i < dd->nnodes; i++)
1436 fprintf(debug, " %d", ma->ncg[i]);
1438 fprintf(debug, "\n");
1442 /* Collect the charge group indices on the master */
1444 ncg_home*sizeof(int), cg,
1445 DDMASTER(dd) ? ma->ibuf : NULL,
1446 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1447 DDMASTER(dd) ? ma->cg : NULL);
1449 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1452 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1455 gmx_domdec_master_t *ma;
1456 int n, i, c, a, nalloc = 0;
1465 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1466 dd->rank, dd->mpi_comm_all);
1471 /* Copy the master coordinates to the global array */
1472 cgs_gl = &dd->comm->cgs_gl;
1474 n = DDMASTERRANK(dd);
1476 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1478 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1480 copy_rvec(lv[a++], v[c]);
1484 for (n = 0; n < dd->nnodes; n++)
1488 if (ma->nat[n] > nalloc)
1490 nalloc = over_alloc_dd(ma->nat[n]);
1491 srenew(buf, nalloc);
1494 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1495 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1498 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1500 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1502 copy_rvec(buf[a++], v[c]);
1511 static void get_commbuffer_counts(gmx_domdec_t *dd,
1512 int **counts, int **disps)
1514 gmx_domdec_master_t *ma;
1519 /* Make the rvec count and displacment arrays */
1521 *disps = ma->ibuf + dd->nnodes;
1522 for (n = 0; n < dd->nnodes; n++)
1524 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1525 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1529 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1532 gmx_domdec_master_t *ma;
1533 int *rcounts = NULL, *disps = NULL;
1542 get_commbuffer_counts(dd, &rcounts, &disps);
1547 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1551 cgs_gl = &dd->comm->cgs_gl;
1554 for (n = 0; n < dd->nnodes; n++)
1556 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1558 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1560 copy_rvec(buf[a++], v[c]);
1567 void dd_collect_vec(gmx_domdec_t *dd,
1568 t_state *state_local, rvec *lv, rvec *v)
1570 dd_collect_cg(dd, state_local);
1572 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1574 dd_collect_vec_sendrecv(dd, lv, v);
1578 dd_collect_vec_gatherv(dd, lv, v);
1583 void dd_collect_state(gmx_domdec_t *dd,
1584 t_state *state_local, t_state *state)
1588 nh = state->nhchainlength;
1592 for (i = 0; i < efptNR; i++)
1594 state->lambda[i] = state_local->lambda[i];
1596 state->fep_state = state_local->fep_state;
1597 state->veta = state_local->veta;
1598 state->vol0 = state_local->vol0;
1599 copy_mat(state_local->box, state->box);
1600 copy_mat(state_local->boxv, state->boxv);
1601 copy_mat(state_local->svir_prev, state->svir_prev);
1602 copy_mat(state_local->fvir_prev, state->fvir_prev);
1603 copy_mat(state_local->pres_prev, state->pres_prev);
1605 for (i = 0; i < state_local->ngtc; i++)
1607 for (j = 0; j < nh; j++)
1609 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1610 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1612 state->therm_integral[i] = state_local->therm_integral[i];
1614 for (i = 0; i < state_local->nnhpres; i++)
1616 for (j = 0; j < nh; j++)
1618 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1619 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1623 for (est = 0; est < estNR; est++)
1625 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1630 dd_collect_vec(dd, state_local, state_local->x, state->x);
1633 dd_collect_vec(dd, state_local, state_local->v, state->v);
1636 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1639 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1641 case estDISRE_INITF:
1642 case estDISRE_RM3TAV:
1643 case estORIRE_INITF:
1647 gmx_incons("Unknown state entry encountered in dd_collect_state");
1653 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1659 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1662 state->nalloc = over_alloc_dd(nalloc);
1664 for (est = 0; est < estNR; est++)
1666 if (EST_DISTR(est) && (state->flags & (1<<est)))
1671 srenew(state->x, state->nalloc);
1674 srenew(state->v, state->nalloc);
1677 srenew(state->sd_X, state->nalloc);
1680 srenew(state->cg_p, state->nalloc);
1682 case estDISRE_INITF:
1683 case estDISRE_RM3TAV:
1684 case estORIRE_INITF:
1686 /* No reallocation required */
1689 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1696 srenew(*f, state->nalloc);
1700 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1703 if (nalloc > fr->cg_nalloc)
1707 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1709 fr->cg_nalloc = over_alloc_dd(nalloc);
1710 srenew(fr->cginfo, fr->cg_nalloc);
1711 if (fr->cutoff_scheme == ecutsGROUP)
1713 srenew(fr->cg_cm, fr->cg_nalloc);
1716 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1718 /* We don't use charge groups, we use x in state to set up
1719 * the atom communication.
1721 dd_realloc_state(state, f, nalloc);
1725 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1728 gmx_domdec_master_t *ma;
1729 int n, i, c, a, nalloc = 0;
1736 for (n = 0; n < dd->nnodes; n++)
1740 if (ma->nat[n] > nalloc)
1742 nalloc = over_alloc_dd(ma->nat[n]);
1743 srenew(buf, nalloc);
1745 /* Use lv as a temporary buffer */
1747 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1749 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1751 copy_rvec(v[c], buf[a++]);
1754 if (a != ma->nat[n])
1756 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1761 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1762 DDRANK(dd, n), n, dd->mpi_comm_all);
1767 n = DDMASTERRANK(dd);
1769 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1771 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1773 copy_rvec(v[c], lv[a++]);
1780 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1781 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1786 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1789 gmx_domdec_master_t *ma;
1790 int *scounts = NULL, *disps = NULL;
1798 get_commbuffer_counts(dd, &scounts, &disps);
1802 for (n = 0; n < dd->nnodes; n++)
1804 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1806 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1808 copy_rvec(v[c], buf[a++]);
1814 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1817 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1819 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1821 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1825 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1829 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1832 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1833 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1834 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1836 if (dfhist->nlambda > 0)
1838 int nlam = dfhist->nlambda;
1839 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1840 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1841 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1842 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1843 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1844 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1846 for (i = 0; i < nlam; i++)
1848 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1849 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1850 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1851 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1852 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1853 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1858 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1859 t_state *state, t_state *state_local,
1864 nh = state->nhchainlength;
1868 for (i = 0; i < efptNR; i++)
1870 state_local->lambda[i] = state->lambda[i];
1872 state_local->fep_state = state->fep_state;
1873 state_local->veta = state->veta;
1874 state_local->vol0 = state->vol0;
1875 copy_mat(state->box, state_local->box);
1876 copy_mat(state->box_rel, state_local->box_rel);
1877 copy_mat(state->boxv, state_local->boxv);
1878 copy_mat(state->svir_prev, state_local->svir_prev);
1879 copy_mat(state->fvir_prev, state_local->fvir_prev);
1880 copy_df_history(&state_local->dfhist, &state->dfhist);
1881 for (i = 0; i < state_local->ngtc; i++)
1883 for (j = 0; j < nh; j++)
1885 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1886 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1888 state_local->therm_integral[i] = state->therm_integral[i];
1890 for (i = 0; i < state_local->nnhpres; i++)
1892 for (j = 0; j < nh; j++)
1894 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1895 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1899 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1900 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1901 dd_bcast(dd, sizeof(real), &state_local->veta);
1902 dd_bcast(dd, sizeof(real), &state_local->vol0);
1903 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1904 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1905 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1906 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1907 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1908 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1909 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1910 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1911 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1912 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1914 /* communicate df_history -- required for restarting from checkpoint */
1915 dd_distribute_dfhist(dd, &state_local->dfhist);
1917 if (dd->nat_home > state_local->nalloc)
1919 dd_realloc_state(state_local, f, dd->nat_home);
1921 for (i = 0; i < estNR; i++)
1923 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1928 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1931 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1934 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1937 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1939 case estDISRE_INITF:
1940 case estDISRE_RM3TAV:
1941 case estORIRE_INITF:
1943 /* Not implemented yet */
1946 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1952 static char dim2char(int dim)
1958 case XX: c = 'X'; break;
1959 case YY: c = 'Y'; break;
1960 case ZZ: c = 'Z'; break;
1961 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1967 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1968 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1970 rvec grid_s[2], *grid_r = NULL, cx, r;
1971 char fname[STRLEN], buf[22];
1973 int a, i, d, z, y, x;
1977 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1978 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1982 snew(grid_r, 2*dd->nnodes);
1985 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1989 for (d = 0; d < DIM; d++)
1991 for (i = 0; i < DIM; i++)
1999 if (d < ddbox->npbcdim && dd->nc[d] > 1)
2001 tric[d][i] = box[i][d]/box[i][i];
2010 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2011 out = gmx_fio_fopen(fname, "w");
2012 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2014 for (i = 0; i < dd->nnodes; i++)
2016 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2017 for (d = 0; d < DIM; d++)
2019 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2021 for (z = 0; z < 2; z++)
2023 for (y = 0; y < 2; y++)
2025 for (x = 0; x < 2; x++)
2027 cx[XX] = grid_r[i*2+x][XX];
2028 cx[YY] = grid_r[i*2+y][YY];
2029 cx[ZZ] = grid_r[i*2+z][ZZ];
2031 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2032 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2036 for (d = 0; d < DIM; d++)
2038 for (x = 0; x < 4; x++)
2042 case 0: y = 1 + i*8 + 2*x; break;
2043 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2044 case 2: y = 1 + i*8 + x; break;
2046 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2050 gmx_fio_fclose(out);
2055 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2056 gmx_mtop_t *mtop, t_commrec *cr,
2057 int natoms, rvec x[], matrix box)
2059 char fname[STRLEN], buf[22];
2061 int i, ii, resnr, c;
2062 char *atomname, *resname;
2069 natoms = dd->comm->nat[ddnatVSITE];
2072 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2074 out = gmx_fio_fopen(fname, "w");
2076 fprintf(out, "TITLE %s\n", title);
2077 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2078 for (i = 0; i < natoms; i++)
2080 ii = dd->gatindex[i];
2081 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2082 if (i < dd->comm->nat[ddnatZONE])
2085 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2091 else if (i < dd->comm->nat[ddnatVSITE])
2093 b = dd->comm->zones.n;
2097 b = dd->comm->zones.n + 1;
2099 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2100 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2102 fprintf(out, "TER\n");
2104 gmx_fio_fclose(out);
2107 real dd_cutoff_multibody(const gmx_domdec_t *dd)
2109 gmx_domdec_comm_t *comm;
2116 if (comm->bInterCGBondeds)
2118 if (comm->cutoff_mbody > 0)
2120 r = comm->cutoff_mbody;
2124 /* cutoff_mbody=0 means we do not have DLB */
2125 r = comm->cellsize_min[dd->dim[0]];
2126 for (di = 1; di < dd->ndim; di++)
2128 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2130 if (comm->bBondComm)
2132 r = std::max(r, comm->cutoff_mbody);
2136 r = std::min(r, comm->cutoff);
2144 real dd_cutoff_twobody(const gmx_domdec_t *dd)
2148 r_mb = dd_cutoff_multibody(dd);
2150 return std::max(dd->comm->cutoff, r_mb);
2154 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2158 nc = dd->nc[dd->comm->cartpmedim];
2159 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2160 copy_ivec(coord, coord_pme);
2161 coord_pme[dd->comm->cartpmedim] =
2162 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2165 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2167 /* Here we assign a PME node to communicate with this DD node
2168 * by assuming that the major index of both is x.
2169 * We add cr->npmenodes/2 to obtain an even distribution.
2171 return (ddindex*npme + npme/2)/ndd;
2174 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2176 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2179 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2181 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2184 static int *dd_pmenodes(t_commrec *cr)
2189 snew(pmenodes, cr->npmenodes);
2191 for (i = 0; i < cr->dd->nnodes; i++)
2193 p0 = cr_ddindex2pmeindex(cr, i);
2194 p1 = cr_ddindex2pmeindex(cr, i+1);
2195 if (i+1 == cr->dd->nnodes || p1 > p0)
2199 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2201 pmenodes[n] = i + 1 + n;
2209 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2217 if (dd->comm->bCartesian) {
2218 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2219 dd_coords2pmecoords(dd,coords,coords_pme);
2220 copy_ivec(dd->ntot,nc);
2221 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2222 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2224 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2226 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2232 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2237 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2239 gmx_domdec_comm_t *comm;
2241 int ddindex, nodeid = -1;
2243 comm = cr->dd->comm;
2248 if (comm->bCartesianPP_PME)
2251 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2256 ddindex = dd_index(cr->dd->nc, coords);
2257 if (comm->bCartesianPP)
2259 nodeid = comm->ddindex2simnodeid[ddindex];
2265 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2277 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2280 gmx_domdec_comm_t *comm;
2287 /* This assumes a uniform x domain decomposition grid cell size */
2288 if (comm->bCartesianPP_PME)
2291 ivec coord, coord_pme;
2292 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2293 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2295 /* This is a PP node */
2296 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2297 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2301 else if (comm->bCartesianPP)
2303 if (sim_nodeid < dd->nnodes)
2305 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2310 /* This assumes DD cells with identical x coordinates
2311 * are numbered sequentially.
2313 if (dd->comm->pmenodes == NULL)
2315 if (sim_nodeid < dd->nnodes)
2317 /* The DD index equals the nodeid */
2318 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2324 while (sim_nodeid > dd->comm->pmenodes[i])
2328 if (sim_nodeid < dd->comm->pmenodes[i])
2330 pmenode = dd->comm->pmenodes[i];
2338 void get_pme_nnodes(const gmx_domdec_t *dd,
2339 int *npmenodes_x, int *npmenodes_y)
2343 *npmenodes_x = dd->comm->npmenodes_x;
2344 *npmenodes_y = dd->comm->npmenodes_y;
2353 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2354 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2358 ivec coord, coord_pme;
2362 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2365 for (x = 0; x < dd->nc[XX]; x++)
2367 for (y = 0; y < dd->nc[YY]; y++)
2369 for (z = 0; z < dd->nc[ZZ]; z++)
2371 if (dd->comm->bCartesianPP_PME)
2376 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2377 if (dd->ci[XX] == coord_pme[XX] &&
2378 dd->ci[YY] == coord_pme[YY] &&
2379 dd->ci[ZZ] == coord_pme[ZZ])
2381 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2386 /* The slab corresponds to the nodeid in the PME group */
2387 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2389 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2396 /* The last PP-only node is the peer node */
2397 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2401 fprintf(debug, "Receive coordinates from PP ranks:");
2402 for (x = 0; x < *nmy_ddnodes; x++)
2404 fprintf(debug, " %d", (*my_ddnodes)[x]);
2406 fprintf(debug, "\n");
2410 static gmx_bool receive_vir_ener(t_commrec *cr)
2412 gmx_domdec_comm_t *comm;
2417 if (cr->npmenodes < cr->dd->nnodes)
2419 comm = cr->dd->comm;
2420 if (comm->bCartesianPP_PME)
2422 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2425 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2426 coords[comm->cartpmedim]++;
2427 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2430 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2431 if (dd_simnode2pmenode(cr, rank) == pmenode)
2433 /* This is not the last PP node for pmenode */
2441 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2442 if (cr->sim_nodeid+1 < cr->nnodes &&
2443 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2445 /* This is not the last PP node for pmenode */
2454 static void set_zones_ncg_home(gmx_domdec_t *dd)
2456 gmx_domdec_zones_t *zones;
2459 zones = &dd->comm->zones;
2461 zones->cg_range[0] = 0;
2462 for (i = 1; i < zones->n+1; i++)
2464 zones->cg_range[i] = dd->ncg_home;
2466 /* zone_ncg1[0] should always be equal to ncg_home */
2467 dd->comm->zone_ncg1[0] = dd->ncg_home;
2470 static void rebuild_cgindex(gmx_domdec_t *dd,
2471 const int *gcgs_index, t_state *state)
2473 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2476 dd_cg_gl = dd->index_gl;
2477 cgindex = dd->cgindex;
2480 for (i = 0; i < state->ncg_gl; i++)
2484 dd_cg_gl[i] = cg_gl;
2485 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2489 dd->ncg_home = state->ncg_gl;
2492 set_zones_ncg_home(dd);
2495 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2497 while (cg >= cginfo_mb->cg_end)
2502 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2505 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2506 t_forcerec *fr, char *bLocalCG)
2508 cginfo_mb_t *cginfo_mb;
2514 cginfo_mb = fr->cginfo_mb;
2515 cginfo = fr->cginfo;
2517 for (cg = cg0; cg < cg1; cg++)
2519 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2523 if (bLocalCG != NULL)
2525 for (cg = cg0; cg < cg1; cg++)
2527 bLocalCG[index_gl[cg]] = TRUE;
2532 static void make_dd_indices(gmx_domdec_t *dd,
2533 const int *gcgs_index, int cg_start)
2535 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2536 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2539 if (dd->nat_tot > dd->gatindex_nalloc)
2541 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2542 srenew(dd->gatindex, dd->gatindex_nalloc);
2545 nzone = dd->comm->zones.n;
2546 zone2cg = dd->comm->zones.cg_range;
2547 zone_ncg1 = dd->comm->zone_ncg1;
2548 index_gl = dd->index_gl;
2549 gatindex = dd->gatindex;
2550 bCGs = dd->comm->bCGs;
2552 if (zone2cg[1] != dd->ncg_home)
2554 gmx_incons("dd->ncg_zone is not up to date");
2557 /* Make the local to global and global to local atom index */
2558 a = dd->cgindex[cg_start];
2559 for (zone = 0; zone < nzone; zone++)
2567 cg0 = zone2cg[zone];
2569 cg1 = zone2cg[zone+1];
2570 cg1_p1 = cg0 + zone_ncg1[zone];
2572 for (cg = cg0; cg < cg1; cg++)
2577 /* Signal that this cg is from more than one pulse away */
2580 cg_gl = index_gl[cg];
2583 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2586 ga2la_set(dd->ga2la, a_gl, a, zone1);
2592 gatindex[a] = cg_gl;
2593 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2600 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2606 if (bLocalCG == NULL)
2610 for (i = 0; i < dd->ncg_tot; i++)
2612 if (!bLocalCG[dd->index_gl[i]])
2615 "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);
2620 for (i = 0; i < ncg_sys; i++)
2627 if (ngl != dd->ncg_tot)
2629 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);
2636 static void check_index_consistency(gmx_domdec_t *dd,
2637 int natoms_sys, int ncg_sys,
2640 int nerr, ngl, i, a, cell;
2645 if (dd->comm->DD_debug > 1)
2647 snew(have, natoms_sys);
2648 for (a = 0; a < dd->nat_tot; a++)
2650 if (have[dd->gatindex[a]] > 0)
2652 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);
2656 have[dd->gatindex[a]] = a + 1;
2662 snew(have, dd->nat_tot);
2665 for (i = 0; i < natoms_sys; i++)
2667 if (ga2la_get(dd->ga2la, i, &a, &cell))
2669 if (a >= dd->nat_tot)
2671 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);
2677 if (dd->gatindex[a] != i)
2679 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);
2686 if (ngl != dd->nat_tot)
2689 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2690 dd->rank, where, ngl, dd->nat_tot);
2692 for (a = 0; a < dd->nat_tot; a++)
2697 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2698 dd->rank, where, a+1, dd->gatindex[a]+1);
2703 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2707 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2708 dd->rank, where, nerr);
2712 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2719 /* Clear the whole list without searching */
2720 ga2la_clear(dd->ga2la);
2724 for (i = a_start; i < dd->nat_tot; i++)
2726 ga2la_del(dd->ga2la, dd->gatindex[i]);
2730 bLocalCG = dd->comm->bLocalCG;
2733 for (i = cg_start; i < dd->ncg_tot; i++)
2735 bLocalCG[dd->index_gl[i]] = FALSE;
2739 dd_clear_local_vsite_indices(dd);
2741 if (dd->constraints)
2743 dd_clear_local_constraint_indices(dd);
2747 /* This function should be used for moving the domain boudaries during DLB,
2748 * for obtaining the minimum cell size. It checks the initially set limit
2749 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2750 * and, possibly, a longer cut-off limit set for PME load balancing.
2752 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2756 cellsize_min = comm->cellsize_min[dim];
2758 if (!comm->bVacDLBNoLimit)
2760 /* The cut-off might have changed, e.g. by PME load balacning,
2761 * from the value used to set comm->cellsize_min, so check it.
2763 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2765 if (comm->bPMELoadBalDLBLimits)
2767 /* Check for the cut-off limit set by the PME load balancing */
2768 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2772 return cellsize_min;
2775 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2778 real grid_jump_limit;
2780 /* The distance between the boundaries of cells at distance
2781 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2782 * and by the fact that cells should not be shifted by more than
2783 * half their size, such that cg's only shift by one cell
2784 * at redecomposition.
2786 grid_jump_limit = comm->cellsize_limit;
2787 if (!comm->bVacDLBNoLimit)
2789 if (comm->bPMELoadBalDLBLimits)
2791 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2793 grid_jump_limit = std::max(grid_jump_limit,
2794 cutoff/comm->cd[dim_ind].np);
2797 return grid_jump_limit;
2800 static gmx_bool check_grid_jump(gmx_int64_t step,
2806 gmx_domdec_comm_t *comm;
2815 for (d = 1; d < dd->ndim; d++)
2818 limit = grid_jump_limit(comm, cutoff, d);
2819 bfac = ddbox->box_size[dim];
2820 if (ddbox->tric_dir[dim])
2822 bfac *= ddbox->skew_fac[dim];
2824 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2825 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2833 /* This error should never be triggered under normal
2834 * circumstances, but you never know ...
2836 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.",
2837 gmx_step_str(step, buf),
2838 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2846 static int dd_load_count(gmx_domdec_comm_t *comm)
2848 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2851 static float dd_force_load(gmx_domdec_comm_t *comm)
2858 if (comm->eFlop > 1)
2860 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2865 load = comm->cycl[ddCyclF];
2866 if (comm->cycl_n[ddCyclF] > 1)
2868 /* Subtract the maximum of the last n cycle counts
2869 * to get rid of possible high counts due to other sources,
2870 * for instance system activity, that would otherwise
2871 * affect the dynamic load balancing.
2873 load -= comm->cycl_max[ddCyclF];
2877 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2879 float gpu_wait, gpu_wait_sum;
2881 gpu_wait = comm->cycl[ddCyclWaitGPU];
2882 if (comm->cycl_n[ddCyclF] > 1)
2884 /* We should remove the WaitGPU time of the same MD step
2885 * as the one with the maximum F time, since the F time
2886 * and the wait time are not independent.
2887 * Furthermore, the step for the max F time should be chosen
2888 * the same on all ranks that share the same GPU.
2889 * But to keep the code simple, we remove the average instead.
2890 * The main reason for artificially long times at some steps
2891 * is spurious CPU activity or MPI time, so we don't expect
2892 * that changes in the GPU wait time matter a lot here.
2894 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2896 /* Sum the wait times over the ranks that share the same GPU */
2897 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2898 comm->mpi_comm_gpu_shared);
2899 /* Replace the wait time by the average over the ranks */
2900 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2908 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2910 gmx_domdec_comm_t *comm;
2915 snew(*dim_f, dd->nc[dim]+1);
2917 for (i = 1; i < dd->nc[dim]; i++)
2919 if (comm->slb_frac[dim])
2921 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2925 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2928 (*dim_f)[dd->nc[dim]] = 1;
2931 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2933 int pmeindex, slab, nso, i;
2936 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2942 ddpme->dim = dimind;
2944 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2946 ddpme->nslab = (ddpme->dim == 0 ?
2947 dd->comm->npmenodes_x :
2948 dd->comm->npmenodes_y);
2950 if (ddpme->nslab <= 1)
2955 nso = dd->comm->npmenodes/ddpme->nslab;
2956 /* Determine for each PME slab the PP location range for dimension dim */
2957 snew(ddpme->pp_min, ddpme->nslab);
2958 snew(ddpme->pp_max, ddpme->nslab);
2959 for (slab = 0; slab < ddpme->nslab; slab++)
2961 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2962 ddpme->pp_max[slab] = 0;
2964 for (i = 0; i < dd->nnodes; i++)
2966 ddindex2xyz(dd->nc, i, xyz);
2967 /* For y only use our y/z slab.
2968 * This assumes that the PME x grid size matches the DD grid size.
2970 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2972 pmeindex = ddindex2pmeindex(dd, i);
2975 slab = pmeindex/nso;
2979 slab = pmeindex % ddpme->nslab;
2981 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2982 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2986 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2989 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2991 if (dd->comm->ddpme[0].dim == XX)
2993 return dd->comm->ddpme[0].maxshift;
3001 int dd_pme_maxshift_y(gmx_domdec_t *dd)
3003 if (dd->comm->ddpme[0].dim == YY)
3005 return dd->comm->ddpme[0].maxshift;
3007 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3009 return dd->comm->ddpme[1].maxshift;
3017 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3018 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3020 gmx_domdec_comm_t *comm;
3023 real range, pme_boundary;
3027 nc = dd->nc[ddpme->dim];
3030 if (!ddpme->dim_match)
3032 /* PP decomposition is not along dim: the worst situation */
3035 else if (ns <= 3 || (bUniform && ns == nc))
3037 /* The optimal situation */
3042 /* We need to check for all pme nodes which nodes they
3043 * could possibly need to communicate with.
3045 xmin = ddpme->pp_min;
3046 xmax = ddpme->pp_max;
3047 /* Allow for atoms to be maximally 2/3 times the cut-off
3048 * out of their DD cell. This is a reasonable balance between
3049 * between performance and support for most charge-group/cut-off
3052 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3053 /* Avoid extra communication when we are exactly at a boundary */
3057 for (s = 0; s < ns; s++)
3059 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3060 pme_boundary = (real)s/ns;
3063 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3065 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3069 pme_boundary = (real)(s+1)/ns;
3072 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3074 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3081 ddpme->maxshift = sh;
3085 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3086 ddpme->dim, ddpme->maxshift);
3090 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3094 for (d = 0; d < dd->ndim; d++)
3097 if (dim < ddbox->nboundeddim &&
3098 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3099 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3101 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",
3102 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3103 dd->nc[dim], dd->comm->cellsize_limit);
3109 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3112 /* Set the domain boundaries. Use for static (or no) load balancing,
3113 * and also for the starting state for dynamic load balancing.
3114 * setmode determine if and where the boundaries are stored, use enum above.
3115 * Returns the number communication pulses in npulse.
3117 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3118 int setmode, ivec npulse)
3120 gmx_domdec_comm_t *comm;
3123 real *cell_x, cell_dx, cellsize;
3127 for (d = 0; d < DIM; d++)
3129 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3131 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3134 cell_dx = ddbox->box_size[d]/dd->nc[d];
3137 case setcellsizeslbMASTER:
3138 for (j = 0; j < dd->nc[d]+1; j++)
3140 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3143 case setcellsizeslbLOCAL:
3144 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3145 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3150 cellsize = cell_dx*ddbox->skew_fac[d];
3151 while (cellsize*npulse[d] < comm->cutoff)
3155 cellsize_min[d] = cellsize;
3159 /* Statically load balanced grid */
3160 /* Also when we are not doing a master distribution we determine
3161 * all cell borders in a loop to obtain identical values
3162 * to the master distribution case and to determine npulse.
3164 if (setmode == setcellsizeslbMASTER)
3166 cell_x = dd->ma->cell_x[d];
3170 snew(cell_x, dd->nc[d]+1);
3172 cell_x[0] = ddbox->box0[d];
3173 for (j = 0; j < dd->nc[d]; j++)
3175 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3176 cell_x[j+1] = cell_x[j] + cell_dx;
3177 cellsize = cell_dx*ddbox->skew_fac[d];
3178 while (cellsize*npulse[d] < comm->cutoff &&
3179 npulse[d] < dd->nc[d]-1)
3183 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3185 if (setmode == setcellsizeslbLOCAL)
3187 comm->cell_x0[d] = cell_x[dd->ci[d]];
3188 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3190 if (setmode != setcellsizeslbMASTER)
3195 /* The following limitation is to avoid that a cell would receive
3196 * some of its own home charge groups back over the periodic boundary.
3197 * Double charge groups cause trouble with the global indices.
3199 if (d < ddbox->npbcdim &&
3200 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3202 char error_string[STRLEN];
3204 sprintf(error_string,
3205 "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",
3206 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3208 dd->nc[d], dd->nc[d],
3209 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3211 if (setmode == setcellsizeslbLOCAL)
3213 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3217 gmx_fatal(FARGS, error_string);
3222 if (!comm->bDynLoadBal)
3224 copy_rvec(cellsize_min, comm->cellsize_min);
3227 for (d = 0; d < comm->npmedecompdim; d++)
3229 set_pme_maxshift(dd, &comm->ddpme[d],
3230 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3231 comm->ddpme[d].slb_dim_f);
3236 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3237 int d, int dim, gmx_domdec_root_t *root,
3239 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3241 gmx_domdec_comm_t *comm;
3242 int ncd, i, j, nmin, nmin_old;
3243 gmx_bool bLimLo, bLimHi;
3245 real fac, halfway, cellsize_limit_f_i, region_size;
3246 gmx_bool bPBC, bLastHi = FALSE;
3247 int nrange[] = {range[0], range[1]};
3249 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3255 bPBC = (dim < ddbox->npbcdim);
3257 cell_size = root->buf_ncd;
3261 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3264 /* First we need to check if the scaling does not make cells
3265 * smaller than the smallest allowed size.
3266 * We need to do this iteratively, since if a cell is too small,
3267 * it needs to be enlarged, which makes all the other cells smaller,
3268 * which could in turn make another cell smaller than allowed.
3270 for (i = range[0]; i < range[1]; i++)
3272 root->bCellMin[i] = FALSE;
3278 /* We need the total for normalization */
3280 for (i = range[0]; i < range[1]; i++)
3282 if (root->bCellMin[i] == FALSE)
3284 fac += cell_size[i];
3287 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3288 /* Determine the cell boundaries */
3289 for (i = range[0]; i < range[1]; i++)
3291 if (root->bCellMin[i] == FALSE)
3293 cell_size[i] *= fac;
3294 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3296 cellsize_limit_f_i = 0;
3300 cellsize_limit_f_i = cellsize_limit_f;
3302 if (cell_size[i] < cellsize_limit_f_i)
3304 root->bCellMin[i] = TRUE;
3305 cell_size[i] = cellsize_limit_f_i;
3309 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3312 while (nmin > nmin_old);
3315 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3316 /* For this check we should not use DD_CELL_MARGIN,
3317 * but a slightly smaller factor,
3318 * since rounding could get use below the limit.
3320 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3323 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",
3324 gmx_step_str(step, buf),
3325 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3326 ncd, comm->cellsize_min[dim]);
3329 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3333 /* Check if the boundary did not displace more than halfway
3334 * each of the cells it bounds, as this could cause problems,
3335 * especially when the differences between cell sizes are large.
3336 * If changes are applied, they will not make cells smaller
3337 * than the cut-off, as we check all the boundaries which
3338 * might be affected by a change and if the old state was ok,
3339 * the cells will at most be shrunk back to their old size.
3341 for (i = range[0]+1; i < range[1]; i++)
3343 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3344 if (root->cell_f[i] < halfway)
3346 root->cell_f[i] = halfway;
3347 /* Check if the change also causes shifts of the next boundaries */
3348 for (j = i+1; j < range[1]; j++)
3350 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3352 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3356 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3357 if (root->cell_f[i] > halfway)
3359 root->cell_f[i] = halfway;
3360 /* Check if the change also causes shifts of the next boundaries */
3361 for (j = i-1; j >= range[0]+1; j--)
3363 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3365 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3372 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3373 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3374 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3375 * for a and b nrange is used */
3378 /* Take care of the staggering of the cell boundaries */
3381 for (i = range[0]; i < range[1]; i++)
3383 root->cell_f_max0[i] = root->cell_f[i];
3384 root->cell_f_min1[i] = root->cell_f[i+1];
3389 for (i = range[0]+1; i < range[1]; i++)
3391 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3392 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3393 if (bLimLo && bLimHi)
3395 /* Both limits violated, try the best we can */
3396 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3397 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3398 nrange[0] = range[0];
3400 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3403 nrange[1] = range[1];
3404 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3410 /* root->cell_f[i] = root->bound_min[i]; */
3411 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3414 else if (bLimHi && !bLastHi)
3417 if (nrange[1] < range[1]) /* found a LimLo before */
3419 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3420 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3421 nrange[0] = nrange[1];
3423 root->cell_f[i] = root->bound_max[i];
3425 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3427 nrange[1] = range[1];
3430 if (nrange[1] < range[1]) /* found last a LimLo */
3432 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3433 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3434 nrange[0] = nrange[1];
3435 nrange[1] = range[1];
3436 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3438 else if (nrange[0] > range[0]) /* found at least one LimHi */
3440 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3447 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3448 int d, int dim, gmx_domdec_root_t *root,
3449 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3450 gmx_bool bUniform, gmx_int64_t step)
3452 gmx_domdec_comm_t *comm;
3453 int ncd, d1, i, pos;
3455 real load_aver, load_i, imbalance, change, change_max, sc;
3456 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3460 int range[] = { 0, 0 };
3464 /* Convert the maximum change from the input percentage to a fraction */
3465 change_limit = comm->dlb_scale_lim*0.01;
3469 bPBC = (dim < ddbox->npbcdim);
3471 cell_size = root->buf_ncd;
3473 /* Store the original boundaries */
3474 for (i = 0; i < ncd+1; i++)
3476 root->old_cell_f[i] = root->cell_f[i];
3480 for (i = 0; i < ncd; i++)
3482 cell_size[i] = 1.0/ncd;
3485 else if (dd_load_count(comm) > 0)
3487 load_aver = comm->load[d].sum_m/ncd;
3489 for (i = 0; i < ncd; i++)
3491 /* Determine the relative imbalance of cell i */
3492 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3493 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3494 /* Determine the change of the cell size using underrelaxation */
3495 change = -relax*imbalance;
3496 change_max = std::max(change_max, std::max(change, -change));
3498 /* Limit the amount of scaling.
3499 * We need to use the same rescaling for all cells in one row,
3500 * otherwise the load balancing might not converge.
3503 if (change_max > change_limit)
3505 sc *= change_limit/change_max;
3507 for (i = 0; i < ncd; i++)
3509 /* Determine the relative imbalance of cell i */
3510 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3511 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3512 /* Determine the change of the cell size using underrelaxation */
3513 change = -sc*imbalance;
3514 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3518 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3519 cellsize_limit_f *= DD_CELL_MARGIN;
3520 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3521 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3522 if (ddbox->tric_dir[dim])
3524 cellsize_limit_f /= ddbox->skew_fac[dim];
3525 dist_min_f /= ddbox->skew_fac[dim];
3527 if (bDynamicBox && d > 0)
3529 dist_min_f *= DD_PRES_SCALE_MARGIN;
3531 if (d > 0 && !bUniform)
3533 /* Make sure that the grid is not shifted too much */
3534 for (i = 1; i < ncd; i++)
3536 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3538 gmx_incons("Inconsistent DD boundary staggering limits!");
3540 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3541 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3544 root->bound_min[i] += 0.5*space;
3546 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3547 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3550 root->bound_max[i] += 0.5*space;
3555 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3557 root->cell_f_max0[i-1] + dist_min_f,
3558 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3559 root->cell_f_min1[i] - dist_min_f);
3564 root->cell_f[0] = 0;
3565 root->cell_f[ncd] = 1;
3566 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3569 /* After the checks above, the cells should obey the cut-off
3570 * restrictions, but it does not hurt to check.
3572 for (i = 0; i < ncd; i++)
3576 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3577 dim, i, root->cell_f[i], root->cell_f[i+1]);
3580 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3581 root->cell_f[i+1] - root->cell_f[i] <
3582 cellsize_limit_f/DD_CELL_MARGIN)
3586 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3587 gmx_step_str(step, buf), dim2char(dim), i,
3588 (root->cell_f[i+1] - root->cell_f[i])
3589 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3594 /* Store the cell boundaries of the lower dimensions at the end */
3595 for (d1 = 0; d1 < d; d1++)
3597 root->cell_f[pos++] = comm->cell_f0[d1];
3598 root->cell_f[pos++] = comm->cell_f1[d1];
3601 if (d < comm->npmedecompdim)
3603 /* The master determines the maximum shift for
3604 * the coordinate communication between separate PME nodes.
3606 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3608 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3611 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3615 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3616 gmx_ddbox_t *ddbox, int dimind)
3618 gmx_domdec_comm_t *comm;
3623 /* Set the cell dimensions */
3624 dim = dd->dim[dimind];
3625 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3626 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3627 if (dim >= ddbox->nboundeddim)
3629 comm->cell_x0[dim] += ddbox->box0[dim];
3630 comm->cell_x1[dim] += ddbox->box0[dim];
3634 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3635 int d, int dim, real *cell_f_row,
3638 gmx_domdec_comm_t *comm;
3644 /* Each node would only need to know two fractions,
3645 * but it is probably cheaper to broadcast the whole array.
3647 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3648 0, comm->mpi_comm_load[d]);
3650 /* Copy the fractions for this dimension from the buffer */
3651 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3652 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3653 /* The whole array was communicated, so set the buffer position */
3654 pos = dd->nc[dim] + 1;
3655 for (d1 = 0; d1 <= d; d1++)
3659 /* Copy the cell fractions of the lower dimensions */
3660 comm->cell_f0[d1] = cell_f_row[pos++];
3661 comm->cell_f1[d1] = cell_f_row[pos++];
3663 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3665 /* Convert the communicated shift from float to int */
3666 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3669 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3673 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3674 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3675 gmx_bool bUniform, gmx_int64_t step)
3677 gmx_domdec_comm_t *comm;
3679 gmx_bool bRowMember, bRowRoot;
3684 for (d = 0; d < dd->ndim; d++)
3689 for (d1 = d; d1 < dd->ndim; d1++)
3691 if (dd->ci[dd->dim[d1]] > 0)
3704 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3705 ddbox, bDynamicBox, bUniform, step);
3706 cell_f_row = comm->root[d]->cell_f;
3710 cell_f_row = comm->cell_f_row;
3712 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3717 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3721 /* This function assumes the box is static and should therefore
3722 * not be called when the box has changed since the last
3723 * call to dd_partition_system.
3725 for (d = 0; d < dd->ndim; d++)
3727 relative_to_absolute_cell_bounds(dd, ddbox, d);
3733 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3734 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3735 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3736 gmx_wallcycle_t wcycle)
3738 gmx_domdec_comm_t *comm;
3745 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3746 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3747 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3749 else if (bDynamicBox)
3751 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3754 /* Set the dimensions for which no DD is used */
3755 for (dim = 0; dim < DIM; dim++)
3757 if (dd->nc[dim] == 1)
3759 comm->cell_x0[dim] = 0;
3760 comm->cell_x1[dim] = ddbox->box_size[dim];
3761 if (dim >= ddbox->nboundeddim)
3763 comm->cell_x0[dim] += ddbox->box0[dim];
3764 comm->cell_x1[dim] += ddbox->box0[dim];
3770 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3773 gmx_domdec_comm_dim_t *cd;
3775 for (d = 0; d < dd->ndim; d++)
3777 cd = &dd->comm->cd[d];
3778 np = npulse[dd->dim[d]];
3779 if (np > cd->np_nalloc)
3783 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3784 dim2char(dd->dim[d]), np);
3786 if (DDMASTER(dd) && cd->np_nalloc > 0)
3788 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3790 srenew(cd->ind, np);
3791 for (i = cd->np_nalloc; i < np; i++)
3793 cd->ind[i].index = NULL;
3794 cd->ind[i].nalloc = 0;
3803 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3804 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3805 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3806 gmx_wallcycle_t wcycle)
3808 gmx_domdec_comm_t *comm;
3814 /* Copy the old cell boundaries for the cg displacement check */
3815 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3816 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3818 if (comm->bDynLoadBal)
3822 check_box_size(dd, ddbox);
3824 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3828 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3829 realloc_comm_ind(dd, npulse);
3834 for (d = 0; d < DIM; d++)
3836 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3837 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3842 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3844 rvec cell_ns_x0, rvec cell_ns_x1,
3847 gmx_domdec_comm_t *comm;
3852 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3854 dim = dd->dim[dim_ind];
3856 /* Without PBC we don't have restrictions on the outer cells */
3857 if (!(dim >= ddbox->npbcdim &&
3858 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3859 comm->bDynLoadBal &&
3860 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3861 comm->cellsize_min[dim])
3864 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",
3865 gmx_step_str(step, buf), dim2char(dim),
3866 comm->cell_x1[dim] - comm->cell_x0[dim],
3867 ddbox->skew_fac[dim],
3868 dd->comm->cellsize_min[dim],
3869 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3873 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3875 /* Communicate the boundaries and update cell_ns_x0/1 */
3876 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3877 if (dd->bGridJump && dd->ndim > 1)
3879 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3884 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3888 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3896 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3897 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3906 static void check_screw_box(matrix box)
3908 /* Mathematical limitation */
3909 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3911 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3914 /* Limitation due to the asymmetry of the eighth shell method */
3915 if (box[ZZ][YY] != 0)
3917 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3921 static void distribute_cg(FILE *fplog,
3922 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3925 gmx_domdec_master_t *ma;
3926 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3927 int i, icg, j, k, k0, k1, d;
3931 real nrcg, inv_ncg, pos_d;
3937 if (tmp_ind == NULL)
3939 snew(tmp_nalloc, dd->nnodes);
3940 snew(tmp_ind, dd->nnodes);
3941 for (i = 0; i < dd->nnodes; i++)
3943 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3944 snew(tmp_ind[i], tmp_nalloc[i]);
3948 /* Clear the count */
3949 for (i = 0; i < dd->nnodes; i++)
3955 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3957 cgindex = cgs->index;
3959 /* Compute the center of geometry for all charge groups */
3960 for (icg = 0; icg < cgs->nr; icg++)
3963 k1 = cgindex[icg+1];
3967 copy_rvec(pos[k0], cg_cm);
3974 for (k = k0; (k < k1); k++)
3976 rvec_inc(cg_cm, pos[k]);
3978 for (d = 0; (d < DIM); d++)
3980 cg_cm[d] *= inv_ncg;
3983 /* Put the charge group in the box and determine the cell index */
3984 for (d = DIM-1; d >= 0; d--)
3987 if (d < dd->npbcdim)
3989 bScrew = (dd->bScrewPBC && d == XX);
3990 if (tric_dir[d] && dd->nc[d] > 1)
3992 /* Use triclinic coordintates for this dimension */
3993 for (j = d+1; j < DIM; j++)
3995 pos_d += cg_cm[j]*tcm[j][d];
3998 while (pos_d >= box[d][d])
4001 rvec_dec(cg_cm, box[d]);
4004 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4005 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4007 for (k = k0; (k < k1); k++)
4009 rvec_dec(pos[k], box[d]);
4012 pos[k][YY] = box[YY][YY] - pos[k][YY];
4013 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4020 rvec_inc(cg_cm, box[d]);
4023 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4024 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4026 for (k = k0; (k < k1); k++)
4028 rvec_inc(pos[k], box[d]);
4031 pos[k][YY] = box[YY][YY] - pos[k][YY];
4032 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4037 /* This could be done more efficiently */
4039 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4044 i = dd_index(dd->nc, ind);
4045 if (ma->ncg[i] == tmp_nalloc[i])
4047 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4048 srenew(tmp_ind[i], tmp_nalloc[i]);
4050 tmp_ind[i][ma->ncg[i]] = icg;
4052 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4056 for (i = 0; i < dd->nnodes; i++)
4059 for (k = 0; k < ma->ncg[i]; k++)
4061 ma->cg[k1++] = tmp_ind[i][k];
4064 ma->index[dd->nnodes] = k1;
4066 for (i = 0; i < dd->nnodes; i++)
4075 /* Here we avoid int overflows due to #atoms^2: use double, dsqr */
4076 int nat_sum, nat_min, nat_max;
4081 nat_min = ma->nat[0];
4082 nat_max = ma->nat[0];
4083 for (i = 0; i < dd->nnodes; i++)
4085 nat_sum += ma->nat[i];
4086 nat2_sum += dsqr(ma->nat[i]);
4087 nat_min = std::min(nat_min, ma->nat[i]);
4088 nat_max = std::max(nat_max, ma->nat[i]);
4090 nat_sum /= dd->nnodes;
4091 nat2_sum /= dd->nnodes;
4093 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
4096 static_cast<int>(sqrt(nat2_sum - dsqr(nat_sum) + 0.5)),
4101 static void get_cg_distribution(FILE *fplog, gmx_domdec_t *dd,
4102 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4105 gmx_domdec_master_t *ma = NULL;
4108 int *ibuf, buf2[2] = { 0, 0 };
4109 gmx_bool bMaster = DDMASTER(dd);
4117 check_screw_box(box);
4120 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4122 distribute_cg(fplog, box, ddbox->tric_dir, cgs, pos, dd);
4123 for (i = 0; i < dd->nnodes; i++)
4125 ma->ibuf[2*i] = ma->ncg[i];
4126 ma->ibuf[2*i+1] = ma->nat[i];
4134 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4136 dd->ncg_home = buf2[0];
4137 dd->nat_home = buf2[1];
4138 dd->ncg_tot = dd->ncg_home;
4139 dd->nat_tot = dd->nat_home;
4140 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4142 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4143 srenew(dd->index_gl, dd->cg_nalloc);
4144 srenew(dd->cgindex, dd->cg_nalloc+1);
4148 for (i = 0; i < dd->nnodes; i++)
4150 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4151 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4156 bMaster ? ma->ibuf : NULL,
4157 bMaster ? ma->ibuf+dd->nnodes : NULL,
4158 bMaster ? ma->cg : NULL,
4159 dd->ncg_home*sizeof(int), dd->index_gl);
4161 /* Determine the home charge group sizes */
4163 for (i = 0; i < dd->ncg_home; i++)
4165 cg_gl = dd->index_gl[i];
4167 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4172 fprintf(debug, "Home charge groups:\n");
4173 for (i = 0; i < dd->ncg_home; i++)
4175 fprintf(debug, " %d", dd->index_gl[i]);
4178 fprintf(debug, "\n");
4181 fprintf(debug, "\n");
4185 static int compact_and_copy_vec_at(int ncg, int *move,
4188 rvec *src, gmx_domdec_comm_t *comm,
4191 int m, icg, i, i0, i1, nrcg;
4197 for (m = 0; m < DIM*2; m++)
4203 for (icg = 0; icg < ncg; icg++)
4205 i1 = cgindex[icg+1];
4211 /* Compact the home array in place */
4212 for (i = i0; i < i1; i++)
4214 copy_rvec(src[i], src[home_pos++]);
4220 /* Copy to the communication buffer */
4222 pos_vec[m] += 1 + vec*nrcg;
4223 for (i = i0; i < i1; i++)
4225 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4227 pos_vec[m] += (nvec - vec - 1)*nrcg;
4231 home_pos += i1 - i0;
4239 static int compact_and_copy_vec_cg(int ncg, int *move,
4241 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4244 int m, icg, i0, i1, nrcg;
4250 for (m = 0; m < DIM*2; m++)
4256 for (icg = 0; icg < ncg; icg++)
4258 i1 = cgindex[icg+1];
4264 /* Compact the home array in place */
4265 copy_rvec(src[icg], src[home_pos++]);
4271 /* Copy to the communication buffer */
4272 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4273 pos_vec[m] += 1 + nrcg*nvec;
4285 static int compact_ind(int ncg, int *move,
4286 int *index_gl, int *cgindex,
4288 gmx_ga2la_t ga2la, char *bLocalCG,
4291 int cg, nat, a0, a1, a, a_gl;
4296 for (cg = 0; cg < ncg; cg++)
4302 /* Compact the home arrays in place.
4303 * Anything that can be done here avoids access to global arrays.
4305 cgindex[home_pos] = nat;
4306 for (a = a0; a < a1; a++)
4309 gatindex[nat] = a_gl;
4310 /* The cell number stays 0, so we don't need to set it */
4311 ga2la_change_la(ga2la, a_gl, nat);
4314 index_gl[home_pos] = index_gl[cg];
4315 cginfo[home_pos] = cginfo[cg];
4316 /* The charge group remains local, so bLocalCG does not change */
4321 /* Clear the global indices */
4322 for (a = a0; a < a1; a++)
4324 ga2la_del(ga2la, gatindex[a]);
4328 bLocalCG[index_gl[cg]] = FALSE;
4332 cgindex[home_pos] = nat;
4337 static void clear_and_mark_ind(int ncg, int *move,
4338 int *index_gl, int *cgindex, int *gatindex,
4339 gmx_ga2la_t ga2la, char *bLocalCG,
4344 for (cg = 0; cg < ncg; cg++)
4350 /* Clear the global indices */
4351 for (a = a0; a < a1; a++)
4353 ga2la_del(ga2la, gatindex[a]);
4357 bLocalCG[index_gl[cg]] = FALSE;
4359 /* Signal that this cg has moved using the ns cell index.
4360 * Here we set it to -1. fill_grid will change it
4361 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4363 cell_index[cg] = -1;
4368 static void print_cg_move(FILE *fplog,
4370 gmx_int64_t step, int cg, int dim, int dir,
4371 gmx_bool bHaveCgcmOld, real limitd,
4372 rvec cm_old, rvec cm_new, real pos_d)
4374 gmx_domdec_comm_t *comm;
4379 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4382 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4383 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4384 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4388 /* We don't have a limiting distance available: don't print it */
4389 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4390 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4391 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4393 fprintf(fplog, "distance out of cell %f\n",
4394 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4397 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4398 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4400 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4401 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4402 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4404 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4405 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4407 comm->cell_x0[dim], comm->cell_x1[dim]);
4410 static void cg_move_error(FILE *fplog,
4412 gmx_int64_t step, int cg, int dim, int dir,
4413 gmx_bool bHaveCgcmOld, real limitd,
4414 rvec cm_old, rvec cm_new, real pos_d)
4418 print_cg_move(fplog, dd, step, cg, dim, dir,
4419 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4421 print_cg_move(stderr, dd, step, cg, dim, dir,
4422 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4424 "%s moved too far between two domain decomposition steps\n"
4425 "This usually means that your system is not well equilibrated",
4426 dd->comm->bCGs ? "A charge group" : "An atom");
4429 static void rotate_state_atom(t_state *state, int a)
4433 for (est = 0; est < estNR; est++)
4435 if (EST_DISTR(est) && (state->flags & (1<<est)))
4440 /* Rotate the complete state; for a rectangular box only */
4441 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4442 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4445 state->v[a][YY] = -state->v[a][YY];
4446 state->v[a][ZZ] = -state->v[a][ZZ];
4449 state->sd_X[a][YY] = -state->sd_X[a][YY];
4450 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4453 state->cg_p[a][YY] = -state->cg_p[a][YY];
4454 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4456 case estDISRE_INITF:
4457 case estDISRE_RM3TAV:
4458 case estORIRE_INITF:
4460 /* These are distances, so not affected by rotation */
4463 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4469 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4471 if (natoms > comm->moved_nalloc)
4473 /* Contents should be preserved here */
4474 comm->moved_nalloc = over_alloc_dd(natoms);
4475 srenew(comm->moved, comm->moved_nalloc);
4481 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4484 ivec tric_dir, matrix tcm,
4485 rvec cell_x0, rvec cell_x1,
4486 rvec limitd, rvec limit0, rvec limit1,
4488 int cg_start, int cg_end,
4493 int cg, k, k0, k1, d, dim, d2;
4498 real inv_ncg, pos_d;
4501 npbcdim = dd->npbcdim;
4503 for (cg = cg_start; cg < cg_end; cg++)
4510 copy_rvec(state->x[k0], cm_new);
4517 for (k = k0; (k < k1); k++)
4519 rvec_inc(cm_new, state->x[k]);
4521 for (d = 0; (d < DIM); d++)
4523 cm_new[d] = inv_ncg*cm_new[d];
4528 /* Do pbc and check DD cell boundary crossings */
4529 for (d = DIM-1; d >= 0; d--)
4533 bScrew = (dd->bScrewPBC && d == XX);
4534 /* Determine the location of this cg in lattice coordinates */
4538 for (d2 = d+1; d2 < DIM; d2++)
4540 pos_d += cm_new[d2]*tcm[d2][d];
4543 /* Put the charge group in the triclinic unit-cell */
4544 if (pos_d >= cell_x1[d])
4546 if (pos_d >= limit1[d])
4548 cg_move_error(fplog, dd, step, cg, d, 1,
4549 cg_cm != state->x, limitd[d],
4550 cg_cm[cg], cm_new, pos_d);
4553 if (dd->ci[d] == dd->nc[d] - 1)
4555 rvec_dec(cm_new, state->box[d]);
4558 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4559 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4561 for (k = k0; (k < k1); k++)
4563 rvec_dec(state->x[k], state->box[d]);
4566 rotate_state_atom(state, k);
4571 else if (pos_d < cell_x0[d])
4573 if (pos_d < limit0[d])
4575 cg_move_error(fplog, dd, step, cg, d, -1,
4576 cg_cm != state->x, limitd[d],
4577 cg_cm[cg], cm_new, pos_d);
4582 rvec_inc(cm_new, state->box[d]);
4585 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4586 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4588 for (k = k0; (k < k1); k++)
4590 rvec_inc(state->x[k], state->box[d]);
4593 rotate_state_atom(state, k);
4599 else if (d < npbcdim)
4601 /* Put the charge group in the rectangular unit-cell */
4602 while (cm_new[d] >= state->box[d][d])
4604 rvec_dec(cm_new, state->box[d]);
4605 for (k = k0; (k < k1); k++)
4607 rvec_dec(state->x[k], state->box[d]);
4610 while (cm_new[d] < 0)
4612 rvec_inc(cm_new, state->box[d]);
4613 for (k = k0; (k < k1); k++)
4615 rvec_inc(state->x[k], state->box[d]);
4621 copy_rvec(cm_new, cg_cm[cg]);
4623 /* Determine where this cg should go */
4626 for (d = 0; d < dd->ndim; d++)
4631 flag |= DD_FLAG_FW(d);
4637 else if (dev[dim] == -1)
4639 flag |= DD_FLAG_BW(d);
4642 if (dd->nc[dim] > 2)
4653 /* Temporarily store the flag in move */
4654 move[cg] = mc + flag;
4658 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4659 gmx_domdec_t *dd, ivec tric_dir,
4660 t_state *state, rvec **f,
4669 int ncg[DIM*2], nat[DIM*2];
4670 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4671 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4672 int sbuf[2], rbuf[2];
4673 int home_pos_cg, home_pos_at, buf_pos;
4675 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4678 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4680 cginfo_mb_t *cginfo_mb;
4681 gmx_domdec_comm_t *comm;
4683 int nthread, thread;
4687 check_screw_box(state->box);
4691 if (fr->cutoff_scheme == ecutsGROUP)
4696 for (i = 0; i < estNR; i++)
4702 case estX: /* Always present */ break;
4703 case estV: bV = (state->flags & (1<<i)); break;
4704 case estSDX: bSDX = (state->flags & (1<<i)); break;
4705 case estCGP: bCGP = (state->flags & (1<<i)); break;
4708 case estDISRE_INITF:
4709 case estDISRE_RM3TAV:
4710 case estORIRE_INITF:
4712 /* No processing required */
4715 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4720 if (dd->ncg_tot > comm->nalloc_int)
4722 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4723 srenew(comm->buf_int, comm->nalloc_int);
4725 move = comm->buf_int;
4727 /* Clear the count */
4728 for (c = 0; c < dd->ndim*2; c++)
4734 npbcdim = dd->npbcdim;
4736 for (d = 0; (d < DIM); d++)
4738 limitd[d] = dd->comm->cellsize_min[d];
4739 if (d >= npbcdim && dd->ci[d] == 0)
4741 cell_x0[d] = -GMX_FLOAT_MAX;
4745 cell_x0[d] = comm->cell_x0[d];
4747 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4749 cell_x1[d] = GMX_FLOAT_MAX;
4753 cell_x1[d] = comm->cell_x1[d];
4757 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4758 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4762 /* We check after communication if a charge group moved
4763 * more than one cell. Set the pre-comm check limit to float_max.
4765 limit0[d] = -GMX_FLOAT_MAX;
4766 limit1[d] = GMX_FLOAT_MAX;
4770 make_tric_corr_matrix(npbcdim, state->box, tcm);
4772 cgindex = dd->cgindex;
4774 nthread = gmx_omp_nthreads_get(emntDomdec);
4776 /* Compute the center of geometry for all home charge groups
4777 * and put them in the box and determine where they should go.
4779 #pragma omp parallel for num_threads(nthread) schedule(static)
4780 for (thread = 0; thread < nthread; thread++)
4782 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4783 cell_x0, cell_x1, limitd, limit0, limit1,
4785 ( thread *dd->ncg_home)/nthread,
4786 ((thread+1)*dd->ncg_home)/nthread,
4787 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4791 for (cg = 0; cg < dd->ncg_home; cg++)
4796 flag = mc & ~DD_FLAG_NRCG;
4797 mc = mc & DD_FLAG_NRCG;
4800 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4802 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4803 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4805 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4806 /* We store the cg size in the lower 16 bits
4807 * and the place where the charge group should go
4808 * in the next 6 bits. This saves some communication volume.
4810 nrcg = cgindex[cg+1] - cgindex[cg];
4811 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4817 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4818 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4821 for (i = 0; i < dd->ndim*2; i++)
4823 *ncg_moved += ncg[i];
4840 /* Make sure the communication buffers are large enough */
4841 for (mc = 0; mc < dd->ndim*2; mc++)
4843 nvr = ncg[mc] + nat[mc]*nvec;
4844 if (nvr > comm->cgcm_state_nalloc[mc])
4846 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4847 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4851 switch (fr->cutoff_scheme)
4854 /* Recalculating cg_cm might be cheaper than communicating,
4855 * but that could give rise to rounding issues.
4858 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4859 nvec, cg_cm, comm, bCompact);
4862 /* Without charge groups we send the moved atom coordinates
4863 * over twice. This is so the code below can be used without
4864 * many conditionals for both for with and without charge groups.
4867 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4868 nvec, state->x, comm, FALSE);
4871 home_pos_cg -= *ncg_moved;
4875 gmx_incons("unimplemented");
4881 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4882 nvec, vec++, state->x, comm, bCompact);
4885 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4886 nvec, vec++, state->v, comm, bCompact);
4890 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4891 nvec, vec++, state->sd_X, comm, bCompact);
4895 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4896 nvec, vec++, state->cg_p, comm, bCompact);
4901 compact_ind(dd->ncg_home, move,
4902 dd->index_gl, dd->cgindex, dd->gatindex,
4903 dd->ga2la, comm->bLocalCG,
4908 if (fr->cutoff_scheme == ecutsVERLET)
4910 moved = get_moved(comm, dd->ncg_home);
4912 for (k = 0; k < dd->ncg_home; k++)
4919 moved = fr->ns.grid->cell_index;
4922 clear_and_mark_ind(dd->ncg_home, move,
4923 dd->index_gl, dd->cgindex, dd->gatindex,
4924 dd->ga2la, comm->bLocalCG,
4928 cginfo_mb = fr->cginfo_mb;
4930 *ncg_stay_home = home_pos_cg;
4931 for (d = 0; d < dd->ndim; d++)
4936 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4939 /* Communicate the cg and atom counts */
4944 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4945 d, dir, sbuf[0], sbuf[1]);
4947 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4949 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4951 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4952 srenew(comm->buf_int, comm->nalloc_int);
4955 /* Communicate the charge group indices, sizes and flags */
4956 dd_sendrecv_int(dd, d, dir,
4957 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4958 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4960 nvs = ncg[cdd] + nat[cdd]*nvec;
4961 i = rbuf[0] + rbuf[1] *nvec;
4962 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4964 /* Communicate cgcm and state */
4965 dd_sendrecv_rvec(dd, d, dir,
4966 comm->cgcm_state[cdd], nvs,
4967 comm->vbuf.v+nvr, i);
4968 ncg_recv += rbuf[0];
4972 /* Process the received charge groups */
4974 for (cg = 0; cg < ncg_recv; cg++)
4976 flag = comm->buf_int[cg*DD_CGIBS+1];
4978 if (dim >= npbcdim && dd->nc[dim] > 2)
4980 /* No pbc in this dim and more than one domain boundary.
4981 * We do a separate check if a charge group didn't move too far.
4983 if (((flag & DD_FLAG_FW(d)) &&
4984 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4985 ((flag & DD_FLAG_BW(d)) &&
4986 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4988 cg_move_error(fplog, dd, step, cg, dim,
4989 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4990 fr->cutoff_scheme == ecutsGROUP, 0,
4991 comm->vbuf.v[buf_pos],
4992 comm->vbuf.v[buf_pos],
4993 comm->vbuf.v[buf_pos][dim]);
5000 /* Check which direction this cg should go */
5001 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
5005 /* The cell boundaries for dimension d2 are not equal
5006 * for each cell row of the lower dimension(s),
5007 * therefore we might need to redetermine where
5008 * this cg should go.
5011 /* If this cg crosses the box boundary in dimension d2
5012 * we can use the communicated flag, so we do not
5013 * have to worry about pbc.
5015 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5016 (flag & DD_FLAG_FW(d2))) ||
5017 (dd->ci[dim2] == 0 &&
5018 (flag & DD_FLAG_BW(d2)))))
5020 /* Clear the two flags for this dimension */
5021 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5022 /* Determine the location of this cg
5023 * in lattice coordinates
5025 pos_d = comm->vbuf.v[buf_pos][dim2];
5028 for (d3 = dim2+1; d3 < DIM; d3++)
5031 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5034 /* Check of we are not at the box edge.
5035 * pbc is only handled in the first step above,
5036 * but this check could move over pbc while
5037 * the first step did not due to different rounding.
5039 if (pos_d >= cell_x1[dim2] &&
5040 dd->ci[dim2] != dd->nc[dim2]-1)
5042 flag |= DD_FLAG_FW(d2);
5044 else if (pos_d < cell_x0[dim2] &&
5047 flag |= DD_FLAG_BW(d2);
5049 comm->buf_int[cg*DD_CGIBS+1] = flag;
5052 /* Set to which neighboring cell this cg should go */
5053 if (flag & DD_FLAG_FW(d2))
5057 else if (flag & DD_FLAG_BW(d2))
5059 if (dd->nc[dd->dim[d2]] > 2)
5071 nrcg = flag & DD_FLAG_NRCG;
5074 if (home_pos_cg+1 > dd->cg_nalloc)
5076 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5077 srenew(dd->index_gl, dd->cg_nalloc);
5078 srenew(dd->cgindex, dd->cg_nalloc+1);
5080 /* Set the global charge group index and size */
5081 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5082 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5083 /* Copy the state from the buffer */
5084 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5085 if (fr->cutoff_scheme == ecutsGROUP)
5088 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5092 /* Set the cginfo */
5093 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5094 dd->index_gl[home_pos_cg]);
5097 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5100 if (home_pos_at+nrcg > state->nalloc)
5102 dd_realloc_state(state, f, home_pos_at+nrcg);
5104 for (i = 0; i < nrcg; i++)
5106 copy_rvec(comm->vbuf.v[buf_pos++],
5107 state->x[home_pos_at+i]);
5111 for (i = 0; i < nrcg; i++)
5113 copy_rvec(comm->vbuf.v[buf_pos++],
5114 state->v[home_pos_at+i]);
5119 for (i = 0; i < nrcg; i++)
5121 copy_rvec(comm->vbuf.v[buf_pos++],
5122 state->sd_X[home_pos_at+i]);
5127 for (i = 0; i < nrcg; i++)
5129 copy_rvec(comm->vbuf.v[buf_pos++],
5130 state->cg_p[home_pos_at+i]);
5134 home_pos_at += nrcg;
5138 /* Reallocate the buffers if necessary */
5139 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5141 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5142 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5144 nvr = ncg[mc] + nat[mc]*nvec;
5145 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5147 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5148 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5150 /* Copy from the receive to the send buffers */
5151 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5152 comm->buf_int + cg*DD_CGIBS,
5153 DD_CGIBS*sizeof(int));
5154 memcpy(comm->cgcm_state[mc][nvr],
5155 comm->vbuf.v[buf_pos],
5156 (1+nrcg*nvec)*sizeof(rvec));
5157 buf_pos += 1 + nrcg*nvec;
5164 /* With sorting (!bCompact) the indices are now only partially up to date
5165 * and ncg_home and nat_home are not the real count, since there are
5166 * "holes" in the arrays for the charge groups that moved to neighbors.
5168 if (fr->cutoff_scheme == ecutsVERLET)
5170 moved = get_moved(comm, home_pos_cg);
5172 for (i = dd->ncg_home; i < home_pos_cg; i++)
5177 dd->ncg_home = home_pos_cg;
5178 dd->nat_home = home_pos_at;
5183 "Finished repartitioning: cgs moved out %d, new home %d\n",
5184 *ncg_moved, dd->ncg_home-*ncg_moved);
5189 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5191 dd->comm->cycl[ddCycl] += cycles;
5192 dd->comm->cycl_n[ddCycl]++;
5193 if (cycles > dd->comm->cycl_max[ddCycl])
5195 dd->comm->cycl_max[ddCycl] = cycles;
5199 static double force_flop_count(t_nrnb *nrnb)
5206 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5208 /* To get closer to the real timings, we half the count
5209 * for the normal loops and again half it for water loops.
5212 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5214 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5218 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5221 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5224 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5226 sum += nrnb->n[i]*cost_nrnb(i);
5229 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5231 sum += nrnb->n[i]*cost_nrnb(i);
5237 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5239 if (dd->comm->eFlop)
5241 dd->comm->flop -= force_flop_count(nrnb);
5244 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5246 if (dd->comm->eFlop)
5248 dd->comm->flop += force_flop_count(nrnb);
5253 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5257 for (i = 0; i < ddCyclNr; i++)
5259 dd->comm->cycl[i] = 0;
5260 dd->comm->cycl_n[i] = 0;
5261 dd->comm->cycl_max[i] = 0;
5264 dd->comm->flop_n = 0;
5267 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5269 gmx_domdec_comm_t *comm;
5270 gmx_domdec_load_t *load;
5271 gmx_domdec_root_t *root = NULL;
5273 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5278 fprintf(debug, "get_load_distribution start\n");
5281 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5285 bSepPME = (dd->pme_nodeid >= 0);
5287 if (dd->ndim == 0 && bSepPME)
5289 /* Without decomposition, but with PME nodes, we need the load */
5290 comm->load[0].mdf = comm->cycl[ddCyclPPduringPME];
5291 comm->load[0].pme = comm->cycl[ddCyclPME];
5294 for (d = dd->ndim-1; d >= 0; d--)
5297 /* Check if we participate in the communication in this dimension */
5298 if (d == dd->ndim-1 ||
5299 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5301 load = &comm->load[d];
5304 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5307 if (d == dd->ndim-1)
5309 sbuf[pos++] = dd_force_load(comm);
5310 sbuf[pos++] = sbuf[0];
5313 sbuf[pos++] = sbuf[0];
5314 sbuf[pos++] = cell_frac;
5317 sbuf[pos++] = comm->cell_f_max0[d];
5318 sbuf[pos++] = comm->cell_f_min1[d];
5323 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5324 sbuf[pos++] = comm->cycl[ddCyclPME];
5329 sbuf[pos++] = comm->load[d+1].sum;
5330 sbuf[pos++] = comm->load[d+1].max;
5333 sbuf[pos++] = comm->load[d+1].sum_m;
5334 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5335 sbuf[pos++] = comm->load[d+1].flags;
5338 sbuf[pos++] = comm->cell_f_max0[d];
5339 sbuf[pos++] = comm->cell_f_min1[d];
5344 sbuf[pos++] = comm->load[d+1].mdf;
5345 sbuf[pos++] = comm->load[d+1].pme;
5349 /* Communicate a row in DD direction d.
5350 * The communicators are setup such that the root always has rank 0.
5353 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5354 load->load, load->nload*sizeof(float), MPI_BYTE,
5355 0, comm->mpi_comm_load[d]);
5357 if (dd->ci[dim] == dd->master_ci[dim])
5359 /* We are the root, process this row */
5360 if (comm->bDynLoadBal)
5362 root = comm->root[d];
5372 for (i = 0; i < dd->nc[dim]; i++)
5374 load->sum += load->load[pos++];
5375 load->max = std::max(load->max, load->load[pos]);
5381 /* This direction could not be load balanced properly,
5382 * therefore we need to use the maximum iso the average load.
5384 load->sum_m = std::max(load->sum_m, load->load[pos]);
5388 load->sum_m += load->load[pos];
5391 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5395 load->flags = (int)(load->load[pos++] + 0.5);
5399 root->cell_f_max0[i] = load->load[pos++];
5400 root->cell_f_min1[i] = load->load[pos++];
5405 load->mdf = std::max(load->mdf, load->load[pos]);
5407 load->pme = std::max(load->pme, load->load[pos]);
5411 if (comm->bDynLoadBal && root->bLimited)
5413 load->sum_m *= dd->nc[dim];
5414 load->flags |= (1<<d);
5422 comm->nload += dd_load_count(comm);
5423 comm->load_step += comm->cycl[ddCyclStep];
5424 comm->load_sum += comm->load[0].sum;
5425 comm->load_max += comm->load[0].max;
5426 if (comm->bDynLoadBal)
5428 for (d = 0; d < dd->ndim; d++)
5430 if (comm->load[0].flags & (1<<d))
5432 comm->load_lim[d]++;
5438 comm->load_mdf += comm->load[0].mdf;
5439 comm->load_pme += comm->load[0].pme;
5443 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5447 fprintf(debug, "get_load_distribution finished\n");
5451 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5453 /* Return the relative performance loss on the total run time
5454 * due to the force calculation load imbalance.
5456 if (dd->comm->nload > 0)
5459 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5460 (dd->comm->load_step*dd->nnodes);
5468 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5471 int npp, npme, nnodes, d, limp;
5472 float imbal, pme_f_ratio, lossf = 0, lossp = 0;
5474 gmx_domdec_comm_t *comm;
5477 if (DDMASTER(dd) && comm->nload > 0)
5480 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5481 nnodes = npp + npme;
5484 imbal = comm->load_max*npp/comm->load_sum - 1;
5485 lossf = dd_force_imb_perf_loss(dd);
5486 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5487 fprintf(fplog, "%s", buf);
5488 fprintf(stderr, "\n");
5489 fprintf(stderr, "%s", buf);
5490 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5491 fprintf(fplog, "%s", buf);
5492 fprintf(stderr, "%s", buf);
5495 if (comm->bDynLoadBal)
5497 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5498 for (d = 0; d < dd->ndim; d++)
5500 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5501 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5507 sprintf(buf+strlen(buf), "\n");
5508 fprintf(fplog, "%s", buf);
5509 fprintf(stderr, "%s", buf);
5513 pme_f_ratio = comm->load_pme/comm->load_mdf;
5514 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5517 lossp *= (float)npme/(float)nnodes;
5521 lossp *= (float)npp/(float)nnodes;
5523 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5524 fprintf(fplog, "%s", buf);
5525 fprintf(stderr, "%s", buf);
5526 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5527 fprintf(fplog, "%s", buf);
5528 fprintf(stderr, "%s", buf);
5530 fprintf(fplog, "\n");
5531 fprintf(stderr, "\n");
5533 if (lossf >= DD_PERF_LOSS_WARN)
5536 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5537 " in the domain decomposition.\n", lossf*100);
5538 if (!comm->bDynLoadBal)
5540 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5544 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5546 fprintf(fplog, "%s\n", buf);
5547 fprintf(stderr, "%s\n", buf);
5549 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5552 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5553 " had %s work to do than the PP ranks.\n"
5554 " You might want to %s the number of PME ranks\n"
5555 " or %s the cut-off and the grid spacing.\n",
5557 (lossp < 0) ? "less" : "more",
5558 (lossp < 0) ? "decrease" : "increase",
5559 (lossp < 0) ? "decrease" : "increase");
5560 fprintf(fplog, "%s\n", buf);
5561 fprintf(stderr, "%s\n", buf);
5566 static float dd_vol_min(gmx_domdec_t *dd)
5568 return dd->comm->load[0].cvol_min*dd->nnodes;
5571 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5573 return dd->comm->load[0].flags;
5576 static float dd_f_imbal(gmx_domdec_t *dd)
5578 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5581 float dd_pme_f_ratio(gmx_domdec_t *dd)
5583 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5585 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5593 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5598 flags = dd_load_flags(dd);
5602 "DD load balancing is limited by minimum cell size in dimension");
5603 for (d = 0; d < dd->ndim; d++)
5607 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5610 fprintf(fplog, "\n");
5612 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5613 if (dd->comm->bDynLoadBal)
5615 fprintf(fplog, " vol min/aver %5.3f%c",
5616 dd_vol_min(dd), flags ? '!' : ' ');
5620 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5622 if (dd->comm->cycl_n[ddCyclPME])
5624 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5626 fprintf(fplog, "\n\n");
5629 static void dd_print_load_verbose(gmx_domdec_t *dd)
5631 if (dd->comm->bDynLoadBal)
5633 fprintf(stderr, "vol %4.2f%c ",
5634 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5638 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5640 if (dd->comm->cycl_n[ddCyclPME])
5642 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5647 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5652 gmx_domdec_root_t *root;
5653 gmx_bool bPartOfGroup = FALSE;
5655 dim = dd->dim[dim_ind];
5656 copy_ivec(loc, loc_c);
5657 for (i = 0; i < dd->nc[dim]; i++)
5660 rank = dd_index(dd->nc, loc_c);
5661 if (rank == dd->rank)
5663 /* This process is part of the group */
5664 bPartOfGroup = TRUE;
5667 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5671 dd->comm->mpi_comm_load[dim_ind] = c_row;
5672 if (dd->comm->eDLB != edlbNO)
5674 if (dd->ci[dim] == dd->master_ci[dim])
5676 /* This is the root process of this row */
5677 snew(dd->comm->root[dim_ind], 1);
5678 root = dd->comm->root[dim_ind];
5679 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5680 snew(root->old_cell_f, dd->nc[dim]+1);
5681 snew(root->bCellMin, dd->nc[dim]);
5684 snew(root->cell_f_max0, dd->nc[dim]);
5685 snew(root->cell_f_min1, dd->nc[dim]);
5686 snew(root->bound_min, dd->nc[dim]);
5687 snew(root->bound_max, dd->nc[dim]);
5689 snew(root->buf_ncd, dd->nc[dim]);
5693 /* This is not a root process, we only need to receive cell_f */
5694 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5697 if (dd->ci[dim] == dd->master_ci[dim])
5699 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5705 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5706 const gmx_hw_info_t gmx_unused *hwinfo,
5707 const gmx_hw_opt_t gmx_unused *hw_opt)
5710 int physicalnode_id_hash;
5713 MPI_Comm mpi_comm_pp_physicalnode;
5715 if (!(cr->duty & DUTY_PP) || hw_opt->gpu_opt.n_dev_use == 0)
5717 /* Only PP nodes (currently) use GPUs.
5718 * If we don't have GPUs, there are no resources to share.
5723 physicalnode_id_hash = gmx_physicalnode_id_hash();
5725 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5731 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5732 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5733 dd->rank, physicalnode_id_hash, gpu_id);
5735 /* Split the PP communicator over the physical nodes */
5736 /* TODO: See if we should store this (before), as it's also used for
5737 * for the nodecomm summution.
5739 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5740 &mpi_comm_pp_physicalnode);
5741 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5742 &dd->comm->mpi_comm_gpu_shared);
5743 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5744 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5748 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5751 /* Note that some ranks could share a GPU, while others don't */
5753 if (dd->comm->nrank_gpu_shared == 1)
5755 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5760 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5763 int dim0, dim1, i, j;
5768 fprintf(debug, "Making load communicators\n");
5771 snew(dd->comm->load, std::max(dd->ndim, 1));
5772 snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
5780 make_load_communicator(dd, 0, loc);
5784 for (i = 0; i < dd->nc[dim0]; i++)
5787 make_load_communicator(dd, 1, loc);
5793 for (i = 0; i < dd->nc[dim0]; i++)
5797 for (j = 0; j < dd->nc[dim1]; j++)
5800 make_load_communicator(dd, 2, loc);
5807 fprintf(debug, "Finished making load communicators\n");
5812 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5814 int d, dim, i, j, m;
5817 ivec dd_zp[DD_MAXIZONE];
5818 gmx_domdec_zones_t *zones;
5819 gmx_domdec_ns_ranges_t *izone;
5821 for (d = 0; d < dd->ndim; d++)
5824 copy_ivec(dd->ci, tmp);
5825 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5826 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5827 copy_ivec(dd->ci, tmp);
5828 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5829 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5832 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5835 dd->neighbor[d][1]);
5841 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5843 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5844 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5851 for (i = 0; i < nzonep; i++)
5853 copy_ivec(dd_zp3[i], dd_zp[i]);
5859 for (i = 0; i < nzonep; i++)
5861 copy_ivec(dd_zp2[i], dd_zp[i]);
5867 for (i = 0; i < nzonep; i++)
5869 copy_ivec(dd_zp1[i], dd_zp[i]);
5875 for (i = 0; i < nzonep; i++)
5877 copy_ivec(dd_zp0[i], dd_zp[i]);
5881 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5886 zones = &dd->comm->zones;
5888 for (i = 0; i < nzone; i++)
5891 clear_ivec(zones->shift[i]);
5892 for (d = 0; d < dd->ndim; d++)
5894 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5899 for (i = 0; i < nzone; i++)
5901 for (d = 0; d < DIM; d++)
5903 s[d] = dd->ci[d] - zones->shift[i][d];
5908 else if (s[d] >= dd->nc[d])
5914 zones->nizone = nzonep;
5915 for (i = 0; i < zones->nizone; i++)
5917 if (dd_zp[i][0] != i)
5919 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5921 izone = &zones->izone[i];
5922 izone->j0 = dd_zp[i][1];
5923 izone->j1 = dd_zp[i][2];
5924 for (dim = 0; dim < DIM; dim++)
5926 if (dd->nc[dim] == 1)
5928 /* All shifts should be allowed */
5929 izone->shift0[dim] = -1;
5930 izone->shift1[dim] = 1;
5935 izone->shift0[d] = 0;
5936 izone->shift1[d] = 0;
5937 for(j=izone->j0; j<izone->j1; j++) {
5938 if (dd->shift[j][d] > dd->shift[i][d])
5939 izone->shift0[d] = -1;
5940 if (dd->shift[j][d] < dd->shift[i][d])
5941 izone->shift1[d] = 1;
5947 /* Assume the shift are not more than 1 cell */
5948 izone->shift0[dim] = 1;
5949 izone->shift1[dim] = -1;
5950 for (j = izone->j0; j < izone->j1; j++)
5952 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5953 if (shift_diff < izone->shift0[dim])
5955 izone->shift0[dim] = shift_diff;
5957 if (shift_diff > izone->shift1[dim])
5959 izone->shift1[dim] = shift_diff;
5966 if (dd->comm->eDLB != edlbNO)
5968 snew(dd->comm->root, dd->ndim);
5971 if (dd->comm->bRecordLoad)
5973 make_load_communicators(dd);
5977 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5983 gmx_domdec_comm_t *comm;
5990 if (comm->bCartesianPP)
5992 /* Set up cartesian communication for the particle-particle part */
5995 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5996 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5999 for (int i = 0; i < DIM; i++)
6003 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
6005 /* We overwrite the old communicator with the new cartesian one */
6006 cr->mpi_comm_mygroup = comm_cart;
6009 dd->mpi_comm_all = cr->mpi_comm_mygroup;
6010 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
6012 if (comm->bCartesianPP_PME)
6014 /* Since we want to use the original cartesian setup for sim,
6015 * and not the one after split, we need to make an index.
6017 snew(comm->ddindex2ddnodeid, dd->nnodes);
6018 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
6019 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
6020 /* Get the rank of the DD master,
6021 * above we made sure that the master node is a PP node.
6031 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
6033 else if (comm->bCartesianPP)
6035 if (cr->npmenodes == 0)
6037 /* The PP communicator is also
6038 * the communicator for this simulation
6040 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6042 cr->nodeid = dd->rank;
6044 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6046 /* We need to make an index to go from the coordinates
6047 * to the nodeid of this simulation.
6049 snew(comm->ddindex2simnodeid, dd->nnodes);
6050 snew(buf, dd->nnodes);
6051 if (cr->duty & DUTY_PP)
6053 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6055 /* Communicate the ddindex to simulation nodeid index */
6056 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6057 cr->mpi_comm_mysim);
6060 /* Determine the master coordinates and rank.
6061 * The DD master should be the same node as the master of this sim.
6063 for (int i = 0; i < dd->nnodes; i++)
6065 if (comm->ddindex2simnodeid[i] == 0)
6067 ddindex2xyz(dd->nc, i, dd->master_ci);
6068 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6073 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6078 /* No Cartesian communicators */
6079 /* We use the rank in dd->comm->all as DD index */
6080 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6081 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6083 clear_ivec(dd->master_ci);
6090 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6091 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6096 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6097 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6101 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6105 gmx_domdec_comm_t *comm;
6110 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6113 snew(comm->ddindex2simnodeid, dd->nnodes);
6114 snew(buf, dd->nnodes);
6115 if (cr->duty & DUTY_PP)
6117 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6119 /* Communicate the ddindex to simulation nodeid index */
6120 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6121 cr->mpi_comm_mysim);
6127 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6128 int ncg, int natoms)
6130 gmx_domdec_master_t *ma;
6135 snew(ma->ncg, dd->nnodes);
6136 snew(ma->index, dd->nnodes+1);
6138 snew(ma->nat, dd->nnodes);
6139 snew(ma->ibuf, dd->nnodes*2);
6140 snew(ma->cell_x, DIM);
6141 for (i = 0; i < DIM; i++)
6143 snew(ma->cell_x[i], dd->nc[i]+1);
6146 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6152 snew(ma->vbuf, natoms);
6158 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6159 int gmx_unused reorder)
6162 gmx_domdec_comm_t *comm;
6172 if (comm->bCartesianPP)
6174 for (i = 1; i < DIM; i++)
6176 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6178 if (bDiv[YY] || bDiv[ZZ])
6180 comm->bCartesianPP_PME = TRUE;
6181 /* If we have 2D PME decomposition, which is always in x+y,
6182 * we stack the PME only nodes in z.
6183 * Otherwise we choose the direction that provides the thinnest slab
6184 * of PME only nodes as this will have the least effect
6185 * on the PP communication.
6186 * But for the PME communication the opposite might be better.
6188 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6190 dd->nc[YY] > dd->nc[ZZ]))
6192 comm->cartpmedim = ZZ;
6196 comm->cartpmedim = YY;
6198 comm->ntot[comm->cartpmedim]
6199 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6203 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]);
6205 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6210 if (comm->bCartesianPP_PME)
6217 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]);
6220 for (i = 0; i < DIM; i++)
6224 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6226 MPI_Comm_rank(comm_cart, &rank);
6227 if (MASTERNODE(cr) && rank != 0)
6229 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6232 /* With this assigment we loose the link to the original communicator
6233 * which will usually be MPI_COMM_WORLD, unless have multisim.
6235 cr->mpi_comm_mysim = comm_cart;
6236 cr->sim_nodeid = rank;
6238 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6242 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6243 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6246 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6250 if (cr->npmenodes == 0 ||
6251 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6253 cr->duty = DUTY_PME;
6256 /* Split the sim communicator into PP and PME only nodes */
6257 MPI_Comm_split(cr->mpi_comm_mysim,
6259 dd_index(comm->ntot, dd->ci),
6260 &cr->mpi_comm_mygroup);
6264 switch (dd_node_order)
6269 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6272 case ddnoINTERLEAVE:
6273 /* Interleave the PP-only and PME-only nodes,
6274 * as on clusters with dual-core machines this will double
6275 * the communication bandwidth of the PME processes
6276 * and thus speed up the PP <-> PME and inter PME communication.
6280 fprintf(fplog, "Interleaving PP and PME ranks\n");
6282 comm->pmenodes = dd_pmenodes(cr);
6287 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6290 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6292 cr->duty = DUTY_PME;
6299 /* Split the sim communicator into PP and PME only nodes */
6300 MPI_Comm_split(cr->mpi_comm_mysim,
6303 &cr->mpi_comm_mygroup);
6304 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6310 fprintf(fplog, "This rank does only %s work.\n\n",
6311 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6315 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6318 gmx_domdec_comm_t *comm;
6324 copy_ivec(dd->nc, comm->ntot);
6326 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6327 comm->bCartesianPP_PME = FALSE;
6329 /* Reorder the nodes by default. This might change the MPI ranks.
6330 * Real reordering is only supported on very few architectures,
6331 * Blue Gene is one of them.
6333 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6335 if (cr->npmenodes > 0)
6337 /* Split the communicator into a PP and PME part */
6338 split_communicator(fplog, cr, dd_node_order, CartReorder);
6339 if (comm->bCartesianPP_PME)
6341 /* We (possibly) reordered the nodes in split_communicator,
6342 * so it is no longer required in make_pp_communicator.
6344 CartReorder = FALSE;
6349 /* All nodes do PP and PME */
6351 /* We do not require separate communicators */
6352 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6356 if (cr->duty & DUTY_PP)
6358 /* Copy or make a new PP communicator */
6359 make_pp_communicator(fplog, cr, CartReorder);
6363 receive_ddindex2simnodeid(cr);
6366 if (!(cr->duty & DUTY_PME))
6368 /* Set up the commnuication to our PME node */
6369 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6370 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6373 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6374 dd->pme_nodeid, dd->pme_receive_vir_ener);
6379 dd->pme_nodeid = -1;
6384 dd->ma = init_gmx_domdec_master_t(dd,
6386 comm->cgs_gl.index[comm->cgs_gl.nr]);
6390 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6392 real *slb_frac, tot;
6397 if (nc > 1 && size_string != NULL)
6401 fprintf(fplog, "Using static load balancing for the %s direction\n",
6406 for (i = 0; i < nc; i++)
6409 sscanf(size_string, "%20lf%n", &dbl, &n);
6412 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6421 fprintf(fplog, "Relative cell sizes:");
6423 for (i = 0; i < nc; i++)
6428 fprintf(fplog, " %5.3f", slb_frac[i]);
6433 fprintf(fplog, "\n");
6440 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6443 gmx_mtop_ilistloop_t iloop;
6447 iloop = gmx_mtop_ilistloop_init(mtop);
6448 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6450 for (ftype = 0; ftype < F_NRE; ftype++)
6452 if ((interaction_function[ftype].flags & IF_BOND) &&
6455 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6463 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6469 val = getenv(env_var);
6472 if (sscanf(val, "%20d", &nst) <= 0)
6478 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6486 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6490 fprintf(stderr, "\n%s\n", warn_string);
6494 fprintf(fplog, "\n%s\n", warn_string);
6498 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6499 t_inputrec *ir, FILE *fplog)
6501 if (ir->ePBC == epbcSCREW &&
6502 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6504 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6507 if (ir->ns_type == ensSIMPLE)
6509 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6512 if (ir->nstlist == 0)
6514 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6517 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6519 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6523 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6528 r = ddbox->box_size[XX];
6529 for (di = 0; di < dd->ndim; di++)
6532 /* Check using the initial average cell size */
6533 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6539 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6540 const char *dlb_opt, gmx_bool bRecordLoad,
6541 unsigned long Flags, t_inputrec *ir)
6548 case 'a': eDLB = edlbAUTO; break;
6549 case 'n': eDLB = edlbNO; break;
6550 case 'y': eDLB = edlbYES; break;
6551 default: gmx_incons("Unknown dlb_opt");
6554 if (Flags & MD_RERUN)
6559 if (!EI_DYNAMICS(ir->eI))
6561 if (eDLB == edlbYES)
6563 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6564 dd_warning(cr, fplog, buf);
6572 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6577 if (Flags & MD_REPRODUCIBLE)
6584 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6588 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6591 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6599 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6604 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6606 /* Decomposition order z,y,x */
6609 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6611 for (dim = DIM-1; dim >= 0; dim--)
6613 if (dd->nc[dim] > 1)
6615 dd->dim[dd->ndim++] = dim;
6621 /* Decomposition order x,y,z */
6622 for (dim = 0; dim < DIM; dim++)
6624 if (dd->nc[dim] > 1)
6626 dd->dim[dd->ndim++] = dim;
6632 static gmx_domdec_comm_t *init_dd_comm()
6634 gmx_domdec_comm_t *comm;
6638 snew(comm->cggl_flag, DIM*2);
6639 snew(comm->cgcm_state, DIM*2);
6640 for (i = 0; i < DIM*2; i++)
6642 comm->cggl_flag_nalloc[i] = 0;
6643 comm->cgcm_state_nalloc[i] = 0;
6646 comm->nalloc_int = 0;
6647 comm->buf_int = NULL;
6649 vec_rvec_init(&comm->vbuf);
6651 comm->n_load_have = 0;
6652 comm->n_load_collect = 0;
6654 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6656 comm->sum_nat[i] = 0;
6660 comm->load_step = 0;
6663 clear_ivec(comm->load_lim);
6670 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6671 unsigned long Flags,
6673 real comm_distance_min, real rconstr,
6674 const char *dlb_opt, real dlb_scale,
6675 const char *sizex, const char *sizey, const char *sizez,
6676 gmx_mtop_t *mtop, t_inputrec *ir,
6677 matrix box, rvec *x,
6679 int *npme_x, int *npme_y)
6682 gmx_domdec_comm_t *comm;
6684 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6687 const real tenPercentMargin = 1.1;
6692 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6697 dd->comm = init_dd_comm();
6699 snew(comm->cggl_flag, DIM*2);
6700 snew(comm->cgcm_state, DIM*2);
6702 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6703 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6705 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6706 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6707 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6708 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6709 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6710 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6711 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6712 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6714 dd->pme_recv_f_alloc = 0;
6715 dd->pme_recv_f_buf = NULL;
6717 if (dd->bSendRecv2 && fplog)
6719 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");
6725 fprintf(fplog, "Will load balance based on FLOP count\n");
6727 if (comm->eFlop > 1)
6729 srand(1+cr->nodeid);
6731 comm->bRecordLoad = TRUE;
6735 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6739 /* Initialize to GPU share count to 0, might change later */
6740 comm->nrank_gpu_shared = 0;
6742 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6743 comm->bDLB_locked = FALSE;
6744 comm->bCheckWhetherToTurnDlbOn = TRUE;
6746 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6749 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6751 dd->bGridJump = comm->bDynLoadBal;
6752 comm->bPMELoadBalDLBLimits = FALSE;
6754 if (comm->nstSortCG)
6758 if (comm->nstSortCG == 1)
6760 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6764 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6768 snew(comm->sort, 1);
6774 fprintf(fplog, "Will not sort the charge groups\n");
6778 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6780 comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
6781 mtop->bIntermolecularInteractions);
6782 if (comm->bInterCGBondeds)
6784 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6788 comm->bInterCGMultiBody = FALSE;
6791 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6792 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6794 if (ir->rlistlong == 0)
6796 /* Set the cut-off to some very large value,
6797 * so we don't need if statements everywhere in the code.
6798 * We use sqrt, since the cut-off is squared in some places.
6800 comm->cutoff = GMX_CUTOFF_INF;
6804 comm->cutoff = ir->rlistlong;
6806 comm->cutoff_mbody = 0;
6808 comm->cellsize_limit = 0;
6809 comm->bBondComm = FALSE;
6811 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6812 * within nstlist steps. Since boundaries are allowed to displace by half
6813 * a cell size, DD cells should be at least the size of the list buffer.
6815 comm->cellsize_limit = std::max(comm->cellsize_limit,
6816 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6818 if (comm->bInterCGBondeds)
6820 if (comm_distance_min > 0)
6822 comm->cutoff_mbody = comm_distance_min;
6823 if (Flags & MD_DDBONDCOMM)
6825 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6829 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6831 r_bonded_limit = comm->cutoff_mbody;
6833 else if (ir->bPeriodicMols)
6835 /* Can not easily determine the required cut-off */
6836 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");
6837 comm->cutoff_mbody = comm->cutoff/2;
6838 r_bonded_limit = comm->cutoff_mbody;
6844 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6845 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6847 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6848 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6850 /* We use an initial margin of 10% for the minimum cell size,
6851 * except when we are just below the non-bonded cut-off.
6853 if (Flags & MD_DDBONDCOMM)
6855 if (std::max(r_2b, r_mb) > comm->cutoff)
6857 r_bonded = std::max(r_2b, r_mb);
6858 r_bonded_limit = tenPercentMargin*r_bonded;
6859 comm->bBondComm = TRUE;
6864 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6866 /* We determine cutoff_mbody later */
6870 /* No special bonded communication,
6871 * simply increase the DD cut-off.
6873 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6874 comm->cutoff_mbody = r_bonded_limit;
6875 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6881 "Minimum cell size due to bonded interactions: %.3f nm\n",
6884 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6887 if (dd->bInterCGcons && rconstr <= 0)
6889 /* There is a cell size limit due to the constraints (P-LINCS) */
6890 rconstr = constr_r_max(fplog, mtop, ir);
6894 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6896 if (rconstr > comm->cellsize_limit)
6898 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6902 else if (rconstr > 0 && fplog)
6904 /* Here we do not check for dd->bInterCGcons,
6905 * because one can also set a cell size limit for virtual sites only
6906 * and at this point we don't know yet if there are intercg v-sites.
6909 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6912 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6914 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6918 copy_ivec(nc, dd->nc);
6919 set_dd_dim(fplog, dd);
6920 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6922 if (cr->npmenodes == -1)
6926 acs = average_cellsize_min(dd, ddbox);
6927 if (acs < comm->cellsize_limit)
6931 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6933 gmx_fatal_collective(FARGS, cr, NULL,
6934 "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",
6935 acs, comm->cellsize_limit);
6940 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6942 /* We need to choose the optimal DD grid and possibly PME nodes */
6943 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6944 comm->eDLB != edlbNO, dlb_scale,
6945 comm->cellsize_limit, comm->cutoff,
6946 comm->bInterCGBondeds);
6948 if (dd->nc[XX] == 0)
6950 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6951 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6952 !bC ? "-rdd" : "-rcon",
6953 comm->eDLB != edlbNO ? " or -dds" : "",
6954 bC ? " or your LINCS settings" : "");
6956 gmx_fatal_collective(FARGS, cr, NULL,
6957 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6959 "Look in the log file for details on the domain decomposition",
6960 cr->nnodes-cr->npmenodes, limit, buf);
6962 set_dd_dim(fplog, dd);
6968 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6969 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6972 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6973 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6975 gmx_fatal_collective(FARGS, cr, NULL,
6976 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6977 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6979 if (cr->npmenodes > dd->nnodes)
6981 gmx_fatal_collective(FARGS, cr, NULL,
6982 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6984 if (cr->npmenodes > 0)
6986 comm->npmenodes = cr->npmenodes;
6990 comm->npmenodes = dd->nnodes;
6993 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6995 /* The following choices should match those
6996 * in comm_cost_est in domdec_setup.c.
6997 * Note that here the checks have to take into account
6998 * that the decomposition might occur in a different order than xyz
6999 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
7000 * in which case they will not match those in comm_cost_est,
7001 * but since that is mainly for testing purposes that's fine.
7003 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
7004 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
7005 getenv("GMX_PMEONEDD") == NULL)
7007 comm->npmedecompdim = 2;
7008 comm->npmenodes_x = dd->nc[XX];
7009 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
7013 /* In case nc is 1 in both x and y we could still choose to
7014 * decompose pme in y instead of x, but we use x for simplicity.
7016 comm->npmedecompdim = 1;
7017 if (dd->dim[0] == YY)
7019 comm->npmenodes_x = 1;
7020 comm->npmenodes_y = comm->npmenodes;
7024 comm->npmenodes_x = comm->npmenodes;
7025 comm->npmenodes_y = 1;
7030 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
7031 comm->npmenodes_x, comm->npmenodes_y, 1);
7036 comm->npmedecompdim = 0;
7037 comm->npmenodes_x = 0;
7038 comm->npmenodes_y = 0;
7041 /* Technically we don't need both of these,
7042 * but it simplifies code not having to recalculate it.
7044 *npme_x = comm->npmenodes_x;
7045 *npme_y = comm->npmenodes_y;
7047 snew(comm->slb_frac, DIM);
7048 if (comm->eDLB == edlbNO)
7050 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7051 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7052 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7055 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7057 if (comm->bBondComm || comm->eDLB != edlbNO)
7059 /* Set the bonded communication distance to halfway
7060 * the minimum and the maximum,
7061 * since the extra communication cost is nearly zero.
7063 acs = average_cellsize_min(dd, ddbox);
7064 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7065 if (comm->eDLB != edlbNO)
7067 /* Check if this does not limit the scaling */
7068 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7070 if (!comm->bBondComm)
7072 /* Without bBondComm do not go beyond the n.b. cut-off */
7073 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7074 if (comm->cellsize_limit >= comm->cutoff)
7076 /* We don't loose a lot of efficieny
7077 * when increasing it to the n.b. cut-off.
7078 * It can even be slightly faster, because we need
7079 * less checks for the communication setup.
7081 comm->cutoff_mbody = comm->cutoff;
7084 /* Check if we did not end up below our original limit */
7085 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7087 if (comm->cutoff_mbody > comm->cellsize_limit)
7089 comm->cellsize_limit = comm->cutoff_mbody;
7092 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7097 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7098 "cellsize limit %f\n",
7099 comm->bBondComm, comm->cellsize_limit);
7104 check_dd_restrictions(cr, dd, ir, fplog);
7107 comm->partition_step = INT_MIN;
7110 clear_dd_cycle_counts(dd);
7115 static void set_dlb_limits(gmx_domdec_t *dd)
7120 for (d = 0; d < dd->ndim; d++)
7122 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7123 dd->comm->cellsize_min[dd->dim[d]] =
7124 dd->comm->cellsize_min_dlb[dd->dim[d]];
7129 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7132 gmx_domdec_comm_t *comm;
7142 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);
7145 cellsize_min = comm->cellsize_min[dd->dim[0]];
7146 for (d = 1; d < dd->ndim; d++)
7148 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7151 if (cellsize_min < comm->cellsize_limit*1.05)
7153 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");
7155 /* Change DLB from "auto" to "no". */
7156 comm->eDLB = edlbNO;
7161 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7162 comm->bDynLoadBal = TRUE;
7163 dd->bGridJump = TRUE;
7167 /* We can set the required cell size info here,
7168 * so we do not need to communicate this.
7169 * The grid is completely uniform.
7171 for (d = 0; d < dd->ndim; d++)
7175 comm->load[d].sum_m = comm->load[d].sum;
7177 nc = dd->nc[dd->dim[d]];
7178 for (i = 0; i < nc; i++)
7180 comm->root[d]->cell_f[i] = i/(real)nc;
7183 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7184 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7187 comm->root[d]->cell_f[nc] = 1.0;
7192 static char *init_bLocalCG(gmx_mtop_t *mtop)
7197 ncg = ncg_mtop(mtop);
7198 snew(bLocalCG, ncg);
7199 for (cg = 0; cg < ncg; cg++)
7201 bLocalCG[cg] = FALSE;
7207 void dd_init_bondeds(FILE *fplog,
7208 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7210 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7212 gmx_domdec_comm_t *comm;
7214 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7218 if (comm->bBondComm)
7220 /* Communicate atoms beyond the cut-off for bonded interactions */
7223 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7225 comm->bLocalCG = init_bLocalCG(mtop);
7229 /* Only communicate atoms based on cut-off */
7230 comm->cglink = NULL;
7231 comm->bLocalCG = NULL;
7235 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7237 gmx_bool bDynLoadBal, real dlb_scale,
7240 gmx_domdec_comm_t *comm;
7255 fprintf(fplog, "The maximum number of communication pulses is:");
7256 for (d = 0; d < dd->ndim; d++)
7258 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7260 fprintf(fplog, "\n");
7261 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7262 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7263 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7264 for (d = 0; d < DIM; d++)
7268 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7275 comm->cellsize_min_dlb[d]/
7276 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7278 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7281 fprintf(fplog, "\n");
7285 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7286 fprintf(fplog, "The initial number of communication pulses is:");
7287 for (d = 0; d < dd->ndim; d++)
7289 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7291 fprintf(fplog, "\n");
7292 fprintf(fplog, "The initial domain decomposition cell size is:");
7293 for (d = 0; d < DIM; d++)
7297 fprintf(fplog, " %c %.2f nm",
7298 dim2char(d), dd->comm->cellsize_min[d]);
7301 fprintf(fplog, "\n\n");
7304 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7306 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7307 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7308 "non-bonded interactions", "", comm->cutoff);
7312 limit = dd->comm->cellsize_limit;
7316 if (dynamic_dd_box(ddbox, ir))
7318 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7320 limit = dd->comm->cellsize_min[XX];
7321 for (d = 1; d < DIM; d++)
7323 limit = std::min(limit, dd->comm->cellsize_min[d]);
7327 if (comm->bInterCGBondeds)
7329 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7330 "two-body bonded interactions", "(-rdd)",
7331 std::max(comm->cutoff, comm->cutoff_mbody));
7332 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7333 "multi-body bonded interactions", "(-rdd)",
7334 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7338 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7339 "virtual site constructions", "(-rcon)", limit);
7341 if (dd->constraint_comm)
7343 sprintf(buf, "atoms separated by up to %d constraints",
7345 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7346 buf, "(-rcon)", limit);
7348 fprintf(fplog, "\n");
7354 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7356 const t_inputrec *ir,
7357 const gmx_ddbox_t *ddbox)
7359 gmx_domdec_comm_t *comm;
7360 int d, dim, npulse, npulse_d_max, npulse_d;
7365 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7367 /* Determine the maximum number of comm. pulses in one dimension */
7369 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7371 /* Determine the maximum required number of grid pulses */
7372 if (comm->cellsize_limit >= comm->cutoff)
7374 /* Only a single pulse is required */
7377 else if (!bNoCutOff && comm->cellsize_limit > 0)
7379 /* We round down slightly here to avoid overhead due to the latency
7380 * of extra communication calls when the cut-off
7381 * would be only slightly longer than the cell size.
7382 * Later cellsize_limit is redetermined,
7383 * so we can not miss interactions due to this rounding.
7385 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7389 /* There is no cell size limit */
7390 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7393 if (!bNoCutOff && npulse > 1)
7395 /* See if we can do with less pulses, based on dlb_scale */
7397 for (d = 0; d < dd->ndim; d++)
7400 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7401 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7402 npulse_d_max = std::max(npulse_d_max, npulse_d);
7404 npulse = std::min(npulse, npulse_d_max);
7407 /* This env var can override npulse */
7408 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7415 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7416 for (d = 0; d < dd->ndim; d++)
7418 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7419 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7420 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7421 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7422 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7424 comm->bVacDLBNoLimit = FALSE;
7428 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7429 if (!comm->bVacDLBNoLimit)
7431 comm->cellsize_limit = std::max(comm->cellsize_limit,
7432 comm->cutoff/comm->maxpulse);
7434 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7435 /* Set the minimum cell size for each DD dimension */
7436 for (d = 0; d < dd->ndim; d++)
7438 if (comm->bVacDLBNoLimit ||
7439 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7441 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7445 comm->cellsize_min_dlb[dd->dim[d]] =
7446 comm->cutoff/comm->cd[d].np_dlb;
7449 if (comm->cutoff_mbody <= 0)
7451 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7453 if (comm->bDynLoadBal)
7459 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7461 /* If each molecule is a single charge group
7462 * or we use domain decomposition for each periodic dimension,
7463 * we do not need to take pbc into account for the bonded interactions.
7465 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7468 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7471 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7472 t_inputrec *ir, gmx_ddbox_t *ddbox)
7474 gmx_domdec_comm_t *comm;
7480 /* Initialize the thread data.
7481 * This can not be done in init_domain_decomposition,
7482 * as the numbers of threads is determined later.
7484 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7487 snew(comm->dth, comm->nth);
7490 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7492 init_ddpme(dd, &comm->ddpme[0], 0);
7493 if (comm->npmedecompdim >= 2)
7495 init_ddpme(dd, &comm->ddpme[1], 1);
7500 comm->npmenodes = 0;
7501 if (dd->pme_nodeid >= 0)
7503 gmx_fatal_collective(FARGS, NULL, dd,
7504 "Can not have separate PME ranks without PME electrostatics");
7510 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7512 if (comm->eDLB != edlbNO)
7514 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7517 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7518 if (comm->eDLB == edlbAUTO)
7522 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7524 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7527 if (ir->ePBC == epbcNONE)
7529 vol_frac = 1 - 1/(double)dd->nnodes;
7534 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7538 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7540 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7542 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7545 static gmx_bool test_dd_cutoff(t_commrec *cr,
7546 t_state *state, t_inputrec *ir,
7557 set_ddbox(dd, FALSE, cr, ir, state->box,
7558 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7562 for (d = 0; d < dd->ndim; d++)
7566 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7567 if (dynamic_dd_box(&ddbox, ir))
7569 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7572 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7574 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7575 dd->comm->cd[d].np_dlb > 0)
7577 if (np > dd->comm->cd[d].np_dlb)
7582 /* If a current local cell size is smaller than the requested
7583 * cut-off, we could still fix it, but this gets very complicated.
7584 * Without fixing here, we might actually need more checks.
7586 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7593 if (dd->comm->eDLB != edlbNO)
7595 /* If DLB is not active yet, we don't need to check the grid jumps.
7596 * Actually we shouldn't, because then the grid jump data is not set.
7598 if (dd->comm->bDynLoadBal &&
7599 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7604 gmx_sumi(1, &LocallyLimited, cr);
7606 if (LocallyLimited > 0)
7615 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7618 gmx_bool bCutoffAllowed;
7620 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7624 cr->dd->comm->cutoff = cutoff_req;
7627 return bCutoffAllowed;
7630 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
7632 gmx_domdec_comm_t *comm;
7634 comm = cr->dd->comm;
7636 /* Turn on the DLB limiting (might have been on already) */
7637 comm->bPMELoadBalDLBLimits = TRUE;
7639 /* Change the cut-off limit */
7640 comm->PMELoadBal_max_cutoff = cutoff;
7644 fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
7645 comm->PMELoadBal_max_cutoff);
7649 /* Sets whether we should later check the load imbalance data, so that
7650 * we can trigger dynamic load balancing if enough imbalance has
7653 * Used after PME load balancing unlocks DLB, so that the check
7654 * whether DLB will be useful can happen immediately.
7656 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
7658 if (dd->comm->eDLB == edlbAUTO && !dd_dlb_is_locked(dd))
7660 dd->comm->bCheckWhetherToTurnDlbOn = bValue;
7664 /* Returns if we should check whether there has been enough load
7665 * imbalance to trigger dynamic load balancing.
7667 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
7669 const int nddp_chk_dlb = 100;
7671 if (dd->comm->eDLB != edlbAUTO)
7676 /* We should check whether we should use DLB directly after
7678 if (dd->comm->bCheckWhetherToTurnDlbOn)
7680 /* This flag was set when the PME load-balancing routines
7681 unlocked DLB, and should now be cleared. */
7682 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
7685 /* We should also check whether we should use DLB every 100
7686 * partitionings (we do not do this every partioning, so that we
7687 * avoid excessive communication). */
7688 if (dd->comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1)
7696 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
7698 return dd->comm->bDynLoadBal;
7701 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7703 return dd->comm->bDLB_locked;
7706 void dd_dlb_lock(gmx_domdec_t *dd)
7708 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7709 if (dd->comm->eDLB == edlbAUTO)
7711 dd->comm->bDLB_locked = TRUE;
7715 void dd_dlb_unlock(gmx_domdec_t *dd)
7717 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7718 if (dd->comm->eDLB == edlbAUTO)
7720 dd->comm->bDLB_locked = FALSE;
7721 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, !dd->comm->bDynLoadBal);
7725 static void merge_cg_buffers(int ncell,
7726 gmx_domdec_comm_dim_t *cd, int pulse,
7728 int *index_gl, int *recv_i,
7729 rvec *cg_cm, rvec *recv_vr,
7731 cginfo_mb_t *cginfo_mb, int *cginfo)
7733 gmx_domdec_ind_t *ind, *ind_p;
7734 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7735 int shift, shift_at;
7737 ind = &cd->ind[pulse];
7739 /* First correct the already stored data */
7740 shift = ind->nrecv[ncell];
7741 for (cell = ncell-1; cell >= 0; cell--)
7743 shift -= ind->nrecv[cell];
7746 /* Move the cg's present from previous grid pulses */
7747 cg0 = ncg_cell[ncell+cell];
7748 cg1 = ncg_cell[ncell+cell+1];
7749 cgindex[cg1+shift] = cgindex[cg1];
7750 for (cg = cg1-1; cg >= cg0; cg--)
7752 index_gl[cg+shift] = index_gl[cg];
7753 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7754 cgindex[cg+shift] = cgindex[cg];
7755 cginfo[cg+shift] = cginfo[cg];
7757 /* Correct the already stored send indices for the shift */
7758 for (p = 1; p <= pulse; p++)
7760 ind_p = &cd->ind[p];
7762 for (c = 0; c < cell; c++)
7764 cg0 += ind_p->nsend[c];
7766 cg1 = cg0 + ind_p->nsend[cell];
7767 for (cg = cg0; cg < cg1; cg++)
7769 ind_p->index[cg] += shift;
7775 /* Merge in the communicated buffers */
7779 for (cell = 0; cell < ncell; cell++)
7781 cg1 = ncg_cell[ncell+cell+1] + shift;
7784 /* Correct the old cg indices */
7785 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7787 cgindex[cg+1] += shift_at;
7790 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7792 /* Copy this charge group from the buffer */
7793 index_gl[cg1] = recv_i[cg0];
7794 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7795 /* Add it to the cgindex */
7796 cg_gl = index_gl[cg1];
7797 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7798 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7799 cgindex[cg1+1] = cgindex[cg1] + nat;
7804 shift += ind->nrecv[cell];
7805 ncg_cell[ncell+cell+1] = cg1;
7809 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7810 int nzone, int cg0, const int *cgindex)
7814 /* Store the atom block boundaries for easy copying of communication buffers
7817 for (zone = 0; zone < nzone; zone++)
7819 for (p = 0; p < cd->np; p++)
7821 cd->ind[p].cell2at0[zone] = cgindex[cg];
7822 cg += cd->ind[p].nrecv[zone];
7823 cd->ind[p].cell2at1[zone] = cgindex[cg];
7828 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7834 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7836 if (!bLocalCG[link->a[i]])
7845 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7847 real c[DIM][4]; /* the corners for the non-bonded communication */
7848 real cr0; /* corner for rounding */
7849 real cr1[4]; /* corners for rounding */
7850 real bc[DIM]; /* corners for bounded communication */
7851 real bcr1; /* corner for rounding for bonded communication */
7854 /* Determine the corners of the domain(s) we are communicating with */
7856 set_dd_corners(const gmx_domdec_t *dd,
7857 int dim0, int dim1, int dim2,
7861 const gmx_domdec_comm_t *comm;
7862 const gmx_domdec_zones_t *zones;
7867 zones = &comm->zones;
7869 /* Keep the compiler happy */
7873 /* The first dimension is equal for all cells */
7874 c->c[0][0] = comm->cell_x0[dim0];
7877 c->bc[0] = c->c[0][0];
7882 /* This cell row is only seen from the first row */
7883 c->c[1][0] = comm->cell_x0[dim1];
7884 /* All rows can see this row */
7885 c->c[1][1] = comm->cell_x0[dim1];
7888 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7891 /* For the multi-body distance we need the maximum */
7892 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7895 /* Set the upper-right corner for rounding */
7896 c->cr0 = comm->cell_x1[dim0];
7901 for (j = 0; j < 4; j++)
7903 c->c[2][j] = comm->cell_x0[dim2];
7907 /* Use the maximum of the i-cells that see a j-cell */
7908 for (i = 0; i < zones->nizone; i++)
7910 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7915 std::max(c->c[2][j-4],
7916 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7922 /* For the multi-body distance we need the maximum */
7923 c->bc[2] = comm->cell_x0[dim2];
7924 for (i = 0; i < 2; i++)
7926 for (j = 0; j < 2; j++)
7928 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7934 /* Set the upper-right corner for rounding */
7935 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7936 * Only cell (0,0,0) can see cell 7 (1,1,1)
7938 c->cr1[0] = comm->cell_x1[dim1];
7939 c->cr1[3] = comm->cell_x1[dim1];
7942 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7945 /* For the multi-body distance we need the maximum */
7946 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7953 /* Determine which cg's we need to send in this pulse from this zone */
7955 get_zone_pulse_cgs(gmx_domdec_t *dd,
7956 int zonei, int zone,
7958 const int *index_gl,
7960 int dim, int dim_ind,
7961 int dim0, int dim1, int dim2,
7962 real r_comm2, real r_bcomm2,
7966 real skew_fac2_d, real skew_fac_01,
7967 rvec *v_d, rvec *v_0, rvec *v_1,
7968 const dd_corners_t *c,
7970 gmx_bool bDistBonded,
7976 gmx_domdec_ind_t *ind,
7977 int **ibuf, int *ibuf_nalloc,
7983 gmx_domdec_comm_t *comm;
7985 gmx_bool bDistMB_pulse;
7987 real r2, rb2, r, tric_sh;
7990 int nsend_z, nsend, nat;
7994 bScrew = (dd->bScrewPBC && dim == XX);
7996 bDistMB_pulse = (bDistMB && bDistBonded);
8002 for (cg = cg0; cg < cg1; cg++)
8006 if (tric_dist[dim_ind] == 0)
8008 /* Rectangular direction, easy */
8009 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
8016 r = cg_cm[cg][dim] - c->bc[dim_ind];
8022 /* Rounding gives at most a 16% reduction
8023 * in communicated atoms
8025 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8027 r = cg_cm[cg][dim0] - c->cr0;
8028 /* This is the first dimension, so always r >= 0 */
8035 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8037 r = cg_cm[cg][dim1] - c->cr1[zone];
8044 r = cg_cm[cg][dim1] - c->bcr1;
8054 /* Triclinic direction, more complicated */
8057 /* Rounding, conservative as the skew_fac multiplication
8058 * will slightly underestimate the distance.
8060 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8062 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
8063 for (i = dim0+1; i < DIM; i++)
8065 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
8067 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
8070 rb[dim0] = rn[dim0];
8073 /* Take care that the cell planes along dim0 might not
8074 * be orthogonal to those along dim1 and dim2.
8076 for (i = 1; i <= dim_ind; i++)
8079 if (normal[dim0][dimd] > 0)
8081 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
8084 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
8089 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8091 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
8093 for (i = dim1+1; i < DIM; i++)
8095 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
8097 rn[dim1] += tric_sh;
8100 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
8101 /* Take care of coupling of the distances
8102 * to the planes along dim0 and dim1 through dim2.
8104 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
8105 /* Take care that the cell planes along dim1
8106 * might not be orthogonal to that along dim2.
8108 if (normal[dim1][dim2] > 0)
8110 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8116 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8119 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8120 /* Take care of coupling of the distances
8121 * to the planes along dim0 and dim1 through dim2.
8123 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8124 /* Take care that the cell planes along dim1
8125 * might not be orthogonal to that along dim2.
8127 if (normal[dim1][dim2] > 0)
8129 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8134 /* The distance along the communication direction */
8135 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8137 for (i = dim+1; i < DIM; i++)
8139 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8144 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8145 /* Take care of coupling of the distances
8146 * to the planes along dim0 and dim1 through dim2.
8148 if (dim_ind == 1 && zonei == 1)
8150 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8156 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8159 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8160 /* Take care of coupling of the distances
8161 * to the planes along dim0 and dim1 through dim2.
8163 if (dim_ind == 1 && zonei == 1)
8165 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8173 ((bDistMB && rb2 < r_bcomm2) ||
8174 (bDist2B && r2 < r_bcomm2)) &&
8176 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8177 missing_link(comm->cglink, index_gl[cg],
8180 /* Make an index to the local charge groups */
8181 if (nsend+1 > ind->nalloc)
8183 ind->nalloc = over_alloc_large(nsend+1);
8184 srenew(ind->index, ind->nalloc);
8186 if (nsend+1 > *ibuf_nalloc)
8188 *ibuf_nalloc = over_alloc_large(nsend+1);
8189 srenew(*ibuf, *ibuf_nalloc);
8191 ind->index[nsend] = cg;
8192 (*ibuf)[nsend] = index_gl[cg];
8194 vec_rvec_check_alloc(vbuf, nsend+1);
8196 if (dd->ci[dim] == 0)
8198 /* Correct cg_cm for pbc */
8199 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8202 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8203 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8208 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8211 nat += cgindex[cg+1] - cgindex[cg];
8217 *nsend_z_ptr = nsend_z;
8220 static void setup_dd_communication(gmx_domdec_t *dd,
8221 matrix box, gmx_ddbox_t *ddbox,
8222 t_forcerec *fr, t_state *state, rvec **f)
8224 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8225 int nzone, nzone_send, zone, zonei, cg0, cg1;
8226 int c, i, cg, cg_gl, nrcg;
8227 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8228 gmx_domdec_comm_t *comm;
8229 gmx_domdec_zones_t *zones;
8230 gmx_domdec_comm_dim_t *cd;
8231 gmx_domdec_ind_t *ind;
8232 cginfo_mb_t *cginfo_mb;
8233 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8234 real r_comm2, r_bcomm2;
8235 dd_corners_t corners;
8237 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8238 real skew_fac2_d, skew_fac_01;
8245 fprintf(debug, "Setting up DD communication\n");
8250 switch (fr->cutoff_scheme)
8259 gmx_incons("unimplemented");
8263 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8265 /* Check if we need to use triclinic distances */
8266 tric_dist[dim_ind] = 0;
8267 for (i = 0; i <= dim_ind; i++)
8269 if (ddbox->tric_dir[dd->dim[i]])
8271 tric_dist[dim_ind] = 1;
8276 bBondComm = comm->bBondComm;
8278 /* Do we need to determine extra distances for multi-body bondeds? */
8279 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8281 /* Do we need to determine extra distances for only two-body bondeds? */
8282 bDist2B = (bBondComm && !bDistMB);
8284 r_comm2 = sqr(comm->cutoff);
8285 r_bcomm2 = sqr(comm->cutoff_mbody);
8289 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8292 zones = &comm->zones;
8295 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8296 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8298 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8300 /* Triclinic stuff */
8301 normal = ddbox->normal;
8305 v_0 = ddbox->v[dim0];
8306 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8308 /* Determine the coupling coefficient for the distances
8309 * to the cell planes along dim0 and dim1 through dim2.
8310 * This is required for correct rounding.
8313 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8316 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8322 v_1 = ddbox->v[dim1];
8325 zone_cg_range = zones->cg_range;
8326 index_gl = dd->index_gl;
8327 cgindex = dd->cgindex;
8328 cginfo_mb = fr->cginfo_mb;
8330 zone_cg_range[0] = 0;
8331 zone_cg_range[1] = dd->ncg_home;
8332 comm->zone_ncg1[0] = dd->ncg_home;
8333 pos_cg = dd->ncg_home;
8335 nat_tot = dd->nat_home;
8337 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8339 dim = dd->dim[dim_ind];
8340 cd = &comm->cd[dim_ind];
8342 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8344 /* No pbc in this dimension, the first node should not comm. */
8352 v_d = ddbox->v[dim];
8353 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8355 cd->bInPlace = TRUE;
8356 for (p = 0; p < cd->np; p++)
8358 /* Only atoms communicated in the first pulse are used
8359 * for multi-body bonded interactions or for bBondComm.
8361 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8366 for (zone = 0; zone < nzone_send; zone++)
8368 if (tric_dist[dim_ind] && dim_ind > 0)
8370 /* Determine slightly more optimized skew_fac's
8372 * This reduces the number of communicated atoms
8373 * by about 10% for 3D DD of rhombic dodecahedra.
8375 for (dimd = 0; dimd < dim; dimd++)
8377 sf2_round[dimd] = 1;
8378 if (ddbox->tric_dir[dimd])
8380 for (i = dd->dim[dimd]+1; i < DIM; i++)
8382 /* If we are shifted in dimension i
8383 * and the cell plane is tilted forward
8384 * in dimension i, skip this coupling.
8386 if (!(zones->shift[nzone+zone][i] &&
8387 ddbox->v[dimd][i][dimd] >= 0))
8390 sqr(ddbox->v[dimd][i][dimd]);
8393 sf2_round[dimd] = 1/sf2_round[dimd];
8398 zonei = zone_perm[dim_ind][zone];
8401 /* Here we permutate the zones to obtain a convenient order
8402 * for neighbor searching
8404 cg0 = zone_cg_range[zonei];
8405 cg1 = zone_cg_range[zonei+1];
8409 /* Look only at the cg's received in the previous grid pulse
8411 cg1 = zone_cg_range[nzone+zone+1];
8412 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8415 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8416 for (th = 0; th < comm->nth; th++)
8418 gmx_domdec_ind_t *ind_p;
8419 int **ibuf_p, *ibuf_nalloc_p;
8421 int *nsend_p, *nat_p;
8427 /* Thread 0 writes in the comm buffers */
8429 ibuf_p = &comm->buf_int;
8430 ibuf_nalloc_p = &comm->nalloc_int;
8431 vbuf_p = &comm->vbuf;
8434 nsend_zone_p = &ind->nsend[zone];
8438 /* Other threads write into temp buffers */
8439 ind_p = &comm->dth[th].ind;
8440 ibuf_p = &comm->dth[th].ibuf;
8441 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8442 vbuf_p = &comm->dth[th].vbuf;
8443 nsend_p = &comm->dth[th].nsend;
8444 nat_p = &comm->dth[th].nat;
8445 nsend_zone_p = &comm->dth[th].nsend_zone;
8447 comm->dth[th].nsend = 0;
8448 comm->dth[th].nat = 0;
8449 comm->dth[th].nsend_zone = 0;
8459 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8460 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8463 /* Get the cg's for this pulse in this zone */
8464 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8466 dim, dim_ind, dim0, dim1, dim2,
8469 normal, skew_fac2_d, skew_fac_01,
8470 v_d, v_0, v_1, &corners, sf2_round,
8471 bDistBonded, bBondComm,
8475 ibuf_p, ibuf_nalloc_p,
8481 /* Append data of threads>=1 to the communication buffers */
8482 for (th = 1; th < comm->nth; th++)
8484 dd_comm_setup_work_t *dth;
8487 dth = &comm->dth[th];
8489 ns1 = nsend + dth->nsend_zone;
8490 if (ns1 > ind->nalloc)
8492 ind->nalloc = over_alloc_dd(ns1);
8493 srenew(ind->index, ind->nalloc);
8495 if (ns1 > comm->nalloc_int)
8497 comm->nalloc_int = over_alloc_dd(ns1);
8498 srenew(comm->buf_int, comm->nalloc_int);
8500 if (ns1 > comm->vbuf.nalloc)
8502 comm->vbuf.nalloc = over_alloc_dd(ns1);
8503 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8506 for (i = 0; i < dth->nsend_zone; i++)
8508 ind->index[nsend] = dth->ind.index[i];
8509 comm->buf_int[nsend] = dth->ibuf[i];
8510 copy_rvec(dth->vbuf.v[i],
8511 comm->vbuf.v[nsend]);
8515 ind->nsend[zone] += dth->nsend_zone;
8518 /* Clear the counts in case we do not have pbc */
8519 for (zone = nzone_send; zone < nzone; zone++)
8521 ind->nsend[zone] = 0;
8523 ind->nsend[nzone] = nsend;
8524 ind->nsend[nzone+1] = nat;
8525 /* Communicate the number of cg's and atoms to receive */
8526 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8527 ind->nsend, nzone+2,
8528 ind->nrecv, nzone+2);
8530 /* The rvec buffer is also required for atom buffers of size nsend
8531 * in dd_move_x and dd_move_f.
8533 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8537 /* We can receive in place if only the last zone is not empty */
8538 for (zone = 0; zone < nzone-1; zone++)
8540 if (ind->nrecv[zone] > 0)
8542 cd->bInPlace = FALSE;
8547 /* The int buffer is only required here for the cg indices */
8548 if (ind->nrecv[nzone] > comm->nalloc_int2)
8550 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8551 srenew(comm->buf_int2, comm->nalloc_int2);
8553 /* The rvec buffer is also required for atom buffers
8554 * of size nrecv in dd_move_x and dd_move_f.
8556 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8557 vec_rvec_check_alloc(&comm->vbuf2, i);
8561 /* Make space for the global cg indices */
8562 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8563 || dd->cg_nalloc == 0)
8565 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8566 srenew(index_gl, dd->cg_nalloc);
8567 srenew(cgindex, dd->cg_nalloc+1);
8569 /* Communicate the global cg indices */
8572 recv_i = index_gl + pos_cg;
8576 recv_i = comm->buf_int2;
8578 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8579 comm->buf_int, nsend,
8580 recv_i, ind->nrecv[nzone]);
8582 /* Make space for cg_cm */
8583 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8584 if (fr->cutoff_scheme == ecutsGROUP)
8592 /* Communicate cg_cm */
8595 recv_vr = cg_cm + pos_cg;
8599 recv_vr = comm->vbuf2.v;
8601 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8602 comm->vbuf.v, nsend,
8603 recv_vr, ind->nrecv[nzone]);
8605 /* Make the charge group index */
8608 zone = (p == 0 ? 0 : nzone - 1);
8609 while (zone < nzone)
8611 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8613 cg_gl = index_gl[pos_cg];
8614 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8615 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8616 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8619 /* Update the charge group presence,
8620 * so we can use it in the next pass of the loop.
8622 comm->bLocalCG[cg_gl] = TRUE;
8628 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8631 zone_cg_range[nzone+zone] = pos_cg;
8636 /* This part of the code is never executed with bBondComm. */
8637 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8638 index_gl, recv_i, cg_cm, recv_vr,
8639 cgindex, fr->cginfo_mb, fr->cginfo);
8640 pos_cg += ind->nrecv[nzone];
8642 nat_tot += ind->nrecv[nzone+1];
8646 /* Store the atom block for easy copying of communication buffers */
8647 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8651 dd->index_gl = index_gl;
8652 dd->cgindex = cgindex;
8654 dd->ncg_tot = zone_cg_range[zones->n];
8655 dd->nat_tot = nat_tot;
8656 comm->nat[ddnatHOME] = dd->nat_home;
8657 for (i = ddnatZONE; i < ddnatNR; i++)
8659 comm->nat[i] = dd->nat_tot;
8664 /* We don't need to update cginfo, since that was alrady done above.
8665 * So we pass NULL for the forcerec.
8667 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8668 NULL, comm->bLocalCG);
8673 fprintf(debug, "Finished setting up DD communication, zones:");
8674 for (c = 0; c < zones->n; c++)
8676 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8678 fprintf(debug, "\n");
8682 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8686 for (c = 0; c < zones->nizone; c++)
8688 zones->izone[c].cg1 = zones->cg_range[c+1];
8689 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8690 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8694 static void set_zones_size(gmx_domdec_t *dd,
8695 matrix box, const gmx_ddbox_t *ddbox,
8696 int zone_start, int zone_end)
8698 gmx_domdec_comm_t *comm;
8699 gmx_domdec_zones_t *zones;
8708 zones = &comm->zones;
8710 /* Do we need to determine extra distances for multi-body bondeds? */
8711 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8713 for (z = zone_start; z < zone_end; z++)
8715 /* Copy cell limits to zone limits.
8716 * Valid for non-DD dims and non-shifted dims.
8718 copy_rvec(comm->cell_x0, zones->size[z].x0);
8719 copy_rvec(comm->cell_x1, zones->size[z].x1);
8722 for (d = 0; d < dd->ndim; d++)
8726 for (z = 0; z < zones->n; z++)
8728 /* With a staggered grid we have different sizes
8729 * for non-shifted dimensions.
8731 if (dd->bGridJump && zones->shift[z][dim] == 0)
8735 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8736 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8740 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8741 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8747 rcmbs = comm->cutoff_mbody;
8748 if (ddbox->tric_dir[dim])
8750 rcs /= ddbox->skew_fac[dim];
8751 rcmbs /= ddbox->skew_fac[dim];
8754 /* Set the lower limit for the shifted zone dimensions */
8755 for (z = zone_start; z < zone_end; z++)
8757 if (zones->shift[z][dim] > 0)
8760 if (!dd->bGridJump || d == 0)
8762 zones->size[z].x0[dim] = comm->cell_x1[dim];
8763 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8767 /* Here we take the lower limit of the zone from
8768 * the lowest domain of the zone below.
8772 zones->size[z].x0[dim] =
8773 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8779 zones->size[z].x0[dim] =
8780 zones->size[zone_perm[2][z-4]].x0[dim];
8784 zones->size[z].x0[dim] =
8785 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8788 /* A temporary limit, is updated below */
8789 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8793 for (zi = 0; zi < zones->nizone; zi++)
8795 if (zones->shift[zi][dim] == 0)
8797 /* This takes the whole zone into account.
8798 * With multiple pulses this will lead
8799 * to a larger zone then strictly necessary.
8801 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8802 zones->size[zi].x1[dim]+rcmbs);
8810 /* Loop over the i-zones to set the upper limit of each
8813 for (zi = 0; zi < zones->nizone; zi++)
8815 if (zones->shift[zi][dim] == 0)
8817 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8819 if (zones->shift[z][dim] > 0)
8821 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8822 zones->size[zi].x1[dim]+rcs);
8829 for (z = zone_start; z < zone_end; z++)
8831 /* Initialization only required to keep the compiler happy */
8832 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8835 /* To determine the bounding box for a zone we need to find
8836 * the extreme corners of 4, 2 or 1 corners.
8838 nc = 1 << (ddbox->nboundeddim - 1);
8840 for (c = 0; c < nc; c++)
8842 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8846 corner[YY] = zones->size[z].x0[YY];
8850 corner[YY] = zones->size[z].x1[YY];
8854 corner[ZZ] = zones->size[z].x0[ZZ];
8858 corner[ZZ] = zones->size[z].x1[ZZ];
8860 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8861 box[ZZ][1 - dd->dim[0]] != 0)
8863 /* With 1D domain decomposition the cg's are not in
8864 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8865 * Shift the corner of the z-vector back to along the box
8866 * vector of dimension d, so it will later end up at 0 along d.
8867 * This can affect the location of this corner along dd->dim[0]
8868 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8870 int d = 1 - dd->dim[0];
8872 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8874 /* Apply the triclinic couplings */
8875 assert(ddbox->npbcdim <= DIM);
8876 for (i = YY; i < ddbox->npbcdim; i++)
8878 for (j = XX; j < i; j++)
8880 corner[j] += corner[i]*box[i][j]/box[i][i];
8885 copy_rvec(corner, corner_min);
8886 copy_rvec(corner, corner_max);
8890 for (i = 0; i < DIM; i++)
8892 corner_min[i] = std::min(corner_min[i], corner[i]);
8893 corner_max[i] = std::max(corner_max[i], corner[i]);
8897 /* Copy the extreme cornes without offset along x */
8898 for (i = 0; i < DIM; i++)
8900 zones->size[z].bb_x0[i] = corner_min[i];
8901 zones->size[z].bb_x1[i] = corner_max[i];
8903 /* Add the offset along x */
8904 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8905 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8908 if (zone_start == 0)
8911 for (dim = 0; dim < DIM; dim++)
8913 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8915 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8920 for (z = zone_start; z < zone_end; z++)
8922 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8924 zones->size[z].x0[XX], zones->size[z].x1[XX],
8925 zones->size[z].x0[YY], zones->size[z].x1[YY],
8926 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8927 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8929 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8930 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8931 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8936 static int comp_cgsort(const void *a, const void *b)
8940 gmx_cgsort_t *cga, *cgb;
8941 cga = (gmx_cgsort_t *)a;
8942 cgb = (gmx_cgsort_t *)b;
8944 comp = cga->nsc - cgb->nsc;
8947 comp = cga->ind_gl - cgb->ind_gl;
8953 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8958 /* Order the data */
8959 for (i = 0; i < n; i++)
8961 buf[i] = a[sort[i].ind];
8964 /* Copy back to the original array */
8965 for (i = 0; i < n; i++)
8971 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8976 /* Order the data */
8977 for (i = 0; i < n; i++)
8979 copy_rvec(v[sort[i].ind], buf[i]);
8982 /* Copy back to the original array */
8983 for (i = 0; i < n; i++)
8985 copy_rvec(buf[i], v[i]);
8989 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8992 int a, atot, cg, cg0, cg1, i;
8994 if (cgindex == NULL)
8996 /* Avoid the useless loop of the atoms within a cg */
8997 order_vec_cg(ncg, sort, v, buf);
9002 /* Order the data */
9004 for (cg = 0; cg < ncg; cg++)
9006 cg0 = cgindex[sort[cg].ind];
9007 cg1 = cgindex[sort[cg].ind+1];
9008 for (i = cg0; i < cg1; i++)
9010 copy_rvec(v[i], buf[a]);
9016 /* Copy back to the original array */
9017 for (a = 0; a < atot; a++)
9019 copy_rvec(buf[a], v[a]);
9023 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
9024 int nsort_new, gmx_cgsort_t *sort_new,
9025 gmx_cgsort_t *sort1)
9029 /* The new indices are not very ordered, so we qsort them */
9030 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
9032 /* sort2 is already ordered, so now we can merge the two arrays */
9036 while (i2 < nsort2 || i_new < nsort_new)
9040 sort1[i1++] = sort_new[i_new++];
9042 else if (i_new == nsort_new)
9044 sort1[i1++] = sort2[i2++];
9046 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
9047 (sort2[i2].nsc == sort_new[i_new].nsc &&
9048 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
9050 sort1[i1++] = sort2[i2++];
9054 sort1[i1++] = sort_new[i_new++];
9059 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
9061 gmx_domdec_sort_t *sort;
9062 gmx_cgsort_t *cgsort, *sort_i;
9063 int ncg_new, nsort2, nsort_new, i, *a, moved;
9065 sort = dd->comm->sort;
9067 a = fr->ns.grid->cell_index;
9069 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
9071 if (ncg_home_old >= 0)
9073 /* The charge groups that remained in the same ns grid cell
9074 * are completely ordered. So we can sort efficiently by sorting
9075 * the charge groups that did move into the stationary list.
9080 for (i = 0; i < dd->ncg_home; i++)
9082 /* Check if this cg did not move to another node */
9085 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
9087 /* This cg is new on this node or moved ns grid cell */
9088 if (nsort_new >= sort->sort_new_nalloc)
9090 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
9091 srenew(sort->sort_new, sort->sort_new_nalloc);
9093 sort_i = &(sort->sort_new[nsort_new++]);
9097 /* This cg did not move */
9098 sort_i = &(sort->sort2[nsort2++]);
9100 /* Sort on the ns grid cell indices
9101 * and the global topology index.
9102 * index_gl is irrelevant with cell ns,
9103 * but we set it here anyhow to avoid a conditional.
9106 sort_i->ind_gl = dd->index_gl[i];
9113 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9116 /* Sort efficiently */
9117 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9122 cgsort = sort->sort;
9124 for (i = 0; i < dd->ncg_home; i++)
9126 /* Sort on the ns grid cell indices
9127 * and the global topology index
9129 cgsort[i].nsc = a[i];
9130 cgsort[i].ind_gl = dd->index_gl[i];
9132 if (cgsort[i].nsc < moved)
9139 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9141 /* Determine the order of the charge groups using qsort */
9142 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9148 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9151 int ncg_new, i, *a, na;
9153 sort = dd->comm->sort->sort;
9155 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9158 for (i = 0; i < na; i++)
9162 sort[ncg_new].ind = a[i];
9170 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9173 gmx_domdec_sort_t *sort;
9174 gmx_cgsort_t *cgsort;
9176 int ncg_new, i, *ibuf, cgsize;
9179 sort = dd->comm->sort;
9181 if (dd->ncg_home > sort->sort_nalloc)
9183 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9184 srenew(sort->sort, sort->sort_nalloc);
9185 srenew(sort->sort2, sort->sort_nalloc);
9187 cgsort = sort->sort;
9189 switch (fr->cutoff_scheme)
9192 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9195 ncg_new = dd_sort_order_nbnxn(dd, fr);
9198 gmx_incons("unimplemented");
9202 /* We alloc with the old size, since cgindex is still old */
9203 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9204 vbuf = dd->comm->vbuf.v;
9208 cgindex = dd->cgindex;
9215 /* Remove the charge groups which are no longer at home here */
9216 dd->ncg_home = ncg_new;
9219 fprintf(debug, "Set the new home charge group count to %d\n",
9223 /* Reorder the state */
9224 for (i = 0; i < estNR; i++)
9226 if (EST_DISTR(i) && (state->flags & (1<<i)))
9231 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9234 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9237 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9240 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9244 case estDISRE_INITF:
9245 case estDISRE_RM3TAV:
9246 case estORIRE_INITF:
9248 /* No ordering required */
9251 gmx_incons("Unknown state entry encountered in dd_sort_state");
9256 if (fr->cutoff_scheme == ecutsGROUP)
9259 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9262 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9264 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9265 srenew(sort->ibuf, sort->ibuf_nalloc);
9268 /* Reorder the global cg index */
9269 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9270 /* Reorder the cginfo */
9271 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9272 /* Rebuild the local cg index */
9276 for (i = 0; i < dd->ncg_home; i++)
9278 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9279 ibuf[i+1] = ibuf[i] + cgsize;
9281 for (i = 0; i < dd->ncg_home+1; i++)
9283 dd->cgindex[i] = ibuf[i];
9288 for (i = 0; i < dd->ncg_home+1; i++)
9293 /* Set the home atom number */
9294 dd->nat_home = dd->cgindex[dd->ncg_home];
9296 if (fr->cutoff_scheme == ecutsVERLET)
9298 /* The atoms are now exactly in grid order, update the grid order */
9299 nbnxn_set_atomorder(fr->nbv->nbs);
9303 /* Copy the sorted ns cell indices back to the ns grid struct */
9304 for (i = 0; i < dd->ncg_home; i++)
9306 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9308 fr->ns.grid->nr = dd->ncg_home;
9312 static void add_dd_statistics(gmx_domdec_t *dd)
9314 gmx_domdec_comm_t *comm;
9319 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9321 comm->sum_nat[ddnat-ddnatZONE] +=
9322 comm->nat[ddnat] - comm->nat[ddnat-1];
9327 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9329 gmx_domdec_comm_t *comm;
9334 /* Reset all the statistics and counters for total run counting */
9335 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9337 comm->sum_nat[ddnat-ddnatZONE] = 0;
9341 comm->load_step = 0;
9344 clear_ivec(comm->load_lim);
9349 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9351 gmx_domdec_comm_t *comm;
9355 comm = cr->dd->comm;
9357 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9364 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");
9366 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9368 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9373 " av. #atoms communicated per step for force: %d x %.1f\n",
9377 if (cr->dd->vsite_comm)
9380 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9381 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9386 if (cr->dd->constraint_comm)
9389 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9390 1 + ir->nLincsIter, av);
9394 gmx_incons(" Unknown type for DD statistics");
9397 fprintf(fplog, "\n");
9399 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9401 print_dd_load_av(fplog, cr->dd);
9405 void dd_partition_system(FILE *fplog,
9408 gmx_bool bMasterState,
9410 t_state *state_global,
9411 gmx_mtop_t *top_global,
9413 t_state *state_local,
9416 gmx_localtop_t *top_local,
9419 gmx_shellfc_t shellfc,
9420 gmx_constr_t constr,
9422 gmx_wallcycle_t wcycle,
9426 gmx_domdec_comm_t *comm;
9427 gmx_ddbox_t ddbox = {0};
9429 gmx_int64_t step_pcoupl;
9430 rvec cell_ns_x0, cell_ns_x1;
9431 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9432 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bTurnOnDLB, bLogLoad;
9433 gmx_bool bRedist, bSortCG, bResortAll;
9434 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9438 wallcycle_start(wcycle, ewcDOMDEC);
9443 bBoxChanged = (bMasterState || DEFORM(*ir));
9444 if (ir->epc != epcNO)
9446 /* With nstpcouple > 1 pressure coupling happens.
9447 * one step after calculating the pressure.
9448 * Box scaling happens at the end of the MD step,
9449 * after the DD partitioning.
9450 * We therefore have to do DLB in the first partitioning
9451 * after an MD step where P-coupling occured.
9452 * We need to determine the last step in which p-coupling occurred.
9453 * MRS -- need to validate this for vv?
9458 step_pcoupl = step - 1;
9462 step_pcoupl = ((step - 1)/n)*n + 1;
9464 if (step_pcoupl >= comm->partition_step)
9470 bNStGlobalComm = (step % nstglobalcomm == 0);
9472 if (!comm->bDynLoadBal)
9478 /* Should we do dynamic load balacing this step?
9479 * Since it requires (possibly expensive) global communication,
9480 * we might want to do DLB less frequently.
9482 if (bBoxChanged || ir->epc != epcNO)
9484 bDoDLB = bBoxChanged;
9488 bDoDLB = bNStGlobalComm;
9492 /* Check if we have recorded loads on the nodes */
9493 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9495 bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
9497 /* Print load every nstlog, first and last step to the log file */
9498 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9499 comm->n_load_collect == 0 ||
9501 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9503 /* Avoid extra communication due to verbose screen output
9504 * when nstglobalcomm is set.
9506 if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
9507 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9509 get_load_distribution(dd, wcycle);
9514 dd_print_load(fplog, dd, step-1);
9518 dd_print_load_verbose(dd);
9521 comm->n_load_collect++;
9523 if (bCheckWhetherToTurnDlbOn)
9525 /* Since the timings are node dependent, the master decides */
9528 /* Here we check if the max PME rank load is more than 0.98
9529 * the max PP force load. If so, PP DLB will not help,
9530 * since we are (almost) limited by PME. Furthermore,
9531 * DLB will cause a significant extra x/f redistribution
9532 * cost on the PME ranks, which will then surely result
9533 * in lower total performance.
9534 * This check might be fragile, since one measurement
9535 * below 0.98 (although only done once every 100 DD part.)
9536 * could turn on DLB for the rest of the run.
9538 if (cr->npmenodes > 0 &&
9539 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9546 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9550 fprintf(debug, "step %s, imb loss %f\n",
9551 gmx_step_str(step, sbuf),
9552 dd_force_imb_perf_loss(dd));
9555 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9558 turn_on_dlb(fplog, cr, step);
9563 comm->n_load_have++;
9566 cgs_gl = &comm->cgs_gl;
9571 /* Clear the old state */
9572 clear_dd_indices(dd, 0, 0);
9575 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9576 TRUE, cgs_gl, state_global->x, &ddbox);
9578 get_cg_distribution(fplog, dd, cgs_gl,
9579 state_global->box, &ddbox, state_global->x);
9581 dd_distribute_state(dd, cgs_gl,
9582 state_global, state_local, f);
9584 dd_make_local_cgs(dd, &top_local->cgs);
9586 /* Ensure that we have space for the new distribution */
9587 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9589 if (fr->cutoff_scheme == ecutsGROUP)
9591 calc_cgcm(fplog, 0, dd->ncg_home,
9592 &top_local->cgs, state_local->x, fr->cg_cm);
9595 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9597 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9599 else if (state_local->ddp_count != dd->ddp_count)
9601 if (state_local->ddp_count > dd->ddp_count)
9603 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9606 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9608 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);
9611 /* Clear the old state */
9612 clear_dd_indices(dd, 0, 0);
9614 /* Build the new indices */
9615 rebuild_cgindex(dd, cgs_gl->index, state_local);
9616 make_dd_indices(dd, cgs_gl->index, 0);
9617 ncgindex_set = dd->ncg_home;
9619 if (fr->cutoff_scheme == ecutsGROUP)
9621 /* Redetermine the cg COMs */
9622 calc_cgcm(fplog, 0, dd->ncg_home,
9623 &top_local->cgs, state_local->x, fr->cg_cm);
9626 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9628 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9630 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9631 TRUE, &top_local->cgs, state_local->x, &ddbox);
9633 bRedist = comm->bDynLoadBal;
9637 /* We have the full state, only redistribute the cgs */
9639 /* Clear the non-home indices */
9640 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9643 /* Avoid global communication for dim's without pbc and -gcom */
9644 if (!bNStGlobalComm)
9646 copy_rvec(comm->box0, ddbox.box0 );
9647 copy_rvec(comm->box_size, ddbox.box_size);
9649 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9650 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9655 /* For dim's without pbc and -gcom */
9656 copy_rvec(ddbox.box0, comm->box0 );
9657 copy_rvec(ddbox.box_size, comm->box_size);
9659 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9662 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9664 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9667 /* Check if we should sort the charge groups */
9668 if (comm->nstSortCG > 0)
9670 bSortCG = (bMasterState ||
9671 (bRedist && (step % comm->nstSortCG == 0)));
9678 ncg_home_old = dd->ncg_home;
9683 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9685 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9687 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9689 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9692 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9694 &comm->cell_x0, &comm->cell_x1,
9695 dd->ncg_home, fr->cg_cm,
9696 cell_ns_x0, cell_ns_x1, &grid_density);
9700 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9703 switch (fr->cutoff_scheme)
9706 copy_ivec(fr->ns.grid->n, ncells_old);
9707 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9708 state_local->box, cell_ns_x0, cell_ns_x1,
9709 fr->rlistlong, grid_density);
9712 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9715 gmx_incons("unimplemented");
9717 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9718 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9722 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9724 /* Sort the state on charge group position.
9725 * This enables exact restarts from this step.
9726 * It also improves performance by about 15% with larger numbers
9727 * of atoms per node.
9730 /* Fill the ns grid with the home cell,
9731 * so we can sort with the indices.
9733 set_zones_ncg_home(dd);
9735 switch (fr->cutoff_scheme)
9738 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9740 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9742 comm->zones.size[0].bb_x0,
9743 comm->zones.size[0].bb_x1,
9745 comm->zones.dens_zone0,
9748 ncg_moved, bRedist ? comm->moved : NULL,
9749 fr->nbv->grp[eintLocal].kernel_type,
9750 fr->nbv->grp[eintLocal].nbat);
9752 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9755 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9756 0, dd->ncg_home, fr->cg_cm);
9758 copy_ivec(fr->ns.grid->n, ncells_new);
9761 gmx_incons("unimplemented");
9764 bResortAll = bMasterState;
9766 /* Check if we can user the old order and ns grid cell indices
9767 * of the charge groups to sort the charge groups efficiently.
9769 if (ncells_new[XX] != ncells_old[XX] ||
9770 ncells_new[YY] != ncells_old[YY] ||
9771 ncells_new[ZZ] != ncells_old[ZZ])
9778 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9779 gmx_step_str(step, sbuf), dd->ncg_home);
9781 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9782 bResortAll ? -1 : ncg_home_old);
9783 /* Rebuild all the indices */
9784 ga2la_clear(dd->ga2la);
9787 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9790 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9792 /* Setup up the communication and communicate the coordinates */
9793 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9795 /* Set the indices */
9796 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9798 /* Set the charge group boundaries for neighbor searching */
9799 set_cg_boundaries(&comm->zones);
9801 if (fr->cutoff_scheme == ecutsVERLET)
9803 set_zones_size(dd, state_local->box, &ddbox,
9804 bSortCG ? 1 : 0, comm->zones.n);
9807 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9810 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9811 -1,state_local->x,state_local->box);
9814 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9816 /* Extract a local topology from the global topology */
9817 for (i = 0; i < dd->ndim; i++)
9819 np[dd->dim[i]] = comm->cd[i].np;
9821 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9822 comm->cellsize_min, np,
9824 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9825 vsite, top_global, top_local);
9827 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9829 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9831 /* Set up the special atom communication */
9832 n = comm->nat[ddnatZONE];
9833 for (i = ddnatZONE+1; i < ddnatNR; i++)
9838 if (vsite && vsite->n_intercg_vsite)
9840 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9844 if (dd->bInterCGcons || dd->bInterCGsettles)
9846 /* Only for inter-cg constraints we need special code */
9847 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9848 constr, ir->nProjOrder,
9849 top_local->idef.il);
9853 gmx_incons("Unknown special atom type setup");
9858 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9860 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9862 /* Make space for the extra coordinates for virtual site
9863 * or constraint communication.
9865 state_local->natoms = comm->nat[ddnatNR-1];
9866 if (state_local->natoms > state_local->nalloc)
9868 dd_realloc_state(state_local, f, state_local->natoms);
9871 if (fr->bF_NoVirSum)
9873 if (vsite && vsite->n_intercg_vsite)
9875 nat_f_novirsum = comm->nat[ddnatVSITE];
9879 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9881 nat_f_novirsum = dd->nat_tot;
9885 nat_f_novirsum = dd->nat_home;
9894 /* Set the number of atoms required for the force calculation.
9895 * Forces need to be constrained when using a twin-range setup
9896 * or with energy minimization. For simple simulations we could
9897 * avoid some allocation, zeroing and copying, but this is
9898 * probably not worth the complications ande checking.
9900 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9901 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9903 /* We make the all mdatoms up to nat_tot_con.
9904 * We could save some work by only setting invmass
9905 * between nat_tot and nat_tot_con.
9907 /* This call also sets the new number of home particles to dd->nat_home */
9908 atoms2md(top_global, ir,
9909 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9911 /* Now we have the charges we can sort the FE interactions */
9912 dd_sort_local_top(dd, mdatoms, top_local);
9916 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9917 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9918 mdatoms, FALSE, vsite);
9923 /* Make the local shell stuff, currently no communication is done */
9924 make_local_shells(cr, mdatoms, shellfc);
9927 if (ir->implicit_solvent)
9929 make_local_gb(cr, fr->born, ir->gb_algorithm);
9932 setup_bonded_threading(fr, &top_local->idef);
9934 if (!(cr->duty & DUTY_PME))
9936 /* Send the charges and/or c6/sigmas to our PME only node */
9937 gmx_pme_send_parameters(cr,
9939 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9940 mdatoms->chargeA, mdatoms->chargeB,
9941 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9942 mdatoms->sigmaA, mdatoms->sigmaB,
9943 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9948 set_constraints(constr, top_local, ir, mdatoms, cr);
9953 /* Update the local pull groups */
9954 dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
9959 /* Update the local rotation groups */
9960 dd_make_local_rotation_groups(dd, ir->rot);
9963 if (ir->eSwapCoords != eswapNO)
9965 /* Update the local groups needed for ion swapping */
9966 dd_make_local_swap_groups(dd, ir->swap);
9969 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9970 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9972 add_dd_statistics(dd);
9974 /* Make sure we only count the cycles for this DD partitioning */
9975 clear_dd_cycle_counts(dd);
9977 /* Because the order of the atoms might have changed since
9978 * the last vsite construction, we need to communicate the constructing
9979 * atom coordinates again (for spreading the forces this MD step).
9981 dd_move_x_vsites(dd, state_local->box, state_local->x);
9983 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9985 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9987 dd_move_x(dd, state_local->box, state_local->x);
9988 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9989 -1, state_local->x, state_local->box);
9992 /* Store the partitioning step */
9993 comm->partition_step = step;
9995 /* Increase the DD partitioning counter */
9997 /* The state currently matches this DD partitioning count, store it */
9998 state_local->ddp_count = dd->ddp_count;
10001 /* The DD master node knows the complete cg distribution,
10002 * store the count so we can possibly skip the cg info communication.
10004 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
10007 if (comm->DD_debug > 0)
10009 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
10010 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
10011 "after partitioning");
10014 wallcycle_stop(wcycle, ewcDOMDEC);