2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
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/imd/imd.h"
56 #include "gromacs/legacyheaders/bonded-threading.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/gpu_utils.h"
64 #include "gromacs/legacyheaders/mdatoms.h"
65 #include "gromacs/legacyheaders/mdrun.h"
66 #include "gromacs/legacyheaders/names.h"
67 #include "gromacs/legacyheaders/network.h"
68 #include "gromacs/legacyheaders/nrnb.h"
69 #include "gromacs/legacyheaders/nsgrid.h"
70 #include "gromacs/legacyheaders/shellfc.h"
71 #include "gromacs/legacyheaders/typedefs.h"
72 #include "gromacs/legacyheaders/vsite.h"
73 #include "gromacs/legacyheaders/types/commrec.h"
74 #include "gromacs/legacyheaders/types/constr.h"
75 #include "gromacs/legacyheaders/types/enums.h"
76 #include "gromacs/legacyheaders/types/forcerec.h"
77 #include "gromacs/legacyheaders/types/hw_info.h"
78 #include "gromacs/legacyheaders/types/ifunc.h"
79 #include "gromacs/legacyheaders/types/inputrec.h"
80 #include "gromacs/legacyheaders/types/mdatom.h"
81 #include "gromacs/legacyheaders/types/nrnb.h"
82 #include "gromacs/legacyheaders/types/ns.h"
83 #include "gromacs/legacyheaders/types/nsgrid.h"
84 #include "gromacs/legacyheaders/types/shellfc.h"
85 #include "gromacs/legacyheaders/types/simple.h"
86 #include "gromacs/legacyheaders/types/state.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 #define DDRANK(dd, rank) (rank)
111 #define DDMASTERRANK(dd) (dd->masterrank)
113 typedef struct gmx_domdec_master
115 /* The cell boundaries */
117 /* The global charge group division */
118 int *ncg; /* Number of home charge groups for each node */
119 int *index; /* Index of nnodes+1 into cg */
120 int *cg; /* Global charge group index */
121 int *nat; /* Number of home atoms for each node. */
122 int *ibuf; /* Buffer for communication */
123 rvec *vbuf; /* Buffer for state scattering and gathering */
124 } gmx_domdec_master_t;
128 /* The numbers of charge groups to send and receive for each cell
129 * that requires communication, the last entry contains the total
130 * number of atoms that needs to be communicated.
132 int nsend[DD_MAXIZONE+2];
133 int nrecv[DD_MAXIZONE+2];
134 /* The charge groups to send */
137 /* The atom range for non-in-place communication */
138 int cell2at0[DD_MAXIZONE];
139 int cell2at1[DD_MAXIZONE];
144 int np; /* Number of grid pulses in this dimension */
145 int np_dlb; /* For dlb, for use with edlbAUTO */
146 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
148 gmx_bool bInPlace; /* Can we communicate in place? */
149 } gmx_domdec_comm_dim_t;
153 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
154 real *cell_f; /* State var.: cell boundaries, box relative */
155 real *old_cell_f; /* Temp. var.: old cell size */
156 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
157 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
158 real *bound_min; /* Temp. var.: lower limit for cell boundary */
159 real *bound_max; /* Temp. var.: upper limit for cell boundary */
160 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
161 real *buf_ncd; /* Temp. var. */
164 #define DD_NLOAD_MAX 9
166 /* Here floats are accurate enough, since these variables
167 * only influence the load balancing, not the actual MD results.
194 gmx_cgsort_t *sort_new;
206 /* This enum determines the order of the coordinates.
207 * ddnatHOME and ddnatZONE should be first and second,
208 * the others can be ordered as wanted.
211 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
215 edlbAUTO, edlbNO, edlbYES, edlbNR
217 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
221 int dim; /* The dimension */
222 gmx_bool dim_match; /* Tells if DD and PME dims match */
223 int nslab; /* The number of PME slabs in this dimension */
224 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
225 int *pp_min; /* The minimum pp node location, size nslab */
226 int *pp_max; /* The maximum pp node location,size nslab */
227 int maxshift; /* The maximum shift for coordinate redistribution in PME */
232 real min0; /* The minimum bottom of this zone */
233 real max1; /* The maximum top of this zone */
234 real min1; /* The minimum top of this zone */
235 real mch0; /* The maximum bottom communicaton height for this zone */
236 real mch1; /* The maximum top communicaton height for this zone */
237 real p1_0; /* The bottom value of the first cell in this zone */
238 real p1_1; /* The top value of the first cell in this zone */
243 gmx_domdec_ind_t ind;
250 } dd_comm_setup_work_t;
252 typedef struct gmx_domdec_comm
254 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
255 * unless stated otherwise.
258 /* The number of decomposition dimensions for PME, 0: no PME */
260 /* The number of nodes doing PME (PP/PME or only PME) */
264 /* The communication setup including the PME only nodes */
265 gmx_bool bCartesianPP_PME;
268 int *pmenodes; /* size npmenodes */
269 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
270 * but with bCartesianPP_PME */
271 gmx_ddpme_t ddpme[2];
273 /* The DD particle-particle nodes only */
274 gmx_bool bCartesianPP;
275 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
277 /* The global charge groups */
280 /* Should we sort the cgs */
282 gmx_domdec_sort_t *sort;
284 /* Are there charge groups? */
287 /* Are there bonded and multi-body interactions between charge groups? */
288 gmx_bool bInterCGBondeds;
289 gmx_bool bInterCGMultiBody;
291 /* Data for the optional bonded interaction atom communication range */
298 /* Is eDLB=edlbAUTO locked such that we currently can't turn it on? */
299 gmx_bool bDLB_locked;
300 /* Are we actually using DLB? */
301 gmx_bool bDynLoadBal;
303 /* Cell sizes for static load balancing, first index cartesian */
306 /* The width of the communicated boundaries */
309 /* The minimum cell size (including triclinic correction) */
311 /* For dlb, for use with edlbAUTO */
312 rvec cellsize_min_dlb;
313 /* The lower limit for the DD cell size with DLB */
315 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
316 gmx_bool bVacDLBNoLimit;
318 /* With PME load balancing we set limits on DLB */
319 gmx_bool bPMELoadBalDLBLimits;
320 /* DLB needs to take into account that we want to allow this maximum
321 * cut-off (for PME load balancing), this could limit cell boundaries.
323 real PMELoadBal_max_cutoff;
325 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
327 /* box0 and box_size are required with dim's without pbc and -gcom */
331 /* The cell boundaries */
335 /* The old location of the cell boundaries, to check cg displacements */
339 /* The communication setup and charge group boundaries for the zones */
340 gmx_domdec_zones_t zones;
342 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
343 * cell boundaries of neighboring cells for dynamic load balancing.
345 gmx_ddzone_t zone_d1[2];
346 gmx_ddzone_t zone_d2[2][2];
348 /* The coordinate/force communication setup and indices */
349 gmx_domdec_comm_dim_t cd[DIM];
350 /* The maximum number of cells to communicate with in one dimension */
353 /* Which cg distribution is stored on the master node */
354 int master_cg_ddp_count;
356 /* The number of cg's received from the direct neighbors */
357 int zone_ncg1[DD_MAXZONE];
359 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
362 /* Array for signalling if atoms have moved to another domain */
366 /* Communication buffer for general use */
370 /* Communication buffer for general use */
373 /* Temporary storage for thread parallel communication setup */
375 dd_comm_setup_work_t *dth;
377 /* Communication buffers only used with multiple grid pulses */
382 /* Communication buffers for local redistribution */
384 int cggl_flag_nalloc[DIM*2];
386 int cgcm_state_nalloc[DIM*2];
388 /* Cell sizes for dynamic load balancing */
389 gmx_domdec_root_t **root;
393 real cell_f_max0[DIM];
394 real cell_f_min1[DIM];
396 /* Stuff for load communication */
397 gmx_bool bRecordLoad;
398 gmx_domdec_load_t *load;
399 int nrank_gpu_shared;
401 MPI_Comm *mpi_comm_load;
402 MPI_Comm mpi_comm_gpu_shared;
405 /* Maximum DLB scaling per load balancing step in percent */
409 float cycl[ddCyclNr];
410 int cycl_n[ddCyclNr];
411 float cycl_max[ddCyclNr];
412 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
416 /* How many times have did we have load measurements */
418 /* How many times have we collected the load measurements */
422 double sum_nat[ddnatNR-ddnatZONE];
432 /* The last partition step */
433 gmx_int64_t partition_step;
441 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
444 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
445 #define DD_FLAG_NRCG 65535
446 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
447 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
449 /* Zone permutation required to obtain consecutive charge groups
450 * for neighbor searching.
452 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
454 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
455 * components see only j zones with that component 0.
458 /* The DD zone order */
459 static const ivec dd_zo[DD_MAXZONE] =
460 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
465 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
470 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
475 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
477 /* Factors used to avoid problems due to rounding issues */
478 #define DD_CELL_MARGIN 1.0001
479 #define DD_CELL_MARGIN2 1.00005
480 /* Factor to account for pressure scaling during nstlist steps */
481 #define DD_PRES_SCALE_MARGIN 1.02
483 /* Turn on DLB when the load imbalance causes this amount of total loss.
484 * There is a bit of overhead with DLB and it's difficult to achieve
485 * a load imbalance of less than 2% with DLB.
487 #define DD_PERF_LOSS_DLB_ON 0.02
489 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
490 #define DD_PERF_LOSS_WARN 0.05
492 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
494 /* Use separate MPI send and receive commands
495 * when nnodes <= GMX_DD_NNODES_SENDRECV.
496 * This saves memory (and some copying for small nnodes).
497 * For high parallelization scatter and gather calls are used.
499 #define GMX_DD_NNODES_SENDRECV 4
503 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
505 static void index2xyz(ivec nc,int ind,ivec xyz)
507 xyz[XX] = ind % nc[XX];
508 xyz[YY] = (ind / nc[XX]) % nc[YY];
509 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
513 /* This order is required to minimize the coordinate communication in PME
514 * which uses decomposition in the x direction.
516 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
518 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
520 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
521 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
522 xyz[ZZ] = ind % nc[ZZ];
525 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
530 ddindex = dd_index(dd->nc, c);
531 if (dd->comm->bCartesianPP_PME)
533 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
535 else if (dd->comm->bCartesianPP)
538 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
549 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
551 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
554 int ddglatnr(gmx_domdec_t *dd, int i)
564 if (i >= dd->comm->nat[ddnatNR-1])
566 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
568 atnr = dd->gatindex[i] + 1;
574 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
576 return &dd->comm->cgs_gl;
579 static void vec_rvec_init(vec_rvec_t *v)
585 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
589 v->nalloc = over_alloc_dd(n);
590 srenew(v->v, v->nalloc);
594 void dd_store_state(gmx_domdec_t *dd, t_state *state)
598 if (state->ddp_count != dd->ddp_count)
600 gmx_incons("The state does not the domain decomposition state");
603 state->ncg_gl = dd->ncg_home;
604 if (state->ncg_gl > state->cg_gl_nalloc)
606 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
607 srenew(state->cg_gl, state->cg_gl_nalloc);
609 for (i = 0; i < state->ncg_gl; i++)
611 state->cg_gl[i] = dd->index_gl[i];
614 state->ddp_count_cg_gl = dd->ddp_count;
617 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
619 return &dd->comm->zones;
622 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
623 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
625 gmx_domdec_zones_t *zones;
628 zones = &dd->comm->zones;
631 while (icg >= zones->izone[izone].cg1)
640 else if (izone < zones->nizone)
642 *jcg0 = zones->izone[izone].jcg0;
646 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
647 icg, izone, zones->nizone);
650 *jcg1 = zones->izone[izone].jcg1;
652 for (d = 0; d < dd->ndim; d++)
655 shift0[dim] = zones->izone[izone].shift0[dim];
656 shift1[dim] = zones->izone[izone].shift1[dim];
657 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
659 /* A conservative approach, this can be optimized */
666 int dd_natoms_vsite(gmx_domdec_t *dd)
668 return dd->comm->nat[ddnatVSITE];
671 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
673 *at_start = dd->comm->nat[ddnatCON-1];
674 *at_end = dd->comm->nat[ddnatCON];
677 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
679 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
680 int *index, *cgindex;
681 gmx_domdec_comm_t *comm;
682 gmx_domdec_comm_dim_t *cd;
683 gmx_domdec_ind_t *ind;
684 rvec shift = {0, 0, 0}, *buf, *rbuf;
685 gmx_bool bPBC, bScrew;
689 cgindex = dd->cgindex;
694 nat_tot = dd->nat_home;
695 for (d = 0; d < dd->ndim; d++)
697 bPBC = (dd->ci[dd->dim[d]] == 0);
698 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
701 copy_rvec(box[dd->dim[d]], shift);
704 for (p = 0; p < cd->np; p++)
711 for (i = 0; i < ind->nsend[nzone]; i++)
713 at0 = cgindex[index[i]];
714 at1 = cgindex[index[i]+1];
715 for (j = at0; j < at1; j++)
717 copy_rvec(x[j], buf[n]);
724 for (i = 0; i < ind->nsend[nzone]; i++)
726 at0 = cgindex[index[i]];
727 at1 = cgindex[index[i]+1];
728 for (j = at0; j < at1; j++)
730 /* We need to shift the coordinates */
731 rvec_add(x[j], shift, buf[n]);
738 for (i = 0; i < ind->nsend[nzone]; i++)
740 at0 = cgindex[index[i]];
741 at1 = cgindex[index[i]+1];
742 for (j = at0; j < at1; j++)
745 buf[n][XX] = x[j][XX] + shift[XX];
747 * This operation requires a special shift force
748 * treatment, which is performed in calc_vir.
750 buf[n][YY] = box[YY][YY] - x[j][YY];
751 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
763 rbuf = comm->vbuf2.v;
765 /* Send and receive the coordinates */
766 dd_sendrecv_rvec(dd, d, dddirBackward,
767 buf, ind->nsend[nzone+1],
768 rbuf, ind->nrecv[nzone+1]);
772 for (zone = 0; zone < nzone; zone++)
774 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
776 copy_rvec(rbuf[j], x[i]);
781 nat_tot += ind->nrecv[nzone+1];
787 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
789 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
790 int *index, *cgindex;
791 gmx_domdec_comm_t *comm;
792 gmx_domdec_comm_dim_t *cd;
793 gmx_domdec_ind_t *ind;
797 gmx_bool bShiftForcesNeedPbc, bScrew;
801 cgindex = dd->cgindex;
805 nzone = comm->zones.n/2;
806 nat_tot = dd->nat_tot;
807 for (d = dd->ndim-1; d >= 0; d--)
809 /* Only forces in domains near the PBC boundaries need to
810 consider PBC in the treatment of fshift */
811 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
812 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
813 if (fshift == NULL && !bScrew)
815 bShiftForcesNeedPbc = FALSE;
817 /* Determine which shift vector we need */
823 for (p = cd->np-1; p >= 0; p--)
826 nat_tot -= ind->nrecv[nzone+1];
833 sbuf = comm->vbuf2.v;
835 for (zone = 0; zone < nzone; zone++)
837 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
839 copy_rvec(f[i], sbuf[j]);
844 /* Communicate the forces */
845 dd_sendrecv_rvec(dd, d, dddirForward,
846 sbuf, ind->nrecv[nzone+1],
847 buf, ind->nsend[nzone+1]);
849 /* Add the received forces */
851 if (!bShiftForcesNeedPbc)
853 for (i = 0; i < ind->nsend[nzone]; i++)
855 at0 = cgindex[index[i]];
856 at1 = cgindex[index[i]+1];
857 for (j = at0; j < at1; j++)
859 rvec_inc(f[j], buf[n]);
866 /* fshift should always be defined if this function is
867 * called when bShiftForcesNeedPbc is true */
868 assert(NULL != fshift);
869 for (i = 0; i < ind->nsend[nzone]; i++)
871 at0 = cgindex[index[i]];
872 at1 = cgindex[index[i]+1];
873 for (j = at0; j < at1; j++)
875 rvec_inc(f[j], buf[n]);
876 /* Add this force to the shift force */
877 rvec_inc(fshift[is], buf[n]);
884 for (i = 0; i < ind->nsend[nzone]; i++)
886 at0 = cgindex[index[i]];
887 at1 = cgindex[index[i]+1];
888 for (j = at0; j < at1; j++)
890 /* Rotate the force */
891 f[j][XX] += buf[n][XX];
892 f[j][YY] -= buf[n][YY];
893 f[j][ZZ] -= buf[n][ZZ];
896 /* Add this force to the shift force */
897 rvec_inc(fshift[is], buf[n]);
908 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
910 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
911 int *index, *cgindex;
912 gmx_domdec_comm_t *comm;
913 gmx_domdec_comm_dim_t *cd;
914 gmx_domdec_ind_t *ind;
919 cgindex = dd->cgindex;
921 buf = &comm->vbuf.v[0][0];
924 nat_tot = dd->nat_home;
925 for (d = 0; d < dd->ndim; d++)
928 for (p = 0; p < cd->np; p++)
933 for (i = 0; i < ind->nsend[nzone]; i++)
935 at0 = cgindex[index[i]];
936 at1 = cgindex[index[i]+1];
937 for (j = at0; j < at1; j++)
950 rbuf = &comm->vbuf2.v[0][0];
952 /* Send and receive the coordinates */
953 dd_sendrecv_real(dd, d, dddirBackward,
954 buf, ind->nsend[nzone+1],
955 rbuf, ind->nrecv[nzone+1]);
959 for (zone = 0; zone < nzone; zone++)
961 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
968 nat_tot += ind->nrecv[nzone+1];
974 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
976 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
977 int *index, *cgindex;
978 gmx_domdec_comm_t *comm;
979 gmx_domdec_comm_dim_t *cd;
980 gmx_domdec_ind_t *ind;
985 cgindex = dd->cgindex;
987 buf = &comm->vbuf.v[0][0];
989 nzone = comm->zones.n/2;
990 nat_tot = dd->nat_tot;
991 for (d = dd->ndim-1; d >= 0; d--)
994 for (p = cd->np-1; p >= 0; p--)
997 nat_tot -= ind->nrecv[nzone+1];
1004 sbuf = &comm->vbuf2.v[0][0];
1006 for (zone = 0; zone < nzone; zone++)
1008 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
1015 /* Communicate the forces */
1016 dd_sendrecv_real(dd, d, dddirForward,
1017 sbuf, ind->nrecv[nzone+1],
1018 buf, ind->nsend[nzone+1]);
1020 /* Add the received forces */
1022 for (i = 0; i < ind->nsend[nzone]; i++)
1024 at0 = cgindex[index[i]];
1025 at1 = cgindex[index[i]+1];
1026 for (j = at0; j < at1; j++)
1037 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1039 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",
1041 zone->min0, zone->max1,
1042 zone->mch0, zone->mch0,
1043 zone->p1_0, zone->p1_1);
1047 #define DDZONECOMM_MAXZONE 5
1048 #define DDZONECOMM_BUFSIZE 3
1050 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1051 int ddimind, int direction,
1052 gmx_ddzone_t *buf_s, int n_s,
1053 gmx_ddzone_t *buf_r, int n_r)
1055 #define ZBS DDZONECOMM_BUFSIZE
1056 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1057 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1060 for (i = 0; i < n_s; i++)
1062 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1063 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1064 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1065 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1066 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1067 vbuf_s[i*ZBS+1][2] = 0;
1068 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1069 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1070 vbuf_s[i*ZBS+2][2] = 0;
1073 dd_sendrecv_rvec(dd, ddimind, direction,
1077 for (i = 0; i < n_r; i++)
1079 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1080 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1081 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1082 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1083 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1084 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1085 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1091 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1092 rvec cell_ns_x0, rvec cell_ns_x1)
1094 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1096 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1097 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1098 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1099 rvec extr_s[2], extr_r[2];
1101 real dist_d, c = 0, det;
1102 gmx_domdec_comm_t *comm;
1103 gmx_bool bPBC, bUse;
1107 for (d = 1; d < dd->ndim; d++)
1110 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1111 zp->min0 = cell_ns_x0[dim];
1112 zp->max1 = cell_ns_x1[dim];
1113 zp->min1 = cell_ns_x1[dim];
1114 zp->mch0 = cell_ns_x0[dim];
1115 zp->mch1 = cell_ns_x1[dim];
1116 zp->p1_0 = cell_ns_x0[dim];
1117 zp->p1_1 = cell_ns_x1[dim];
1120 for (d = dd->ndim-2; d >= 0; d--)
1123 bPBC = (dim < ddbox->npbcdim);
1125 /* Use an rvec to store two reals */
1126 extr_s[d][0] = comm->cell_f0[d+1];
1127 extr_s[d][1] = comm->cell_f1[d+1];
1128 extr_s[d][2] = comm->cell_f1[d+1];
1131 /* Store the extremes in the backward sending buffer,
1132 * so the get updated separately from the forward communication.
1134 for (d1 = d; d1 < dd->ndim-1; d1++)
1136 /* We invert the order to be able to use the same loop for buf_e */
1137 buf_s[pos].min0 = extr_s[d1][1];
1138 buf_s[pos].max1 = extr_s[d1][0];
1139 buf_s[pos].min1 = extr_s[d1][2];
1140 buf_s[pos].mch0 = 0;
1141 buf_s[pos].mch1 = 0;
1142 /* Store the cell corner of the dimension we communicate along */
1143 buf_s[pos].p1_0 = comm->cell_x0[dim];
1144 buf_s[pos].p1_1 = 0;
1148 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1151 if (dd->ndim == 3 && d == 0)
1153 buf_s[pos] = comm->zone_d2[0][1];
1155 buf_s[pos] = comm->zone_d1[0];
1159 /* We only need to communicate the extremes
1160 * in the forward direction
1162 npulse = comm->cd[d].np;
1165 /* Take the minimum to avoid double communication */
1166 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1170 /* Without PBC we should really not communicate over
1171 * the boundaries, but implementing that complicates
1172 * the communication setup and therefore we simply
1173 * do all communication, but ignore some data.
1175 npulse_min = npulse;
1177 for (p = 0; p < npulse_min; p++)
1179 /* Communicate the extremes forward */
1180 bUse = (bPBC || dd->ci[dim] > 0);
1182 dd_sendrecv_rvec(dd, d, dddirForward,
1183 extr_s+d, dd->ndim-d-1,
1184 extr_r+d, dd->ndim-d-1);
1188 for (d1 = d; d1 < dd->ndim-1; d1++)
1190 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1191 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1192 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1198 for (p = 0; p < npulse; p++)
1200 /* Communicate all the zone information backward */
1201 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1203 dd_sendrecv_ddzone(dd, d, dddirBackward,
1210 for (d1 = d+1; d1 < dd->ndim; d1++)
1212 /* Determine the decrease of maximum required
1213 * communication height along d1 due to the distance along d,
1214 * this avoids a lot of useless atom communication.
1216 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1218 if (ddbox->tric_dir[dim])
1220 /* c is the off-diagonal coupling between the cell planes
1221 * along directions d and d1.
1223 c = ddbox->v[dim][dd->dim[d1]][dim];
1229 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1232 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1236 /* A negative value signals out of range */
1242 /* Accumulate the extremes over all pulses */
1243 for (i = 0; i < buf_size; i++)
1247 buf_e[i] = buf_r[i];
1253 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1254 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1255 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1258 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1266 if (bUse && dh[d1] >= 0)
1268 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1269 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1272 /* Copy the received buffer to the send buffer,
1273 * to pass the data through with the next pulse.
1275 buf_s[i] = buf_r[i];
1277 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1278 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1280 /* Store the extremes */
1283 for (d1 = d; d1 < dd->ndim-1; d1++)
1285 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1286 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1287 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1291 if (d == 1 || (d == 0 && dd->ndim == 3))
1293 for (i = d; i < 2; i++)
1295 comm->zone_d2[1-d][i] = buf_e[pos];
1301 comm->zone_d1[1] = buf_e[pos];
1311 for (i = 0; i < 2; i++)
1315 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1317 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1318 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1324 for (i = 0; i < 2; i++)
1326 for (j = 0; j < 2; j++)
1330 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1332 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1333 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1337 for (d = 1; d < dd->ndim; d++)
1339 comm->cell_f_max0[d] = extr_s[d-1][0];
1340 comm->cell_f_min1[d] = extr_s[d-1][1];
1343 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1344 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1349 static void dd_collect_cg(gmx_domdec_t *dd,
1350 t_state *state_local)
1352 gmx_domdec_master_t *ma = NULL;
1353 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1355 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1357 /* The master has the correct distribution */
1361 if (state_local->ddp_count == dd->ddp_count)
1363 /* The local state and DD are in sync, use the DD indices */
1364 ncg_home = dd->ncg_home;
1366 nat_home = dd->nat_home;
1368 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1370 /* The DD is out of sync with the local state, but we have stored
1371 * the cg indices with the local state, so we can use those.
1375 cgs_gl = &dd->comm->cgs_gl;
1377 ncg_home = state_local->ncg_gl;
1378 cg = state_local->cg_gl;
1380 for (i = 0; i < ncg_home; i++)
1382 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1387 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1401 /* Collect the charge group and atom counts on the master */
1402 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1407 for (i = 0; i < dd->nnodes; i++)
1409 ma->ncg[i] = ma->ibuf[2*i];
1410 ma->nat[i] = ma->ibuf[2*i+1];
1411 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1414 /* Make byte counts and indices */
1415 for (i = 0; i < dd->nnodes; i++)
1417 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1418 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1422 fprintf(debug, "Initial charge group distribution: ");
1423 for (i = 0; i < dd->nnodes; i++)
1425 fprintf(debug, " %d", ma->ncg[i]);
1427 fprintf(debug, "\n");
1431 /* Collect the charge group indices on the master */
1433 ncg_home*sizeof(int), cg,
1434 DDMASTER(dd) ? ma->ibuf : NULL,
1435 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1436 DDMASTER(dd) ? ma->cg : NULL);
1438 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1441 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1444 gmx_domdec_master_t *ma;
1445 int n, i, c, a, nalloc = 0;
1454 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1455 dd->rank, dd->mpi_comm_all);
1460 /* Copy the master coordinates to the global array */
1461 cgs_gl = &dd->comm->cgs_gl;
1463 n = DDMASTERRANK(dd);
1465 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1467 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1469 copy_rvec(lv[a++], v[c]);
1473 for (n = 0; n < dd->nnodes; n++)
1477 if (ma->nat[n] > nalloc)
1479 nalloc = over_alloc_dd(ma->nat[n]);
1480 srenew(buf, nalloc);
1483 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1484 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1487 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1489 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1491 copy_rvec(buf[a++], v[c]);
1500 static void get_commbuffer_counts(gmx_domdec_t *dd,
1501 int **counts, int **disps)
1503 gmx_domdec_master_t *ma;
1508 /* Make the rvec count and displacment arrays */
1510 *disps = ma->ibuf + dd->nnodes;
1511 for (n = 0; n < dd->nnodes; n++)
1513 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1514 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1518 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1521 gmx_domdec_master_t *ma;
1522 int *rcounts = NULL, *disps = NULL;
1531 get_commbuffer_counts(dd, &rcounts, &disps);
1536 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1540 cgs_gl = &dd->comm->cgs_gl;
1543 for (n = 0; n < dd->nnodes; n++)
1545 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1547 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1549 copy_rvec(buf[a++], v[c]);
1556 void dd_collect_vec(gmx_domdec_t *dd,
1557 t_state *state_local, rvec *lv, rvec *v)
1559 dd_collect_cg(dd, state_local);
1561 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1563 dd_collect_vec_sendrecv(dd, lv, v);
1567 dd_collect_vec_gatherv(dd, lv, v);
1572 void dd_collect_state(gmx_domdec_t *dd,
1573 t_state *state_local, t_state *state)
1577 nh = state->nhchainlength;
1581 for (i = 0; i < efptNR; i++)
1583 state->lambda[i] = state_local->lambda[i];
1585 state->fep_state = state_local->fep_state;
1586 state->veta = state_local->veta;
1587 state->vol0 = state_local->vol0;
1588 copy_mat(state_local->box, state->box);
1589 copy_mat(state_local->boxv, state->boxv);
1590 copy_mat(state_local->svir_prev, state->svir_prev);
1591 copy_mat(state_local->fvir_prev, state->fvir_prev);
1592 copy_mat(state_local->pres_prev, state->pres_prev);
1594 for (i = 0; i < state_local->ngtc; i++)
1596 for (j = 0; j < nh; j++)
1598 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1599 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1601 state->therm_integral[i] = state_local->therm_integral[i];
1603 for (i = 0; i < state_local->nnhpres; i++)
1605 for (j = 0; j < nh; j++)
1607 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1608 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1612 for (est = 0; est < estNR; est++)
1614 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1619 dd_collect_vec(dd, state_local, state_local->x, state->x);
1622 dd_collect_vec(dd, state_local, state_local->v, state->v);
1625 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1628 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1630 case estDISRE_INITF:
1631 case estDISRE_RM3TAV:
1632 case estORIRE_INITF:
1636 gmx_incons("Unknown state entry encountered in dd_collect_state");
1642 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1648 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1651 state->nalloc = over_alloc_dd(nalloc);
1653 for (est = 0; est < estNR; est++)
1655 if (EST_DISTR(est) && (state->flags & (1<<est)))
1660 srenew(state->x, state->nalloc);
1663 srenew(state->v, state->nalloc);
1666 srenew(state->sd_X, state->nalloc);
1669 srenew(state->cg_p, state->nalloc);
1671 case estDISRE_INITF:
1672 case estDISRE_RM3TAV:
1673 case estORIRE_INITF:
1675 /* No reallocation required */
1678 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1685 srenew(*f, state->nalloc);
1689 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1692 if (nalloc > fr->cg_nalloc)
1696 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1698 fr->cg_nalloc = over_alloc_dd(nalloc);
1699 srenew(fr->cginfo, fr->cg_nalloc);
1700 if (fr->cutoff_scheme == ecutsGROUP)
1702 srenew(fr->cg_cm, fr->cg_nalloc);
1705 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1707 /* We don't use charge groups, we use x in state to set up
1708 * the atom communication.
1710 dd_realloc_state(state, f, nalloc);
1714 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1717 gmx_domdec_master_t *ma;
1718 int n, i, c, a, nalloc = 0;
1725 for (n = 0; n < dd->nnodes; n++)
1729 if (ma->nat[n] > nalloc)
1731 nalloc = over_alloc_dd(ma->nat[n]);
1732 srenew(buf, nalloc);
1734 /* Use lv as a temporary buffer */
1736 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1738 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1740 copy_rvec(v[c], buf[a++]);
1743 if (a != ma->nat[n])
1745 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1750 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1751 DDRANK(dd, n), n, dd->mpi_comm_all);
1756 n = DDMASTERRANK(dd);
1758 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1760 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1762 copy_rvec(v[c], lv[a++]);
1769 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1770 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1775 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1778 gmx_domdec_master_t *ma;
1779 int *scounts = NULL, *disps = NULL;
1787 get_commbuffer_counts(dd, &scounts, &disps);
1791 for (n = 0; n < dd->nnodes; n++)
1793 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1795 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1797 copy_rvec(v[c], buf[a++]);
1803 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1806 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1808 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1810 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1814 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1818 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1821 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1822 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1823 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1825 if (dfhist->nlambda > 0)
1827 int nlam = dfhist->nlambda;
1828 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1829 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1830 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1831 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1832 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1833 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1835 for (i = 0; i < nlam; i++)
1837 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1838 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1839 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1840 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1841 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1842 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1847 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1848 t_state *state, t_state *state_local,
1853 nh = state->nhchainlength;
1857 for (i = 0; i < efptNR; i++)
1859 state_local->lambda[i] = state->lambda[i];
1861 state_local->fep_state = state->fep_state;
1862 state_local->veta = state->veta;
1863 state_local->vol0 = state->vol0;
1864 copy_mat(state->box, state_local->box);
1865 copy_mat(state->box_rel, state_local->box_rel);
1866 copy_mat(state->boxv, state_local->boxv);
1867 copy_mat(state->svir_prev, state_local->svir_prev);
1868 copy_mat(state->fvir_prev, state_local->fvir_prev);
1869 copy_df_history(&state_local->dfhist, &state->dfhist);
1870 for (i = 0; i < state_local->ngtc; i++)
1872 for (j = 0; j < nh; j++)
1874 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1875 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1877 state_local->therm_integral[i] = state->therm_integral[i];
1879 for (i = 0; i < state_local->nnhpres; i++)
1881 for (j = 0; j < nh; j++)
1883 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1884 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1888 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1889 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1890 dd_bcast(dd, sizeof(real), &state_local->veta);
1891 dd_bcast(dd, sizeof(real), &state_local->vol0);
1892 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1893 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1894 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1895 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1896 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1897 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1898 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1899 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1900 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1901 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1903 /* communicate df_history -- required for restarting from checkpoint */
1904 dd_distribute_dfhist(dd, &state_local->dfhist);
1906 if (dd->nat_home > state_local->nalloc)
1908 dd_realloc_state(state_local, f, dd->nat_home);
1910 for (i = 0; i < estNR; i++)
1912 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1917 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1920 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1923 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1926 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1928 case estDISRE_INITF:
1929 case estDISRE_RM3TAV:
1930 case estORIRE_INITF:
1932 /* Not implemented yet */
1935 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1941 static char dim2char(int dim)
1947 case XX: c = 'X'; break;
1948 case YY: c = 'Y'; break;
1949 case ZZ: c = 'Z'; break;
1950 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1956 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1957 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1959 rvec grid_s[2], *grid_r = NULL, cx, r;
1960 char fname[STRLEN], buf[22];
1962 int a, i, d, z, y, x;
1966 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1967 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1971 snew(grid_r, 2*dd->nnodes);
1974 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1978 for (d = 0; d < DIM; d++)
1980 for (i = 0; i < DIM; i++)
1988 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1990 tric[d][i] = box[i][d]/box[i][i];
1999 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2000 out = gmx_fio_fopen(fname, "w");
2001 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2003 for (i = 0; i < dd->nnodes; i++)
2005 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2006 for (d = 0; d < DIM; d++)
2008 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2010 for (z = 0; z < 2; z++)
2012 for (y = 0; y < 2; y++)
2014 for (x = 0; x < 2; x++)
2016 cx[XX] = grid_r[i*2+x][XX];
2017 cx[YY] = grid_r[i*2+y][YY];
2018 cx[ZZ] = grid_r[i*2+z][ZZ];
2020 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2021 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2025 for (d = 0; d < DIM; d++)
2027 for (x = 0; x < 4; x++)
2031 case 0: y = 1 + i*8 + 2*x; break;
2032 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2033 case 2: y = 1 + i*8 + x; break;
2035 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2039 gmx_fio_fclose(out);
2044 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2045 gmx_mtop_t *mtop, t_commrec *cr,
2046 int natoms, rvec x[], matrix box)
2048 char fname[STRLEN], buf[22];
2050 int i, ii, resnr, c;
2051 char *atomname, *resname;
2058 natoms = dd->comm->nat[ddnatVSITE];
2061 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2063 out = gmx_fio_fopen(fname, "w");
2065 fprintf(out, "TITLE %s\n", title);
2066 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2067 for (i = 0; i < natoms; i++)
2069 ii = dd->gatindex[i];
2070 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2071 if (i < dd->comm->nat[ddnatZONE])
2074 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2080 else if (i < dd->comm->nat[ddnatVSITE])
2082 b = dd->comm->zones.n;
2086 b = dd->comm->zones.n + 1;
2088 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2089 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2091 fprintf(out, "TER\n");
2093 gmx_fio_fclose(out);
2096 real dd_cutoff_mbody(gmx_domdec_t *dd)
2098 gmx_domdec_comm_t *comm;
2105 if (comm->bInterCGBondeds)
2107 if (comm->cutoff_mbody > 0)
2109 r = comm->cutoff_mbody;
2113 /* cutoff_mbody=0 means we do not have DLB */
2114 r = comm->cellsize_min[dd->dim[0]];
2115 for (di = 1; di < dd->ndim; di++)
2117 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2119 if (comm->bBondComm)
2121 r = std::max(r, comm->cutoff_mbody);
2125 r = std::min(r, comm->cutoff);
2133 real dd_cutoff_twobody(gmx_domdec_t *dd)
2137 r_mb = dd_cutoff_mbody(dd);
2139 return std::max(dd->comm->cutoff, r_mb);
2143 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2147 nc = dd->nc[dd->comm->cartpmedim];
2148 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2149 copy_ivec(coord, coord_pme);
2150 coord_pme[dd->comm->cartpmedim] =
2151 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2154 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2156 /* Here we assign a PME node to communicate with this DD node
2157 * by assuming that the major index of both is x.
2158 * We add cr->npmenodes/2 to obtain an even distribution.
2160 return (ddindex*npme + npme/2)/ndd;
2163 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2165 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2168 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2170 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2173 static int *dd_pmenodes(t_commrec *cr)
2178 snew(pmenodes, cr->npmenodes);
2180 for (i = 0; i < cr->dd->nnodes; i++)
2182 p0 = cr_ddindex2pmeindex(cr, i);
2183 p1 = cr_ddindex2pmeindex(cr, i+1);
2184 if (i+1 == cr->dd->nnodes || p1 > p0)
2188 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2190 pmenodes[n] = i + 1 + n;
2198 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2206 if (dd->comm->bCartesian) {
2207 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2208 dd_coords2pmecoords(dd,coords,coords_pme);
2209 copy_ivec(dd->ntot,nc);
2210 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2211 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2213 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2215 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2221 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2226 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2228 gmx_domdec_comm_t *comm;
2230 int ddindex, nodeid = -1;
2232 comm = cr->dd->comm;
2237 if (comm->bCartesianPP_PME)
2240 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2245 ddindex = dd_index(cr->dd->nc, coords);
2246 if (comm->bCartesianPP)
2248 nodeid = comm->ddindex2simnodeid[ddindex];
2254 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2266 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2269 gmx_domdec_comm_t *comm;
2276 /* This assumes a uniform x domain decomposition grid cell size */
2277 if (comm->bCartesianPP_PME)
2280 ivec coord, coord_pme;
2281 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2282 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2284 /* This is a PP node */
2285 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2286 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2290 else if (comm->bCartesianPP)
2292 if (sim_nodeid < dd->nnodes)
2294 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2299 /* This assumes DD cells with identical x coordinates
2300 * are numbered sequentially.
2302 if (dd->comm->pmenodes == NULL)
2304 if (sim_nodeid < dd->nnodes)
2306 /* The DD index equals the nodeid */
2307 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2313 while (sim_nodeid > dd->comm->pmenodes[i])
2317 if (sim_nodeid < dd->comm->pmenodes[i])
2319 pmenode = dd->comm->pmenodes[i];
2327 void get_pme_nnodes(const gmx_domdec_t *dd,
2328 int *npmenodes_x, int *npmenodes_y)
2332 *npmenodes_x = dd->comm->npmenodes_x;
2333 *npmenodes_y = dd->comm->npmenodes_y;
2342 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2344 gmx_bool bPMEOnlyNode;
2346 if (DOMAINDECOMP(cr))
2348 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2352 bPMEOnlyNode = FALSE;
2355 return bPMEOnlyNode;
2358 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2359 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2363 ivec coord, coord_pme;
2367 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2370 for (x = 0; x < dd->nc[XX]; x++)
2372 for (y = 0; y < dd->nc[YY]; y++)
2374 for (z = 0; z < dd->nc[ZZ]; z++)
2376 if (dd->comm->bCartesianPP_PME)
2381 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2382 if (dd->ci[XX] == coord_pme[XX] &&
2383 dd->ci[YY] == coord_pme[YY] &&
2384 dd->ci[ZZ] == coord_pme[ZZ])
2386 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2391 /* The slab corresponds to the nodeid in the PME group */
2392 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2394 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2401 /* The last PP-only node is the peer node */
2402 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2406 fprintf(debug, "Receive coordinates from PP ranks:");
2407 for (x = 0; x < *nmy_ddnodes; x++)
2409 fprintf(debug, " %d", (*my_ddnodes)[x]);
2411 fprintf(debug, "\n");
2415 static gmx_bool receive_vir_ener(t_commrec *cr)
2417 gmx_domdec_comm_t *comm;
2422 if (cr->npmenodes < cr->dd->nnodes)
2424 comm = cr->dd->comm;
2425 if (comm->bCartesianPP_PME)
2427 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2430 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2431 coords[comm->cartpmedim]++;
2432 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2435 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2436 if (dd_simnode2pmenode(cr, rank) == pmenode)
2438 /* This is not the last PP node for pmenode */
2446 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2447 if (cr->sim_nodeid+1 < cr->nnodes &&
2448 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2450 /* This is not the last PP node for pmenode */
2459 static void set_zones_ncg_home(gmx_domdec_t *dd)
2461 gmx_domdec_zones_t *zones;
2464 zones = &dd->comm->zones;
2466 zones->cg_range[0] = 0;
2467 for (i = 1; i < zones->n+1; i++)
2469 zones->cg_range[i] = dd->ncg_home;
2471 /* zone_ncg1[0] should always be equal to ncg_home */
2472 dd->comm->zone_ncg1[0] = dd->ncg_home;
2475 static void rebuild_cgindex(gmx_domdec_t *dd,
2476 const int *gcgs_index, t_state *state)
2478 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2481 dd_cg_gl = dd->index_gl;
2482 cgindex = dd->cgindex;
2485 for (i = 0; i < state->ncg_gl; i++)
2489 dd_cg_gl[i] = cg_gl;
2490 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2494 dd->ncg_home = state->ncg_gl;
2497 set_zones_ncg_home(dd);
2500 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2502 while (cg >= cginfo_mb->cg_end)
2507 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2510 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2511 t_forcerec *fr, char *bLocalCG)
2513 cginfo_mb_t *cginfo_mb;
2519 cginfo_mb = fr->cginfo_mb;
2520 cginfo = fr->cginfo;
2522 for (cg = cg0; cg < cg1; cg++)
2524 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2528 if (bLocalCG != NULL)
2530 for (cg = cg0; cg < cg1; cg++)
2532 bLocalCG[index_gl[cg]] = TRUE;
2537 static void make_dd_indices(gmx_domdec_t *dd,
2538 const int *gcgs_index, int cg_start)
2540 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2541 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2544 if (dd->nat_tot > dd->gatindex_nalloc)
2546 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2547 srenew(dd->gatindex, dd->gatindex_nalloc);
2550 nzone = dd->comm->zones.n;
2551 zone2cg = dd->comm->zones.cg_range;
2552 zone_ncg1 = dd->comm->zone_ncg1;
2553 index_gl = dd->index_gl;
2554 gatindex = dd->gatindex;
2555 bCGs = dd->comm->bCGs;
2557 if (zone2cg[1] != dd->ncg_home)
2559 gmx_incons("dd->ncg_zone is not up to date");
2562 /* Make the local to global and global to local atom index */
2563 a = dd->cgindex[cg_start];
2564 for (zone = 0; zone < nzone; zone++)
2572 cg0 = zone2cg[zone];
2574 cg1 = zone2cg[zone+1];
2575 cg1_p1 = cg0 + zone_ncg1[zone];
2577 for (cg = cg0; cg < cg1; cg++)
2582 /* Signal that this cg is from more than one pulse away */
2585 cg_gl = index_gl[cg];
2588 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2591 ga2la_set(dd->ga2la, a_gl, a, zone1);
2597 gatindex[a] = cg_gl;
2598 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2605 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2611 if (bLocalCG == NULL)
2615 for (i = 0; i < dd->ncg_tot; i++)
2617 if (!bLocalCG[dd->index_gl[i]])
2620 "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);
2625 for (i = 0; i < ncg_sys; i++)
2632 if (ngl != dd->ncg_tot)
2634 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);
2641 static void check_index_consistency(gmx_domdec_t *dd,
2642 int natoms_sys, int ncg_sys,
2645 int nerr, ngl, i, a, cell;
2650 if (dd->comm->DD_debug > 1)
2652 snew(have, natoms_sys);
2653 for (a = 0; a < dd->nat_tot; a++)
2655 if (have[dd->gatindex[a]] > 0)
2657 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);
2661 have[dd->gatindex[a]] = a + 1;
2667 snew(have, dd->nat_tot);
2670 for (i = 0; i < natoms_sys; i++)
2672 if (ga2la_get(dd->ga2la, i, &a, &cell))
2674 if (a >= dd->nat_tot)
2676 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);
2682 if (dd->gatindex[a] != i)
2684 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);
2691 if (ngl != dd->nat_tot)
2694 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2695 dd->rank, where, ngl, dd->nat_tot);
2697 for (a = 0; a < dd->nat_tot; a++)
2702 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2703 dd->rank, where, a+1, dd->gatindex[a]+1);
2708 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2712 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2713 dd->rank, where, nerr);
2717 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2724 /* Clear the whole list without searching */
2725 ga2la_clear(dd->ga2la);
2729 for (i = a_start; i < dd->nat_tot; i++)
2731 ga2la_del(dd->ga2la, dd->gatindex[i]);
2735 bLocalCG = dd->comm->bLocalCG;
2738 for (i = cg_start; i < dd->ncg_tot; i++)
2740 bLocalCG[dd->index_gl[i]] = FALSE;
2744 dd_clear_local_vsite_indices(dd);
2746 if (dd->constraints)
2748 dd_clear_local_constraint_indices(dd);
2752 /* This function should be used for moving the domain boudaries during DLB,
2753 * for obtaining the minimum cell size. It checks the initially set limit
2754 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2755 * and, possibly, a longer cut-off limit set for PME load balancing.
2757 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2761 cellsize_min = comm->cellsize_min[dim];
2763 if (!comm->bVacDLBNoLimit)
2765 /* The cut-off might have changed, e.g. by PME load balacning,
2766 * from the value used to set comm->cellsize_min, so check it.
2768 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2770 if (comm->bPMELoadBalDLBLimits)
2772 /* Check for the cut-off limit set by the PME load balancing */
2773 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2777 return cellsize_min;
2780 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2783 real grid_jump_limit;
2785 /* The distance between the boundaries of cells at distance
2786 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2787 * and by the fact that cells should not be shifted by more than
2788 * half their size, such that cg's only shift by one cell
2789 * at redecomposition.
2791 grid_jump_limit = comm->cellsize_limit;
2792 if (!comm->bVacDLBNoLimit)
2794 if (comm->bPMELoadBalDLBLimits)
2796 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2798 grid_jump_limit = std::max(grid_jump_limit,
2799 cutoff/comm->cd[dim_ind].np);
2802 return grid_jump_limit;
2805 static gmx_bool check_grid_jump(gmx_int64_t step,
2811 gmx_domdec_comm_t *comm;
2820 for (d = 1; d < dd->ndim; d++)
2823 limit = grid_jump_limit(comm, cutoff, d);
2824 bfac = ddbox->box_size[dim];
2825 if (ddbox->tric_dir[dim])
2827 bfac *= ddbox->skew_fac[dim];
2829 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2830 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2838 /* This error should never be triggered under normal
2839 * circumstances, but you never know ...
2841 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.",
2842 gmx_step_str(step, buf),
2843 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2851 static int dd_load_count(gmx_domdec_comm_t *comm)
2853 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2856 static float dd_force_load(gmx_domdec_comm_t *comm)
2863 if (comm->eFlop > 1)
2865 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2870 load = comm->cycl[ddCyclF];
2871 if (comm->cycl_n[ddCyclF] > 1)
2873 /* Subtract the maximum of the last n cycle counts
2874 * to get rid of possible high counts due to other sources,
2875 * for instance system activity, that would otherwise
2876 * affect the dynamic load balancing.
2878 load -= comm->cycl_max[ddCyclF];
2882 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2884 float gpu_wait, gpu_wait_sum;
2886 gpu_wait = comm->cycl[ddCyclWaitGPU];
2887 if (comm->cycl_n[ddCyclF] > 1)
2889 /* We should remove the WaitGPU time of the same MD step
2890 * as the one with the maximum F time, since the F time
2891 * and the wait time are not independent.
2892 * Furthermore, the step for the max F time should be chosen
2893 * the same on all ranks that share the same GPU.
2894 * But to keep the code simple, we remove the average instead.
2895 * The main reason for artificially long times at some steps
2896 * is spurious CPU activity or MPI time, so we don't expect
2897 * that changes in the GPU wait time matter a lot here.
2899 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2901 /* Sum the wait times over the ranks that share the same GPU */
2902 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2903 comm->mpi_comm_gpu_shared);
2904 /* Replace the wait time by the average over the ranks */
2905 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2913 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2915 gmx_domdec_comm_t *comm;
2920 snew(*dim_f, dd->nc[dim]+1);
2922 for (i = 1; i < dd->nc[dim]; i++)
2924 if (comm->slb_frac[dim])
2926 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2930 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2933 (*dim_f)[dd->nc[dim]] = 1;
2936 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2938 int pmeindex, slab, nso, i;
2941 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2947 ddpme->dim = dimind;
2949 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2951 ddpme->nslab = (ddpme->dim == 0 ?
2952 dd->comm->npmenodes_x :
2953 dd->comm->npmenodes_y);
2955 if (ddpme->nslab <= 1)
2960 nso = dd->comm->npmenodes/ddpme->nslab;
2961 /* Determine for each PME slab the PP location range for dimension dim */
2962 snew(ddpme->pp_min, ddpme->nslab);
2963 snew(ddpme->pp_max, ddpme->nslab);
2964 for (slab = 0; slab < ddpme->nslab; slab++)
2966 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2967 ddpme->pp_max[slab] = 0;
2969 for (i = 0; i < dd->nnodes; i++)
2971 ddindex2xyz(dd->nc, i, xyz);
2972 /* For y only use our y/z slab.
2973 * This assumes that the PME x grid size matches the DD grid size.
2975 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2977 pmeindex = ddindex2pmeindex(dd, i);
2980 slab = pmeindex/nso;
2984 slab = pmeindex % ddpme->nslab;
2986 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2987 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2991 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2994 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2996 if (dd->comm->ddpme[0].dim == XX)
2998 return dd->comm->ddpme[0].maxshift;
3006 int dd_pme_maxshift_y(gmx_domdec_t *dd)
3008 if (dd->comm->ddpme[0].dim == YY)
3010 return dd->comm->ddpme[0].maxshift;
3012 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3014 return dd->comm->ddpme[1].maxshift;
3022 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3023 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3025 gmx_domdec_comm_t *comm;
3028 real range, pme_boundary;
3032 nc = dd->nc[ddpme->dim];
3035 if (!ddpme->dim_match)
3037 /* PP decomposition is not along dim: the worst situation */
3040 else if (ns <= 3 || (bUniform && ns == nc))
3042 /* The optimal situation */
3047 /* We need to check for all pme nodes which nodes they
3048 * could possibly need to communicate with.
3050 xmin = ddpme->pp_min;
3051 xmax = ddpme->pp_max;
3052 /* Allow for atoms to be maximally 2/3 times the cut-off
3053 * out of their DD cell. This is a reasonable balance between
3054 * between performance and support for most charge-group/cut-off
3057 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3058 /* Avoid extra communication when we are exactly at a boundary */
3062 for (s = 0; s < ns; s++)
3064 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3065 pme_boundary = (real)s/ns;
3068 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3070 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3074 pme_boundary = (real)(s+1)/ns;
3077 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3079 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3086 ddpme->maxshift = sh;
3090 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3091 ddpme->dim, ddpme->maxshift);
3095 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3099 for (d = 0; d < dd->ndim; d++)
3102 if (dim < ddbox->nboundeddim &&
3103 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3104 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3106 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",
3107 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3108 dd->nc[dim], dd->comm->cellsize_limit);
3114 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3117 /* Set the domain boundaries. Use for static (or no) load balancing,
3118 * and also for the starting state for dynamic load balancing.
3119 * setmode determine if and where the boundaries are stored, use enum above.
3120 * Returns the number communication pulses in npulse.
3122 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3123 int setmode, ivec npulse)
3125 gmx_domdec_comm_t *comm;
3128 real *cell_x, cell_dx, cellsize;
3132 for (d = 0; d < DIM; d++)
3134 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3136 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3139 cell_dx = ddbox->box_size[d]/dd->nc[d];
3142 case setcellsizeslbMASTER:
3143 for (j = 0; j < dd->nc[d]+1; j++)
3145 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3148 case setcellsizeslbLOCAL:
3149 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3150 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3155 cellsize = cell_dx*ddbox->skew_fac[d];
3156 while (cellsize*npulse[d] < comm->cutoff)
3160 cellsize_min[d] = cellsize;
3164 /* Statically load balanced grid */
3165 /* Also when we are not doing a master distribution we determine
3166 * all cell borders in a loop to obtain identical values
3167 * to the master distribution case and to determine npulse.
3169 if (setmode == setcellsizeslbMASTER)
3171 cell_x = dd->ma->cell_x[d];
3175 snew(cell_x, dd->nc[d]+1);
3177 cell_x[0] = ddbox->box0[d];
3178 for (j = 0; j < dd->nc[d]; j++)
3180 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3181 cell_x[j+1] = cell_x[j] + cell_dx;
3182 cellsize = cell_dx*ddbox->skew_fac[d];
3183 while (cellsize*npulse[d] < comm->cutoff &&
3184 npulse[d] < dd->nc[d]-1)
3188 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3190 if (setmode == setcellsizeslbLOCAL)
3192 comm->cell_x0[d] = cell_x[dd->ci[d]];
3193 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3195 if (setmode != setcellsizeslbMASTER)
3200 /* The following limitation is to avoid that a cell would receive
3201 * some of its own home charge groups back over the periodic boundary.
3202 * Double charge groups cause trouble with the global indices.
3204 if (d < ddbox->npbcdim &&
3205 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3207 char error_string[STRLEN];
3209 sprintf(error_string,
3210 "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",
3211 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3213 dd->nc[d], dd->nc[d],
3214 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3216 if (setmode == setcellsizeslbLOCAL)
3218 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3222 gmx_fatal(FARGS, error_string);
3227 if (!comm->bDynLoadBal)
3229 copy_rvec(cellsize_min, comm->cellsize_min);
3232 for (d = 0; d < comm->npmedecompdim; d++)
3234 set_pme_maxshift(dd, &comm->ddpme[d],
3235 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3236 comm->ddpme[d].slb_dim_f);
3241 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3242 int d, int dim, gmx_domdec_root_t *root,
3244 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3246 gmx_domdec_comm_t *comm;
3247 int ncd, i, j, nmin, nmin_old;
3248 gmx_bool bLimLo, bLimHi;
3250 real fac, halfway, cellsize_limit_f_i, region_size;
3251 gmx_bool bPBC, bLastHi = FALSE;
3252 int nrange[] = {range[0], range[1]};
3254 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3260 bPBC = (dim < ddbox->npbcdim);
3262 cell_size = root->buf_ncd;
3266 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3269 /* First we need to check if the scaling does not make cells
3270 * smaller than the smallest allowed size.
3271 * We need to do this iteratively, since if a cell is too small,
3272 * it needs to be enlarged, which makes all the other cells smaller,
3273 * which could in turn make another cell smaller than allowed.
3275 for (i = range[0]; i < range[1]; i++)
3277 root->bCellMin[i] = FALSE;
3283 /* We need the total for normalization */
3285 for (i = range[0]; i < range[1]; i++)
3287 if (root->bCellMin[i] == FALSE)
3289 fac += cell_size[i];
3292 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3293 /* Determine the cell boundaries */
3294 for (i = range[0]; i < range[1]; i++)
3296 if (root->bCellMin[i] == FALSE)
3298 cell_size[i] *= fac;
3299 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3301 cellsize_limit_f_i = 0;
3305 cellsize_limit_f_i = cellsize_limit_f;
3307 if (cell_size[i] < cellsize_limit_f_i)
3309 root->bCellMin[i] = TRUE;
3310 cell_size[i] = cellsize_limit_f_i;
3314 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3317 while (nmin > nmin_old);
3320 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3321 /* For this check we should not use DD_CELL_MARGIN,
3322 * but a slightly smaller factor,
3323 * since rounding could get use below the limit.
3325 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3328 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",
3329 gmx_step_str(step, buf),
3330 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3331 ncd, comm->cellsize_min[dim]);
3334 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3338 /* Check if the boundary did not displace more than halfway
3339 * each of the cells it bounds, as this could cause problems,
3340 * especially when the differences between cell sizes are large.
3341 * If changes are applied, they will not make cells smaller
3342 * than the cut-off, as we check all the boundaries which
3343 * might be affected by a change and if the old state was ok,
3344 * the cells will at most be shrunk back to their old size.
3346 for (i = range[0]+1; i < range[1]; i++)
3348 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3349 if (root->cell_f[i] < halfway)
3351 root->cell_f[i] = halfway;
3352 /* Check if the change also causes shifts of the next boundaries */
3353 for (j = i+1; j < range[1]; j++)
3355 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3357 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3361 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3362 if (root->cell_f[i] > halfway)
3364 root->cell_f[i] = halfway;
3365 /* Check if the change also causes shifts of the next boundaries */
3366 for (j = i-1; j >= range[0]+1; j--)
3368 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3370 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3377 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3378 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3379 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3380 * for a and b nrange is used */
3383 /* Take care of the staggering of the cell boundaries */
3386 for (i = range[0]; i < range[1]; i++)
3388 root->cell_f_max0[i] = root->cell_f[i];
3389 root->cell_f_min1[i] = root->cell_f[i+1];
3394 for (i = range[0]+1; i < range[1]; i++)
3396 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3397 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3398 if (bLimLo && bLimHi)
3400 /* Both limits violated, try the best we can */
3401 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3402 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3403 nrange[0] = range[0];
3405 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3408 nrange[1] = range[1];
3409 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3415 /* root->cell_f[i] = root->bound_min[i]; */
3416 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3419 else if (bLimHi && !bLastHi)
3422 if (nrange[1] < range[1]) /* found a LimLo before */
3424 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3425 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3426 nrange[0] = nrange[1];
3428 root->cell_f[i] = root->bound_max[i];
3430 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3432 nrange[1] = range[1];
3435 if (nrange[1] < range[1]) /* found last a LimLo */
3437 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3438 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3439 nrange[0] = nrange[1];
3440 nrange[1] = range[1];
3441 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3443 else if (nrange[0] > range[0]) /* found at least one LimHi */
3445 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3452 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3453 int d, int dim, gmx_domdec_root_t *root,
3454 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3455 gmx_bool bUniform, gmx_int64_t step)
3457 gmx_domdec_comm_t *comm;
3458 int ncd, d1, i, pos;
3460 real load_aver, load_i, imbalance, change, change_max, sc;
3461 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3465 int range[] = { 0, 0 };
3469 /* Convert the maximum change from the input percentage to a fraction */
3470 change_limit = comm->dlb_scale_lim*0.01;
3474 bPBC = (dim < ddbox->npbcdim);
3476 cell_size = root->buf_ncd;
3478 /* Store the original boundaries */
3479 for (i = 0; i < ncd+1; i++)
3481 root->old_cell_f[i] = root->cell_f[i];
3485 for (i = 0; i < ncd; i++)
3487 cell_size[i] = 1.0/ncd;
3490 else if (dd_load_count(comm) > 0)
3492 load_aver = comm->load[d].sum_m/ncd;
3494 for (i = 0; i < ncd; i++)
3496 /* Determine the relative imbalance of cell i */
3497 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3498 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3499 /* Determine the change of the cell size using underrelaxation */
3500 change = -relax*imbalance;
3501 change_max = std::max(change_max, std::max(change, -change));
3503 /* Limit the amount of scaling.
3504 * We need to use the same rescaling for all cells in one row,
3505 * otherwise the load balancing might not converge.
3508 if (change_max > change_limit)
3510 sc *= change_limit/change_max;
3512 for (i = 0; i < ncd; i++)
3514 /* Determine the relative imbalance of cell i */
3515 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3516 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3517 /* Determine the change of the cell size using underrelaxation */
3518 change = -sc*imbalance;
3519 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3523 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3524 cellsize_limit_f *= DD_CELL_MARGIN;
3525 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3526 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3527 if (ddbox->tric_dir[dim])
3529 cellsize_limit_f /= ddbox->skew_fac[dim];
3530 dist_min_f /= ddbox->skew_fac[dim];
3532 if (bDynamicBox && d > 0)
3534 dist_min_f *= DD_PRES_SCALE_MARGIN;
3536 if (d > 0 && !bUniform)
3538 /* Make sure that the grid is not shifted too much */
3539 for (i = 1; i < ncd; i++)
3541 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3543 gmx_incons("Inconsistent DD boundary staggering limits!");
3545 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3546 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3549 root->bound_min[i] += 0.5*space;
3551 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3552 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3555 root->bound_max[i] += 0.5*space;
3560 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3562 root->cell_f_max0[i-1] + dist_min_f,
3563 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3564 root->cell_f_min1[i] - dist_min_f);
3569 root->cell_f[0] = 0;
3570 root->cell_f[ncd] = 1;
3571 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3574 /* After the checks above, the cells should obey the cut-off
3575 * restrictions, but it does not hurt to check.
3577 for (i = 0; i < ncd; i++)
3581 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3582 dim, i, root->cell_f[i], root->cell_f[i+1]);
3585 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3586 root->cell_f[i+1] - root->cell_f[i] <
3587 cellsize_limit_f/DD_CELL_MARGIN)
3591 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3592 gmx_step_str(step, buf), dim2char(dim), i,
3593 (root->cell_f[i+1] - root->cell_f[i])
3594 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3599 /* Store the cell boundaries of the lower dimensions at the end */
3600 for (d1 = 0; d1 < d; d1++)
3602 root->cell_f[pos++] = comm->cell_f0[d1];
3603 root->cell_f[pos++] = comm->cell_f1[d1];
3606 if (d < comm->npmedecompdim)
3608 /* The master determines the maximum shift for
3609 * the coordinate communication between separate PME nodes.
3611 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3613 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3616 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3620 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3621 gmx_ddbox_t *ddbox, int dimind)
3623 gmx_domdec_comm_t *comm;
3628 /* Set the cell dimensions */
3629 dim = dd->dim[dimind];
3630 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3631 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3632 if (dim >= ddbox->nboundeddim)
3634 comm->cell_x0[dim] += ddbox->box0[dim];
3635 comm->cell_x1[dim] += ddbox->box0[dim];
3639 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3640 int d, int dim, real *cell_f_row,
3643 gmx_domdec_comm_t *comm;
3649 /* Each node would only need to know two fractions,
3650 * but it is probably cheaper to broadcast the whole array.
3652 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3653 0, comm->mpi_comm_load[d]);
3655 /* Copy the fractions for this dimension from the buffer */
3656 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3657 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3658 /* The whole array was communicated, so set the buffer position */
3659 pos = dd->nc[dim] + 1;
3660 for (d1 = 0; d1 <= d; d1++)
3664 /* Copy the cell fractions of the lower dimensions */
3665 comm->cell_f0[d1] = cell_f_row[pos++];
3666 comm->cell_f1[d1] = cell_f_row[pos++];
3668 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3670 /* Convert the communicated shift from float to int */
3671 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3674 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3678 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3679 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3680 gmx_bool bUniform, gmx_int64_t step)
3682 gmx_domdec_comm_t *comm;
3684 gmx_bool bRowMember, bRowRoot;
3689 for (d = 0; d < dd->ndim; d++)
3694 for (d1 = d; d1 < dd->ndim; d1++)
3696 if (dd->ci[dd->dim[d1]] > 0)
3709 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3710 ddbox, bDynamicBox, bUniform, step);
3711 cell_f_row = comm->root[d]->cell_f;
3715 cell_f_row = comm->cell_f_row;
3717 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3722 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3726 /* This function assumes the box is static and should therefore
3727 * not be called when the box has changed since the last
3728 * call to dd_partition_system.
3730 for (d = 0; d < dd->ndim; d++)
3732 relative_to_absolute_cell_bounds(dd, ddbox, d);
3738 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3739 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3740 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3741 gmx_wallcycle_t wcycle)
3743 gmx_domdec_comm_t *comm;
3750 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3751 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3752 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3754 else if (bDynamicBox)
3756 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3759 /* Set the dimensions for which no DD is used */
3760 for (dim = 0; dim < DIM; dim++)
3762 if (dd->nc[dim] == 1)
3764 comm->cell_x0[dim] = 0;
3765 comm->cell_x1[dim] = ddbox->box_size[dim];
3766 if (dim >= ddbox->nboundeddim)
3768 comm->cell_x0[dim] += ddbox->box0[dim];
3769 comm->cell_x1[dim] += ddbox->box0[dim];
3775 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3778 gmx_domdec_comm_dim_t *cd;
3780 for (d = 0; d < dd->ndim; d++)
3782 cd = &dd->comm->cd[d];
3783 np = npulse[dd->dim[d]];
3784 if (np > cd->np_nalloc)
3788 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3789 dim2char(dd->dim[d]), np);
3791 if (DDMASTER(dd) && cd->np_nalloc > 0)
3793 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3795 srenew(cd->ind, np);
3796 for (i = cd->np_nalloc; i < np; i++)
3798 cd->ind[i].index = NULL;
3799 cd->ind[i].nalloc = 0;
3808 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3809 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3810 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3811 gmx_wallcycle_t wcycle)
3813 gmx_domdec_comm_t *comm;
3819 /* Copy the old cell boundaries for the cg displacement check */
3820 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3821 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3823 if (comm->bDynLoadBal)
3827 check_box_size(dd, ddbox);
3829 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3833 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3834 realloc_comm_ind(dd, npulse);
3839 for (d = 0; d < DIM; d++)
3841 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3842 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3847 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3849 rvec cell_ns_x0, rvec cell_ns_x1,
3852 gmx_domdec_comm_t *comm;
3857 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3859 dim = dd->dim[dim_ind];
3861 /* Without PBC we don't have restrictions on the outer cells */
3862 if (!(dim >= ddbox->npbcdim &&
3863 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3864 comm->bDynLoadBal &&
3865 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3866 comm->cellsize_min[dim])
3869 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",
3870 gmx_step_str(step, buf), dim2char(dim),
3871 comm->cell_x1[dim] - comm->cell_x0[dim],
3872 ddbox->skew_fac[dim],
3873 dd->comm->cellsize_min[dim],
3874 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3878 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3880 /* Communicate the boundaries and update cell_ns_x0/1 */
3881 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3882 if (dd->bGridJump && dd->ndim > 1)
3884 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3889 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3893 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3901 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3902 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3911 static void check_screw_box(matrix box)
3913 /* Mathematical limitation */
3914 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3916 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3919 /* Limitation due to the asymmetry of the eighth shell method */
3920 if (box[ZZ][YY] != 0)
3922 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3926 static void distribute_cg(FILE *fplog, gmx_int64_t step,
3927 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3930 gmx_domdec_master_t *ma;
3931 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3932 int i, icg, j, k, k0, k1, d;
3936 real nrcg, inv_ncg, pos_d;
3942 if (tmp_ind == NULL)
3944 snew(tmp_nalloc, dd->nnodes);
3945 snew(tmp_ind, dd->nnodes);
3946 for (i = 0; i < dd->nnodes; i++)
3948 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3949 snew(tmp_ind[i], tmp_nalloc[i]);
3953 /* Clear the count */
3954 for (i = 0; i < dd->nnodes; i++)
3960 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3962 cgindex = cgs->index;
3964 /* Compute the center of geometry for all charge groups */
3965 for (icg = 0; icg < cgs->nr; icg++)
3968 k1 = cgindex[icg+1];
3972 copy_rvec(pos[k0], cg_cm);
3979 for (k = k0; (k < k1); k++)
3981 rvec_inc(cg_cm, pos[k]);
3983 for (d = 0; (d < DIM); d++)
3985 cg_cm[d] *= inv_ncg;
3988 /* Put the charge group in the box and determine the cell index */
3989 for (d = DIM-1; d >= 0; d--)
3992 if (d < dd->npbcdim)
3994 bScrew = (dd->bScrewPBC && d == XX);
3995 if (tric_dir[d] && dd->nc[d] > 1)
3997 /* Use triclinic coordintates for this dimension */
3998 for (j = d+1; j < DIM; j++)
4000 pos_d += cg_cm[j]*tcm[j][d];
4003 while (pos_d >= box[d][d])
4006 rvec_dec(cg_cm, box[d]);
4009 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4010 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4012 for (k = k0; (k < k1); k++)
4014 rvec_dec(pos[k], box[d]);
4017 pos[k][YY] = box[YY][YY] - pos[k][YY];
4018 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4025 rvec_inc(cg_cm, box[d]);
4028 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4029 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4031 for (k = k0; (k < k1); k++)
4033 rvec_inc(pos[k], box[d]);
4036 pos[k][YY] = box[YY][YY] - pos[k][YY];
4037 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4042 /* This could be done more efficiently */
4044 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4049 i = dd_index(dd->nc, ind);
4050 if (ma->ncg[i] == tmp_nalloc[i])
4052 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4053 srenew(tmp_ind[i], tmp_nalloc[i]);
4055 tmp_ind[i][ma->ncg[i]] = icg;
4057 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4061 for (i = 0; i < dd->nnodes; i++)
4064 for (k = 0; k < ma->ncg[i]; k++)
4066 ma->cg[k1++] = tmp_ind[i][k];
4069 ma->index[dd->nnodes] = k1;
4071 for (i = 0; i < dd->nnodes; i++)
4081 fprintf(fplog, "Charge group distribution at step %s:",
4082 gmx_step_str(step, buf));
4083 for (i = 0; i < dd->nnodes; i++)
4085 fprintf(fplog, " %d", ma->ncg[i]);
4087 fprintf(fplog, "\n");
4091 static void get_cg_distribution(FILE *fplog, gmx_int64_t step, gmx_domdec_t *dd,
4092 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4095 gmx_domdec_master_t *ma = NULL;
4098 int *ibuf, buf2[2] = { 0, 0 };
4099 gmx_bool bMaster = DDMASTER(dd);
4107 check_screw_box(box);
4110 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4112 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4113 for (i = 0; i < dd->nnodes; i++)
4115 ma->ibuf[2*i] = ma->ncg[i];
4116 ma->ibuf[2*i+1] = ma->nat[i];
4124 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4126 dd->ncg_home = buf2[0];
4127 dd->nat_home = buf2[1];
4128 dd->ncg_tot = dd->ncg_home;
4129 dd->nat_tot = dd->nat_home;
4130 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4132 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4133 srenew(dd->index_gl, dd->cg_nalloc);
4134 srenew(dd->cgindex, dd->cg_nalloc+1);
4138 for (i = 0; i < dd->nnodes; i++)
4140 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4141 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4146 bMaster ? ma->ibuf : NULL,
4147 bMaster ? ma->ibuf+dd->nnodes : NULL,
4148 bMaster ? ma->cg : NULL,
4149 dd->ncg_home*sizeof(int), dd->index_gl);
4151 /* Determine the home charge group sizes */
4153 for (i = 0; i < dd->ncg_home; i++)
4155 cg_gl = dd->index_gl[i];
4157 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4162 fprintf(debug, "Home charge groups:\n");
4163 for (i = 0; i < dd->ncg_home; i++)
4165 fprintf(debug, " %d", dd->index_gl[i]);
4168 fprintf(debug, "\n");
4171 fprintf(debug, "\n");
4175 static int compact_and_copy_vec_at(int ncg, int *move,
4178 rvec *src, gmx_domdec_comm_t *comm,
4181 int m, icg, i, i0, i1, nrcg;
4187 for (m = 0; m < DIM*2; m++)
4193 for (icg = 0; icg < ncg; icg++)
4195 i1 = cgindex[icg+1];
4201 /* Compact the home array in place */
4202 for (i = i0; i < i1; i++)
4204 copy_rvec(src[i], src[home_pos++]);
4210 /* Copy to the communication buffer */
4212 pos_vec[m] += 1 + vec*nrcg;
4213 for (i = i0; i < i1; i++)
4215 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4217 pos_vec[m] += (nvec - vec - 1)*nrcg;
4221 home_pos += i1 - i0;
4229 static int compact_and_copy_vec_cg(int ncg, int *move,
4231 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4234 int m, icg, i0, i1, nrcg;
4240 for (m = 0; m < DIM*2; m++)
4246 for (icg = 0; icg < ncg; icg++)
4248 i1 = cgindex[icg+1];
4254 /* Compact the home array in place */
4255 copy_rvec(src[icg], src[home_pos++]);
4261 /* Copy to the communication buffer */
4262 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4263 pos_vec[m] += 1 + nrcg*nvec;
4275 static int compact_ind(int ncg, int *move,
4276 int *index_gl, int *cgindex,
4278 gmx_ga2la_t ga2la, char *bLocalCG,
4281 int cg, nat, a0, a1, a, a_gl;
4286 for (cg = 0; cg < ncg; cg++)
4292 /* Compact the home arrays in place.
4293 * Anything that can be done here avoids access to global arrays.
4295 cgindex[home_pos] = nat;
4296 for (a = a0; a < a1; a++)
4299 gatindex[nat] = a_gl;
4300 /* The cell number stays 0, so we don't need to set it */
4301 ga2la_change_la(ga2la, a_gl, nat);
4304 index_gl[home_pos] = index_gl[cg];
4305 cginfo[home_pos] = cginfo[cg];
4306 /* The charge group remains local, so bLocalCG does not change */
4311 /* Clear the global indices */
4312 for (a = a0; a < a1; a++)
4314 ga2la_del(ga2la, gatindex[a]);
4318 bLocalCG[index_gl[cg]] = FALSE;
4322 cgindex[home_pos] = nat;
4327 static void clear_and_mark_ind(int ncg, int *move,
4328 int *index_gl, int *cgindex, int *gatindex,
4329 gmx_ga2la_t ga2la, char *bLocalCG,
4334 for (cg = 0; cg < ncg; cg++)
4340 /* Clear the global indices */
4341 for (a = a0; a < a1; a++)
4343 ga2la_del(ga2la, gatindex[a]);
4347 bLocalCG[index_gl[cg]] = FALSE;
4349 /* Signal that this cg has moved using the ns cell index.
4350 * Here we set it to -1. fill_grid will change it
4351 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4353 cell_index[cg] = -1;
4358 static void print_cg_move(FILE *fplog,
4360 gmx_int64_t step, int cg, int dim, int dir,
4361 gmx_bool bHaveCgcmOld, real limitd,
4362 rvec cm_old, rvec cm_new, real pos_d)
4364 gmx_domdec_comm_t *comm;
4369 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4372 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4373 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4374 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4378 /* We don't have a limiting distance available: don't print it */
4379 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4380 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4381 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4383 fprintf(fplog, "distance out of cell %f\n",
4384 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4387 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4388 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4390 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4391 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4392 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4394 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4395 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4397 comm->cell_x0[dim], comm->cell_x1[dim]);
4400 static void cg_move_error(FILE *fplog,
4402 gmx_int64_t step, int cg, int dim, int dir,
4403 gmx_bool bHaveCgcmOld, real limitd,
4404 rvec cm_old, rvec cm_new, real pos_d)
4408 print_cg_move(fplog, dd, step, cg, dim, dir,
4409 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4411 print_cg_move(stderr, dd, step, cg, dim, dir,
4412 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4414 "%s moved too far between two domain decomposition steps\n"
4415 "This usually means that your system is not well equilibrated",
4416 dd->comm->bCGs ? "A charge group" : "An atom");
4419 static void rotate_state_atom(t_state *state, int a)
4423 for (est = 0; est < estNR; est++)
4425 if (EST_DISTR(est) && (state->flags & (1<<est)))
4430 /* Rotate the complete state; for a rectangular box only */
4431 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4432 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4435 state->v[a][YY] = -state->v[a][YY];
4436 state->v[a][ZZ] = -state->v[a][ZZ];
4439 state->sd_X[a][YY] = -state->sd_X[a][YY];
4440 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4443 state->cg_p[a][YY] = -state->cg_p[a][YY];
4444 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4446 case estDISRE_INITF:
4447 case estDISRE_RM3TAV:
4448 case estORIRE_INITF:
4450 /* These are distances, so not affected by rotation */
4453 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4459 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4461 if (natoms > comm->moved_nalloc)
4463 /* Contents should be preserved here */
4464 comm->moved_nalloc = over_alloc_dd(natoms);
4465 srenew(comm->moved, comm->moved_nalloc);
4471 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4474 ivec tric_dir, matrix tcm,
4475 rvec cell_x0, rvec cell_x1,
4476 rvec limitd, rvec limit0, rvec limit1,
4478 int cg_start, int cg_end,
4483 int cg, k, k0, k1, d, dim, d2;
4488 real inv_ncg, pos_d;
4491 npbcdim = dd->npbcdim;
4493 for (cg = cg_start; cg < cg_end; cg++)
4500 copy_rvec(state->x[k0], cm_new);
4507 for (k = k0; (k < k1); k++)
4509 rvec_inc(cm_new, state->x[k]);
4511 for (d = 0; (d < DIM); d++)
4513 cm_new[d] = inv_ncg*cm_new[d];
4518 /* Do pbc and check DD cell boundary crossings */
4519 for (d = DIM-1; d >= 0; d--)
4523 bScrew = (dd->bScrewPBC && d == XX);
4524 /* Determine the location of this cg in lattice coordinates */
4528 for (d2 = d+1; d2 < DIM; d2++)
4530 pos_d += cm_new[d2]*tcm[d2][d];
4533 /* Put the charge group in the triclinic unit-cell */
4534 if (pos_d >= cell_x1[d])
4536 if (pos_d >= limit1[d])
4538 cg_move_error(fplog, dd, step, cg, d, 1,
4539 cg_cm != state->x, limitd[d],
4540 cg_cm[cg], cm_new, pos_d);
4543 if (dd->ci[d] == dd->nc[d] - 1)
4545 rvec_dec(cm_new, state->box[d]);
4548 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4549 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4551 for (k = k0; (k < k1); k++)
4553 rvec_dec(state->x[k], state->box[d]);
4556 rotate_state_atom(state, k);
4561 else if (pos_d < cell_x0[d])
4563 if (pos_d < limit0[d])
4565 cg_move_error(fplog, dd, step, cg, d, -1,
4566 cg_cm != state->x, limitd[d],
4567 cg_cm[cg], cm_new, pos_d);
4572 rvec_inc(cm_new, state->box[d]);
4575 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4576 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4578 for (k = k0; (k < k1); k++)
4580 rvec_inc(state->x[k], state->box[d]);
4583 rotate_state_atom(state, k);
4589 else if (d < npbcdim)
4591 /* Put the charge group in the rectangular unit-cell */
4592 while (cm_new[d] >= state->box[d][d])
4594 rvec_dec(cm_new, state->box[d]);
4595 for (k = k0; (k < k1); k++)
4597 rvec_dec(state->x[k], state->box[d]);
4600 while (cm_new[d] < 0)
4602 rvec_inc(cm_new, state->box[d]);
4603 for (k = k0; (k < k1); k++)
4605 rvec_inc(state->x[k], state->box[d]);
4611 copy_rvec(cm_new, cg_cm[cg]);
4613 /* Determine where this cg should go */
4616 for (d = 0; d < dd->ndim; d++)
4621 flag |= DD_FLAG_FW(d);
4627 else if (dev[dim] == -1)
4629 flag |= DD_FLAG_BW(d);
4632 if (dd->nc[dim] > 2)
4643 /* Temporarily store the flag in move */
4644 move[cg] = mc + flag;
4648 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4649 gmx_domdec_t *dd, ivec tric_dir,
4650 t_state *state, rvec **f,
4659 int ncg[DIM*2], nat[DIM*2];
4660 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4661 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4662 int sbuf[2], rbuf[2];
4663 int home_pos_cg, home_pos_at, buf_pos;
4665 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4668 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4670 cginfo_mb_t *cginfo_mb;
4671 gmx_domdec_comm_t *comm;
4673 int nthread, thread;
4677 check_screw_box(state->box);
4681 if (fr->cutoff_scheme == ecutsGROUP)
4686 for (i = 0; i < estNR; i++)
4692 case estX: /* Always present */ break;
4693 case estV: bV = (state->flags & (1<<i)); break;
4694 case estSDX: bSDX = (state->flags & (1<<i)); break;
4695 case estCGP: bCGP = (state->flags & (1<<i)); break;
4698 case estDISRE_INITF:
4699 case estDISRE_RM3TAV:
4700 case estORIRE_INITF:
4702 /* No processing required */
4705 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4710 if (dd->ncg_tot > comm->nalloc_int)
4712 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4713 srenew(comm->buf_int, comm->nalloc_int);
4715 move = comm->buf_int;
4717 /* Clear the count */
4718 for (c = 0; c < dd->ndim*2; c++)
4724 npbcdim = dd->npbcdim;
4726 for (d = 0; (d < DIM); d++)
4728 limitd[d] = dd->comm->cellsize_min[d];
4729 if (d >= npbcdim && dd->ci[d] == 0)
4731 cell_x0[d] = -GMX_FLOAT_MAX;
4735 cell_x0[d] = comm->cell_x0[d];
4737 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4739 cell_x1[d] = GMX_FLOAT_MAX;
4743 cell_x1[d] = comm->cell_x1[d];
4747 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4748 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4752 /* We check after communication if a charge group moved
4753 * more than one cell. Set the pre-comm check limit to float_max.
4755 limit0[d] = -GMX_FLOAT_MAX;
4756 limit1[d] = GMX_FLOAT_MAX;
4760 make_tric_corr_matrix(npbcdim, state->box, tcm);
4762 cgindex = dd->cgindex;
4764 nthread = gmx_omp_nthreads_get(emntDomdec);
4766 /* Compute the center of geometry for all home charge groups
4767 * and put them in the box and determine where they should go.
4769 #pragma omp parallel for num_threads(nthread) schedule(static)
4770 for (thread = 0; thread < nthread; thread++)
4772 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4773 cell_x0, cell_x1, limitd, limit0, limit1,
4775 ( thread *dd->ncg_home)/nthread,
4776 ((thread+1)*dd->ncg_home)/nthread,
4777 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4781 for (cg = 0; cg < dd->ncg_home; cg++)
4786 flag = mc & ~DD_FLAG_NRCG;
4787 mc = mc & DD_FLAG_NRCG;
4790 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4792 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4793 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4795 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4796 /* We store the cg size in the lower 16 bits
4797 * and the place where the charge group should go
4798 * in the next 6 bits. This saves some communication volume.
4800 nrcg = cgindex[cg+1] - cgindex[cg];
4801 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4807 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4808 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4811 for (i = 0; i < dd->ndim*2; i++)
4813 *ncg_moved += ncg[i];
4830 /* Make sure the communication buffers are large enough */
4831 for (mc = 0; mc < dd->ndim*2; mc++)
4833 nvr = ncg[mc] + nat[mc]*nvec;
4834 if (nvr > comm->cgcm_state_nalloc[mc])
4836 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4837 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4841 switch (fr->cutoff_scheme)
4844 /* Recalculating cg_cm might be cheaper than communicating,
4845 * but that could give rise to rounding issues.
4848 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4849 nvec, cg_cm, comm, bCompact);
4852 /* Without charge groups we send the moved atom coordinates
4853 * over twice. This is so the code below can be used without
4854 * many conditionals for both for with and without charge groups.
4857 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4858 nvec, state->x, comm, FALSE);
4861 home_pos_cg -= *ncg_moved;
4865 gmx_incons("unimplemented");
4871 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4872 nvec, vec++, state->x, comm, bCompact);
4875 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4876 nvec, vec++, state->v, comm, bCompact);
4880 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4881 nvec, vec++, state->sd_X, comm, bCompact);
4885 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4886 nvec, vec++, state->cg_p, comm, bCompact);
4891 compact_ind(dd->ncg_home, move,
4892 dd->index_gl, dd->cgindex, dd->gatindex,
4893 dd->ga2la, comm->bLocalCG,
4898 if (fr->cutoff_scheme == ecutsVERLET)
4900 moved = get_moved(comm, dd->ncg_home);
4902 for (k = 0; k < dd->ncg_home; k++)
4909 moved = fr->ns.grid->cell_index;
4912 clear_and_mark_ind(dd->ncg_home, move,
4913 dd->index_gl, dd->cgindex, dd->gatindex,
4914 dd->ga2la, comm->bLocalCG,
4918 cginfo_mb = fr->cginfo_mb;
4920 *ncg_stay_home = home_pos_cg;
4921 for (d = 0; d < dd->ndim; d++)
4926 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4929 /* Communicate the cg and atom counts */
4934 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4935 d, dir, sbuf[0], sbuf[1]);
4937 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4939 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4941 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4942 srenew(comm->buf_int, comm->nalloc_int);
4945 /* Communicate the charge group indices, sizes and flags */
4946 dd_sendrecv_int(dd, d, dir,
4947 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4948 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4950 nvs = ncg[cdd] + nat[cdd]*nvec;
4951 i = rbuf[0] + rbuf[1] *nvec;
4952 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4954 /* Communicate cgcm and state */
4955 dd_sendrecv_rvec(dd, d, dir,
4956 comm->cgcm_state[cdd], nvs,
4957 comm->vbuf.v+nvr, i);
4958 ncg_recv += rbuf[0];
4962 /* Process the received charge groups */
4964 for (cg = 0; cg < ncg_recv; cg++)
4966 flag = comm->buf_int[cg*DD_CGIBS+1];
4968 if (dim >= npbcdim && dd->nc[dim] > 2)
4970 /* No pbc in this dim and more than one domain boundary.
4971 * We do a separate check if a charge group didn't move too far.
4973 if (((flag & DD_FLAG_FW(d)) &&
4974 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4975 ((flag & DD_FLAG_BW(d)) &&
4976 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4978 cg_move_error(fplog, dd, step, cg, dim,
4979 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4980 fr->cutoff_scheme == ecutsGROUP, 0,
4981 comm->vbuf.v[buf_pos],
4982 comm->vbuf.v[buf_pos],
4983 comm->vbuf.v[buf_pos][dim]);
4990 /* Check which direction this cg should go */
4991 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4995 /* The cell boundaries for dimension d2 are not equal
4996 * for each cell row of the lower dimension(s),
4997 * therefore we might need to redetermine where
4998 * this cg should go.
5001 /* If this cg crosses the box boundary in dimension d2
5002 * we can use the communicated flag, so we do not
5003 * have to worry about pbc.
5005 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5006 (flag & DD_FLAG_FW(d2))) ||
5007 (dd->ci[dim2] == 0 &&
5008 (flag & DD_FLAG_BW(d2)))))
5010 /* Clear the two flags for this dimension */
5011 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5012 /* Determine the location of this cg
5013 * in lattice coordinates
5015 pos_d = comm->vbuf.v[buf_pos][dim2];
5018 for (d3 = dim2+1; d3 < DIM; d3++)
5021 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5024 /* Check of we are not at the box edge.
5025 * pbc is only handled in the first step above,
5026 * but this check could move over pbc while
5027 * the first step did not due to different rounding.
5029 if (pos_d >= cell_x1[dim2] &&
5030 dd->ci[dim2] != dd->nc[dim2]-1)
5032 flag |= DD_FLAG_FW(d2);
5034 else if (pos_d < cell_x0[dim2] &&
5037 flag |= DD_FLAG_BW(d2);
5039 comm->buf_int[cg*DD_CGIBS+1] = flag;
5042 /* Set to which neighboring cell this cg should go */
5043 if (flag & DD_FLAG_FW(d2))
5047 else if (flag & DD_FLAG_BW(d2))
5049 if (dd->nc[dd->dim[d2]] > 2)
5061 nrcg = flag & DD_FLAG_NRCG;
5064 if (home_pos_cg+1 > dd->cg_nalloc)
5066 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5067 srenew(dd->index_gl, dd->cg_nalloc);
5068 srenew(dd->cgindex, dd->cg_nalloc+1);
5070 /* Set the global charge group index and size */
5071 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5072 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5073 /* Copy the state from the buffer */
5074 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5075 if (fr->cutoff_scheme == ecutsGROUP)
5078 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5082 /* Set the cginfo */
5083 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5084 dd->index_gl[home_pos_cg]);
5087 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5090 if (home_pos_at+nrcg > state->nalloc)
5092 dd_realloc_state(state, f, home_pos_at+nrcg);
5094 for (i = 0; i < nrcg; i++)
5096 copy_rvec(comm->vbuf.v[buf_pos++],
5097 state->x[home_pos_at+i]);
5101 for (i = 0; i < nrcg; i++)
5103 copy_rvec(comm->vbuf.v[buf_pos++],
5104 state->v[home_pos_at+i]);
5109 for (i = 0; i < nrcg; i++)
5111 copy_rvec(comm->vbuf.v[buf_pos++],
5112 state->sd_X[home_pos_at+i]);
5117 for (i = 0; i < nrcg; i++)
5119 copy_rvec(comm->vbuf.v[buf_pos++],
5120 state->cg_p[home_pos_at+i]);
5124 home_pos_at += nrcg;
5128 /* Reallocate the buffers if necessary */
5129 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5131 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5132 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5134 nvr = ncg[mc] + nat[mc]*nvec;
5135 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5137 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5138 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5140 /* Copy from the receive to the send buffers */
5141 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5142 comm->buf_int + cg*DD_CGIBS,
5143 DD_CGIBS*sizeof(int));
5144 memcpy(comm->cgcm_state[mc][nvr],
5145 comm->vbuf.v[buf_pos],
5146 (1+nrcg*nvec)*sizeof(rvec));
5147 buf_pos += 1 + nrcg*nvec;
5154 /* With sorting (!bCompact) the indices are now only partially up to date
5155 * and ncg_home and nat_home are not the real count, since there are
5156 * "holes" in the arrays for the charge groups that moved to neighbors.
5158 if (fr->cutoff_scheme == ecutsVERLET)
5160 moved = get_moved(comm, home_pos_cg);
5162 for (i = dd->ncg_home; i < home_pos_cg; i++)
5167 dd->ncg_home = home_pos_cg;
5168 dd->nat_home = home_pos_at;
5173 "Finished repartitioning: cgs moved out %d, new home %d\n",
5174 *ncg_moved, dd->ncg_home-*ncg_moved);
5179 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5181 dd->comm->cycl[ddCycl] += cycles;
5182 dd->comm->cycl_n[ddCycl]++;
5183 if (cycles > dd->comm->cycl_max[ddCycl])
5185 dd->comm->cycl_max[ddCycl] = cycles;
5189 static double force_flop_count(t_nrnb *nrnb)
5196 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5198 /* To get closer to the real timings, we half the count
5199 * for the normal loops and again half it for water loops.
5202 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5204 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5208 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5211 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5214 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5216 sum += nrnb->n[i]*cost_nrnb(i);
5219 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5221 sum += nrnb->n[i]*cost_nrnb(i);
5227 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5229 if (dd->comm->eFlop)
5231 dd->comm->flop -= force_flop_count(nrnb);
5234 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5236 if (dd->comm->eFlop)
5238 dd->comm->flop += force_flop_count(nrnb);
5243 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5247 for (i = 0; i < ddCyclNr; i++)
5249 dd->comm->cycl[i] = 0;
5250 dd->comm->cycl_n[i] = 0;
5251 dd->comm->cycl_max[i] = 0;
5254 dd->comm->flop_n = 0;
5257 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5259 gmx_domdec_comm_t *comm;
5260 gmx_domdec_load_t *load;
5261 gmx_domdec_root_t *root = NULL;
5263 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5268 fprintf(debug, "get_load_distribution start\n");
5271 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5275 bSepPME = (dd->pme_nodeid >= 0);
5277 for (d = dd->ndim-1; d >= 0; d--)
5280 /* Check if we participate in the communication in this dimension */
5281 if (d == dd->ndim-1 ||
5282 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5284 load = &comm->load[d];
5287 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5290 if (d == dd->ndim-1)
5292 sbuf[pos++] = dd_force_load(comm);
5293 sbuf[pos++] = sbuf[0];
5296 sbuf[pos++] = sbuf[0];
5297 sbuf[pos++] = cell_frac;
5300 sbuf[pos++] = comm->cell_f_max0[d];
5301 sbuf[pos++] = comm->cell_f_min1[d];
5306 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5307 sbuf[pos++] = comm->cycl[ddCyclPME];
5312 sbuf[pos++] = comm->load[d+1].sum;
5313 sbuf[pos++] = comm->load[d+1].max;
5316 sbuf[pos++] = comm->load[d+1].sum_m;
5317 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5318 sbuf[pos++] = comm->load[d+1].flags;
5321 sbuf[pos++] = comm->cell_f_max0[d];
5322 sbuf[pos++] = comm->cell_f_min1[d];
5327 sbuf[pos++] = comm->load[d+1].mdf;
5328 sbuf[pos++] = comm->load[d+1].pme;
5332 /* Communicate a row in DD direction d.
5333 * The communicators are setup such that the root always has rank 0.
5336 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5337 load->load, load->nload*sizeof(float), MPI_BYTE,
5338 0, comm->mpi_comm_load[d]);
5340 if (dd->ci[dim] == dd->master_ci[dim])
5342 /* We are the root, process this row */
5343 if (comm->bDynLoadBal)
5345 root = comm->root[d];
5355 for (i = 0; i < dd->nc[dim]; i++)
5357 load->sum += load->load[pos++];
5358 load->max = std::max(load->max, load->load[pos]);
5364 /* This direction could not be load balanced properly,
5365 * therefore we need to use the maximum iso the average load.
5367 load->sum_m = std::max(load->sum_m, load->load[pos]);
5371 load->sum_m += load->load[pos];
5374 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5378 load->flags = (int)(load->load[pos++] + 0.5);
5382 root->cell_f_max0[i] = load->load[pos++];
5383 root->cell_f_min1[i] = load->load[pos++];
5388 load->mdf = std::max(load->mdf, load->load[pos]);
5390 load->pme = std::max(load->pme, load->load[pos]);
5394 if (comm->bDynLoadBal && root->bLimited)
5396 load->sum_m *= dd->nc[dim];
5397 load->flags |= (1<<d);
5405 comm->nload += dd_load_count(comm);
5406 comm->load_step += comm->cycl[ddCyclStep];
5407 comm->load_sum += comm->load[0].sum;
5408 comm->load_max += comm->load[0].max;
5409 if (comm->bDynLoadBal)
5411 for (d = 0; d < dd->ndim; d++)
5413 if (comm->load[0].flags & (1<<d))
5415 comm->load_lim[d]++;
5421 comm->load_mdf += comm->load[0].mdf;
5422 comm->load_pme += comm->load[0].pme;
5426 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5430 fprintf(debug, "get_load_distribution finished\n");
5434 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5436 /* Return the relative performance loss on the total run time
5437 * due to the force calculation load imbalance.
5439 if (dd->comm->nload > 0)
5442 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5443 (dd->comm->load_step*dd->nnodes);
5451 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5454 int npp, npme, nnodes, d, limp;
5455 float imbal, pme_f_ratio, lossf, lossp = 0;
5457 gmx_domdec_comm_t *comm;
5460 if (DDMASTER(dd) && comm->nload > 0)
5463 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5464 nnodes = npp + npme;
5465 imbal = comm->load_max*npp/comm->load_sum - 1;
5466 lossf = dd_force_imb_perf_loss(dd);
5467 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5468 fprintf(fplog, "%s", buf);
5469 fprintf(stderr, "\n");
5470 fprintf(stderr, "%s", buf);
5471 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5472 fprintf(fplog, "%s", buf);
5473 fprintf(stderr, "%s", buf);
5475 if (comm->bDynLoadBal)
5477 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5478 for (d = 0; d < dd->ndim; d++)
5480 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5481 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5487 sprintf(buf+strlen(buf), "\n");
5488 fprintf(fplog, "%s", buf);
5489 fprintf(stderr, "%s", buf);
5493 pme_f_ratio = comm->load_pme/comm->load_mdf;
5494 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5497 lossp *= (float)npme/(float)nnodes;
5501 lossp *= (float)npp/(float)nnodes;
5503 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5504 fprintf(fplog, "%s", buf);
5505 fprintf(stderr, "%s", buf);
5506 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5507 fprintf(fplog, "%s", buf);
5508 fprintf(stderr, "%s", buf);
5510 fprintf(fplog, "\n");
5511 fprintf(stderr, "\n");
5513 if (lossf >= DD_PERF_LOSS_WARN)
5516 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5517 " in the domain decomposition.\n", lossf*100);
5518 if (!comm->bDynLoadBal)
5520 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5524 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5526 fprintf(fplog, "%s\n", buf);
5527 fprintf(stderr, "%s\n", buf);
5529 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5532 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5533 " had %s work to do than the PP ranks.\n"
5534 " You might want to %s the number of PME ranks\n"
5535 " or %s the cut-off and the grid spacing.\n",
5537 (lossp < 0) ? "less" : "more",
5538 (lossp < 0) ? "decrease" : "increase",
5539 (lossp < 0) ? "decrease" : "increase");
5540 fprintf(fplog, "%s\n", buf);
5541 fprintf(stderr, "%s\n", buf);
5546 static float dd_vol_min(gmx_domdec_t *dd)
5548 return dd->comm->load[0].cvol_min*dd->nnodes;
5551 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5553 return dd->comm->load[0].flags;
5556 static float dd_f_imbal(gmx_domdec_t *dd)
5558 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5561 float dd_pme_f_ratio(gmx_domdec_t *dd)
5563 if (dd->comm->cycl_n[ddCyclPME] > 0)
5565 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5573 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5578 flags = dd_load_flags(dd);
5582 "DD load balancing is limited by minimum cell size in dimension");
5583 for (d = 0; d < dd->ndim; d++)
5587 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5590 fprintf(fplog, "\n");
5592 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5593 if (dd->comm->bDynLoadBal)
5595 fprintf(fplog, " vol min/aver %5.3f%c",
5596 dd_vol_min(dd), flags ? '!' : ' ');
5598 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5599 if (dd->comm->cycl_n[ddCyclPME])
5601 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5603 fprintf(fplog, "\n\n");
5606 static void dd_print_load_verbose(gmx_domdec_t *dd)
5608 if (dd->comm->bDynLoadBal)
5610 fprintf(stderr, "vol %4.2f%c ",
5611 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5613 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5614 if (dd->comm->cycl_n[ddCyclPME])
5616 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5621 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5626 gmx_domdec_root_t *root;
5627 gmx_bool bPartOfGroup = FALSE;
5629 dim = dd->dim[dim_ind];
5630 copy_ivec(loc, loc_c);
5631 for (i = 0; i < dd->nc[dim]; i++)
5634 rank = dd_index(dd->nc, loc_c);
5635 if (rank == dd->rank)
5637 /* This process is part of the group */
5638 bPartOfGroup = TRUE;
5641 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5645 dd->comm->mpi_comm_load[dim_ind] = c_row;
5646 if (dd->comm->eDLB != edlbNO)
5648 if (dd->ci[dim] == dd->master_ci[dim])
5650 /* This is the root process of this row */
5651 snew(dd->comm->root[dim_ind], 1);
5652 root = dd->comm->root[dim_ind];
5653 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5654 snew(root->old_cell_f, dd->nc[dim]+1);
5655 snew(root->bCellMin, dd->nc[dim]);
5658 snew(root->cell_f_max0, dd->nc[dim]);
5659 snew(root->cell_f_min1, dd->nc[dim]);
5660 snew(root->bound_min, dd->nc[dim]);
5661 snew(root->bound_max, dd->nc[dim]);
5663 snew(root->buf_ncd, dd->nc[dim]);
5667 /* This is not a root process, we only need to receive cell_f */
5668 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5671 if (dd->ci[dim] == dd->master_ci[dim])
5673 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5679 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5680 const gmx_hw_info_t gmx_unused *hwinfo,
5681 const gmx_hw_opt_t gmx_unused *hw_opt)
5684 int physicalnode_id_hash;
5687 MPI_Comm mpi_comm_pp_physicalnode;
5689 if (!(cr->duty & DUTY_PP) ||
5690 hw_opt->gpu_opt.ncuda_dev_use == 0)
5692 /* Only PP nodes (currently) use GPUs.
5693 * If we don't have GPUs, there are no resources to share.
5698 physicalnode_id_hash = gmx_physicalnode_id_hash();
5700 gpu_id = get_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5706 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5707 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5708 dd->rank, physicalnode_id_hash, gpu_id);
5710 /* Split the PP communicator over the physical nodes */
5711 /* TODO: See if we should store this (before), as it's also used for
5712 * for the nodecomm summution.
5714 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5715 &mpi_comm_pp_physicalnode);
5716 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5717 &dd->comm->mpi_comm_gpu_shared);
5718 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5719 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5723 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5726 /* Note that some ranks could share a GPU, while others don't */
5728 if (dd->comm->nrank_gpu_shared == 1)
5730 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5735 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5738 int dim0, dim1, i, j;
5743 fprintf(debug, "Making load communicators\n");
5746 snew(dd->comm->load, dd->ndim);
5747 snew(dd->comm->mpi_comm_load, dd->ndim);
5750 make_load_communicator(dd, 0, loc);
5754 for (i = 0; i < dd->nc[dim0]; i++)
5757 make_load_communicator(dd, 1, loc);
5763 for (i = 0; i < dd->nc[dim0]; i++)
5767 for (j = 0; j < dd->nc[dim1]; j++)
5770 make_load_communicator(dd, 2, loc);
5777 fprintf(debug, "Finished making load communicators\n");
5782 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5784 int d, dim, i, j, m;
5787 ivec dd_zp[DD_MAXIZONE];
5788 gmx_domdec_zones_t *zones;
5789 gmx_domdec_ns_ranges_t *izone;
5791 for (d = 0; d < dd->ndim; d++)
5794 copy_ivec(dd->ci, tmp);
5795 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5796 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5797 copy_ivec(dd->ci, tmp);
5798 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5799 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5802 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5805 dd->neighbor[d][1]);
5811 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5813 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5814 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5821 for (i = 0; i < nzonep; i++)
5823 copy_ivec(dd_zp3[i], dd_zp[i]);
5829 for (i = 0; i < nzonep; i++)
5831 copy_ivec(dd_zp2[i], dd_zp[i]);
5837 for (i = 0; i < nzonep; i++)
5839 copy_ivec(dd_zp1[i], dd_zp[i]);
5843 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5848 zones = &dd->comm->zones;
5850 for (i = 0; i < nzone; i++)
5853 clear_ivec(zones->shift[i]);
5854 for (d = 0; d < dd->ndim; d++)
5856 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5861 for (i = 0; i < nzone; i++)
5863 for (d = 0; d < DIM; d++)
5865 s[d] = dd->ci[d] - zones->shift[i][d];
5870 else if (s[d] >= dd->nc[d])
5876 zones->nizone = nzonep;
5877 for (i = 0; i < zones->nizone; i++)
5879 if (dd_zp[i][0] != i)
5881 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5883 izone = &zones->izone[i];
5884 izone->j0 = dd_zp[i][1];
5885 izone->j1 = dd_zp[i][2];
5886 for (dim = 0; dim < DIM; dim++)
5888 if (dd->nc[dim] == 1)
5890 /* All shifts should be allowed */
5891 izone->shift0[dim] = -1;
5892 izone->shift1[dim] = 1;
5897 izone->shift0[d] = 0;
5898 izone->shift1[d] = 0;
5899 for(j=izone->j0; j<izone->j1; j++) {
5900 if (dd->shift[j][d] > dd->shift[i][d])
5901 izone->shift0[d] = -1;
5902 if (dd->shift[j][d] < dd->shift[i][d])
5903 izone->shift1[d] = 1;
5909 /* Assume the shift are not more than 1 cell */
5910 izone->shift0[dim] = 1;
5911 izone->shift1[dim] = -1;
5912 for (j = izone->j0; j < izone->j1; j++)
5914 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5915 if (shift_diff < izone->shift0[dim])
5917 izone->shift0[dim] = shift_diff;
5919 if (shift_diff > izone->shift1[dim])
5921 izone->shift1[dim] = shift_diff;
5928 if (dd->comm->eDLB != edlbNO)
5930 snew(dd->comm->root, dd->ndim);
5933 if (dd->comm->bRecordLoad)
5935 make_load_communicators(dd);
5939 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5945 gmx_domdec_comm_t *comm;
5952 if (comm->bCartesianPP)
5954 /* Set up cartesian communication for the particle-particle part */
5957 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5958 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5961 for (int i = 0; i < DIM; i++)
5965 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5967 /* We overwrite the old communicator with the new cartesian one */
5968 cr->mpi_comm_mygroup = comm_cart;
5971 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5972 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5974 if (comm->bCartesianPP_PME)
5976 /* Since we want to use the original cartesian setup for sim,
5977 * and not the one after split, we need to make an index.
5979 snew(comm->ddindex2ddnodeid, dd->nnodes);
5980 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5981 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5982 /* Get the rank of the DD master,
5983 * above we made sure that the master node is a PP node.
5993 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5995 else if (comm->bCartesianPP)
5997 if (cr->npmenodes == 0)
5999 /* The PP communicator is also
6000 * the communicator for this simulation
6002 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6004 cr->nodeid = dd->rank;
6006 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6008 /* We need to make an index to go from the coordinates
6009 * to the nodeid of this simulation.
6011 snew(comm->ddindex2simnodeid, dd->nnodes);
6012 snew(buf, dd->nnodes);
6013 if (cr->duty & DUTY_PP)
6015 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6017 /* Communicate the ddindex to simulation nodeid index */
6018 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6019 cr->mpi_comm_mysim);
6022 /* Determine the master coordinates and rank.
6023 * The DD master should be the same node as the master of this sim.
6025 for (int i = 0; i < dd->nnodes; i++)
6027 if (comm->ddindex2simnodeid[i] == 0)
6029 ddindex2xyz(dd->nc, i, dd->master_ci);
6030 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6035 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6040 /* No Cartesian communicators */
6041 /* We use the rank in dd->comm->all as DD index */
6042 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6043 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6045 clear_ivec(dd->master_ci);
6052 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6053 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6058 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6059 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6063 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6067 gmx_domdec_comm_t *comm;
6072 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6075 snew(comm->ddindex2simnodeid, dd->nnodes);
6076 snew(buf, dd->nnodes);
6077 if (cr->duty & DUTY_PP)
6079 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6081 /* Communicate the ddindex to simulation nodeid index */
6082 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6083 cr->mpi_comm_mysim);
6089 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6090 int ncg, int natoms)
6092 gmx_domdec_master_t *ma;
6097 snew(ma->ncg, dd->nnodes);
6098 snew(ma->index, dd->nnodes+1);
6100 snew(ma->nat, dd->nnodes);
6101 snew(ma->ibuf, dd->nnodes*2);
6102 snew(ma->cell_x, DIM);
6103 for (i = 0; i < DIM; i++)
6105 snew(ma->cell_x[i], dd->nc[i]+1);
6108 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6114 snew(ma->vbuf, natoms);
6120 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6121 int gmx_unused reorder)
6124 gmx_domdec_comm_t *comm;
6134 if (comm->bCartesianPP)
6136 for (i = 1; i < DIM; i++)
6138 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6140 if (bDiv[YY] || bDiv[ZZ])
6142 comm->bCartesianPP_PME = TRUE;
6143 /* If we have 2D PME decomposition, which is always in x+y,
6144 * we stack the PME only nodes in z.
6145 * Otherwise we choose the direction that provides the thinnest slab
6146 * of PME only nodes as this will have the least effect
6147 * on the PP communication.
6148 * But for the PME communication the opposite might be better.
6150 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6152 dd->nc[YY] > dd->nc[ZZ]))
6154 comm->cartpmedim = ZZ;
6158 comm->cartpmedim = YY;
6160 comm->ntot[comm->cartpmedim]
6161 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6165 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]);
6167 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6172 if (comm->bCartesianPP_PME)
6179 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]);
6182 for (i = 0; i < DIM; i++)
6186 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6188 MPI_Comm_rank(comm_cart, &rank);
6189 if (MASTERNODE(cr) && rank != 0)
6191 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6194 /* With this assigment we loose the link to the original communicator
6195 * which will usually be MPI_COMM_WORLD, unless have multisim.
6197 cr->mpi_comm_mysim = comm_cart;
6198 cr->sim_nodeid = rank;
6200 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6204 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6205 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6208 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6212 if (cr->npmenodes == 0 ||
6213 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6215 cr->duty = DUTY_PME;
6218 /* Split the sim communicator into PP and PME only nodes */
6219 MPI_Comm_split(cr->mpi_comm_mysim,
6221 dd_index(comm->ntot, dd->ci),
6222 &cr->mpi_comm_mygroup);
6226 switch (dd_node_order)
6231 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6234 case ddnoINTERLEAVE:
6235 /* Interleave the PP-only and PME-only nodes,
6236 * as on clusters with dual-core machines this will double
6237 * the communication bandwidth of the PME processes
6238 * and thus speed up the PP <-> PME and inter PME communication.
6242 fprintf(fplog, "Interleaving PP and PME ranks\n");
6244 comm->pmenodes = dd_pmenodes(cr);
6249 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6252 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6254 cr->duty = DUTY_PME;
6261 /* Split the sim communicator into PP and PME only nodes */
6262 MPI_Comm_split(cr->mpi_comm_mysim,
6265 &cr->mpi_comm_mygroup);
6266 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6272 fprintf(fplog, "This rank does only %s work.\n\n",
6273 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6277 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6280 gmx_domdec_comm_t *comm;
6286 copy_ivec(dd->nc, comm->ntot);
6288 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6289 comm->bCartesianPP_PME = FALSE;
6291 /* Reorder the nodes by default. This might change the MPI ranks.
6292 * Real reordering is only supported on very few architectures,
6293 * Blue Gene is one of them.
6295 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6297 if (cr->npmenodes > 0)
6299 /* Split the communicator into a PP and PME part */
6300 split_communicator(fplog, cr, dd_node_order, CartReorder);
6301 if (comm->bCartesianPP_PME)
6303 /* We (possibly) reordered the nodes in split_communicator,
6304 * so it is no longer required in make_pp_communicator.
6306 CartReorder = FALSE;
6311 /* All nodes do PP and PME */
6313 /* We do not require separate communicators */
6314 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6318 if (cr->duty & DUTY_PP)
6320 /* Copy or make a new PP communicator */
6321 make_pp_communicator(fplog, cr, CartReorder);
6325 receive_ddindex2simnodeid(cr);
6328 if (!(cr->duty & DUTY_PME))
6330 /* Set up the commnuication to our PME node */
6331 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6332 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6335 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6336 dd->pme_nodeid, dd->pme_receive_vir_ener);
6341 dd->pme_nodeid = -1;
6346 dd->ma = init_gmx_domdec_master_t(dd,
6348 comm->cgs_gl.index[comm->cgs_gl.nr]);
6352 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6354 real *slb_frac, tot;
6359 if (nc > 1 && size_string != NULL)
6363 fprintf(fplog, "Using static load balancing for the %s direction\n",
6368 for (i = 0; i < nc; i++)
6371 sscanf(size_string, "%20lf%n", &dbl, &n);
6374 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6383 fprintf(fplog, "Relative cell sizes:");
6385 for (i = 0; i < nc; i++)
6390 fprintf(fplog, " %5.3f", slb_frac[i]);
6395 fprintf(fplog, "\n");
6402 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6405 gmx_mtop_ilistloop_t iloop;
6409 iloop = gmx_mtop_ilistloop_init(mtop);
6410 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6412 for (ftype = 0; ftype < F_NRE; ftype++)
6414 if ((interaction_function[ftype].flags & IF_BOND) &&
6417 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6425 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6431 val = getenv(env_var);
6434 if (sscanf(val, "%20d", &nst) <= 0)
6440 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6448 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6452 fprintf(stderr, "\n%s\n", warn_string);
6456 fprintf(fplog, "\n%s\n", warn_string);
6460 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6461 t_inputrec *ir, FILE *fplog)
6463 if (ir->ePBC == epbcSCREW &&
6464 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6466 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6469 if (ir->ns_type == ensSIMPLE)
6471 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6474 if (ir->nstlist == 0)
6476 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6479 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6481 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6485 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6490 r = ddbox->box_size[XX];
6491 for (di = 0; di < dd->ndim; di++)
6494 /* Check using the initial average cell size */
6495 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6501 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6502 const char *dlb_opt, gmx_bool bRecordLoad,
6503 unsigned long Flags, t_inputrec *ir)
6510 case 'a': eDLB = edlbAUTO; break;
6511 case 'n': eDLB = edlbNO; break;
6512 case 'y': eDLB = edlbYES; break;
6513 default: gmx_incons("Unknown dlb_opt");
6516 if (Flags & MD_RERUN)
6521 if (!EI_DYNAMICS(ir->eI))
6523 if (eDLB == edlbYES)
6525 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6526 dd_warning(cr, fplog, buf);
6534 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6539 if (Flags & MD_REPRODUCIBLE)
6546 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6550 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6553 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6561 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6566 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6568 /* Decomposition order z,y,x */
6571 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6573 for (dim = DIM-1; dim >= 0; dim--)
6575 if (dd->nc[dim] > 1)
6577 dd->dim[dd->ndim++] = dim;
6583 /* Decomposition order x,y,z */
6584 for (dim = 0; dim < DIM; dim++)
6586 if (dd->nc[dim] > 1)
6588 dd->dim[dd->ndim++] = dim;
6594 static gmx_domdec_comm_t *init_dd_comm()
6596 gmx_domdec_comm_t *comm;
6600 snew(comm->cggl_flag, DIM*2);
6601 snew(comm->cgcm_state, DIM*2);
6602 for (i = 0; i < DIM*2; i++)
6604 comm->cggl_flag_nalloc[i] = 0;
6605 comm->cgcm_state_nalloc[i] = 0;
6608 comm->nalloc_int = 0;
6609 comm->buf_int = NULL;
6611 vec_rvec_init(&comm->vbuf);
6613 comm->n_load_have = 0;
6614 comm->n_load_collect = 0;
6616 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6618 comm->sum_nat[i] = 0;
6622 comm->load_step = 0;
6625 clear_ivec(comm->load_lim);
6632 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6633 unsigned long Flags,
6635 real comm_distance_min, real rconstr,
6636 const char *dlb_opt, real dlb_scale,
6637 const char *sizex, const char *sizey, const char *sizez,
6638 gmx_mtop_t *mtop, t_inputrec *ir,
6639 matrix box, rvec *x,
6641 int *npme_x, int *npme_y)
6644 gmx_domdec_comm_t *comm;
6646 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6649 const real tenPercentMargin = 1.1;
6654 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6659 dd->comm = init_dd_comm();
6661 snew(comm->cggl_flag, DIM*2);
6662 snew(comm->cgcm_state, DIM*2);
6664 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6665 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6667 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6668 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6669 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6670 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6671 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6672 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6673 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6674 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6676 dd->pme_recv_f_alloc = 0;
6677 dd->pme_recv_f_buf = NULL;
6679 if (dd->bSendRecv2 && fplog)
6681 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");
6687 fprintf(fplog, "Will load balance based on FLOP count\n");
6689 if (comm->eFlop > 1)
6691 srand(1+cr->nodeid);
6693 comm->bRecordLoad = TRUE;
6697 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6701 /* Initialize to GPU share count to 0, might change later */
6702 comm->nrank_gpu_shared = 0;
6704 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6705 comm->bDLB_locked = FALSE;
6707 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6710 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6712 dd->bGridJump = comm->bDynLoadBal;
6713 comm->bPMELoadBalDLBLimits = FALSE;
6715 if (comm->nstSortCG)
6719 if (comm->nstSortCG == 1)
6721 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6725 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6729 snew(comm->sort, 1);
6735 fprintf(fplog, "Will not sort the charge groups\n");
6739 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6741 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6742 if (comm->bInterCGBondeds)
6744 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6748 comm->bInterCGMultiBody = FALSE;
6751 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6752 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6754 if (ir->rlistlong == 0)
6756 /* Set the cut-off to some very large value,
6757 * so we don't need if statements everywhere in the code.
6758 * We use sqrt, since the cut-off is squared in some places.
6760 comm->cutoff = GMX_CUTOFF_INF;
6764 comm->cutoff = ir->rlistlong;
6766 comm->cutoff_mbody = 0;
6768 comm->cellsize_limit = 0;
6769 comm->bBondComm = FALSE;
6771 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6772 * within nstlist steps. Since boundaries are allowed to displace by half
6773 * a cell size, DD cells should be at least the size of the list buffer.
6775 comm->cellsize_limit = std::max(comm->cellsize_limit,
6776 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6778 if (comm->bInterCGBondeds)
6780 if (comm_distance_min > 0)
6782 comm->cutoff_mbody = comm_distance_min;
6783 if (Flags & MD_DDBONDCOMM)
6785 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6789 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6791 r_bonded_limit = comm->cutoff_mbody;
6793 else if (ir->bPeriodicMols)
6795 /* Can not easily determine the required cut-off */
6796 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");
6797 comm->cutoff_mbody = comm->cutoff/2;
6798 r_bonded_limit = comm->cutoff_mbody;
6804 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6805 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6807 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6808 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6810 /* We use an initial margin of 10% for the minimum cell size,
6811 * except when we are just below the non-bonded cut-off.
6813 if (Flags & MD_DDBONDCOMM)
6815 if (std::max(r_2b, r_mb) > comm->cutoff)
6817 r_bonded = std::max(r_2b, r_mb);
6818 r_bonded_limit = tenPercentMargin*r_bonded;
6819 comm->bBondComm = TRUE;
6824 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6826 /* We determine cutoff_mbody later */
6830 /* No special bonded communication,
6831 * simply increase the DD cut-off.
6833 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6834 comm->cutoff_mbody = r_bonded_limit;
6835 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6841 "Minimum cell size due to bonded interactions: %.3f nm\n",
6844 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6847 if (dd->bInterCGcons && rconstr <= 0)
6849 /* There is a cell size limit due to the constraints (P-LINCS) */
6850 rconstr = constr_r_max(fplog, mtop, ir);
6854 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6856 if (rconstr > comm->cellsize_limit)
6858 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6862 else if (rconstr > 0 && fplog)
6864 /* Here we do not check for dd->bInterCGcons,
6865 * because one can also set a cell size limit for virtual sites only
6866 * and at this point we don't know yet if there are intercg v-sites.
6869 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6872 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6874 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6878 copy_ivec(nc, dd->nc);
6879 set_dd_dim(fplog, dd);
6880 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6882 if (cr->npmenodes == -1)
6886 acs = average_cellsize_min(dd, ddbox);
6887 if (acs < comm->cellsize_limit)
6891 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6893 gmx_fatal_collective(FARGS, cr, NULL,
6894 "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",
6895 acs, comm->cellsize_limit);
6900 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6902 /* We need to choose the optimal DD grid and possibly PME nodes */
6903 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6904 comm->eDLB != edlbNO, dlb_scale,
6905 comm->cellsize_limit, comm->cutoff,
6906 comm->bInterCGBondeds);
6908 if (dd->nc[XX] == 0)
6910 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6911 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6912 !bC ? "-rdd" : "-rcon",
6913 comm->eDLB != edlbNO ? " or -dds" : "",
6914 bC ? " or your LINCS settings" : "");
6916 gmx_fatal_collective(FARGS, cr, NULL,
6917 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6919 "Look in the log file for details on the domain decomposition",
6920 cr->nnodes-cr->npmenodes, limit, buf);
6922 set_dd_dim(fplog, dd);
6928 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6929 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6932 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6933 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6935 gmx_fatal_collective(FARGS, cr, NULL,
6936 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6937 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6939 if (cr->npmenodes > dd->nnodes)
6941 gmx_fatal_collective(FARGS, cr, NULL,
6942 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6944 if (cr->npmenodes > 0)
6946 comm->npmenodes = cr->npmenodes;
6950 comm->npmenodes = dd->nnodes;
6953 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6955 /* The following choices should match those
6956 * in comm_cost_est in domdec_setup.c.
6957 * Note that here the checks have to take into account
6958 * that the decomposition might occur in a different order than xyz
6959 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6960 * in which case they will not match those in comm_cost_est,
6961 * but since that is mainly for testing purposes that's fine.
6963 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6964 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6965 getenv("GMX_PMEONEDD") == NULL)
6967 comm->npmedecompdim = 2;
6968 comm->npmenodes_x = dd->nc[XX];
6969 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6973 /* In case nc is 1 in both x and y we could still choose to
6974 * decompose pme in y instead of x, but we use x for simplicity.
6976 comm->npmedecompdim = 1;
6977 if (dd->dim[0] == YY)
6979 comm->npmenodes_x = 1;
6980 comm->npmenodes_y = comm->npmenodes;
6984 comm->npmenodes_x = comm->npmenodes;
6985 comm->npmenodes_y = 1;
6990 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6991 comm->npmenodes_x, comm->npmenodes_y, 1);
6996 comm->npmedecompdim = 0;
6997 comm->npmenodes_x = 0;
6998 comm->npmenodes_y = 0;
7001 /* Technically we don't need both of these,
7002 * but it simplifies code not having to recalculate it.
7004 *npme_x = comm->npmenodes_x;
7005 *npme_y = comm->npmenodes_y;
7007 snew(comm->slb_frac, DIM);
7008 if (comm->eDLB == edlbNO)
7010 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7011 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7012 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7015 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7017 if (comm->bBondComm || comm->eDLB != edlbNO)
7019 /* Set the bonded communication distance to halfway
7020 * the minimum and the maximum,
7021 * since the extra communication cost is nearly zero.
7023 acs = average_cellsize_min(dd, ddbox);
7024 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7025 if (comm->eDLB != edlbNO)
7027 /* Check if this does not limit the scaling */
7028 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7030 if (!comm->bBondComm)
7032 /* Without bBondComm do not go beyond the n.b. cut-off */
7033 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7034 if (comm->cellsize_limit >= comm->cutoff)
7036 /* We don't loose a lot of efficieny
7037 * when increasing it to the n.b. cut-off.
7038 * It can even be slightly faster, because we need
7039 * less checks for the communication setup.
7041 comm->cutoff_mbody = comm->cutoff;
7044 /* Check if we did not end up below our original limit */
7045 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7047 if (comm->cutoff_mbody > comm->cellsize_limit)
7049 comm->cellsize_limit = comm->cutoff_mbody;
7052 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7057 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7058 "cellsize limit %f\n",
7059 comm->bBondComm, comm->cellsize_limit);
7064 check_dd_restrictions(cr, dd, ir, fplog);
7067 comm->partition_step = INT_MIN;
7070 clear_dd_cycle_counts(dd);
7075 static void set_dlb_limits(gmx_domdec_t *dd)
7080 for (d = 0; d < dd->ndim; d++)
7082 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7083 dd->comm->cellsize_min[dd->dim[d]] =
7084 dd->comm->cellsize_min_dlb[dd->dim[d]];
7089 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7092 gmx_domdec_comm_t *comm;
7102 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);
7105 cellsize_min = comm->cellsize_min[dd->dim[0]];
7106 for (d = 1; d < dd->ndim; d++)
7108 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7111 if (cellsize_min < comm->cellsize_limit*1.05)
7113 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");
7115 /* Change DLB from "auto" to "no". */
7116 comm->eDLB = edlbNO;
7121 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7122 comm->bDynLoadBal = TRUE;
7123 dd->bGridJump = TRUE;
7127 /* We can set the required cell size info here,
7128 * so we do not need to communicate this.
7129 * The grid is completely uniform.
7131 for (d = 0; d < dd->ndim; d++)
7135 comm->load[d].sum_m = comm->load[d].sum;
7137 nc = dd->nc[dd->dim[d]];
7138 for (i = 0; i < nc; i++)
7140 comm->root[d]->cell_f[i] = i/(real)nc;
7143 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7144 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7147 comm->root[d]->cell_f[nc] = 1.0;
7152 static char *init_bLocalCG(gmx_mtop_t *mtop)
7157 ncg = ncg_mtop(mtop);
7158 snew(bLocalCG, ncg);
7159 for (cg = 0; cg < ncg; cg++)
7161 bLocalCG[cg] = FALSE;
7167 void dd_init_bondeds(FILE *fplog,
7168 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7170 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7172 gmx_domdec_comm_t *comm;
7174 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7178 if (comm->bBondComm)
7180 /* Communicate atoms beyond the cut-off for bonded interactions */
7183 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7185 comm->bLocalCG = init_bLocalCG(mtop);
7189 /* Only communicate atoms based on cut-off */
7190 comm->cglink = NULL;
7191 comm->bLocalCG = NULL;
7195 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7197 gmx_bool bDynLoadBal, real dlb_scale,
7200 gmx_domdec_comm_t *comm;
7215 fprintf(fplog, "The maximum number of communication pulses is:");
7216 for (d = 0; d < dd->ndim; d++)
7218 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7220 fprintf(fplog, "\n");
7221 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7222 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7223 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7224 for (d = 0; d < DIM; d++)
7228 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7235 comm->cellsize_min_dlb[d]/
7236 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7238 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7241 fprintf(fplog, "\n");
7245 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7246 fprintf(fplog, "The initial number of communication pulses is:");
7247 for (d = 0; d < dd->ndim; d++)
7249 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7251 fprintf(fplog, "\n");
7252 fprintf(fplog, "The initial domain decomposition cell size is:");
7253 for (d = 0; d < DIM; d++)
7257 fprintf(fplog, " %c %.2f nm",
7258 dim2char(d), dd->comm->cellsize_min[d]);
7261 fprintf(fplog, "\n\n");
7264 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7266 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7267 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7268 "non-bonded interactions", "", comm->cutoff);
7272 limit = dd->comm->cellsize_limit;
7276 if (dynamic_dd_box(ddbox, ir))
7278 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7280 limit = dd->comm->cellsize_min[XX];
7281 for (d = 1; d < DIM; d++)
7283 limit = std::min(limit, dd->comm->cellsize_min[d]);
7287 if (comm->bInterCGBondeds)
7289 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7290 "two-body bonded interactions", "(-rdd)",
7291 std::max(comm->cutoff, comm->cutoff_mbody));
7292 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7293 "multi-body bonded interactions", "(-rdd)",
7294 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7298 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7299 "virtual site constructions", "(-rcon)", limit);
7301 if (dd->constraint_comm)
7303 sprintf(buf, "atoms separated by up to %d constraints",
7305 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7306 buf, "(-rcon)", limit);
7308 fprintf(fplog, "\n");
7314 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7316 const t_inputrec *ir,
7317 const gmx_ddbox_t *ddbox)
7319 gmx_domdec_comm_t *comm;
7320 int d, dim, npulse, npulse_d_max, npulse_d;
7325 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7327 /* Determine the maximum number of comm. pulses in one dimension */
7329 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7331 /* Determine the maximum required number of grid pulses */
7332 if (comm->cellsize_limit >= comm->cutoff)
7334 /* Only a single pulse is required */
7337 else if (!bNoCutOff && comm->cellsize_limit > 0)
7339 /* We round down slightly here to avoid overhead due to the latency
7340 * of extra communication calls when the cut-off
7341 * would be only slightly longer than the cell size.
7342 * Later cellsize_limit is redetermined,
7343 * so we can not miss interactions due to this rounding.
7345 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7349 /* There is no cell size limit */
7350 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7353 if (!bNoCutOff && npulse > 1)
7355 /* See if we can do with less pulses, based on dlb_scale */
7357 for (d = 0; d < dd->ndim; d++)
7360 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7361 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7362 npulse_d_max = std::max(npulse_d_max, npulse_d);
7364 npulse = std::min(npulse, npulse_d_max);
7367 /* This env var can override npulse */
7368 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7375 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7376 for (d = 0; d < dd->ndim; d++)
7378 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7379 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7380 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7381 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7382 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7384 comm->bVacDLBNoLimit = FALSE;
7388 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7389 if (!comm->bVacDLBNoLimit)
7391 comm->cellsize_limit = std::max(comm->cellsize_limit,
7392 comm->cutoff/comm->maxpulse);
7394 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7395 /* Set the minimum cell size for each DD dimension */
7396 for (d = 0; d < dd->ndim; d++)
7398 if (comm->bVacDLBNoLimit ||
7399 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7401 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7405 comm->cellsize_min_dlb[dd->dim[d]] =
7406 comm->cutoff/comm->cd[d].np_dlb;
7409 if (comm->cutoff_mbody <= 0)
7411 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7413 if (comm->bDynLoadBal)
7419 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7421 /* If each molecule is a single charge group
7422 * or we use domain decomposition for each periodic dimension,
7423 * we do not need to take pbc into account for the bonded interactions.
7425 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7428 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7431 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7432 t_inputrec *ir, gmx_ddbox_t *ddbox)
7434 gmx_domdec_comm_t *comm;
7440 /* Initialize the thread data.
7441 * This can not be done in init_domain_decomposition,
7442 * as the numbers of threads is determined later.
7444 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7447 snew(comm->dth, comm->nth);
7450 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7452 init_ddpme(dd, &comm->ddpme[0], 0);
7453 if (comm->npmedecompdim >= 2)
7455 init_ddpme(dd, &comm->ddpme[1], 1);
7460 comm->npmenodes = 0;
7461 if (dd->pme_nodeid >= 0)
7463 gmx_fatal_collective(FARGS, NULL, dd,
7464 "Can not have separate PME ranks without PME electrostatics");
7470 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7472 if (comm->eDLB != edlbNO)
7474 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7477 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7478 if (comm->eDLB == edlbAUTO)
7482 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7484 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7487 if (ir->ePBC == epbcNONE)
7489 vol_frac = 1 - 1/(double)dd->nnodes;
7494 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7498 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7500 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7502 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7505 static gmx_bool test_dd_cutoff(t_commrec *cr,
7506 t_state *state, t_inputrec *ir,
7517 set_ddbox(dd, FALSE, cr, ir, state->box,
7518 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7522 for (d = 0; d < dd->ndim; d++)
7526 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7527 if (dynamic_dd_box(&ddbox, ir))
7529 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7532 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7534 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7535 dd->comm->cd[d].np_dlb > 0)
7537 if (np > dd->comm->cd[d].np_dlb)
7542 /* If a current local cell size is smaller than the requested
7543 * cut-off, we could still fix it, but this gets very complicated.
7544 * Without fixing here, we might actually need more checks.
7546 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7553 if (dd->comm->eDLB != edlbNO)
7555 /* If DLB is not active yet, we don't need to check the grid jumps.
7556 * Actually we shouldn't, because then the grid jump data is not set.
7558 if (dd->comm->bDynLoadBal &&
7559 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7564 gmx_sumi(1, &LocallyLimited, cr);
7566 if (LocallyLimited > 0)
7575 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7578 gmx_bool bCutoffAllowed;
7580 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7584 cr->dd->comm->cutoff = cutoff_req;
7587 return bCutoffAllowed;
7590 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7592 gmx_domdec_comm_t *comm;
7594 comm = cr->dd->comm;
7596 /* Turn on the DLB limiting (might have been on already) */
7597 comm->bPMELoadBalDLBLimits = TRUE;
7599 /* Change the cut-off limit */
7600 comm->PMELoadBal_max_cutoff = comm->cutoff;
7603 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7605 return dd->comm->bDLB_locked;
7608 void dd_dlb_set_lock(gmx_domdec_t *dd, gmx_bool bValue)
7610 /* We can only lock the DLB when it is set to auto, otherwise don't lock */
7611 if (dd->comm->eDLB == edlbAUTO)
7613 dd->comm->bDLB_locked = bValue;
7617 static void merge_cg_buffers(int ncell,
7618 gmx_domdec_comm_dim_t *cd, int pulse,
7620 int *index_gl, int *recv_i,
7621 rvec *cg_cm, rvec *recv_vr,
7623 cginfo_mb_t *cginfo_mb, int *cginfo)
7625 gmx_domdec_ind_t *ind, *ind_p;
7626 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7627 int shift, shift_at;
7629 ind = &cd->ind[pulse];
7631 /* First correct the already stored data */
7632 shift = ind->nrecv[ncell];
7633 for (cell = ncell-1; cell >= 0; cell--)
7635 shift -= ind->nrecv[cell];
7638 /* Move the cg's present from previous grid pulses */
7639 cg0 = ncg_cell[ncell+cell];
7640 cg1 = ncg_cell[ncell+cell+1];
7641 cgindex[cg1+shift] = cgindex[cg1];
7642 for (cg = cg1-1; cg >= cg0; cg--)
7644 index_gl[cg+shift] = index_gl[cg];
7645 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7646 cgindex[cg+shift] = cgindex[cg];
7647 cginfo[cg+shift] = cginfo[cg];
7649 /* Correct the already stored send indices for the shift */
7650 for (p = 1; p <= pulse; p++)
7652 ind_p = &cd->ind[p];
7654 for (c = 0; c < cell; c++)
7656 cg0 += ind_p->nsend[c];
7658 cg1 = cg0 + ind_p->nsend[cell];
7659 for (cg = cg0; cg < cg1; cg++)
7661 ind_p->index[cg] += shift;
7667 /* Merge in the communicated buffers */
7671 for (cell = 0; cell < ncell; cell++)
7673 cg1 = ncg_cell[ncell+cell+1] + shift;
7676 /* Correct the old cg indices */
7677 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7679 cgindex[cg+1] += shift_at;
7682 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7684 /* Copy this charge group from the buffer */
7685 index_gl[cg1] = recv_i[cg0];
7686 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7687 /* Add it to the cgindex */
7688 cg_gl = index_gl[cg1];
7689 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7690 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7691 cgindex[cg1+1] = cgindex[cg1] + nat;
7696 shift += ind->nrecv[cell];
7697 ncg_cell[ncell+cell+1] = cg1;
7701 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7702 int nzone, int cg0, const int *cgindex)
7706 /* Store the atom block boundaries for easy copying of communication buffers
7709 for (zone = 0; zone < nzone; zone++)
7711 for (p = 0; p < cd->np; p++)
7713 cd->ind[p].cell2at0[zone] = cgindex[cg];
7714 cg += cd->ind[p].nrecv[zone];
7715 cd->ind[p].cell2at1[zone] = cgindex[cg];
7720 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7726 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7728 if (!bLocalCG[link->a[i]])
7737 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7739 real c[DIM][4]; /* the corners for the non-bonded communication */
7740 real cr0; /* corner for rounding */
7741 real cr1[4]; /* corners for rounding */
7742 real bc[DIM]; /* corners for bounded communication */
7743 real bcr1; /* corner for rounding for bonded communication */
7746 /* Determine the corners of the domain(s) we are communicating with */
7748 set_dd_corners(const gmx_domdec_t *dd,
7749 int dim0, int dim1, int dim2,
7753 const gmx_domdec_comm_t *comm;
7754 const gmx_domdec_zones_t *zones;
7759 zones = &comm->zones;
7761 /* Keep the compiler happy */
7765 /* The first dimension is equal for all cells */
7766 c->c[0][0] = comm->cell_x0[dim0];
7769 c->bc[0] = c->c[0][0];
7774 /* This cell row is only seen from the first row */
7775 c->c[1][0] = comm->cell_x0[dim1];
7776 /* All rows can see this row */
7777 c->c[1][1] = comm->cell_x0[dim1];
7780 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7783 /* For the multi-body distance we need the maximum */
7784 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7787 /* Set the upper-right corner for rounding */
7788 c->cr0 = comm->cell_x1[dim0];
7793 for (j = 0; j < 4; j++)
7795 c->c[2][j] = comm->cell_x0[dim2];
7799 /* Use the maximum of the i-cells that see a j-cell */
7800 for (i = 0; i < zones->nizone; i++)
7802 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7807 std::max(c->c[2][j-4],
7808 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7814 /* For the multi-body distance we need the maximum */
7815 c->bc[2] = comm->cell_x0[dim2];
7816 for (i = 0; i < 2; i++)
7818 for (j = 0; j < 2; j++)
7820 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7826 /* Set the upper-right corner for rounding */
7827 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7828 * Only cell (0,0,0) can see cell 7 (1,1,1)
7830 c->cr1[0] = comm->cell_x1[dim1];
7831 c->cr1[3] = comm->cell_x1[dim1];
7834 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7837 /* For the multi-body distance we need the maximum */
7838 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7845 /* Determine which cg's we need to send in this pulse from this zone */
7847 get_zone_pulse_cgs(gmx_domdec_t *dd,
7848 int zonei, int zone,
7850 const int *index_gl,
7852 int dim, int dim_ind,
7853 int dim0, int dim1, int dim2,
7854 real r_comm2, real r_bcomm2,
7858 real skew_fac2_d, real skew_fac_01,
7859 rvec *v_d, rvec *v_0, rvec *v_1,
7860 const dd_corners_t *c,
7862 gmx_bool bDistBonded,
7868 gmx_domdec_ind_t *ind,
7869 int **ibuf, int *ibuf_nalloc,
7875 gmx_domdec_comm_t *comm;
7877 gmx_bool bDistMB_pulse;
7879 real r2, rb2, r, tric_sh;
7882 int nsend_z, nsend, nat;
7886 bScrew = (dd->bScrewPBC && dim == XX);
7888 bDistMB_pulse = (bDistMB && bDistBonded);
7894 for (cg = cg0; cg < cg1; cg++)
7898 if (tric_dist[dim_ind] == 0)
7900 /* Rectangular direction, easy */
7901 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7908 r = cg_cm[cg][dim] - c->bc[dim_ind];
7914 /* Rounding gives at most a 16% reduction
7915 * in communicated atoms
7917 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7919 r = cg_cm[cg][dim0] - c->cr0;
7920 /* This is the first dimension, so always r >= 0 */
7927 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7929 r = cg_cm[cg][dim1] - c->cr1[zone];
7936 r = cg_cm[cg][dim1] - c->bcr1;
7946 /* Triclinic direction, more complicated */
7949 /* Rounding, conservative as the skew_fac multiplication
7950 * will slightly underestimate the distance.
7952 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7954 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7955 for (i = dim0+1; i < DIM; i++)
7957 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7959 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7962 rb[dim0] = rn[dim0];
7965 /* Take care that the cell planes along dim0 might not
7966 * be orthogonal to those along dim1 and dim2.
7968 for (i = 1; i <= dim_ind; i++)
7971 if (normal[dim0][dimd] > 0)
7973 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7976 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7981 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7983 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7985 for (i = dim1+1; i < DIM; i++)
7987 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7989 rn[dim1] += tric_sh;
7992 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7993 /* Take care of coupling of the distances
7994 * to the planes along dim0 and dim1 through dim2.
7996 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7997 /* Take care that the cell planes along dim1
7998 * might not be orthogonal to that along dim2.
8000 if (normal[dim1][dim2] > 0)
8002 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8008 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8011 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8012 /* Take care of coupling of the distances
8013 * to the planes along dim0 and dim1 through dim2.
8015 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8016 /* Take care that the cell planes along dim1
8017 * might not be orthogonal to that along dim2.
8019 if (normal[dim1][dim2] > 0)
8021 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8026 /* The distance along the communication direction */
8027 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8029 for (i = dim+1; i < DIM; i++)
8031 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8036 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8037 /* Take care of coupling of the distances
8038 * to the planes along dim0 and dim1 through dim2.
8040 if (dim_ind == 1 && zonei == 1)
8042 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8048 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8051 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8052 /* Take care of coupling of the distances
8053 * to the planes along dim0 and dim1 through dim2.
8055 if (dim_ind == 1 && zonei == 1)
8057 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8065 ((bDistMB && rb2 < r_bcomm2) ||
8066 (bDist2B && r2 < r_bcomm2)) &&
8068 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8069 missing_link(comm->cglink, index_gl[cg],
8072 /* Make an index to the local charge groups */
8073 if (nsend+1 > ind->nalloc)
8075 ind->nalloc = over_alloc_large(nsend+1);
8076 srenew(ind->index, ind->nalloc);
8078 if (nsend+1 > *ibuf_nalloc)
8080 *ibuf_nalloc = over_alloc_large(nsend+1);
8081 srenew(*ibuf, *ibuf_nalloc);
8083 ind->index[nsend] = cg;
8084 (*ibuf)[nsend] = index_gl[cg];
8086 vec_rvec_check_alloc(vbuf, nsend+1);
8088 if (dd->ci[dim] == 0)
8090 /* Correct cg_cm for pbc */
8091 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8094 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8095 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8100 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8103 nat += cgindex[cg+1] - cgindex[cg];
8109 *nsend_z_ptr = nsend_z;
8112 static void setup_dd_communication(gmx_domdec_t *dd,
8113 matrix box, gmx_ddbox_t *ddbox,
8114 t_forcerec *fr, t_state *state, rvec **f)
8116 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8117 int nzone, nzone_send, zone, zonei, cg0, cg1;
8118 int c, i, cg, cg_gl, nrcg;
8119 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8120 gmx_domdec_comm_t *comm;
8121 gmx_domdec_zones_t *zones;
8122 gmx_domdec_comm_dim_t *cd;
8123 gmx_domdec_ind_t *ind;
8124 cginfo_mb_t *cginfo_mb;
8125 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8126 real r_comm2, r_bcomm2;
8127 dd_corners_t corners;
8129 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8130 real skew_fac2_d, skew_fac_01;
8137 fprintf(debug, "Setting up DD communication\n");
8142 switch (fr->cutoff_scheme)
8151 gmx_incons("unimplemented");
8155 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8157 /* Check if we need to use triclinic distances */
8158 tric_dist[dim_ind] = 0;
8159 for (i = 0; i <= dim_ind; i++)
8161 if (ddbox->tric_dir[dd->dim[i]])
8163 tric_dist[dim_ind] = 1;
8168 bBondComm = comm->bBondComm;
8170 /* Do we need to determine extra distances for multi-body bondeds? */
8171 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8173 /* Do we need to determine extra distances for only two-body bondeds? */
8174 bDist2B = (bBondComm && !bDistMB);
8176 r_comm2 = sqr(comm->cutoff);
8177 r_bcomm2 = sqr(comm->cutoff_mbody);
8181 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8184 zones = &comm->zones;
8187 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8188 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8190 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8192 /* Triclinic stuff */
8193 normal = ddbox->normal;
8197 v_0 = ddbox->v[dim0];
8198 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8200 /* Determine the coupling coefficient for the distances
8201 * to the cell planes along dim0 and dim1 through dim2.
8202 * This is required for correct rounding.
8205 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8208 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8214 v_1 = ddbox->v[dim1];
8217 zone_cg_range = zones->cg_range;
8218 index_gl = dd->index_gl;
8219 cgindex = dd->cgindex;
8220 cginfo_mb = fr->cginfo_mb;
8222 zone_cg_range[0] = 0;
8223 zone_cg_range[1] = dd->ncg_home;
8224 comm->zone_ncg1[0] = dd->ncg_home;
8225 pos_cg = dd->ncg_home;
8227 nat_tot = dd->nat_home;
8229 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8231 dim = dd->dim[dim_ind];
8232 cd = &comm->cd[dim_ind];
8234 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8236 /* No pbc in this dimension, the first node should not comm. */
8244 v_d = ddbox->v[dim];
8245 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8247 cd->bInPlace = TRUE;
8248 for (p = 0; p < cd->np; p++)
8250 /* Only atoms communicated in the first pulse are used
8251 * for multi-body bonded interactions or for bBondComm.
8253 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8258 for (zone = 0; zone < nzone_send; zone++)
8260 if (tric_dist[dim_ind] && dim_ind > 0)
8262 /* Determine slightly more optimized skew_fac's
8264 * This reduces the number of communicated atoms
8265 * by about 10% for 3D DD of rhombic dodecahedra.
8267 for (dimd = 0; dimd < dim; dimd++)
8269 sf2_round[dimd] = 1;
8270 if (ddbox->tric_dir[dimd])
8272 for (i = dd->dim[dimd]+1; i < DIM; i++)
8274 /* If we are shifted in dimension i
8275 * and the cell plane is tilted forward
8276 * in dimension i, skip this coupling.
8278 if (!(zones->shift[nzone+zone][i] &&
8279 ddbox->v[dimd][i][dimd] >= 0))
8282 sqr(ddbox->v[dimd][i][dimd]);
8285 sf2_round[dimd] = 1/sf2_round[dimd];
8290 zonei = zone_perm[dim_ind][zone];
8293 /* Here we permutate the zones to obtain a convenient order
8294 * for neighbor searching
8296 cg0 = zone_cg_range[zonei];
8297 cg1 = zone_cg_range[zonei+1];
8301 /* Look only at the cg's received in the previous grid pulse
8303 cg1 = zone_cg_range[nzone+zone+1];
8304 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8307 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8308 for (th = 0; th < comm->nth; th++)
8310 gmx_domdec_ind_t *ind_p;
8311 int **ibuf_p, *ibuf_nalloc_p;
8313 int *nsend_p, *nat_p;
8319 /* Thread 0 writes in the comm buffers */
8321 ibuf_p = &comm->buf_int;
8322 ibuf_nalloc_p = &comm->nalloc_int;
8323 vbuf_p = &comm->vbuf;
8326 nsend_zone_p = &ind->nsend[zone];
8330 /* Other threads write into temp buffers */
8331 ind_p = &comm->dth[th].ind;
8332 ibuf_p = &comm->dth[th].ibuf;
8333 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8334 vbuf_p = &comm->dth[th].vbuf;
8335 nsend_p = &comm->dth[th].nsend;
8336 nat_p = &comm->dth[th].nat;
8337 nsend_zone_p = &comm->dth[th].nsend_zone;
8339 comm->dth[th].nsend = 0;
8340 comm->dth[th].nat = 0;
8341 comm->dth[th].nsend_zone = 0;
8351 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8352 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8355 /* Get the cg's for this pulse in this zone */
8356 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8358 dim, dim_ind, dim0, dim1, dim2,
8361 normal, skew_fac2_d, skew_fac_01,
8362 v_d, v_0, v_1, &corners, sf2_round,
8363 bDistBonded, bBondComm,
8367 ibuf_p, ibuf_nalloc_p,
8373 /* Append data of threads>=1 to the communication buffers */
8374 for (th = 1; th < comm->nth; th++)
8376 dd_comm_setup_work_t *dth;
8379 dth = &comm->dth[th];
8381 ns1 = nsend + dth->nsend_zone;
8382 if (ns1 > ind->nalloc)
8384 ind->nalloc = over_alloc_dd(ns1);
8385 srenew(ind->index, ind->nalloc);
8387 if (ns1 > comm->nalloc_int)
8389 comm->nalloc_int = over_alloc_dd(ns1);
8390 srenew(comm->buf_int, comm->nalloc_int);
8392 if (ns1 > comm->vbuf.nalloc)
8394 comm->vbuf.nalloc = over_alloc_dd(ns1);
8395 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8398 for (i = 0; i < dth->nsend_zone; i++)
8400 ind->index[nsend] = dth->ind.index[i];
8401 comm->buf_int[nsend] = dth->ibuf[i];
8402 copy_rvec(dth->vbuf.v[i],
8403 comm->vbuf.v[nsend]);
8407 ind->nsend[zone] += dth->nsend_zone;
8410 /* Clear the counts in case we do not have pbc */
8411 for (zone = nzone_send; zone < nzone; zone++)
8413 ind->nsend[zone] = 0;
8415 ind->nsend[nzone] = nsend;
8416 ind->nsend[nzone+1] = nat;
8417 /* Communicate the number of cg's and atoms to receive */
8418 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8419 ind->nsend, nzone+2,
8420 ind->nrecv, nzone+2);
8422 /* The rvec buffer is also required for atom buffers of size nsend
8423 * in dd_move_x and dd_move_f.
8425 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8429 /* We can receive in place if only the last zone is not empty */
8430 for (zone = 0; zone < nzone-1; zone++)
8432 if (ind->nrecv[zone] > 0)
8434 cd->bInPlace = FALSE;
8439 /* The int buffer is only required here for the cg indices */
8440 if (ind->nrecv[nzone] > comm->nalloc_int2)
8442 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8443 srenew(comm->buf_int2, comm->nalloc_int2);
8445 /* The rvec buffer is also required for atom buffers
8446 * of size nrecv in dd_move_x and dd_move_f.
8448 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8449 vec_rvec_check_alloc(&comm->vbuf2, i);
8453 /* Make space for the global cg indices */
8454 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8455 || dd->cg_nalloc == 0)
8457 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8458 srenew(index_gl, dd->cg_nalloc);
8459 srenew(cgindex, dd->cg_nalloc+1);
8461 /* Communicate the global cg indices */
8464 recv_i = index_gl + pos_cg;
8468 recv_i = comm->buf_int2;
8470 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8471 comm->buf_int, nsend,
8472 recv_i, ind->nrecv[nzone]);
8474 /* Make space for cg_cm */
8475 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8476 if (fr->cutoff_scheme == ecutsGROUP)
8484 /* Communicate cg_cm */
8487 recv_vr = cg_cm + pos_cg;
8491 recv_vr = comm->vbuf2.v;
8493 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8494 comm->vbuf.v, nsend,
8495 recv_vr, ind->nrecv[nzone]);
8497 /* Make the charge group index */
8500 zone = (p == 0 ? 0 : nzone - 1);
8501 while (zone < nzone)
8503 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8505 cg_gl = index_gl[pos_cg];
8506 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8507 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8508 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8511 /* Update the charge group presence,
8512 * so we can use it in the next pass of the loop.
8514 comm->bLocalCG[cg_gl] = TRUE;
8520 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8523 zone_cg_range[nzone+zone] = pos_cg;
8528 /* This part of the code is never executed with bBondComm. */
8529 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8530 index_gl, recv_i, cg_cm, recv_vr,
8531 cgindex, fr->cginfo_mb, fr->cginfo);
8532 pos_cg += ind->nrecv[nzone];
8534 nat_tot += ind->nrecv[nzone+1];
8538 /* Store the atom block for easy copying of communication buffers */
8539 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8543 dd->index_gl = index_gl;
8544 dd->cgindex = cgindex;
8546 dd->ncg_tot = zone_cg_range[zones->n];
8547 dd->nat_tot = nat_tot;
8548 comm->nat[ddnatHOME] = dd->nat_home;
8549 for (i = ddnatZONE; i < ddnatNR; i++)
8551 comm->nat[i] = dd->nat_tot;
8556 /* We don't need to update cginfo, since that was alrady done above.
8557 * So we pass NULL for the forcerec.
8559 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8560 NULL, comm->bLocalCG);
8565 fprintf(debug, "Finished setting up DD communication, zones:");
8566 for (c = 0; c < zones->n; c++)
8568 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8570 fprintf(debug, "\n");
8574 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8578 for (c = 0; c < zones->nizone; c++)
8580 zones->izone[c].cg1 = zones->cg_range[c+1];
8581 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8582 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8586 static void set_zones_size(gmx_domdec_t *dd,
8587 matrix box, const gmx_ddbox_t *ddbox,
8588 int zone_start, int zone_end)
8590 gmx_domdec_comm_t *comm;
8591 gmx_domdec_zones_t *zones;
8600 zones = &comm->zones;
8602 /* Do we need to determine extra distances for multi-body bondeds? */
8603 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8605 for (z = zone_start; z < zone_end; z++)
8607 /* Copy cell limits to zone limits.
8608 * Valid for non-DD dims and non-shifted dims.
8610 copy_rvec(comm->cell_x0, zones->size[z].x0);
8611 copy_rvec(comm->cell_x1, zones->size[z].x1);
8614 for (d = 0; d < dd->ndim; d++)
8618 for (z = 0; z < zones->n; z++)
8620 /* With a staggered grid we have different sizes
8621 * for non-shifted dimensions.
8623 if (dd->bGridJump && zones->shift[z][dim] == 0)
8627 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8628 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8632 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8633 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8639 rcmbs = comm->cutoff_mbody;
8640 if (ddbox->tric_dir[dim])
8642 rcs /= ddbox->skew_fac[dim];
8643 rcmbs /= ddbox->skew_fac[dim];
8646 /* Set the lower limit for the shifted zone dimensions */
8647 for (z = zone_start; z < zone_end; z++)
8649 if (zones->shift[z][dim] > 0)
8652 if (!dd->bGridJump || d == 0)
8654 zones->size[z].x0[dim] = comm->cell_x1[dim];
8655 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8659 /* Here we take the lower limit of the zone from
8660 * the lowest domain of the zone below.
8664 zones->size[z].x0[dim] =
8665 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8671 zones->size[z].x0[dim] =
8672 zones->size[zone_perm[2][z-4]].x0[dim];
8676 zones->size[z].x0[dim] =
8677 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8680 /* A temporary limit, is updated below */
8681 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8685 for (zi = 0; zi < zones->nizone; zi++)
8687 if (zones->shift[zi][dim] == 0)
8689 /* This takes the whole zone into account.
8690 * With multiple pulses this will lead
8691 * to a larger zone then strictly necessary.
8693 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8694 zones->size[zi].x1[dim]+rcmbs);
8702 /* Loop over the i-zones to set the upper limit of each
8705 for (zi = 0; zi < zones->nizone; zi++)
8707 if (zones->shift[zi][dim] == 0)
8709 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8711 if (zones->shift[z][dim] > 0)
8713 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8714 zones->size[zi].x1[dim]+rcs);
8721 for (z = zone_start; z < zone_end; z++)
8723 /* Initialization only required to keep the compiler happy */
8724 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8727 /* To determine the bounding box for a zone we need to find
8728 * the extreme corners of 4, 2 or 1 corners.
8730 nc = 1 << (ddbox->npbcdim - 1);
8732 for (c = 0; c < nc; c++)
8734 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8738 corner[YY] = zones->size[z].x0[YY];
8742 corner[YY] = zones->size[z].x1[YY];
8746 corner[ZZ] = zones->size[z].x0[ZZ];
8750 corner[ZZ] = zones->size[z].x1[ZZ];
8752 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8754 /* With 1D domain decomposition the cg's are not in
8755 * the triclinic box, but triclinic x-y and rectangular y-z.
8756 * Shift y back, so it will later end up at 0.
8758 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8760 /* Apply the triclinic couplings */
8761 assert(ddbox->npbcdim <= DIM);
8762 for (i = YY; i < ddbox->npbcdim; i++)
8764 for (j = XX; j < i; j++)
8766 corner[j] += corner[i]*box[i][j]/box[i][i];
8771 copy_rvec(corner, corner_min);
8772 copy_rvec(corner, corner_max);
8776 for (i = 0; i < DIM; i++)
8778 corner_min[i] = std::min(corner_min[i], corner[i]);
8779 corner_max[i] = std::max(corner_max[i], corner[i]);
8783 /* Copy the extreme cornes without offset along x */
8784 for (i = 0; i < DIM; i++)
8786 zones->size[z].bb_x0[i] = corner_min[i];
8787 zones->size[z].bb_x1[i] = corner_max[i];
8789 /* Add the offset along x */
8790 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8791 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8794 if (zone_start == 0)
8797 for (dim = 0; dim < DIM; dim++)
8799 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8801 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8806 for (z = zone_start; z < zone_end; z++)
8808 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8810 zones->size[z].x0[XX], zones->size[z].x1[XX],
8811 zones->size[z].x0[YY], zones->size[z].x1[YY],
8812 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8813 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8815 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8816 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8817 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8822 static int comp_cgsort(const void *a, const void *b)
8826 gmx_cgsort_t *cga, *cgb;
8827 cga = (gmx_cgsort_t *)a;
8828 cgb = (gmx_cgsort_t *)b;
8830 comp = cga->nsc - cgb->nsc;
8833 comp = cga->ind_gl - cgb->ind_gl;
8839 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8844 /* Order the data */
8845 for (i = 0; i < n; i++)
8847 buf[i] = a[sort[i].ind];
8850 /* Copy back to the original array */
8851 for (i = 0; i < n; i++)
8857 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8862 /* Order the data */
8863 for (i = 0; i < n; i++)
8865 copy_rvec(v[sort[i].ind], buf[i]);
8868 /* Copy back to the original array */
8869 for (i = 0; i < n; i++)
8871 copy_rvec(buf[i], v[i]);
8875 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8878 int a, atot, cg, cg0, cg1, i;
8880 if (cgindex == NULL)
8882 /* Avoid the useless loop of the atoms within a cg */
8883 order_vec_cg(ncg, sort, v, buf);
8888 /* Order the data */
8890 for (cg = 0; cg < ncg; cg++)
8892 cg0 = cgindex[sort[cg].ind];
8893 cg1 = cgindex[sort[cg].ind+1];
8894 for (i = cg0; i < cg1; i++)
8896 copy_rvec(v[i], buf[a]);
8902 /* Copy back to the original array */
8903 for (a = 0; a < atot; a++)
8905 copy_rvec(buf[a], v[a]);
8909 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8910 int nsort_new, gmx_cgsort_t *sort_new,
8911 gmx_cgsort_t *sort1)
8915 /* The new indices are not very ordered, so we qsort them */
8916 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8918 /* sort2 is already ordered, so now we can merge the two arrays */
8922 while (i2 < nsort2 || i_new < nsort_new)
8926 sort1[i1++] = sort_new[i_new++];
8928 else if (i_new == nsort_new)
8930 sort1[i1++] = sort2[i2++];
8932 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8933 (sort2[i2].nsc == sort_new[i_new].nsc &&
8934 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8936 sort1[i1++] = sort2[i2++];
8940 sort1[i1++] = sort_new[i_new++];
8945 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8947 gmx_domdec_sort_t *sort;
8948 gmx_cgsort_t *cgsort, *sort_i;
8949 int ncg_new, nsort2, nsort_new, i, *a, moved;
8951 sort = dd->comm->sort;
8953 a = fr->ns.grid->cell_index;
8955 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8957 if (ncg_home_old >= 0)
8959 /* The charge groups that remained in the same ns grid cell
8960 * are completely ordered. So we can sort efficiently by sorting
8961 * the charge groups that did move into the stationary list.
8966 for (i = 0; i < dd->ncg_home; i++)
8968 /* Check if this cg did not move to another node */
8971 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8973 /* This cg is new on this node or moved ns grid cell */
8974 if (nsort_new >= sort->sort_new_nalloc)
8976 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8977 srenew(sort->sort_new, sort->sort_new_nalloc);
8979 sort_i = &(sort->sort_new[nsort_new++]);
8983 /* This cg did not move */
8984 sort_i = &(sort->sort2[nsort2++]);
8986 /* Sort on the ns grid cell indices
8987 * and the global topology index.
8988 * index_gl is irrelevant with cell ns,
8989 * but we set it here anyhow to avoid a conditional.
8992 sort_i->ind_gl = dd->index_gl[i];
8999 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9002 /* Sort efficiently */
9003 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9008 cgsort = sort->sort;
9010 for (i = 0; i < dd->ncg_home; i++)
9012 /* Sort on the ns grid cell indices
9013 * and the global topology index
9015 cgsort[i].nsc = a[i];
9016 cgsort[i].ind_gl = dd->index_gl[i];
9018 if (cgsort[i].nsc < moved)
9025 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9027 /* Determine the order of the charge groups using qsort */
9028 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9034 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9037 int ncg_new, i, *a, na;
9039 sort = dd->comm->sort->sort;
9041 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9044 for (i = 0; i < na; i++)
9048 sort[ncg_new].ind = a[i];
9056 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9059 gmx_domdec_sort_t *sort;
9060 gmx_cgsort_t *cgsort;
9062 int ncg_new, i, *ibuf, cgsize;
9065 sort = dd->comm->sort;
9067 if (dd->ncg_home > sort->sort_nalloc)
9069 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9070 srenew(sort->sort, sort->sort_nalloc);
9071 srenew(sort->sort2, sort->sort_nalloc);
9073 cgsort = sort->sort;
9075 switch (fr->cutoff_scheme)
9078 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9081 ncg_new = dd_sort_order_nbnxn(dd, fr);
9084 gmx_incons("unimplemented");
9088 /* We alloc with the old size, since cgindex is still old */
9089 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9090 vbuf = dd->comm->vbuf.v;
9094 cgindex = dd->cgindex;
9101 /* Remove the charge groups which are no longer at home here */
9102 dd->ncg_home = ncg_new;
9105 fprintf(debug, "Set the new home charge group count to %d\n",
9109 /* Reorder the state */
9110 for (i = 0; i < estNR; i++)
9112 if (EST_DISTR(i) && (state->flags & (1<<i)))
9117 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9120 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9123 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9126 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9130 case estDISRE_INITF:
9131 case estDISRE_RM3TAV:
9132 case estORIRE_INITF:
9134 /* No ordering required */
9137 gmx_incons("Unknown state entry encountered in dd_sort_state");
9142 if (fr->cutoff_scheme == ecutsGROUP)
9145 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9148 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9150 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9151 srenew(sort->ibuf, sort->ibuf_nalloc);
9154 /* Reorder the global cg index */
9155 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9156 /* Reorder the cginfo */
9157 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9158 /* Rebuild the local cg index */
9162 for (i = 0; i < dd->ncg_home; i++)
9164 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9165 ibuf[i+1] = ibuf[i] + cgsize;
9167 for (i = 0; i < dd->ncg_home+1; i++)
9169 dd->cgindex[i] = ibuf[i];
9174 for (i = 0; i < dd->ncg_home+1; i++)
9179 /* Set the home atom number */
9180 dd->nat_home = dd->cgindex[dd->ncg_home];
9182 if (fr->cutoff_scheme == ecutsVERLET)
9184 /* The atoms are now exactly in grid order, update the grid order */
9185 nbnxn_set_atomorder(fr->nbv->nbs);
9189 /* Copy the sorted ns cell indices back to the ns grid struct */
9190 for (i = 0; i < dd->ncg_home; i++)
9192 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9194 fr->ns.grid->nr = dd->ncg_home;
9198 static void add_dd_statistics(gmx_domdec_t *dd)
9200 gmx_domdec_comm_t *comm;
9205 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9207 comm->sum_nat[ddnat-ddnatZONE] +=
9208 comm->nat[ddnat] - comm->nat[ddnat-1];
9213 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9215 gmx_domdec_comm_t *comm;
9220 /* Reset all the statistics and counters for total run counting */
9221 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9223 comm->sum_nat[ddnat-ddnatZONE] = 0;
9227 comm->load_step = 0;
9230 clear_ivec(comm->load_lim);
9235 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9237 gmx_domdec_comm_t *comm;
9241 comm = cr->dd->comm;
9243 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9250 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");
9252 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9254 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9259 " av. #atoms communicated per step for force: %d x %.1f\n",
9263 if (cr->dd->vsite_comm)
9266 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9267 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9272 if (cr->dd->constraint_comm)
9275 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9276 1 + ir->nLincsIter, av);
9280 gmx_incons(" Unknown type for DD statistics");
9283 fprintf(fplog, "\n");
9285 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9287 print_dd_load_av(fplog, cr->dd);
9291 void dd_partition_system(FILE *fplog,
9294 gmx_bool bMasterState,
9296 t_state *state_global,
9297 gmx_mtop_t *top_global,
9299 t_state *state_local,
9302 gmx_localtop_t *top_local,
9305 gmx_shellfc_t shellfc,
9306 gmx_constr_t constr,
9308 gmx_wallcycle_t wcycle,
9312 gmx_domdec_comm_t *comm;
9313 gmx_ddbox_t ddbox = {0};
9315 gmx_int64_t step_pcoupl;
9316 rvec cell_ns_x0, cell_ns_x1;
9317 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9318 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9319 gmx_bool bRedist, bSortCG, bResortAll;
9320 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9327 bBoxChanged = (bMasterState || DEFORM(*ir));
9328 if (ir->epc != epcNO)
9330 /* With nstpcouple > 1 pressure coupling happens.
9331 * one step after calculating the pressure.
9332 * Box scaling happens at the end of the MD step,
9333 * after the DD partitioning.
9334 * We therefore have to do DLB in the first partitioning
9335 * after an MD step where P-coupling occured.
9336 * We need to determine the last step in which p-coupling occurred.
9337 * MRS -- need to validate this for vv?
9342 step_pcoupl = step - 1;
9346 step_pcoupl = ((step - 1)/n)*n + 1;
9348 if (step_pcoupl >= comm->partition_step)
9354 bNStGlobalComm = (step % nstglobalcomm == 0);
9356 if (!comm->bDynLoadBal)
9362 /* Should we do dynamic load balacing this step?
9363 * Since it requires (possibly expensive) global communication,
9364 * we might want to do DLB less frequently.
9366 if (bBoxChanged || ir->epc != epcNO)
9368 bDoDLB = bBoxChanged;
9372 bDoDLB = bNStGlobalComm;
9376 /* Check if we have recorded loads on the nodes */
9377 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9379 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal && !dd_dlb_is_locked(dd))
9381 /* Check if we should use DLB at the second partitioning
9382 * and every 100 partitionings,
9383 * so the extra communication cost is negligible.
9385 const int nddp_chk_dlb = 100;
9386 bCheckDLB = (comm->n_load_collect == 0 ||
9387 comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
9394 /* Print load every nstlog, first and last step to the log file */
9395 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9396 comm->n_load_collect == 0 ||
9398 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9400 /* Avoid extra communication due to verbose screen output
9401 * when nstglobalcomm is set.
9403 if (bDoDLB || bLogLoad || bCheckDLB ||
9404 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9406 get_load_distribution(dd, wcycle);
9411 dd_print_load(fplog, dd, step-1);
9415 dd_print_load_verbose(dd);
9418 comm->n_load_collect++;
9422 /* Since the timings are node dependent, the master decides */
9425 /* Here we check if the max PME rank load is more than 0.98
9426 * the max PP force load. If so, PP DLB will not help,
9427 * since we are (almost) limited by PME. Furthermore,
9428 * DLB will cause a significant extra x/f redistribution
9429 * cost on the PME ranks, which will then surely result
9430 * in lower total performance.
9431 * This check might be fragile, since one measurement
9432 * below 0.98 (although only done once every 100 DD part.)
9433 * could turn on DLB for the rest of the run.
9435 if (cr->npmenodes > 0 &&
9436 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9443 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9447 fprintf(debug, "step %s, imb loss %f\n",
9448 gmx_step_str(step, sbuf),
9449 dd_force_imb_perf_loss(dd));
9452 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9455 turn_on_dlb(fplog, cr, step);
9460 comm->n_load_have++;
9463 cgs_gl = &comm->cgs_gl;
9468 /* Clear the old state */
9469 clear_dd_indices(dd, 0, 0);
9472 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9473 TRUE, cgs_gl, state_global->x, &ddbox);
9475 get_cg_distribution(fplog, step, dd, cgs_gl,
9476 state_global->box, &ddbox, state_global->x);
9478 dd_distribute_state(dd, cgs_gl,
9479 state_global, state_local, f);
9481 dd_make_local_cgs(dd, &top_local->cgs);
9483 /* Ensure that we have space for the new distribution */
9484 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9486 if (fr->cutoff_scheme == ecutsGROUP)
9488 calc_cgcm(fplog, 0, dd->ncg_home,
9489 &top_local->cgs, state_local->x, fr->cg_cm);
9492 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9494 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9496 else if (state_local->ddp_count != dd->ddp_count)
9498 if (state_local->ddp_count > dd->ddp_count)
9500 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9503 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9505 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);
9508 /* Clear the old state */
9509 clear_dd_indices(dd, 0, 0);
9511 /* Build the new indices */
9512 rebuild_cgindex(dd, cgs_gl->index, state_local);
9513 make_dd_indices(dd, cgs_gl->index, 0);
9514 ncgindex_set = dd->ncg_home;
9516 if (fr->cutoff_scheme == ecutsGROUP)
9518 /* Redetermine the cg COMs */
9519 calc_cgcm(fplog, 0, dd->ncg_home,
9520 &top_local->cgs, state_local->x, fr->cg_cm);
9523 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9525 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9527 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9528 TRUE, &top_local->cgs, state_local->x, &ddbox);
9530 bRedist = comm->bDynLoadBal;
9534 /* We have the full state, only redistribute the cgs */
9536 /* Clear the non-home indices */
9537 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9540 /* Avoid global communication for dim's without pbc and -gcom */
9541 if (!bNStGlobalComm)
9543 copy_rvec(comm->box0, ddbox.box0 );
9544 copy_rvec(comm->box_size, ddbox.box_size);
9546 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9547 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9552 /* For dim's without pbc and -gcom */
9553 copy_rvec(ddbox.box0, comm->box0 );
9554 copy_rvec(ddbox.box_size, comm->box_size);
9556 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9559 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9561 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9564 /* Check if we should sort the charge groups */
9565 if (comm->nstSortCG > 0)
9567 bSortCG = (bMasterState ||
9568 (bRedist && (step % comm->nstSortCG == 0)));
9575 ncg_home_old = dd->ncg_home;
9580 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9582 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9584 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9586 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9589 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9591 &comm->cell_x0, &comm->cell_x1,
9592 dd->ncg_home, fr->cg_cm,
9593 cell_ns_x0, cell_ns_x1, &grid_density);
9597 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9600 switch (fr->cutoff_scheme)
9603 copy_ivec(fr->ns.grid->n, ncells_old);
9604 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9605 state_local->box, cell_ns_x0, cell_ns_x1,
9606 fr->rlistlong, grid_density);
9609 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9612 gmx_incons("unimplemented");
9614 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9615 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9619 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9621 /* Sort the state on charge group position.
9622 * This enables exact restarts from this step.
9623 * It also improves performance by about 15% with larger numbers
9624 * of atoms per node.
9627 /* Fill the ns grid with the home cell,
9628 * so we can sort with the indices.
9630 set_zones_ncg_home(dd);
9632 switch (fr->cutoff_scheme)
9635 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9637 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9639 comm->zones.size[0].bb_x0,
9640 comm->zones.size[0].bb_x1,
9642 comm->zones.dens_zone0,
9645 ncg_moved, bRedist ? comm->moved : NULL,
9646 fr->nbv->grp[eintLocal].kernel_type,
9647 fr->nbv->grp[eintLocal].nbat);
9649 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9652 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9653 0, dd->ncg_home, fr->cg_cm);
9655 copy_ivec(fr->ns.grid->n, ncells_new);
9658 gmx_incons("unimplemented");
9661 bResortAll = bMasterState;
9663 /* Check if we can user the old order and ns grid cell indices
9664 * of the charge groups to sort the charge groups efficiently.
9666 if (ncells_new[XX] != ncells_old[XX] ||
9667 ncells_new[YY] != ncells_old[YY] ||
9668 ncells_new[ZZ] != ncells_old[ZZ])
9675 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9676 gmx_step_str(step, sbuf), dd->ncg_home);
9678 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9679 bResortAll ? -1 : ncg_home_old);
9680 /* Rebuild all the indices */
9681 ga2la_clear(dd->ga2la);
9684 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9687 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9689 /* Setup up the communication and communicate the coordinates */
9690 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9692 /* Set the indices */
9693 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9695 /* Set the charge group boundaries for neighbor searching */
9696 set_cg_boundaries(&comm->zones);
9698 if (fr->cutoff_scheme == ecutsVERLET)
9700 set_zones_size(dd, state_local->box, &ddbox,
9701 bSortCG ? 1 : 0, comm->zones.n);
9704 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9707 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9708 -1,state_local->x,state_local->box);
9711 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9713 /* Extract a local topology from the global topology */
9714 for (i = 0; i < dd->ndim; i++)
9716 np[dd->dim[i]] = comm->cd[i].np;
9718 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9719 comm->cellsize_min, np,
9721 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9722 vsite, top_global, top_local);
9724 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9726 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9728 /* Set up the special atom communication */
9729 n = comm->nat[ddnatZONE];
9730 for (i = ddnatZONE+1; i < ddnatNR; i++)
9735 if (vsite && vsite->n_intercg_vsite)
9737 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9741 if (dd->bInterCGcons || dd->bInterCGsettles)
9743 /* Only for inter-cg constraints we need special code */
9744 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9745 constr, ir->nProjOrder,
9746 top_local->idef.il);
9750 gmx_incons("Unknown special atom type setup");
9755 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9757 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9759 /* Make space for the extra coordinates for virtual site
9760 * or constraint communication.
9762 state_local->natoms = comm->nat[ddnatNR-1];
9763 if (state_local->natoms > state_local->nalloc)
9765 dd_realloc_state(state_local, f, state_local->natoms);
9768 if (fr->bF_NoVirSum)
9770 if (vsite && vsite->n_intercg_vsite)
9772 nat_f_novirsum = comm->nat[ddnatVSITE];
9776 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9778 nat_f_novirsum = dd->nat_tot;
9782 nat_f_novirsum = dd->nat_home;
9791 /* Set the number of atoms required for the force calculation.
9792 * Forces need to be constrained when using a twin-range setup
9793 * or with energy minimization. For simple simulations we could
9794 * avoid some allocation, zeroing and copying, but this is
9795 * probably not worth the complications ande checking.
9797 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9798 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9800 /* We make the all mdatoms up to nat_tot_con.
9801 * We could save some work by only setting invmass
9802 * between nat_tot and nat_tot_con.
9804 /* This call also sets the new number of home particles to dd->nat_home */
9805 atoms2md(top_global, ir,
9806 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9808 /* Now we have the charges we can sort the FE interactions */
9809 dd_sort_local_top(dd, mdatoms, top_local);
9813 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9814 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9815 mdatoms, FALSE, vsite);
9820 /* Make the local shell stuff, currently no communication is done */
9821 make_local_shells(cr, mdatoms, shellfc);
9824 if (ir->implicit_solvent)
9826 make_local_gb(cr, fr->born, ir->gb_algorithm);
9829 setup_bonded_threading(fr, &top_local->idef);
9831 if (!(cr->duty & DUTY_PME))
9833 /* Send the charges and/or c6/sigmas to our PME only node */
9834 gmx_pme_send_parameters(cr,
9836 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9837 mdatoms->chargeA, mdatoms->chargeB,
9838 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9839 mdatoms->sigmaA, mdatoms->sigmaB,
9840 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9845 set_constraints(constr, top_local, ir, mdatoms, cr);
9848 if (ir->ePull != epullNO)
9850 /* Update the local pull groups */
9851 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9856 /* Update the local rotation groups */
9857 dd_make_local_rotation_groups(dd, ir->rot);
9860 if (ir->eSwapCoords != eswapNO)
9862 /* Update the local groups needed for ion swapping */
9863 dd_make_local_swap_groups(dd, ir->swap);
9866 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9867 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9869 add_dd_statistics(dd);
9871 /* Make sure we only count the cycles for this DD partitioning */
9872 clear_dd_cycle_counts(dd);
9874 /* Because the order of the atoms might have changed since
9875 * the last vsite construction, we need to communicate the constructing
9876 * atom coordinates again (for spreading the forces this MD step).
9878 dd_move_x_vsites(dd, state_local->box, state_local->x);
9880 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9882 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9884 dd_move_x(dd, state_local->box, state_local->x);
9885 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9886 -1, state_local->x, state_local->box);
9889 /* Store the partitioning step */
9890 comm->partition_step = step;
9892 /* Increase the DD partitioning counter */
9894 /* The state currently matches this DD partitioning count, store it */
9895 state_local->ddp_count = dd->ddp_count;
9898 /* The DD master node knows the complete cg distribution,
9899 * store the count so we can possibly skip the cg info communication.
9901 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9904 if (comm->DD_debug > 0)
9906 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9907 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9908 "after partitioning");