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 /* Should only be called on the DD master rank */
5584 assert(DDMASTER(dd));
5586 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5588 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5596 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5601 flags = dd_load_flags(dd);
5605 "DD load balancing is limited by minimum cell size in dimension");
5606 for (d = 0; d < dd->ndim; d++)
5610 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5613 fprintf(fplog, "\n");
5615 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5616 if (dd->comm->bDynLoadBal)
5618 fprintf(fplog, " vol min/aver %5.3f%c",
5619 dd_vol_min(dd), flags ? '!' : ' ');
5623 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5625 if (dd->comm->cycl_n[ddCyclPME])
5627 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5629 fprintf(fplog, "\n\n");
5632 static void dd_print_load_verbose(gmx_domdec_t *dd)
5634 if (dd->comm->bDynLoadBal)
5636 fprintf(stderr, "vol %4.2f%c ",
5637 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5641 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5643 if (dd->comm->cycl_n[ddCyclPME])
5645 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5650 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5655 gmx_domdec_root_t *root;
5656 gmx_bool bPartOfGroup = FALSE;
5658 dim = dd->dim[dim_ind];
5659 copy_ivec(loc, loc_c);
5660 for (i = 0; i < dd->nc[dim]; i++)
5663 rank = dd_index(dd->nc, loc_c);
5664 if (rank == dd->rank)
5666 /* This process is part of the group */
5667 bPartOfGroup = TRUE;
5670 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5674 dd->comm->mpi_comm_load[dim_ind] = c_row;
5675 if (dd->comm->eDLB != edlbNO)
5677 if (dd->ci[dim] == dd->master_ci[dim])
5679 /* This is the root process of this row */
5680 snew(dd->comm->root[dim_ind], 1);
5681 root = dd->comm->root[dim_ind];
5682 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5683 snew(root->old_cell_f, dd->nc[dim]+1);
5684 snew(root->bCellMin, dd->nc[dim]);
5687 snew(root->cell_f_max0, dd->nc[dim]);
5688 snew(root->cell_f_min1, dd->nc[dim]);
5689 snew(root->bound_min, dd->nc[dim]);
5690 snew(root->bound_max, dd->nc[dim]);
5692 snew(root->buf_ncd, dd->nc[dim]);
5696 /* This is not a root process, we only need to receive cell_f */
5697 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5700 if (dd->ci[dim] == dd->master_ci[dim])
5702 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5708 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5709 const gmx_hw_info_t gmx_unused *hwinfo,
5710 const gmx_hw_opt_t gmx_unused *hw_opt)
5713 int physicalnode_id_hash;
5716 MPI_Comm mpi_comm_pp_physicalnode;
5718 if (!(cr->duty & DUTY_PP) || hw_opt->gpu_opt.n_dev_use == 0)
5720 /* Only PP nodes (currently) use GPUs.
5721 * If we don't have GPUs, there are no resources to share.
5726 physicalnode_id_hash = gmx_physicalnode_id_hash();
5728 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5734 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5735 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5736 dd->rank, physicalnode_id_hash, gpu_id);
5738 /* Split the PP communicator over the physical nodes */
5739 /* TODO: See if we should store this (before), as it's also used for
5740 * for the nodecomm summution.
5742 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5743 &mpi_comm_pp_physicalnode);
5744 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5745 &dd->comm->mpi_comm_gpu_shared);
5746 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5747 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5751 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5754 /* Note that some ranks could share a GPU, while others don't */
5756 if (dd->comm->nrank_gpu_shared == 1)
5758 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5763 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5766 int dim0, dim1, i, j;
5771 fprintf(debug, "Making load communicators\n");
5774 snew(dd->comm->load, std::max(dd->ndim, 1));
5775 snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
5783 make_load_communicator(dd, 0, loc);
5787 for (i = 0; i < dd->nc[dim0]; i++)
5790 make_load_communicator(dd, 1, loc);
5796 for (i = 0; i < dd->nc[dim0]; i++)
5800 for (j = 0; j < dd->nc[dim1]; j++)
5803 make_load_communicator(dd, 2, loc);
5810 fprintf(debug, "Finished making load communicators\n");
5815 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5817 int d, dim, i, j, m;
5820 ivec dd_zp[DD_MAXIZONE];
5821 gmx_domdec_zones_t *zones;
5822 gmx_domdec_ns_ranges_t *izone;
5824 for (d = 0; d < dd->ndim; d++)
5827 copy_ivec(dd->ci, tmp);
5828 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5829 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5830 copy_ivec(dd->ci, tmp);
5831 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5832 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5835 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5838 dd->neighbor[d][1]);
5844 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5846 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5847 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5854 for (i = 0; i < nzonep; i++)
5856 copy_ivec(dd_zp3[i], dd_zp[i]);
5862 for (i = 0; i < nzonep; i++)
5864 copy_ivec(dd_zp2[i], dd_zp[i]);
5870 for (i = 0; i < nzonep; i++)
5872 copy_ivec(dd_zp1[i], dd_zp[i]);
5878 for (i = 0; i < nzonep; i++)
5880 copy_ivec(dd_zp0[i], dd_zp[i]);
5884 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5889 zones = &dd->comm->zones;
5891 for (i = 0; i < nzone; i++)
5894 clear_ivec(zones->shift[i]);
5895 for (d = 0; d < dd->ndim; d++)
5897 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5902 for (i = 0; i < nzone; i++)
5904 for (d = 0; d < DIM; d++)
5906 s[d] = dd->ci[d] - zones->shift[i][d];
5911 else if (s[d] >= dd->nc[d])
5917 zones->nizone = nzonep;
5918 for (i = 0; i < zones->nizone; i++)
5920 if (dd_zp[i][0] != i)
5922 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5924 izone = &zones->izone[i];
5925 izone->j0 = dd_zp[i][1];
5926 izone->j1 = dd_zp[i][2];
5927 for (dim = 0; dim < DIM; dim++)
5929 if (dd->nc[dim] == 1)
5931 /* All shifts should be allowed */
5932 izone->shift0[dim] = -1;
5933 izone->shift1[dim] = 1;
5938 izone->shift0[d] = 0;
5939 izone->shift1[d] = 0;
5940 for(j=izone->j0; j<izone->j1; j++) {
5941 if (dd->shift[j][d] > dd->shift[i][d])
5942 izone->shift0[d] = -1;
5943 if (dd->shift[j][d] < dd->shift[i][d])
5944 izone->shift1[d] = 1;
5950 /* Assume the shift are not more than 1 cell */
5951 izone->shift0[dim] = 1;
5952 izone->shift1[dim] = -1;
5953 for (j = izone->j0; j < izone->j1; j++)
5955 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5956 if (shift_diff < izone->shift0[dim])
5958 izone->shift0[dim] = shift_diff;
5960 if (shift_diff > izone->shift1[dim])
5962 izone->shift1[dim] = shift_diff;
5969 if (dd->comm->eDLB != edlbNO)
5971 snew(dd->comm->root, dd->ndim);
5974 if (dd->comm->bRecordLoad)
5976 make_load_communicators(dd);
5980 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5986 gmx_domdec_comm_t *comm;
5993 if (comm->bCartesianPP)
5995 /* Set up cartesian communication for the particle-particle part */
5998 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5999 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
6002 for (int i = 0; i < DIM; i++)
6006 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
6008 /* We overwrite the old communicator with the new cartesian one */
6009 cr->mpi_comm_mygroup = comm_cart;
6012 dd->mpi_comm_all = cr->mpi_comm_mygroup;
6013 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
6015 if (comm->bCartesianPP_PME)
6017 /* Since we want to use the original cartesian setup for sim,
6018 * and not the one after split, we need to make an index.
6020 snew(comm->ddindex2ddnodeid, dd->nnodes);
6021 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
6022 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
6023 /* Get the rank of the DD master,
6024 * above we made sure that the master node is a PP node.
6034 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
6036 else if (comm->bCartesianPP)
6038 if (cr->npmenodes == 0)
6040 /* The PP communicator is also
6041 * the communicator for this simulation
6043 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6045 cr->nodeid = dd->rank;
6047 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6049 /* We need to make an index to go from the coordinates
6050 * to the nodeid of this simulation.
6052 snew(comm->ddindex2simnodeid, dd->nnodes);
6053 snew(buf, dd->nnodes);
6054 if (cr->duty & DUTY_PP)
6056 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6058 /* Communicate the ddindex to simulation nodeid index */
6059 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6060 cr->mpi_comm_mysim);
6063 /* Determine the master coordinates and rank.
6064 * The DD master should be the same node as the master of this sim.
6066 for (int i = 0; i < dd->nnodes; i++)
6068 if (comm->ddindex2simnodeid[i] == 0)
6070 ddindex2xyz(dd->nc, i, dd->master_ci);
6071 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6076 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6081 /* No Cartesian communicators */
6082 /* We use the rank in dd->comm->all as DD index */
6083 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6084 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6086 clear_ivec(dd->master_ci);
6093 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6094 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6099 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6100 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6104 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6108 gmx_domdec_comm_t *comm;
6113 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6116 snew(comm->ddindex2simnodeid, dd->nnodes);
6117 snew(buf, dd->nnodes);
6118 if (cr->duty & DUTY_PP)
6120 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6122 /* Communicate the ddindex to simulation nodeid index */
6123 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6124 cr->mpi_comm_mysim);
6130 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6131 int ncg, int natoms)
6133 gmx_domdec_master_t *ma;
6138 snew(ma->ncg, dd->nnodes);
6139 snew(ma->index, dd->nnodes+1);
6141 snew(ma->nat, dd->nnodes);
6142 snew(ma->ibuf, dd->nnodes*2);
6143 snew(ma->cell_x, DIM);
6144 for (i = 0; i < DIM; i++)
6146 snew(ma->cell_x[i], dd->nc[i]+1);
6149 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6155 snew(ma->vbuf, natoms);
6161 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6162 int gmx_unused reorder)
6165 gmx_domdec_comm_t *comm;
6175 if (comm->bCartesianPP)
6177 for (i = 1; i < DIM; i++)
6179 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6181 if (bDiv[YY] || bDiv[ZZ])
6183 comm->bCartesianPP_PME = TRUE;
6184 /* If we have 2D PME decomposition, which is always in x+y,
6185 * we stack the PME only nodes in z.
6186 * Otherwise we choose the direction that provides the thinnest slab
6187 * of PME only nodes as this will have the least effect
6188 * on the PP communication.
6189 * But for the PME communication the opposite might be better.
6191 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6193 dd->nc[YY] > dd->nc[ZZ]))
6195 comm->cartpmedim = ZZ;
6199 comm->cartpmedim = YY;
6201 comm->ntot[comm->cartpmedim]
6202 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6206 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]);
6208 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6213 if (comm->bCartesianPP_PME)
6220 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]);
6223 for (i = 0; i < DIM; i++)
6227 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6229 MPI_Comm_rank(comm_cart, &rank);
6230 if (MASTERNODE(cr) && rank != 0)
6232 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6235 /* With this assigment we loose the link to the original communicator
6236 * which will usually be MPI_COMM_WORLD, unless have multisim.
6238 cr->mpi_comm_mysim = comm_cart;
6239 cr->sim_nodeid = rank;
6241 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6245 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6246 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6249 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6253 if (cr->npmenodes == 0 ||
6254 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6256 cr->duty = DUTY_PME;
6259 /* Split the sim communicator into PP and PME only nodes */
6260 MPI_Comm_split(cr->mpi_comm_mysim,
6262 dd_index(comm->ntot, dd->ci),
6263 &cr->mpi_comm_mygroup);
6267 switch (dd_node_order)
6272 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6275 case ddnoINTERLEAVE:
6276 /* Interleave the PP-only and PME-only nodes,
6277 * as on clusters with dual-core machines this will double
6278 * the communication bandwidth of the PME processes
6279 * and thus speed up the PP <-> PME and inter PME communication.
6283 fprintf(fplog, "Interleaving PP and PME ranks\n");
6285 comm->pmenodes = dd_pmenodes(cr);
6290 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6293 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6295 cr->duty = DUTY_PME;
6302 /* Split the sim communicator into PP and PME only nodes */
6303 MPI_Comm_split(cr->mpi_comm_mysim,
6306 &cr->mpi_comm_mygroup);
6307 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6313 fprintf(fplog, "This rank does only %s work.\n\n",
6314 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6318 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6321 gmx_domdec_comm_t *comm;
6327 copy_ivec(dd->nc, comm->ntot);
6329 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6330 comm->bCartesianPP_PME = FALSE;
6332 /* Reorder the nodes by default. This might change the MPI ranks.
6333 * Real reordering is only supported on very few architectures,
6334 * Blue Gene is one of them.
6336 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6338 if (cr->npmenodes > 0)
6340 /* Split the communicator into a PP and PME part */
6341 split_communicator(fplog, cr, dd_node_order, CartReorder);
6342 if (comm->bCartesianPP_PME)
6344 /* We (possibly) reordered the nodes in split_communicator,
6345 * so it is no longer required in make_pp_communicator.
6347 CartReorder = FALSE;
6352 /* All nodes do PP and PME */
6354 /* We do not require separate communicators */
6355 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6359 if (cr->duty & DUTY_PP)
6361 /* Copy or make a new PP communicator */
6362 make_pp_communicator(fplog, cr, CartReorder);
6366 receive_ddindex2simnodeid(cr);
6369 if (!(cr->duty & DUTY_PME))
6371 /* Set up the commnuication to our PME node */
6372 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6373 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6376 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6377 dd->pme_nodeid, dd->pme_receive_vir_ener);
6382 dd->pme_nodeid = -1;
6387 dd->ma = init_gmx_domdec_master_t(dd,
6389 comm->cgs_gl.index[comm->cgs_gl.nr]);
6393 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6395 real *slb_frac, tot;
6400 if (nc > 1 && size_string != NULL)
6404 fprintf(fplog, "Using static load balancing for the %s direction\n",
6409 for (i = 0; i < nc; i++)
6412 sscanf(size_string, "%20lf%n", &dbl, &n);
6415 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6424 fprintf(fplog, "Relative cell sizes:");
6426 for (i = 0; i < nc; i++)
6431 fprintf(fplog, " %5.3f", slb_frac[i]);
6436 fprintf(fplog, "\n");
6443 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6446 gmx_mtop_ilistloop_t iloop;
6450 iloop = gmx_mtop_ilistloop_init(mtop);
6451 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6453 for (ftype = 0; ftype < F_NRE; ftype++)
6455 if ((interaction_function[ftype].flags & IF_BOND) &&
6458 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6466 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6472 val = getenv(env_var);
6475 if (sscanf(val, "%20d", &nst) <= 0)
6481 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6489 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6493 fprintf(stderr, "\n%s\n", warn_string);
6497 fprintf(fplog, "\n%s\n", warn_string);
6501 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6502 t_inputrec *ir, FILE *fplog)
6504 if (ir->ePBC == epbcSCREW &&
6505 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6507 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6510 if (ir->ns_type == ensSIMPLE)
6512 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6515 if (ir->nstlist == 0)
6517 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6520 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6522 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6526 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6531 r = ddbox->box_size[XX];
6532 for (di = 0; di < dd->ndim; di++)
6535 /* Check using the initial average cell size */
6536 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6542 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6543 const char *dlb_opt, gmx_bool bRecordLoad,
6544 unsigned long Flags, t_inputrec *ir)
6551 case 'a': eDLB = edlbAUTO; break;
6552 case 'n': eDLB = edlbNO; break;
6553 case 'y': eDLB = edlbYES; break;
6554 default: gmx_incons("Unknown dlb_opt");
6557 if (Flags & MD_RERUN)
6562 if (!EI_DYNAMICS(ir->eI))
6564 if (eDLB == edlbYES)
6566 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6567 dd_warning(cr, fplog, buf);
6575 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6580 if (Flags & MD_REPRODUCIBLE)
6587 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6591 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6594 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6602 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6607 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6609 /* Decomposition order z,y,x */
6612 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6614 for (dim = DIM-1; dim >= 0; dim--)
6616 if (dd->nc[dim] > 1)
6618 dd->dim[dd->ndim++] = dim;
6624 /* Decomposition order x,y,z */
6625 for (dim = 0; dim < DIM; dim++)
6627 if (dd->nc[dim] > 1)
6629 dd->dim[dd->ndim++] = dim;
6635 static gmx_domdec_comm_t *init_dd_comm()
6637 gmx_domdec_comm_t *comm;
6641 snew(comm->cggl_flag, DIM*2);
6642 snew(comm->cgcm_state, DIM*2);
6643 for (i = 0; i < DIM*2; i++)
6645 comm->cggl_flag_nalloc[i] = 0;
6646 comm->cgcm_state_nalloc[i] = 0;
6649 comm->nalloc_int = 0;
6650 comm->buf_int = NULL;
6652 vec_rvec_init(&comm->vbuf);
6654 comm->n_load_have = 0;
6655 comm->n_load_collect = 0;
6657 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6659 comm->sum_nat[i] = 0;
6663 comm->load_step = 0;
6666 clear_ivec(comm->load_lim);
6673 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6674 unsigned long Flags,
6676 real comm_distance_min, real rconstr,
6677 const char *dlb_opt, real dlb_scale,
6678 const char *sizex, const char *sizey, const char *sizez,
6679 gmx_mtop_t *mtop, t_inputrec *ir,
6680 matrix box, rvec *x,
6682 int *npme_x, int *npme_y)
6685 gmx_domdec_comm_t *comm;
6687 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6690 const real tenPercentMargin = 1.1;
6695 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6700 dd->comm = init_dd_comm();
6702 snew(comm->cggl_flag, DIM*2);
6703 snew(comm->cgcm_state, DIM*2);
6705 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6706 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6708 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6709 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6710 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6711 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6712 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6713 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6714 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6715 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6717 dd->pme_recv_f_alloc = 0;
6718 dd->pme_recv_f_buf = NULL;
6720 if (dd->bSendRecv2 && fplog)
6722 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");
6728 fprintf(fplog, "Will load balance based on FLOP count\n");
6730 if (comm->eFlop > 1)
6732 srand(1+cr->nodeid);
6734 comm->bRecordLoad = TRUE;
6738 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6742 /* Initialize to GPU share count to 0, might change later */
6743 comm->nrank_gpu_shared = 0;
6745 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6746 comm->bDLB_locked = FALSE;
6747 comm->bCheckWhetherToTurnDlbOn = TRUE;
6749 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6752 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6754 dd->bGridJump = comm->bDynLoadBal;
6755 comm->bPMELoadBalDLBLimits = FALSE;
6757 if (comm->nstSortCG)
6761 if (comm->nstSortCG == 1)
6763 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6767 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6771 snew(comm->sort, 1);
6777 fprintf(fplog, "Will not sort the charge groups\n");
6781 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6783 comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
6784 mtop->bIntermolecularInteractions);
6785 if (comm->bInterCGBondeds)
6787 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6791 comm->bInterCGMultiBody = FALSE;
6794 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6795 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6797 if (ir->rlistlong == 0)
6799 /* Set the cut-off to some very large value,
6800 * so we don't need if statements everywhere in the code.
6801 * We use sqrt, since the cut-off is squared in some places.
6803 comm->cutoff = GMX_CUTOFF_INF;
6807 comm->cutoff = ir->rlistlong;
6809 comm->cutoff_mbody = 0;
6811 comm->cellsize_limit = 0;
6812 comm->bBondComm = FALSE;
6814 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6815 * within nstlist steps. Since boundaries are allowed to displace by half
6816 * a cell size, DD cells should be at least the size of the list buffer.
6818 comm->cellsize_limit = std::max(comm->cellsize_limit,
6819 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6821 if (comm->bInterCGBondeds)
6823 if (comm_distance_min > 0)
6825 comm->cutoff_mbody = comm_distance_min;
6826 if (Flags & MD_DDBONDCOMM)
6828 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6832 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6834 r_bonded_limit = comm->cutoff_mbody;
6836 else if (ir->bPeriodicMols)
6838 /* Can not easily determine the required cut-off */
6839 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");
6840 comm->cutoff_mbody = comm->cutoff/2;
6841 r_bonded_limit = comm->cutoff_mbody;
6847 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6848 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6850 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6851 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6853 /* We use an initial margin of 10% for the minimum cell size,
6854 * except when we are just below the non-bonded cut-off.
6856 if (Flags & MD_DDBONDCOMM)
6858 if (std::max(r_2b, r_mb) > comm->cutoff)
6860 r_bonded = std::max(r_2b, r_mb);
6861 r_bonded_limit = tenPercentMargin*r_bonded;
6862 comm->bBondComm = TRUE;
6867 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6869 /* We determine cutoff_mbody later */
6873 /* No special bonded communication,
6874 * simply increase the DD cut-off.
6876 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6877 comm->cutoff_mbody = r_bonded_limit;
6878 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6884 "Minimum cell size due to bonded interactions: %.3f nm\n",
6887 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6890 if (dd->bInterCGcons && rconstr <= 0)
6892 /* There is a cell size limit due to the constraints (P-LINCS) */
6893 rconstr = constr_r_max(fplog, mtop, ir);
6897 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6899 if (rconstr > comm->cellsize_limit)
6901 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6905 else if (rconstr > 0 && fplog)
6907 /* Here we do not check for dd->bInterCGcons,
6908 * because one can also set a cell size limit for virtual sites only
6909 * and at this point we don't know yet if there are intercg v-sites.
6912 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6915 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6917 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6921 copy_ivec(nc, dd->nc);
6922 set_dd_dim(fplog, dd);
6923 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6925 if (cr->npmenodes == -1)
6929 acs = average_cellsize_min(dd, ddbox);
6930 if (acs < comm->cellsize_limit)
6934 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6936 gmx_fatal_collective(FARGS, cr, NULL,
6937 "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",
6938 acs, comm->cellsize_limit);
6943 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6945 /* We need to choose the optimal DD grid and possibly PME nodes */
6946 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6947 comm->eDLB != edlbNO, dlb_scale,
6948 comm->cellsize_limit, comm->cutoff,
6949 comm->bInterCGBondeds);
6951 if (dd->nc[XX] == 0)
6953 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6954 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6955 !bC ? "-rdd" : "-rcon",
6956 comm->eDLB != edlbNO ? " or -dds" : "",
6957 bC ? " or your LINCS settings" : "");
6959 gmx_fatal_collective(FARGS, cr, NULL,
6960 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6962 "Look in the log file for details on the domain decomposition",
6963 cr->nnodes-cr->npmenodes, limit, buf);
6965 set_dd_dim(fplog, dd);
6971 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6972 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6975 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6976 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6978 gmx_fatal_collective(FARGS, cr, NULL,
6979 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6980 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6982 if (cr->npmenodes > dd->nnodes)
6984 gmx_fatal_collective(FARGS, cr, NULL,
6985 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6987 if (cr->npmenodes > 0)
6989 comm->npmenodes = cr->npmenodes;
6993 comm->npmenodes = dd->nnodes;
6996 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6998 /* The following choices should match those
6999 * in comm_cost_est in domdec_setup.c.
7000 * Note that here the checks have to take into account
7001 * that the decomposition might occur in a different order than xyz
7002 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
7003 * in which case they will not match those in comm_cost_est,
7004 * but since that is mainly for testing purposes that's fine.
7006 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
7007 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
7008 getenv("GMX_PMEONEDD") == NULL)
7010 comm->npmedecompdim = 2;
7011 comm->npmenodes_x = dd->nc[XX];
7012 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
7016 /* In case nc is 1 in both x and y we could still choose to
7017 * decompose pme in y instead of x, but we use x for simplicity.
7019 comm->npmedecompdim = 1;
7020 if (dd->dim[0] == YY)
7022 comm->npmenodes_x = 1;
7023 comm->npmenodes_y = comm->npmenodes;
7027 comm->npmenodes_x = comm->npmenodes;
7028 comm->npmenodes_y = 1;
7033 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
7034 comm->npmenodes_x, comm->npmenodes_y, 1);
7039 comm->npmedecompdim = 0;
7040 comm->npmenodes_x = 0;
7041 comm->npmenodes_y = 0;
7044 /* Technically we don't need both of these,
7045 * but it simplifies code not having to recalculate it.
7047 *npme_x = comm->npmenodes_x;
7048 *npme_y = comm->npmenodes_y;
7050 snew(comm->slb_frac, DIM);
7051 if (comm->eDLB == edlbNO)
7053 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7054 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7055 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7058 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7060 if (comm->bBondComm || comm->eDLB != edlbNO)
7062 /* Set the bonded communication distance to halfway
7063 * the minimum and the maximum,
7064 * since the extra communication cost is nearly zero.
7066 acs = average_cellsize_min(dd, ddbox);
7067 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7068 if (comm->eDLB != edlbNO)
7070 /* Check if this does not limit the scaling */
7071 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7073 if (!comm->bBondComm)
7075 /* Without bBondComm do not go beyond the n.b. cut-off */
7076 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7077 if (comm->cellsize_limit >= comm->cutoff)
7079 /* We don't loose a lot of efficieny
7080 * when increasing it to the n.b. cut-off.
7081 * It can even be slightly faster, because we need
7082 * less checks for the communication setup.
7084 comm->cutoff_mbody = comm->cutoff;
7087 /* Check if we did not end up below our original limit */
7088 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7090 if (comm->cutoff_mbody > comm->cellsize_limit)
7092 comm->cellsize_limit = comm->cutoff_mbody;
7095 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7100 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7101 "cellsize limit %f\n",
7102 comm->bBondComm, comm->cellsize_limit);
7107 check_dd_restrictions(cr, dd, ir, fplog);
7110 comm->partition_step = INT_MIN;
7113 clear_dd_cycle_counts(dd);
7118 static void set_dlb_limits(gmx_domdec_t *dd)
7123 for (d = 0; d < dd->ndim; d++)
7125 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7126 dd->comm->cellsize_min[dd->dim[d]] =
7127 dd->comm->cellsize_min_dlb[dd->dim[d]];
7132 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7135 gmx_domdec_comm_t *comm;
7145 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);
7148 cellsize_min = comm->cellsize_min[dd->dim[0]];
7149 for (d = 1; d < dd->ndim; d++)
7151 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7154 if (cellsize_min < comm->cellsize_limit*1.05)
7156 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");
7158 /* Change DLB from "auto" to "no". */
7159 comm->eDLB = edlbNO;
7164 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7165 comm->bDynLoadBal = TRUE;
7166 dd->bGridJump = TRUE;
7170 /* We can set the required cell size info here,
7171 * so we do not need to communicate this.
7172 * The grid is completely uniform.
7174 for (d = 0; d < dd->ndim; d++)
7178 comm->load[d].sum_m = comm->load[d].sum;
7180 nc = dd->nc[dd->dim[d]];
7181 for (i = 0; i < nc; i++)
7183 comm->root[d]->cell_f[i] = i/(real)nc;
7186 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7187 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7190 comm->root[d]->cell_f[nc] = 1.0;
7195 static char *init_bLocalCG(gmx_mtop_t *mtop)
7200 ncg = ncg_mtop(mtop);
7201 snew(bLocalCG, ncg);
7202 for (cg = 0; cg < ncg; cg++)
7204 bLocalCG[cg] = FALSE;
7210 void dd_init_bondeds(FILE *fplog,
7211 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7213 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7215 gmx_domdec_comm_t *comm;
7217 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7221 if (comm->bBondComm)
7223 /* Communicate atoms beyond the cut-off for bonded interactions */
7226 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7228 comm->bLocalCG = init_bLocalCG(mtop);
7232 /* Only communicate atoms based on cut-off */
7233 comm->cglink = NULL;
7234 comm->bLocalCG = NULL;
7238 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7240 gmx_bool bDynLoadBal, real dlb_scale,
7243 gmx_domdec_comm_t *comm;
7258 fprintf(fplog, "The maximum number of communication pulses is:");
7259 for (d = 0; d < dd->ndim; d++)
7261 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7263 fprintf(fplog, "\n");
7264 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7265 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7266 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7267 for (d = 0; d < DIM; d++)
7271 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7278 comm->cellsize_min_dlb[d]/
7279 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7281 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7284 fprintf(fplog, "\n");
7288 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7289 fprintf(fplog, "The initial number of communication pulses is:");
7290 for (d = 0; d < dd->ndim; d++)
7292 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7294 fprintf(fplog, "\n");
7295 fprintf(fplog, "The initial domain decomposition cell size is:");
7296 for (d = 0; d < DIM; d++)
7300 fprintf(fplog, " %c %.2f nm",
7301 dim2char(d), dd->comm->cellsize_min[d]);
7304 fprintf(fplog, "\n\n");
7307 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7309 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7310 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7311 "non-bonded interactions", "", comm->cutoff);
7315 limit = dd->comm->cellsize_limit;
7319 if (dynamic_dd_box(ddbox, ir))
7321 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7323 limit = dd->comm->cellsize_min[XX];
7324 for (d = 1; d < DIM; d++)
7326 limit = std::min(limit, dd->comm->cellsize_min[d]);
7330 if (comm->bInterCGBondeds)
7332 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7333 "two-body bonded interactions", "(-rdd)",
7334 std::max(comm->cutoff, comm->cutoff_mbody));
7335 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7336 "multi-body bonded interactions", "(-rdd)",
7337 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7341 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7342 "virtual site constructions", "(-rcon)", limit);
7344 if (dd->constraint_comm)
7346 sprintf(buf, "atoms separated by up to %d constraints",
7348 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7349 buf, "(-rcon)", limit);
7351 fprintf(fplog, "\n");
7357 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7359 const t_inputrec *ir,
7360 const gmx_ddbox_t *ddbox)
7362 gmx_domdec_comm_t *comm;
7363 int d, dim, npulse, npulse_d_max, npulse_d;
7368 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7370 /* Determine the maximum number of comm. pulses in one dimension */
7372 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7374 /* Determine the maximum required number of grid pulses */
7375 if (comm->cellsize_limit >= comm->cutoff)
7377 /* Only a single pulse is required */
7380 else if (!bNoCutOff && comm->cellsize_limit > 0)
7382 /* We round down slightly here to avoid overhead due to the latency
7383 * of extra communication calls when the cut-off
7384 * would be only slightly longer than the cell size.
7385 * Later cellsize_limit is redetermined,
7386 * so we can not miss interactions due to this rounding.
7388 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7392 /* There is no cell size limit */
7393 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7396 if (!bNoCutOff && npulse > 1)
7398 /* See if we can do with less pulses, based on dlb_scale */
7400 for (d = 0; d < dd->ndim; d++)
7403 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7404 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7405 npulse_d_max = std::max(npulse_d_max, npulse_d);
7407 npulse = std::min(npulse, npulse_d_max);
7410 /* This env var can override npulse */
7411 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7418 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7419 for (d = 0; d < dd->ndim; d++)
7421 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7422 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7423 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7424 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7425 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7427 comm->bVacDLBNoLimit = FALSE;
7431 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7432 if (!comm->bVacDLBNoLimit)
7434 comm->cellsize_limit = std::max(comm->cellsize_limit,
7435 comm->cutoff/comm->maxpulse);
7437 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7438 /* Set the minimum cell size for each DD dimension */
7439 for (d = 0; d < dd->ndim; d++)
7441 if (comm->bVacDLBNoLimit ||
7442 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7444 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7448 comm->cellsize_min_dlb[dd->dim[d]] =
7449 comm->cutoff/comm->cd[d].np_dlb;
7452 if (comm->cutoff_mbody <= 0)
7454 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7456 if (comm->bDynLoadBal)
7462 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7464 /* If each molecule is a single charge group
7465 * or we use domain decomposition for each periodic dimension,
7466 * we do not need to take pbc into account for the bonded interactions.
7468 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7471 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7474 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7475 t_inputrec *ir, gmx_ddbox_t *ddbox)
7477 gmx_domdec_comm_t *comm;
7483 /* Initialize the thread data.
7484 * This can not be done in init_domain_decomposition,
7485 * as the numbers of threads is determined later.
7487 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7490 snew(comm->dth, comm->nth);
7493 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7495 init_ddpme(dd, &comm->ddpme[0], 0);
7496 if (comm->npmedecompdim >= 2)
7498 init_ddpme(dd, &comm->ddpme[1], 1);
7503 comm->npmenodes = 0;
7504 if (dd->pme_nodeid >= 0)
7506 gmx_fatal_collective(FARGS, NULL, dd,
7507 "Can not have separate PME ranks without PME electrostatics");
7513 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7515 if (comm->eDLB != edlbNO)
7517 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7520 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7521 if (comm->eDLB == edlbAUTO)
7525 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7527 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7530 if (ir->ePBC == epbcNONE)
7532 vol_frac = 1 - 1/(double)dd->nnodes;
7537 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7541 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7543 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7545 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7548 static gmx_bool test_dd_cutoff(t_commrec *cr,
7549 t_state *state, t_inputrec *ir,
7560 set_ddbox(dd, FALSE, cr, ir, state->box,
7561 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7565 for (d = 0; d < dd->ndim; d++)
7569 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7570 if (dynamic_dd_box(&ddbox, ir))
7572 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7575 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7577 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7578 dd->comm->cd[d].np_dlb > 0)
7580 if (np > dd->comm->cd[d].np_dlb)
7585 /* If a current local cell size is smaller than the requested
7586 * cut-off, we could still fix it, but this gets very complicated.
7587 * Without fixing here, we might actually need more checks.
7589 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7596 if (dd->comm->eDLB != edlbNO)
7598 /* If DLB is not active yet, we don't need to check the grid jumps.
7599 * Actually we shouldn't, because then the grid jump data is not set.
7601 if (dd->comm->bDynLoadBal &&
7602 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7607 gmx_sumi(1, &LocallyLimited, cr);
7609 if (LocallyLimited > 0)
7618 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7621 gmx_bool bCutoffAllowed;
7623 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7627 cr->dd->comm->cutoff = cutoff_req;
7630 return bCutoffAllowed;
7633 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
7635 gmx_domdec_comm_t *comm;
7637 comm = cr->dd->comm;
7639 /* Turn on the DLB limiting (might have been on already) */
7640 comm->bPMELoadBalDLBLimits = TRUE;
7642 /* Change the cut-off limit */
7643 comm->PMELoadBal_max_cutoff = cutoff;
7647 fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
7648 comm->PMELoadBal_max_cutoff);
7652 /* Sets whether we should later check the load imbalance data, so that
7653 * we can trigger dynamic load balancing if enough imbalance has
7656 * Used after PME load balancing unlocks DLB, so that the check
7657 * whether DLB will be useful can happen immediately.
7659 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
7661 if (dd->comm->eDLB == edlbAUTO)
7663 assert(!dd_dlb_is_locked(dd));
7665 dd->comm->bCheckWhetherToTurnDlbOn = bValue;
7669 /* Returns if we should check whether there has been enough load
7670 * imbalance to trigger dynamic load balancing.
7672 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
7674 const int nddp_chk_dlb = 100;
7676 if (dd->comm->eDLB != edlbAUTO || dd_dlb_is_locked(dd))
7681 /* We should check whether we should use DLB directly after
7683 if (dd->comm->bCheckWhetherToTurnDlbOn)
7685 /* This flag was set when the PME load-balancing routines
7686 unlocked DLB, and should now be cleared. */
7687 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
7690 /* We should also check whether we should use DLB every 100
7691 * partitionings (we do not do this every partioning, so that we
7692 * avoid excessive communication). */
7693 if (dd->comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1)
7701 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
7703 return dd->comm->bDynLoadBal;
7706 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7708 return dd->comm->bDLB_locked;
7711 void dd_dlb_lock(gmx_domdec_t *dd)
7713 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7714 if (dd->comm->eDLB == edlbAUTO)
7716 dd->comm->bDLB_locked = TRUE;
7720 void dd_dlb_unlock(gmx_domdec_t *dd)
7722 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7723 if (dd->comm->eDLB == edlbAUTO)
7725 dd->comm->bDLB_locked = FALSE;
7726 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, !dd->comm->bDynLoadBal);
7730 static void merge_cg_buffers(int ncell,
7731 gmx_domdec_comm_dim_t *cd, int pulse,
7733 int *index_gl, int *recv_i,
7734 rvec *cg_cm, rvec *recv_vr,
7736 cginfo_mb_t *cginfo_mb, int *cginfo)
7738 gmx_domdec_ind_t *ind, *ind_p;
7739 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7740 int shift, shift_at;
7742 ind = &cd->ind[pulse];
7744 /* First correct the already stored data */
7745 shift = ind->nrecv[ncell];
7746 for (cell = ncell-1; cell >= 0; cell--)
7748 shift -= ind->nrecv[cell];
7751 /* Move the cg's present from previous grid pulses */
7752 cg0 = ncg_cell[ncell+cell];
7753 cg1 = ncg_cell[ncell+cell+1];
7754 cgindex[cg1+shift] = cgindex[cg1];
7755 for (cg = cg1-1; cg >= cg0; cg--)
7757 index_gl[cg+shift] = index_gl[cg];
7758 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7759 cgindex[cg+shift] = cgindex[cg];
7760 cginfo[cg+shift] = cginfo[cg];
7762 /* Correct the already stored send indices for the shift */
7763 for (p = 1; p <= pulse; p++)
7765 ind_p = &cd->ind[p];
7767 for (c = 0; c < cell; c++)
7769 cg0 += ind_p->nsend[c];
7771 cg1 = cg0 + ind_p->nsend[cell];
7772 for (cg = cg0; cg < cg1; cg++)
7774 ind_p->index[cg] += shift;
7780 /* Merge in the communicated buffers */
7784 for (cell = 0; cell < ncell; cell++)
7786 cg1 = ncg_cell[ncell+cell+1] + shift;
7789 /* Correct the old cg indices */
7790 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7792 cgindex[cg+1] += shift_at;
7795 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7797 /* Copy this charge group from the buffer */
7798 index_gl[cg1] = recv_i[cg0];
7799 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7800 /* Add it to the cgindex */
7801 cg_gl = index_gl[cg1];
7802 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7803 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7804 cgindex[cg1+1] = cgindex[cg1] + nat;
7809 shift += ind->nrecv[cell];
7810 ncg_cell[ncell+cell+1] = cg1;
7814 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7815 int nzone, int cg0, const int *cgindex)
7819 /* Store the atom block boundaries for easy copying of communication buffers
7822 for (zone = 0; zone < nzone; zone++)
7824 for (p = 0; p < cd->np; p++)
7826 cd->ind[p].cell2at0[zone] = cgindex[cg];
7827 cg += cd->ind[p].nrecv[zone];
7828 cd->ind[p].cell2at1[zone] = cgindex[cg];
7833 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7839 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7841 if (!bLocalCG[link->a[i]])
7850 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7852 real c[DIM][4]; /* the corners for the non-bonded communication */
7853 real cr0; /* corner for rounding */
7854 real cr1[4]; /* corners for rounding */
7855 real bc[DIM]; /* corners for bounded communication */
7856 real bcr1; /* corner for rounding for bonded communication */
7859 /* Determine the corners of the domain(s) we are communicating with */
7861 set_dd_corners(const gmx_domdec_t *dd,
7862 int dim0, int dim1, int dim2,
7866 const gmx_domdec_comm_t *comm;
7867 const gmx_domdec_zones_t *zones;
7872 zones = &comm->zones;
7874 /* Keep the compiler happy */
7878 /* The first dimension is equal for all cells */
7879 c->c[0][0] = comm->cell_x0[dim0];
7882 c->bc[0] = c->c[0][0];
7887 /* This cell row is only seen from the first row */
7888 c->c[1][0] = comm->cell_x0[dim1];
7889 /* All rows can see this row */
7890 c->c[1][1] = comm->cell_x0[dim1];
7893 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7896 /* For the multi-body distance we need the maximum */
7897 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7900 /* Set the upper-right corner for rounding */
7901 c->cr0 = comm->cell_x1[dim0];
7906 for (j = 0; j < 4; j++)
7908 c->c[2][j] = comm->cell_x0[dim2];
7912 /* Use the maximum of the i-cells that see a j-cell */
7913 for (i = 0; i < zones->nizone; i++)
7915 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7920 std::max(c->c[2][j-4],
7921 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7927 /* For the multi-body distance we need the maximum */
7928 c->bc[2] = comm->cell_x0[dim2];
7929 for (i = 0; i < 2; i++)
7931 for (j = 0; j < 2; j++)
7933 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7939 /* Set the upper-right corner for rounding */
7940 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7941 * Only cell (0,0,0) can see cell 7 (1,1,1)
7943 c->cr1[0] = comm->cell_x1[dim1];
7944 c->cr1[3] = comm->cell_x1[dim1];
7947 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7950 /* For the multi-body distance we need the maximum */
7951 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7958 /* Determine which cg's we need to send in this pulse from this zone */
7960 get_zone_pulse_cgs(gmx_domdec_t *dd,
7961 int zonei, int zone,
7963 const int *index_gl,
7965 int dim, int dim_ind,
7966 int dim0, int dim1, int dim2,
7967 real r_comm2, real r_bcomm2,
7971 real skew_fac2_d, real skew_fac_01,
7972 rvec *v_d, rvec *v_0, rvec *v_1,
7973 const dd_corners_t *c,
7975 gmx_bool bDistBonded,
7981 gmx_domdec_ind_t *ind,
7982 int **ibuf, int *ibuf_nalloc,
7988 gmx_domdec_comm_t *comm;
7990 gmx_bool bDistMB_pulse;
7992 real r2, rb2, r, tric_sh;
7995 int nsend_z, nsend, nat;
7999 bScrew = (dd->bScrewPBC && dim == XX);
8001 bDistMB_pulse = (bDistMB && bDistBonded);
8007 for (cg = cg0; cg < cg1; cg++)
8011 if (tric_dist[dim_ind] == 0)
8013 /* Rectangular direction, easy */
8014 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
8021 r = cg_cm[cg][dim] - c->bc[dim_ind];
8027 /* Rounding gives at most a 16% reduction
8028 * in communicated atoms
8030 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8032 r = cg_cm[cg][dim0] - c->cr0;
8033 /* This is the first dimension, so always r >= 0 */
8040 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8042 r = cg_cm[cg][dim1] - c->cr1[zone];
8049 r = cg_cm[cg][dim1] - c->bcr1;
8059 /* Triclinic direction, more complicated */
8062 /* Rounding, conservative as the skew_fac multiplication
8063 * will slightly underestimate the distance.
8065 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8067 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
8068 for (i = dim0+1; i < DIM; i++)
8070 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
8072 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
8075 rb[dim0] = rn[dim0];
8078 /* Take care that the cell planes along dim0 might not
8079 * be orthogonal to those along dim1 and dim2.
8081 for (i = 1; i <= dim_ind; i++)
8084 if (normal[dim0][dimd] > 0)
8086 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
8089 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
8094 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8096 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
8098 for (i = dim1+1; i < DIM; i++)
8100 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
8102 rn[dim1] += tric_sh;
8105 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
8106 /* Take care of coupling of the distances
8107 * to the planes along dim0 and dim1 through dim2.
8109 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
8110 /* Take care that the cell planes along dim1
8111 * might not be orthogonal to that along dim2.
8113 if (normal[dim1][dim2] > 0)
8115 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8121 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8124 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8125 /* Take care of coupling of the distances
8126 * to the planes along dim0 and dim1 through dim2.
8128 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8129 /* Take care that the cell planes along dim1
8130 * might not be orthogonal to that along dim2.
8132 if (normal[dim1][dim2] > 0)
8134 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8139 /* The distance along the communication direction */
8140 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8142 for (i = dim+1; i < DIM; i++)
8144 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8149 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8150 /* Take care of coupling of the distances
8151 * to the planes along dim0 and dim1 through dim2.
8153 if (dim_ind == 1 && zonei == 1)
8155 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8161 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8164 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8165 /* Take care of coupling of the distances
8166 * to the planes along dim0 and dim1 through dim2.
8168 if (dim_ind == 1 && zonei == 1)
8170 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8178 ((bDistMB && rb2 < r_bcomm2) ||
8179 (bDist2B && r2 < r_bcomm2)) &&
8181 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8182 missing_link(comm->cglink, index_gl[cg],
8185 /* Make an index to the local charge groups */
8186 if (nsend+1 > ind->nalloc)
8188 ind->nalloc = over_alloc_large(nsend+1);
8189 srenew(ind->index, ind->nalloc);
8191 if (nsend+1 > *ibuf_nalloc)
8193 *ibuf_nalloc = over_alloc_large(nsend+1);
8194 srenew(*ibuf, *ibuf_nalloc);
8196 ind->index[nsend] = cg;
8197 (*ibuf)[nsend] = index_gl[cg];
8199 vec_rvec_check_alloc(vbuf, nsend+1);
8201 if (dd->ci[dim] == 0)
8203 /* Correct cg_cm for pbc */
8204 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8207 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8208 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8213 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8216 nat += cgindex[cg+1] - cgindex[cg];
8222 *nsend_z_ptr = nsend_z;
8225 static void setup_dd_communication(gmx_domdec_t *dd,
8226 matrix box, gmx_ddbox_t *ddbox,
8227 t_forcerec *fr, t_state *state, rvec **f)
8229 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8230 int nzone, nzone_send, zone, zonei, cg0, cg1;
8231 int c, i, cg, cg_gl, nrcg;
8232 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8233 gmx_domdec_comm_t *comm;
8234 gmx_domdec_zones_t *zones;
8235 gmx_domdec_comm_dim_t *cd;
8236 gmx_domdec_ind_t *ind;
8237 cginfo_mb_t *cginfo_mb;
8238 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8239 real r_comm2, r_bcomm2;
8240 dd_corners_t corners;
8242 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8243 real skew_fac2_d, skew_fac_01;
8250 fprintf(debug, "Setting up DD communication\n");
8255 switch (fr->cutoff_scheme)
8264 gmx_incons("unimplemented");
8268 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8270 /* Check if we need to use triclinic distances */
8271 tric_dist[dim_ind] = 0;
8272 for (i = 0; i <= dim_ind; i++)
8274 if (ddbox->tric_dir[dd->dim[i]])
8276 tric_dist[dim_ind] = 1;
8281 bBondComm = comm->bBondComm;
8283 /* Do we need to determine extra distances for multi-body bondeds? */
8284 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8286 /* Do we need to determine extra distances for only two-body bondeds? */
8287 bDist2B = (bBondComm && !bDistMB);
8289 r_comm2 = sqr(comm->cutoff);
8290 r_bcomm2 = sqr(comm->cutoff_mbody);
8294 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8297 zones = &comm->zones;
8300 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8301 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8303 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8305 /* Triclinic stuff */
8306 normal = ddbox->normal;
8310 v_0 = ddbox->v[dim0];
8311 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8313 /* Determine the coupling coefficient for the distances
8314 * to the cell planes along dim0 and dim1 through dim2.
8315 * This is required for correct rounding.
8318 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8321 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8327 v_1 = ddbox->v[dim1];
8330 zone_cg_range = zones->cg_range;
8331 index_gl = dd->index_gl;
8332 cgindex = dd->cgindex;
8333 cginfo_mb = fr->cginfo_mb;
8335 zone_cg_range[0] = 0;
8336 zone_cg_range[1] = dd->ncg_home;
8337 comm->zone_ncg1[0] = dd->ncg_home;
8338 pos_cg = dd->ncg_home;
8340 nat_tot = dd->nat_home;
8342 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8344 dim = dd->dim[dim_ind];
8345 cd = &comm->cd[dim_ind];
8347 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8349 /* No pbc in this dimension, the first node should not comm. */
8357 v_d = ddbox->v[dim];
8358 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8360 cd->bInPlace = TRUE;
8361 for (p = 0; p < cd->np; p++)
8363 /* Only atoms communicated in the first pulse are used
8364 * for multi-body bonded interactions or for bBondComm.
8366 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8371 for (zone = 0; zone < nzone_send; zone++)
8373 if (tric_dist[dim_ind] && dim_ind > 0)
8375 /* Determine slightly more optimized skew_fac's
8377 * This reduces the number of communicated atoms
8378 * by about 10% for 3D DD of rhombic dodecahedra.
8380 for (dimd = 0; dimd < dim; dimd++)
8382 sf2_round[dimd] = 1;
8383 if (ddbox->tric_dir[dimd])
8385 for (i = dd->dim[dimd]+1; i < DIM; i++)
8387 /* If we are shifted in dimension i
8388 * and the cell plane is tilted forward
8389 * in dimension i, skip this coupling.
8391 if (!(zones->shift[nzone+zone][i] &&
8392 ddbox->v[dimd][i][dimd] >= 0))
8395 sqr(ddbox->v[dimd][i][dimd]);
8398 sf2_round[dimd] = 1/sf2_round[dimd];
8403 zonei = zone_perm[dim_ind][zone];
8406 /* Here we permutate the zones to obtain a convenient order
8407 * for neighbor searching
8409 cg0 = zone_cg_range[zonei];
8410 cg1 = zone_cg_range[zonei+1];
8414 /* Look only at the cg's received in the previous grid pulse
8416 cg1 = zone_cg_range[nzone+zone+1];
8417 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8420 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8421 for (th = 0; th < comm->nth; th++)
8423 gmx_domdec_ind_t *ind_p;
8424 int **ibuf_p, *ibuf_nalloc_p;
8426 int *nsend_p, *nat_p;
8432 /* Thread 0 writes in the comm buffers */
8434 ibuf_p = &comm->buf_int;
8435 ibuf_nalloc_p = &comm->nalloc_int;
8436 vbuf_p = &comm->vbuf;
8439 nsend_zone_p = &ind->nsend[zone];
8443 /* Other threads write into temp buffers */
8444 ind_p = &comm->dth[th].ind;
8445 ibuf_p = &comm->dth[th].ibuf;
8446 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8447 vbuf_p = &comm->dth[th].vbuf;
8448 nsend_p = &comm->dth[th].nsend;
8449 nat_p = &comm->dth[th].nat;
8450 nsend_zone_p = &comm->dth[th].nsend_zone;
8452 comm->dth[th].nsend = 0;
8453 comm->dth[th].nat = 0;
8454 comm->dth[th].nsend_zone = 0;
8464 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8465 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8468 /* Get the cg's for this pulse in this zone */
8469 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8471 dim, dim_ind, dim0, dim1, dim2,
8474 normal, skew_fac2_d, skew_fac_01,
8475 v_d, v_0, v_1, &corners, sf2_round,
8476 bDistBonded, bBondComm,
8480 ibuf_p, ibuf_nalloc_p,
8486 /* Append data of threads>=1 to the communication buffers */
8487 for (th = 1; th < comm->nth; th++)
8489 dd_comm_setup_work_t *dth;
8492 dth = &comm->dth[th];
8494 ns1 = nsend + dth->nsend_zone;
8495 if (ns1 > ind->nalloc)
8497 ind->nalloc = over_alloc_dd(ns1);
8498 srenew(ind->index, ind->nalloc);
8500 if (ns1 > comm->nalloc_int)
8502 comm->nalloc_int = over_alloc_dd(ns1);
8503 srenew(comm->buf_int, comm->nalloc_int);
8505 if (ns1 > comm->vbuf.nalloc)
8507 comm->vbuf.nalloc = over_alloc_dd(ns1);
8508 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8511 for (i = 0; i < dth->nsend_zone; i++)
8513 ind->index[nsend] = dth->ind.index[i];
8514 comm->buf_int[nsend] = dth->ibuf[i];
8515 copy_rvec(dth->vbuf.v[i],
8516 comm->vbuf.v[nsend]);
8520 ind->nsend[zone] += dth->nsend_zone;
8523 /* Clear the counts in case we do not have pbc */
8524 for (zone = nzone_send; zone < nzone; zone++)
8526 ind->nsend[zone] = 0;
8528 ind->nsend[nzone] = nsend;
8529 ind->nsend[nzone+1] = nat;
8530 /* Communicate the number of cg's and atoms to receive */
8531 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8532 ind->nsend, nzone+2,
8533 ind->nrecv, nzone+2);
8535 /* The rvec buffer is also required for atom buffers of size nsend
8536 * in dd_move_x and dd_move_f.
8538 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8542 /* We can receive in place if only the last zone is not empty */
8543 for (zone = 0; zone < nzone-1; zone++)
8545 if (ind->nrecv[zone] > 0)
8547 cd->bInPlace = FALSE;
8552 /* The int buffer is only required here for the cg indices */
8553 if (ind->nrecv[nzone] > comm->nalloc_int2)
8555 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8556 srenew(comm->buf_int2, comm->nalloc_int2);
8558 /* The rvec buffer is also required for atom buffers
8559 * of size nrecv in dd_move_x and dd_move_f.
8561 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8562 vec_rvec_check_alloc(&comm->vbuf2, i);
8566 /* Make space for the global cg indices */
8567 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8568 || dd->cg_nalloc == 0)
8570 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8571 srenew(index_gl, dd->cg_nalloc);
8572 srenew(cgindex, dd->cg_nalloc+1);
8574 /* Communicate the global cg indices */
8577 recv_i = index_gl + pos_cg;
8581 recv_i = comm->buf_int2;
8583 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8584 comm->buf_int, nsend,
8585 recv_i, ind->nrecv[nzone]);
8587 /* Make space for cg_cm */
8588 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8589 if (fr->cutoff_scheme == ecutsGROUP)
8597 /* Communicate cg_cm */
8600 recv_vr = cg_cm + pos_cg;
8604 recv_vr = comm->vbuf2.v;
8606 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8607 comm->vbuf.v, nsend,
8608 recv_vr, ind->nrecv[nzone]);
8610 /* Make the charge group index */
8613 zone = (p == 0 ? 0 : nzone - 1);
8614 while (zone < nzone)
8616 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8618 cg_gl = index_gl[pos_cg];
8619 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8620 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8621 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8624 /* Update the charge group presence,
8625 * so we can use it in the next pass of the loop.
8627 comm->bLocalCG[cg_gl] = TRUE;
8633 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8636 zone_cg_range[nzone+zone] = pos_cg;
8641 /* This part of the code is never executed with bBondComm. */
8642 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8643 index_gl, recv_i, cg_cm, recv_vr,
8644 cgindex, fr->cginfo_mb, fr->cginfo);
8645 pos_cg += ind->nrecv[nzone];
8647 nat_tot += ind->nrecv[nzone+1];
8651 /* Store the atom block for easy copying of communication buffers */
8652 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8656 dd->index_gl = index_gl;
8657 dd->cgindex = cgindex;
8659 dd->ncg_tot = zone_cg_range[zones->n];
8660 dd->nat_tot = nat_tot;
8661 comm->nat[ddnatHOME] = dd->nat_home;
8662 for (i = ddnatZONE; i < ddnatNR; i++)
8664 comm->nat[i] = dd->nat_tot;
8669 /* We don't need to update cginfo, since that was alrady done above.
8670 * So we pass NULL for the forcerec.
8672 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8673 NULL, comm->bLocalCG);
8678 fprintf(debug, "Finished setting up DD communication, zones:");
8679 for (c = 0; c < zones->n; c++)
8681 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8683 fprintf(debug, "\n");
8687 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8691 for (c = 0; c < zones->nizone; c++)
8693 zones->izone[c].cg1 = zones->cg_range[c+1];
8694 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8695 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8699 static void set_zones_size(gmx_domdec_t *dd,
8700 matrix box, const gmx_ddbox_t *ddbox,
8701 int zone_start, int zone_end)
8703 gmx_domdec_comm_t *comm;
8704 gmx_domdec_zones_t *zones;
8713 zones = &comm->zones;
8715 /* Do we need to determine extra distances for multi-body bondeds? */
8716 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8718 for (z = zone_start; z < zone_end; z++)
8720 /* Copy cell limits to zone limits.
8721 * Valid for non-DD dims and non-shifted dims.
8723 copy_rvec(comm->cell_x0, zones->size[z].x0);
8724 copy_rvec(comm->cell_x1, zones->size[z].x1);
8727 for (d = 0; d < dd->ndim; d++)
8731 for (z = 0; z < zones->n; z++)
8733 /* With a staggered grid we have different sizes
8734 * for non-shifted dimensions.
8736 if (dd->bGridJump && zones->shift[z][dim] == 0)
8740 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8741 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8745 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8746 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8752 rcmbs = comm->cutoff_mbody;
8753 if (ddbox->tric_dir[dim])
8755 rcs /= ddbox->skew_fac[dim];
8756 rcmbs /= ddbox->skew_fac[dim];
8759 /* Set the lower limit for the shifted zone dimensions */
8760 for (z = zone_start; z < zone_end; z++)
8762 if (zones->shift[z][dim] > 0)
8765 if (!dd->bGridJump || d == 0)
8767 zones->size[z].x0[dim] = comm->cell_x1[dim];
8768 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8772 /* Here we take the lower limit of the zone from
8773 * the lowest domain of the zone below.
8777 zones->size[z].x0[dim] =
8778 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8784 zones->size[z].x0[dim] =
8785 zones->size[zone_perm[2][z-4]].x0[dim];
8789 zones->size[z].x0[dim] =
8790 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8793 /* A temporary limit, is updated below */
8794 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8798 for (zi = 0; zi < zones->nizone; zi++)
8800 if (zones->shift[zi][dim] == 0)
8802 /* This takes the whole zone into account.
8803 * With multiple pulses this will lead
8804 * to a larger zone then strictly necessary.
8806 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8807 zones->size[zi].x1[dim]+rcmbs);
8815 /* Loop over the i-zones to set the upper limit of each
8818 for (zi = 0; zi < zones->nizone; zi++)
8820 if (zones->shift[zi][dim] == 0)
8822 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8824 if (zones->shift[z][dim] > 0)
8826 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8827 zones->size[zi].x1[dim]+rcs);
8834 for (z = zone_start; z < zone_end; z++)
8836 /* Initialization only required to keep the compiler happy */
8837 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8840 /* To determine the bounding box for a zone we need to find
8841 * the extreme corners of 4, 2 or 1 corners.
8843 nc = 1 << (ddbox->nboundeddim - 1);
8845 for (c = 0; c < nc; c++)
8847 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8851 corner[YY] = zones->size[z].x0[YY];
8855 corner[YY] = zones->size[z].x1[YY];
8859 corner[ZZ] = zones->size[z].x0[ZZ];
8863 corner[ZZ] = zones->size[z].x1[ZZ];
8865 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8866 box[ZZ][1 - dd->dim[0]] != 0)
8868 /* With 1D domain decomposition the cg's are not in
8869 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8870 * Shift the corner of the z-vector back to along the box
8871 * vector of dimension d, so it will later end up at 0 along d.
8872 * This can affect the location of this corner along dd->dim[0]
8873 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8875 int d = 1 - dd->dim[0];
8877 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8879 /* Apply the triclinic couplings */
8880 assert(ddbox->npbcdim <= DIM);
8881 for (i = YY; i < ddbox->npbcdim; i++)
8883 for (j = XX; j < i; j++)
8885 corner[j] += corner[i]*box[i][j]/box[i][i];
8890 copy_rvec(corner, corner_min);
8891 copy_rvec(corner, corner_max);
8895 for (i = 0; i < DIM; i++)
8897 corner_min[i] = std::min(corner_min[i], corner[i]);
8898 corner_max[i] = std::max(corner_max[i], corner[i]);
8902 /* Copy the extreme cornes without offset along x */
8903 for (i = 0; i < DIM; i++)
8905 zones->size[z].bb_x0[i] = corner_min[i];
8906 zones->size[z].bb_x1[i] = corner_max[i];
8908 /* Add the offset along x */
8909 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8910 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8913 if (zone_start == 0)
8916 for (dim = 0; dim < DIM; dim++)
8918 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8920 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8925 for (z = zone_start; z < zone_end; z++)
8927 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8929 zones->size[z].x0[XX], zones->size[z].x1[XX],
8930 zones->size[z].x0[YY], zones->size[z].x1[YY],
8931 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8932 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8934 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8935 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8936 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8941 static int comp_cgsort(const void *a, const void *b)
8945 gmx_cgsort_t *cga, *cgb;
8946 cga = (gmx_cgsort_t *)a;
8947 cgb = (gmx_cgsort_t *)b;
8949 comp = cga->nsc - cgb->nsc;
8952 comp = cga->ind_gl - cgb->ind_gl;
8958 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8963 /* Order the data */
8964 for (i = 0; i < n; i++)
8966 buf[i] = a[sort[i].ind];
8969 /* Copy back to the original array */
8970 for (i = 0; i < n; i++)
8976 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8981 /* Order the data */
8982 for (i = 0; i < n; i++)
8984 copy_rvec(v[sort[i].ind], buf[i]);
8987 /* Copy back to the original array */
8988 for (i = 0; i < n; i++)
8990 copy_rvec(buf[i], v[i]);
8994 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8997 int a, atot, cg, cg0, cg1, i;
8999 if (cgindex == NULL)
9001 /* Avoid the useless loop of the atoms within a cg */
9002 order_vec_cg(ncg, sort, v, buf);
9007 /* Order the data */
9009 for (cg = 0; cg < ncg; cg++)
9011 cg0 = cgindex[sort[cg].ind];
9012 cg1 = cgindex[sort[cg].ind+1];
9013 for (i = cg0; i < cg1; i++)
9015 copy_rvec(v[i], buf[a]);
9021 /* Copy back to the original array */
9022 for (a = 0; a < atot; a++)
9024 copy_rvec(buf[a], v[a]);
9028 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
9029 int nsort_new, gmx_cgsort_t *sort_new,
9030 gmx_cgsort_t *sort1)
9034 /* The new indices are not very ordered, so we qsort them */
9035 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
9037 /* sort2 is already ordered, so now we can merge the two arrays */
9041 while (i2 < nsort2 || i_new < nsort_new)
9045 sort1[i1++] = sort_new[i_new++];
9047 else if (i_new == nsort_new)
9049 sort1[i1++] = sort2[i2++];
9051 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
9052 (sort2[i2].nsc == sort_new[i_new].nsc &&
9053 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
9055 sort1[i1++] = sort2[i2++];
9059 sort1[i1++] = sort_new[i_new++];
9064 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
9066 gmx_domdec_sort_t *sort;
9067 gmx_cgsort_t *cgsort, *sort_i;
9068 int ncg_new, nsort2, nsort_new, i, *a, moved;
9070 sort = dd->comm->sort;
9072 a = fr->ns.grid->cell_index;
9074 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
9076 if (ncg_home_old >= 0)
9078 /* The charge groups that remained in the same ns grid cell
9079 * are completely ordered. So we can sort efficiently by sorting
9080 * the charge groups that did move into the stationary list.
9085 for (i = 0; i < dd->ncg_home; i++)
9087 /* Check if this cg did not move to another node */
9090 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
9092 /* This cg is new on this node or moved ns grid cell */
9093 if (nsort_new >= sort->sort_new_nalloc)
9095 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
9096 srenew(sort->sort_new, sort->sort_new_nalloc);
9098 sort_i = &(sort->sort_new[nsort_new++]);
9102 /* This cg did not move */
9103 sort_i = &(sort->sort2[nsort2++]);
9105 /* Sort on the ns grid cell indices
9106 * and the global topology index.
9107 * index_gl is irrelevant with cell ns,
9108 * but we set it here anyhow to avoid a conditional.
9111 sort_i->ind_gl = dd->index_gl[i];
9118 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9121 /* Sort efficiently */
9122 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9127 cgsort = sort->sort;
9129 for (i = 0; i < dd->ncg_home; i++)
9131 /* Sort on the ns grid cell indices
9132 * and the global topology index
9134 cgsort[i].nsc = a[i];
9135 cgsort[i].ind_gl = dd->index_gl[i];
9137 if (cgsort[i].nsc < moved)
9144 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9146 /* Determine the order of the charge groups using qsort */
9147 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9153 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9156 int ncg_new, i, *a, na;
9158 sort = dd->comm->sort->sort;
9160 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9163 for (i = 0; i < na; i++)
9167 sort[ncg_new].ind = a[i];
9175 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9178 gmx_domdec_sort_t *sort;
9179 gmx_cgsort_t *cgsort;
9181 int ncg_new, i, *ibuf, cgsize;
9184 sort = dd->comm->sort;
9186 if (dd->ncg_home > sort->sort_nalloc)
9188 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9189 srenew(sort->sort, sort->sort_nalloc);
9190 srenew(sort->sort2, sort->sort_nalloc);
9192 cgsort = sort->sort;
9194 switch (fr->cutoff_scheme)
9197 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9200 ncg_new = dd_sort_order_nbnxn(dd, fr);
9203 gmx_incons("unimplemented");
9207 /* We alloc with the old size, since cgindex is still old */
9208 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9209 vbuf = dd->comm->vbuf.v;
9213 cgindex = dd->cgindex;
9220 /* Remove the charge groups which are no longer at home here */
9221 dd->ncg_home = ncg_new;
9224 fprintf(debug, "Set the new home charge group count to %d\n",
9228 /* Reorder the state */
9229 for (i = 0; i < estNR; i++)
9231 if (EST_DISTR(i) && (state->flags & (1<<i)))
9236 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9239 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9242 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9245 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9249 case estDISRE_INITF:
9250 case estDISRE_RM3TAV:
9251 case estORIRE_INITF:
9253 /* No ordering required */
9256 gmx_incons("Unknown state entry encountered in dd_sort_state");
9261 if (fr->cutoff_scheme == ecutsGROUP)
9264 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9267 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9269 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9270 srenew(sort->ibuf, sort->ibuf_nalloc);
9273 /* Reorder the global cg index */
9274 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9275 /* Reorder the cginfo */
9276 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9277 /* Rebuild the local cg index */
9281 for (i = 0; i < dd->ncg_home; i++)
9283 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9284 ibuf[i+1] = ibuf[i] + cgsize;
9286 for (i = 0; i < dd->ncg_home+1; i++)
9288 dd->cgindex[i] = ibuf[i];
9293 for (i = 0; i < dd->ncg_home+1; i++)
9298 /* Set the home atom number */
9299 dd->nat_home = dd->cgindex[dd->ncg_home];
9301 if (fr->cutoff_scheme == ecutsVERLET)
9303 /* The atoms are now exactly in grid order, update the grid order */
9304 nbnxn_set_atomorder(fr->nbv->nbs);
9308 /* Copy the sorted ns cell indices back to the ns grid struct */
9309 for (i = 0; i < dd->ncg_home; i++)
9311 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9313 fr->ns.grid->nr = dd->ncg_home;
9317 static void add_dd_statistics(gmx_domdec_t *dd)
9319 gmx_domdec_comm_t *comm;
9324 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9326 comm->sum_nat[ddnat-ddnatZONE] +=
9327 comm->nat[ddnat] - comm->nat[ddnat-1];
9332 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9334 gmx_domdec_comm_t *comm;
9339 /* Reset all the statistics and counters for total run counting */
9340 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9342 comm->sum_nat[ddnat-ddnatZONE] = 0;
9346 comm->load_step = 0;
9349 clear_ivec(comm->load_lim);
9354 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9356 gmx_domdec_comm_t *comm;
9360 comm = cr->dd->comm;
9362 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9369 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");
9371 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9373 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9378 " av. #atoms communicated per step for force: %d x %.1f\n",
9382 if (cr->dd->vsite_comm)
9385 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9386 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9391 if (cr->dd->constraint_comm)
9394 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9395 1 + ir->nLincsIter, av);
9399 gmx_incons(" Unknown type for DD statistics");
9402 fprintf(fplog, "\n");
9404 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9406 print_dd_load_av(fplog, cr->dd);
9410 void dd_partition_system(FILE *fplog,
9413 gmx_bool bMasterState,
9415 t_state *state_global,
9416 gmx_mtop_t *top_global,
9418 t_state *state_local,
9421 gmx_localtop_t *top_local,
9424 gmx_shellfc_t shellfc,
9425 gmx_constr_t constr,
9427 gmx_wallcycle_t wcycle,
9431 gmx_domdec_comm_t *comm;
9432 gmx_ddbox_t ddbox = {0};
9434 gmx_int64_t step_pcoupl;
9435 rvec cell_ns_x0, cell_ns_x1;
9436 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9437 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bTurnOnDLB, bLogLoad;
9438 gmx_bool bRedist, bSortCG, bResortAll;
9439 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9443 wallcycle_start(wcycle, ewcDOMDEC);
9448 bBoxChanged = (bMasterState || DEFORM(*ir));
9449 if (ir->epc != epcNO)
9451 /* With nstpcouple > 1 pressure coupling happens.
9452 * one step after calculating the pressure.
9453 * Box scaling happens at the end of the MD step,
9454 * after the DD partitioning.
9455 * We therefore have to do DLB in the first partitioning
9456 * after an MD step where P-coupling occured.
9457 * We need to determine the last step in which p-coupling occurred.
9458 * MRS -- need to validate this for vv?
9463 step_pcoupl = step - 1;
9467 step_pcoupl = ((step - 1)/n)*n + 1;
9469 if (step_pcoupl >= comm->partition_step)
9475 bNStGlobalComm = (step % nstglobalcomm == 0);
9477 if (!comm->bDynLoadBal)
9483 /* Should we do dynamic load balacing this step?
9484 * Since it requires (possibly expensive) global communication,
9485 * we might want to do DLB less frequently.
9487 if (bBoxChanged || ir->epc != epcNO)
9489 bDoDLB = bBoxChanged;
9493 bDoDLB = bNStGlobalComm;
9497 /* Check if we have recorded loads on the nodes */
9498 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9500 bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
9502 /* Print load every nstlog, first and last step to the log file */
9503 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9504 comm->n_load_collect == 0 ||
9506 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9508 /* Avoid extra communication due to verbose screen output
9509 * when nstglobalcomm is set.
9511 if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
9512 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9514 get_load_distribution(dd, wcycle);
9519 dd_print_load(fplog, dd, step-1);
9523 dd_print_load_verbose(dd);
9526 comm->n_load_collect++;
9528 if (bCheckWhetherToTurnDlbOn)
9530 /* Since the timings are node dependent, the master decides */
9533 /* Here we check if the max PME rank load is more than 0.98
9534 * the max PP force load. If so, PP DLB will not help,
9535 * since we are (almost) limited by PME. Furthermore,
9536 * DLB will cause a significant extra x/f redistribution
9537 * cost on the PME ranks, which will then surely result
9538 * in lower total performance.
9539 * This check might be fragile, since one measurement
9540 * below 0.98 (although only done once every 100 DD part.)
9541 * could turn on DLB for the rest of the run.
9543 if (cr->npmenodes > 0 &&
9544 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9551 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9555 fprintf(debug, "step %s, imb loss %f\n",
9556 gmx_step_str(step, sbuf),
9557 dd_force_imb_perf_loss(dd));
9560 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9563 turn_on_dlb(fplog, cr, step);
9568 comm->n_load_have++;
9571 cgs_gl = &comm->cgs_gl;
9576 /* Clear the old state */
9577 clear_dd_indices(dd, 0, 0);
9580 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9581 TRUE, cgs_gl, state_global->x, &ddbox);
9583 get_cg_distribution(fplog, dd, cgs_gl,
9584 state_global->box, &ddbox, state_global->x);
9586 dd_distribute_state(dd, cgs_gl,
9587 state_global, state_local, f);
9589 dd_make_local_cgs(dd, &top_local->cgs);
9591 /* Ensure that we have space for the new distribution */
9592 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9594 if (fr->cutoff_scheme == ecutsGROUP)
9596 calc_cgcm(fplog, 0, dd->ncg_home,
9597 &top_local->cgs, state_local->x, fr->cg_cm);
9600 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9602 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9604 else if (state_local->ddp_count != dd->ddp_count)
9606 if (state_local->ddp_count > dd->ddp_count)
9608 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9611 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9613 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);
9616 /* Clear the old state */
9617 clear_dd_indices(dd, 0, 0);
9619 /* Build the new indices */
9620 rebuild_cgindex(dd, cgs_gl->index, state_local);
9621 make_dd_indices(dd, cgs_gl->index, 0);
9622 ncgindex_set = dd->ncg_home;
9624 if (fr->cutoff_scheme == ecutsGROUP)
9626 /* Redetermine the cg COMs */
9627 calc_cgcm(fplog, 0, dd->ncg_home,
9628 &top_local->cgs, state_local->x, fr->cg_cm);
9631 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9633 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9635 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9636 TRUE, &top_local->cgs, state_local->x, &ddbox);
9638 bRedist = comm->bDynLoadBal;
9642 /* We have the full state, only redistribute the cgs */
9644 /* Clear the non-home indices */
9645 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9648 /* Avoid global communication for dim's without pbc and -gcom */
9649 if (!bNStGlobalComm)
9651 copy_rvec(comm->box0, ddbox.box0 );
9652 copy_rvec(comm->box_size, ddbox.box_size);
9654 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9655 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9660 /* For dim's without pbc and -gcom */
9661 copy_rvec(ddbox.box0, comm->box0 );
9662 copy_rvec(ddbox.box_size, comm->box_size);
9664 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9667 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9669 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9672 /* Check if we should sort the charge groups */
9673 if (comm->nstSortCG > 0)
9675 bSortCG = (bMasterState ||
9676 (bRedist && (step % comm->nstSortCG == 0)));
9683 ncg_home_old = dd->ncg_home;
9688 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9690 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9692 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9694 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9697 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9699 &comm->cell_x0, &comm->cell_x1,
9700 dd->ncg_home, fr->cg_cm,
9701 cell_ns_x0, cell_ns_x1, &grid_density);
9705 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9708 switch (fr->cutoff_scheme)
9711 copy_ivec(fr->ns.grid->n, ncells_old);
9712 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9713 state_local->box, cell_ns_x0, cell_ns_x1,
9714 fr->rlistlong, grid_density);
9717 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9720 gmx_incons("unimplemented");
9722 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9723 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9727 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9729 /* Sort the state on charge group position.
9730 * This enables exact restarts from this step.
9731 * It also improves performance by about 15% with larger numbers
9732 * of atoms per node.
9735 /* Fill the ns grid with the home cell,
9736 * so we can sort with the indices.
9738 set_zones_ncg_home(dd);
9740 switch (fr->cutoff_scheme)
9743 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9745 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9747 comm->zones.size[0].bb_x0,
9748 comm->zones.size[0].bb_x1,
9750 comm->zones.dens_zone0,
9753 ncg_moved, bRedist ? comm->moved : NULL,
9754 fr->nbv->grp[eintLocal].kernel_type,
9755 fr->nbv->grp[eintLocal].nbat);
9757 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9760 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9761 0, dd->ncg_home, fr->cg_cm);
9763 copy_ivec(fr->ns.grid->n, ncells_new);
9766 gmx_incons("unimplemented");
9769 bResortAll = bMasterState;
9771 /* Check if we can user the old order and ns grid cell indices
9772 * of the charge groups to sort the charge groups efficiently.
9774 if (ncells_new[XX] != ncells_old[XX] ||
9775 ncells_new[YY] != ncells_old[YY] ||
9776 ncells_new[ZZ] != ncells_old[ZZ])
9783 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9784 gmx_step_str(step, sbuf), dd->ncg_home);
9786 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9787 bResortAll ? -1 : ncg_home_old);
9788 /* Rebuild all the indices */
9789 ga2la_clear(dd->ga2la);
9792 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9795 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9797 /* Setup up the communication and communicate the coordinates */
9798 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9800 /* Set the indices */
9801 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9803 /* Set the charge group boundaries for neighbor searching */
9804 set_cg_boundaries(&comm->zones);
9806 if (fr->cutoff_scheme == ecutsVERLET)
9808 set_zones_size(dd, state_local->box, &ddbox,
9809 bSortCG ? 1 : 0, comm->zones.n);
9812 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9815 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9816 -1,state_local->x,state_local->box);
9819 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9821 /* Extract a local topology from the global topology */
9822 for (i = 0; i < dd->ndim; i++)
9824 np[dd->dim[i]] = comm->cd[i].np;
9826 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9827 comm->cellsize_min, np,
9829 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9830 vsite, top_global, top_local);
9832 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9834 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9836 /* Set up the special atom communication */
9837 n = comm->nat[ddnatZONE];
9838 for (i = ddnatZONE+1; i < ddnatNR; i++)
9843 if (vsite && vsite->n_intercg_vsite)
9845 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9849 if (dd->bInterCGcons || dd->bInterCGsettles)
9851 /* Only for inter-cg constraints we need special code */
9852 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9853 constr, ir->nProjOrder,
9854 top_local->idef.il);
9858 gmx_incons("Unknown special atom type setup");
9863 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9865 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9867 /* Make space for the extra coordinates for virtual site
9868 * or constraint communication.
9870 state_local->natoms = comm->nat[ddnatNR-1];
9871 if (state_local->natoms > state_local->nalloc)
9873 dd_realloc_state(state_local, f, state_local->natoms);
9876 if (fr->bF_NoVirSum)
9878 if (vsite && vsite->n_intercg_vsite)
9880 nat_f_novirsum = comm->nat[ddnatVSITE];
9884 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9886 nat_f_novirsum = dd->nat_tot;
9890 nat_f_novirsum = dd->nat_home;
9899 /* Set the number of atoms required for the force calculation.
9900 * Forces need to be constrained when using a twin-range setup
9901 * or with energy minimization. For simple simulations we could
9902 * avoid some allocation, zeroing and copying, but this is
9903 * probably not worth the complications ande checking.
9905 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9906 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9908 /* We make the all mdatoms up to nat_tot_con.
9909 * We could save some work by only setting invmass
9910 * between nat_tot and nat_tot_con.
9912 /* This call also sets the new number of home particles to dd->nat_home */
9913 atoms2md(top_global, ir,
9914 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9916 /* Now we have the charges we can sort the FE interactions */
9917 dd_sort_local_top(dd, mdatoms, top_local);
9921 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9922 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9923 mdatoms, FALSE, vsite);
9928 /* Make the local shell stuff, currently no communication is done */
9929 make_local_shells(cr, mdatoms, shellfc);
9932 if (ir->implicit_solvent)
9934 make_local_gb(cr, fr->born, ir->gb_algorithm);
9937 setup_bonded_threading(fr, &top_local->idef);
9939 if (!(cr->duty & DUTY_PME))
9941 /* Send the charges and/or c6/sigmas to our PME only node */
9942 gmx_pme_send_parameters(cr,
9944 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9945 mdatoms->chargeA, mdatoms->chargeB,
9946 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9947 mdatoms->sigmaA, mdatoms->sigmaB,
9948 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9953 set_constraints(constr, top_local, ir, mdatoms, cr);
9958 /* Update the local pull groups */
9959 dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
9964 /* Update the local rotation groups */
9965 dd_make_local_rotation_groups(dd, ir->rot);
9968 if (ir->eSwapCoords != eswapNO)
9970 /* Update the local groups needed for ion swapping */
9971 dd_make_local_swap_groups(dd, ir->swap);
9974 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9975 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9977 add_dd_statistics(dd);
9979 /* Make sure we only count the cycles for this DD partitioning */
9980 clear_dd_cycle_counts(dd);
9982 /* Because the order of the atoms might have changed since
9983 * the last vsite construction, we need to communicate the constructing
9984 * atom coordinates again (for spreading the forces this MD step).
9986 dd_move_x_vsites(dd, state_local->box, state_local->x);
9988 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9990 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9992 dd_move_x(dd, state_local->box, state_local->x);
9993 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9994 -1, state_local->x, state_local->box);
9997 /* Store the partitioning step */
9998 comm->partition_step = step;
10000 /* Increase the DD partitioning counter */
10002 /* The state currently matches this DD partitioning count, store it */
10003 state_local->ddp_count = dd->ddp_count;
10006 /* The DD master node knows the complete cg distribution,
10007 * store the count so we can possibly skip the cg info communication.
10009 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
10012 if (comm->DD_debug > 0)
10014 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
10015 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
10016 "after partitioning");
10019 wallcycle_stop(wcycle, ewcDOMDEC);