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 edlbsOffForever, /* DLB is off and will never be turned on */
220 edlbsOffCanTurnOn, /* DLB is off and will turn on with imbalance */
221 edlbsOffTemporarilyLocked, /* DLB is off and temporarily can not turn on */
222 edlbsOn, /* DLB is on and will stay on forever */
225 /* Allowed DLB state transitions:
226 * edlbsOffCanTurnOn -> edlbsOn
227 * edlbsOffCanTurnOn -> edlbsOffForever
228 * edlbsOffCanTurnOn -> edlbsOffTemporarilyLocked
229 * edlbsOffTemporarilyLocked -> edlbsOffCanTurnOn
232 const char *edlbs_names[edlbsNR] = { "off", "auto", "locked", "on" };
236 int dim; /* The dimension */
237 gmx_bool dim_match; /* Tells if DD and PME dims match */
238 int nslab; /* The number of PME slabs in this dimension */
239 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
240 int *pp_min; /* The minimum pp node location, size nslab */
241 int *pp_max; /* The maximum pp node location,size nslab */
242 int maxshift; /* The maximum shift for coordinate redistribution in PME */
247 real min0; /* The minimum bottom of this zone */
248 real max1; /* The maximum top of this zone */
249 real min1; /* The minimum top of this zone */
250 real mch0; /* The maximum bottom communicaton height for this zone */
251 real mch1; /* The maximum top communicaton height for this zone */
252 real p1_0; /* The bottom value of the first cell in this zone */
253 real p1_1; /* The top value of the first cell in this zone */
258 gmx_domdec_ind_t ind;
265 } dd_comm_setup_work_t;
267 typedef struct gmx_domdec_comm
269 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
270 * unless stated otherwise.
273 /* The number of decomposition dimensions for PME, 0: no PME */
275 /* The number of nodes doing PME (PP/PME or only PME) */
279 /* The communication setup including the PME only nodes */
280 gmx_bool bCartesianPP_PME;
283 int *pmenodes; /* size npmenodes */
284 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
285 * but with bCartesianPP_PME */
286 gmx_ddpme_t ddpme[2];
288 /* The DD particle-particle nodes only */
289 gmx_bool bCartesianPP;
290 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
292 /* The global charge groups */
295 /* Should we sort the cgs */
297 gmx_domdec_sort_t *sort;
299 /* Are there charge groups? */
302 /* Are there bonded and multi-body interactions between charge groups? */
303 gmx_bool bInterCGBondeds;
304 gmx_bool bInterCGMultiBody;
306 /* Data for the optional bonded interaction atom communication range */
311 /* The DLB state, possible values are defined above */
313 /* With dlbState=edlbsOffCanTurnOn, should we check if to DLB on at the next DD? */
314 gmx_bool bCheckWhetherToTurnDlbOn;
316 /* Cell sizes for static load balancing, first index cartesian */
319 /* The width of the communicated boundaries */
322 /* The minimum cell size (including triclinic correction) */
324 /* For dlb, for use with edlbAUTO */
325 rvec cellsize_min_dlb;
326 /* The lower limit for the DD cell size with DLB */
328 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
329 gmx_bool bVacDLBNoLimit;
331 /* With PME load balancing we set limits on DLB */
332 gmx_bool bPMELoadBalDLBLimits;
333 /* DLB needs to take into account that we want to allow this maximum
334 * cut-off (for PME load balancing), this could limit cell boundaries.
336 real PMELoadBal_max_cutoff;
338 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
340 /* box0 and box_size are required with dim's without pbc and -gcom */
344 /* The cell boundaries */
348 /* The old location of the cell boundaries, to check cg displacements */
352 /* The communication setup and charge group boundaries for the zones */
353 gmx_domdec_zones_t zones;
355 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
356 * cell boundaries of neighboring cells for dynamic load balancing.
358 gmx_ddzone_t zone_d1[2];
359 gmx_ddzone_t zone_d2[2][2];
361 /* The coordinate/force communication setup and indices */
362 gmx_domdec_comm_dim_t cd[DIM];
363 /* The maximum number of cells to communicate with in one dimension */
366 /* Which cg distribution is stored on the master node */
367 int master_cg_ddp_count;
369 /* The number of cg's received from the direct neighbors */
370 int zone_ncg1[DD_MAXZONE];
372 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
375 /* Array for signalling if atoms have moved to another domain */
379 /* Communication buffer for general use */
383 /* Communication buffer for general use */
386 /* Temporary storage for thread parallel communication setup */
388 dd_comm_setup_work_t *dth;
390 /* Communication buffers only used with multiple grid pulses */
395 /* Communication buffers for local redistribution */
397 int cggl_flag_nalloc[DIM*2];
399 int cgcm_state_nalloc[DIM*2];
401 /* Cell sizes for dynamic load balancing */
402 gmx_domdec_root_t **root;
406 real cell_f_max0[DIM];
407 real cell_f_min1[DIM];
409 /* Stuff for load communication */
410 gmx_bool bRecordLoad;
411 gmx_domdec_load_t *load;
412 int nrank_gpu_shared;
414 MPI_Comm *mpi_comm_load;
415 MPI_Comm mpi_comm_gpu_shared;
418 /* Maximum DLB scaling per load balancing step in percent */
422 float cycl[ddCyclNr];
423 int cycl_n[ddCyclNr];
424 float cycl_max[ddCyclNr];
425 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
429 /* How many times have did we have load measurements */
431 /* How many times have we collected the load measurements */
435 double sum_nat[ddnatNR-ddnatZONE];
445 /* The last partition step */
446 gmx_int64_t partition_step;
454 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
457 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
458 #define DD_FLAG_NRCG 65535
459 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
460 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
462 /* Zone permutation required to obtain consecutive charge groups
463 * for neighbor searching.
465 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
467 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
468 * components see only j zones with that component 0.
471 /* The DD zone order */
472 static const ivec dd_zo[DD_MAXZONE] =
473 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
478 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
483 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
488 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
493 static const ivec dd_zp0[dd_zp0n] = {{0, 0, 1}};
495 /* Factors used to avoid problems due to rounding issues */
496 #define DD_CELL_MARGIN 1.0001
497 #define DD_CELL_MARGIN2 1.00005
498 /* Factor to account for pressure scaling during nstlist steps */
499 #define DD_PRES_SCALE_MARGIN 1.02
501 /* Turn on DLB when the load imbalance causes this amount of total loss.
502 * There is a bit of overhead with DLB and it's difficult to achieve
503 * a load imbalance of less than 2% with DLB.
505 #define DD_PERF_LOSS_DLB_ON 0.02
507 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
508 #define DD_PERF_LOSS_WARN 0.05
510 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
512 /* Use separate MPI send and receive commands
513 * when nnodes <= GMX_DD_NNODES_SENDRECV.
514 * This saves memory (and some copying for small nnodes).
515 * For high parallelization scatter and gather calls are used.
517 #define GMX_DD_NNODES_SENDRECV 4
521 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
523 static void index2xyz(ivec nc,int ind,ivec xyz)
525 xyz[XX] = ind % nc[XX];
526 xyz[YY] = (ind / nc[XX]) % nc[YY];
527 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
531 /* This order is required to minimize the coordinate communication in PME
532 * which uses decomposition in the x direction.
534 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
536 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
538 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
539 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
540 xyz[ZZ] = ind % nc[ZZ];
543 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
548 ddindex = dd_index(dd->nc, c);
549 if (dd->comm->bCartesianPP_PME)
551 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
553 else if (dd->comm->bCartesianPP)
556 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
567 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
569 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
572 int ddglatnr(gmx_domdec_t *dd, int i)
582 if (i >= dd->comm->nat[ddnatNR-1])
584 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
586 atnr = dd->gatindex[i] + 1;
592 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
594 return &dd->comm->cgs_gl;
597 static bool dlbIsOn(const gmx_domdec_comm_t *comm)
599 return (comm->dlbState == edlbsOn);
602 static void vec_rvec_init(vec_rvec_t *v)
608 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
612 v->nalloc = over_alloc_dd(n);
613 srenew(v->v, v->nalloc);
617 void dd_store_state(gmx_domdec_t *dd, t_state *state)
621 if (state->ddp_count != dd->ddp_count)
623 gmx_incons("The state does not the domain decomposition state");
626 state->ncg_gl = dd->ncg_home;
627 if (state->ncg_gl > state->cg_gl_nalloc)
629 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
630 srenew(state->cg_gl, state->cg_gl_nalloc);
632 for (i = 0; i < state->ncg_gl; i++)
634 state->cg_gl[i] = dd->index_gl[i];
637 state->ddp_count_cg_gl = dd->ddp_count;
640 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
642 return &dd->comm->zones;
645 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
646 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
648 gmx_domdec_zones_t *zones;
651 zones = &dd->comm->zones;
654 while (icg >= zones->izone[izone].cg1)
663 else if (izone < zones->nizone)
665 *jcg0 = zones->izone[izone].jcg0;
669 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
670 icg, izone, zones->nizone);
673 *jcg1 = zones->izone[izone].jcg1;
675 for (d = 0; d < dd->ndim; d++)
678 shift0[dim] = zones->izone[izone].shift0[dim];
679 shift1[dim] = zones->izone[izone].shift1[dim];
680 if (dd->comm->tric_dir[dim] || (dlbIsOn(dd->comm) && d > 0))
682 /* A conservative approach, this can be optimized */
689 int dd_natoms_vsite(gmx_domdec_t *dd)
691 return dd->comm->nat[ddnatVSITE];
694 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
696 *at_start = dd->comm->nat[ddnatCON-1];
697 *at_end = dd->comm->nat[ddnatCON];
700 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
702 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
703 int *index, *cgindex;
704 gmx_domdec_comm_t *comm;
705 gmx_domdec_comm_dim_t *cd;
706 gmx_domdec_ind_t *ind;
707 rvec shift = {0, 0, 0}, *buf, *rbuf;
708 gmx_bool bPBC, bScrew;
712 cgindex = dd->cgindex;
717 nat_tot = dd->nat_home;
718 for (d = 0; d < dd->ndim; d++)
720 bPBC = (dd->ci[dd->dim[d]] == 0);
721 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
724 copy_rvec(box[dd->dim[d]], shift);
727 for (p = 0; p < cd->np; p++)
734 for (i = 0; i < ind->nsend[nzone]; i++)
736 at0 = cgindex[index[i]];
737 at1 = cgindex[index[i]+1];
738 for (j = at0; j < at1; j++)
740 copy_rvec(x[j], buf[n]);
747 for (i = 0; i < ind->nsend[nzone]; i++)
749 at0 = cgindex[index[i]];
750 at1 = cgindex[index[i]+1];
751 for (j = at0; j < at1; j++)
753 /* We need to shift the coordinates */
754 rvec_add(x[j], shift, buf[n]);
761 for (i = 0; i < ind->nsend[nzone]; i++)
763 at0 = cgindex[index[i]];
764 at1 = cgindex[index[i]+1];
765 for (j = at0; j < at1; j++)
768 buf[n][XX] = x[j][XX] + shift[XX];
770 * This operation requires a special shift force
771 * treatment, which is performed in calc_vir.
773 buf[n][YY] = box[YY][YY] - x[j][YY];
774 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
786 rbuf = comm->vbuf2.v;
788 /* Send and receive the coordinates */
789 dd_sendrecv_rvec(dd, d, dddirBackward,
790 buf, ind->nsend[nzone+1],
791 rbuf, ind->nrecv[nzone+1]);
795 for (zone = 0; zone < nzone; zone++)
797 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
799 copy_rvec(rbuf[j], x[i]);
804 nat_tot += ind->nrecv[nzone+1];
810 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
812 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
813 int *index, *cgindex;
814 gmx_domdec_comm_t *comm;
815 gmx_domdec_comm_dim_t *cd;
816 gmx_domdec_ind_t *ind;
820 gmx_bool bShiftForcesNeedPbc, bScrew;
824 cgindex = dd->cgindex;
828 nzone = comm->zones.n/2;
829 nat_tot = dd->nat_tot;
830 for (d = dd->ndim-1; d >= 0; d--)
832 /* Only forces in domains near the PBC boundaries need to
833 consider PBC in the treatment of fshift */
834 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
835 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
836 if (fshift == NULL && !bScrew)
838 bShiftForcesNeedPbc = FALSE;
840 /* Determine which shift vector we need */
846 for (p = cd->np-1; p >= 0; p--)
849 nat_tot -= ind->nrecv[nzone+1];
856 sbuf = comm->vbuf2.v;
858 for (zone = 0; zone < nzone; zone++)
860 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
862 copy_rvec(f[i], sbuf[j]);
867 /* Communicate the forces */
868 dd_sendrecv_rvec(dd, d, dddirForward,
869 sbuf, ind->nrecv[nzone+1],
870 buf, ind->nsend[nzone+1]);
872 /* Add the received forces */
874 if (!bShiftForcesNeedPbc)
876 for (i = 0; i < ind->nsend[nzone]; i++)
878 at0 = cgindex[index[i]];
879 at1 = cgindex[index[i]+1];
880 for (j = at0; j < at1; j++)
882 rvec_inc(f[j], buf[n]);
889 /* fshift should always be defined if this function is
890 * called when bShiftForcesNeedPbc is true */
891 assert(NULL != fshift);
892 for (i = 0; i < ind->nsend[nzone]; i++)
894 at0 = cgindex[index[i]];
895 at1 = cgindex[index[i]+1];
896 for (j = at0; j < at1; j++)
898 rvec_inc(f[j], buf[n]);
899 /* Add this force to the shift force */
900 rvec_inc(fshift[is], buf[n]);
907 for (i = 0; i < ind->nsend[nzone]; i++)
909 at0 = cgindex[index[i]];
910 at1 = cgindex[index[i]+1];
911 for (j = at0; j < at1; j++)
913 /* Rotate the force */
914 f[j][XX] += buf[n][XX];
915 f[j][YY] -= buf[n][YY];
916 f[j][ZZ] -= buf[n][ZZ];
919 /* Add this force to the shift force */
920 rvec_inc(fshift[is], buf[n]);
931 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
933 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
934 int *index, *cgindex;
935 gmx_domdec_comm_t *comm;
936 gmx_domdec_comm_dim_t *cd;
937 gmx_domdec_ind_t *ind;
942 cgindex = dd->cgindex;
944 buf = &comm->vbuf.v[0][0];
947 nat_tot = dd->nat_home;
948 for (d = 0; d < dd->ndim; d++)
951 for (p = 0; p < cd->np; p++)
956 for (i = 0; i < ind->nsend[nzone]; i++)
958 at0 = cgindex[index[i]];
959 at1 = cgindex[index[i]+1];
960 for (j = at0; j < at1; j++)
973 rbuf = &comm->vbuf2.v[0][0];
975 /* Send and receive the coordinates */
976 dd_sendrecv_real(dd, d, dddirBackward,
977 buf, ind->nsend[nzone+1],
978 rbuf, ind->nrecv[nzone+1]);
982 for (zone = 0; zone < nzone; zone++)
984 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
991 nat_tot += ind->nrecv[nzone+1];
997 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
999 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
1000 int *index, *cgindex;
1001 gmx_domdec_comm_t *comm;
1002 gmx_domdec_comm_dim_t *cd;
1003 gmx_domdec_ind_t *ind;
1008 cgindex = dd->cgindex;
1010 buf = &comm->vbuf.v[0][0];
1012 nzone = comm->zones.n/2;
1013 nat_tot = dd->nat_tot;
1014 for (d = dd->ndim-1; d >= 0; d--)
1017 for (p = cd->np-1; p >= 0; p--)
1020 nat_tot -= ind->nrecv[nzone+1];
1027 sbuf = &comm->vbuf2.v[0][0];
1029 for (zone = 0; zone < nzone; zone++)
1031 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
1038 /* Communicate the forces */
1039 dd_sendrecv_real(dd, d, dddirForward,
1040 sbuf, ind->nrecv[nzone+1],
1041 buf, ind->nsend[nzone+1]);
1043 /* Add the received forces */
1045 for (i = 0; i < ind->nsend[nzone]; i++)
1047 at0 = cgindex[index[i]];
1048 at1 = cgindex[index[i]+1];
1049 for (j = at0; j < at1; j++)
1060 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1062 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",
1064 zone->min0, zone->max1,
1065 zone->mch0, zone->mch0,
1066 zone->p1_0, zone->p1_1);
1070 #define DDZONECOMM_MAXZONE 5
1071 #define DDZONECOMM_BUFSIZE 3
1073 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1074 int ddimind, int direction,
1075 gmx_ddzone_t *buf_s, int n_s,
1076 gmx_ddzone_t *buf_r, int n_r)
1078 #define ZBS DDZONECOMM_BUFSIZE
1079 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1080 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1083 for (i = 0; i < n_s; i++)
1085 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1086 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1087 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1088 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1089 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1090 vbuf_s[i*ZBS+1][2] = 0;
1091 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1092 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1093 vbuf_s[i*ZBS+2][2] = 0;
1096 dd_sendrecv_rvec(dd, ddimind, direction,
1100 for (i = 0; i < n_r; i++)
1102 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1103 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1104 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1105 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1106 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1107 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1108 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1114 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1115 rvec cell_ns_x0, rvec cell_ns_x1)
1117 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1119 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1120 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1121 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1122 rvec extr_s[2], extr_r[2];
1124 real dist_d, c = 0, det;
1125 gmx_domdec_comm_t *comm;
1126 gmx_bool bPBC, bUse;
1130 for (d = 1; d < dd->ndim; d++)
1133 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1134 zp->min0 = cell_ns_x0[dim];
1135 zp->max1 = cell_ns_x1[dim];
1136 zp->min1 = cell_ns_x1[dim];
1137 zp->mch0 = cell_ns_x0[dim];
1138 zp->mch1 = cell_ns_x1[dim];
1139 zp->p1_0 = cell_ns_x0[dim];
1140 zp->p1_1 = cell_ns_x1[dim];
1143 for (d = dd->ndim-2; d >= 0; d--)
1146 bPBC = (dim < ddbox->npbcdim);
1148 /* Use an rvec to store two reals */
1149 extr_s[d][0] = comm->cell_f0[d+1];
1150 extr_s[d][1] = comm->cell_f1[d+1];
1151 extr_s[d][2] = comm->cell_f1[d+1];
1154 /* Store the extremes in the backward sending buffer,
1155 * so the get updated separately from the forward communication.
1157 for (d1 = d; d1 < dd->ndim-1; d1++)
1159 /* We invert the order to be able to use the same loop for buf_e */
1160 buf_s[pos].min0 = extr_s[d1][1];
1161 buf_s[pos].max1 = extr_s[d1][0];
1162 buf_s[pos].min1 = extr_s[d1][2];
1163 buf_s[pos].mch0 = 0;
1164 buf_s[pos].mch1 = 0;
1165 /* Store the cell corner of the dimension we communicate along */
1166 buf_s[pos].p1_0 = comm->cell_x0[dim];
1167 buf_s[pos].p1_1 = 0;
1171 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1174 if (dd->ndim == 3 && d == 0)
1176 buf_s[pos] = comm->zone_d2[0][1];
1178 buf_s[pos] = comm->zone_d1[0];
1182 /* We only need to communicate the extremes
1183 * in the forward direction
1185 npulse = comm->cd[d].np;
1188 /* Take the minimum to avoid double communication */
1189 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1193 /* Without PBC we should really not communicate over
1194 * the boundaries, but implementing that complicates
1195 * the communication setup and therefore we simply
1196 * do all communication, but ignore some data.
1198 npulse_min = npulse;
1200 for (p = 0; p < npulse_min; p++)
1202 /* Communicate the extremes forward */
1203 bUse = (bPBC || dd->ci[dim] > 0);
1205 dd_sendrecv_rvec(dd, d, dddirForward,
1206 extr_s+d, dd->ndim-d-1,
1207 extr_r+d, dd->ndim-d-1);
1211 for (d1 = d; d1 < dd->ndim-1; d1++)
1213 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1214 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1215 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1221 for (p = 0; p < npulse; p++)
1223 /* Communicate all the zone information backward */
1224 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1226 dd_sendrecv_ddzone(dd, d, dddirBackward,
1233 for (d1 = d+1; d1 < dd->ndim; d1++)
1235 /* Determine the decrease of maximum required
1236 * communication height along d1 due to the distance along d,
1237 * this avoids a lot of useless atom communication.
1239 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1241 if (ddbox->tric_dir[dim])
1243 /* c is the off-diagonal coupling between the cell planes
1244 * along directions d and d1.
1246 c = ddbox->v[dim][dd->dim[d1]][dim];
1252 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1255 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1259 /* A negative value signals out of range */
1265 /* Accumulate the extremes over all pulses */
1266 for (i = 0; i < buf_size; i++)
1270 buf_e[i] = buf_r[i];
1276 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1277 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1278 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1281 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1289 if (bUse && dh[d1] >= 0)
1291 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1292 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1295 /* Copy the received buffer to the send buffer,
1296 * to pass the data through with the next pulse.
1298 buf_s[i] = buf_r[i];
1300 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1301 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1303 /* Store the extremes */
1306 for (d1 = d; d1 < dd->ndim-1; d1++)
1308 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1309 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1310 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1314 if (d == 1 || (d == 0 && dd->ndim == 3))
1316 for (i = d; i < 2; i++)
1318 comm->zone_d2[1-d][i] = buf_e[pos];
1324 comm->zone_d1[1] = buf_e[pos];
1334 for (i = 0; i < 2; i++)
1338 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1340 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1341 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1347 for (i = 0; i < 2; i++)
1349 for (j = 0; j < 2; j++)
1353 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1355 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1356 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1360 for (d = 1; d < dd->ndim; d++)
1362 comm->cell_f_max0[d] = extr_s[d-1][0];
1363 comm->cell_f_min1[d] = extr_s[d-1][1];
1366 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1367 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1372 static void dd_collect_cg(gmx_domdec_t *dd,
1373 t_state *state_local)
1375 gmx_domdec_master_t *ma = NULL;
1376 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1378 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1380 /* The master has the correct distribution */
1384 if (state_local->ddp_count == dd->ddp_count)
1386 /* The local state and DD are in sync, use the DD indices */
1387 ncg_home = dd->ncg_home;
1389 nat_home = dd->nat_home;
1391 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1393 /* The DD is out of sync with the local state, but we have stored
1394 * the cg indices with the local state, so we can use those.
1398 cgs_gl = &dd->comm->cgs_gl;
1400 ncg_home = state_local->ncg_gl;
1401 cg = state_local->cg_gl;
1403 for (i = 0; i < ncg_home; i++)
1405 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1410 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1424 /* Collect the charge group and atom counts on the master */
1425 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1430 for (i = 0; i < dd->nnodes; i++)
1432 ma->ncg[i] = ma->ibuf[2*i];
1433 ma->nat[i] = ma->ibuf[2*i+1];
1434 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1437 /* Make byte counts and indices */
1438 for (i = 0; i < dd->nnodes; i++)
1440 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1441 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1445 fprintf(debug, "Initial charge group distribution: ");
1446 for (i = 0; i < dd->nnodes; i++)
1448 fprintf(debug, " %d", ma->ncg[i]);
1450 fprintf(debug, "\n");
1454 /* Collect the charge group indices on the master */
1456 ncg_home*sizeof(int), cg,
1457 DDMASTER(dd) ? ma->ibuf : NULL,
1458 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1459 DDMASTER(dd) ? ma->cg : NULL);
1461 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1464 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1467 gmx_domdec_master_t *ma;
1468 int n, i, c, a, nalloc = 0;
1477 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1478 dd->rank, dd->mpi_comm_all);
1483 /* Copy the master coordinates to the global array */
1484 cgs_gl = &dd->comm->cgs_gl;
1486 n = DDMASTERRANK(dd);
1488 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1490 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1492 copy_rvec(lv[a++], v[c]);
1496 for (n = 0; n < dd->nnodes; n++)
1500 if (ma->nat[n] > nalloc)
1502 nalloc = over_alloc_dd(ma->nat[n]);
1503 srenew(buf, nalloc);
1506 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1507 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1510 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1512 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1514 copy_rvec(buf[a++], v[c]);
1523 static void get_commbuffer_counts(gmx_domdec_t *dd,
1524 int **counts, int **disps)
1526 gmx_domdec_master_t *ma;
1531 /* Make the rvec count and displacment arrays */
1533 *disps = ma->ibuf + dd->nnodes;
1534 for (n = 0; n < dd->nnodes; n++)
1536 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1537 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1541 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1544 gmx_domdec_master_t *ma;
1545 int *rcounts = NULL, *disps = NULL;
1554 get_commbuffer_counts(dd, &rcounts, &disps);
1559 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1563 cgs_gl = &dd->comm->cgs_gl;
1566 for (n = 0; n < dd->nnodes; n++)
1568 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1570 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1572 copy_rvec(buf[a++], v[c]);
1579 void dd_collect_vec(gmx_domdec_t *dd,
1580 t_state *state_local, rvec *lv, rvec *v)
1582 dd_collect_cg(dd, state_local);
1584 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1586 dd_collect_vec_sendrecv(dd, lv, v);
1590 dd_collect_vec_gatherv(dd, lv, v);
1595 void dd_collect_state(gmx_domdec_t *dd,
1596 t_state *state_local, t_state *state)
1600 nh = state->nhchainlength;
1604 for (i = 0; i < efptNR; i++)
1606 state->lambda[i] = state_local->lambda[i];
1608 state->fep_state = state_local->fep_state;
1609 state->veta = state_local->veta;
1610 state->vol0 = state_local->vol0;
1611 copy_mat(state_local->box, state->box);
1612 copy_mat(state_local->boxv, state->boxv);
1613 copy_mat(state_local->svir_prev, state->svir_prev);
1614 copy_mat(state_local->fvir_prev, state->fvir_prev);
1615 copy_mat(state_local->pres_prev, state->pres_prev);
1617 for (i = 0; i < state_local->ngtc; i++)
1619 for (j = 0; j < nh; j++)
1621 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1622 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1624 state->therm_integral[i] = state_local->therm_integral[i];
1626 for (i = 0; i < state_local->nnhpres; i++)
1628 for (j = 0; j < nh; j++)
1630 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1631 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1635 for (est = 0; est < estNR; est++)
1637 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1642 dd_collect_vec(dd, state_local, state_local->x, state->x);
1645 dd_collect_vec(dd, state_local, state_local->v, state->v);
1648 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1651 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1653 case estDISRE_INITF:
1654 case estDISRE_RM3TAV:
1655 case estORIRE_INITF:
1659 gmx_incons("Unknown state entry encountered in dd_collect_state");
1665 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1671 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1674 state->nalloc = over_alloc_dd(nalloc);
1676 for (est = 0; est < estNR; est++)
1678 if (EST_DISTR(est) && (state->flags & (1<<est)))
1683 srenew(state->x, state->nalloc);
1686 srenew(state->v, state->nalloc);
1689 srenew(state->sd_X, state->nalloc);
1692 srenew(state->cg_p, state->nalloc);
1694 case estDISRE_INITF:
1695 case estDISRE_RM3TAV:
1696 case estORIRE_INITF:
1698 /* No reallocation required */
1701 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1708 srenew(*f, state->nalloc);
1712 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1715 if (nalloc > fr->cg_nalloc)
1719 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1721 fr->cg_nalloc = over_alloc_dd(nalloc);
1722 srenew(fr->cginfo, fr->cg_nalloc);
1723 if (fr->cutoff_scheme == ecutsGROUP)
1725 srenew(fr->cg_cm, fr->cg_nalloc);
1728 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1730 /* We don't use charge groups, we use x in state to set up
1731 * the atom communication.
1733 dd_realloc_state(state, f, nalloc);
1737 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1740 gmx_domdec_master_t *ma;
1741 int n, i, c, a, nalloc = 0;
1748 for (n = 0; n < dd->nnodes; n++)
1752 if (ma->nat[n] > nalloc)
1754 nalloc = over_alloc_dd(ma->nat[n]);
1755 srenew(buf, nalloc);
1757 /* Use lv as a temporary buffer */
1759 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1761 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1763 copy_rvec(v[c], buf[a++]);
1766 if (a != ma->nat[n])
1768 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1773 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1774 DDRANK(dd, n), n, dd->mpi_comm_all);
1779 n = DDMASTERRANK(dd);
1781 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1783 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1785 copy_rvec(v[c], lv[a++]);
1792 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1793 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1798 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1801 gmx_domdec_master_t *ma;
1802 int *scounts = NULL, *disps = NULL;
1810 get_commbuffer_counts(dd, &scounts, &disps);
1814 for (n = 0; n < dd->nnodes; n++)
1816 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1818 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1820 copy_rvec(v[c], buf[a++]);
1826 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1829 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1831 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1833 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1837 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1841 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1844 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1845 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1846 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1848 if (dfhist->nlambda > 0)
1850 int nlam = dfhist->nlambda;
1851 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1852 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1853 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1854 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1855 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1856 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1858 for (i = 0; i < nlam; i++)
1860 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1861 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1862 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1863 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1864 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1865 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1870 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1871 t_state *state, t_state *state_local,
1876 nh = state->nhchainlength;
1880 for (i = 0; i < efptNR; i++)
1882 state_local->lambda[i] = state->lambda[i];
1884 state_local->fep_state = state->fep_state;
1885 state_local->veta = state->veta;
1886 state_local->vol0 = state->vol0;
1887 copy_mat(state->box, state_local->box);
1888 copy_mat(state->box_rel, state_local->box_rel);
1889 copy_mat(state->boxv, state_local->boxv);
1890 copy_mat(state->svir_prev, state_local->svir_prev);
1891 copy_mat(state->fvir_prev, state_local->fvir_prev);
1892 copy_df_history(&state_local->dfhist, &state->dfhist);
1893 for (i = 0; i < state_local->ngtc; i++)
1895 for (j = 0; j < nh; j++)
1897 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1898 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1900 state_local->therm_integral[i] = state->therm_integral[i];
1902 for (i = 0; i < state_local->nnhpres; i++)
1904 for (j = 0; j < nh; j++)
1906 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1907 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1911 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1912 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1913 dd_bcast(dd, sizeof(real), &state_local->veta);
1914 dd_bcast(dd, sizeof(real), &state_local->vol0);
1915 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1916 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1917 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1918 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1919 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1920 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1921 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1922 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1923 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1924 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1926 /* communicate df_history -- required for restarting from checkpoint */
1927 dd_distribute_dfhist(dd, &state_local->dfhist);
1929 if (dd->nat_home > state_local->nalloc)
1931 dd_realloc_state(state_local, f, dd->nat_home);
1933 for (i = 0; i < estNR; i++)
1935 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1940 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1943 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1946 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1949 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1951 case estDISRE_INITF:
1952 case estDISRE_RM3TAV:
1953 case estORIRE_INITF:
1955 /* Not implemented yet */
1958 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1964 static char dim2char(int dim)
1970 case XX: c = 'X'; break;
1971 case YY: c = 'Y'; break;
1972 case ZZ: c = 'Z'; break;
1973 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1979 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1980 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1982 rvec grid_s[2], *grid_r = NULL, cx, r;
1983 char fname[STRLEN], buf[22];
1985 int a, i, d, z, y, x;
1989 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1990 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1994 snew(grid_r, 2*dd->nnodes);
1997 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
2001 for (d = 0; d < DIM; d++)
2003 for (i = 0; i < DIM; i++)
2011 if (d < ddbox->npbcdim && dd->nc[d] > 1)
2013 tric[d][i] = box[i][d]/box[i][i];
2022 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2023 out = gmx_fio_fopen(fname, "w");
2024 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2026 for (i = 0; i < dd->nnodes; i++)
2028 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2029 for (d = 0; d < DIM; d++)
2031 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2033 for (z = 0; z < 2; z++)
2035 for (y = 0; y < 2; y++)
2037 for (x = 0; x < 2; x++)
2039 cx[XX] = grid_r[i*2+x][XX];
2040 cx[YY] = grid_r[i*2+y][YY];
2041 cx[ZZ] = grid_r[i*2+z][ZZ];
2043 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2044 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2048 for (d = 0; d < DIM; d++)
2050 for (x = 0; x < 4; x++)
2054 case 0: y = 1 + i*8 + 2*x; break;
2055 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2056 case 2: y = 1 + i*8 + x; break;
2058 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2062 gmx_fio_fclose(out);
2067 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2068 gmx_mtop_t *mtop, t_commrec *cr,
2069 int natoms, rvec x[], matrix box)
2071 char fname[STRLEN], buf[22];
2073 int i, ii, resnr, c;
2074 char *atomname, *resname;
2081 natoms = dd->comm->nat[ddnatVSITE];
2084 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2086 out = gmx_fio_fopen(fname, "w");
2088 fprintf(out, "TITLE %s\n", title);
2089 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2090 for (i = 0; i < natoms; i++)
2092 ii = dd->gatindex[i];
2093 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2094 if (i < dd->comm->nat[ddnatZONE])
2097 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2103 else if (i < dd->comm->nat[ddnatVSITE])
2105 b = dd->comm->zones.n;
2109 b = dd->comm->zones.n + 1;
2111 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2112 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2114 fprintf(out, "TER\n");
2116 gmx_fio_fclose(out);
2119 real dd_cutoff_multibody(const gmx_domdec_t *dd)
2121 gmx_domdec_comm_t *comm;
2128 if (comm->bInterCGBondeds)
2130 if (comm->cutoff_mbody > 0)
2132 r = comm->cutoff_mbody;
2136 /* cutoff_mbody=0 means we do not have DLB */
2137 r = comm->cellsize_min[dd->dim[0]];
2138 for (di = 1; di < dd->ndim; di++)
2140 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2142 if (comm->bBondComm)
2144 r = std::max(r, comm->cutoff_mbody);
2148 r = std::min(r, comm->cutoff);
2156 real dd_cutoff_twobody(const gmx_domdec_t *dd)
2160 r_mb = dd_cutoff_multibody(dd);
2162 return std::max(dd->comm->cutoff, r_mb);
2166 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2170 nc = dd->nc[dd->comm->cartpmedim];
2171 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2172 copy_ivec(coord, coord_pme);
2173 coord_pme[dd->comm->cartpmedim] =
2174 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2177 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2179 /* Here we assign a PME node to communicate with this DD node
2180 * by assuming that the major index of both is x.
2181 * We add cr->npmenodes/2 to obtain an even distribution.
2183 return (ddindex*npme + npme/2)/ndd;
2186 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2188 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2191 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2193 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2196 static int *dd_pmenodes(t_commrec *cr)
2201 snew(pmenodes, cr->npmenodes);
2203 for (i = 0; i < cr->dd->nnodes; i++)
2205 p0 = cr_ddindex2pmeindex(cr, i);
2206 p1 = cr_ddindex2pmeindex(cr, i+1);
2207 if (i+1 == cr->dd->nnodes || p1 > p0)
2211 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2213 pmenodes[n] = i + 1 + n;
2221 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2229 if (dd->comm->bCartesian) {
2230 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2231 dd_coords2pmecoords(dd,coords,coords_pme);
2232 copy_ivec(dd->ntot,nc);
2233 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2234 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2236 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2238 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2244 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2249 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2251 gmx_domdec_comm_t *comm;
2253 int ddindex, nodeid = -1;
2255 comm = cr->dd->comm;
2260 if (comm->bCartesianPP_PME)
2263 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2268 ddindex = dd_index(cr->dd->nc, coords);
2269 if (comm->bCartesianPP)
2271 nodeid = comm->ddindex2simnodeid[ddindex];
2277 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2289 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2292 gmx_domdec_comm_t *comm;
2299 /* This assumes a uniform x domain decomposition grid cell size */
2300 if (comm->bCartesianPP_PME)
2303 ivec coord, coord_pme;
2304 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2305 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2307 /* This is a PP node */
2308 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2309 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2313 else if (comm->bCartesianPP)
2315 if (sim_nodeid < dd->nnodes)
2317 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2322 /* This assumes DD cells with identical x coordinates
2323 * are numbered sequentially.
2325 if (dd->comm->pmenodes == NULL)
2327 if (sim_nodeid < dd->nnodes)
2329 /* The DD index equals the nodeid */
2330 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2336 while (sim_nodeid > dd->comm->pmenodes[i])
2340 if (sim_nodeid < dd->comm->pmenodes[i])
2342 pmenode = dd->comm->pmenodes[i];
2350 void get_pme_nnodes(const gmx_domdec_t *dd,
2351 int *npmenodes_x, int *npmenodes_y)
2355 *npmenodes_x = dd->comm->npmenodes_x;
2356 *npmenodes_y = dd->comm->npmenodes_y;
2365 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2366 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2370 ivec coord, coord_pme;
2374 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2377 for (x = 0; x < dd->nc[XX]; x++)
2379 for (y = 0; y < dd->nc[YY]; y++)
2381 for (z = 0; z < dd->nc[ZZ]; z++)
2383 if (dd->comm->bCartesianPP_PME)
2388 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2389 if (dd->ci[XX] == coord_pme[XX] &&
2390 dd->ci[YY] == coord_pme[YY] &&
2391 dd->ci[ZZ] == coord_pme[ZZ])
2393 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2398 /* The slab corresponds to the nodeid in the PME group */
2399 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2401 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2408 /* The last PP-only node is the peer node */
2409 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2413 fprintf(debug, "Receive coordinates from PP ranks:");
2414 for (x = 0; x < *nmy_ddnodes; x++)
2416 fprintf(debug, " %d", (*my_ddnodes)[x]);
2418 fprintf(debug, "\n");
2422 static gmx_bool receive_vir_ener(t_commrec *cr)
2424 gmx_domdec_comm_t *comm;
2429 if (cr->npmenodes < cr->dd->nnodes)
2431 comm = cr->dd->comm;
2432 if (comm->bCartesianPP_PME)
2434 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2437 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2438 coords[comm->cartpmedim]++;
2439 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2442 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2443 if (dd_simnode2pmenode(cr, rank) == pmenode)
2445 /* This is not the last PP node for pmenode */
2453 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2454 if (cr->sim_nodeid+1 < cr->nnodes &&
2455 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2457 /* This is not the last PP node for pmenode */
2466 static void set_zones_ncg_home(gmx_domdec_t *dd)
2468 gmx_domdec_zones_t *zones;
2471 zones = &dd->comm->zones;
2473 zones->cg_range[0] = 0;
2474 for (i = 1; i < zones->n+1; i++)
2476 zones->cg_range[i] = dd->ncg_home;
2478 /* zone_ncg1[0] should always be equal to ncg_home */
2479 dd->comm->zone_ncg1[0] = dd->ncg_home;
2482 static void rebuild_cgindex(gmx_domdec_t *dd,
2483 const int *gcgs_index, t_state *state)
2485 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2488 dd_cg_gl = dd->index_gl;
2489 cgindex = dd->cgindex;
2492 for (i = 0; i < state->ncg_gl; i++)
2496 dd_cg_gl[i] = cg_gl;
2497 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2501 dd->ncg_home = state->ncg_gl;
2504 set_zones_ncg_home(dd);
2507 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2509 while (cg >= cginfo_mb->cg_end)
2514 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2517 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2518 t_forcerec *fr, char *bLocalCG)
2520 cginfo_mb_t *cginfo_mb;
2526 cginfo_mb = fr->cginfo_mb;
2527 cginfo = fr->cginfo;
2529 for (cg = cg0; cg < cg1; cg++)
2531 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2535 if (bLocalCG != NULL)
2537 for (cg = cg0; cg < cg1; cg++)
2539 bLocalCG[index_gl[cg]] = TRUE;
2544 static void make_dd_indices(gmx_domdec_t *dd,
2545 const int *gcgs_index, int cg_start)
2547 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2548 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2551 if (dd->nat_tot > dd->gatindex_nalloc)
2553 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2554 srenew(dd->gatindex, dd->gatindex_nalloc);
2557 nzone = dd->comm->zones.n;
2558 zone2cg = dd->comm->zones.cg_range;
2559 zone_ncg1 = dd->comm->zone_ncg1;
2560 index_gl = dd->index_gl;
2561 gatindex = dd->gatindex;
2562 bCGs = dd->comm->bCGs;
2564 if (zone2cg[1] != dd->ncg_home)
2566 gmx_incons("dd->ncg_zone is not up to date");
2569 /* Make the local to global and global to local atom index */
2570 a = dd->cgindex[cg_start];
2571 for (zone = 0; zone < nzone; zone++)
2579 cg0 = zone2cg[zone];
2581 cg1 = zone2cg[zone+1];
2582 cg1_p1 = cg0 + zone_ncg1[zone];
2584 for (cg = cg0; cg < cg1; cg++)
2589 /* Signal that this cg is from more than one pulse away */
2592 cg_gl = index_gl[cg];
2595 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2598 ga2la_set(dd->ga2la, a_gl, a, zone1);
2604 gatindex[a] = cg_gl;
2605 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2612 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2618 if (bLocalCG == NULL)
2622 for (i = 0; i < dd->ncg_tot; i++)
2624 if (!bLocalCG[dd->index_gl[i]])
2627 "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);
2632 for (i = 0; i < ncg_sys; i++)
2639 if (ngl != dd->ncg_tot)
2641 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);
2648 static void check_index_consistency(gmx_domdec_t *dd,
2649 int natoms_sys, int ncg_sys,
2652 int nerr, ngl, i, a, cell;
2657 if (dd->comm->DD_debug > 1)
2659 snew(have, natoms_sys);
2660 for (a = 0; a < dd->nat_tot; a++)
2662 if (have[dd->gatindex[a]] > 0)
2664 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);
2668 have[dd->gatindex[a]] = a + 1;
2674 snew(have, dd->nat_tot);
2677 for (i = 0; i < natoms_sys; i++)
2679 if (ga2la_get(dd->ga2la, i, &a, &cell))
2681 if (a >= dd->nat_tot)
2683 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);
2689 if (dd->gatindex[a] != i)
2691 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);
2698 if (ngl != dd->nat_tot)
2701 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2702 dd->rank, where, ngl, dd->nat_tot);
2704 for (a = 0; a < dd->nat_tot; a++)
2709 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2710 dd->rank, where, a+1, dd->gatindex[a]+1);
2715 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2719 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2720 dd->rank, where, nerr);
2724 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2731 /* Clear the whole list without searching */
2732 ga2la_clear(dd->ga2la);
2736 for (i = a_start; i < dd->nat_tot; i++)
2738 ga2la_del(dd->ga2la, dd->gatindex[i]);
2742 bLocalCG = dd->comm->bLocalCG;
2745 for (i = cg_start; i < dd->ncg_tot; i++)
2747 bLocalCG[dd->index_gl[i]] = FALSE;
2751 dd_clear_local_vsite_indices(dd);
2753 if (dd->constraints)
2755 dd_clear_local_constraint_indices(dd);
2759 /* This function should be used for moving the domain boudaries during DLB,
2760 * for obtaining the minimum cell size. It checks the initially set limit
2761 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2762 * and, possibly, a longer cut-off limit set for PME load balancing.
2764 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2768 cellsize_min = comm->cellsize_min[dim];
2770 if (!comm->bVacDLBNoLimit)
2772 /* The cut-off might have changed, e.g. by PME load balacning,
2773 * from the value used to set comm->cellsize_min, so check it.
2775 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2777 if (comm->bPMELoadBalDLBLimits)
2779 /* Check for the cut-off limit set by the PME load balancing */
2780 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2784 return cellsize_min;
2787 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2790 real grid_jump_limit;
2792 /* The distance between the boundaries of cells at distance
2793 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2794 * and by the fact that cells should not be shifted by more than
2795 * half their size, such that cg's only shift by one cell
2796 * at redecomposition.
2798 grid_jump_limit = comm->cellsize_limit;
2799 if (!comm->bVacDLBNoLimit)
2801 if (comm->bPMELoadBalDLBLimits)
2803 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2805 grid_jump_limit = std::max(grid_jump_limit,
2806 cutoff/comm->cd[dim_ind].np);
2809 return grid_jump_limit;
2812 static gmx_bool check_grid_jump(gmx_int64_t step,
2818 gmx_domdec_comm_t *comm;
2827 for (d = 1; d < dd->ndim; d++)
2830 limit = grid_jump_limit(comm, cutoff, d);
2831 bfac = ddbox->box_size[dim];
2832 if (ddbox->tric_dir[dim])
2834 bfac *= ddbox->skew_fac[dim];
2836 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2837 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2845 /* This error should never be triggered under normal
2846 * circumstances, but you never know ...
2848 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.",
2849 gmx_step_str(step, buf),
2850 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2858 static int dd_load_count(gmx_domdec_comm_t *comm)
2860 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2863 static float dd_force_load(gmx_domdec_comm_t *comm)
2870 if (comm->eFlop > 1)
2872 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2877 load = comm->cycl[ddCyclF];
2878 if (comm->cycl_n[ddCyclF] > 1)
2880 /* Subtract the maximum of the last n cycle counts
2881 * to get rid of possible high counts due to other sources,
2882 * for instance system activity, that would otherwise
2883 * affect the dynamic load balancing.
2885 load -= comm->cycl_max[ddCyclF];
2889 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2891 float gpu_wait, gpu_wait_sum;
2893 gpu_wait = comm->cycl[ddCyclWaitGPU];
2894 if (comm->cycl_n[ddCyclF] > 1)
2896 /* We should remove the WaitGPU time of the same MD step
2897 * as the one with the maximum F time, since the F time
2898 * and the wait time are not independent.
2899 * Furthermore, the step for the max F time should be chosen
2900 * the same on all ranks that share the same GPU.
2901 * But to keep the code simple, we remove the average instead.
2902 * The main reason for artificially long times at some steps
2903 * is spurious CPU activity or MPI time, so we don't expect
2904 * that changes in the GPU wait time matter a lot here.
2906 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2908 /* Sum the wait times over the ranks that share the same GPU */
2909 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2910 comm->mpi_comm_gpu_shared);
2911 /* Replace the wait time by the average over the ranks */
2912 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2920 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2922 gmx_domdec_comm_t *comm;
2927 snew(*dim_f, dd->nc[dim]+1);
2929 for (i = 1; i < dd->nc[dim]; i++)
2931 if (comm->slb_frac[dim])
2933 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2937 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2940 (*dim_f)[dd->nc[dim]] = 1;
2943 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2945 int pmeindex, slab, nso, i;
2948 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2954 ddpme->dim = dimind;
2956 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2958 ddpme->nslab = (ddpme->dim == 0 ?
2959 dd->comm->npmenodes_x :
2960 dd->comm->npmenodes_y);
2962 if (ddpme->nslab <= 1)
2967 nso = dd->comm->npmenodes/ddpme->nslab;
2968 /* Determine for each PME slab the PP location range for dimension dim */
2969 snew(ddpme->pp_min, ddpme->nslab);
2970 snew(ddpme->pp_max, ddpme->nslab);
2971 for (slab = 0; slab < ddpme->nslab; slab++)
2973 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2974 ddpme->pp_max[slab] = 0;
2976 for (i = 0; i < dd->nnodes; i++)
2978 ddindex2xyz(dd->nc, i, xyz);
2979 /* For y only use our y/z slab.
2980 * This assumes that the PME x grid size matches the DD grid size.
2982 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2984 pmeindex = ddindex2pmeindex(dd, i);
2987 slab = pmeindex/nso;
2991 slab = pmeindex % ddpme->nslab;
2993 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2994 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2998 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
3001 int dd_pme_maxshift_x(gmx_domdec_t *dd)
3003 if (dd->comm->ddpme[0].dim == XX)
3005 return dd->comm->ddpme[0].maxshift;
3013 int dd_pme_maxshift_y(gmx_domdec_t *dd)
3015 if (dd->comm->ddpme[0].dim == YY)
3017 return dd->comm->ddpme[0].maxshift;
3019 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3021 return dd->comm->ddpme[1].maxshift;
3029 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3030 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3032 gmx_domdec_comm_t *comm;
3035 real range, pme_boundary;
3039 nc = dd->nc[ddpme->dim];
3042 if (!ddpme->dim_match)
3044 /* PP decomposition is not along dim: the worst situation */
3047 else if (ns <= 3 || (bUniform && ns == nc))
3049 /* The optimal situation */
3054 /* We need to check for all pme nodes which nodes they
3055 * could possibly need to communicate with.
3057 xmin = ddpme->pp_min;
3058 xmax = ddpme->pp_max;
3059 /* Allow for atoms to be maximally 2/3 times the cut-off
3060 * out of their DD cell. This is a reasonable balance between
3061 * between performance and support for most charge-group/cut-off
3064 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3065 /* Avoid extra communication when we are exactly at a boundary */
3069 for (s = 0; s < ns; s++)
3071 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3072 pme_boundary = (real)s/ns;
3075 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3077 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3081 pme_boundary = (real)(s+1)/ns;
3084 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3086 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3093 ddpme->maxshift = sh;
3097 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3098 ddpme->dim, ddpme->maxshift);
3102 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3106 for (d = 0; d < dd->ndim; d++)
3109 if (dim < ddbox->nboundeddim &&
3110 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3111 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3113 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",
3114 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3115 dd->nc[dim], dd->comm->cellsize_limit);
3121 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3124 /* Set the domain boundaries. Use for static (or no) load balancing,
3125 * and also for the starting state for dynamic load balancing.
3126 * setmode determine if and where the boundaries are stored, use enum above.
3127 * Returns the number communication pulses in npulse.
3129 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3130 int setmode, ivec npulse)
3132 gmx_domdec_comm_t *comm;
3135 real *cell_x, cell_dx, cellsize;
3139 for (d = 0; d < DIM; d++)
3141 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3143 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3146 cell_dx = ddbox->box_size[d]/dd->nc[d];
3149 case setcellsizeslbMASTER:
3150 for (j = 0; j < dd->nc[d]+1; j++)
3152 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3155 case setcellsizeslbLOCAL:
3156 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3157 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3162 cellsize = cell_dx*ddbox->skew_fac[d];
3163 while (cellsize*npulse[d] < comm->cutoff)
3167 cellsize_min[d] = cellsize;
3171 /* Statically load balanced grid */
3172 /* Also when we are not doing a master distribution we determine
3173 * all cell borders in a loop to obtain identical values
3174 * to the master distribution case and to determine npulse.
3176 if (setmode == setcellsizeslbMASTER)
3178 cell_x = dd->ma->cell_x[d];
3182 snew(cell_x, dd->nc[d]+1);
3184 cell_x[0] = ddbox->box0[d];
3185 for (j = 0; j < dd->nc[d]; j++)
3187 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3188 cell_x[j+1] = cell_x[j] + cell_dx;
3189 cellsize = cell_dx*ddbox->skew_fac[d];
3190 while (cellsize*npulse[d] < comm->cutoff &&
3191 npulse[d] < dd->nc[d]-1)
3195 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3197 if (setmode == setcellsizeslbLOCAL)
3199 comm->cell_x0[d] = cell_x[dd->ci[d]];
3200 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3202 if (setmode != setcellsizeslbMASTER)
3207 /* The following limitation is to avoid that a cell would receive
3208 * some of its own home charge groups back over the periodic boundary.
3209 * Double charge groups cause trouble with the global indices.
3211 if (d < ddbox->npbcdim &&
3212 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3214 char error_string[STRLEN];
3216 sprintf(error_string,
3217 "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",
3218 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3220 dd->nc[d], dd->nc[d],
3221 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3223 if (setmode == setcellsizeslbLOCAL)
3225 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3229 gmx_fatal(FARGS, error_string);
3236 copy_rvec(cellsize_min, comm->cellsize_min);
3239 for (d = 0; d < comm->npmedecompdim; d++)
3241 set_pme_maxshift(dd, &comm->ddpme[d],
3242 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3243 comm->ddpme[d].slb_dim_f);
3248 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3249 int d, int dim, gmx_domdec_root_t *root,
3251 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3253 gmx_domdec_comm_t *comm;
3254 int ncd, i, j, nmin, nmin_old;
3255 gmx_bool bLimLo, bLimHi;
3257 real fac, halfway, cellsize_limit_f_i, region_size;
3258 gmx_bool bPBC, bLastHi = FALSE;
3259 int nrange[] = {range[0], range[1]};
3261 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3267 bPBC = (dim < ddbox->npbcdim);
3269 cell_size = root->buf_ncd;
3273 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3276 /* First we need to check if the scaling does not make cells
3277 * smaller than the smallest allowed size.
3278 * We need to do this iteratively, since if a cell is too small,
3279 * it needs to be enlarged, which makes all the other cells smaller,
3280 * which could in turn make another cell smaller than allowed.
3282 for (i = range[0]; i < range[1]; i++)
3284 root->bCellMin[i] = FALSE;
3290 /* We need the total for normalization */
3292 for (i = range[0]; i < range[1]; i++)
3294 if (root->bCellMin[i] == FALSE)
3296 fac += cell_size[i];
3299 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3300 /* Determine the cell boundaries */
3301 for (i = range[0]; i < range[1]; i++)
3303 if (root->bCellMin[i] == FALSE)
3305 cell_size[i] *= fac;
3306 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3308 cellsize_limit_f_i = 0;
3312 cellsize_limit_f_i = cellsize_limit_f;
3314 if (cell_size[i] < cellsize_limit_f_i)
3316 root->bCellMin[i] = TRUE;
3317 cell_size[i] = cellsize_limit_f_i;
3321 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3324 while (nmin > nmin_old);
3327 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3328 /* For this check we should not use DD_CELL_MARGIN,
3329 * but a slightly smaller factor,
3330 * since rounding could get use below the limit.
3332 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3335 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",
3336 gmx_step_str(step, buf),
3337 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3338 ncd, comm->cellsize_min[dim]);
3341 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3345 /* Check if the boundary did not displace more than halfway
3346 * each of the cells it bounds, as this could cause problems,
3347 * especially when the differences between cell sizes are large.
3348 * If changes are applied, they will not make cells smaller
3349 * than the cut-off, as we check all the boundaries which
3350 * might be affected by a change and if the old state was ok,
3351 * the cells will at most be shrunk back to their old size.
3353 for (i = range[0]+1; i < range[1]; i++)
3355 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3356 if (root->cell_f[i] < halfway)
3358 root->cell_f[i] = halfway;
3359 /* Check if the change also causes shifts of the next boundaries */
3360 for (j = i+1; j < range[1]; j++)
3362 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3364 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3368 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3369 if (root->cell_f[i] > halfway)
3371 root->cell_f[i] = halfway;
3372 /* Check if the change also causes shifts of the next boundaries */
3373 for (j = i-1; j >= range[0]+1; j--)
3375 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3377 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3384 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3385 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3386 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3387 * for a and b nrange is used */
3390 /* Take care of the staggering of the cell boundaries */
3393 for (i = range[0]; i < range[1]; i++)
3395 root->cell_f_max0[i] = root->cell_f[i];
3396 root->cell_f_min1[i] = root->cell_f[i+1];
3401 for (i = range[0]+1; i < range[1]; i++)
3403 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3404 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3405 if (bLimLo && bLimHi)
3407 /* Both limits violated, try the best we can */
3408 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3409 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3410 nrange[0] = range[0];
3412 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3415 nrange[1] = range[1];
3416 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3422 /* root->cell_f[i] = root->bound_min[i]; */
3423 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3426 else if (bLimHi && !bLastHi)
3429 if (nrange[1] < range[1]) /* found a LimLo before */
3431 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3432 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3433 nrange[0] = nrange[1];
3435 root->cell_f[i] = root->bound_max[i];
3437 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3439 nrange[1] = range[1];
3442 if (nrange[1] < range[1]) /* found last a LimLo */
3444 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3445 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3446 nrange[0] = nrange[1];
3447 nrange[1] = range[1];
3448 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3450 else if (nrange[0] > range[0]) /* found at least one LimHi */
3452 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3459 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3460 int d, int dim, gmx_domdec_root_t *root,
3461 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3462 gmx_bool bUniform, gmx_int64_t step)
3464 gmx_domdec_comm_t *comm;
3465 int ncd, d1, i, pos;
3467 real load_aver, load_i, imbalance, change, change_max, sc;
3468 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3472 int range[] = { 0, 0 };
3476 /* Convert the maximum change from the input percentage to a fraction */
3477 change_limit = comm->dlb_scale_lim*0.01;
3481 bPBC = (dim < ddbox->npbcdim);
3483 cell_size = root->buf_ncd;
3485 /* Store the original boundaries */
3486 for (i = 0; i < ncd+1; i++)
3488 root->old_cell_f[i] = root->cell_f[i];
3492 for (i = 0; i < ncd; i++)
3494 cell_size[i] = 1.0/ncd;
3497 else if (dd_load_count(comm) > 0)
3499 load_aver = comm->load[d].sum_m/ncd;
3501 for (i = 0; i < ncd; i++)
3503 /* Determine the relative imbalance of cell i */
3504 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3505 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3506 /* Determine the change of the cell size using underrelaxation */
3507 change = -relax*imbalance;
3508 change_max = std::max(change_max, std::max(change, -change));
3510 /* Limit the amount of scaling.
3511 * We need to use the same rescaling for all cells in one row,
3512 * otherwise the load balancing might not converge.
3515 if (change_max > change_limit)
3517 sc *= change_limit/change_max;
3519 for (i = 0; i < ncd; i++)
3521 /* Determine the relative imbalance of cell i */
3522 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3523 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3524 /* Determine the change of the cell size using underrelaxation */
3525 change = -sc*imbalance;
3526 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3530 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3531 cellsize_limit_f *= DD_CELL_MARGIN;
3532 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3533 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3534 if (ddbox->tric_dir[dim])
3536 cellsize_limit_f /= ddbox->skew_fac[dim];
3537 dist_min_f /= ddbox->skew_fac[dim];
3539 if (bDynamicBox && d > 0)
3541 dist_min_f *= DD_PRES_SCALE_MARGIN;
3543 if (d > 0 && !bUniform)
3545 /* Make sure that the grid is not shifted too much */
3546 for (i = 1; i < ncd; i++)
3548 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3550 gmx_incons("Inconsistent DD boundary staggering limits!");
3552 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3553 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3556 root->bound_min[i] += 0.5*space;
3558 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3559 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3562 root->bound_max[i] += 0.5*space;
3567 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3569 root->cell_f_max0[i-1] + dist_min_f,
3570 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3571 root->cell_f_min1[i] - dist_min_f);
3576 root->cell_f[0] = 0;
3577 root->cell_f[ncd] = 1;
3578 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3581 /* After the checks above, the cells should obey the cut-off
3582 * restrictions, but it does not hurt to check.
3584 for (i = 0; i < ncd; i++)
3588 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3589 dim, i, root->cell_f[i], root->cell_f[i+1]);
3592 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3593 root->cell_f[i+1] - root->cell_f[i] <
3594 cellsize_limit_f/DD_CELL_MARGIN)
3598 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3599 gmx_step_str(step, buf), dim2char(dim), i,
3600 (root->cell_f[i+1] - root->cell_f[i])
3601 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3606 /* Store the cell boundaries of the lower dimensions at the end */
3607 for (d1 = 0; d1 < d; d1++)
3609 root->cell_f[pos++] = comm->cell_f0[d1];
3610 root->cell_f[pos++] = comm->cell_f1[d1];
3613 if (d < comm->npmedecompdim)
3615 /* The master determines the maximum shift for
3616 * the coordinate communication between separate PME nodes.
3618 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3620 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3623 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3627 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3628 gmx_ddbox_t *ddbox, int dimind)
3630 gmx_domdec_comm_t *comm;
3635 /* Set the cell dimensions */
3636 dim = dd->dim[dimind];
3637 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3638 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3639 if (dim >= ddbox->nboundeddim)
3641 comm->cell_x0[dim] += ddbox->box0[dim];
3642 comm->cell_x1[dim] += ddbox->box0[dim];
3646 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3647 int d, int dim, real *cell_f_row,
3650 gmx_domdec_comm_t *comm;
3656 /* Each node would only need to know two fractions,
3657 * but it is probably cheaper to broadcast the whole array.
3659 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3660 0, comm->mpi_comm_load[d]);
3662 /* Copy the fractions for this dimension from the buffer */
3663 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3664 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3665 /* The whole array was communicated, so set the buffer position */
3666 pos = dd->nc[dim] + 1;
3667 for (d1 = 0; d1 <= d; d1++)
3671 /* Copy the cell fractions of the lower dimensions */
3672 comm->cell_f0[d1] = cell_f_row[pos++];
3673 comm->cell_f1[d1] = cell_f_row[pos++];
3675 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3677 /* Convert the communicated shift from float to int */
3678 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3681 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3685 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3686 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3687 gmx_bool bUniform, gmx_int64_t step)
3689 gmx_domdec_comm_t *comm;
3691 gmx_bool bRowMember, bRowRoot;
3696 for (d = 0; d < dd->ndim; d++)
3701 for (d1 = d; d1 < dd->ndim; d1++)
3703 if (dd->ci[dd->dim[d1]] > 0)
3716 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3717 ddbox, bDynamicBox, bUniform, step);
3718 cell_f_row = comm->root[d]->cell_f;
3722 cell_f_row = comm->cell_f_row;
3724 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3729 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3733 /* This function assumes the box is static and should therefore
3734 * not be called when the box has changed since the last
3735 * call to dd_partition_system.
3737 for (d = 0; d < dd->ndim; d++)
3739 relative_to_absolute_cell_bounds(dd, ddbox, d);
3745 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3746 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3747 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3748 gmx_wallcycle_t wcycle)
3750 gmx_domdec_comm_t *comm;
3757 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3758 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3759 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3761 else if (bDynamicBox)
3763 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3766 /* Set the dimensions for which no DD is used */
3767 for (dim = 0; dim < DIM; dim++)
3769 if (dd->nc[dim] == 1)
3771 comm->cell_x0[dim] = 0;
3772 comm->cell_x1[dim] = ddbox->box_size[dim];
3773 if (dim >= ddbox->nboundeddim)
3775 comm->cell_x0[dim] += ddbox->box0[dim];
3776 comm->cell_x1[dim] += ddbox->box0[dim];
3782 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3785 gmx_domdec_comm_dim_t *cd;
3787 for (d = 0; d < dd->ndim; d++)
3789 cd = &dd->comm->cd[d];
3790 np = npulse[dd->dim[d]];
3791 if (np > cd->np_nalloc)
3795 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3796 dim2char(dd->dim[d]), np);
3798 if (DDMASTER(dd) && cd->np_nalloc > 0)
3800 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3802 srenew(cd->ind, np);
3803 for (i = cd->np_nalloc; i < np; i++)
3805 cd->ind[i].index = NULL;
3806 cd->ind[i].nalloc = 0;
3815 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3816 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3817 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3818 gmx_wallcycle_t wcycle)
3820 gmx_domdec_comm_t *comm;
3826 /* Copy the old cell boundaries for the cg displacement check */
3827 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3828 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3834 check_box_size(dd, ddbox);
3836 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3840 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3841 realloc_comm_ind(dd, npulse);
3846 for (d = 0; d < DIM; d++)
3848 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3849 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3854 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3856 rvec cell_ns_x0, rvec cell_ns_x1,
3859 gmx_domdec_comm_t *comm;
3864 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3866 dim = dd->dim[dim_ind];
3868 /* Without PBC we don't have restrictions on the outer cells */
3869 if (!(dim >= ddbox->npbcdim &&
3870 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3872 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3873 comm->cellsize_min[dim])
3876 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",
3877 gmx_step_str(step, buf), dim2char(dim),
3878 comm->cell_x1[dim] - comm->cell_x0[dim],
3879 ddbox->skew_fac[dim],
3880 dd->comm->cellsize_min[dim],
3881 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3885 if ((dlbIsOn(dd->comm) && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3887 /* Communicate the boundaries and update cell_ns_x0/1 */
3888 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3889 if (dlbIsOn(dd->comm) && dd->ndim > 1)
3891 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3896 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3900 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3908 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3909 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3918 static void check_screw_box(matrix box)
3920 /* Mathematical limitation */
3921 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3923 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3926 /* Limitation due to the asymmetry of the eighth shell method */
3927 if (box[ZZ][YY] != 0)
3929 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3933 static void distribute_cg(FILE *fplog,
3934 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3937 gmx_domdec_master_t *ma;
3938 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3939 int i, icg, j, k, k0, k1, d;
3943 real nrcg, inv_ncg, pos_d;
3949 if (tmp_ind == NULL)
3951 snew(tmp_nalloc, dd->nnodes);
3952 snew(tmp_ind, dd->nnodes);
3953 for (i = 0; i < dd->nnodes; i++)
3955 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3956 snew(tmp_ind[i], tmp_nalloc[i]);
3960 /* Clear the count */
3961 for (i = 0; i < dd->nnodes; i++)
3967 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3969 cgindex = cgs->index;
3971 /* Compute the center of geometry for all charge groups */
3972 for (icg = 0; icg < cgs->nr; icg++)
3975 k1 = cgindex[icg+1];
3979 copy_rvec(pos[k0], cg_cm);
3986 for (k = k0; (k < k1); k++)
3988 rvec_inc(cg_cm, pos[k]);
3990 for (d = 0; (d < DIM); d++)
3992 cg_cm[d] *= inv_ncg;
3995 /* Put the charge group in the box and determine the cell index */
3996 for (d = DIM-1; d >= 0; d--)
3999 if (d < dd->npbcdim)
4001 bScrew = (dd->bScrewPBC && d == XX);
4002 if (tric_dir[d] && dd->nc[d] > 1)
4004 /* Use triclinic coordintates for this dimension */
4005 for (j = d+1; j < DIM; j++)
4007 pos_d += cg_cm[j]*tcm[j][d];
4010 while (pos_d >= box[d][d])
4013 rvec_dec(cg_cm, box[d]);
4016 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4017 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4019 for (k = k0; (k < k1); k++)
4021 rvec_dec(pos[k], box[d]);
4024 pos[k][YY] = box[YY][YY] - pos[k][YY];
4025 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4032 rvec_inc(cg_cm, box[d]);
4035 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4036 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4038 for (k = k0; (k < k1); k++)
4040 rvec_inc(pos[k], box[d]);
4043 pos[k][YY] = box[YY][YY] - pos[k][YY];
4044 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4049 /* This could be done more efficiently */
4051 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4056 i = dd_index(dd->nc, ind);
4057 if (ma->ncg[i] == tmp_nalloc[i])
4059 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4060 srenew(tmp_ind[i], tmp_nalloc[i]);
4062 tmp_ind[i][ma->ncg[i]] = icg;
4064 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4068 for (i = 0; i < dd->nnodes; i++)
4071 for (k = 0; k < ma->ncg[i]; k++)
4073 ma->cg[k1++] = tmp_ind[i][k];
4076 ma->index[dd->nnodes] = k1;
4078 for (i = 0; i < dd->nnodes; i++)
4087 /* Here we avoid int overflows due to #atoms^2: use double, dsqr */
4088 int nat_sum, nat_min, nat_max;
4093 nat_min = ma->nat[0];
4094 nat_max = ma->nat[0];
4095 for (i = 0; i < dd->nnodes; i++)
4097 nat_sum += ma->nat[i];
4098 nat2_sum += dsqr(ma->nat[i]);
4099 nat_min = std::min(nat_min, ma->nat[i]);
4100 nat_max = std::max(nat_max, ma->nat[i]);
4102 nat_sum /= dd->nnodes;
4103 nat2_sum /= dd->nnodes;
4105 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
4108 static_cast<int>(sqrt(nat2_sum - dsqr(nat_sum) + 0.5)),
4113 static void get_cg_distribution(FILE *fplog, gmx_domdec_t *dd,
4114 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4117 gmx_domdec_master_t *ma = NULL;
4120 int *ibuf, buf2[2] = { 0, 0 };
4121 gmx_bool bMaster = DDMASTER(dd);
4129 check_screw_box(box);
4132 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4134 distribute_cg(fplog, box, ddbox->tric_dir, cgs, pos, dd);
4135 for (i = 0; i < dd->nnodes; i++)
4137 ma->ibuf[2*i] = ma->ncg[i];
4138 ma->ibuf[2*i+1] = ma->nat[i];
4146 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4148 dd->ncg_home = buf2[0];
4149 dd->nat_home = buf2[1];
4150 dd->ncg_tot = dd->ncg_home;
4151 dd->nat_tot = dd->nat_home;
4152 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4154 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4155 srenew(dd->index_gl, dd->cg_nalloc);
4156 srenew(dd->cgindex, dd->cg_nalloc+1);
4160 for (i = 0; i < dd->nnodes; i++)
4162 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4163 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4168 bMaster ? ma->ibuf : NULL,
4169 bMaster ? ma->ibuf+dd->nnodes : NULL,
4170 bMaster ? ma->cg : NULL,
4171 dd->ncg_home*sizeof(int), dd->index_gl);
4173 /* Determine the home charge group sizes */
4175 for (i = 0; i < dd->ncg_home; i++)
4177 cg_gl = dd->index_gl[i];
4179 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4184 fprintf(debug, "Home charge groups:\n");
4185 for (i = 0; i < dd->ncg_home; i++)
4187 fprintf(debug, " %d", dd->index_gl[i]);
4190 fprintf(debug, "\n");
4193 fprintf(debug, "\n");
4197 static int compact_and_copy_vec_at(int ncg, int *move,
4200 rvec *src, gmx_domdec_comm_t *comm,
4203 int m, icg, i, i0, i1, nrcg;
4209 for (m = 0; m < DIM*2; m++)
4215 for (icg = 0; icg < ncg; icg++)
4217 i1 = cgindex[icg+1];
4223 /* Compact the home array in place */
4224 for (i = i0; i < i1; i++)
4226 copy_rvec(src[i], src[home_pos++]);
4232 /* Copy to the communication buffer */
4234 pos_vec[m] += 1 + vec*nrcg;
4235 for (i = i0; i < i1; i++)
4237 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4239 pos_vec[m] += (nvec - vec - 1)*nrcg;
4243 home_pos += i1 - i0;
4251 static int compact_and_copy_vec_cg(int ncg, int *move,
4253 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4256 int m, icg, i0, i1, nrcg;
4262 for (m = 0; m < DIM*2; m++)
4268 for (icg = 0; icg < ncg; icg++)
4270 i1 = cgindex[icg+1];
4276 /* Compact the home array in place */
4277 copy_rvec(src[icg], src[home_pos++]);
4283 /* Copy to the communication buffer */
4284 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4285 pos_vec[m] += 1 + nrcg*nvec;
4297 static int compact_ind(int ncg, int *move,
4298 int *index_gl, int *cgindex,
4300 gmx_ga2la_t ga2la, char *bLocalCG,
4303 int cg, nat, a0, a1, a, a_gl;
4308 for (cg = 0; cg < ncg; cg++)
4314 /* Compact the home arrays in place.
4315 * Anything that can be done here avoids access to global arrays.
4317 cgindex[home_pos] = nat;
4318 for (a = a0; a < a1; a++)
4321 gatindex[nat] = a_gl;
4322 /* The cell number stays 0, so we don't need to set it */
4323 ga2la_change_la(ga2la, a_gl, nat);
4326 index_gl[home_pos] = index_gl[cg];
4327 cginfo[home_pos] = cginfo[cg];
4328 /* The charge group remains local, so bLocalCG does not change */
4333 /* Clear the global indices */
4334 for (a = a0; a < a1; a++)
4336 ga2la_del(ga2la, gatindex[a]);
4340 bLocalCG[index_gl[cg]] = FALSE;
4344 cgindex[home_pos] = nat;
4349 static void clear_and_mark_ind(int ncg, int *move,
4350 int *index_gl, int *cgindex, int *gatindex,
4351 gmx_ga2la_t ga2la, char *bLocalCG,
4356 for (cg = 0; cg < ncg; cg++)
4362 /* Clear the global indices */
4363 for (a = a0; a < a1; a++)
4365 ga2la_del(ga2la, gatindex[a]);
4369 bLocalCG[index_gl[cg]] = FALSE;
4371 /* Signal that this cg has moved using the ns cell index.
4372 * Here we set it to -1. fill_grid will change it
4373 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4375 cell_index[cg] = -1;
4380 static void print_cg_move(FILE *fplog,
4382 gmx_int64_t step, int cg, int dim, int dir,
4383 gmx_bool bHaveCgcmOld, real limitd,
4384 rvec cm_old, rvec cm_new, real pos_d)
4386 gmx_domdec_comm_t *comm;
4391 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4394 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4395 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4396 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4400 /* We don't have a limiting distance available: don't print it */
4401 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4402 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4403 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4405 fprintf(fplog, "distance out of cell %f\n",
4406 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4409 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4410 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4412 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4413 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4414 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4416 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4417 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4419 comm->cell_x0[dim], comm->cell_x1[dim]);
4422 static void cg_move_error(FILE *fplog,
4424 gmx_int64_t step, int cg, int dim, int dir,
4425 gmx_bool bHaveCgcmOld, real limitd,
4426 rvec cm_old, rvec cm_new, real pos_d)
4430 print_cg_move(fplog, dd, step, cg, dim, dir,
4431 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4433 print_cg_move(stderr, dd, step, cg, dim, dir,
4434 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4436 "%s moved too far between two domain decomposition steps\n"
4437 "This usually means that your system is not well equilibrated",
4438 dd->comm->bCGs ? "A charge group" : "An atom");
4441 static void rotate_state_atom(t_state *state, int a)
4445 for (est = 0; est < estNR; est++)
4447 if (EST_DISTR(est) && (state->flags & (1<<est)))
4452 /* Rotate the complete state; for a rectangular box only */
4453 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4454 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4457 state->v[a][YY] = -state->v[a][YY];
4458 state->v[a][ZZ] = -state->v[a][ZZ];
4461 state->sd_X[a][YY] = -state->sd_X[a][YY];
4462 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4465 state->cg_p[a][YY] = -state->cg_p[a][YY];
4466 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4468 case estDISRE_INITF:
4469 case estDISRE_RM3TAV:
4470 case estORIRE_INITF:
4472 /* These are distances, so not affected by rotation */
4475 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4481 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4483 if (natoms > comm->moved_nalloc)
4485 /* Contents should be preserved here */
4486 comm->moved_nalloc = over_alloc_dd(natoms);
4487 srenew(comm->moved, comm->moved_nalloc);
4493 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4496 ivec tric_dir, matrix tcm,
4497 rvec cell_x0, rvec cell_x1,
4498 rvec limitd, rvec limit0, rvec limit1,
4500 int cg_start, int cg_end,
4505 int cg, k, k0, k1, d, dim, d2;
4510 real inv_ncg, pos_d;
4513 npbcdim = dd->npbcdim;
4515 for (cg = cg_start; cg < cg_end; cg++)
4522 copy_rvec(state->x[k0], cm_new);
4529 for (k = k0; (k < k1); k++)
4531 rvec_inc(cm_new, state->x[k]);
4533 for (d = 0; (d < DIM); d++)
4535 cm_new[d] = inv_ncg*cm_new[d];
4540 /* Do pbc and check DD cell boundary crossings */
4541 for (d = DIM-1; d >= 0; d--)
4545 bScrew = (dd->bScrewPBC && d == XX);
4546 /* Determine the location of this cg in lattice coordinates */
4550 for (d2 = d+1; d2 < DIM; d2++)
4552 pos_d += cm_new[d2]*tcm[d2][d];
4555 /* Put the charge group in the triclinic unit-cell */
4556 if (pos_d >= cell_x1[d])
4558 if (pos_d >= limit1[d])
4560 cg_move_error(fplog, dd, step, cg, d, 1,
4561 cg_cm != state->x, limitd[d],
4562 cg_cm[cg], cm_new, pos_d);
4565 if (dd->ci[d] == dd->nc[d] - 1)
4567 rvec_dec(cm_new, state->box[d]);
4570 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4571 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4573 for (k = k0; (k < k1); k++)
4575 rvec_dec(state->x[k], state->box[d]);
4578 rotate_state_atom(state, k);
4583 else if (pos_d < cell_x0[d])
4585 if (pos_d < limit0[d])
4587 cg_move_error(fplog, dd, step, cg, d, -1,
4588 cg_cm != state->x, limitd[d],
4589 cg_cm[cg], cm_new, pos_d);
4594 rvec_inc(cm_new, state->box[d]);
4597 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4598 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4600 for (k = k0; (k < k1); k++)
4602 rvec_inc(state->x[k], state->box[d]);
4605 rotate_state_atom(state, k);
4611 else if (d < npbcdim)
4613 /* Put the charge group in the rectangular unit-cell */
4614 while (cm_new[d] >= state->box[d][d])
4616 rvec_dec(cm_new, state->box[d]);
4617 for (k = k0; (k < k1); k++)
4619 rvec_dec(state->x[k], state->box[d]);
4622 while (cm_new[d] < 0)
4624 rvec_inc(cm_new, state->box[d]);
4625 for (k = k0; (k < k1); k++)
4627 rvec_inc(state->x[k], state->box[d]);
4633 copy_rvec(cm_new, cg_cm[cg]);
4635 /* Determine where this cg should go */
4638 for (d = 0; d < dd->ndim; d++)
4643 flag |= DD_FLAG_FW(d);
4649 else if (dev[dim] == -1)
4651 flag |= DD_FLAG_BW(d);
4654 if (dd->nc[dim] > 2)
4665 /* Temporarily store the flag in move */
4666 move[cg] = mc + flag;
4670 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4671 gmx_domdec_t *dd, ivec tric_dir,
4672 t_state *state, rvec **f,
4681 int ncg[DIM*2], nat[DIM*2];
4682 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4683 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4684 int sbuf[2], rbuf[2];
4685 int home_pos_cg, home_pos_at, buf_pos;
4687 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4690 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4692 cginfo_mb_t *cginfo_mb;
4693 gmx_domdec_comm_t *comm;
4695 int nthread, thread;
4699 check_screw_box(state->box);
4703 if (fr->cutoff_scheme == ecutsGROUP)
4708 for (i = 0; i < estNR; i++)
4714 case estX: /* Always present */ break;
4715 case estV: bV = (state->flags & (1<<i)); break;
4716 case estSDX: bSDX = (state->flags & (1<<i)); break;
4717 case estCGP: bCGP = (state->flags & (1<<i)); break;
4720 case estDISRE_INITF:
4721 case estDISRE_RM3TAV:
4722 case estORIRE_INITF:
4724 /* No processing required */
4727 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4732 if (dd->ncg_tot > comm->nalloc_int)
4734 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4735 srenew(comm->buf_int, comm->nalloc_int);
4737 move = comm->buf_int;
4739 /* Clear the count */
4740 for (c = 0; c < dd->ndim*2; c++)
4746 npbcdim = dd->npbcdim;
4748 for (d = 0; (d < DIM); d++)
4750 limitd[d] = dd->comm->cellsize_min[d];
4751 if (d >= npbcdim && dd->ci[d] == 0)
4753 cell_x0[d] = -GMX_FLOAT_MAX;
4757 cell_x0[d] = comm->cell_x0[d];
4759 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4761 cell_x1[d] = GMX_FLOAT_MAX;
4765 cell_x1[d] = comm->cell_x1[d];
4769 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4770 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4774 /* We check after communication if a charge group moved
4775 * more than one cell. Set the pre-comm check limit to float_max.
4777 limit0[d] = -GMX_FLOAT_MAX;
4778 limit1[d] = GMX_FLOAT_MAX;
4782 make_tric_corr_matrix(npbcdim, state->box, tcm);
4784 cgindex = dd->cgindex;
4786 nthread = gmx_omp_nthreads_get(emntDomdec);
4788 /* Compute the center of geometry for all home charge groups
4789 * and put them in the box and determine where they should go.
4791 #pragma omp parallel for num_threads(nthread) schedule(static)
4792 for (thread = 0; thread < nthread; thread++)
4794 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4795 cell_x0, cell_x1, limitd, limit0, limit1,
4797 ( thread *dd->ncg_home)/nthread,
4798 ((thread+1)*dd->ncg_home)/nthread,
4799 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4803 for (cg = 0; cg < dd->ncg_home; cg++)
4808 flag = mc & ~DD_FLAG_NRCG;
4809 mc = mc & DD_FLAG_NRCG;
4812 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4814 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4815 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4817 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4818 /* We store the cg size in the lower 16 bits
4819 * and the place where the charge group should go
4820 * in the next 6 bits. This saves some communication volume.
4822 nrcg = cgindex[cg+1] - cgindex[cg];
4823 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4829 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4830 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4833 for (i = 0; i < dd->ndim*2; i++)
4835 *ncg_moved += ncg[i];
4852 /* Make sure the communication buffers are large enough */
4853 for (mc = 0; mc < dd->ndim*2; mc++)
4855 nvr = ncg[mc] + nat[mc]*nvec;
4856 if (nvr > comm->cgcm_state_nalloc[mc])
4858 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4859 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4863 switch (fr->cutoff_scheme)
4866 /* Recalculating cg_cm might be cheaper than communicating,
4867 * but that could give rise to rounding issues.
4870 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4871 nvec, cg_cm, comm, bCompact);
4874 /* Without charge groups we send the moved atom coordinates
4875 * over twice. This is so the code below can be used without
4876 * many conditionals for both for with and without charge groups.
4879 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4880 nvec, state->x, comm, FALSE);
4883 home_pos_cg -= *ncg_moved;
4887 gmx_incons("unimplemented");
4893 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4894 nvec, vec++, state->x, comm, bCompact);
4897 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4898 nvec, vec++, state->v, comm, bCompact);
4902 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4903 nvec, vec++, state->sd_X, comm, bCompact);
4907 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4908 nvec, vec++, state->cg_p, comm, bCompact);
4913 compact_ind(dd->ncg_home, move,
4914 dd->index_gl, dd->cgindex, dd->gatindex,
4915 dd->ga2la, comm->bLocalCG,
4920 if (fr->cutoff_scheme == ecutsVERLET)
4922 moved = get_moved(comm, dd->ncg_home);
4924 for (k = 0; k < dd->ncg_home; k++)
4931 moved = fr->ns.grid->cell_index;
4934 clear_and_mark_ind(dd->ncg_home, move,
4935 dd->index_gl, dd->cgindex, dd->gatindex,
4936 dd->ga2la, comm->bLocalCG,
4940 cginfo_mb = fr->cginfo_mb;
4942 *ncg_stay_home = home_pos_cg;
4943 for (d = 0; d < dd->ndim; d++)
4948 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4951 /* Communicate the cg and atom counts */
4956 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4957 d, dir, sbuf[0], sbuf[1]);
4959 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4961 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4963 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4964 srenew(comm->buf_int, comm->nalloc_int);
4967 /* Communicate the charge group indices, sizes and flags */
4968 dd_sendrecv_int(dd, d, dir,
4969 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4970 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4972 nvs = ncg[cdd] + nat[cdd]*nvec;
4973 i = rbuf[0] + rbuf[1] *nvec;
4974 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4976 /* Communicate cgcm and state */
4977 dd_sendrecv_rvec(dd, d, dir,
4978 comm->cgcm_state[cdd], nvs,
4979 comm->vbuf.v+nvr, i);
4980 ncg_recv += rbuf[0];
4984 /* Process the received charge groups */
4986 for (cg = 0; cg < ncg_recv; cg++)
4988 flag = comm->buf_int[cg*DD_CGIBS+1];
4990 if (dim >= npbcdim && dd->nc[dim] > 2)
4992 /* No pbc in this dim and more than one domain boundary.
4993 * We do a separate check if a charge group didn't move too far.
4995 if (((flag & DD_FLAG_FW(d)) &&
4996 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4997 ((flag & DD_FLAG_BW(d)) &&
4998 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
5000 cg_move_error(fplog, dd, step, cg, dim,
5001 (flag & DD_FLAG_FW(d)) ? 1 : 0,
5002 fr->cutoff_scheme == ecutsGROUP, 0,
5003 comm->vbuf.v[buf_pos],
5004 comm->vbuf.v[buf_pos],
5005 comm->vbuf.v[buf_pos][dim]);
5012 /* Check which direction this cg should go */
5013 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
5015 if (dlbIsOn(dd->comm))
5017 /* The cell boundaries for dimension d2 are not equal
5018 * for each cell row of the lower dimension(s),
5019 * therefore we might need to redetermine where
5020 * this cg should go.
5023 /* If this cg crosses the box boundary in dimension d2
5024 * we can use the communicated flag, so we do not
5025 * have to worry about pbc.
5027 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5028 (flag & DD_FLAG_FW(d2))) ||
5029 (dd->ci[dim2] == 0 &&
5030 (flag & DD_FLAG_BW(d2)))))
5032 /* Clear the two flags for this dimension */
5033 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5034 /* Determine the location of this cg
5035 * in lattice coordinates
5037 pos_d = comm->vbuf.v[buf_pos][dim2];
5040 for (d3 = dim2+1; d3 < DIM; d3++)
5043 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5046 /* Check of we are not at the box edge.
5047 * pbc is only handled in the first step above,
5048 * but this check could move over pbc while
5049 * the first step did not due to different rounding.
5051 if (pos_d >= cell_x1[dim2] &&
5052 dd->ci[dim2] != dd->nc[dim2]-1)
5054 flag |= DD_FLAG_FW(d2);
5056 else if (pos_d < cell_x0[dim2] &&
5059 flag |= DD_FLAG_BW(d2);
5061 comm->buf_int[cg*DD_CGIBS+1] = flag;
5064 /* Set to which neighboring cell this cg should go */
5065 if (flag & DD_FLAG_FW(d2))
5069 else if (flag & DD_FLAG_BW(d2))
5071 if (dd->nc[dd->dim[d2]] > 2)
5083 nrcg = flag & DD_FLAG_NRCG;
5086 if (home_pos_cg+1 > dd->cg_nalloc)
5088 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5089 srenew(dd->index_gl, dd->cg_nalloc);
5090 srenew(dd->cgindex, dd->cg_nalloc+1);
5092 /* Set the global charge group index and size */
5093 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5094 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5095 /* Copy the state from the buffer */
5096 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5097 if (fr->cutoff_scheme == ecutsGROUP)
5100 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5104 /* Set the cginfo */
5105 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5106 dd->index_gl[home_pos_cg]);
5109 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5112 if (home_pos_at+nrcg > state->nalloc)
5114 dd_realloc_state(state, f, home_pos_at+nrcg);
5116 for (i = 0; i < nrcg; i++)
5118 copy_rvec(comm->vbuf.v[buf_pos++],
5119 state->x[home_pos_at+i]);
5123 for (i = 0; i < nrcg; i++)
5125 copy_rvec(comm->vbuf.v[buf_pos++],
5126 state->v[home_pos_at+i]);
5131 for (i = 0; i < nrcg; i++)
5133 copy_rvec(comm->vbuf.v[buf_pos++],
5134 state->sd_X[home_pos_at+i]);
5139 for (i = 0; i < nrcg; i++)
5141 copy_rvec(comm->vbuf.v[buf_pos++],
5142 state->cg_p[home_pos_at+i]);
5146 home_pos_at += nrcg;
5150 /* Reallocate the buffers if necessary */
5151 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5153 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5154 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5156 nvr = ncg[mc] + nat[mc]*nvec;
5157 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5159 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5160 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5162 /* Copy from the receive to the send buffers */
5163 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5164 comm->buf_int + cg*DD_CGIBS,
5165 DD_CGIBS*sizeof(int));
5166 memcpy(comm->cgcm_state[mc][nvr],
5167 comm->vbuf.v[buf_pos],
5168 (1+nrcg*nvec)*sizeof(rvec));
5169 buf_pos += 1 + nrcg*nvec;
5176 /* With sorting (!bCompact) the indices are now only partially up to date
5177 * and ncg_home and nat_home are not the real count, since there are
5178 * "holes" in the arrays for the charge groups that moved to neighbors.
5180 if (fr->cutoff_scheme == ecutsVERLET)
5182 moved = get_moved(comm, home_pos_cg);
5184 for (i = dd->ncg_home; i < home_pos_cg; i++)
5189 dd->ncg_home = home_pos_cg;
5190 dd->nat_home = home_pos_at;
5195 "Finished repartitioning: cgs moved out %d, new home %d\n",
5196 *ncg_moved, dd->ncg_home-*ncg_moved);
5201 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5203 dd->comm->cycl[ddCycl] += cycles;
5204 dd->comm->cycl_n[ddCycl]++;
5205 if (cycles > dd->comm->cycl_max[ddCycl])
5207 dd->comm->cycl_max[ddCycl] = cycles;
5211 static double force_flop_count(t_nrnb *nrnb)
5218 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5220 /* To get closer to the real timings, we half the count
5221 * for the normal loops and again half it for water loops.
5224 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5226 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5230 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5233 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5236 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5238 sum += nrnb->n[i]*cost_nrnb(i);
5241 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5243 sum += nrnb->n[i]*cost_nrnb(i);
5249 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5251 if (dd->comm->eFlop)
5253 dd->comm->flop -= force_flop_count(nrnb);
5256 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5258 if (dd->comm->eFlop)
5260 dd->comm->flop += force_flop_count(nrnb);
5265 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5269 for (i = 0; i < ddCyclNr; i++)
5271 dd->comm->cycl[i] = 0;
5272 dd->comm->cycl_n[i] = 0;
5273 dd->comm->cycl_max[i] = 0;
5276 dd->comm->flop_n = 0;
5279 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5281 gmx_domdec_comm_t *comm;
5282 gmx_domdec_load_t *load;
5283 gmx_domdec_root_t *root = NULL;
5285 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5290 fprintf(debug, "get_load_distribution start\n");
5293 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5297 bSepPME = (dd->pme_nodeid >= 0);
5299 if (dd->ndim == 0 && bSepPME)
5301 /* Without decomposition, but with PME nodes, we need the load */
5302 comm->load[0].mdf = comm->cycl[ddCyclPPduringPME];
5303 comm->load[0].pme = comm->cycl[ddCyclPME];
5306 for (d = dd->ndim-1; d >= 0; d--)
5309 /* Check if we participate in the communication in this dimension */
5310 if (d == dd->ndim-1 ||
5311 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5313 load = &comm->load[d];
5314 if (dlbIsOn(dd->comm))
5316 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5319 if (d == dd->ndim-1)
5321 sbuf[pos++] = dd_force_load(comm);
5322 sbuf[pos++] = sbuf[0];
5323 if (dlbIsOn(dd->comm))
5325 sbuf[pos++] = sbuf[0];
5326 sbuf[pos++] = cell_frac;
5329 sbuf[pos++] = comm->cell_f_max0[d];
5330 sbuf[pos++] = comm->cell_f_min1[d];
5335 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5336 sbuf[pos++] = comm->cycl[ddCyclPME];
5341 sbuf[pos++] = comm->load[d+1].sum;
5342 sbuf[pos++] = comm->load[d+1].max;
5343 if (dlbIsOn(dd->comm))
5345 sbuf[pos++] = comm->load[d+1].sum_m;
5346 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5347 sbuf[pos++] = comm->load[d+1].flags;
5350 sbuf[pos++] = comm->cell_f_max0[d];
5351 sbuf[pos++] = comm->cell_f_min1[d];
5356 sbuf[pos++] = comm->load[d+1].mdf;
5357 sbuf[pos++] = comm->load[d+1].pme;
5361 /* Communicate a row in DD direction d.
5362 * The communicators are setup such that the root always has rank 0.
5365 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5366 load->load, load->nload*sizeof(float), MPI_BYTE,
5367 0, comm->mpi_comm_load[d]);
5369 if (dd->ci[dim] == dd->master_ci[dim])
5371 /* We are the root, process this row */
5374 root = comm->root[d];
5384 for (i = 0; i < dd->nc[dim]; i++)
5386 load->sum += load->load[pos++];
5387 load->max = std::max(load->max, load->load[pos]);
5389 if (dlbIsOn(dd->comm))
5393 /* This direction could not be load balanced properly,
5394 * therefore we need to use the maximum iso the average load.
5396 load->sum_m = std::max(load->sum_m, load->load[pos]);
5400 load->sum_m += load->load[pos];
5403 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5407 load->flags = (int)(load->load[pos++] + 0.5);
5411 root->cell_f_max0[i] = load->load[pos++];
5412 root->cell_f_min1[i] = load->load[pos++];
5417 load->mdf = std::max(load->mdf, load->load[pos]);
5419 load->pme = std::max(load->pme, load->load[pos]);
5423 if (dlbIsOn(comm) && root->bLimited)
5425 load->sum_m *= dd->nc[dim];
5426 load->flags |= (1<<d);
5434 comm->nload += dd_load_count(comm);
5435 comm->load_step += comm->cycl[ddCyclStep];
5436 comm->load_sum += comm->load[0].sum;
5437 comm->load_max += comm->load[0].max;
5440 for (d = 0; d < dd->ndim; d++)
5442 if (comm->load[0].flags & (1<<d))
5444 comm->load_lim[d]++;
5450 comm->load_mdf += comm->load[0].mdf;
5451 comm->load_pme += comm->load[0].pme;
5455 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5459 fprintf(debug, "get_load_distribution finished\n");
5463 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5465 /* Return the relative performance loss on the total run time
5466 * due to the force calculation load imbalance.
5468 if (dd->comm->nload > 0)
5471 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5472 (dd->comm->load_step*dd->nnodes);
5480 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5483 int npp, npme, nnodes, d, limp;
5484 float imbal, pme_f_ratio, lossf = 0, lossp = 0;
5486 gmx_domdec_comm_t *comm;
5489 if (DDMASTER(dd) && comm->nload > 0)
5492 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5493 nnodes = npp + npme;
5496 imbal = comm->load_max*npp/comm->load_sum - 1;
5497 lossf = dd_force_imb_perf_loss(dd);
5498 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5499 fprintf(fplog, "%s", buf);
5500 fprintf(stderr, "\n");
5501 fprintf(stderr, "%s", buf);
5502 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5503 fprintf(fplog, "%s", buf);
5504 fprintf(stderr, "%s", buf);
5509 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5510 for (d = 0; d < dd->ndim; d++)
5512 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5513 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5519 sprintf(buf+strlen(buf), "\n");
5520 fprintf(fplog, "%s", buf);
5521 fprintf(stderr, "%s", buf);
5525 pme_f_ratio = comm->load_pme/comm->load_mdf;
5526 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5529 lossp *= (float)npme/(float)nnodes;
5533 lossp *= (float)npp/(float)nnodes;
5535 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5536 fprintf(fplog, "%s", buf);
5537 fprintf(stderr, "%s", buf);
5538 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5539 fprintf(fplog, "%s", buf);
5540 fprintf(stderr, "%s", buf);
5542 fprintf(fplog, "\n");
5543 fprintf(stderr, "\n");
5545 if (lossf >= DD_PERF_LOSS_WARN)
5548 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5549 " in the domain decomposition.\n", lossf*100);
5552 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5556 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5558 fprintf(fplog, "%s\n", buf);
5559 fprintf(stderr, "%s\n", buf);
5561 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5564 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5565 " had %s work to do than the PP ranks.\n"
5566 " You might want to %s the number of PME ranks\n"
5567 " or %s the cut-off and the grid spacing.\n",
5569 (lossp < 0) ? "less" : "more",
5570 (lossp < 0) ? "decrease" : "increase",
5571 (lossp < 0) ? "decrease" : "increase");
5572 fprintf(fplog, "%s\n", buf);
5573 fprintf(stderr, "%s\n", buf);
5578 static float dd_vol_min(gmx_domdec_t *dd)
5580 return dd->comm->load[0].cvol_min*dd->nnodes;
5583 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5585 return dd->comm->load[0].flags;
5588 static float dd_f_imbal(gmx_domdec_t *dd)
5590 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5593 float dd_pme_f_ratio(gmx_domdec_t *dd)
5595 /* Should only be called on the DD master rank */
5596 assert(DDMASTER(dd));
5598 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5600 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5608 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5613 flags = dd_load_flags(dd);
5617 "DD load balancing is limited by minimum cell size in dimension");
5618 for (d = 0; d < dd->ndim; d++)
5622 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5625 fprintf(fplog, "\n");
5627 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5628 if (dlbIsOn(dd->comm))
5630 fprintf(fplog, " vol min/aver %5.3f%c",
5631 dd_vol_min(dd), flags ? '!' : ' ');
5635 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5637 if (dd->comm->cycl_n[ddCyclPME])
5639 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5641 fprintf(fplog, "\n\n");
5644 static void dd_print_load_verbose(gmx_domdec_t *dd)
5646 if (dlbIsOn(dd->comm))
5648 fprintf(stderr, "vol %4.2f%c ",
5649 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5653 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5655 if (dd->comm->cycl_n[ddCyclPME])
5657 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5662 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5667 gmx_domdec_root_t *root;
5668 gmx_bool bPartOfGroup = FALSE;
5670 dim = dd->dim[dim_ind];
5671 copy_ivec(loc, loc_c);
5672 for (i = 0; i < dd->nc[dim]; i++)
5675 rank = dd_index(dd->nc, loc_c);
5676 if (rank == dd->rank)
5678 /* This process is part of the group */
5679 bPartOfGroup = TRUE;
5682 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5686 dd->comm->mpi_comm_load[dim_ind] = c_row;
5687 if (dd->comm->dlbState != edlbsOffForever)
5689 if (dd->ci[dim] == dd->master_ci[dim])
5691 /* This is the root process of this row */
5692 snew(dd->comm->root[dim_ind], 1);
5693 root = dd->comm->root[dim_ind];
5694 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5695 snew(root->old_cell_f, dd->nc[dim]+1);
5696 snew(root->bCellMin, dd->nc[dim]);
5699 snew(root->cell_f_max0, dd->nc[dim]);
5700 snew(root->cell_f_min1, dd->nc[dim]);
5701 snew(root->bound_min, dd->nc[dim]);
5702 snew(root->bound_max, dd->nc[dim]);
5704 snew(root->buf_ncd, dd->nc[dim]);
5708 /* This is not a root process, we only need to receive cell_f */
5709 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5712 if (dd->ci[dim] == dd->master_ci[dim])
5714 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5720 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5721 const gmx_hw_info_t gmx_unused *hwinfo,
5722 const gmx_hw_opt_t gmx_unused *hw_opt)
5725 int physicalnode_id_hash;
5728 MPI_Comm mpi_comm_pp_physicalnode;
5730 if (!(cr->duty & DUTY_PP) || hw_opt->gpu_opt.n_dev_use == 0)
5732 /* Only PP nodes (currently) use GPUs.
5733 * If we don't have GPUs, there are no resources to share.
5738 physicalnode_id_hash = gmx_physicalnode_id_hash();
5740 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5746 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5747 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5748 dd->rank, physicalnode_id_hash, gpu_id);
5750 /* Split the PP communicator over the physical nodes */
5751 /* TODO: See if we should store this (before), as it's also used for
5752 * for the nodecomm summution.
5754 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5755 &mpi_comm_pp_physicalnode);
5756 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5757 &dd->comm->mpi_comm_gpu_shared);
5758 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5759 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5763 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5766 /* Note that some ranks could share a GPU, while others don't */
5768 if (dd->comm->nrank_gpu_shared == 1)
5770 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5775 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5778 int dim0, dim1, i, j;
5783 fprintf(debug, "Making load communicators\n");
5786 snew(dd->comm->load, std::max(dd->ndim, 1));
5787 snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
5795 make_load_communicator(dd, 0, loc);
5799 for (i = 0; i < dd->nc[dim0]; i++)
5802 make_load_communicator(dd, 1, loc);
5808 for (i = 0; i < dd->nc[dim0]; i++)
5812 for (j = 0; j < dd->nc[dim1]; j++)
5815 make_load_communicator(dd, 2, loc);
5822 fprintf(debug, "Finished making load communicators\n");
5827 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5829 int d, dim, i, j, m;
5832 ivec dd_zp[DD_MAXIZONE];
5833 gmx_domdec_zones_t *zones;
5834 gmx_domdec_ns_ranges_t *izone;
5836 for (d = 0; d < dd->ndim; d++)
5839 copy_ivec(dd->ci, tmp);
5840 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5841 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5842 copy_ivec(dd->ci, tmp);
5843 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5844 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5847 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5850 dd->neighbor[d][1]);
5856 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5858 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5859 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5866 for (i = 0; i < nzonep; i++)
5868 copy_ivec(dd_zp3[i], dd_zp[i]);
5874 for (i = 0; i < nzonep; i++)
5876 copy_ivec(dd_zp2[i], dd_zp[i]);
5882 for (i = 0; i < nzonep; i++)
5884 copy_ivec(dd_zp1[i], dd_zp[i]);
5890 for (i = 0; i < nzonep; i++)
5892 copy_ivec(dd_zp0[i], dd_zp[i]);
5896 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5901 zones = &dd->comm->zones;
5903 for (i = 0; i < nzone; i++)
5906 clear_ivec(zones->shift[i]);
5907 for (d = 0; d < dd->ndim; d++)
5909 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5914 for (i = 0; i < nzone; i++)
5916 for (d = 0; d < DIM; d++)
5918 s[d] = dd->ci[d] - zones->shift[i][d];
5923 else if (s[d] >= dd->nc[d])
5929 zones->nizone = nzonep;
5930 for (i = 0; i < zones->nizone; i++)
5932 if (dd_zp[i][0] != i)
5934 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5936 izone = &zones->izone[i];
5937 izone->j0 = dd_zp[i][1];
5938 izone->j1 = dd_zp[i][2];
5939 for (dim = 0; dim < DIM; dim++)
5941 if (dd->nc[dim] == 1)
5943 /* All shifts should be allowed */
5944 izone->shift0[dim] = -1;
5945 izone->shift1[dim] = 1;
5950 izone->shift0[d] = 0;
5951 izone->shift1[d] = 0;
5952 for(j=izone->j0; j<izone->j1; j++) {
5953 if (dd->shift[j][d] > dd->shift[i][d])
5954 izone->shift0[d] = -1;
5955 if (dd->shift[j][d] < dd->shift[i][d])
5956 izone->shift1[d] = 1;
5962 /* Assume the shift are not more than 1 cell */
5963 izone->shift0[dim] = 1;
5964 izone->shift1[dim] = -1;
5965 for (j = izone->j0; j < izone->j1; j++)
5967 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5968 if (shift_diff < izone->shift0[dim])
5970 izone->shift0[dim] = shift_diff;
5972 if (shift_diff > izone->shift1[dim])
5974 izone->shift1[dim] = shift_diff;
5981 if (dd->comm->dlbState != edlbsOffForever)
5983 snew(dd->comm->root, dd->ndim);
5986 if (dd->comm->bRecordLoad)
5988 make_load_communicators(dd);
5992 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5998 gmx_domdec_comm_t *comm;
6005 if (comm->bCartesianPP)
6007 /* Set up cartesian communication for the particle-particle part */
6010 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
6011 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
6014 for (int i = 0; i < DIM; i++)
6018 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
6020 /* We overwrite the old communicator with the new cartesian one */
6021 cr->mpi_comm_mygroup = comm_cart;
6024 dd->mpi_comm_all = cr->mpi_comm_mygroup;
6025 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
6027 if (comm->bCartesianPP_PME)
6029 /* Since we want to use the original cartesian setup for sim,
6030 * and not the one after split, we need to make an index.
6032 snew(comm->ddindex2ddnodeid, dd->nnodes);
6033 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
6034 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
6035 /* Get the rank of the DD master,
6036 * above we made sure that the master node is a PP node.
6046 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
6048 else if (comm->bCartesianPP)
6050 if (cr->npmenodes == 0)
6052 /* The PP communicator is also
6053 * the communicator for this simulation
6055 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6057 cr->nodeid = dd->rank;
6059 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6061 /* We need to make an index to go from the coordinates
6062 * to the nodeid of this simulation.
6064 snew(comm->ddindex2simnodeid, dd->nnodes);
6065 snew(buf, dd->nnodes);
6066 if (cr->duty & DUTY_PP)
6068 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6070 /* Communicate the ddindex to simulation nodeid index */
6071 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6072 cr->mpi_comm_mysim);
6075 /* Determine the master coordinates and rank.
6076 * The DD master should be the same node as the master of this sim.
6078 for (int i = 0; i < dd->nnodes; i++)
6080 if (comm->ddindex2simnodeid[i] == 0)
6082 ddindex2xyz(dd->nc, i, dd->master_ci);
6083 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6088 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6093 /* No Cartesian communicators */
6094 /* We use the rank in dd->comm->all as DD index */
6095 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6096 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6098 clear_ivec(dd->master_ci);
6105 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6106 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6111 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6112 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6116 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6120 gmx_domdec_comm_t *comm;
6125 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6128 snew(comm->ddindex2simnodeid, dd->nnodes);
6129 snew(buf, dd->nnodes);
6130 if (cr->duty & DUTY_PP)
6132 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6134 /* Communicate the ddindex to simulation nodeid index */
6135 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6136 cr->mpi_comm_mysim);
6142 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6143 int ncg, int natoms)
6145 gmx_domdec_master_t *ma;
6150 snew(ma->ncg, dd->nnodes);
6151 snew(ma->index, dd->nnodes+1);
6153 snew(ma->nat, dd->nnodes);
6154 snew(ma->ibuf, dd->nnodes*2);
6155 snew(ma->cell_x, DIM);
6156 for (i = 0; i < DIM; i++)
6158 snew(ma->cell_x[i], dd->nc[i]+1);
6161 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6167 snew(ma->vbuf, natoms);
6173 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6174 int gmx_unused reorder)
6177 gmx_domdec_comm_t *comm;
6187 if (comm->bCartesianPP)
6189 for (i = 1; i < DIM; i++)
6191 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6193 if (bDiv[YY] || bDiv[ZZ])
6195 comm->bCartesianPP_PME = TRUE;
6196 /* If we have 2D PME decomposition, which is always in x+y,
6197 * we stack the PME only nodes in z.
6198 * Otherwise we choose the direction that provides the thinnest slab
6199 * of PME only nodes as this will have the least effect
6200 * on the PP communication.
6201 * But for the PME communication the opposite might be better.
6203 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6205 dd->nc[YY] > dd->nc[ZZ]))
6207 comm->cartpmedim = ZZ;
6211 comm->cartpmedim = YY;
6213 comm->ntot[comm->cartpmedim]
6214 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6218 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]);
6220 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6225 if (comm->bCartesianPP_PME)
6232 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]);
6235 for (i = 0; i < DIM; i++)
6239 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6241 MPI_Comm_rank(comm_cart, &rank);
6242 if (MASTERNODE(cr) && rank != 0)
6244 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6247 /* With this assigment we loose the link to the original communicator
6248 * which will usually be MPI_COMM_WORLD, unless have multisim.
6250 cr->mpi_comm_mysim = comm_cart;
6251 cr->sim_nodeid = rank;
6253 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6257 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6258 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6261 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6265 if (cr->npmenodes == 0 ||
6266 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6268 cr->duty = DUTY_PME;
6271 /* Split the sim communicator into PP and PME only nodes */
6272 MPI_Comm_split(cr->mpi_comm_mysim,
6274 dd_index(comm->ntot, dd->ci),
6275 &cr->mpi_comm_mygroup);
6279 switch (dd_node_order)
6284 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6287 case ddnoINTERLEAVE:
6288 /* Interleave the PP-only and PME-only nodes,
6289 * as on clusters with dual-core machines this will double
6290 * the communication bandwidth of the PME processes
6291 * and thus speed up the PP <-> PME and inter PME communication.
6295 fprintf(fplog, "Interleaving PP and PME ranks\n");
6297 comm->pmenodes = dd_pmenodes(cr);
6302 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6305 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6307 cr->duty = DUTY_PME;
6314 /* Split the sim communicator into PP and PME only nodes */
6315 MPI_Comm_split(cr->mpi_comm_mysim,
6318 &cr->mpi_comm_mygroup);
6319 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6325 fprintf(fplog, "This rank does only %s work.\n\n",
6326 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6330 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6333 gmx_domdec_comm_t *comm;
6339 copy_ivec(dd->nc, comm->ntot);
6341 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6342 comm->bCartesianPP_PME = FALSE;
6344 /* Reorder the nodes by default. This might change the MPI ranks.
6345 * Real reordering is only supported on very few architectures,
6346 * Blue Gene is one of them.
6348 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6350 if (cr->npmenodes > 0)
6352 /* Split the communicator into a PP and PME part */
6353 split_communicator(fplog, cr, dd_node_order, CartReorder);
6354 if (comm->bCartesianPP_PME)
6356 /* We (possibly) reordered the nodes in split_communicator,
6357 * so it is no longer required in make_pp_communicator.
6359 CartReorder = FALSE;
6364 /* All nodes do PP and PME */
6366 /* We do not require separate communicators */
6367 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6371 if (cr->duty & DUTY_PP)
6373 /* Copy or make a new PP communicator */
6374 make_pp_communicator(fplog, cr, CartReorder);
6378 receive_ddindex2simnodeid(cr);
6381 if (!(cr->duty & DUTY_PME))
6383 /* Set up the commnuication to our PME node */
6384 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6385 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6388 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6389 dd->pme_nodeid, dd->pme_receive_vir_ener);
6394 dd->pme_nodeid = -1;
6399 dd->ma = init_gmx_domdec_master_t(dd,
6401 comm->cgs_gl.index[comm->cgs_gl.nr]);
6405 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6407 real *slb_frac, tot;
6412 if (nc > 1 && size_string != NULL)
6416 fprintf(fplog, "Using static load balancing for the %s direction\n",
6421 for (i = 0; i < nc; i++)
6424 sscanf(size_string, "%20lf%n", &dbl, &n);
6427 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6436 fprintf(fplog, "Relative cell sizes:");
6438 for (i = 0; i < nc; i++)
6443 fprintf(fplog, " %5.3f", slb_frac[i]);
6448 fprintf(fplog, "\n");
6455 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6458 gmx_mtop_ilistloop_t iloop;
6462 iloop = gmx_mtop_ilistloop_init(mtop);
6463 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6465 for (ftype = 0; ftype < F_NRE; ftype++)
6467 if ((interaction_function[ftype].flags & IF_BOND) &&
6470 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6478 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6484 val = getenv(env_var);
6487 if (sscanf(val, "%20d", &nst) <= 0)
6493 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6501 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6505 fprintf(stderr, "\n%s\n", warn_string);
6509 fprintf(fplog, "\n%s\n", warn_string);
6513 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6514 t_inputrec *ir, FILE *fplog)
6516 if (ir->ePBC == epbcSCREW &&
6517 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6519 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6522 if (ir->ns_type == ensSIMPLE)
6524 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6527 if (ir->nstlist == 0)
6529 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6532 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6534 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6538 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6543 r = ddbox->box_size[XX];
6544 for (di = 0; di < dd->ndim; di++)
6547 /* Check using the initial average cell size */
6548 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6554 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6555 const char *dlb_opt, gmx_bool bRecordLoad,
6556 unsigned long Flags, t_inputrec *ir)
6563 case 'a': dlbState = edlbsOffCanTurnOn; break;
6564 case 'n': dlbState = edlbsOffForever; break;
6565 case 'y': dlbState = edlbsOn; break;
6566 default: gmx_incons("Unknown dlb_opt");
6569 if (Flags & MD_RERUN)
6571 return edlbsOffForever;
6574 if (!EI_DYNAMICS(ir->eI))
6576 if (dlbState == edlbsOn)
6578 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6579 dd_warning(cr, fplog, buf);
6582 return edlbsOffForever;
6587 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6589 return edlbsOffForever;
6592 if (Flags & MD_REPRODUCIBLE)
6596 case edlbsOffForever:
6598 case edlbsOffCanTurnOn:
6599 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6600 dlbState = edlbsOffForever;
6603 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6606 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", dlbState);
6614 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6619 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6621 /* Decomposition order z,y,x */
6624 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6626 for (dim = DIM-1; dim >= 0; dim--)
6628 if (dd->nc[dim] > 1)
6630 dd->dim[dd->ndim++] = dim;
6636 /* Decomposition order x,y,z */
6637 for (dim = 0; dim < DIM; dim++)
6639 if (dd->nc[dim] > 1)
6641 dd->dim[dd->ndim++] = dim;
6647 static gmx_domdec_comm_t *init_dd_comm()
6649 gmx_domdec_comm_t *comm;
6653 snew(comm->cggl_flag, DIM*2);
6654 snew(comm->cgcm_state, DIM*2);
6655 for (i = 0; i < DIM*2; i++)
6657 comm->cggl_flag_nalloc[i] = 0;
6658 comm->cgcm_state_nalloc[i] = 0;
6661 comm->nalloc_int = 0;
6662 comm->buf_int = NULL;
6664 vec_rvec_init(&comm->vbuf);
6666 comm->n_load_have = 0;
6667 comm->n_load_collect = 0;
6669 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6671 comm->sum_nat[i] = 0;
6675 comm->load_step = 0;
6678 clear_ivec(comm->load_lim);
6685 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6686 unsigned long Flags,
6688 real comm_distance_min, real rconstr,
6689 const char *dlb_opt, real dlb_scale,
6690 const char *sizex, const char *sizey, const char *sizez,
6691 gmx_mtop_t *mtop, t_inputrec *ir,
6692 matrix box, rvec *x,
6694 int *npme_x, int *npme_y)
6697 gmx_domdec_comm_t *comm;
6699 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6702 const real tenPercentMargin = 1.1;
6707 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6712 dd->comm = init_dd_comm();
6714 snew(comm->cggl_flag, DIM*2);
6715 snew(comm->cgcm_state, DIM*2);
6717 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6718 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6720 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6721 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6722 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6723 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6724 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6725 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6726 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6727 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6729 dd->pme_recv_f_alloc = 0;
6730 dd->pme_recv_f_buf = NULL;
6732 if (dd->bSendRecv2 && fplog)
6734 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");
6740 fprintf(fplog, "Will load balance based on FLOP count\n");
6742 if (comm->eFlop > 1)
6744 srand(1+cr->nodeid);
6746 comm->bRecordLoad = TRUE;
6750 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6754 /* Initialize to GPU share count to 0, might change later */
6755 comm->nrank_gpu_shared = 0;
6757 comm->dlbState = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6758 comm->bCheckWhetherToTurnDlbOn = TRUE;
6762 fprintf(fplog, "Dynamic load balancing: %s\n",
6763 edlbs_names[comm->dlbState]);
6765 comm->bPMELoadBalDLBLimits = FALSE;
6767 if (comm->nstSortCG)
6771 if (comm->nstSortCG == 1)
6773 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6777 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6781 snew(comm->sort, 1);
6787 fprintf(fplog, "Will not sort the charge groups\n");
6791 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6793 comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
6794 mtop->bIntermolecularInteractions);
6795 if (comm->bInterCGBondeds)
6797 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6801 comm->bInterCGMultiBody = FALSE;
6804 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6805 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6807 if (ir->rlistlong == 0)
6809 /* Set the cut-off to some very large value,
6810 * so we don't need if statements everywhere in the code.
6811 * We use sqrt, since the cut-off is squared in some places.
6813 comm->cutoff = GMX_CUTOFF_INF;
6817 comm->cutoff = ir->rlistlong;
6819 comm->cutoff_mbody = 0;
6821 comm->cellsize_limit = 0;
6822 comm->bBondComm = FALSE;
6824 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6825 * within nstlist steps. Since boundaries are allowed to displace by half
6826 * a cell size, DD cells should be at least the size of the list buffer.
6828 comm->cellsize_limit = std::max(comm->cellsize_limit,
6829 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6831 if (comm->bInterCGBondeds)
6833 if (comm_distance_min > 0)
6835 comm->cutoff_mbody = comm_distance_min;
6836 if (Flags & MD_DDBONDCOMM)
6838 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6842 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6844 r_bonded_limit = comm->cutoff_mbody;
6846 else if (ir->bPeriodicMols)
6848 /* Can not easily determine the required cut-off */
6849 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");
6850 comm->cutoff_mbody = comm->cutoff/2;
6851 r_bonded_limit = comm->cutoff_mbody;
6857 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6858 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6860 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6861 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6863 /* We use an initial margin of 10% for the minimum cell size,
6864 * except when we are just below the non-bonded cut-off.
6866 if (Flags & MD_DDBONDCOMM)
6868 if (std::max(r_2b, r_mb) > comm->cutoff)
6870 r_bonded = std::max(r_2b, r_mb);
6871 r_bonded_limit = tenPercentMargin*r_bonded;
6872 comm->bBondComm = TRUE;
6877 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6879 /* We determine cutoff_mbody later */
6883 /* No special bonded communication,
6884 * simply increase the DD cut-off.
6886 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6887 comm->cutoff_mbody = r_bonded_limit;
6888 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6894 "Minimum cell size due to bonded interactions: %.3f nm\n",
6897 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6900 if (dd->bInterCGcons && rconstr <= 0)
6902 /* There is a cell size limit due to the constraints (P-LINCS) */
6903 rconstr = constr_r_max(fplog, mtop, ir);
6907 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6909 if (rconstr > comm->cellsize_limit)
6911 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6915 else if (rconstr > 0 && fplog)
6917 /* Here we do not check for dd->bInterCGcons,
6918 * because one can also set a cell size limit for virtual sites only
6919 * and at this point we don't know yet if there are intercg v-sites.
6922 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6925 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6927 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6931 copy_ivec(nc, dd->nc);
6932 set_dd_dim(fplog, dd);
6933 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6935 if (cr->npmenodes == -1)
6939 acs = average_cellsize_min(dd, ddbox);
6940 if (acs < comm->cellsize_limit)
6944 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6946 gmx_fatal_collective(FARGS, cr, NULL,
6947 "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",
6948 acs, comm->cellsize_limit);
6953 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6955 /* We need to choose the optimal DD grid and possibly PME nodes */
6956 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6957 comm->dlbState != edlbsOffForever, dlb_scale,
6958 comm->cellsize_limit, comm->cutoff,
6959 comm->bInterCGBondeds);
6961 if (dd->nc[XX] == 0)
6963 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6964 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6965 !bC ? "-rdd" : "-rcon",
6966 comm->dlbState != edlbsOffForever ? " or -dds" : "",
6967 bC ? " or your LINCS settings" : "");
6969 gmx_fatal_collective(FARGS, cr, NULL,
6970 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6972 "Look in the log file for details on the domain decomposition",
6973 cr->nnodes-cr->npmenodes, limit, buf);
6975 set_dd_dim(fplog, dd);
6981 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6982 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6985 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6986 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6988 gmx_fatal_collective(FARGS, cr, NULL,
6989 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6990 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6992 if (cr->npmenodes > dd->nnodes)
6994 gmx_fatal_collective(FARGS, cr, NULL,
6995 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6997 if (cr->npmenodes > 0)
6999 comm->npmenodes = cr->npmenodes;
7003 comm->npmenodes = dd->nnodes;
7006 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7008 /* The following choices should match those
7009 * in comm_cost_est in domdec_setup.c.
7010 * Note that here the checks have to take into account
7011 * that the decomposition might occur in a different order than xyz
7012 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
7013 * in which case they will not match those in comm_cost_est,
7014 * but since that is mainly for testing purposes that's fine.
7016 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
7017 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
7018 getenv("GMX_PMEONEDD") == NULL)
7020 comm->npmedecompdim = 2;
7021 comm->npmenodes_x = dd->nc[XX];
7022 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
7026 /* In case nc is 1 in both x and y we could still choose to
7027 * decompose pme in y instead of x, but we use x for simplicity.
7029 comm->npmedecompdim = 1;
7030 if (dd->dim[0] == YY)
7032 comm->npmenodes_x = 1;
7033 comm->npmenodes_y = comm->npmenodes;
7037 comm->npmenodes_x = comm->npmenodes;
7038 comm->npmenodes_y = 1;
7043 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
7044 comm->npmenodes_x, comm->npmenodes_y, 1);
7049 comm->npmedecompdim = 0;
7050 comm->npmenodes_x = 0;
7051 comm->npmenodes_y = 0;
7054 /* Technically we don't need both of these,
7055 * but it simplifies code not having to recalculate it.
7057 *npme_x = comm->npmenodes_x;
7058 *npme_y = comm->npmenodes_y;
7060 snew(comm->slb_frac, DIM);
7061 if (comm->dlbState == edlbsOffForever)
7063 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7064 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7065 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7068 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7070 if (comm->bBondComm || comm->dlbState != edlbsOffForever)
7072 /* Set the bonded communication distance to halfway
7073 * the minimum and the maximum,
7074 * since the extra communication cost is nearly zero.
7076 acs = average_cellsize_min(dd, ddbox);
7077 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7078 if (comm->dlbState != edlbsOffForever)
7080 /* Check if this does not limit the scaling */
7081 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7083 if (!comm->bBondComm)
7085 /* Without bBondComm do not go beyond the n.b. cut-off */
7086 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7087 if (comm->cellsize_limit >= comm->cutoff)
7089 /* We don't loose a lot of efficieny
7090 * when increasing it to the n.b. cut-off.
7091 * It can even be slightly faster, because we need
7092 * less checks for the communication setup.
7094 comm->cutoff_mbody = comm->cutoff;
7097 /* Check if we did not end up below our original limit */
7098 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7100 if (comm->cutoff_mbody > comm->cellsize_limit)
7102 comm->cellsize_limit = comm->cutoff_mbody;
7105 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7110 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7111 "cellsize limit %f\n",
7112 comm->bBondComm, comm->cellsize_limit);
7117 check_dd_restrictions(cr, dd, ir, fplog);
7120 comm->partition_step = INT_MIN;
7123 clear_dd_cycle_counts(dd);
7128 static void set_dlb_limits(gmx_domdec_t *dd)
7133 for (d = 0; d < dd->ndim; d++)
7135 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7136 dd->comm->cellsize_min[dd->dim[d]] =
7137 dd->comm->cellsize_min_dlb[dd->dim[d]];
7142 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7145 gmx_domdec_comm_t *comm;
7155 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);
7158 cellsize_min = comm->cellsize_min[dd->dim[0]];
7159 for (d = 1; d < dd->ndim; d++)
7161 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7164 if (cellsize_min < comm->cellsize_limit*1.05)
7166 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");
7168 /* Change DLB from "auto" to "no". */
7169 comm->dlbState = edlbsOffForever;
7174 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7175 comm->dlbState = edlbsOn;
7179 /* We can set the required cell size info here,
7180 * so we do not need to communicate this.
7181 * The grid is completely uniform.
7183 for (d = 0; d < dd->ndim; d++)
7187 comm->load[d].sum_m = comm->load[d].sum;
7189 nc = dd->nc[dd->dim[d]];
7190 for (i = 0; i < nc; i++)
7192 comm->root[d]->cell_f[i] = i/(real)nc;
7195 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7196 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7199 comm->root[d]->cell_f[nc] = 1.0;
7204 static char *init_bLocalCG(gmx_mtop_t *mtop)
7209 ncg = ncg_mtop(mtop);
7210 snew(bLocalCG, ncg);
7211 for (cg = 0; cg < ncg; cg++)
7213 bLocalCG[cg] = FALSE;
7219 void dd_init_bondeds(FILE *fplog,
7220 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7222 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7224 gmx_domdec_comm_t *comm;
7226 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7230 if (comm->bBondComm)
7232 /* Communicate atoms beyond the cut-off for bonded interactions */
7235 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7237 comm->bLocalCG = init_bLocalCG(mtop);
7241 /* Only communicate atoms based on cut-off */
7242 comm->cglink = NULL;
7243 comm->bLocalCG = NULL;
7247 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7249 gmx_bool bDynLoadBal, real dlb_scale,
7252 gmx_domdec_comm_t *comm;
7267 fprintf(fplog, "The maximum number of communication pulses is:");
7268 for (d = 0; d < dd->ndim; d++)
7270 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7272 fprintf(fplog, "\n");
7273 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7274 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7275 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7276 for (d = 0; d < DIM; d++)
7280 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7287 comm->cellsize_min_dlb[d]/
7288 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7290 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7293 fprintf(fplog, "\n");
7297 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7298 fprintf(fplog, "The initial number of communication pulses is:");
7299 for (d = 0; d < dd->ndim; d++)
7301 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7303 fprintf(fplog, "\n");
7304 fprintf(fplog, "The initial domain decomposition cell size is:");
7305 for (d = 0; d < DIM; d++)
7309 fprintf(fplog, " %c %.2f nm",
7310 dim2char(d), dd->comm->cellsize_min[d]);
7313 fprintf(fplog, "\n\n");
7316 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7318 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7319 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7320 "non-bonded interactions", "", comm->cutoff);
7324 limit = dd->comm->cellsize_limit;
7328 if (dynamic_dd_box(ddbox, ir))
7330 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7332 limit = dd->comm->cellsize_min[XX];
7333 for (d = 1; d < DIM; d++)
7335 limit = std::min(limit, dd->comm->cellsize_min[d]);
7339 if (comm->bInterCGBondeds)
7341 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7342 "two-body bonded interactions", "(-rdd)",
7343 std::max(comm->cutoff, comm->cutoff_mbody));
7344 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7345 "multi-body bonded interactions", "(-rdd)",
7346 (comm->bBondComm || dlbIsOn(dd->comm)) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7350 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7351 "virtual site constructions", "(-rcon)", limit);
7353 if (dd->constraint_comm)
7355 sprintf(buf, "atoms separated by up to %d constraints",
7357 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7358 buf, "(-rcon)", limit);
7360 fprintf(fplog, "\n");
7366 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7368 const t_inputrec *ir,
7369 const gmx_ddbox_t *ddbox)
7371 gmx_domdec_comm_t *comm;
7372 int d, dim, npulse, npulse_d_max, npulse_d;
7377 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7379 /* Determine the maximum number of comm. pulses in one dimension */
7381 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7383 /* Determine the maximum required number of grid pulses */
7384 if (comm->cellsize_limit >= comm->cutoff)
7386 /* Only a single pulse is required */
7389 else if (!bNoCutOff && comm->cellsize_limit > 0)
7391 /* We round down slightly here to avoid overhead due to the latency
7392 * of extra communication calls when the cut-off
7393 * would be only slightly longer than the cell size.
7394 * Later cellsize_limit is redetermined,
7395 * so we can not miss interactions due to this rounding.
7397 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7401 /* There is no cell size limit */
7402 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7405 if (!bNoCutOff && npulse > 1)
7407 /* See if we can do with less pulses, based on dlb_scale */
7409 for (d = 0; d < dd->ndim; d++)
7412 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7413 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7414 npulse_d_max = std::max(npulse_d_max, npulse_d);
7416 npulse = std::min(npulse, npulse_d_max);
7419 /* This env var can override npulse */
7420 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7427 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7428 for (d = 0; d < dd->ndim; d++)
7430 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7431 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7432 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7433 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7434 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7436 comm->bVacDLBNoLimit = FALSE;
7440 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7441 if (!comm->bVacDLBNoLimit)
7443 comm->cellsize_limit = std::max(comm->cellsize_limit,
7444 comm->cutoff/comm->maxpulse);
7446 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7447 /* Set the minimum cell size for each DD dimension */
7448 for (d = 0; d < dd->ndim; d++)
7450 if (comm->bVacDLBNoLimit ||
7451 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7453 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7457 comm->cellsize_min_dlb[dd->dim[d]] =
7458 comm->cutoff/comm->cd[d].np_dlb;
7461 if (comm->cutoff_mbody <= 0)
7463 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7471 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7473 /* If each molecule is a single charge group
7474 * or we use domain decomposition for each periodic dimension,
7475 * we do not need to take pbc into account for the bonded interactions.
7477 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7480 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7483 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7484 t_inputrec *ir, gmx_ddbox_t *ddbox)
7486 gmx_domdec_comm_t *comm;
7492 /* Initialize the thread data.
7493 * This can not be done in init_domain_decomposition,
7494 * as the numbers of threads is determined later.
7496 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7499 snew(comm->dth, comm->nth);
7502 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7504 init_ddpme(dd, &comm->ddpme[0], 0);
7505 if (comm->npmedecompdim >= 2)
7507 init_ddpme(dd, &comm->ddpme[1], 1);
7512 comm->npmenodes = 0;
7513 if (dd->pme_nodeid >= 0)
7515 gmx_fatal_collective(FARGS, NULL, dd,
7516 "Can not have separate PME ranks without PME electrostatics");
7522 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7524 if (comm->dlbState != edlbsOffForever)
7526 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7529 print_dd_settings(fplog, dd, ir, dlbIsOn(comm), dlb_scale, ddbox);
7530 if (comm->dlbState == edlbsOffCanTurnOn)
7534 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7536 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7539 if (ir->ePBC == epbcNONE)
7541 vol_frac = 1 - 1/(double)dd->nnodes;
7546 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7550 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7552 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7554 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7557 static gmx_bool test_dd_cutoff(t_commrec *cr,
7558 t_state *state, t_inputrec *ir,
7569 set_ddbox(dd, FALSE, cr, ir, state->box,
7570 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7574 for (d = 0; d < dd->ndim; d++)
7578 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7579 if (dynamic_dd_box(&ddbox, ir))
7581 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7584 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7586 if (dd->comm->dlbState != edlbsOffForever && dim < ddbox.npbcdim &&
7587 dd->comm->cd[d].np_dlb > 0)
7589 if (np > dd->comm->cd[d].np_dlb)
7594 /* If a current local cell size is smaller than the requested
7595 * cut-off, we could still fix it, but this gets very complicated.
7596 * Without fixing here, we might actually need more checks.
7598 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7605 if (dd->comm->dlbState != edlbsOffForever)
7607 /* If DLB is not active yet, we don't need to check the grid jumps.
7608 * Actually we shouldn't, because then the grid jump data is not set.
7610 if (dlbIsOn(dd->comm) &&
7611 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7616 gmx_sumi(1, &LocallyLimited, cr);
7618 if (LocallyLimited > 0)
7627 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7630 gmx_bool bCutoffAllowed;
7632 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7636 cr->dd->comm->cutoff = cutoff_req;
7639 return bCutoffAllowed;
7642 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
7644 gmx_domdec_comm_t *comm;
7646 comm = cr->dd->comm;
7648 /* Turn on the DLB limiting (might have been on already) */
7649 comm->bPMELoadBalDLBLimits = TRUE;
7651 /* Change the cut-off limit */
7652 comm->PMELoadBal_max_cutoff = cutoff;
7656 fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
7657 comm->PMELoadBal_max_cutoff);
7661 /* Sets whether we should later check the load imbalance data, so that
7662 * we can trigger dynamic load balancing if enough imbalance has
7665 * Used after PME load balancing unlocks DLB, so that the check
7666 * whether DLB will be useful can happen immediately.
7668 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
7670 if (dd->comm->dlbState == edlbsOffCanTurnOn)
7672 dd->comm->bCheckWhetherToTurnDlbOn = bValue;
7676 /* Returns if we should check whether there has been enough load
7677 * imbalance to trigger dynamic load balancing.
7679 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
7681 const int nddp_chk_dlb = 100;
7683 if (dd->comm->dlbState != edlbsOffCanTurnOn)
7688 /* We should check whether we should use DLB directly after
7690 if (dd->comm->bCheckWhetherToTurnDlbOn)
7692 /* This flag was set when the PME load-balancing routines
7693 unlocked DLB, and should now be cleared. */
7694 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
7697 /* We should also check whether we should use DLB every 100
7698 * partitionings (we do not do this every partioning, so that we
7699 * avoid excessive communication). */
7700 if (dd->comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1)
7708 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
7710 return (dd->comm->dlbState == edlbsOn);
7713 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7715 return (dd->comm->dlbState == edlbsOffTemporarilyLocked);
7718 void dd_dlb_lock(gmx_domdec_t *dd)
7720 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7721 if (dd->comm->dlbState == edlbsOffCanTurnOn)
7723 dd->comm->dlbState = edlbsOffTemporarilyLocked;
7727 void dd_dlb_unlock(gmx_domdec_t *dd)
7729 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7730 if (dd->comm->dlbState == edlbsOffTemporarilyLocked)
7732 dd->comm->dlbState = edlbsOffCanTurnOn;
7733 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
7737 static void merge_cg_buffers(int ncell,
7738 gmx_domdec_comm_dim_t *cd, int pulse,
7740 int *index_gl, int *recv_i,
7741 rvec *cg_cm, rvec *recv_vr,
7743 cginfo_mb_t *cginfo_mb, int *cginfo)
7745 gmx_domdec_ind_t *ind, *ind_p;
7746 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7747 int shift, shift_at;
7749 ind = &cd->ind[pulse];
7751 /* First correct the already stored data */
7752 shift = ind->nrecv[ncell];
7753 for (cell = ncell-1; cell >= 0; cell--)
7755 shift -= ind->nrecv[cell];
7758 /* Move the cg's present from previous grid pulses */
7759 cg0 = ncg_cell[ncell+cell];
7760 cg1 = ncg_cell[ncell+cell+1];
7761 cgindex[cg1+shift] = cgindex[cg1];
7762 for (cg = cg1-1; cg >= cg0; cg--)
7764 index_gl[cg+shift] = index_gl[cg];
7765 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7766 cgindex[cg+shift] = cgindex[cg];
7767 cginfo[cg+shift] = cginfo[cg];
7769 /* Correct the already stored send indices for the shift */
7770 for (p = 1; p <= pulse; p++)
7772 ind_p = &cd->ind[p];
7774 for (c = 0; c < cell; c++)
7776 cg0 += ind_p->nsend[c];
7778 cg1 = cg0 + ind_p->nsend[cell];
7779 for (cg = cg0; cg < cg1; cg++)
7781 ind_p->index[cg] += shift;
7787 /* Merge in the communicated buffers */
7791 for (cell = 0; cell < ncell; cell++)
7793 cg1 = ncg_cell[ncell+cell+1] + shift;
7796 /* Correct the old cg indices */
7797 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7799 cgindex[cg+1] += shift_at;
7802 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7804 /* Copy this charge group from the buffer */
7805 index_gl[cg1] = recv_i[cg0];
7806 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7807 /* Add it to the cgindex */
7808 cg_gl = index_gl[cg1];
7809 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7810 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7811 cgindex[cg1+1] = cgindex[cg1] + nat;
7816 shift += ind->nrecv[cell];
7817 ncg_cell[ncell+cell+1] = cg1;
7821 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7822 int nzone, int cg0, const int *cgindex)
7826 /* Store the atom block boundaries for easy copying of communication buffers
7829 for (zone = 0; zone < nzone; zone++)
7831 for (p = 0; p < cd->np; p++)
7833 cd->ind[p].cell2at0[zone] = cgindex[cg];
7834 cg += cd->ind[p].nrecv[zone];
7835 cd->ind[p].cell2at1[zone] = cgindex[cg];
7840 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7846 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7848 if (!bLocalCG[link->a[i]])
7857 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7859 real c[DIM][4]; /* the corners for the non-bonded communication */
7860 real cr0; /* corner for rounding */
7861 real cr1[4]; /* corners for rounding */
7862 real bc[DIM]; /* corners for bounded communication */
7863 real bcr1; /* corner for rounding for bonded communication */
7866 /* Determine the corners of the domain(s) we are communicating with */
7868 set_dd_corners(const gmx_domdec_t *dd,
7869 int dim0, int dim1, int dim2,
7873 const gmx_domdec_comm_t *comm;
7874 const gmx_domdec_zones_t *zones;
7879 zones = &comm->zones;
7881 /* Keep the compiler happy */
7885 /* The first dimension is equal for all cells */
7886 c->c[0][0] = comm->cell_x0[dim0];
7889 c->bc[0] = c->c[0][0];
7894 /* This cell row is only seen from the first row */
7895 c->c[1][0] = comm->cell_x0[dim1];
7896 /* All rows can see this row */
7897 c->c[1][1] = comm->cell_x0[dim1];
7898 if (dlbIsOn(dd->comm))
7900 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7903 /* For the multi-body distance we need the maximum */
7904 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7907 /* Set the upper-right corner for rounding */
7908 c->cr0 = comm->cell_x1[dim0];
7913 for (j = 0; j < 4; j++)
7915 c->c[2][j] = comm->cell_x0[dim2];
7917 if (dlbIsOn(dd->comm))
7919 /* Use the maximum of the i-cells that see a j-cell */
7920 for (i = 0; i < zones->nizone; i++)
7922 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7927 std::max(c->c[2][j-4],
7928 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7934 /* For the multi-body distance we need the maximum */
7935 c->bc[2] = comm->cell_x0[dim2];
7936 for (i = 0; i < 2; i++)
7938 for (j = 0; j < 2; j++)
7940 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7946 /* Set the upper-right corner for rounding */
7947 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7948 * Only cell (0,0,0) can see cell 7 (1,1,1)
7950 c->cr1[0] = comm->cell_x1[dim1];
7951 c->cr1[3] = comm->cell_x1[dim1];
7952 if (dlbIsOn(dd->comm))
7954 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7957 /* For the multi-body distance we need the maximum */
7958 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7965 /* Determine which cg's we need to send in this pulse from this zone */
7967 get_zone_pulse_cgs(gmx_domdec_t *dd,
7968 int zonei, int zone,
7970 const int *index_gl,
7972 int dim, int dim_ind,
7973 int dim0, int dim1, int dim2,
7974 real r_comm2, real r_bcomm2,
7978 real skew_fac2_d, real skew_fac_01,
7979 rvec *v_d, rvec *v_0, rvec *v_1,
7980 const dd_corners_t *c,
7982 gmx_bool bDistBonded,
7988 gmx_domdec_ind_t *ind,
7989 int **ibuf, int *ibuf_nalloc,
7995 gmx_domdec_comm_t *comm;
7997 gmx_bool bDistMB_pulse;
7999 real r2, rb2, r, tric_sh;
8002 int nsend_z, nsend, nat;
8006 bScrew = (dd->bScrewPBC && dim == XX);
8008 bDistMB_pulse = (bDistMB && bDistBonded);
8014 for (cg = cg0; cg < cg1; cg++)
8018 if (tric_dist[dim_ind] == 0)
8020 /* Rectangular direction, easy */
8021 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
8028 r = cg_cm[cg][dim] - c->bc[dim_ind];
8034 /* Rounding gives at most a 16% reduction
8035 * in communicated atoms
8037 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8039 r = cg_cm[cg][dim0] - c->cr0;
8040 /* This is the first dimension, so always r >= 0 */
8047 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8049 r = cg_cm[cg][dim1] - c->cr1[zone];
8056 r = cg_cm[cg][dim1] - c->bcr1;
8066 /* Triclinic direction, more complicated */
8069 /* Rounding, conservative as the skew_fac multiplication
8070 * will slightly underestimate the distance.
8072 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8074 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
8075 for (i = dim0+1; i < DIM; i++)
8077 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
8079 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
8082 rb[dim0] = rn[dim0];
8085 /* Take care that the cell planes along dim0 might not
8086 * be orthogonal to those along dim1 and dim2.
8088 for (i = 1; i <= dim_ind; i++)
8091 if (normal[dim0][dimd] > 0)
8093 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
8096 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
8101 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8103 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
8105 for (i = dim1+1; i < DIM; i++)
8107 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
8109 rn[dim1] += tric_sh;
8112 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
8113 /* Take care of coupling of the distances
8114 * to the planes along dim0 and dim1 through dim2.
8116 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
8117 /* Take care that the cell planes along dim1
8118 * might not be orthogonal to that along dim2.
8120 if (normal[dim1][dim2] > 0)
8122 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8128 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8131 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8132 /* Take care of coupling of the distances
8133 * to the planes along dim0 and dim1 through dim2.
8135 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8136 /* Take care that the cell planes along dim1
8137 * might not be orthogonal to that along dim2.
8139 if (normal[dim1][dim2] > 0)
8141 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8146 /* The distance along the communication direction */
8147 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8149 for (i = dim+1; i < DIM; i++)
8151 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8156 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8157 /* Take care of coupling of the distances
8158 * to the planes along dim0 and dim1 through dim2.
8160 if (dim_ind == 1 && zonei == 1)
8162 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8168 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8171 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8172 /* Take care of coupling of the distances
8173 * to the planes along dim0 and dim1 through dim2.
8175 if (dim_ind == 1 && zonei == 1)
8177 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8185 ((bDistMB && rb2 < r_bcomm2) ||
8186 (bDist2B && r2 < r_bcomm2)) &&
8188 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8189 missing_link(comm->cglink, index_gl[cg],
8192 /* Make an index to the local charge groups */
8193 if (nsend+1 > ind->nalloc)
8195 ind->nalloc = over_alloc_large(nsend+1);
8196 srenew(ind->index, ind->nalloc);
8198 if (nsend+1 > *ibuf_nalloc)
8200 *ibuf_nalloc = over_alloc_large(nsend+1);
8201 srenew(*ibuf, *ibuf_nalloc);
8203 ind->index[nsend] = cg;
8204 (*ibuf)[nsend] = index_gl[cg];
8206 vec_rvec_check_alloc(vbuf, nsend+1);
8208 if (dd->ci[dim] == 0)
8210 /* Correct cg_cm for pbc */
8211 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8214 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8215 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8220 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8223 nat += cgindex[cg+1] - cgindex[cg];
8229 *nsend_z_ptr = nsend_z;
8232 static void setup_dd_communication(gmx_domdec_t *dd,
8233 matrix box, gmx_ddbox_t *ddbox,
8234 t_forcerec *fr, t_state *state, rvec **f)
8236 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8237 int nzone, nzone_send, zone, zonei, cg0, cg1;
8238 int c, i, cg, cg_gl, nrcg;
8239 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8240 gmx_domdec_comm_t *comm;
8241 gmx_domdec_zones_t *zones;
8242 gmx_domdec_comm_dim_t *cd;
8243 gmx_domdec_ind_t *ind;
8244 cginfo_mb_t *cginfo_mb;
8245 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8246 real r_comm2, r_bcomm2;
8247 dd_corners_t corners;
8249 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8250 real skew_fac2_d, skew_fac_01;
8257 fprintf(debug, "Setting up DD communication\n");
8262 switch (fr->cutoff_scheme)
8271 gmx_incons("unimplemented");
8275 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8277 /* Check if we need to use triclinic distances */
8278 tric_dist[dim_ind] = 0;
8279 for (i = 0; i <= dim_ind; i++)
8281 if (ddbox->tric_dir[dd->dim[i]])
8283 tric_dist[dim_ind] = 1;
8288 bBondComm = comm->bBondComm;
8290 /* Do we need to determine extra distances for multi-body bondeds? */
8291 bDistMB = (comm->bInterCGMultiBody && dlbIsOn(dd->comm) && dd->ndim > 1);
8293 /* Do we need to determine extra distances for only two-body bondeds? */
8294 bDist2B = (bBondComm && !bDistMB);
8296 r_comm2 = sqr(comm->cutoff);
8297 r_bcomm2 = sqr(comm->cutoff_mbody);
8301 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8304 zones = &comm->zones;
8307 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8308 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8310 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8312 /* Triclinic stuff */
8313 normal = ddbox->normal;
8317 v_0 = ddbox->v[dim0];
8318 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8320 /* Determine the coupling coefficient for the distances
8321 * to the cell planes along dim0 and dim1 through dim2.
8322 * This is required for correct rounding.
8325 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8328 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8334 v_1 = ddbox->v[dim1];
8337 zone_cg_range = zones->cg_range;
8338 index_gl = dd->index_gl;
8339 cgindex = dd->cgindex;
8340 cginfo_mb = fr->cginfo_mb;
8342 zone_cg_range[0] = 0;
8343 zone_cg_range[1] = dd->ncg_home;
8344 comm->zone_ncg1[0] = dd->ncg_home;
8345 pos_cg = dd->ncg_home;
8347 nat_tot = dd->nat_home;
8349 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8351 dim = dd->dim[dim_ind];
8352 cd = &comm->cd[dim_ind];
8354 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8356 /* No pbc in this dimension, the first node should not comm. */
8364 v_d = ddbox->v[dim];
8365 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8367 cd->bInPlace = TRUE;
8368 for (p = 0; p < cd->np; p++)
8370 /* Only atoms communicated in the first pulse are used
8371 * for multi-body bonded interactions or for bBondComm.
8373 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8378 for (zone = 0; zone < nzone_send; zone++)
8380 if (tric_dist[dim_ind] && dim_ind > 0)
8382 /* Determine slightly more optimized skew_fac's
8384 * This reduces the number of communicated atoms
8385 * by about 10% for 3D DD of rhombic dodecahedra.
8387 for (dimd = 0; dimd < dim; dimd++)
8389 sf2_round[dimd] = 1;
8390 if (ddbox->tric_dir[dimd])
8392 for (i = dd->dim[dimd]+1; i < DIM; i++)
8394 /* If we are shifted in dimension i
8395 * and the cell plane is tilted forward
8396 * in dimension i, skip this coupling.
8398 if (!(zones->shift[nzone+zone][i] &&
8399 ddbox->v[dimd][i][dimd] >= 0))
8402 sqr(ddbox->v[dimd][i][dimd]);
8405 sf2_round[dimd] = 1/sf2_round[dimd];
8410 zonei = zone_perm[dim_ind][zone];
8413 /* Here we permutate the zones to obtain a convenient order
8414 * for neighbor searching
8416 cg0 = zone_cg_range[zonei];
8417 cg1 = zone_cg_range[zonei+1];
8421 /* Look only at the cg's received in the previous grid pulse
8423 cg1 = zone_cg_range[nzone+zone+1];
8424 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8427 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8428 for (th = 0; th < comm->nth; th++)
8430 gmx_domdec_ind_t *ind_p;
8431 int **ibuf_p, *ibuf_nalloc_p;
8433 int *nsend_p, *nat_p;
8439 /* Thread 0 writes in the comm buffers */
8441 ibuf_p = &comm->buf_int;
8442 ibuf_nalloc_p = &comm->nalloc_int;
8443 vbuf_p = &comm->vbuf;
8446 nsend_zone_p = &ind->nsend[zone];
8450 /* Other threads write into temp buffers */
8451 ind_p = &comm->dth[th].ind;
8452 ibuf_p = &comm->dth[th].ibuf;
8453 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8454 vbuf_p = &comm->dth[th].vbuf;
8455 nsend_p = &comm->dth[th].nsend;
8456 nat_p = &comm->dth[th].nat;
8457 nsend_zone_p = &comm->dth[th].nsend_zone;
8459 comm->dth[th].nsend = 0;
8460 comm->dth[th].nat = 0;
8461 comm->dth[th].nsend_zone = 0;
8471 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8472 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8475 /* Get the cg's for this pulse in this zone */
8476 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8478 dim, dim_ind, dim0, dim1, dim2,
8481 normal, skew_fac2_d, skew_fac_01,
8482 v_d, v_0, v_1, &corners, sf2_round,
8483 bDistBonded, bBondComm,
8487 ibuf_p, ibuf_nalloc_p,
8493 /* Append data of threads>=1 to the communication buffers */
8494 for (th = 1; th < comm->nth; th++)
8496 dd_comm_setup_work_t *dth;
8499 dth = &comm->dth[th];
8501 ns1 = nsend + dth->nsend_zone;
8502 if (ns1 > ind->nalloc)
8504 ind->nalloc = over_alloc_dd(ns1);
8505 srenew(ind->index, ind->nalloc);
8507 if (ns1 > comm->nalloc_int)
8509 comm->nalloc_int = over_alloc_dd(ns1);
8510 srenew(comm->buf_int, comm->nalloc_int);
8512 if (ns1 > comm->vbuf.nalloc)
8514 comm->vbuf.nalloc = over_alloc_dd(ns1);
8515 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8518 for (i = 0; i < dth->nsend_zone; i++)
8520 ind->index[nsend] = dth->ind.index[i];
8521 comm->buf_int[nsend] = dth->ibuf[i];
8522 copy_rvec(dth->vbuf.v[i],
8523 comm->vbuf.v[nsend]);
8527 ind->nsend[zone] += dth->nsend_zone;
8530 /* Clear the counts in case we do not have pbc */
8531 for (zone = nzone_send; zone < nzone; zone++)
8533 ind->nsend[zone] = 0;
8535 ind->nsend[nzone] = nsend;
8536 ind->nsend[nzone+1] = nat;
8537 /* Communicate the number of cg's and atoms to receive */
8538 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8539 ind->nsend, nzone+2,
8540 ind->nrecv, nzone+2);
8542 /* The rvec buffer is also required for atom buffers of size nsend
8543 * in dd_move_x and dd_move_f.
8545 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8549 /* We can receive in place if only the last zone is not empty */
8550 for (zone = 0; zone < nzone-1; zone++)
8552 if (ind->nrecv[zone] > 0)
8554 cd->bInPlace = FALSE;
8559 /* The int buffer is only required here for the cg indices */
8560 if (ind->nrecv[nzone] > comm->nalloc_int2)
8562 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8563 srenew(comm->buf_int2, comm->nalloc_int2);
8565 /* The rvec buffer is also required for atom buffers
8566 * of size nrecv in dd_move_x and dd_move_f.
8568 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8569 vec_rvec_check_alloc(&comm->vbuf2, i);
8573 /* Make space for the global cg indices */
8574 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8575 || dd->cg_nalloc == 0)
8577 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8578 srenew(index_gl, dd->cg_nalloc);
8579 srenew(cgindex, dd->cg_nalloc+1);
8581 /* Communicate the global cg indices */
8584 recv_i = index_gl + pos_cg;
8588 recv_i = comm->buf_int2;
8590 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8591 comm->buf_int, nsend,
8592 recv_i, ind->nrecv[nzone]);
8594 /* Make space for cg_cm */
8595 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8596 if (fr->cutoff_scheme == ecutsGROUP)
8604 /* Communicate cg_cm */
8607 recv_vr = cg_cm + pos_cg;
8611 recv_vr = comm->vbuf2.v;
8613 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8614 comm->vbuf.v, nsend,
8615 recv_vr, ind->nrecv[nzone]);
8617 /* Make the charge group index */
8620 zone = (p == 0 ? 0 : nzone - 1);
8621 while (zone < nzone)
8623 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8625 cg_gl = index_gl[pos_cg];
8626 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8627 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8628 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8631 /* Update the charge group presence,
8632 * so we can use it in the next pass of the loop.
8634 comm->bLocalCG[cg_gl] = TRUE;
8640 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8643 zone_cg_range[nzone+zone] = pos_cg;
8648 /* This part of the code is never executed with bBondComm. */
8649 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8650 index_gl, recv_i, cg_cm, recv_vr,
8651 cgindex, fr->cginfo_mb, fr->cginfo);
8652 pos_cg += ind->nrecv[nzone];
8654 nat_tot += ind->nrecv[nzone+1];
8658 /* Store the atom block for easy copying of communication buffers */
8659 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8663 dd->index_gl = index_gl;
8664 dd->cgindex = cgindex;
8666 dd->ncg_tot = zone_cg_range[zones->n];
8667 dd->nat_tot = nat_tot;
8668 comm->nat[ddnatHOME] = dd->nat_home;
8669 for (i = ddnatZONE; i < ddnatNR; i++)
8671 comm->nat[i] = dd->nat_tot;
8676 /* We don't need to update cginfo, since that was alrady done above.
8677 * So we pass NULL for the forcerec.
8679 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8680 NULL, comm->bLocalCG);
8685 fprintf(debug, "Finished setting up DD communication, zones:");
8686 for (c = 0; c < zones->n; c++)
8688 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8690 fprintf(debug, "\n");
8694 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8698 for (c = 0; c < zones->nizone; c++)
8700 zones->izone[c].cg1 = zones->cg_range[c+1];
8701 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8702 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8706 static void set_zones_size(gmx_domdec_t *dd,
8707 matrix box, const gmx_ddbox_t *ddbox,
8708 int zone_start, int zone_end)
8710 gmx_domdec_comm_t *comm;
8711 gmx_domdec_zones_t *zones;
8720 zones = &comm->zones;
8722 /* Do we need to determine extra distances for multi-body bondeds? */
8723 bDistMB = (comm->bInterCGMultiBody && dlbIsOn(dd->comm) && dd->ndim > 1);
8725 for (z = zone_start; z < zone_end; z++)
8727 /* Copy cell limits to zone limits.
8728 * Valid for non-DD dims and non-shifted dims.
8730 copy_rvec(comm->cell_x0, zones->size[z].x0);
8731 copy_rvec(comm->cell_x1, zones->size[z].x1);
8734 for (d = 0; d < dd->ndim; d++)
8738 for (z = 0; z < zones->n; z++)
8740 /* With a staggered grid we have different sizes
8741 * for non-shifted dimensions.
8743 if (dlbIsOn(dd->comm) && zones->shift[z][dim] == 0)
8747 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8748 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8752 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8753 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8759 rcmbs = comm->cutoff_mbody;
8760 if (ddbox->tric_dir[dim])
8762 rcs /= ddbox->skew_fac[dim];
8763 rcmbs /= ddbox->skew_fac[dim];
8766 /* Set the lower limit for the shifted zone dimensions */
8767 for (z = zone_start; z < zone_end; z++)
8769 if (zones->shift[z][dim] > 0)
8772 if (!dlbIsOn(dd->comm) || d == 0)
8774 zones->size[z].x0[dim] = comm->cell_x1[dim];
8775 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8779 /* Here we take the lower limit of the zone from
8780 * the lowest domain of the zone below.
8784 zones->size[z].x0[dim] =
8785 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8791 zones->size[z].x0[dim] =
8792 zones->size[zone_perm[2][z-4]].x0[dim];
8796 zones->size[z].x0[dim] =
8797 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8800 /* A temporary limit, is updated below */
8801 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8805 for (zi = 0; zi < zones->nizone; zi++)
8807 if (zones->shift[zi][dim] == 0)
8809 /* This takes the whole zone into account.
8810 * With multiple pulses this will lead
8811 * to a larger zone then strictly necessary.
8813 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8814 zones->size[zi].x1[dim]+rcmbs);
8822 /* Loop over the i-zones to set the upper limit of each
8825 for (zi = 0; zi < zones->nizone; zi++)
8827 if (zones->shift[zi][dim] == 0)
8829 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8831 if (zones->shift[z][dim] > 0)
8833 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8834 zones->size[zi].x1[dim]+rcs);
8841 for (z = zone_start; z < zone_end; z++)
8843 /* Initialization only required to keep the compiler happy */
8844 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8847 /* To determine the bounding box for a zone we need to find
8848 * the extreme corners of 4, 2 or 1 corners.
8850 nc = 1 << (ddbox->nboundeddim - 1);
8852 for (c = 0; c < nc; c++)
8854 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8858 corner[YY] = zones->size[z].x0[YY];
8862 corner[YY] = zones->size[z].x1[YY];
8866 corner[ZZ] = zones->size[z].x0[ZZ];
8870 corner[ZZ] = zones->size[z].x1[ZZ];
8872 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8873 box[ZZ][1 - dd->dim[0]] != 0)
8875 /* With 1D domain decomposition the cg's are not in
8876 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8877 * Shift the corner of the z-vector back to along the box
8878 * vector of dimension d, so it will later end up at 0 along d.
8879 * This can affect the location of this corner along dd->dim[0]
8880 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8882 int d = 1 - dd->dim[0];
8884 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8886 /* Apply the triclinic couplings */
8887 assert(ddbox->npbcdim <= DIM);
8888 for (i = YY; i < ddbox->npbcdim; i++)
8890 for (j = XX; j < i; j++)
8892 corner[j] += corner[i]*box[i][j]/box[i][i];
8897 copy_rvec(corner, corner_min);
8898 copy_rvec(corner, corner_max);
8902 for (i = 0; i < DIM; i++)
8904 corner_min[i] = std::min(corner_min[i], corner[i]);
8905 corner_max[i] = std::max(corner_max[i], corner[i]);
8909 /* Copy the extreme cornes without offset along x */
8910 for (i = 0; i < DIM; i++)
8912 zones->size[z].bb_x0[i] = corner_min[i];
8913 zones->size[z].bb_x1[i] = corner_max[i];
8915 /* Add the offset along x */
8916 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8917 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8920 if (zone_start == 0)
8923 for (dim = 0; dim < DIM; dim++)
8925 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8927 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8932 for (z = zone_start; z < zone_end; z++)
8934 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8936 zones->size[z].x0[XX], zones->size[z].x1[XX],
8937 zones->size[z].x0[YY], zones->size[z].x1[YY],
8938 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8939 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8941 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8942 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8943 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8948 static int comp_cgsort(const void *a, const void *b)
8952 gmx_cgsort_t *cga, *cgb;
8953 cga = (gmx_cgsort_t *)a;
8954 cgb = (gmx_cgsort_t *)b;
8956 comp = cga->nsc - cgb->nsc;
8959 comp = cga->ind_gl - cgb->ind_gl;
8965 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8970 /* Order the data */
8971 for (i = 0; i < n; i++)
8973 buf[i] = a[sort[i].ind];
8976 /* Copy back to the original array */
8977 for (i = 0; i < n; i++)
8983 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8988 /* Order the data */
8989 for (i = 0; i < n; i++)
8991 copy_rvec(v[sort[i].ind], buf[i]);
8994 /* Copy back to the original array */
8995 for (i = 0; i < n; i++)
8997 copy_rvec(buf[i], v[i]);
9001 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
9004 int a, atot, cg, cg0, cg1, i;
9006 if (cgindex == NULL)
9008 /* Avoid the useless loop of the atoms within a cg */
9009 order_vec_cg(ncg, sort, v, buf);
9014 /* Order the data */
9016 for (cg = 0; cg < ncg; cg++)
9018 cg0 = cgindex[sort[cg].ind];
9019 cg1 = cgindex[sort[cg].ind+1];
9020 for (i = cg0; i < cg1; i++)
9022 copy_rvec(v[i], buf[a]);
9028 /* Copy back to the original array */
9029 for (a = 0; a < atot; a++)
9031 copy_rvec(buf[a], v[a]);
9035 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
9036 int nsort_new, gmx_cgsort_t *sort_new,
9037 gmx_cgsort_t *sort1)
9041 /* The new indices are not very ordered, so we qsort them */
9042 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
9044 /* sort2 is already ordered, so now we can merge the two arrays */
9048 while (i2 < nsort2 || i_new < nsort_new)
9052 sort1[i1++] = sort_new[i_new++];
9054 else if (i_new == nsort_new)
9056 sort1[i1++] = sort2[i2++];
9058 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
9059 (sort2[i2].nsc == sort_new[i_new].nsc &&
9060 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
9062 sort1[i1++] = sort2[i2++];
9066 sort1[i1++] = sort_new[i_new++];
9071 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
9073 gmx_domdec_sort_t *sort;
9074 gmx_cgsort_t *cgsort, *sort_i;
9075 int ncg_new, nsort2, nsort_new, i, *a, moved;
9077 sort = dd->comm->sort;
9079 a = fr->ns.grid->cell_index;
9081 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
9083 if (ncg_home_old >= 0)
9085 /* The charge groups that remained in the same ns grid cell
9086 * are completely ordered. So we can sort efficiently by sorting
9087 * the charge groups that did move into the stationary list.
9092 for (i = 0; i < dd->ncg_home; i++)
9094 /* Check if this cg did not move to another node */
9097 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
9099 /* This cg is new on this node or moved ns grid cell */
9100 if (nsort_new >= sort->sort_new_nalloc)
9102 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
9103 srenew(sort->sort_new, sort->sort_new_nalloc);
9105 sort_i = &(sort->sort_new[nsort_new++]);
9109 /* This cg did not move */
9110 sort_i = &(sort->sort2[nsort2++]);
9112 /* Sort on the ns grid cell indices
9113 * and the global topology index.
9114 * index_gl is irrelevant with cell ns,
9115 * but we set it here anyhow to avoid a conditional.
9118 sort_i->ind_gl = dd->index_gl[i];
9125 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9128 /* Sort efficiently */
9129 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9134 cgsort = sort->sort;
9136 for (i = 0; i < dd->ncg_home; i++)
9138 /* Sort on the ns grid cell indices
9139 * and the global topology index
9141 cgsort[i].nsc = a[i];
9142 cgsort[i].ind_gl = dd->index_gl[i];
9144 if (cgsort[i].nsc < moved)
9151 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9153 /* Determine the order of the charge groups using qsort */
9154 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9160 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9163 int ncg_new, i, *a, na;
9165 sort = dd->comm->sort->sort;
9167 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9170 for (i = 0; i < na; i++)
9174 sort[ncg_new].ind = a[i];
9182 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9185 gmx_domdec_sort_t *sort;
9186 gmx_cgsort_t *cgsort;
9188 int ncg_new, i, *ibuf, cgsize;
9191 sort = dd->comm->sort;
9193 if (dd->ncg_home > sort->sort_nalloc)
9195 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9196 srenew(sort->sort, sort->sort_nalloc);
9197 srenew(sort->sort2, sort->sort_nalloc);
9199 cgsort = sort->sort;
9201 switch (fr->cutoff_scheme)
9204 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9207 ncg_new = dd_sort_order_nbnxn(dd, fr);
9210 gmx_incons("unimplemented");
9214 /* We alloc with the old size, since cgindex is still old */
9215 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9216 vbuf = dd->comm->vbuf.v;
9220 cgindex = dd->cgindex;
9227 /* Remove the charge groups which are no longer at home here */
9228 dd->ncg_home = ncg_new;
9231 fprintf(debug, "Set the new home charge group count to %d\n",
9235 /* Reorder the state */
9236 for (i = 0; i < estNR; i++)
9238 if (EST_DISTR(i) && (state->flags & (1<<i)))
9243 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9246 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9249 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9252 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9256 case estDISRE_INITF:
9257 case estDISRE_RM3TAV:
9258 case estORIRE_INITF:
9260 /* No ordering required */
9263 gmx_incons("Unknown state entry encountered in dd_sort_state");
9268 if (fr->cutoff_scheme == ecutsGROUP)
9271 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9274 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9276 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9277 srenew(sort->ibuf, sort->ibuf_nalloc);
9280 /* Reorder the global cg index */
9281 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9282 /* Reorder the cginfo */
9283 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9284 /* Rebuild the local cg index */
9288 for (i = 0; i < dd->ncg_home; i++)
9290 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9291 ibuf[i+1] = ibuf[i] + cgsize;
9293 for (i = 0; i < dd->ncg_home+1; i++)
9295 dd->cgindex[i] = ibuf[i];
9300 for (i = 0; i < dd->ncg_home+1; i++)
9305 /* Set the home atom number */
9306 dd->nat_home = dd->cgindex[dd->ncg_home];
9308 if (fr->cutoff_scheme == ecutsVERLET)
9310 /* The atoms are now exactly in grid order, update the grid order */
9311 nbnxn_set_atomorder(fr->nbv->nbs);
9315 /* Copy the sorted ns cell indices back to the ns grid struct */
9316 for (i = 0; i < dd->ncg_home; i++)
9318 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9320 fr->ns.grid->nr = dd->ncg_home;
9324 static void add_dd_statistics(gmx_domdec_t *dd)
9326 gmx_domdec_comm_t *comm;
9331 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9333 comm->sum_nat[ddnat-ddnatZONE] +=
9334 comm->nat[ddnat] - comm->nat[ddnat-1];
9339 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9341 gmx_domdec_comm_t *comm;
9346 /* Reset all the statistics and counters for total run counting */
9347 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9349 comm->sum_nat[ddnat-ddnatZONE] = 0;
9353 comm->load_step = 0;
9356 clear_ivec(comm->load_lim);
9361 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9363 gmx_domdec_comm_t *comm;
9367 comm = cr->dd->comm;
9369 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9376 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");
9378 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9380 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9385 " av. #atoms communicated per step for force: %d x %.1f\n",
9389 if (cr->dd->vsite_comm)
9392 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9393 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9398 if (cr->dd->constraint_comm)
9401 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9402 1 + ir->nLincsIter, av);
9406 gmx_incons(" Unknown type for DD statistics");
9409 fprintf(fplog, "\n");
9411 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9413 print_dd_load_av(fplog, cr->dd);
9417 void dd_partition_system(FILE *fplog,
9420 gmx_bool bMasterState,
9422 t_state *state_global,
9423 gmx_mtop_t *top_global,
9425 t_state *state_local,
9428 gmx_localtop_t *top_local,
9431 gmx_shellfc_t shellfc,
9432 gmx_constr_t constr,
9434 gmx_wallcycle_t wcycle,
9438 gmx_domdec_comm_t *comm;
9439 gmx_ddbox_t ddbox = {0};
9441 gmx_int64_t step_pcoupl;
9442 rvec cell_ns_x0, cell_ns_x1;
9443 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9444 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bTurnOnDLB, bLogLoad;
9445 gmx_bool bRedist, bSortCG, bResortAll;
9446 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9450 wallcycle_start(wcycle, ewcDOMDEC);
9455 bBoxChanged = (bMasterState || DEFORM(*ir));
9456 if (ir->epc != epcNO)
9458 /* With nstpcouple > 1 pressure coupling happens.
9459 * one step after calculating the pressure.
9460 * Box scaling happens at the end of the MD step,
9461 * after the DD partitioning.
9462 * We therefore have to do DLB in the first partitioning
9463 * after an MD step where P-coupling occured.
9464 * We need to determine the last step in which p-coupling occurred.
9465 * MRS -- need to validate this for vv?
9470 step_pcoupl = step - 1;
9474 step_pcoupl = ((step - 1)/n)*n + 1;
9476 if (step_pcoupl >= comm->partition_step)
9482 bNStGlobalComm = (step % nstglobalcomm == 0);
9490 /* Should we do dynamic load balacing this step?
9491 * Since it requires (possibly expensive) global communication,
9492 * we might want to do DLB less frequently.
9494 if (bBoxChanged || ir->epc != epcNO)
9496 bDoDLB = bBoxChanged;
9500 bDoDLB = bNStGlobalComm;
9504 /* Check if we have recorded loads on the nodes */
9505 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9507 bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
9509 /* Print load every nstlog, first and last step to the log file */
9510 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9511 comm->n_load_collect == 0 ||
9513 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9515 /* Avoid extra communication due to verbose screen output
9516 * when nstglobalcomm is set.
9518 if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
9519 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9521 get_load_distribution(dd, wcycle);
9526 dd_print_load(fplog, dd, step-1);
9530 dd_print_load_verbose(dd);
9533 comm->n_load_collect++;
9535 if (bCheckWhetherToTurnDlbOn)
9537 /* Since the timings are node dependent, the master decides */
9540 /* Here we check if the max PME rank load is more than 0.98
9541 * the max PP force load. If so, PP DLB will not help,
9542 * since we are (almost) limited by PME. Furthermore,
9543 * DLB will cause a significant extra x/f redistribution
9544 * cost on the PME ranks, which will then surely result
9545 * in lower total performance.
9546 * This check might be fragile, since one measurement
9547 * below 0.98 (although only done once every 100 DD part.)
9548 * could turn on DLB for the rest of the run.
9550 if (cr->npmenodes > 0 &&
9551 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9558 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9562 fprintf(debug, "step %s, imb loss %f\n",
9563 gmx_step_str(step, sbuf),
9564 dd_force_imb_perf_loss(dd));
9567 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9570 turn_on_dlb(fplog, cr, step);
9575 comm->n_load_have++;
9578 cgs_gl = &comm->cgs_gl;
9583 /* Clear the old state */
9584 clear_dd_indices(dd, 0, 0);
9587 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9588 TRUE, cgs_gl, state_global->x, &ddbox);
9590 get_cg_distribution(fplog, dd, cgs_gl,
9591 state_global->box, &ddbox, state_global->x);
9593 dd_distribute_state(dd, cgs_gl,
9594 state_global, state_local, f);
9596 dd_make_local_cgs(dd, &top_local->cgs);
9598 /* Ensure that we have space for the new distribution */
9599 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9601 if (fr->cutoff_scheme == ecutsGROUP)
9603 calc_cgcm(fplog, 0, dd->ncg_home,
9604 &top_local->cgs, state_local->x, fr->cg_cm);
9607 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9609 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9611 else if (state_local->ddp_count != dd->ddp_count)
9613 if (state_local->ddp_count > dd->ddp_count)
9615 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9618 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9620 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);
9623 /* Clear the old state */
9624 clear_dd_indices(dd, 0, 0);
9626 /* Build the new indices */
9627 rebuild_cgindex(dd, cgs_gl->index, state_local);
9628 make_dd_indices(dd, cgs_gl->index, 0);
9629 ncgindex_set = dd->ncg_home;
9631 if (fr->cutoff_scheme == ecutsGROUP)
9633 /* Redetermine the cg COMs */
9634 calc_cgcm(fplog, 0, dd->ncg_home,
9635 &top_local->cgs, state_local->x, fr->cg_cm);
9638 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9640 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9642 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9643 TRUE, &top_local->cgs, state_local->x, &ddbox);
9645 bRedist = dlbIsOn(comm);
9649 /* We have the full state, only redistribute the cgs */
9651 /* Clear the non-home indices */
9652 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9655 /* Avoid global communication for dim's without pbc and -gcom */
9656 if (!bNStGlobalComm)
9658 copy_rvec(comm->box0, ddbox.box0 );
9659 copy_rvec(comm->box_size, ddbox.box_size);
9661 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9662 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9667 /* For dim's without pbc and -gcom */
9668 copy_rvec(ddbox.box0, comm->box0 );
9669 copy_rvec(ddbox.box_size, comm->box_size);
9671 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9674 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9676 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9679 /* Check if we should sort the charge groups */
9680 if (comm->nstSortCG > 0)
9682 bSortCG = (bMasterState ||
9683 (bRedist && (step % comm->nstSortCG == 0)));
9690 ncg_home_old = dd->ncg_home;
9695 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9697 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9699 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9701 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9704 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9706 &comm->cell_x0, &comm->cell_x1,
9707 dd->ncg_home, fr->cg_cm,
9708 cell_ns_x0, cell_ns_x1, &grid_density);
9712 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9715 switch (fr->cutoff_scheme)
9718 copy_ivec(fr->ns.grid->n, ncells_old);
9719 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9720 state_local->box, cell_ns_x0, cell_ns_x1,
9721 fr->rlistlong, grid_density);
9724 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9727 gmx_incons("unimplemented");
9729 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9730 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9734 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9736 /* Sort the state on charge group position.
9737 * This enables exact restarts from this step.
9738 * It also improves performance by about 15% with larger numbers
9739 * of atoms per node.
9742 /* Fill the ns grid with the home cell,
9743 * so we can sort with the indices.
9745 set_zones_ncg_home(dd);
9747 switch (fr->cutoff_scheme)
9750 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9752 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9754 comm->zones.size[0].bb_x0,
9755 comm->zones.size[0].bb_x1,
9757 comm->zones.dens_zone0,
9760 ncg_moved, bRedist ? comm->moved : NULL,
9761 fr->nbv->grp[eintLocal].kernel_type,
9762 fr->nbv->grp[eintLocal].nbat);
9764 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9767 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9768 0, dd->ncg_home, fr->cg_cm);
9770 copy_ivec(fr->ns.grid->n, ncells_new);
9773 gmx_incons("unimplemented");
9776 bResortAll = bMasterState;
9778 /* Check if we can user the old order and ns grid cell indices
9779 * of the charge groups to sort the charge groups efficiently.
9781 if (ncells_new[XX] != ncells_old[XX] ||
9782 ncells_new[YY] != ncells_old[YY] ||
9783 ncells_new[ZZ] != ncells_old[ZZ])
9790 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9791 gmx_step_str(step, sbuf), dd->ncg_home);
9793 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9794 bResortAll ? -1 : ncg_home_old);
9795 /* Rebuild all the indices */
9796 ga2la_clear(dd->ga2la);
9799 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9802 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9804 /* Setup up the communication and communicate the coordinates */
9805 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9807 /* Set the indices */
9808 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9810 /* Set the charge group boundaries for neighbor searching */
9811 set_cg_boundaries(&comm->zones);
9813 if (fr->cutoff_scheme == ecutsVERLET)
9815 set_zones_size(dd, state_local->box, &ddbox,
9816 bSortCG ? 1 : 0, comm->zones.n);
9819 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9822 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9823 -1,state_local->x,state_local->box);
9826 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9828 /* Extract a local topology from the global topology */
9829 for (i = 0; i < dd->ndim; i++)
9831 np[dd->dim[i]] = comm->cd[i].np;
9833 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9834 comm->cellsize_min, np,
9836 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9837 vsite, top_global, top_local);
9839 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9841 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9843 /* Set up the special atom communication */
9844 n = comm->nat[ddnatZONE];
9845 for (i = ddnatZONE+1; i < ddnatNR; i++)
9850 if (vsite && vsite->n_intercg_vsite)
9852 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9856 if (dd->bInterCGcons || dd->bInterCGsettles)
9858 /* Only for inter-cg constraints we need special code */
9859 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9860 constr, ir->nProjOrder,
9861 top_local->idef.il);
9865 gmx_incons("Unknown special atom type setup");
9870 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9872 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9874 /* Make space for the extra coordinates for virtual site
9875 * or constraint communication.
9877 state_local->natoms = comm->nat[ddnatNR-1];
9878 if (state_local->natoms > state_local->nalloc)
9880 dd_realloc_state(state_local, f, state_local->natoms);
9883 if (fr->bF_NoVirSum)
9885 if (vsite && vsite->n_intercg_vsite)
9887 nat_f_novirsum = comm->nat[ddnatVSITE];
9891 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9893 nat_f_novirsum = dd->nat_tot;
9897 nat_f_novirsum = dd->nat_home;
9906 /* Set the number of atoms required for the force calculation.
9907 * Forces need to be constrained when using a twin-range setup
9908 * or with energy minimization. For simple simulations we could
9909 * avoid some allocation, zeroing and copying, but this is
9910 * probably not worth the complications ande checking.
9912 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9913 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9915 /* We make the all mdatoms up to nat_tot_con.
9916 * We could save some work by only setting invmass
9917 * between nat_tot and nat_tot_con.
9919 /* This call also sets the new number of home particles to dd->nat_home */
9920 atoms2md(top_global, ir,
9921 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9923 /* Now we have the charges we can sort the FE interactions */
9924 dd_sort_local_top(dd, mdatoms, top_local);
9928 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9929 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9930 mdatoms, FALSE, vsite);
9935 /* Make the local shell stuff, currently no communication is done */
9936 make_local_shells(cr, mdatoms, shellfc);
9939 if (ir->implicit_solvent)
9941 make_local_gb(cr, fr->born, ir->gb_algorithm);
9944 setup_bonded_threading(fr, &top_local->idef);
9946 if (!(cr->duty & DUTY_PME))
9948 /* Send the charges and/or c6/sigmas to our PME only node */
9949 gmx_pme_send_parameters(cr,
9951 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9952 mdatoms->chargeA, mdatoms->chargeB,
9953 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9954 mdatoms->sigmaA, mdatoms->sigmaB,
9955 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9960 set_constraints(constr, top_local, ir, mdatoms, cr);
9965 /* Update the local pull groups */
9966 dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
9971 /* Update the local rotation groups */
9972 dd_make_local_rotation_groups(dd, ir->rot);
9975 if (ir->eSwapCoords != eswapNO)
9977 /* Update the local groups needed for ion swapping */
9978 dd_make_local_swap_groups(dd, ir->swap);
9981 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9982 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9984 add_dd_statistics(dd);
9986 /* Make sure we only count the cycles for this DD partitioning */
9987 clear_dd_cycle_counts(dd);
9989 /* Because the order of the atoms might have changed since
9990 * the last vsite construction, we need to communicate the constructing
9991 * atom coordinates again (for spreading the forces this MD step).
9993 dd_move_x_vsites(dd, state_local->box, state_local->x);
9995 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9997 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9999 dd_move_x(dd, state_local->box, state_local->x);
10000 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
10001 -1, state_local->x, state_local->box);
10004 /* Store the partitioning step */
10005 comm->partition_step = step;
10007 /* Increase the DD partitioning counter */
10009 /* The state currently matches this DD partitioning count, store it */
10010 state_local->ddp_count = dd->ddp_count;
10013 /* The DD master node knows the complete cg distribution,
10014 * store the count so we can possibly skip the cg info communication.
10016 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
10019 if (comm->DD_debug > 0)
10021 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
10022 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
10023 "after partitioning");
10026 wallcycle_stop(wcycle, ewcDOMDEC);