2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gromacs/domdec/domdec_network.h"
52 #include "gromacs/ewald/pme.h"
53 #include "gromacs/fileio/gmxfio.h"
54 #include "gromacs/fileio/pdbio.h"
55 #include "gromacs/gmxlib/gpu_utils/gpu_utils.h"
56 #include "gromacs/imd/imd.h"
57 #include "gromacs/legacyheaders/chargegroup.h"
58 #include "gromacs/legacyheaders/constr.h"
59 #include "gromacs/legacyheaders/force.h"
60 #include "gromacs/legacyheaders/genborn.h"
61 #include "gromacs/legacyheaders/gmx_ga2la.h"
62 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
63 #include "gromacs/legacyheaders/mdatoms.h"
64 #include "gromacs/legacyheaders/mdrun.h"
65 #include "gromacs/legacyheaders/names.h"
66 #include "gromacs/legacyheaders/network.h"
67 #include "gromacs/legacyheaders/nrnb.h"
68 #include "gromacs/legacyheaders/nsgrid.h"
69 #include "gromacs/legacyheaders/shellfc.h"
70 #include "gromacs/legacyheaders/typedefs.h"
71 #include "gromacs/legacyheaders/vsite.h"
72 #include "gromacs/legacyheaders/types/commrec.h"
73 #include "gromacs/legacyheaders/types/constr.h"
74 #include "gromacs/legacyheaders/types/enums.h"
75 #include "gromacs/legacyheaders/types/forcerec.h"
76 #include "gromacs/legacyheaders/types/hw_info.h"
77 #include "gromacs/legacyheaders/types/ifunc.h"
78 #include "gromacs/legacyheaders/types/inputrec.h"
79 #include "gromacs/legacyheaders/types/mdatom.h"
80 #include "gromacs/legacyheaders/types/nrnb.h"
81 #include "gromacs/legacyheaders/types/ns.h"
82 #include "gromacs/legacyheaders/types/nsgrid.h"
83 #include "gromacs/legacyheaders/types/shellfc.h"
84 #include "gromacs/legacyheaders/types/simple.h"
85 #include "gromacs/legacyheaders/types/state.h"
86 #include "gromacs/listed-forces/manage-threading.h"
87 #include "gromacs/math/vec.h"
88 #include "gromacs/math/vectypes.h"
89 #include "gromacs/mdlib/nb_verlet.h"
90 #include "gromacs/mdlib/nbnxn_search.h"
91 #include "gromacs/pbcutil/ishift.h"
92 #include "gromacs/pbcutil/pbc.h"
93 #include "gromacs/pulling/pull.h"
94 #include "gromacs/pulling/pull_rotation.h"
95 #include "gromacs/swap/swapcoords.h"
96 #include "gromacs/timing/wallcycle.h"
97 #include "gromacs/topology/block.h"
98 #include "gromacs/topology/idef.h"
99 #include "gromacs/topology/mtop_util.h"
100 #include "gromacs/topology/topology.h"
101 #include "gromacs/utility/basedefinitions.h"
102 #include "gromacs/utility/basenetwork.h"
103 #include "gromacs/utility/cstringutil.h"
104 #include "gromacs/utility/fatalerror.h"
105 #include "gromacs/utility/gmxmpi.h"
106 #include "gromacs/utility/qsort_threadsafe.h"
107 #include "gromacs/utility/real.h"
108 #include "gromacs/utility/smalloc.h"
110 #include "domdec_constraints.h"
111 #include "domdec_internal.h"
112 #include "domdec_vsite.h"
114 #define DDRANK(dd, rank) (rank)
115 #define DDMASTERRANK(dd) (dd->masterrank)
117 typedef struct gmx_domdec_master
119 /* The cell boundaries */
121 /* The global charge group division */
122 int *ncg; /* Number of home charge groups for each node */
123 int *index; /* Index of nnodes+1 into cg */
124 int *cg; /* Global charge group index */
125 int *nat; /* Number of home atoms for each node. */
126 int *ibuf; /* Buffer for communication */
127 rvec *vbuf; /* Buffer for state scattering and gathering */
128 } gmx_domdec_master_t;
132 /* The numbers of charge groups to send and receive for each cell
133 * that requires communication, the last entry contains the total
134 * number of atoms that needs to be communicated.
136 int nsend[DD_MAXIZONE+2];
137 int nrecv[DD_MAXIZONE+2];
138 /* The charge groups to send */
141 /* The atom range for non-in-place communication */
142 int cell2at0[DD_MAXIZONE];
143 int cell2at1[DD_MAXIZONE];
148 int np; /* Number of grid pulses in this dimension */
149 int np_dlb; /* For dlb, for use with edlbAUTO */
150 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
152 gmx_bool bInPlace; /* Can we communicate in place? */
153 } gmx_domdec_comm_dim_t;
157 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
158 real *cell_f; /* State var.: cell boundaries, box relative */
159 real *old_cell_f; /* Temp. var.: old cell size */
160 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
161 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
162 real *bound_min; /* Temp. var.: lower limit for cell boundary */
163 real *bound_max; /* Temp. var.: upper limit for cell boundary */
164 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
165 real *buf_ncd; /* Temp. var. */
168 #define DD_NLOAD_MAX 9
170 /* Here floats are accurate enough, since these variables
171 * only influence the load balancing, not the actual MD results.
198 gmx_cgsort_t *sort_new;
210 /* This enum determines the order of the coordinates.
211 * ddnatHOME and ddnatZONE should be first and second,
212 * the others can be ordered as wanted.
215 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
219 edlbAUTO, edlbNO, edlbYES, edlbNR
221 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
225 int dim; /* The dimension */
226 gmx_bool dim_match; /* Tells if DD and PME dims match */
227 int nslab; /* The number of PME slabs in this dimension */
228 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
229 int *pp_min; /* The minimum pp node location, size nslab */
230 int *pp_max; /* The maximum pp node location,size nslab */
231 int maxshift; /* The maximum shift for coordinate redistribution in PME */
236 real min0; /* The minimum bottom of this zone */
237 real max1; /* The maximum top of this zone */
238 real min1; /* The minimum top of this zone */
239 real mch0; /* The maximum bottom communicaton height for this zone */
240 real mch1; /* The maximum top communicaton height for this zone */
241 real p1_0; /* The bottom value of the first cell in this zone */
242 real p1_1; /* The top value of the first cell in this zone */
247 gmx_domdec_ind_t ind;
254 } dd_comm_setup_work_t;
256 typedef struct gmx_domdec_comm
258 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
259 * unless stated otherwise.
262 /* The number of decomposition dimensions for PME, 0: no PME */
264 /* The number of nodes doing PME (PP/PME or only PME) */
268 /* The communication setup including the PME only nodes */
269 gmx_bool bCartesianPP_PME;
272 int *pmenodes; /* size npmenodes */
273 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
274 * but with bCartesianPP_PME */
275 gmx_ddpme_t ddpme[2];
277 /* The DD particle-particle nodes only */
278 gmx_bool bCartesianPP;
279 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
281 /* The global charge groups */
284 /* Should we sort the cgs */
286 gmx_domdec_sort_t *sort;
288 /* Are there charge groups? */
291 /* Are there bonded and multi-body interactions between charge groups? */
292 gmx_bool bInterCGBondeds;
293 gmx_bool bInterCGMultiBody;
295 /* Data for the optional bonded interaction atom communication range */
302 /* Is eDLB=edlbAUTO locked such that we currently can't turn it on? */
303 gmx_bool bDLB_locked;
304 /* With eDLB=edlbAUTO, should we check if to DLB on at the next DD? */
305 gmx_bool bCheckWhetherToTurnDlbOn;
306 /* Are we actually using DLB? */
307 gmx_bool bDynLoadBal;
309 /* Cell sizes for static load balancing, first index cartesian */
312 /* The width of the communicated boundaries */
315 /* The minimum cell size (including triclinic correction) */
317 /* For dlb, for use with edlbAUTO */
318 rvec cellsize_min_dlb;
319 /* The lower limit for the DD cell size with DLB */
321 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
322 gmx_bool bVacDLBNoLimit;
324 /* With PME load balancing we set limits on DLB */
325 gmx_bool bPMELoadBalDLBLimits;
326 /* DLB needs to take into account that we want to allow this maximum
327 * cut-off (for PME load balancing), this could limit cell boundaries.
329 real PMELoadBal_max_cutoff;
331 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
333 /* box0 and box_size are required with dim's without pbc and -gcom */
337 /* The cell boundaries */
341 /* The old location of the cell boundaries, to check cg displacements */
345 /* The communication setup and charge group boundaries for the zones */
346 gmx_domdec_zones_t zones;
348 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
349 * cell boundaries of neighboring cells for dynamic load balancing.
351 gmx_ddzone_t zone_d1[2];
352 gmx_ddzone_t zone_d2[2][2];
354 /* The coordinate/force communication setup and indices */
355 gmx_domdec_comm_dim_t cd[DIM];
356 /* The maximum number of cells to communicate with in one dimension */
359 /* Which cg distribution is stored on the master node */
360 int master_cg_ddp_count;
362 /* The number of cg's received from the direct neighbors */
363 int zone_ncg1[DD_MAXZONE];
365 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
368 /* Array for signalling if atoms have moved to another domain */
372 /* Communication buffer for general use */
376 /* Communication buffer for general use */
379 /* Temporary storage for thread parallel communication setup */
381 dd_comm_setup_work_t *dth;
383 /* Communication buffers only used with multiple grid pulses */
388 /* Communication buffers for local redistribution */
390 int cggl_flag_nalloc[DIM*2];
392 int cgcm_state_nalloc[DIM*2];
394 /* Cell sizes for dynamic load balancing */
395 gmx_domdec_root_t **root;
399 real cell_f_max0[DIM];
400 real cell_f_min1[DIM];
402 /* Stuff for load communication */
403 gmx_bool bRecordLoad;
404 gmx_domdec_load_t *load;
405 int nrank_gpu_shared;
407 MPI_Comm *mpi_comm_load;
408 MPI_Comm mpi_comm_gpu_shared;
411 /* Maximum DLB scaling per load balancing step in percent */
415 float cycl[ddCyclNr];
416 int cycl_n[ddCyclNr];
417 float cycl_max[ddCyclNr];
418 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
422 /* How many times have did we have load measurements */
424 /* How many times have we collected the load measurements */
428 double sum_nat[ddnatNR-ddnatZONE];
438 /* The last partition step */
439 gmx_int64_t partition_step;
447 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
450 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
451 #define DD_FLAG_NRCG 65535
452 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
453 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
455 /* Zone permutation required to obtain consecutive charge groups
456 * for neighbor searching.
458 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
460 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
461 * components see only j zones with that component 0.
464 /* The DD zone order */
465 static const ivec dd_zo[DD_MAXZONE] =
466 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
471 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
476 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
481 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
483 /* Factors used to avoid problems due to rounding issues */
484 #define DD_CELL_MARGIN 1.0001
485 #define DD_CELL_MARGIN2 1.00005
486 /* Factor to account for pressure scaling during nstlist steps */
487 #define DD_PRES_SCALE_MARGIN 1.02
489 /* Turn on DLB when the load imbalance causes this amount of total loss.
490 * There is a bit of overhead with DLB and it's difficult to achieve
491 * a load imbalance of less than 2% with DLB.
493 #define DD_PERF_LOSS_DLB_ON 0.02
495 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
496 #define DD_PERF_LOSS_WARN 0.05
498 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
500 /* Use separate MPI send and receive commands
501 * when nnodes <= GMX_DD_NNODES_SENDRECV.
502 * This saves memory (and some copying for small nnodes).
503 * For high parallelization scatter and gather calls are used.
505 #define GMX_DD_NNODES_SENDRECV 4
509 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
511 static void index2xyz(ivec nc,int ind,ivec xyz)
513 xyz[XX] = ind % nc[XX];
514 xyz[YY] = (ind / nc[XX]) % nc[YY];
515 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
519 /* This order is required to minimize the coordinate communication in PME
520 * which uses decomposition in the x direction.
522 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
524 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
526 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
527 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
528 xyz[ZZ] = ind % nc[ZZ];
531 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
536 ddindex = dd_index(dd->nc, c);
537 if (dd->comm->bCartesianPP_PME)
539 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
541 else if (dd->comm->bCartesianPP)
544 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
555 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
557 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
560 int ddglatnr(gmx_domdec_t *dd, int i)
570 if (i >= dd->comm->nat[ddnatNR-1])
572 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
574 atnr = dd->gatindex[i] + 1;
580 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
582 return &dd->comm->cgs_gl;
585 static void vec_rvec_init(vec_rvec_t *v)
591 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
595 v->nalloc = over_alloc_dd(n);
596 srenew(v->v, v->nalloc);
600 void dd_store_state(gmx_domdec_t *dd, t_state *state)
604 if (state->ddp_count != dd->ddp_count)
606 gmx_incons("The state does not the domain decomposition state");
609 state->ncg_gl = dd->ncg_home;
610 if (state->ncg_gl > state->cg_gl_nalloc)
612 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
613 srenew(state->cg_gl, state->cg_gl_nalloc);
615 for (i = 0; i < state->ncg_gl; i++)
617 state->cg_gl[i] = dd->index_gl[i];
620 state->ddp_count_cg_gl = dd->ddp_count;
623 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
625 return &dd->comm->zones;
628 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
629 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
631 gmx_domdec_zones_t *zones;
634 zones = &dd->comm->zones;
637 while (icg >= zones->izone[izone].cg1)
646 else if (izone < zones->nizone)
648 *jcg0 = zones->izone[izone].jcg0;
652 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
653 icg, izone, zones->nizone);
656 *jcg1 = zones->izone[izone].jcg1;
658 for (d = 0; d < dd->ndim; d++)
661 shift0[dim] = zones->izone[izone].shift0[dim];
662 shift1[dim] = zones->izone[izone].shift1[dim];
663 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
665 /* A conservative approach, this can be optimized */
672 int dd_natoms_vsite(gmx_domdec_t *dd)
674 return dd->comm->nat[ddnatVSITE];
677 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
679 *at_start = dd->comm->nat[ddnatCON-1];
680 *at_end = dd->comm->nat[ddnatCON];
683 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
685 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
686 int *index, *cgindex;
687 gmx_domdec_comm_t *comm;
688 gmx_domdec_comm_dim_t *cd;
689 gmx_domdec_ind_t *ind;
690 rvec shift = {0, 0, 0}, *buf, *rbuf;
691 gmx_bool bPBC, bScrew;
695 cgindex = dd->cgindex;
700 nat_tot = dd->nat_home;
701 for (d = 0; d < dd->ndim; d++)
703 bPBC = (dd->ci[dd->dim[d]] == 0);
704 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
707 copy_rvec(box[dd->dim[d]], shift);
710 for (p = 0; p < cd->np; p++)
717 for (i = 0; i < ind->nsend[nzone]; i++)
719 at0 = cgindex[index[i]];
720 at1 = cgindex[index[i]+1];
721 for (j = at0; j < at1; j++)
723 copy_rvec(x[j], buf[n]);
730 for (i = 0; i < ind->nsend[nzone]; i++)
732 at0 = cgindex[index[i]];
733 at1 = cgindex[index[i]+1];
734 for (j = at0; j < at1; j++)
736 /* We need to shift the coordinates */
737 rvec_add(x[j], shift, buf[n]);
744 for (i = 0; i < ind->nsend[nzone]; i++)
746 at0 = cgindex[index[i]];
747 at1 = cgindex[index[i]+1];
748 for (j = at0; j < at1; j++)
751 buf[n][XX] = x[j][XX] + shift[XX];
753 * This operation requires a special shift force
754 * treatment, which is performed in calc_vir.
756 buf[n][YY] = box[YY][YY] - x[j][YY];
757 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
769 rbuf = comm->vbuf2.v;
771 /* Send and receive the coordinates */
772 dd_sendrecv_rvec(dd, d, dddirBackward,
773 buf, ind->nsend[nzone+1],
774 rbuf, ind->nrecv[nzone+1]);
778 for (zone = 0; zone < nzone; zone++)
780 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
782 copy_rvec(rbuf[j], x[i]);
787 nat_tot += ind->nrecv[nzone+1];
793 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
795 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
796 int *index, *cgindex;
797 gmx_domdec_comm_t *comm;
798 gmx_domdec_comm_dim_t *cd;
799 gmx_domdec_ind_t *ind;
803 gmx_bool bShiftForcesNeedPbc, bScrew;
807 cgindex = dd->cgindex;
811 nzone = comm->zones.n/2;
812 nat_tot = dd->nat_tot;
813 for (d = dd->ndim-1; d >= 0; d--)
815 /* Only forces in domains near the PBC boundaries need to
816 consider PBC in the treatment of fshift */
817 bShiftForcesNeedPbc = (dd->ci[dd->dim[d]] == 0);
818 bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
819 if (fshift == NULL && !bScrew)
821 bShiftForcesNeedPbc = FALSE;
823 /* Determine which shift vector we need */
829 for (p = cd->np-1; p >= 0; p--)
832 nat_tot -= ind->nrecv[nzone+1];
839 sbuf = comm->vbuf2.v;
841 for (zone = 0; zone < nzone; zone++)
843 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
845 copy_rvec(f[i], sbuf[j]);
850 /* Communicate the forces */
851 dd_sendrecv_rvec(dd, d, dddirForward,
852 sbuf, ind->nrecv[nzone+1],
853 buf, ind->nsend[nzone+1]);
855 /* Add the received forces */
857 if (!bShiftForcesNeedPbc)
859 for (i = 0; i < ind->nsend[nzone]; i++)
861 at0 = cgindex[index[i]];
862 at1 = cgindex[index[i]+1];
863 for (j = at0; j < at1; j++)
865 rvec_inc(f[j], buf[n]);
872 /* fshift should always be defined if this function is
873 * called when bShiftForcesNeedPbc is true */
874 assert(NULL != fshift);
875 for (i = 0; i < ind->nsend[nzone]; i++)
877 at0 = cgindex[index[i]];
878 at1 = cgindex[index[i]+1];
879 for (j = at0; j < at1; j++)
881 rvec_inc(f[j], buf[n]);
882 /* Add this force to the shift force */
883 rvec_inc(fshift[is], buf[n]);
890 for (i = 0; i < ind->nsend[nzone]; i++)
892 at0 = cgindex[index[i]];
893 at1 = cgindex[index[i]+1];
894 for (j = at0; j < at1; j++)
896 /* Rotate the force */
897 f[j][XX] += buf[n][XX];
898 f[j][YY] -= buf[n][YY];
899 f[j][ZZ] -= buf[n][ZZ];
902 /* Add this force to the shift force */
903 rvec_inc(fshift[is], buf[n]);
914 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
916 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
917 int *index, *cgindex;
918 gmx_domdec_comm_t *comm;
919 gmx_domdec_comm_dim_t *cd;
920 gmx_domdec_ind_t *ind;
925 cgindex = dd->cgindex;
927 buf = &comm->vbuf.v[0][0];
930 nat_tot = dd->nat_home;
931 for (d = 0; d < dd->ndim; d++)
934 for (p = 0; p < cd->np; p++)
939 for (i = 0; i < ind->nsend[nzone]; i++)
941 at0 = cgindex[index[i]];
942 at1 = cgindex[index[i]+1];
943 for (j = at0; j < at1; j++)
956 rbuf = &comm->vbuf2.v[0][0];
958 /* Send and receive the coordinates */
959 dd_sendrecv_real(dd, d, dddirBackward,
960 buf, ind->nsend[nzone+1],
961 rbuf, ind->nrecv[nzone+1]);
965 for (zone = 0; zone < nzone; zone++)
967 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
974 nat_tot += ind->nrecv[nzone+1];
980 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
982 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
983 int *index, *cgindex;
984 gmx_domdec_comm_t *comm;
985 gmx_domdec_comm_dim_t *cd;
986 gmx_domdec_ind_t *ind;
991 cgindex = dd->cgindex;
993 buf = &comm->vbuf.v[0][0];
995 nzone = comm->zones.n/2;
996 nat_tot = dd->nat_tot;
997 for (d = dd->ndim-1; d >= 0; d--)
1000 for (p = cd->np-1; p >= 0; p--)
1003 nat_tot -= ind->nrecv[nzone+1];
1010 sbuf = &comm->vbuf2.v[0][0];
1012 for (zone = 0; zone < nzone; zone++)
1014 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
1021 /* Communicate the forces */
1022 dd_sendrecv_real(dd, d, dddirForward,
1023 sbuf, ind->nrecv[nzone+1],
1024 buf, ind->nsend[nzone+1]);
1026 /* Add the received forces */
1028 for (i = 0; i < ind->nsend[nzone]; i++)
1030 at0 = cgindex[index[i]];
1031 at1 = cgindex[index[i]+1];
1032 for (j = at0; j < at1; j++)
1043 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1045 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",
1047 zone->min0, zone->max1,
1048 zone->mch0, zone->mch0,
1049 zone->p1_0, zone->p1_1);
1053 #define DDZONECOMM_MAXZONE 5
1054 #define DDZONECOMM_BUFSIZE 3
1056 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1057 int ddimind, int direction,
1058 gmx_ddzone_t *buf_s, int n_s,
1059 gmx_ddzone_t *buf_r, int n_r)
1061 #define ZBS DDZONECOMM_BUFSIZE
1062 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1063 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1066 for (i = 0; i < n_s; i++)
1068 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1069 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1070 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1071 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1072 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1073 vbuf_s[i*ZBS+1][2] = 0;
1074 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1075 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1076 vbuf_s[i*ZBS+2][2] = 0;
1079 dd_sendrecv_rvec(dd, ddimind, direction,
1083 for (i = 0; i < n_r; i++)
1085 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1086 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1087 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1088 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1089 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1090 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1091 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1097 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1098 rvec cell_ns_x0, rvec cell_ns_x1)
1100 int d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
1102 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1103 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1104 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1105 rvec extr_s[2], extr_r[2];
1107 real dist_d, c = 0, det;
1108 gmx_domdec_comm_t *comm;
1109 gmx_bool bPBC, bUse;
1113 for (d = 1; d < dd->ndim; d++)
1116 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1117 zp->min0 = cell_ns_x0[dim];
1118 zp->max1 = cell_ns_x1[dim];
1119 zp->min1 = cell_ns_x1[dim];
1120 zp->mch0 = cell_ns_x0[dim];
1121 zp->mch1 = cell_ns_x1[dim];
1122 zp->p1_0 = cell_ns_x0[dim];
1123 zp->p1_1 = cell_ns_x1[dim];
1126 for (d = dd->ndim-2; d >= 0; d--)
1129 bPBC = (dim < ddbox->npbcdim);
1131 /* Use an rvec to store two reals */
1132 extr_s[d][0] = comm->cell_f0[d+1];
1133 extr_s[d][1] = comm->cell_f1[d+1];
1134 extr_s[d][2] = comm->cell_f1[d+1];
1137 /* Store the extremes in the backward sending buffer,
1138 * so the get updated separately from the forward communication.
1140 for (d1 = d; d1 < dd->ndim-1; d1++)
1142 /* We invert the order to be able to use the same loop for buf_e */
1143 buf_s[pos].min0 = extr_s[d1][1];
1144 buf_s[pos].max1 = extr_s[d1][0];
1145 buf_s[pos].min1 = extr_s[d1][2];
1146 buf_s[pos].mch0 = 0;
1147 buf_s[pos].mch1 = 0;
1148 /* Store the cell corner of the dimension we communicate along */
1149 buf_s[pos].p1_0 = comm->cell_x0[dim];
1150 buf_s[pos].p1_1 = 0;
1154 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1157 if (dd->ndim == 3 && d == 0)
1159 buf_s[pos] = comm->zone_d2[0][1];
1161 buf_s[pos] = comm->zone_d1[0];
1165 /* We only need to communicate the extremes
1166 * in the forward direction
1168 npulse = comm->cd[d].np;
1171 /* Take the minimum to avoid double communication */
1172 npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
1176 /* Without PBC we should really not communicate over
1177 * the boundaries, but implementing that complicates
1178 * the communication setup and therefore we simply
1179 * do all communication, but ignore some data.
1181 npulse_min = npulse;
1183 for (p = 0; p < npulse_min; p++)
1185 /* Communicate the extremes forward */
1186 bUse = (bPBC || dd->ci[dim] > 0);
1188 dd_sendrecv_rvec(dd, d, dddirForward,
1189 extr_s+d, dd->ndim-d-1,
1190 extr_r+d, dd->ndim-d-1);
1194 for (d1 = d; d1 < dd->ndim-1; d1++)
1196 extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
1197 extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
1198 extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
1204 for (p = 0; p < npulse; p++)
1206 /* Communicate all the zone information backward */
1207 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1209 dd_sendrecv_ddzone(dd, d, dddirBackward,
1216 for (d1 = d+1; d1 < dd->ndim; d1++)
1218 /* Determine the decrease of maximum required
1219 * communication height along d1 due to the distance along d,
1220 * this avoids a lot of useless atom communication.
1222 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1224 if (ddbox->tric_dir[dim])
1226 /* c is the off-diagonal coupling between the cell planes
1227 * along directions d and d1.
1229 c = ddbox->v[dim][dd->dim[d1]][dim];
1235 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1238 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1242 /* A negative value signals out of range */
1248 /* Accumulate the extremes over all pulses */
1249 for (i = 0; i < buf_size; i++)
1253 buf_e[i] = buf_r[i];
1259 buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
1260 buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
1261 buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
1264 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1272 if (bUse && dh[d1] >= 0)
1274 buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1275 buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1278 /* Copy the received buffer to the send buffer,
1279 * to pass the data through with the next pulse.
1281 buf_s[i] = buf_r[i];
1283 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1284 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1286 /* Store the extremes */
1289 for (d1 = d; d1 < dd->ndim-1; d1++)
1291 extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
1292 extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
1293 extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
1297 if (d == 1 || (d == 0 && dd->ndim == 3))
1299 for (i = d; i < 2; i++)
1301 comm->zone_d2[1-d][i] = buf_e[pos];
1307 comm->zone_d1[1] = buf_e[pos];
1317 for (i = 0; i < 2; i++)
1321 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1323 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1324 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1330 for (i = 0; i < 2; i++)
1332 for (j = 0; j < 2; j++)
1336 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1338 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1339 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1343 for (d = 1; d < dd->ndim; d++)
1345 comm->cell_f_max0[d] = extr_s[d-1][0];
1346 comm->cell_f_min1[d] = extr_s[d-1][1];
1349 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1350 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1355 static void dd_collect_cg(gmx_domdec_t *dd,
1356 t_state *state_local)
1358 gmx_domdec_master_t *ma = NULL;
1359 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1361 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1363 /* The master has the correct distribution */
1367 if (state_local->ddp_count == dd->ddp_count)
1369 /* The local state and DD are in sync, use the DD indices */
1370 ncg_home = dd->ncg_home;
1372 nat_home = dd->nat_home;
1374 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1376 /* The DD is out of sync with the local state, but we have stored
1377 * the cg indices with the local state, so we can use those.
1381 cgs_gl = &dd->comm->cgs_gl;
1383 ncg_home = state_local->ncg_gl;
1384 cg = state_local->cg_gl;
1386 for (i = 0; i < ncg_home; i++)
1388 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1393 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1407 /* Collect the charge group and atom counts on the master */
1408 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1413 for (i = 0; i < dd->nnodes; i++)
1415 ma->ncg[i] = ma->ibuf[2*i];
1416 ma->nat[i] = ma->ibuf[2*i+1];
1417 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1420 /* Make byte counts and indices */
1421 for (i = 0; i < dd->nnodes; i++)
1423 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1424 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1428 fprintf(debug, "Initial charge group distribution: ");
1429 for (i = 0; i < dd->nnodes; i++)
1431 fprintf(debug, " %d", ma->ncg[i]);
1433 fprintf(debug, "\n");
1437 /* Collect the charge group indices on the master */
1439 ncg_home*sizeof(int), cg,
1440 DDMASTER(dd) ? ma->ibuf : NULL,
1441 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1442 DDMASTER(dd) ? ma->cg : NULL);
1444 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1447 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1450 gmx_domdec_master_t *ma;
1451 int n, i, c, a, nalloc = 0;
1460 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1461 dd->rank, dd->mpi_comm_all);
1466 /* Copy the master coordinates to the global array */
1467 cgs_gl = &dd->comm->cgs_gl;
1469 n = DDMASTERRANK(dd);
1471 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1473 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1475 copy_rvec(lv[a++], v[c]);
1479 for (n = 0; n < dd->nnodes; n++)
1483 if (ma->nat[n] > nalloc)
1485 nalloc = over_alloc_dd(ma->nat[n]);
1486 srenew(buf, nalloc);
1489 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1490 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1493 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1495 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1497 copy_rvec(buf[a++], v[c]);
1506 static void get_commbuffer_counts(gmx_domdec_t *dd,
1507 int **counts, int **disps)
1509 gmx_domdec_master_t *ma;
1514 /* Make the rvec count and displacment arrays */
1516 *disps = ma->ibuf + dd->nnodes;
1517 for (n = 0; n < dd->nnodes; n++)
1519 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1520 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1524 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1527 gmx_domdec_master_t *ma;
1528 int *rcounts = NULL, *disps = NULL;
1537 get_commbuffer_counts(dd, &rcounts, &disps);
1542 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1546 cgs_gl = &dd->comm->cgs_gl;
1549 for (n = 0; n < dd->nnodes; n++)
1551 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1553 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1555 copy_rvec(buf[a++], v[c]);
1562 void dd_collect_vec(gmx_domdec_t *dd,
1563 t_state *state_local, rvec *lv, rvec *v)
1565 dd_collect_cg(dd, state_local);
1567 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1569 dd_collect_vec_sendrecv(dd, lv, v);
1573 dd_collect_vec_gatherv(dd, lv, v);
1578 void dd_collect_state(gmx_domdec_t *dd,
1579 t_state *state_local, t_state *state)
1583 nh = state->nhchainlength;
1587 for (i = 0; i < efptNR; i++)
1589 state->lambda[i] = state_local->lambda[i];
1591 state->fep_state = state_local->fep_state;
1592 state->veta = state_local->veta;
1593 state->vol0 = state_local->vol0;
1594 copy_mat(state_local->box, state->box);
1595 copy_mat(state_local->boxv, state->boxv);
1596 copy_mat(state_local->svir_prev, state->svir_prev);
1597 copy_mat(state_local->fvir_prev, state->fvir_prev);
1598 copy_mat(state_local->pres_prev, state->pres_prev);
1600 for (i = 0; i < state_local->ngtc; i++)
1602 for (j = 0; j < nh; j++)
1604 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1605 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1607 state->therm_integral[i] = state_local->therm_integral[i];
1609 for (i = 0; i < state_local->nnhpres; i++)
1611 for (j = 0; j < nh; j++)
1613 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1614 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1618 for (est = 0; est < estNR; est++)
1620 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1625 dd_collect_vec(dd, state_local, state_local->x, state->x);
1628 dd_collect_vec(dd, state_local, state_local->v, state->v);
1631 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1634 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1636 case estDISRE_INITF:
1637 case estDISRE_RM3TAV:
1638 case estORIRE_INITF:
1642 gmx_incons("Unknown state entry encountered in dd_collect_state");
1648 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1654 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1657 state->nalloc = over_alloc_dd(nalloc);
1659 for (est = 0; est < estNR; est++)
1661 if (EST_DISTR(est) && (state->flags & (1<<est)))
1666 srenew(state->x, state->nalloc);
1669 srenew(state->v, state->nalloc);
1672 srenew(state->sd_X, state->nalloc);
1675 srenew(state->cg_p, state->nalloc);
1677 case estDISRE_INITF:
1678 case estDISRE_RM3TAV:
1679 case estORIRE_INITF:
1681 /* No reallocation required */
1684 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1691 srenew(*f, state->nalloc);
1695 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1698 if (nalloc > fr->cg_nalloc)
1702 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1704 fr->cg_nalloc = over_alloc_dd(nalloc);
1705 srenew(fr->cginfo, fr->cg_nalloc);
1706 if (fr->cutoff_scheme == ecutsGROUP)
1708 srenew(fr->cg_cm, fr->cg_nalloc);
1711 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1713 /* We don't use charge groups, we use x in state to set up
1714 * the atom communication.
1716 dd_realloc_state(state, f, nalloc);
1720 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1723 gmx_domdec_master_t *ma;
1724 int n, i, c, a, nalloc = 0;
1731 for (n = 0; n < dd->nnodes; n++)
1735 if (ma->nat[n] > nalloc)
1737 nalloc = over_alloc_dd(ma->nat[n]);
1738 srenew(buf, nalloc);
1740 /* Use lv as a temporary buffer */
1742 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1744 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1746 copy_rvec(v[c], buf[a++]);
1749 if (a != ma->nat[n])
1751 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1756 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1757 DDRANK(dd, n), n, dd->mpi_comm_all);
1762 n = DDMASTERRANK(dd);
1764 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1766 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1768 copy_rvec(v[c], lv[a++]);
1775 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1776 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1781 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1784 gmx_domdec_master_t *ma;
1785 int *scounts = NULL, *disps = NULL;
1793 get_commbuffer_counts(dd, &scounts, &disps);
1797 for (n = 0; n < dd->nnodes; n++)
1799 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1801 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1803 copy_rvec(v[c], buf[a++]);
1809 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1812 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1814 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1816 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1820 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1824 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1827 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1828 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1829 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1831 if (dfhist->nlambda > 0)
1833 int nlam = dfhist->nlambda;
1834 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1835 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1836 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1837 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1838 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1839 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1841 for (i = 0; i < nlam; i++)
1843 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1844 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1845 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1846 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1847 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1848 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1853 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1854 t_state *state, t_state *state_local,
1859 nh = state->nhchainlength;
1863 for (i = 0; i < efptNR; i++)
1865 state_local->lambda[i] = state->lambda[i];
1867 state_local->fep_state = state->fep_state;
1868 state_local->veta = state->veta;
1869 state_local->vol0 = state->vol0;
1870 copy_mat(state->box, state_local->box);
1871 copy_mat(state->box_rel, state_local->box_rel);
1872 copy_mat(state->boxv, state_local->boxv);
1873 copy_mat(state->svir_prev, state_local->svir_prev);
1874 copy_mat(state->fvir_prev, state_local->fvir_prev);
1875 copy_df_history(&state_local->dfhist, &state->dfhist);
1876 for (i = 0; i < state_local->ngtc; i++)
1878 for (j = 0; j < nh; j++)
1880 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1881 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1883 state_local->therm_integral[i] = state->therm_integral[i];
1885 for (i = 0; i < state_local->nnhpres; i++)
1887 for (j = 0; j < nh; j++)
1889 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1890 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1894 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1895 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1896 dd_bcast(dd, sizeof(real), &state_local->veta);
1897 dd_bcast(dd, sizeof(real), &state_local->vol0);
1898 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1899 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1900 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1901 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1902 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1903 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1904 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1905 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1906 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1907 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1909 /* communicate df_history -- required for restarting from checkpoint */
1910 dd_distribute_dfhist(dd, &state_local->dfhist);
1912 if (dd->nat_home > state_local->nalloc)
1914 dd_realloc_state(state_local, f, dd->nat_home);
1916 for (i = 0; i < estNR; i++)
1918 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1923 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1926 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1929 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1932 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1934 case estDISRE_INITF:
1935 case estDISRE_RM3TAV:
1936 case estORIRE_INITF:
1938 /* Not implemented yet */
1941 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1947 static char dim2char(int dim)
1953 case XX: c = 'X'; break;
1954 case YY: c = 'Y'; break;
1955 case ZZ: c = 'Z'; break;
1956 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1962 static void write_dd_grid_pdb(const char *fn, gmx_int64_t step,
1963 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1965 rvec grid_s[2], *grid_r = NULL, cx, r;
1966 char fname[STRLEN], buf[22];
1968 int a, i, d, z, y, x;
1972 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1973 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1977 snew(grid_r, 2*dd->nnodes);
1980 dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
1984 for (d = 0; d < DIM; d++)
1986 for (i = 0; i < DIM; i++)
1994 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1996 tric[d][i] = box[i][d]/box[i][i];
2005 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2006 out = gmx_fio_fopen(fname, "w");
2007 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2009 for (i = 0; i < dd->nnodes; i++)
2011 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2012 for (d = 0; d < DIM; d++)
2014 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2016 for (z = 0; z < 2; z++)
2018 for (y = 0; y < 2; y++)
2020 for (x = 0; x < 2; x++)
2022 cx[XX] = grid_r[i*2+x][XX];
2023 cx[YY] = grid_r[i*2+y][YY];
2024 cx[ZZ] = grid_r[i*2+z][ZZ];
2026 gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
2027 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
2031 for (d = 0; d < DIM; d++)
2033 for (x = 0; x < 4; x++)
2037 case 0: y = 1 + i*8 + 2*x; break;
2038 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2039 case 2: y = 1 + i*8 + x; break;
2041 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2045 gmx_fio_fclose(out);
2050 void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
2051 gmx_mtop_t *mtop, t_commrec *cr,
2052 int natoms, rvec x[], matrix box)
2054 char fname[STRLEN], buf[22];
2056 int i, ii, resnr, c;
2057 char *atomname, *resname;
2064 natoms = dd->comm->nat[ddnatVSITE];
2067 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2069 out = gmx_fio_fopen(fname, "w");
2071 fprintf(out, "TITLE %s\n", title);
2072 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2073 for (i = 0; i < natoms; i++)
2075 ii = dd->gatindex[i];
2076 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2077 if (i < dd->comm->nat[ddnatZONE])
2080 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2086 else if (i < dd->comm->nat[ddnatVSITE])
2088 b = dd->comm->zones.n;
2092 b = dd->comm->zones.n + 1;
2094 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
2095 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
2097 fprintf(out, "TER\n");
2099 gmx_fio_fclose(out);
2102 real dd_cutoff_multibody(const gmx_domdec_t *dd)
2104 gmx_domdec_comm_t *comm;
2111 if (comm->bInterCGBondeds)
2113 if (comm->cutoff_mbody > 0)
2115 r = comm->cutoff_mbody;
2119 /* cutoff_mbody=0 means we do not have DLB */
2120 r = comm->cellsize_min[dd->dim[0]];
2121 for (di = 1; di < dd->ndim; di++)
2123 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
2125 if (comm->bBondComm)
2127 r = std::max(r, comm->cutoff_mbody);
2131 r = std::min(r, comm->cutoff);
2139 real dd_cutoff_twobody(const gmx_domdec_t *dd)
2143 r_mb = dd_cutoff_multibody(dd);
2145 return std::max(dd->comm->cutoff, r_mb);
2149 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2153 nc = dd->nc[dd->comm->cartpmedim];
2154 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2155 copy_ivec(coord, coord_pme);
2156 coord_pme[dd->comm->cartpmedim] =
2157 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2160 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2162 /* Here we assign a PME node to communicate with this DD node
2163 * by assuming that the major index of both is x.
2164 * We add cr->npmenodes/2 to obtain an even distribution.
2166 return (ddindex*npme + npme/2)/ndd;
2169 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2171 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2174 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2176 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2179 static int *dd_pmenodes(t_commrec *cr)
2184 snew(pmenodes, cr->npmenodes);
2186 for (i = 0; i < cr->dd->nnodes; i++)
2188 p0 = cr_ddindex2pmeindex(cr, i);
2189 p1 = cr_ddindex2pmeindex(cr, i+1);
2190 if (i+1 == cr->dd->nnodes || p1 > p0)
2194 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2196 pmenodes[n] = i + 1 + n;
2204 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2212 if (dd->comm->bCartesian) {
2213 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2214 dd_coords2pmecoords(dd,coords,coords_pme);
2215 copy_ivec(dd->ntot,nc);
2216 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2217 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2219 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2221 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2227 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2232 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2234 gmx_domdec_comm_t *comm;
2236 int ddindex, nodeid = -1;
2238 comm = cr->dd->comm;
2243 if (comm->bCartesianPP_PME)
2246 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2251 ddindex = dd_index(cr->dd->nc, coords);
2252 if (comm->bCartesianPP)
2254 nodeid = comm->ddindex2simnodeid[ddindex];
2260 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2272 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2275 gmx_domdec_comm_t *comm;
2282 /* This assumes a uniform x domain decomposition grid cell size */
2283 if (comm->bCartesianPP_PME)
2286 ivec coord, coord_pme;
2287 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2288 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2290 /* This is a PP node */
2291 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2292 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2296 else if (comm->bCartesianPP)
2298 if (sim_nodeid < dd->nnodes)
2300 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2305 /* This assumes DD cells with identical x coordinates
2306 * are numbered sequentially.
2308 if (dd->comm->pmenodes == NULL)
2310 if (sim_nodeid < dd->nnodes)
2312 /* The DD index equals the nodeid */
2313 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2319 while (sim_nodeid > dd->comm->pmenodes[i])
2323 if (sim_nodeid < dd->comm->pmenodes[i])
2325 pmenode = dd->comm->pmenodes[i];
2333 void get_pme_nnodes(const gmx_domdec_t *dd,
2334 int *npmenodes_x, int *npmenodes_y)
2338 *npmenodes_x = dd->comm->npmenodes_x;
2339 *npmenodes_y = dd->comm->npmenodes_y;
2348 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2349 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2353 ivec coord, coord_pme;
2357 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2360 for (x = 0; x < dd->nc[XX]; x++)
2362 for (y = 0; y < dd->nc[YY]; y++)
2364 for (z = 0; z < dd->nc[ZZ]; z++)
2366 if (dd->comm->bCartesianPP_PME)
2371 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2372 if (dd->ci[XX] == coord_pme[XX] &&
2373 dd->ci[YY] == coord_pme[YY] &&
2374 dd->ci[ZZ] == coord_pme[ZZ])
2376 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2381 /* The slab corresponds to the nodeid in the PME group */
2382 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2384 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2391 /* The last PP-only node is the peer node */
2392 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2396 fprintf(debug, "Receive coordinates from PP ranks:");
2397 for (x = 0; x < *nmy_ddnodes; x++)
2399 fprintf(debug, " %d", (*my_ddnodes)[x]);
2401 fprintf(debug, "\n");
2405 static gmx_bool receive_vir_ener(t_commrec *cr)
2407 gmx_domdec_comm_t *comm;
2412 if (cr->npmenodes < cr->dd->nnodes)
2414 comm = cr->dd->comm;
2415 if (comm->bCartesianPP_PME)
2417 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2420 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2421 coords[comm->cartpmedim]++;
2422 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2425 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2426 if (dd_simnode2pmenode(cr, rank) == pmenode)
2428 /* This is not the last PP node for pmenode */
2436 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2437 if (cr->sim_nodeid+1 < cr->nnodes &&
2438 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2440 /* This is not the last PP node for pmenode */
2449 static void set_zones_ncg_home(gmx_domdec_t *dd)
2451 gmx_domdec_zones_t *zones;
2454 zones = &dd->comm->zones;
2456 zones->cg_range[0] = 0;
2457 for (i = 1; i < zones->n+1; i++)
2459 zones->cg_range[i] = dd->ncg_home;
2461 /* zone_ncg1[0] should always be equal to ncg_home */
2462 dd->comm->zone_ncg1[0] = dd->ncg_home;
2465 static void rebuild_cgindex(gmx_domdec_t *dd,
2466 const int *gcgs_index, t_state *state)
2468 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2471 dd_cg_gl = dd->index_gl;
2472 cgindex = dd->cgindex;
2475 for (i = 0; i < state->ncg_gl; i++)
2479 dd_cg_gl[i] = cg_gl;
2480 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2484 dd->ncg_home = state->ncg_gl;
2487 set_zones_ncg_home(dd);
2490 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2492 while (cg >= cginfo_mb->cg_end)
2497 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2500 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2501 t_forcerec *fr, char *bLocalCG)
2503 cginfo_mb_t *cginfo_mb;
2509 cginfo_mb = fr->cginfo_mb;
2510 cginfo = fr->cginfo;
2512 for (cg = cg0; cg < cg1; cg++)
2514 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2518 if (bLocalCG != NULL)
2520 for (cg = cg0; cg < cg1; cg++)
2522 bLocalCG[index_gl[cg]] = TRUE;
2527 static void make_dd_indices(gmx_domdec_t *dd,
2528 const int *gcgs_index, int cg_start)
2530 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2531 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2534 if (dd->nat_tot > dd->gatindex_nalloc)
2536 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2537 srenew(dd->gatindex, dd->gatindex_nalloc);
2540 nzone = dd->comm->zones.n;
2541 zone2cg = dd->comm->zones.cg_range;
2542 zone_ncg1 = dd->comm->zone_ncg1;
2543 index_gl = dd->index_gl;
2544 gatindex = dd->gatindex;
2545 bCGs = dd->comm->bCGs;
2547 if (zone2cg[1] != dd->ncg_home)
2549 gmx_incons("dd->ncg_zone is not up to date");
2552 /* Make the local to global and global to local atom index */
2553 a = dd->cgindex[cg_start];
2554 for (zone = 0; zone < nzone; zone++)
2562 cg0 = zone2cg[zone];
2564 cg1 = zone2cg[zone+1];
2565 cg1_p1 = cg0 + zone_ncg1[zone];
2567 for (cg = cg0; cg < cg1; cg++)
2572 /* Signal that this cg is from more than one pulse away */
2575 cg_gl = index_gl[cg];
2578 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2581 ga2la_set(dd->ga2la, a_gl, a, zone1);
2587 gatindex[a] = cg_gl;
2588 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2595 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2601 if (bLocalCG == NULL)
2605 for (i = 0; i < dd->ncg_tot; i++)
2607 if (!bLocalCG[dd->index_gl[i]])
2610 "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);
2615 for (i = 0; i < ncg_sys; i++)
2622 if (ngl != dd->ncg_tot)
2624 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);
2631 static void check_index_consistency(gmx_domdec_t *dd,
2632 int natoms_sys, int ncg_sys,
2635 int nerr, ngl, i, a, cell;
2640 if (dd->comm->DD_debug > 1)
2642 snew(have, natoms_sys);
2643 for (a = 0; a < dd->nat_tot; a++)
2645 if (have[dd->gatindex[a]] > 0)
2647 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);
2651 have[dd->gatindex[a]] = a + 1;
2657 snew(have, dd->nat_tot);
2660 for (i = 0; i < natoms_sys; i++)
2662 if (ga2la_get(dd->ga2la, i, &a, &cell))
2664 if (a >= dd->nat_tot)
2666 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);
2672 if (dd->gatindex[a] != i)
2674 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);
2681 if (ngl != dd->nat_tot)
2684 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
2685 dd->rank, where, ngl, dd->nat_tot);
2687 for (a = 0; a < dd->nat_tot; a++)
2692 "DD rank %d, %s: local atom %d, global %d has no global index\n",
2693 dd->rank, where, a+1, dd->gatindex[a]+1);
2698 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2702 gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
2703 dd->rank, where, nerr);
2707 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2714 /* Clear the whole list without searching */
2715 ga2la_clear(dd->ga2la);
2719 for (i = a_start; i < dd->nat_tot; i++)
2721 ga2la_del(dd->ga2la, dd->gatindex[i]);
2725 bLocalCG = dd->comm->bLocalCG;
2728 for (i = cg_start; i < dd->ncg_tot; i++)
2730 bLocalCG[dd->index_gl[i]] = FALSE;
2734 dd_clear_local_vsite_indices(dd);
2736 if (dd->constraints)
2738 dd_clear_local_constraint_indices(dd);
2742 /* This function should be used for moving the domain boudaries during DLB,
2743 * for obtaining the minimum cell size. It checks the initially set limit
2744 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2745 * and, possibly, a longer cut-off limit set for PME load balancing.
2747 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2751 cellsize_min = comm->cellsize_min[dim];
2753 if (!comm->bVacDLBNoLimit)
2755 /* The cut-off might have changed, e.g. by PME load balacning,
2756 * from the value used to set comm->cellsize_min, so check it.
2758 cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2760 if (comm->bPMELoadBalDLBLimits)
2762 /* Check for the cut-off limit set by the PME load balancing */
2763 cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2767 return cellsize_min;
2770 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2773 real grid_jump_limit;
2775 /* The distance between the boundaries of cells at distance
2776 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2777 * and by the fact that cells should not be shifted by more than
2778 * half their size, such that cg's only shift by one cell
2779 * at redecomposition.
2781 grid_jump_limit = comm->cellsize_limit;
2782 if (!comm->bVacDLBNoLimit)
2784 if (comm->bPMELoadBalDLBLimits)
2786 cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
2788 grid_jump_limit = std::max(grid_jump_limit,
2789 cutoff/comm->cd[dim_ind].np);
2792 return grid_jump_limit;
2795 static gmx_bool check_grid_jump(gmx_int64_t step,
2801 gmx_domdec_comm_t *comm;
2810 for (d = 1; d < dd->ndim; d++)
2813 limit = grid_jump_limit(comm, cutoff, d);
2814 bfac = ddbox->box_size[dim];
2815 if (ddbox->tric_dir[dim])
2817 bfac *= ddbox->skew_fac[dim];
2819 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2820 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2828 /* This error should never be triggered under normal
2829 * circumstances, but you never know ...
2831 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.",
2832 gmx_step_str(step, buf),
2833 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2841 static int dd_load_count(gmx_domdec_comm_t *comm)
2843 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2846 static float dd_force_load(gmx_domdec_comm_t *comm)
2853 if (comm->eFlop > 1)
2855 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2860 load = comm->cycl[ddCyclF];
2861 if (comm->cycl_n[ddCyclF] > 1)
2863 /* Subtract the maximum of the last n cycle counts
2864 * to get rid of possible high counts due to other sources,
2865 * for instance system activity, that would otherwise
2866 * affect the dynamic load balancing.
2868 load -= comm->cycl_max[ddCyclF];
2872 if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
2874 float gpu_wait, gpu_wait_sum;
2876 gpu_wait = comm->cycl[ddCyclWaitGPU];
2877 if (comm->cycl_n[ddCyclF] > 1)
2879 /* We should remove the WaitGPU time of the same MD step
2880 * as the one with the maximum F time, since the F time
2881 * and the wait time are not independent.
2882 * Furthermore, the step for the max F time should be chosen
2883 * the same on all ranks that share the same GPU.
2884 * But to keep the code simple, we remove the average instead.
2885 * The main reason for artificially long times at some steps
2886 * is spurious CPU activity or MPI time, so we don't expect
2887 * that changes in the GPU wait time matter a lot here.
2889 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
2891 /* Sum the wait times over the ranks that share the same GPU */
2892 MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
2893 comm->mpi_comm_gpu_shared);
2894 /* Replace the wait time by the average over the ranks */
2895 load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
2903 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2905 gmx_domdec_comm_t *comm;
2910 snew(*dim_f, dd->nc[dim]+1);
2912 for (i = 1; i < dd->nc[dim]; i++)
2914 if (comm->slb_frac[dim])
2916 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2920 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2923 (*dim_f)[dd->nc[dim]] = 1;
2926 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2928 int pmeindex, slab, nso, i;
2931 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2937 ddpme->dim = dimind;
2939 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2941 ddpme->nslab = (ddpme->dim == 0 ?
2942 dd->comm->npmenodes_x :
2943 dd->comm->npmenodes_y);
2945 if (ddpme->nslab <= 1)
2950 nso = dd->comm->npmenodes/ddpme->nslab;
2951 /* Determine for each PME slab the PP location range for dimension dim */
2952 snew(ddpme->pp_min, ddpme->nslab);
2953 snew(ddpme->pp_max, ddpme->nslab);
2954 for (slab = 0; slab < ddpme->nslab; slab++)
2956 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2957 ddpme->pp_max[slab] = 0;
2959 for (i = 0; i < dd->nnodes; i++)
2961 ddindex2xyz(dd->nc, i, xyz);
2962 /* For y only use our y/z slab.
2963 * This assumes that the PME x grid size matches the DD grid size.
2965 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2967 pmeindex = ddindex2pmeindex(dd, i);
2970 slab = pmeindex/nso;
2974 slab = pmeindex % ddpme->nslab;
2976 ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
2977 ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
2981 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2984 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2986 if (dd->comm->ddpme[0].dim == XX)
2988 return dd->comm->ddpme[0].maxshift;
2996 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2998 if (dd->comm->ddpme[0].dim == YY)
3000 return dd->comm->ddpme[0].maxshift;
3002 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3004 return dd->comm->ddpme[1].maxshift;
3012 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3013 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3015 gmx_domdec_comm_t *comm;
3018 real range, pme_boundary;
3022 nc = dd->nc[ddpme->dim];
3025 if (!ddpme->dim_match)
3027 /* PP decomposition is not along dim: the worst situation */
3030 else if (ns <= 3 || (bUniform && ns == nc))
3032 /* The optimal situation */
3037 /* We need to check for all pme nodes which nodes they
3038 * could possibly need to communicate with.
3040 xmin = ddpme->pp_min;
3041 xmax = ddpme->pp_max;
3042 /* Allow for atoms to be maximally 2/3 times the cut-off
3043 * out of their DD cell. This is a reasonable balance between
3044 * between performance and support for most charge-group/cut-off
3047 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3048 /* Avoid extra communication when we are exactly at a boundary */
3052 for (s = 0; s < ns; s++)
3054 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3055 pme_boundary = (real)s/ns;
3058 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3060 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3064 pme_boundary = (real)(s+1)/ns;
3067 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3069 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3076 ddpme->maxshift = sh;
3080 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3081 ddpme->dim, ddpme->maxshift);
3085 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3089 for (d = 0; d < dd->ndim; d++)
3092 if (dim < ddbox->nboundeddim &&
3093 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3094 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3096 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",
3097 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3098 dd->nc[dim], dd->comm->cellsize_limit);
3104 setcellsizeslbLOCAL, setcellsizeslbMASTER, setcellsizeslbPULSE_ONLY
3107 /* Set the domain boundaries. Use for static (or no) load balancing,
3108 * and also for the starting state for dynamic load balancing.
3109 * setmode determine if and where the boundaries are stored, use enum above.
3110 * Returns the number communication pulses in npulse.
3112 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3113 int setmode, ivec npulse)
3115 gmx_domdec_comm_t *comm;
3118 real *cell_x, cell_dx, cellsize;
3122 for (d = 0; d < DIM; d++)
3124 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3126 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3129 cell_dx = ddbox->box_size[d]/dd->nc[d];
3132 case setcellsizeslbMASTER:
3133 for (j = 0; j < dd->nc[d]+1; j++)
3135 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3138 case setcellsizeslbLOCAL:
3139 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3140 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3145 cellsize = cell_dx*ddbox->skew_fac[d];
3146 while (cellsize*npulse[d] < comm->cutoff)
3150 cellsize_min[d] = cellsize;
3154 /* Statically load balanced grid */
3155 /* Also when we are not doing a master distribution we determine
3156 * all cell borders in a loop to obtain identical values
3157 * to the master distribution case and to determine npulse.
3159 if (setmode == setcellsizeslbMASTER)
3161 cell_x = dd->ma->cell_x[d];
3165 snew(cell_x, dd->nc[d]+1);
3167 cell_x[0] = ddbox->box0[d];
3168 for (j = 0; j < dd->nc[d]; j++)
3170 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3171 cell_x[j+1] = cell_x[j] + cell_dx;
3172 cellsize = cell_dx*ddbox->skew_fac[d];
3173 while (cellsize*npulse[d] < comm->cutoff &&
3174 npulse[d] < dd->nc[d]-1)
3178 cellsize_min[d] = std::min(cellsize_min[d], cellsize);
3180 if (setmode == setcellsizeslbLOCAL)
3182 comm->cell_x0[d] = cell_x[dd->ci[d]];
3183 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3185 if (setmode != setcellsizeslbMASTER)
3190 /* The following limitation is to avoid that a cell would receive
3191 * some of its own home charge groups back over the periodic boundary.
3192 * Double charge groups cause trouble with the global indices.
3194 if (d < ddbox->npbcdim &&
3195 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3197 char error_string[STRLEN];
3199 sprintf(error_string,
3200 "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",
3201 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3203 dd->nc[d], dd->nc[d],
3204 dd->nnodes > dd->nc[d] ? "cells" : "ranks");
3206 if (setmode == setcellsizeslbLOCAL)
3208 gmx_fatal_collective(FARGS, NULL, dd, error_string);
3212 gmx_fatal(FARGS, error_string);
3217 if (!comm->bDynLoadBal)
3219 copy_rvec(cellsize_min, comm->cellsize_min);
3222 for (d = 0; d < comm->npmedecompdim; d++)
3224 set_pme_maxshift(dd, &comm->ddpme[d],
3225 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3226 comm->ddpme[d].slb_dim_f);
3231 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3232 int d, int dim, gmx_domdec_root_t *root,
3234 gmx_bool bUniform, gmx_int64_t step, real cellsize_limit_f, int range[])
3236 gmx_domdec_comm_t *comm;
3237 int ncd, i, j, nmin, nmin_old;
3238 gmx_bool bLimLo, bLimHi;
3240 real fac, halfway, cellsize_limit_f_i, region_size;
3241 gmx_bool bPBC, bLastHi = FALSE;
3242 int nrange[] = {range[0], range[1]};
3244 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3250 bPBC = (dim < ddbox->npbcdim);
3252 cell_size = root->buf_ncd;
3256 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3259 /* First we need to check if the scaling does not make cells
3260 * smaller than the smallest allowed size.
3261 * We need to do this iteratively, since if a cell is too small,
3262 * it needs to be enlarged, which makes all the other cells smaller,
3263 * which could in turn make another cell smaller than allowed.
3265 for (i = range[0]; i < range[1]; i++)
3267 root->bCellMin[i] = FALSE;
3273 /* We need the total for normalization */
3275 for (i = range[0]; i < range[1]; i++)
3277 if (root->bCellMin[i] == FALSE)
3279 fac += cell_size[i];
3282 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3283 /* Determine the cell boundaries */
3284 for (i = range[0]; i < range[1]; i++)
3286 if (root->bCellMin[i] == FALSE)
3288 cell_size[i] *= fac;
3289 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3291 cellsize_limit_f_i = 0;
3295 cellsize_limit_f_i = cellsize_limit_f;
3297 if (cell_size[i] < cellsize_limit_f_i)
3299 root->bCellMin[i] = TRUE;
3300 cell_size[i] = cellsize_limit_f_i;
3304 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3307 while (nmin > nmin_old);
3310 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3311 /* For this check we should not use DD_CELL_MARGIN,
3312 * but a slightly smaller factor,
3313 * since rounding could get use below the limit.
3315 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3318 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",
3319 gmx_step_str(step, buf),
3320 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3321 ncd, comm->cellsize_min[dim]);
3324 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3328 /* Check if the boundary did not displace more than halfway
3329 * each of the cells it bounds, as this could cause problems,
3330 * especially when the differences between cell sizes are large.
3331 * If changes are applied, they will not make cells smaller
3332 * than the cut-off, as we check all the boundaries which
3333 * might be affected by a change and if the old state was ok,
3334 * the cells will at most be shrunk back to their old size.
3336 for (i = range[0]+1; i < range[1]; i++)
3338 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3339 if (root->cell_f[i] < halfway)
3341 root->cell_f[i] = halfway;
3342 /* Check if the change also causes shifts of the next boundaries */
3343 for (j = i+1; j < range[1]; j++)
3345 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3347 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3351 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3352 if (root->cell_f[i] > halfway)
3354 root->cell_f[i] = halfway;
3355 /* Check if the change also causes shifts of the next boundaries */
3356 for (j = i-1; j >= range[0]+1; j--)
3358 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3360 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3367 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3368 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3369 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3370 * for a and b nrange is used */
3373 /* Take care of the staggering of the cell boundaries */
3376 for (i = range[0]; i < range[1]; i++)
3378 root->cell_f_max0[i] = root->cell_f[i];
3379 root->cell_f_min1[i] = root->cell_f[i+1];
3384 for (i = range[0]+1; i < range[1]; i++)
3386 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3387 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3388 if (bLimLo && bLimHi)
3390 /* Both limits violated, try the best we can */
3391 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3392 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3393 nrange[0] = range[0];
3395 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3398 nrange[1] = range[1];
3399 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3405 /* root->cell_f[i] = root->bound_min[i]; */
3406 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3409 else if (bLimHi && !bLastHi)
3412 if (nrange[1] < range[1]) /* found a LimLo before */
3414 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3415 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3416 nrange[0] = nrange[1];
3418 root->cell_f[i] = root->bound_max[i];
3420 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3422 nrange[1] = range[1];
3425 if (nrange[1] < range[1]) /* found last a LimLo */
3427 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3428 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3429 nrange[0] = nrange[1];
3430 nrange[1] = range[1];
3431 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3433 else if (nrange[0] > range[0]) /* found at least one LimHi */
3435 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3442 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3443 int d, int dim, gmx_domdec_root_t *root,
3444 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3445 gmx_bool bUniform, gmx_int64_t step)
3447 gmx_domdec_comm_t *comm;
3448 int ncd, d1, i, pos;
3450 real load_aver, load_i, imbalance, change, change_max, sc;
3451 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3455 int range[] = { 0, 0 };
3459 /* Convert the maximum change from the input percentage to a fraction */
3460 change_limit = comm->dlb_scale_lim*0.01;
3464 bPBC = (dim < ddbox->npbcdim);
3466 cell_size = root->buf_ncd;
3468 /* Store the original boundaries */
3469 for (i = 0; i < ncd+1; i++)
3471 root->old_cell_f[i] = root->cell_f[i];
3475 for (i = 0; i < ncd; i++)
3477 cell_size[i] = 1.0/ncd;
3480 else if (dd_load_count(comm) > 0)
3482 load_aver = comm->load[d].sum_m/ncd;
3484 for (i = 0; i < ncd; i++)
3486 /* Determine the relative imbalance of cell i */
3487 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3488 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3489 /* Determine the change of the cell size using underrelaxation */
3490 change = -relax*imbalance;
3491 change_max = std::max(change_max, std::max(change, -change));
3493 /* Limit the amount of scaling.
3494 * We need to use the same rescaling for all cells in one row,
3495 * otherwise the load balancing might not converge.
3498 if (change_max > change_limit)
3500 sc *= change_limit/change_max;
3502 for (i = 0; i < ncd; i++)
3504 /* Determine the relative imbalance of cell i */
3505 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3506 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3507 /* Determine the change of the cell size using underrelaxation */
3508 change = -sc*imbalance;
3509 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3513 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3514 cellsize_limit_f *= DD_CELL_MARGIN;
3515 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3516 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3517 if (ddbox->tric_dir[dim])
3519 cellsize_limit_f /= ddbox->skew_fac[dim];
3520 dist_min_f /= ddbox->skew_fac[dim];
3522 if (bDynamicBox && d > 0)
3524 dist_min_f *= DD_PRES_SCALE_MARGIN;
3526 if (d > 0 && !bUniform)
3528 /* Make sure that the grid is not shifted too much */
3529 for (i = 1; i < ncd; i++)
3531 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3533 gmx_incons("Inconsistent DD boundary staggering limits!");
3535 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3536 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3539 root->bound_min[i] += 0.5*space;
3541 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3542 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3545 root->bound_max[i] += 0.5*space;
3550 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3552 root->cell_f_max0[i-1] + dist_min_f,
3553 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3554 root->cell_f_min1[i] - dist_min_f);
3559 root->cell_f[0] = 0;
3560 root->cell_f[ncd] = 1;
3561 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3564 /* After the checks above, the cells should obey the cut-off
3565 * restrictions, but it does not hurt to check.
3567 for (i = 0; i < ncd; i++)
3571 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3572 dim, i, root->cell_f[i], root->cell_f[i+1]);
3575 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3576 root->cell_f[i+1] - root->cell_f[i] <
3577 cellsize_limit_f/DD_CELL_MARGIN)
3581 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3582 gmx_step_str(step, buf), dim2char(dim), i,
3583 (root->cell_f[i+1] - root->cell_f[i])
3584 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3589 /* Store the cell boundaries of the lower dimensions at the end */
3590 for (d1 = 0; d1 < d; d1++)
3592 root->cell_f[pos++] = comm->cell_f0[d1];
3593 root->cell_f[pos++] = comm->cell_f1[d1];
3596 if (d < comm->npmedecompdim)
3598 /* The master determines the maximum shift for
3599 * the coordinate communication between separate PME nodes.
3601 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3603 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3606 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3610 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3611 gmx_ddbox_t *ddbox, int dimind)
3613 gmx_domdec_comm_t *comm;
3618 /* Set the cell dimensions */
3619 dim = dd->dim[dimind];
3620 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3621 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3622 if (dim >= ddbox->nboundeddim)
3624 comm->cell_x0[dim] += ddbox->box0[dim];
3625 comm->cell_x1[dim] += ddbox->box0[dim];
3629 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3630 int d, int dim, real *cell_f_row,
3633 gmx_domdec_comm_t *comm;
3639 /* Each node would only need to know two fractions,
3640 * but it is probably cheaper to broadcast the whole array.
3642 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3643 0, comm->mpi_comm_load[d]);
3645 /* Copy the fractions for this dimension from the buffer */
3646 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3647 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3648 /* The whole array was communicated, so set the buffer position */
3649 pos = dd->nc[dim] + 1;
3650 for (d1 = 0; d1 <= d; d1++)
3654 /* Copy the cell fractions of the lower dimensions */
3655 comm->cell_f0[d1] = cell_f_row[pos++];
3656 comm->cell_f1[d1] = cell_f_row[pos++];
3658 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3660 /* Convert the communicated shift from float to int */
3661 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3664 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3668 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3669 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3670 gmx_bool bUniform, gmx_int64_t step)
3672 gmx_domdec_comm_t *comm;
3674 gmx_bool bRowMember, bRowRoot;
3679 for (d = 0; d < dd->ndim; d++)
3684 for (d1 = d; d1 < dd->ndim; d1++)
3686 if (dd->ci[dd->dim[d1]] > 0)
3699 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3700 ddbox, bDynamicBox, bUniform, step);
3701 cell_f_row = comm->root[d]->cell_f;
3705 cell_f_row = comm->cell_f_row;
3707 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3712 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3716 /* This function assumes the box is static and should therefore
3717 * not be called when the box has changed since the last
3718 * call to dd_partition_system.
3720 for (d = 0; d < dd->ndim; d++)
3722 relative_to_absolute_cell_bounds(dd, ddbox, d);
3728 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3729 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3730 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3731 gmx_wallcycle_t wcycle)
3733 gmx_domdec_comm_t *comm;
3740 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3741 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3742 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3744 else if (bDynamicBox)
3746 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3749 /* Set the dimensions for which no DD is used */
3750 for (dim = 0; dim < DIM; dim++)
3752 if (dd->nc[dim] == 1)
3754 comm->cell_x0[dim] = 0;
3755 comm->cell_x1[dim] = ddbox->box_size[dim];
3756 if (dim >= ddbox->nboundeddim)
3758 comm->cell_x0[dim] += ddbox->box0[dim];
3759 comm->cell_x1[dim] += ddbox->box0[dim];
3765 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3768 gmx_domdec_comm_dim_t *cd;
3770 for (d = 0; d < dd->ndim; d++)
3772 cd = &dd->comm->cd[d];
3773 np = npulse[dd->dim[d]];
3774 if (np > cd->np_nalloc)
3778 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3779 dim2char(dd->dim[d]), np);
3781 if (DDMASTER(dd) && cd->np_nalloc > 0)
3783 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3785 srenew(cd->ind, np);
3786 for (i = cd->np_nalloc; i < np; i++)
3788 cd->ind[i].index = NULL;
3789 cd->ind[i].nalloc = 0;
3798 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3799 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3800 gmx_bool bUniform, gmx_bool bDoDLB, gmx_int64_t step,
3801 gmx_wallcycle_t wcycle)
3803 gmx_domdec_comm_t *comm;
3809 /* Copy the old cell boundaries for the cg displacement check */
3810 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3811 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3813 if (comm->bDynLoadBal)
3817 check_box_size(dd, ddbox);
3819 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3823 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbLOCAL, npulse);
3824 realloc_comm_ind(dd, npulse);
3829 for (d = 0; d < DIM; d++)
3831 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3832 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3837 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3839 rvec cell_ns_x0, rvec cell_ns_x1,
3842 gmx_domdec_comm_t *comm;
3847 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3849 dim = dd->dim[dim_ind];
3851 /* Without PBC we don't have restrictions on the outer cells */
3852 if (!(dim >= ddbox->npbcdim &&
3853 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3854 comm->bDynLoadBal &&
3855 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3856 comm->cellsize_min[dim])
3859 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",
3860 gmx_step_str(step, buf), dim2char(dim),
3861 comm->cell_x1[dim] - comm->cell_x0[dim],
3862 ddbox->skew_fac[dim],
3863 dd->comm->cellsize_min[dim],
3864 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3868 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3870 /* Communicate the boundaries and update cell_ns_x0/1 */
3871 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3872 if (dd->bGridJump && dd->ndim > 1)
3874 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3879 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3883 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3891 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3892 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3901 static void check_screw_box(matrix box)
3903 /* Mathematical limitation */
3904 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3906 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3909 /* Limitation due to the asymmetry of the eighth shell method */
3910 if (box[ZZ][YY] != 0)
3912 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3916 static void distribute_cg(FILE *fplog,
3917 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3920 gmx_domdec_master_t *ma;
3921 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3922 int i, icg, j, k, k0, k1, d;
3926 real nrcg, inv_ncg, pos_d;
3932 if (tmp_ind == NULL)
3934 snew(tmp_nalloc, dd->nnodes);
3935 snew(tmp_ind, dd->nnodes);
3936 for (i = 0; i < dd->nnodes; i++)
3938 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3939 snew(tmp_ind[i], tmp_nalloc[i]);
3943 /* Clear the count */
3944 for (i = 0; i < dd->nnodes; i++)
3950 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3952 cgindex = cgs->index;
3954 /* Compute the center of geometry for all charge groups */
3955 for (icg = 0; icg < cgs->nr; icg++)
3958 k1 = cgindex[icg+1];
3962 copy_rvec(pos[k0], cg_cm);
3969 for (k = k0; (k < k1); k++)
3971 rvec_inc(cg_cm, pos[k]);
3973 for (d = 0; (d < DIM); d++)
3975 cg_cm[d] *= inv_ncg;
3978 /* Put the charge group in the box and determine the cell index */
3979 for (d = DIM-1; d >= 0; d--)
3982 if (d < dd->npbcdim)
3984 bScrew = (dd->bScrewPBC && d == XX);
3985 if (tric_dir[d] && dd->nc[d] > 1)
3987 /* Use triclinic coordintates for this dimension */
3988 for (j = d+1; j < DIM; j++)
3990 pos_d += cg_cm[j]*tcm[j][d];
3993 while (pos_d >= box[d][d])
3996 rvec_dec(cg_cm, box[d]);
3999 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4000 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4002 for (k = k0; (k < k1); k++)
4004 rvec_dec(pos[k], box[d]);
4007 pos[k][YY] = box[YY][YY] - pos[k][YY];
4008 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4015 rvec_inc(cg_cm, box[d]);
4018 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4019 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4021 for (k = k0; (k < k1); k++)
4023 rvec_inc(pos[k], box[d]);
4026 pos[k][YY] = box[YY][YY] - pos[k][YY];
4027 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4032 /* This could be done more efficiently */
4034 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4039 i = dd_index(dd->nc, ind);
4040 if (ma->ncg[i] == tmp_nalloc[i])
4042 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4043 srenew(tmp_ind[i], tmp_nalloc[i]);
4045 tmp_ind[i][ma->ncg[i]] = icg;
4047 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4051 for (i = 0; i < dd->nnodes; i++)
4054 for (k = 0; k < ma->ncg[i]; k++)
4056 ma->cg[k1++] = tmp_ind[i][k];
4059 ma->index[dd->nnodes] = k1;
4061 for (i = 0; i < dd->nnodes; i++)
4070 /* Here we avoid int overflows due to #atoms^2: use double, dsqr */
4071 int nat_sum, nat_min, nat_max;
4076 nat_min = ma->nat[0];
4077 nat_max = ma->nat[0];
4078 for (i = 0; i < dd->nnodes; i++)
4080 nat_sum += ma->nat[i];
4081 nat2_sum += dsqr(ma->nat[i]);
4082 nat_min = std::min(nat_min, ma->nat[i]);
4083 nat_max = std::max(nat_max, ma->nat[i]);
4085 nat_sum /= dd->nnodes;
4086 nat2_sum /= dd->nnodes;
4088 fprintf(fplog, "Atom distribution over %d domains: av %d stddev %d min %d max %d\n",
4091 static_cast<int>(sqrt(nat2_sum - dsqr(nat_sum) + 0.5)),
4096 static void get_cg_distribution(FILE *fplog, gmx_domdec_t *dd,
4097 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4100 gmx_domdec_master_t *ma = NULL;
4103 int *ibuf, buf2[2] = { 0, 0 };
4104 gmx_bool bMaster = DDMASTER(dd);
4112 check_screw_box(box);
4115 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbMASTER, npulse);
4117 distribute_cg(fplog, box, ddbox->tric_dir, cgs, pos, dd);
4118 for (i = 0; i < dd->nnodes; i++)
4120 ma->ibuf[2*i] = ma->ncg[i];
4121 ma->ibuf[2*i+1] = ma->nat[i];
4129 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4131 dd->ncg_home = buf2[0];
4132 dd->nat_home = buf2[1];
4133 dd->ncg_tot = dd->ncg_home;
4134 dd->nat_tot = dd->nat_home;
4135 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4137 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4138 srenew(dd->index_gl, dd->cg_nalloc);
4139 srenew(dd->cgindex, dd->cg_nalloc+1);
4143 for (i = 0; i < dd->nnodes; i++)
4145 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4146 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4151 bMaster ? ma->ibuf : NULL,
4152 bMaster ? ma->ibuf+dd->nnodes : NULL,
4153 bMaster ? ma->cg : NULL,
4154 dd->ncg_home*sizeof(int), dd->index_gl);
4156 /* Determine the home charge group sizes */
4158 for (i = 0; i < dd->ncg_home; i++)
4160 cg_gl = dd->index_gl[i];
4162 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4167 fprintf(debug, "Home charge groups:\n");
4168 for (i = 0; i < dd->ncg_home; i++)
4170 fprintf(debug, " %d", dd->index_gl[i]);
4173 fprintf(debug, "\n");
4176 fprintf(debug, "\n");
4180 static int compact_and_copy_vec_at(int ncg, int *move,
4183 rvec *src, gmx_domdec_comm_t *comm,
4186 int m, icg, i, i0, i1, nrcg;
4192 for (m = 0; m < DIM*2; m++)
4198 for (icg = 0; icg < ncg; icg++)
4200 i1 = cgindex[icg+1];
4206 /* Compact the home array in place */
4207 for (i = i0; i < i1; i++)
4209 copy_rvec(src[i], src[home_pos++]);
4215 /* Copy to the communication buffer */
4217 pos_vec[m] += 1 + vec*nrcg;
4218 for (i = i0; i < i1; i++)
4220 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4222 pos_vec[m] += (nvec - vec - 1)*nrcg;
4226 home_pos += i1 - i0;
4234 static int compact_and_copy_vec_cg(int ncg, int *move,
4236 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4239 int m, icg, i0, i1, nrcg;
4245 for (m = 0; m < DIM*2; m++)
4251 for (icg = 0; icg < ncg; icg++)
4253 i1 = cgindex[icg+1];
4259 /* Compact the home array in place */
4260 copy_rvec(src[icg], src[home_pos++]);
4266 /* Copy to the communication buffer */
4267 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4268 pos_vec[m] += 1 + nrcg*nvec;
4280 static int compact_ind(int ncg, int *move,
4281 int *index_gl, int *cgindex,
4283 gmx_ga2la_t ga2la, char *bLocalCG,
4286 int cg, nat, a0, a1, a, a_gl;
4291 for (cg = 0; cg < ncg; cg++)
4297 /* Compact the home arrays in place.
4298 * Anything that can be done here avoids access to global arrays.
4300 cgindex[home_pos] = nat;
4301 for (a = a0; a < a1; a++)
4304 gatindex[nat] = a_gl;
4305 /* The cell number stays 0, so we don't need to set it */
4306 ga2la_change_la(ga2la, a_gl, nat);
4309 index_gl[home_pos] = index_gl[cg];
4310 cginfo[home_pos] = cginfo[cg];
4311 /* The charge group remains local, so bLocalCG does not change */
4316 /* Clear the global indices */
4317 for (a = a0; a < a1; a++)
4319 ga2la_del(ga2la, gatindex[a]);
4323 bLocalCG[index_gl[cg]] = FALSE;
4327 cgindex[home_pos] = nat;
4332 static void clear_and_mark_ind(int ncg, int *move,
4333 int *index_gl, int *cgindex, int *gatindex,
4334 gmx_ga2la_t ga2la, char *bLocalCG,
4339 for (cg = 0; cg < ncg; cg++)
4345 /* Clear the global indices */
4346 for (a = a0; a < a1; a++)
4348 ga2la_del(ga2la, gatindex[a]);
4352 bLocalCG[index_gl[cg]] = FALSE;
4354 /* Signal that this cg has moved using the ns cell index.
4355 * Here we set it to -1. fill_grid will change it
4356 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4358 cell_index[cg] = -1;
4363 static void print_cg_move(FILE *fplog,
4365 gmx_int64_t step, int cg, int dim, int dir,
4366 gmx_bool bHaveCgcmOld, real limitd,
4367 rvec cm_old, rvec cm_new, real pos_d)
4369 gmx_domdec_comm_t *comm;
4374 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4377 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4378 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4379 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4383 /* We don't have a limiting distance available: don't print it */
4384 fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
4385 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
4386 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4388 fprintf(fplog, "distance out of cell %f\n",
4389 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4392 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4393 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4395 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4396 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4397 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4399 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4400 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4402 comm->cell_x0[dim], comm->cell_x1[dim]);
4405 static void cg_move_error(FILE *fplog,
4407 gmx_int64_t step, int cg, int dim, int dir,
4408 gmx_bool bHaveCgcmOld, real limitd,
4409 rvec cm_old, rvec cm_new, real pos_d)
4413 print_cg_move(fplog, dd, step, cg, dim, dir,
4414 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4416 print_cg_move(stderr, dd, step, cg, dim, dir,
4417 bHaveCgcmOld, limitd, cm_old, cm_new, pos_d);
4419 "%s moved too far between two domain decomposition steps\n"
4420 "This usually means that your system is not well equilibrated",
4421 dd->comm->bCGs ? "A charge group" : "An atom");
4424 static void rotate_state_atom(t_state *state, int a)
4428 for (est = 0; est < estNR; est++)
4430 if (EST_DISTR(est) && (state->flags & (1<<est)))
4435 /* Rotate the complete state; for a rectangular box only */
4436 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4437 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4440 state->v[a][YY] = -state->v[a][YY];
4441 state->v[a][ZZ] = -state->v[a][ZZ];
4444 state->sd_X[a][YY] = -state->sd_X[a][YY];
4445 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4448 state->cg_p[a][YY] = -state->cg_p[a][YY];
4449 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4451 case estDISRE_INITF:
4452 case estDISRE_RM3TAV:
4453 case estORIRE_INITF:
4455 /* These are distances, so not affected by rotation */
4458 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4464 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4466 if (natoms > comm->moved_nalloc)
4468 /* Contents should be preserved here */
4469 comm->moved_nalloc = over_alloc_dd(natoms);
4470 srenew(comm->moved, comm->moved_nalloc);
4476 static void calc_cg_move(FILE *fplog, gmx_int64_t step,
4479 ivec tric_dir, matrix tcm,
4480 rvec cell_x0, rvec cell_x1,
4481 rvec limitd, rvec limit0, rvec limit1,
4483 int cg_start, int cg_end,
4488 int cg, k, k0, k1, d, dim, d2;
4493 real inv_ncg, pos_d;
4496 npbcdim = dd->npbcdim;
4498 for (cg = cg_start; cg < cg_end; cg++)
4505 copy_rvec(state->x[k0], cm_new);
4512 for (k = k0; (k < k1); k++)
4514 rvec_inc(cm_new, state->x[k]);
4516 for (d = 0; (d < DIM); d++)
4518 cm_new[d] = inv_ncg*cm_new[d];
4523 /* Do pbc and check DD cell boundary crossings */
4524 for (d = DIM-1; d >= 0; d--)
4528 bScrew = (dd->bScrewPBC && d == XX);
4529 /* Determine the location of this cg in lattice coordinates */
4533 for (d2 = d+1; d2 < DIM; d2++)
4535 pos_d += cm_new[d2]*tcm[d2][d];
4538 /* Put the charge group in the triclinic unit-cell */
4539 if (pos_d >= cell_x1[d])
4541 if (pos_d >= limit1[d])
4543 cg_move_error(fplog, dd, step, cg, d, 1,
4544 cg_cm != state->x, limitd[d],
4545 cg_cm[cg], cm_new, pos_d);
4548 if (dd->ci[d] == dd->nc[d] - 1)
4550 rvec_dec(cm_new, state->box[d]);
4553 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4554 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4556 for (k = k0; (k < k1); k++)
4558 rvec_dec(state->x[k], state->box[d]);
4561 rotate_state_atom(state, k);
4566 else if (pos_d < cell_x0[d])
4568 if (pos_d < limit0[d])
4570 cg_move_error(fplog, dd, step, cg, d, -1,
4571 cg_cm != state->x, limitd[d],
4572 cg_cm[cg], cm_new, pos_d);
4577 rvec_inc(cm_new, state->box[d]);
4580 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4581 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4583 for (k = k0; (k < k1); k++)
4585 rvec_inc(state->x[k], state->box[d]);
4588 rotate_state_atom(state, k);
4594 else if (d < npbcdim)
4596 /* Put the charge group in the rectangular unit-cell */
4597 while (cm_new[d] >= state->box[d][d])
4599 rvec_dec(cm_new, state->box[d]);
4600 for (k = k0; (k < k1); k++)
4602 rvec_dec(state->x[k], state->box[d]);
4605 while (cm_new[d] < 0)
4607 rvec_inc(cm_new, state->box[d]);
4608 for (k = k0; (k < k1); k++)
4610 rvec_inc(state->x[k], state->box[d]);
4616 copy_rvec(cm_new, cg_cm[cg]);
4618 /* Determine where this cg should go */
4621 for (d = 0; d < dd->ndim; d++)
4626 flag |= DD_FLAG_FW(d);
4632 else if (dev[dim] == -1)
4634 flag |= DD_FLAG_BW(d);
4637 if (dd->nc[dim] > 2)
4648 /* Temporarily store the flag in move */
4649 move[cg] = mc + flag;
4653 static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
4654 gmx_domdec_t *dd, ivec tric_dir,
4655 t_state *state, rvec **f,
4664 int ncg[DIM*2], nat[DIM*2];
4665 int c, i, cg, k, d, dim, dim2, dir, d2, d3;
4666 int mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
4667 int sbuf[2], rbuf[2];
4668 int home_pos_cg, home_pos_at, buf_pos;
4670 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4673 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
4675 cginfo_mb_t *cginfo_mb;
4676 gmx_domdec_comm_t *comm;
4678 int nthread, thread;
4682 check_screw_box(state->box);
4686 if (fr->cutoff_scheme == ecutsGROUP)
4691 for (i = 0; i < estNR; i++)
4697 case estX: /* Always present */ break;
4698 case estV: bV = (state->flags & (1<<i)); break;
4699 case estSDX: bSDX = (state->flags & (1<<i)); break;
4700 case estCGP: bCGP = (state->flags & (1<<i)); break;
4703 case estDISRE_INITF:
4704 case estDISRE_RM3TAV:
4705 case estORIRE_INITF:
4707 /* No processing required */
4710 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4715 if (dd->ncg_tot > comm->nalloc_int)
4717 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4718 srenew(comm->buf_int, comm->nalloc_int);
4720 move = comm->buf_int;
4722 /* Clear the count */
4723 for (c = 0; c < dd->ndim*2; c++)
4729 npbcdim = dd->npbcdim;
4731 for (d = 0; (d < DIM); d++)
4733 limitd[d] = dd->comm->cellsize_min[d];
4734 if (d >= npbcdim && dd->ci[d] == 0)
4736 cell_x0[d] = -GMX_FLOAT_MAX;
4740 cell_x0[d] = comm->cell_x0[d];
4742 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4744 cell_x1[d] = GMX_FLOAT_MAX;
4748 cell_x1[d] = comm->cell_x1[d];
4752 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4753 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4757 /* We check after communication if a charge group moved
4758 * more than one cell. Set the pre-comm check limit to float_max.
4760 limit0[d] = -GMX_FLOAT_MAX;
4761 limit1[d] = GMX_FLOAT_MAX;
4765 make_tric_corr_matrix(npbcdim, state->box, tcm);
4767 cgindex = dd->cgindex;
4769 nthread = gmx_omp_nthreads_get(emntDomdec);
4771 /* Compute the center of geometry for all home charge groups
4772 * and put them in the box and determine where they should go.
4774 #pragma omp parallel for num_threads(nthread) schedule(static)
4775 for (thread = 0; thread < nthread; thread++)
4777 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4778 cell_x0, cell_x1, limitd, limit0, limit1,
4780 ( thread *dd->ncg_home)/nthread,
4781 ((thread+1)*dd->ncg_home)/nthread,
4782 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4786 for (cg = 0; cg < dd->ncg_home; cg++)
4791 flag = mc & ~DD_FLAG_NRCG;
4792 mc = mc & DD_FLAG_NRCG;
4795 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4797 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4798 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4800 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4801 /* We store the cg size in the lower 16 bits
4802 * and the place where the charge group should go
4803 * in the next 6 bits. This saves some communication volume.
4805 nrcg = cgindex[cg+1] - cgindex[cg];
4806 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4812 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4813 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4816 for (i = 0; i < dd->ndim*2; i++)
4818 *ncg_moved += ncg[i];
4835 /* Make sure the communication buffers are large enough */
4836 for (mc = 0; mc < dd->ndim*2; mc++)
4838 nvr = ncg[mc] + nat[mc]*nvec;
4839 if (nvr > comm->cgcm_state_nalloc[mc])
4841 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4842 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4846 switch (fr->cutoff_scheme)
4849 /* Recalculating cg_cm might be cheaper than communicating,
4850 * but that could give rise to rounding issues.
4853 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4854 nvec, cg_cm, comm, bCompact);
4857 /* Without charge groups we send the moved atom coordinates
4858 * over twice. This is so the code below can be used without
4859 * many conditionals for both for with and without charge groups.
4862 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4863 nvec, state->x, comm, FALSE);
4866 home_pos_cg -= *ncg_moved;
4870 gmx_incons("unimplemented");
4876 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4877 nvec, vec++, state->x, comm, bCompact);
4880 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4881 nvec, vec++, state->v, comm, bCompact);
4885 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4886 nvec, vec++, state->sd_X, comm, bCompact);
4890 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4891 nvec, vec++, state->cg_p, comm, bCompact);
4896 compact_ind(dd->ncg_home, move,
4897 dd->index_gl, dd->cgindex, dd->gatindex,
4898 dd->ga2la, comm->bLocalCG,
4903 if (fr->cutoff_scheme == ecutsVERLET)
4905 moved = get_moved(comm, dd->ncg_home);
4907 for (k = 0; k < dd->ncg_home; k++)
4914 moved = fr->ns.grid->cell_index;
4917 clear_and_mark_ind(dd->ncg_home, move,
4918 dd->index_gl, dd->cgindex, dd->gatindex,
4919 dd->ga2la, comm->bLocalCG,
4923 cginfo_mb = fr->cginfo_mb;
4925 *ncg_stay_home = home_pos_cg;
4926 for (d = 0; d < dd->ndim; d++)
4931 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4934 /* Communicate the cg and atom counts */
4939 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4940 d, dir, sbuf[0], sbuf[1]);
4942 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4944 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4946 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4947 srenew(comm->buf_int, comm->nalloc_int);
4950 /* Communicate the charge group indices, sizes and flags */
4951 dd_sendrecv_int(dd, d, dir,
4952 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4953 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4955 nvs = ncg[cdd] + nat[cdd]*nvec;
4956 i = rbuf[0] + rbuf[1] *nvec;
4957 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4959 /* Communicate cgcm and state */
4960 dd_sendrecv_rvec(dd, d, dir,
4961 comm->cgcm_state[cdd], nvs,
4962 comm->vbuf.v+nvr, i);
4963 ncg_recv += rbuf[0];
4967 /* Process the received charge groups */
4969 for (cg = 0; cg < ncg_recv; cg++)
4971 flag = comm->buf_int[cg*DD_CGIBS+1];
4973 if (dim >= npbcdim && dd->nc[dim] > 2)
4975 /* No pbc in this dim and more than one domain boundary.
4976 * We do a separate check if a charge group didn't move too far.
4978 if (((flag & DD_FLAG_FW(d)) &&
4979 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4980 ((flag & DD_FLAG_BW(d)) &&
4981 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4983 cg_move_error(fplog, dd, step, cg, dim,
4984 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4985 fr->cutoff_scheme == ecutsGROUP, 0,
4986 comm->vbuf.v[buf_pos],
4987 comm->vbuf.v[buf_pos],
4988 comm->vbuf.v[buf_pos][dim]);
4995 /* Check which direction this cg should go */
4996 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
5000 /* The cell boundaries for dimension d2 are not equal
5001 * for each cell row of the lower dimension(s),
5002 * therefore we might need to redetermine where
5003 * this cg should go.
5006 /* If this cg crosses the box boundary in dimension d2
5007 * we can use the communicated flag, so we do not
5008 * have to worry about pbc.
5010 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
5011 (flag & DD_FLAG_FW(d2))) ||
5012 (dd->ci[dim2] == 0 &&
5013 (flag & DD_FLAG_BW(d2)))))
5015 /* Clear the two flags for this dimension */
5016 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
5017 /* Determine the location of this cg
5018 * in lattice coordinates
5020 pos_d = comm->vbuf.v[buf_pos][dim2];
5023 for (d3 = dim2+1; d3 < DIM; d3++)
5026 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
5029 /* Check of we are not at the box edge.
5030 * pbc is only handled in the first step above,
5031 * but this check could move over pbc while
5032 * the first step did not due to different rounding.
5034 if (pos_d >= cell_x1[dim2] &&
5035 dd->ci[dim2] != dd->nc[dim2]-1)
5037 flag |= DD_FLAG_FW(d2);
5039 else if (pos_d < cell_x0[dim2] &&
5042 flag |= DD_FLAG_BW(d2);
5044 comm->buf_int[cg*DD_CGIBS+1] = flag;
5047 /* Set to which neighboring cell this cg should go */
5048 if (flag & DD_FLAG_FW(d2))
5052 else if (flag & DD_FLAG_BW(d2))
5054 if (dd->nc[dd->dim[d2]] > 2)
5066 nrcg = flag & DD_FLAG_NRCG;
5069 if (home_pos_cg+1 > dd->cg_nalloc)
5071 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5072 srenew(dd->index_gl, dd->cg_nalloc);
5073 srenew(dd->cgindex, dd->cg_nalloc+1);
5075 /* Set the global charge group index and size */
5076 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5077 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5078 /* Copy the state from the buffer */
5079 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5080 if (fr->cutoff_scheme == ecutsGROUP)
5083 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5087 /* Set the cginfo */
5088 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5089 dd->index_gl[home_pos_cg]);
5092 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5095 if (home_pos_at+nrcg > state->nalloc)
5097 dd_realloc_state(state, f, home_pos_at+nrcg);
5099 for (i = 0; i < nrcg; i++)
5101 copy_rvec(comm->vbuf.v[buf_pos++],
5102 state->x[home_pos_at+i]);
5106 for (i = 0; i < nrcg; i++)
5108 copy_rvec(comm->vbuf.v[buf_pos++],
5109 state->v[home_pos_at+i]);
5114 for (i = 0; i < nrcg; i++)
5116 copy_rvec(comm->vbuf.v[buf_pos++],
5117 state->sd_X[home_pos_at+i]);
5122 for (i = 0; i < nrcg; i++)
5124 copy_rvec(comm->vbuf.v[buf_pos++],
5125 state->cg_p[home_pos_at+i]);
5129 home_pos_at += nrcg;
5133 /* Reallocate the buffers if necessary */
5134 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5136 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5137 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5139 nvr = ncg[mc] + nat[mc]*nvec;
5140 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5142 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5143 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5145 /* Copy from the receive to the send buffers */
5146 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5147 comm->buf_int + cg*DD_CGIBS,
5148 DD_CGIBS*sizeof(int));
5149 memcpy(comm->cgcm_state[mc][nvr],
5150 comm->vbuf.v[buf_pos],
5151 (1+nrcg*nvec)*sizeof(rvec));
5152 buf_pos += 1 + nrcg*nvec;
5159 /* With sorting (!bCompact) the indices are now only partially up to date
5160 * and ncg_home and nat_home are not the real count, since there are
5161 * "holes" in the arrays for the charge groups that moved to neighbors.
5163 if (fr->cutoff_scheme == ecutsVERLET)
5165 moved = get_moved(comm, home_pos_cg);
5167 for (i = dd->ncg_home; i < home_pos_cg; i++)
5172 dd->ncg_home = home_pos_cg;
5173 dd->nat_home = home_pos_at;
5178 "Finished repartitioning: cgs moved out %d, new home %d\n",
5179 *ncg_moved, dd->ncg_home-*ncg_moved);
5184 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5186 dd->comm->cycl[ddCycl] += cycles;
5187 dd->comm->cycl_n[ddCycl]++;
5188 if (cycles > dd->comm->cycl_max[ddCycl])
5190 dd->comm->cycl_max[ddCycl] = cycles;
5194 static double force_flop_count(t_nrnb *nrnb)
5201 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5203 /* To get closer to the real timings, we half the count
5204 * for the normal loops and again half it for water loops.
5207 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5209 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5213 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5216 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5219 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5221 sum += nrnb->n[i]*cost_nrnb(i);
5224 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5226 sum += nrnb->n[i]*cost_nrnb(i);
5232 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5234 if (dd->comm->eFlop)
5236 dd->comm->flop -= force_flop_count(nrnb);
5239 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5241 if (dd->comm->eFlop)
5243 dd->comm->flop += force_flop_count(nrnb);
5248 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5252 for (i = 0; i < ddCyclNr; i++)
5254 dd->comm->cycl[i] = 0;
5255 dd->comm->cycl_n[i] = 0;
5256 dd->comm->cycl_max[i] = 0;
5259 dd->comm->flop_n = 0;
5262 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5264 gmx_domdec_comm_t *comm;
5265 gmx_domdec_load_t *load;
5266 gmx_domdec_root_t *root = NULL;
5268 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5273 fprintf(debug, "get_load_distribution start\n");
5276 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5280 bSepPME = (dd->pme_nodeid >= 0);
5282 for (d = dd->ndim-1; d >= 0; d--)
5285 /* Check if we participate in the communication in this dimension */
5286 if (d == dd->ndim-1 ||
5287 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5289 load = &comm->load[d];
5292 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5295 if (d == dd->ndim-1)
5297 sbuf[pos++] = dd_force_load(comm);
5298 sbuf[pos++] = sbuf[0];
5301 sbuf[pos++] = sbuf[0];
5302 sbuf[pos++] = cell_frac;
5305 sbuf[pos++] = comm->cell_f_max0[d];
5306 sbuf[pos++] = comm->cell_f_min1[d];
5311 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5312 sbuf[pos++] = comm->cycl[ddCyclPME];
5317 sbuf[pos++] = comm->load[d+1].sum;
5318 sbuf[pos++] = comm->load[d+1].max;
5321 sbuf[pos++] = comm->load[d+1].sum_m;
5322 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5323 sbuf[pos++] = comm->load[d+1].flags;
5326 sbuf[pos++] = comm->cell_f_max0[d];
5327 sbuf[pos++] = comm->cell_f_min1[d];
5332 sbuf[pos++] = comm->load[d+1].mdf;
5333 sbuf[pos++] = comm->load[d+1].pme;
5337 /* Communicate a row in DD direction d.
5338 * The communicators are setup such that the root always has rank 0.
5341 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5342 load->load, load->nload*sizeof(float), MPI_BYTE,
5343 0, comm->mpi_comm_load[d]);
5345 if (dd->ci[dim] == dd->master_ci[dim])
5347 /* We are the root, process this row */
5348 if (comm->bDynLoadBal)
5350 root = comm->root[d];
5360 for (i = 0; i < dd->nc[dim]; i++)
5362 load->sum += load->load[pos++];
5363 load->max = std::max(load->max, load->load[pos]);
5369 /* This direction could not be load balanced properly,
5370 * therefore we need to use the maximum iso the average load.
5372 load->sum_m = std::max(load->sum_m, load->load[pos]);
5376 load->sum_m += load->load[pos];
5379 load->cvol_min = std::min(load->cvol_min, load->load[pos]);
5383 load->flags = (int)(load->load[pos++] + 0.5);
5387 root->cell_f_max0[i] = load->load[pos++];
5388 root->cell_f_min1[i] = load->load[pos++];
5393 load->mdf = std::max(load->mdf, load->load[pos]);
5395 load->pme = std::max(load->pme, load->load[pos]);
5399 if (comm->bDynLoadBal && root->bLimited)
5401 load->sum_m *= dd->nc[dim];
5402 load->flags |= (1<<d);
5410 comm->nload += dd_load_count(comm);
5411 comm->load_step += comm->cycl[ddCyclStep];
5412 comm->load_sum += comm->load[0].sum;
5413 comm->load_max += comm->load[0].max;
5414 if (comm->bDynLoadBal)
5416 for (d = 0; d < dd->ndim; d++)
5418 if (comm->load[0].flags & (1<<d))
5420 comm->load_lim[d]++;
5426 comm->load_mdf += comm->load[0].mdf;
5427 comm->load_pme += comm->load[0].pme;
5431 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5435 fprintf(debug, "get_load_distribution finished\n");
5439 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5441 /* Return the relative performance loss on the total run time
5442 * due to the force calculation load imbalance.
5444 if (dd->comm->nload > 0)
5447 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5448 (dd->comm->load_step*dd->nnodes);
5456 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5459 int npp, npme, nnodes, d, limp;
5460 float imbal, pme_f_ratio, lossf, lossp = 0;
5462 gmx_domdec_comm_t *comm;
5465 if (DDMASTER(dd) && comm->nload > 0)
5468 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5469 nnodes = npp + npme;
5470 imbal = comm->load_max*npp/comm->load_sum - 1;
5471 lossf = dd_force_imb_perf_loss(dd);
5472 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5473 fprintf(fplog, "%s", buf);
5474 fprintf(stderr, "\n");
5475 fprintf(stderr, "%s", buf);
5476 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5477 fprintf(fplog, "%s", buf);
5478 fprintf(stderr, "%s", buf);
5480 if (comm->bDynLoadBal)
5482 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5483 for (d = 0; d < dd->ndim; d++)
5485 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5486 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5492 sprintf(buf+strlen(buf), "\n");
5493 fprintf(fplog, "%s", buf);
5494 fprintf(stderr, "%s", buf);
5498 pme_f_ratio = comm->load_pme/comm->load_mdf;
5499 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5502 lossp *= (float)npme/(float)nnodes;
5506 lossp *= (float)npp/(float)nnodes;
5508 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5509 fprintf(fplog, "%s", buf);
5510 fprintf(stderr, "%s", buf);
5511 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5512 fprintf(fplog, "%s", buf);
5513 fprintf(stderr, "%s", buf);
5515 fprintf(fplog, "\n");
5516 fprintf(stderr, "\n");
5518 if (lossf >= DD_PERF_LOSS_WARN)
5521 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5522 " in the domain decomposition.\n", lossf*100);
5523 if (!comm->bDynLoadBal)
5525 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5529 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5531 fprintf(fplog, "%s\n", buf);
5532 fprintf(stderr, "%s\n", buf);
5534 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
5537 "NOTE: %.1f %% performance was lost because the PME ranks\n"
5538 " had %s work to do than the PP ranks.\n"
5539 " You might want to %s the number of PME ranks\n"
5540 " or %s the cut-off and the grid spacing.\n",
5542 (lossp < 0) ? "less" : "more",
5543 (lossp < 0) ? "decrease" : "increase",
5544 (lossp < 0) ? "decrease" : "increase");
5545 fprintf(fplog, "%s\n", buf);
5546 fprintf(stderr, "%s\n", buf);
5551 static float dd_vol_min(gmx_domdec_t *dd)
5553 return dd->comm->load[0].cvol_min*dd->nnodes;
5556 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5558 return dd->comm->load[0].flags;
5561 static float dd_f_imbal(gmx_domdec_t *dd)
5563 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5566 float dd_pme_f_ratio(gmx_domdec_t *dd)
5568 if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
5570 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5578 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_int64_t step)
5583 flags = dd_load_flags(dd);
5587 "DD load balancing is limited by minimum cell size in dimension");
5588 for (d = 0; d < dd->ndim; d++)
5592 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5595 fprintf(fplog, "\n");
5597 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5598 if (dd->comm->bDynLoadBal)
5600 fprintf(fplog, " vol min/aver %5.3f%c",
5601 dd_vol_min(dd), flags ? '!' : ' ');
5603 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5604 if (dd->comm->cycl_n[ddCyclPME])
5606 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5608 fprintf(fplog, "\n\n");
5611 static void dd_print_load_verbose(gmx_domdec_t *dd)
5613 if (dd->comm->bDynLoadBal)
5615 fprintf(stderr, "vol %4.2f%c ",
5616 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5618 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5619 if (dd->comm->cycl_n[ddCyclPME])
5621 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5626 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5631 gmx_domdec_root_t *root;
5632 gmx_bool bPartOfGroup = FALSE;
5634 dim = dd->dim[dim_ind];
5635 copy_ivec(loc, loc_c);
5636 for (i = 0; i < dd->nc[dim]; i++)
5639 rank = dd_index(dd->nc, loc_c);
5640 if (rank == dd->rank)
5642 /* This process is part of the group */
5643 bPartOfGroup = TRUE;
5646 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5650 dd->comm->mpi_comm_load[dim_ind] = c_row;
5651 if (dd->comm->eDLB != edlbNO)
5653 if (dd->ci[dim] == dd->master_ci[dim])
5655 /* This is the root process of this row */
5656 snew(dd->comm->root[dim_ind], 1);
5657 root = dd->comm->root[dim_ind];
5658 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5659 snew(root->old_cell_f, dd->nc[dim]+1);
5660 snew(root->bCellMin, dd->nc[dim]);
5663 snew(root->cell_f_max0, dd->nc[dim]);
5664 snew(root->cell_f_min1, dd->nc[dim]);
5665 snew(root->bound_min, dd->nc[dim]);
5666 snew(root->bound_max, dd->nc[dim]);
5668 snew(root->buf_ncd, dd->nc[dim]);
5672 /* This is not a root process, we only need to receive cell_f */
5673 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5676 if (dd->ci[dim] == dd->master_ci[dim])
5678 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5684 void dd_setup_dlb_resource_sharing(t_commrec gmx_unused *cr,
5685 const gmx_hw_info_t gmx_unused *hwinfo,
5686 const gmx_hw_opt_t gmx_unused *hw_opt)
5689 int physicalnode_id_hash;
5692 MPI_Comm mpi_comm_pp_physicalnode;
5694 if (!(cr->duty & DUTY_PP) || hw_opt->gpu_opt.n_dev_use == 0)
5696 /* Only PP nodes (currently) use GPUs.
5697 * If we don't have GPUs, there are no resources to share.
5702 physicalnode_id_hash = gmx_physicalnode_id_hash();
5704 gpu_id = get_cuda_gpu_device_id(&hwinfo->gpu_info, &hw_opt->gpu_opt, cr->rank_pp_intranode);
5710 fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
5711 fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
5712 dd->rank, physicalnode_id_hash, gpu_id);
5714 /* Split the PP communicator over the physical nodes */
5715 /* TODO: See if we should store this (before), as it's also used for
5716 * for the nodecomm summution.
5718 MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
5719 &mpi_comm_pp_physicalnode);
5720 MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
5721 &dd->comm->mpi_comm_gpu_shared);
5722 MPI_Comm_free(&mpi_comm_pp_physicalnode);
5723 MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
5727 fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
5730 /* Note that some ranks could share a GPU, while others don't */
5732 if (dd->comm->nrank_gpu_shared == 1)
5734 MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
5739 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
5742 int dim0, dim1, i, j;
5747 fprintf(debug, "Making load communicators\n");
5750 snew(dd->comm->load, dd->ndim);
5751 snew(dd->comm->mpi_comm_load, dd->ndim);
5754 make_load_communicator(dd, 0, loc);
5758 for (i = 0; i < dd->nc[dim0]; i++)
5761 make_load_communicator(dd, 1, loc);
5767 for (i = 0; i < dd->nc[dim0]; i++)
5771 for (j = 0; j < dd->nc[dim1]; j++)
5774 make_load_communicator(dd, 2, loc);
5781 fprintf(debug, "Finished making load communicators\n");
5786 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5788 int d, dim, i, j, m;
5791 ivec dd_zp[DD_MAXIZONE];
5792 gmx_domdec_zones_t *zones;
5793 gmx_domdec_ns_ranges_t *izone;
5795 for (d = 0; d < dd->ndim; d++)
5798 copy_ivec(dd->ci, tmp);
5799 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5800 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5801 copy_ivec(dd->ci, tmp);
5802 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5803 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5806 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5809 dd->neighbor[d][1]);
5815 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5817 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5818 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5825 for (i = 0; i < nzonep; i++)
5827 copy_ivec(dd_zp3[i], dd_zp[i]);
5833 for (i = 0; i < nzonep; i++)
5835 copy_ivec(dd_zp2[i], dd_zp[i]);
5841 for (i = 0; i < nzonep; i++)
5843 copy_ivec(dd_zp1[i], dd_zp[i]);
5847 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5852 zones = &dd->comm->zones;
5854 for (i = 0; i < nzone; i++)
5857 clear_ivec(zones->shift[i]);
5858 for (d = 0; d < dd->ndim; d++)
5860 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5865 for (i = 0; i < nzone; i++)
5867 for (d = 0; d < DIM; d++)
5869 s[d] = dd->ci[d] - zones->shift[i][d];
5874 else if (s[d] >= dd->nc[d])
5880 zones->nizone = nzonep;
5881 for (i = 0; i < zones->nizone; i++)
5883 if (dd_zp[i][0] != i)
5885 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5887 izone = &zones->izone[i];
5888 izone->j0 = dd_zp[i][1];
5889 izone->j1 = dd_zp[i][2];
5890 for (dim = 0; dim < DIM; dim++)
5892 if (dd->nc[dim] == 1)
5894 /* All shifts should be allowed */
5895 izone->shift0[dim] = -1;
5896 izone->shift1[dim] = 1;
5901 izone->shift0[d] = 0;
5902 izone->shift1[d] = 0;
5903 for(j=izone->j0; j<izone->j1; j++) {
5904 if (dd->shift[j][d] > dd->shift[i][d])
5905 izone->shift0[d] = -1;
5906 if (dd->shift[j][d] < dd->shift[i][d])
5907 izone->shift1[d] = 1;
5913 /* Assume the shift are not more than 1 cell */
5914 izone->shift0[dim] = 1;
5915 izone->shift1[dim] = -1;
5916 for (j = izone->j0; j < izone->j1; j++)
5918 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5919 if (shift_diff < izone->shift0[dim])
5921 izone->shift0[dim] = shift_diff;
5923 if (shift_diff > izone->shift1[dim])
5925 izone->shift1[dim] = shift_diff;
5932 if (dd->comm->eDLB != edlbNO)
5934 snew(dd->comm->root, dd->ndim);
5937 if (dd->comm->bRecordLoad)
5939 make_load_communicators(dd);
5943 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
5949 gmx_domdec_comm_t *comm;
5956 if (comm->bCartesianPP)
5958 /* Set up cartesian communication for the particle-particle part */
5961 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5962 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5965 for (int i = 0; i < DIM; i++)
5969 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5971 /* We overwrite the old communicator with the new cartesian one */
5972 cr->mpi_comm_mygroup = comm_cart;
5975 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5976 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5978 if (comm->bCartesianPP_PME)
5980 /* Since we want to use the original cartesian setup for sim,
5981 * and not the one after split, we need to make an index.
5983 snew(comm->ddindex2ddnodeid, dd->nnodes);
5984 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5985 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5986 /* Get the rank of the DD master,
5987 * above we made sure that the master node is a PP node.
5997 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5999 else if (comm->bCartesianPP)
6001 if (cr->npmenodes == 0)
6003 /* The PP communicator is also
6004 * the communicator for this simulation
6006 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
6008 cr->nodeid = dd->rank;
6010 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
6012 /* We need to make an index to go from the coordinates
6013 * to the nodeid of this simulation.
6015 snew(comm->ddindex2simnodeid, dd->nnodes);
6016 snew(buf, dd->nnodes);
6017 if (cr->duty & DUTY_PP)
6019 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6021 /* Communicate the ddindex to simulation nodeid index */
6022 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6023 cr->mpi_comm_mysim);
6026 /* Determine the master coordinates and rank.
6027 * The DD master should be the same node as the master of this sim.
6029 for (int i = 0; i < dd->nnodes; i++)
6031 if (comm->ddindex2simnodeid[i] == 0)
6033 ddindex2xyz(dd->nc, i, dd->master_ci);
6034 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
6039 fprintf(debug, "The master rank is %d\n", dd->masterrank);
6044 /* No Cartesian communicators */
6045 /* We use the rank in dd->comm->all as DD index */
6046 ddindex2xyz(dd->nc, dd->rank, dd->ci);
6047 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
6049 clear_ivec(dd->master_ci);
6056 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6057 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6062 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
6063 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6067 static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
6071 gmx_domdec_comm_t *comm;
6076 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
6079 snew(comm->ddindex2simnodeid, dd->nnodes);
6080 snew(buf, dd->nnodes);
6081 if (cr->duty & DUTY_PP)
6083 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
6085 /* Communicate the ddindex to simulation nodeid index */
6086 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6087 cr->mpi_comm_mysim);
6093 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6094 int ncg, int natoms)
6096 gmx_domdec_master_t *ma;
6101 snew(ma->ncg, dd->nnodes);
6102 snew(ma->index, dd->nnodes+1);
6104 snew(ma->nat, dd->nnodes);
6105 snew(ma->ibuf, dd->nnodes*2);
6106 snew(ma->cell_x, DIM);
6107 for (i = 0; i < DIM; i++)
6109 snew(ma->cell_x[i], dd->nc[i]+1);
6112 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6118 snew(ma->vbuf, natoms);
6124 static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_node_order,
6125 int gmx_unused reorder)
6128 gmx_domdec_comm_t *comm;
6138 if (comm->bCartesianPP)
6140 for (i = 1; i < DIM; i++)
6142 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6144 if (bDiv[YY] || bDiv[ZZ])
6146 comm->bCartesianPP_PME = TRUE;
6147 /* If we have 2D PME decomposition, which is always in x+y,
6148 * we stack the PME only nodes in z.
6149 * Otherwise we choose the direction that provides the thinnest slab
6150 * of PME only nodes as this will have the least effect
6151 * on the PP communication.
6152 * But for the PME communication the opposite might be better.
6154 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6156 dd->nc[YY] > dd->nc[ZZ]))
6158 comm->cartpmedim = ZZ;
6162 comm->cartpmedim = YY;
6164 comm->ntot[comm->cartpmedim]
6165 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6169 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]);
6171 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6176 if (comm->bCartesianPP_PME)
6183 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]);
6186 for (i = 0; i < DIM; i++)
6190 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6192 MPI_Comm_rank(comm_cart, &rank);
6193 if (MASTERNODE(cr) && rank != 0)
6195 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6198 /* With this assigment we loose the link to the original communicator
6199 * which will usually be MPI_COMM_WORLD, unless have multisim.
6201 cr->mpi_comm_mysim = comm_cart;
6202 cr->sim_nodeid = rank;
6204 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6208 fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
6209 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6212 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6216 if (cr->npmenodes == 0 ||
6217 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6219 cr->duty = DUTY_PME;
6222 /* Split the sim communicator into PP and PME only nodes */
6223 MPI_Comm_split(cr->mpi_comm_mysim,
6225 dd_index(comm->ntot, dd->ci),
6226 &cr->mpi_comm_mygroup);
6230 switch (dd_node_order)
6235 fprintf(fplog, "Order of the ranks: PP first, PME last\n");
6238 case ddnoINTERLEAVE:
6239 /* Interleave the PP-only and PME-only nodes,
6240 * as on clusters with dual-core machines this will double
6241 * the communication bandwidth of the PME processes
6242 * and thus speed up the PP <-> PME and inter PME communication.
6246 fprintf(fplog, "Interleaving PP and PME ranks\n");
6248 comm->pmenodes = dd_pmenodes(cr);
6253 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6256 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6258 cr->duty = DUTY_PME;
6265 /* Split the sim communicator into PP and PME only nodes */
6266 MPI_Comm_split(cr->mpi_comm_mysim,
6269 &cr->mpi_comm_mygroup);
6270 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6276 fprintf(fplog, "This rank does only %s work.\n\n",
6277 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6281 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6284 gmx_domdec_comm_t *comm;
6290 copy_ivec(dd->nc, comm->ntot);
6292 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6293 comm->bCartesianPP_PME = FALSE;
6295 /* Reorder the nodes by default. This might change the MPI ranks.
6296 * Real reordering is only supported on very few architectures,
6297 * Blue Gene is one of them.
6299 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6301 if (cr->npmenodes > 0)
6303 /* Split the communicator into a PP and PME part */
6304 split_communicator(fplog, cr, dd_node_order, CartReorder);
6305 if (comm->bCartesianPP_PME)
6307 /* We (possibly) reordered the nodes in split_communicator,
6308 * so it is no longer required in make_pp_communicator.
6310 CartReorder = FALSE;
6315 /* All nodes do PP and PME */
6317 /* We do not require separate communicators */
6318 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6322 if (cr->duty & DUTY_PP)
6324 /* Copy or make a new PP communicator */
6325 make_pp_communicator(fplog, cr, CartReorder);
6329 receive_ddindex2simnodeid(cr);
6332 if (!(cr->duty & DUTY_PME))
6334 /* Set up the commnuication to our PME node */
6335 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6336 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6339 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6340 dd->pme_nodeid, dd->pme_receive_vir_ener);
6345 dd->pme_nodeid = -1;
6350 dd->ma = init_gmx_domdec_master_t(dd,
6352 comm->cgs_gl.index[comm->cgs_gl.nr]);
6356 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6358 real *slb_frac, tot;
6363 if (nc > 1 && size_string != NULL)
6367 fprintf(fplog, "Using static load balancing for the %s direction\n",
6372 for (i = 0; i < nc; i++)
6375 sscanf(size_string, "%20lf%n", &dbl, &n);
6378 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6387 fprintf(fplog, "Relative cell sizes:");
6389 for (i = 0; i < nc; i++)
6394 fprintf(fplog, " %5.3f", slb_frac[i]);
6399 fprintf(fplog, "\n");
6406 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6409 gmx_mtop_ilistloop_t iloop;
6413 iloop = gmx_mtop_ilistloop_init(mtop);
6414 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6416 for (ftype = 0; ftype < F_NRE; ftype++)
6418 if ((interaction_function[ftype].flags & IF_BOND) &&
6421 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6429 static int dd_getenv(FILE *fplog, const char *env_var, int def)
6435 val = getenv(env_var);
6438 if (sscanf(val, "%20d", &nst) <= 0)
6444 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6452 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6456 fprintf(stderr, "\n%s\n", warn_string);
6460 fprintf(fplog, "\n%s\n", warn_string);
6464 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6465 t_inputrec *ir, FILE *fplog)
6467 if (ir->ePBC == epbcSCREW &&
6468 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6470 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6473 if (ir->ns_type == ensSIMPLE)
6475 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
6478 if (ir->nstlist == 0)
6480 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6483 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6485 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6489 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6494 r = ddbox->box_size[XX];
6495 for (di = 0; di < dd->ndim; di++)
6498 /* Check using the initial average cell size */
6499 r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6505 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6506 const char *dlb_opt, gmx_bool bRecordLoad,
6507 unsigned long Flags, t_inputrec *ir)
6514 case 'a': eDLB = edlbAUTO; break;
6515 case 'n': eDLB = edlbNO; break;
6516 case 'y': eDLB = edlbYES; break;
6517 default: gmx_incons("Unknown dlb_opt");
6520 if (Flags & MD_RERUN)
6525 if (!EI_DYNAMICS(ir->eI))
6527 if (eDLB == edlbYES)
6529 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6530 dd_warning(cr, fplog, buf);
6538 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6543 if (Flags & MD_REPRODUCIBLE)
6550 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6554 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6557 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6565 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6570 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6572 /* Decomposition order z,y,x */
6575 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6577 for (dim = DIM-1; dim >= 0; dim--)
6579 if (dd->nc[dim] > 1)
6581 dd->dim[dd->ndim++] = dim;
6587 /* Decomposition order x,y,z */
6588 for (dim = 0; dim < DIM; dim++)
6590 if (dd->nc[dim] > 1)
6592 dd->dim[dd->ndim++] = dim;
6598 static gmx_domdec_comm_t *init_dd_comm()
6600 gmx_domdec_comm_t *comm;
6604 snew(comm->cggl_flag, DIM*2);
6605 snew(comm->cgcm_state, DIM*2);
6606 for (i = 0; i < DIM*2; i++)
6608 comm->cggl_flag_nalloc[i] = 0;
6609 comm->cgcm_state_nalloc[i] = 0;
6612 comm->nalloc_int = 0;
6613 comm->buf_int = NULL;
6615 vec_rvec_init(&comm->vbuf);
6617 comm->n_load_have = 0;
6618 comm->n_load_collect = 0;
6620 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6622 comm->sum_nat[i] = 0;
6626 comm->load_step = 0;
6629 clear_ivec(comm->load_lim);
6636 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6637 unsigned long Flags,
6639 real comm_distance_min, real rconstr,
6640 const char *dlb_opt, real dlb_scale,
6641 const char *sizex, const char *sizey, const char *sizez,
6642 gmx_mtop_t *mtop, t_inputrec *ir,
6643 matrix box, rvec *x,
6645 int *npme_x, int *npme_y)
6648 gmx_domdec_comm_t *comm;
6650 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6653 const real tenPercentMargin = 1.1;
6658 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
6663 dd->comm = init_dd_comm();
6665 snew(comm->cggl_flag, DIM*2);
6666 snew(comm->cgcm_state, DIM*2);
6668 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6669 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6671 dd->bSendRecv2 = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
6672 comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
6673 comm->eFlop = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
6674 recload = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
6675 comm->nstSortCG = dd_getenv(fplog, "GMX_DD_NST_SORT_CHARGE_GROUPS", 1);
6676 comm->nstDDDump = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
6677 comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
6678 comm->DD_debug = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
6680 dd->pme_recv_f_alloc = 0;
6681 dd->pme_recv_f_buf = NULL;
6683 if (dd->bSendRecv2 && fplog)
6685 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");
6691 fprintf(fplog, "Will load balance based on FLOP count\n");
6693 if (comm->eFlop > 1)
6695 srand(1+cr->nodeid);
6697 comm->bRecordLoad = TRUE;
6701 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6705 /* Initialize to GPU share count to 0, might change later */
6706 comm->nrank_gpu_shared = 0;
6708 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6709 comm->bDLB_locked = FALSE;
6710 comm->bCheckWhetherToTurnDlbOn = TRUE;
6712 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6715 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6717 dd->bGridJump = comm->bDynLoadBal;
6718 comm->bPMELoadBalDLBLimits = FALSE;
6720 if (comm->nstSortCG)
6724 if (comm->nstSortCG == 1)
6726 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6730 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6734 snew(comm->sort, 1);
6740 fprintf(fplog, "Will not sort the charge groups\n");
6744 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6746 comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
6747 mtop->bIntermolecularInteractions);
6748 if (comm->bInterCGBondeds)
6750 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6754 comm->bInterCGMultiBody = FALSE;
6757 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6758 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6760 if (ir->rlistlong == 0)
6762 /* Set the cut-off to some very large value,
6763 * so we don't need if statements everywhere in the code.
6764 * We use sqrt, since the cut-off is squared in some places.
6766 comm->cutoff = GMX_CUTOFF_INF;
6770 comm->cutoff = ir->rlistlong;
6772 comm->cutoff_mbody = 0;
6774 comm->cellsize_limit = 0;
6775 comm->bBondComm = FALSE;
6777 /* Atoms should be able to move by up to half the list buffer size (if > 0)
6778 * within nstlist steps. Since boundaries are allowed to displace by half
6779 * a cell size, DD cells should be at least the size of the list buffer.
6781 comm->cellsize_limit = std::max(comm->cellsize_limit,
6782 ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
6784 if (comm->bInterCGBondeds)
6786 if (comm_distance_min > 0)
6788 comm->cutoff_mbody = comm_distance_min;
6789 if (Flags & MD_DDBONDCOMM)
6791 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6795 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6797 r_bonded_limit = comm->cutoff_mbody;
6799 else if (ir->bPeriodicMols)
6801 /* Can not easily determine the required cut-off */
6802 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");
6803 comm->cutoff_mbody = comm->cutoff/2;
6804 r_bonded_limit = comm->cutoff_mbody;
6810 dd_bonded_cg_distance(fplog, mtop, ir, x, box,
6811 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6813 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6814 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6816 /* We use an initial margin of 10% for the minimum cell size,
6817 * except when we are just below the non-bonded cut-off.
6819 if (Flags & MD_DDBONDCOMM)
6821 if (std::max(r_2b, r_mb) > comm->cutoff)
6823 r_bonded = std::max(r_2b, r_mb);
6824 r_bonded_limit = tenPercentMargin*r_bonded;
6825 comm->bBondComm = TRUE;
6830 r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
6832 /* We determine cutoff_mbody later */
6836 /* No special bonded communication,
6837 * simply increase the DD cut-off.
6839 r_bonded_limit = tenPercentMargin*std::max(r_2b, r_mb);
6840 comm->cutoff_mbody = r_bonded_limit;
6841 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
6847 "Minimum cell size due to bonded interactions: %.3f nm\n",
6850 comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
6853 if (dd->bInterCGcons && rconstr <= 0)
6855 /* There is a cell size limit due to the constraints (P-LINCS) */
6856 rconstr = constr_r_max(fplog, mtop, ir);
6860 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6862 if (rconstr > comm->cellsize_limit)
6864 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6868 else if (rconstr > 0 && fplog)
6870 /* Here we do not check for dd->bInterCGcons,
6871 * because one can also set a cell size limit for virtual sites only
6872 * and at this point we don't know yet if there are intercg v-sites.
6875 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6878 comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
6880 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6884 copy_ivec(nc, dd->nc);
6885 set_dd_dim(fplog, dd);
6886 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6888 if (cr->npmenodes == -1)
6892 acs = average_cellsize_min(dd, ddbox);
6893 if (acs < comm->cellsize_limit)
6897 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6899 gmx_fatal_collective(FARGS, cr, NULL,
6900 "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",
6901 acs, comm->cellsize_limit);
6906 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6908 /* We need to choose the optimal DD grid and possibly PME nodes */
6909 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6910 comm->eDLB != edlbNO, dlb_scale,
6911 comm->cellsize_limit, comm->cutoff,
6912 comm->bInterCGBondeds);
6914 if (dd->nc[XX] == 0)
6916 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6917 sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
6918 !bC ? "-rdd" : "-rcon",
6919 comm->eDLB != edlbNO ? " or -dds" : "",
6920 bC ? " or your LINCS settings" : "");
6922 gmx_fatal_collective(FARGS, cr, NULL,
6923 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
6925 "Look in the log file for details on the domain decomposition",
6926 cr->nnodes-cr->npmenodes, limit, buf);
6928 set_dd_dim(fplog, dd);
6934 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
6935 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6938 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6939 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6941 gmx_fatal_collective(FARGS, cr, NULL,
6942 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
6943 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6945 if (cr->npmenodes > dd->nnodes)
6947 gmx_fatal_collective(FARGS, cr, NULL,
6948 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6950 if (cr->npmenodes > 0)
6952 comm->npmenodes = cr->npmenodes;
6956 comm->npmenodes = dd->nnodes;
6959 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
6961 /* The following choices should match those
6962 * in comm_cost_est in domdec_setup.c.
6963 * Note that here the checks have to take into account
6964 * that the decomposition might occur in a different order than xyz
6965 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6966 * in which case they will not match those in comm_cost_est,
6967 * but since that is mainly for testing purposes that's fine.
6969 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6970 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6971 getenv("GMX_PMEONEDD") == NULL)
6973 comm->npmedecompdim = 2;
6974 comm->npmenodes_x = dd->nc[XX];
6975 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6979 /* In case nc is 1 in both x and y we could still choose to
6980 * decompose pme in y instead of x, but we use x for simplicity.
6982 comm->npmedecompdim = 1;
6983 if (dd->dim[0] == YY)
6985 comm->npmenodes_x = 1;
6986 comm->npmenodes_y = comm->npmenodes;
6990 comm->npmenodes_x = comm->npmenodes;
6991 comm->npmenodes_y = 1;
6996 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6997 comm->npmenodes_x, comm->npmenodes_y, 1);
7002 comm->npmedecompdim = 0;
7003 comm->npmenodes_x = 0;
7004 comm->npmenodes_y = 0;
7007 /* Technically we don't need both of these,
7008 * but it simplifies code not having to recalculate it.
7010 *npme_x = comm->npmenodes_x;
7011 *npme_y = comm->npmenodes_y;
7013 snew(comm->slb_frac, DIM);
7014 if (comm->eDLB == edlbNO)
7016 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
7017 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
7018 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
7021 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
7023 if (comm->bBondComm || comm->eDLB != edlbNO)
7025 /* Set the bonded communication distance to halfway
7026 * the minimum and the maximum,
7027 * since the extra communication cost is nearly zero.
7029 acs = average_cellsize_min(dd, ddbox);
7030 comm->cutoff_mbody = 0.5*(r_bonded + acs);
7031 if (comm->eDLB != edlbNO)
7033 /* Check if this does not limit the scaling */
7034 comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
7036 if (!comm->bBondComm)
7038 /* Without bBondComm do not go beyond the n.b. cut-off */
7039 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
7040 if (comm->cellsize_limit >= comm->cutoff)
7042 /* We don't loose a lot of efficieny
7043 * when increasing it to the n.b. cut-off.
7044 * It can even be slightly faster, because we need
7045 * less checks for the communication setup.
7047 comm->cutoff_mbody = comm->cutoff;
7050 /* Check if we did not end up below our original limit */
7051 comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
7053 if (comm->cutoff_mbody > comm->cellsize_limit)
7055 comm->cellsize_limit = comm->cutoff_mbody;
7058 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
7063 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
7064 "cellsize limit %f\n",
7065 comm->bBondComm, comm->cellsize_limit);
7070 check_dd_restrictions(cr, dd, ir, fplog);
7073 comm->partition_step = INT_MIN;
7076 clear_dd_cycle_counts(dd);
7081 static void set_dlb_limits(gmx_domdec_t *dd)
7086 for (d = 0; d < dd->ndim; d++)
7088 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
7089 dd->comm->cellsize_min[dd->dim[d]] =
7090 dd->comm->cellsize_min_dlb[dd->dim[d]];
7095 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
7098 gmx_domdec_comm_t *comm;
7108 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);
7111 cellsize_min = comm->cellsize_min[dd->dim[0]];
7112 for (d = 1; d < dd->ndim; d++)
7114 cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7117 if (cellsize_min < comm->cellsize_limit*1.05)
7119 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");
7121 /* Change DLB from "auto" to "no". */
7122 comm->eDLB = edlbNO;
7127 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7128 comm->bDynLoadBal = TRUE;
7129 dd->bGridJump = TRUE;
7133 /* We can set the required cell size info here,
7134 * so we do not need to communicate this.
7135 * The grid is completely uniform.
7137 for (d = 0; d < dd->ndim; d++)
7141 comm->load[d].sum_m = comm->load[d].sum;
7143 nc = dd->nc[dd->dim[d]];
7144 for (i = 0; i < nc; i++)
7146 comm->root[d]->cell_f[i] = i/(real)nc;
7149 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7150 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7153 comm->root[d]->cell_f[nc] = 1.0;
7158 static char *init_bLocalCG(gmx_mtop_t *mtop)
7163 ncg = ncg_mtop(mtop);
7164 snew(bLocalCG, ncg);
7165 for (cg = 0; cg < ncg; cg++)
7167 bLocalCG[cg] = FALSE;
7173 void dd_init_bondeds(FILE *fplog,
7174 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7176 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7178 gmx_domdec_comm_t *comm;
7180 dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
7184 if (comm->bBondComm)
7186 /* Communicate atoms beyond the cut-off for bonded interactions */
7189 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7191 comm->bLocalCG = init_bLocalCG(mtop);
7195 /* Only communicate atoms based on cut-off */
7196 comm->cglink = NULL;
7197 comm->bLocalCG = NULL;
7201 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7203 gmx_bool bDynLoadBal, real dlb_scale,
7206 gmx_domdec_comm_t *comm;
7221 fprintf(fplog, "The maximum number of communication pulses is:");
7222 for (d = 0; d < dd->ndim; d++)
7224 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7226 fprintf(fplog, "\n");
7227 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7228 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7229 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7230 for (d = 0; d < DIM; d++)
7234 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7241 comm->cellsize_min_dlb[d]/
7242 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7244 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7247 fprintf(fplog, "\n");
7251 set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
7252 fprintf(fplog, "The initial number of communication pulses is:");
7253 for (d = 0; d < dd->ndim; d++)
7255 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7257 fprintf(fplog, "\n");
7258 fprintf(fplog, "The initial domain decomposition cell size is:");
7259 for (d = 0; d < DIM; d++)
7263 fprintf(fplog, " %c %.2f nm",
7264 dim2char(d), dd->comm->cellsize_min[d]);
7267 fprintf(fplog, "\n\n");
7270 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7272 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7273 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7274 "non-bonded interactions", "", comm->cutoff);
7278 limit = dd->comm->cellsize_limit;
7282 if (dynamic_dd_box(ddbox, ir))
7284 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7286 limit = dd->comm->cellsize_min[XX];
7287 for (d = 1; d < DIM; d++)
7289 limit = std::min(limit, dd->comm->cellsize_min[d]);
7293 if (comm->bInterCGBondeds)
7295 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7296 "two-body bonded interactions", "(-rdd)",
7297 std::max(comm->cutoff, comm->cutoff_mbody));
7298 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7299 "multi-body bonded interactions", "(-rdd)",
7300 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
7304 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7305 "virtual site constructions", "(-rcon)", limit);
7307 if (dd->constraint_comm)
7309 sprintf(buf, "atoms separated by up to %d constraints",
7311 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7312 buf, "(-rcon)", limit);
7314 fprintf(fplog, "\n");
7320 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7322 const t_inputrec *ir,
7323 const gmx_ddbox_t *ddbox)
7325 gmx_domdec_comm_t *comm;
7326 int d, dim, npulse, npulse_d_max, npulse_d;
7331 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7333 /* Determine the maximum number of comm. pulses in one dimension */
7335 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7337 /* Determine the maximum required number of grid pulses */
7338 if (comm->cellsize_limit >= comm->cutoff)
7340 /* Only a single pulse is required */
7343 else if (!bNoCutOff && comm->cellsize_limit > 0)
7345 /* We round down slightly here to avoid overhead due to the latency
7346 * of extra communication calls when the cut-off
7347 * would be only slightly longer than the cell size.
7348 * Later cellsize_limit is redetermined,
7349 * so we can not miss interactions due to this rounding.
7351 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7355 /* There is no cell size limit */
7356 npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7359 if (!bNoCutOff && npulse > 1)
7361 /* See if we can do with less pulses, based on dlb_scale */
7363 for (d = 0; d < dd->ndim; d++)
7366 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7367 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7368 npulse_d_max = std::max(npulse_d_max, npulse_d);
7370 npulse = std::min(npulse, npulse_d_max);
7373 /* This env var can override npulse */
7374 d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
7381 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7382 for (d = 0; d < dd->ndim; d++)
7384 comm->cd[d].np_dlb = std::min(npulse, dd->nc[dd->dim[d]]-1);
7385 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7386 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7387 comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
7388 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7390 comm->bVacDLBNoLimit = FALSE;
7394 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7395 if (!comm->bVacDLBNoLimit)
7397 comm->cellsize_limit = std::max(comm->cellsize_limit,
7398 comm->cutoff/comm->maxpulse);
7400 comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
7401 /* Set the minimum cell size for each DD dimension */
7402 for (d = 0; d < dd->ndim; d++)
7404 if (comm->bVacDLBNoLimit ||
7405 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7407 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7411 comm->cellsize_min_dlb[dd->dim[d]] =
7412 comm->cutoff/comm->cd[d].np_dlb;
7415 if (comm->cutoff_mbody <= 0)
7417 comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
7419 if (comm->bDynLoadBal)
7425 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7427 /* If each molecule is a single charge group
7428 * or we use domain decomposition for each periodic dimension,
7429 * we do not need to take pbc into account for the bonded interactions.
7431 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7434 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7437 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7438 t_inputrec *ir, gmx_ddbox_t *ddbox)
7440 gmx_domdec_comm_t *comm;
7446 /* Initialize the thread data.
7447 * This can not be done in init_domain_decomposition,
7448 * as the numbers of threads is determined later.
7450 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7453 snew(comm->dth, comm->nth);
7456 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
7458 init_ddpme(dd, &comm->ddpme[0], 0);
7459 if (comm->npmedecompdim >= 2)
7461 init_ddpme(dd, &comm->ddpme[1], 1);
7466 comm->npmenodes = 0;
7467 if (dd->pme_nodeid >= 0)
7469 gmx_fatal_collective(FARGS, NULL, dd,
7470 "Can not have separate PME ranks without PME electrostatics");
7476 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7478 if (comm->eDLB != edlbNO)
7480 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7483 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7484 if (comm->eDLB == edlbAUTO)
7488 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7490 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7493 if (ir->ePBC == epbcNONE)
7495 vol_frac = 1 - 1/(double)dd->nnodes;
7500 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7504 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7506 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7508 dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
7511 static gmx_bool test_dd_cutoff(t_commrec *cr,
7512 t_state *state, t_inputrec *ir,
7523 set_ddbox(dd, FALSE, cr, ir, state->box,
7524 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7528 for (d = 0; d < dd->ndim; d++)
7532 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7533 if (dynamic_dd_box(&ddbox, ir))
7535 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7538 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7540 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7541 dd->comm->cd[d].np_dlb > 0)
7543 if (np > dd->comm->cd[d].np_dlb)
7548 /* If a current local cell size is smaller than the requested
7549 * cut-off, we could still fix it, but this gets very complicated.
7550 * Without fixing here, we might actually need more checks.
7552 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7559 if (dd->comm->eDLB != edlbNO)
7561 /* If DLB is not active yet, we don't need to check the grid jumps.
7562 * Actually we shouldn't, because then the grid jump data is not set.
7564 if (dd->comm->bDynLoadBal &&
7565 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7570 gmx_sumi(1, &LocallyLimited, cr);
7572 if (LocallyLimited > 0)
7581 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7584 gmx_bool bCutoffAllowed;
7586 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7590 cr->dd->comm->cutoff = cutoff_req;
7593 return bCutoffAllowed;
7596 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
7598 gmx_domdec_comm_t *comm;
7600 comm = cr->dd->comm;
7602 /* Turn on the DLB limiting (might have been on already) */
7603 comm->bPMELoadBalDLBLimits = TRUE;
7605 /* Change the cut-off limit */
7606 comm->PMELoadBal_max_cutoff = cutoff;
7610 fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
7611 comm->PMELoadBal_max_cutoff);
7615 /* Sets whether we should later check the load imbalance data, so that
7616 * we can trigger dynamic load balancing if enough imbalance has
7619 * Used after PME load balancing unlocks DLB, so that the check
7620 * whether DLB will be useful can happen immediately.
7622 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
7624 if (dd->comm->eDLB == edlbAUTO && !dd_dlb_is_locked(dd))
7626 dd->comm->bCheckWhetherToTurnDlbOn = bValue;
7630 /* Returns if we should check whether there has been enough load
7631 * imbalance to trigger dynamic load balancing.
7633 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
7635 const int nddp_chk_dlb = 100;
7637 if (dd->comm->eDLB != edlbAUTO)
7642 /* We should check whether we should use DLB directly after
7644 if (dd->comm->bCheckWhetherToTurnDlbOn)
7646 /* This flag was set when the PME load-balancing routines
7647 unlocked DLB, and should now be cleared. */
7648 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
7651 /* We should also check whether we should use DLB every 100
7652 * partitionings (we do not do this every partioning, so that we
7653 * avoid excessive communication). */
7654 if (dd->comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1)
7662 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
7664 return dd->comm->bDynLoadBal;
7667 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
7669 return dd->comm->bDLB_locked;
7672 void dd_dlb_lock(gmx_domdec_t *dd)
7674 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7675 if (dd->comm->eDLB == edlbAUTO)
7677 dd->comm->bDLB_locked = TRUE;
7681 void dd_dlb_unlock(gmx_domdec_t *dd)
7683 /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
7684 if (dd->comm->eDLB == edlbAUTO)
7686 dd->comm->bDLB_locked = FALSE;
7687 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, !dd->comm->bDynLoadBal);
7691 static void merge_cg_buffers(int ncell,
7692 gmx_domdec_comm_dim_t *cd, int pulse,
7694 int *index_gl, int *recv_i,
7695 rvec *cg_cm, rvec *recv_vr,
7697 cginfo_mb_t *cginfo_mb, int *cginfo)
7699 gmx_domdec_ind_t *ind, *ind_p;
7700 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7701 int shift, shift_at;
7703 ind = &cd->ind[pulse];
7705 /* First correct the already stored data */
7706 shift = ind->nrecv[ncell];
7707 for (cell = ncell-1; cell >= 0; cell--)
7709 shift -= ind->nrecv[cell];
7712 /* Move the cg's present from previous grid pulses */
7713 cg0 = ncg_cell[ncell+cell];
7714 cg1 = ncg_cell[ncell+cell+1];
7715 cgindex[cg1+shift] = cgindex[cg1];
7716 for (cg = cg1-1; cg >= cg0; cg--)
7718 index_gl[cg+shift] = index_gl[cg];
7719 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7720 cgindex[cg+shift] = cgindex[cg];
7721 cginfo[cg+shift] = cginfo[cg];
7723 /* Correct the already stored send indices for the shift */
7724 for (p = 1; p <= pulse; p++)
7726 ind_p = &cd->ind[p];
7728 for (c = 0; c < cell; c++)
7730 cg0 += ind_p->nsend[c];
7732 cg1 = cg0 + ind_p->nsend[cell];
7733 for (cg = cg0; cg < cg1; cg++)
7735 ind_p->index[cg] += shift;
7741 /* Merge in the communicated buffers */
7745 for (cell = 0; cell < ncell; cell++)
7747 cg1 = ncg_cell[ncell+cell+1] + shift;
7750 /* Correct the old cg indices */
7751 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7753 cgindex[cg+1] += shift_at;
7756 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7758 /* Copy this charge group from the buffer */
7759 index_gl[cg1] = recv_i[cg0];
7760 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7761 /* Add it to the cgindex */
7762 cg_gl = index_gl[cg1];
7763 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7764 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7765 cgindex[cg1+1] = cgindex[cg1] + nat;
7770 shift += ind->nrecv[cell];
7771 ncg_cell[ncell+cell+1] = cg1;
7775 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7776 int nzone, int cg0, const int *cgindex)
7780 /* Store the atom block boundaries for easy copying of communication buffers
7783 for (zone = 0; zone < nzone; zone++)
7785 for (p = 0; p < cd->np; p++)
7787 cd->ind[p].cell2at0[zone] = cgindex[cg];
7788 cg += cd->ind[p].nrecv[zone];
7789 cd->ind[p].cell2at1[zone] = cgindex[cg];
7794 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7800 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7802 if (!bLocalCG[link->a[i]])
7811 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7813 real c[DIM][4]; /* the corners for the non-bonded communication */
7814 real cr0; /* corner for rounding */
7815 real cr1[4]; /* corners for rounding */
7816 real bc[DIM]; /* corners for bounded communication */
7817 real bcr1; /* corner for rounding for bonded communication */
7820 /* Determine the corners of the domain(s) we are communicating with */
7822 set_dd_corners(const gmx_domdec_t *dd,
7823 int dim0, int dim1, int dim2,
7827 const gmx_domdec_comm_t *comm;
7828 const gmx_domdec_zones_t *zones;
7833 zones = &comm->zones;
7835 /* Keep the compiler happy */
7839 /* The first dimension is equal for all cells */
7840 c->c[0][0] = comm->cell_x0[dim0];
7843 c->bc[0] = c->c[0][0];
7848 /* This cell row is only seen from the first row */
7849 c->c[1][0] = comm->cell_x0[dim1];
7850 /* All rows can see this row */
7851 c->c[1][1] = comm->cell_x0[dim1];
7854 c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7857 /* For the multi-body distance we need the maximum */
7858 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7861 /* Set the upper-right corner for rounding */
7862 c->cr0 = comm->cell_x1[dim0];
7867 for (j = 0; j < 4; j++)
7869 c->c[2][j] = comm->cell_x0[dim2];
7873 /* Use the maximum of the i-cells that see a j-cell */
7874 for (i = 0; i < zones->nizone; i++)
7876 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7881 std::max(c->c[2][j-4],
7882 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7888 /* For the multi-body distance we need the maximum */
7889 c->bc[2] = comm->cell_x0[dim2];
7890 for (i = 0; i < 2; i++)
7892 for (j = 0; j < 2; j++)
7894 c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
7900 /* Set the upper-right corner for rounding */
7901 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7902 * Only cell (0,0,0) can see cell 7 (1,1,1)
7904 c->cr1[0] = comm->cell_x1[dim1];
7905 c->cr1[3] = comm->cell_x1[dim1];
7908 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7911 /* For the multi-body distance we need the maximum */
7912 c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7919 /* Determine which cg's we need to send in this pulse from this zone */
7921 get_zone_pulse_cgs(gmx_domdec_t *dd,
7922 int zonei, int zone,
7924 const int *index_gl,
7926 int dim, int dim_ind,
7927 int dim0, int dim1, int dim2,
7928 real r_comm2, real r_bcomm2,
7932 real skew_fac2_d, real skew_fac_01,
7933 rvec *v_d, rvec *v_0, rvec *v_1,
7934 const dd_corners_t *c,
7936 gmx_bool bDistBonded,
7942 gmx_domdec_ind_t *ind,
7943 int **ibuf, int *ibuf_nalloc,
7949 gmx_domdec_comm_t *comm;
7951 gmx_bool bDistMB_pulse;
7953 real r2, rb2, r, tric_sh;
7956 int nsend_z, nsend, nat;
7960 bScrew = (dd->bScrewPBC && dim == XX);
7962 bDistMB_pulse = (bDistMB && bDistBonded);
7968 for (cg = cg0; cg < cg1; cg++)
7972 if (tric_dist[dim_ind] == 0)
7974 /* Rectangular direction, easy */
7975 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7982 r = cg_cm[cg][dim] - c->bc[dim_ind];
7988 /* Rounding gives at most a 16% reduction
7989 * in communicated atoms
7991 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7993 r = cg_cm[cg][dim0] - c->cr0;
7994 /* This is the first dimension, so always r >= 0 */
8001 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8003 r = cg_cm[cg][dim1] - c->cr1[zone];
8010 r = cg_cm[cg][dim1] - c->bcr1;
8020 /* Triclinic direction, more complicated */
8023 /* Rounding, conservative as the skew_fac multiplication
8024 * will slightly underestimate the distance.
8026 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
8028 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
8029 for (i = dim0+1; i < DIM; i++)
8031 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
8033 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
8036 rb[dim0] = rn[dim0];
8039 /* Take care that the cell planes along dim0 might not
8040 * be orthogonal to those along dim1 and dim2.
8042 for (i = 1; i <= dim_ind; i++)
8045 if (normal[dim0][dimd] > 0)
8047 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
8050 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
8055 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
8057 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
8059 for (i = dim1+1; i < DIM; i++)
8061 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
8063 rn[dim1] += tric_sh;
8066 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
8067 /* Take care of coupling of the distances
8068 * to the planes along dim0 and dim1 through dim2.
8070 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
8071 /* Take care that the cell planes along dim1
8072 * might not be orthogonal to that along dim2.
8074 if (normal[dim1][dim2] > 0)
8076 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
8082 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
8085 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
8086 /* Take care of coupling of the distances
8087 * to the planes along dim0 and dim1 through dim2.
8089 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
8090 /* Take care that the cell planes along dim1
8091 * might not be orthogonal to that along dim2.
8093 if (normal[dim1][dim2] > 0)
8095 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
8100 /* The distance along the communication direction */
8101 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
8103 for (i = dim+1; i < DIM; i++)
8105 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
8110 r2 += rn[dim]*rn[dim]*skew_fac2_d;
8111 /* Take care of coupling of the distances
8112 * to the planes along dim0 and dim1 through dim2.
8114 if (dim_ind == 1 && zonei == 1)
8116 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
8122 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
8125 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
8126 /* Take care of coupling of the distances
8127 * to the planes along dim0 and dim1 through dim2.
8129 if (dim_ind == 1 && zonei == 1)
8131 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
8139 ((bDistMB && rb2 < r_bcomm2) ||
8140 (bDist2B && r2 < r_bcomm2)) &&
8142 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
8143 missing_link(comm->cglink, index_gl[cg],
8146 /* Make an index to the local charge groups */
8147 if (nsend+1 > ind->nalloc)
8149 ind->nalloc = over_alloc_large(nsend+1);
8150 srenew(ind->index, ind->nalloc);
8152 if (nsend+1 > *ibuf_nalloc)
8154 *ibuf_nalloc = over_alloc_large(nsend+1);
8155 srenew(*ibuf, *ibuf_nalloc);
8157 ind->index[nsend] = cg;
8158 (*ibuf)[nsend] = index_gl[cg];
8160 vec_rvec_check_alloc(vbuf, nsend+1);
8162 if (dd->ci[dim] == 0)
8164 /* Correct cg_cm for pbc */
8165 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
8168 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
8169 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
8174 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
8177 nat += cgindex[cg+1] - cgindex[cg];
8183 *nsend_z_ptr = nsend_z;
8186 static void setup_dd_communication(gmx_domdec_t *dd,
8187 matrix box, gmx_ddbox_t *ddbox,
8188 t_forcerec *fr, t_state *state, rvec **f)
8190 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8191 int nzone, nzone_send, zone, zonei, cg0, cg1;
8192 int c, i, cg, cg_gl, nrcg;
8193 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8194 gmx_domdec_comm_t *comm;
8195 gmx_domdec_zones_t *zones;
8196 gmx_domdec_comm_dim_t *cd;
8197 gmx_domdec_ind_t *ind;
8198 cginfo_mb_t *cginfo_mb;
8199 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8200 real r_comm2, r_bcomm2;
8201 dd_corners_t corners;
8203 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8204 real skew_fac2_d, skew_fac_01;
8211 fprintf(debug, "Setting up DD communication\n");
8216 switch (fr->cutoff_scheme)
8225 gmx_incons("unimplemented");
8229 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8231 /* Check if we need to use triclinic distances */
8232 tric_dist[dim_ind] = 0;
8233 for (i = 0; i <= dim_ind; i++)
8235 if (ddbox->tric_dir[dd->dim[i]])
8237 tric_dist[dim_ind] = 1;
8242 bBondComm = comm->bBondComm;
8244 /* Do we need to determine extra distances for multi-body bondeds? */
8245 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8247 /* Do we need to determine extra distances for only two-body bondeds? */
8248 bDist2B = (bBondComm && !bDistMB);
8250 r_comm2 = sqr(comm->cutoff);
8251 r_bcomm2 = sqr(comm->cutoff_mbody);
8255 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8258 zones = &comm->zones;
8261 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8262 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8264 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8266 /* Triclinic stuff */
8267 normal = ddbox->normal;
8271 v_0 = ddbox->v[dim0];
8272 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8274 /* Determine the coupling coefficient for the distances
8275 * to the cell planes along dim0 and dim1 through dim2.
8276 * This is required for correct rounding.
8279 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8282 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8288 v_1 = ddbox->v[dim1];
8291 zone_cg_range = zones->cg_range;
8292 index_gl = dd->index_gl;
8293 cgindex = dd->cgindex;
8294 cginfo_mb = fr->cginfo_mb;
8296 zone_cg_range[0] = 0;
8297 zone_cg_range[1] = dd->ncg_home;
8298 comm->zone_ncg1[0] = dd->ncg_home;
8299 pos_cg = dd->ncg_home;
8301 nat_tot = dd->nat_home;
8303 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8305 dim = dd->dim[dim_ind];
8306 cd = &comm->cd[dim_ind];
8308 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8310 /* No pbc in this dimension, the first node should not comm. */
8318 v_d = ddbox->v[dim];
8319 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8321 cd->bInPlace = TRUE;
8322 for (p = 0; p < cd->np; p++)
8324 /* Only atoms communicated in the first pulse are used
8325 * for multi-body bonded interactions or for bBondComm.
8327 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8332 for (zone = 0; zone < nzone_send; zone++)
8334 if (tric_dist[dim_ind] && dim_ind > 0)
8336 /* Determine slightly more optimized skew_fac's
8338 * This reduces the number of communicated atoms
8339 * by about 10% for 3D DD of rhombic dodecahedra.
8341 for (dimd = 0; dimd < dim; dimd++)
8343 sf2_round[dimd] = 1;
8344 if (ddbox->tric_dir[dimd])
8346 for (i = dd->dim[dimd]+1; i < DIM; i++)
8348 /* If we are shifted in dimension i
8349 * and the cell plane is tilted forward
8350 * in dimension i, skip this coupling.
8352 if (!(zones->shift[nzone+zone][i] &&
8353 ddbox->v[dimd][i][dimd] >= 0))
8356 sqr(ddbox->v[dimd][i][dimd]);
8359 sf2_round[dimd] = 1/sf2_round[dimd];
8364 zonei = zone_perm[dim_ind][zone];
8367 /* Here we permutate the zones to obtain a convenient order
8368 * for neighbor searching
8370 cg0 = zone_cg_range[zonei];
8371 cg1 = zone_cg_range[zonei+1];
8375 /* Look only at the cg's received in the previous grid pulse
8377 cg1 = zone_cg_range[nzone+zone+1];
8378 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8381 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8382 for (th = 0; th < comm->nth; th++)
8384 gmx_domdec_ind_t *ind_p;
8385 int **ibuf_p, *ibuf_nalloc_p;
8387 int *nsend_p, *nat_p;
8393 /* Thread 0 writes in the comm buffers */
8395 ibuf_p = &comm->buf_int;
8396 ibuf_nalloc_p = &comm->nalloc_int;
8397 vbuf_p = &comm->vbuf;
8400 nsend_zone_p = &ind->nsend[zone];
8404 /* Other threads write into temp buffers */
8405 ind_p = &comm->dth[th].ind;
8406 ibuf_p = &comm->dth[th].ibuf;
8407 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8408 vbuf_p = &comm->dth[th].vbuf;
8409 nsend_p = &comm->dth[th].nsend;
8410 nat_p = &comm->dth[th].nat;
8411 nsend_zone_p = &comm->dth[th].nsend_zone;
8413 comm->dth[th].nsend = 0;
8414 comm->dth[th].nat = 0;
8415 comm->dth[th].nsend_zone = 0;
8425 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8426 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8429 /* Get the cg's for this pulse in this zone */
8430 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8432 dim, dim_ind, dim0, dim1, dim2,
8435 normal, skew_fac2_d, skew_fac_01,
8436 v_d, v_0, v_1, &corners, sf2_round,
8437 bDistBonded, bBondComm,
8441 ibuf_p, ibuf_nalloc_p,
8447 /* Append data of threads>=1 to the communication buffers */
8448 for (th = 1; th < comm->nth; th++)
8450 dd_comm_setup_work_t *dth;
8453 dth = &comm->dth[th];
8455 ns1 = nsend + dth->nsend_zone;
8456 if (ns1 > ind->nalloc)
8458 ind->nalloc = over_alloc_dd(ns1);
8459 srenew(ind->index, ind->nalloc);
8461 if (ns1 > comm->nalloc_int)
8463 comm->nalloc_int = over_alloc_dd(ns1);
8464 srenew(comm->buf_int, comm->nalloc_int);
8466 if (ns1 > comm->vbuf.nalloc)
8468 comm->vbuf.nalloc = over_alloc_dd(ns1);
8469 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8472 for (i = 0; i < dth->nsend_zone; i++)
8474 ind->index[nsend] = dth->ind.index[i];
8475 comm->buf_int[nsend] = dth->ibuf[i];
8476 copy_rvec(dth->vbuf.v[i],
8477 comm->vbuf.v[nsend]);
8481 ind->nsend[zone] += dth->nsend_zone;
8484 /* Clear the counts in case we do not have pbc */
8485 for (zone = nzone_send; zone < nzone; zone++)
8487 ind->nsend[zone] = 0;
8489 ind->nsend[nzone] = nsend;
8490 ind->nsend[nzone+1] = nat;
8491 /* Communicate the number of cg's and atoms to receive */
8492 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8493 ind->nsend, nzone+2,
8494 ind->nrecv, nzone+2);
8496 /* The rvec buffer is also required for atom buffers of size nsend
8497 * in dd_move_x and dd_move_f.
8499 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8503 /* We can receive in place if only the last zone is not empty */
8504 for (zone = 0; zone < nzone-1; zone++)
8506 if (ind->nrecv[zone] > 0)
8508 cd->bInPlace = FALSE;
8513 /* The int buffer is only required here for the cg indices */
8514 if (ind->nrecv[nzone] > comm->nalloc_int2)
8516 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8517 srenew(comm->buf_int2, comm->nalloc_int2);
8519 /* The rvec buffer is also required for atom buffers
8520 * of size nrecv in dd_move_x and dd_move_f.
8522 i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8523 vec_rvec_check_alloc(&comm->vbuf2, i);
8527 /* Make space for the global cg indices */
8528 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8529 || dd->cg_nalloc == 0)
8531 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8532 srenew(index_gl, dd->cg_nalloc);
8533 srenew(cgindex, dd->cg_nalloc+1);
8535 /* Communicate the global cg indices */
8538 recv_i = index_gl + pos_cg;
8542 recv_i = comm->buf_int2;
8544 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8545 comm->buf_int, nsend,
8546 recv_i, ind->nrecv[nzone]);
8548 /* Make space for cg_cm */
8549 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8550 if (fr->cutoff_scheme == ecutsGROUP)
8558 /* Communicate cg_cm */
8561 recv_vr = cg_cm + pos_cg;
8565 recv_vr = comm->vbuf2.v;
8567 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8568 comm->vbuf.v, nsend,
8569 recv_vr, ind->nrecv[nzone]);
8571 /* Make the charge group index */
8574 zone = (p == 0 ? 0 : nzone - 1);
8575 while (zone < nzone)
8577 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8579 cg_gl = index_gl[pos_cg];
8580 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8581 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8582 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8585 /* Update the charge group presence,
8586 * so we can use it in the next pass of the loop.
8588 comm->bLocalCG[cg_gl] = TRUE;
8594 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8597 zone_cg_range[nzone+zone] = pos_cg;
8602 /* This part of the code is never executed with bBondComm. */
8603 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8604 index_gl, recv_i, cg_cm, recv_vr,
8605 cgindex, fr->cginfo_mb, fr->cginfo);
8606 pos_cg += ind->nrecv[nzone];
8608 nat_tot += ind->nrecv[nzone+1];
8612 /* Store the atom block for easy copying of communication buffers */
8613 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8617 dd->index_gl = index_gl;
8618 dd->cgindex = cgindex;
8620 dd->ncg_tot = zone_cg_range[zones->n];
8621 dd->nat_tot = nat_tot;
8622 comm->nat[ddnatHOME] = dd->nat_home;
8623 for (i = ddnatZONE; i < ddnatNR; i++)
8625 comm->nat[i] = dd->nat_tot;
8630 /* We don't need to update cginfo, since that was alrady done above.
8631 * So we pass NULL for the forcerec.
8633 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8634 NULL, comm->bLocalCG);
8639 fprintf(debug, "Finished setting up DD communication, zones:");
8640 for (c = 0; c < zones->n; c++)
8642 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8644 fprintf(debug, "\n");
8648 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8652 for (c = 0; c < zones->nizone; c++)
8654 zones->izone[c].cg1 = zones->cg_range[c+1];
8655 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8656 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8660 static void set_zones_size(gmx_domdec_t *dd,
8661 matrix box, const gmx_ddbox_t *ddbox,
8662 int zone_start, int zone_end)
8664 gmx_domdec_comm_t *comm;
8665 gmx_domdec_zones_t *zones;
8674 zones = &comm->zones;
8676 /* Do we need to determine extra distances for multi-body bondeds? */
8677 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8679 for (z = zone_start; z < zone_end; z++)
8681 /* Copy cell limits to zone limits.
8682 * Valid for non-DD dims and non-shifted dims.
8684 copy_rvec(comm->cell_x0, zones->size[z].x0);
8685 copy_rvec(comm->cell_x1, zones->size[z].x1);
8688 for (d = 0; d < dd->ndim; d++)
8692 for (z = 0; z < zones->n; z++)
8694 /* With a staggered grid we have different sizes
8695 * for non-shifted dimensions.
8697 if (dd->bGridJump && zones->shift[z][dim] == 0)
8701 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8702 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8706 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8707 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8713 rcmbs = comm->cutoff_mbody;
8714 if (ddbox->tric_dir[dim])
8716 rcs /= ddbox->skew_fac[dim];
8717 rcmbs /= ddbox->skew_fac[dim];
8720 /* Set the lower limit for the shifted zone dimensions */
8721 for (z = zone_start; z < zone_end; z++)
8723 if (zones->shift[z][dim] > 0)
8726 if (!dd->bGridJump || d == 0)
8728 zones->size[z].x0[dim] = comm->cell_x1[dim];
8729 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8733 /* Here we take the lower limit of the zone from
8734 * the lowest domain of the zone below.
8738 zones->size[z].x0[dim] =
8739 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8745 zones->size[z].x0[dim] =
8746 zones->size[zone_perm[2][z-4]].x0[dim];
8750 zones->size[z].x0[dim] =
8751 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8754 /* A temporary limit, is updated below */
8755 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8759 for (zi = 0; zi < zones->nizone; zi++)
8761 if (zones->shift[zi][dim] == 0)
8763 /* This takes the whole zone into account.
8764 * With multiple pulses this will lead
8765 * to a larger zone then strictly necessary.
8767 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8768 zones->size[zi].x1[dim]+rcmbs);
8776 /* Loop over the i-zones to set the upper limit of each
8779 for (zi = 0; zi < zones->nizone; zi++)
8781 if (zones->shift[zi][dim] == 0)
8783 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8785 if (zones->shift[z][dim] > 0)
8787 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
8788 zones->size[zi].x1[dim]+rcs);
8795 for (z = zone_start; z < zone_end; z++)
8797 /* Initialization only required to keep the compiler happy */
8798 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8801 /* To determine the bounding box for a zone we need to find
8802 * the extreme corners of 4, 2 or 1 corners.
8804 nc = 1 << (ddbox->nboundeddim - 1);
8806 for (c = 0; c < nc; c++)
8808 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8812 corner[YY] = zones->size[z].x0[YY];
8816 corner[YY] = zones->size[z].x1[YY];
8820 corner[ZZ] = zones->size[z].x0[ZZ];
8824 corner[ZZ] = zones->size[z].x1[ZZ];
8826 if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
8827 box[ZZ][1 - dd->dim[0]] != 0)
8829 /* With 1D domain decomposition the cg's are not in
8830 * the triclinic box, but triclinic x-y and rectangular y/x-z.
8831 * Shift the corner of the z-vector back to along the box
8832 * vector of dimension d, so it will later end up at 0 along d.
8833 * This can affect the location of this corner along dd->dim[0]
8834 * through the matrix operation below if box[d][dd->dim[0]]!=0.
8836 int d = 1 - dd->dim[0];
8838 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
8840 /* Apply the triclinic couplings */
8841 assert(ddbox->npbcdim <= DIM);
8842 for (i = YY; i < ddbox->npbcdim; i++)
8844 for (j = XX; j < i; j++)
8846 corner[j] += corner[i]*box[i][j]/box[i][i];
8851 copy_rvec(corner, corner_min);
8852 copy_rvec(corner, corner_max);
8856 for (i = 0; i < DIM; i++)
8858 corner_min[i] = std::min(corner_min[i], corner[i]);
8859 corner_max[i] = std::max(corner_max[i], corner[i]);
8863 /* Copy the extreme cornes without offset along x */
8864 for (i = 0; i < DIM; i++)
8866 zones->size[z].bb_x0[i] = corner_min[i];
8867 zones->size[z].bb_x1[i] = corner_max[i];
8869 /* Add the offset along x */
8870 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8871 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8874 if (zone_start == 0)
8877 for (dim = 0; dim < DIM; dim++)
8879 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8881 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8886 for (z = zone_start; z < zone_end; z++)
8888 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8890 zones->size[z].x0[XX], zones->size[z].x1[XX],
8891 zones->size[z].x0[YY], zones->size[z].x1[YY],
8892 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8893 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8895 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8896 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8897 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8902 static int comp_cgsort(const void *a, const void *b)
8906 gmx_cgsort_t *cga, *cgb;
8907 cga = (gmx_cgsort_t *)a;
8908 cgb = (gmx_cgsort_t *)b;
8910 comp = cga->nsc - cgb->nsc;
8913 comp = cga->ind_gl - cgb->ind_gl;
8919 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8924 /* Order the data */
8925 for (i = 0; i < n; i++)
8927 buf[i] = a[sort[i].ind];
8930 /* Copy back to the original array */
8931 for (i = 0; i < n; i++)
8937 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8942 /* Order the data */
8943 for (i = 0; i < n; i++)
8945 copy_rvec(v[sort[i].ind], buf[i]);
8948 /* Copy back to the original array */
8949 for (i = 0; i < n; i++)
8951 copy_rvec(buf[i], v[i]);
8955 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8958 int a, atot, cg, cg0, cg1, i;
8960 if (cgindex == NULL)
8962 /* Avoid the useless loop of the atoms within a cg */
8963 order_vec_cg(ncg, sort, v, buf);
8968 /* Order the data */
8970 for (cg = 0; cg < ncg; cg++)
8972 cg0 = cgindex[sort[cg].ind];
8973 cg1 = cgindex[sort[cg].ind+1];
8974 for (i = cg0; i < cg1; i++)
8976 copy_rvec(v[i], buf[a]);
8982 /* Copy back to the original array */
8983 for (a = 0; a < atot; a++)
8985 copy_rvec(buf[a], v[a]);
8989 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8990 int nsort_new, gmx_cgsort_t *sort_new,
8991 gmx_cgsort_t *sort1)
8995 /* The new indices are not very ordered, so we qsort them */
8996 gmx_qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8998 /* sort2 is already ordered, so now we can merge the two arrays */
9002 while (i2 < nsort2 || i_new < nsort_new)
9006 sort1[i1++] = sort_new[i_new++];
9008 else if (i_new == nsort_new)
9010 sort1[i1++] = sort2[i2++];
9012 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
9013 (sort2[i2].nsc == sort_new[i_new].nsc &&
9014 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
9016 sort1[i1++] = sort2[i2++];
9020 sort1[i1++] = sort_new[i_new++];
9025 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
9027 gmx_domdec_sort_t *sort;
9028 gmx_cgsort_t *cgsort, *sort_i;
9029 int ncg_new, nsort2, nsort_new, i, *a, moved;
9031 sort = dd->comm->sort;
9033 a = fr->ns.grid->cell_index;
9035 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
9037 if (ncg_home_old >= 0)
9039 /* The charge groups that remained in the same ns grid cell
9040 * are completely ordered. So we can sort efficiently by sorting
9041 * the charge groups that did move into the stationary list.
9046 for (i = 0; i < dd->ncg_home; i++)
9048 /* Check if this cg did not move to another node */
9051 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
9053 /* This cg is new on this node or moved ns grid cell */
9054 if (nsort_new >= sort->sort_new_nalloc)
9056 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
9057 srenew(sort->sort_new, sort->sort_new_nalloc);
9059 sort_i = &(sort->sort_new[nsort_new++]);
9063 /* This cg did not move */
9064 sort_i = &(sort->sort2[nsort2++]);
9066 /* Sort on the ns grid cell indices
9067 * and the global topology index.
9068 * index_gl is irrelevant with cell ns,
9069 * but we set it here anyhow to avoid a conditional.
9072 sort_i->ind_gl = dd->index_gl[i];
9079 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
9082 /* Sort efficiently */
9083 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
9088 cgsort = sort->sort;
9090 for (i = 0; i < dd->ncg_home; i++)
9092 /* Sort on the ns grid cell indices
9093 * and the global topology index
9095 cgsort[i].nsc = a[i];
9096 cgsort[i].ind_gl = dd->index_gl[i];
9098 if (cgsort[i].nsc < moved)
9105 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
9107 /* Determine the order of the charge groups using qsort */
9108 gmx_qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
9114 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
9117 int ncg_new, i, *a, na;
9119 sort = dd->comm->sort->sort;
9121 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
9124 for (i = 0; i < na; i++)
9128 sort[ncg_new].ind = a[i];
9136 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
9139 gmx_domdec_sort_t *sort;
9140 gmx_cgsort_t *cgsort;
9142 int ncg_new, i, *ibuf, cgsize;
9145 sort = dd->comm->sort;
9147 if (dd->ncg_home > sort->sort_nalloc)
9149 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
9150 srenew(sort->sort, sort->sort_nalloc);
9151 srenew(sort->sort2, sort->sort_nalloc);
9153 cgsort = sort->sort;
9155 switch (fr->cutoff_scheme)
9158 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
9161 ncg_new = dd_sort_order_nbnxn(dd, fr);
9164 gmx_incons("unimplemented");
9168 /* We alloc with the old size, since cgindex is still old */
9169 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
9170 vbuf = dd->comm->vbuf.v;
9174 cgindex = dd->cgindex;
9181 /* Remove the charge groups which are no longer at home here */
9182 dd->ncg_home = ncg_new;
9185 fprintf(debug, "Set the new home charge group count to %d\n",
9189 /* Reorder the state */
9190 for (i = 0; i < estNR; i++)
9192 if (EST_DISTR(i) && (state->flags & (1<<i)))
9197 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9200 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9203 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9206 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9210 case estDISRE_INITF:
9211 case estDISRE_RM3TAV:
9212 case estORIRE_INITF:
9214 /* No ordering required */
9217 gmx_incons("Unknown state entry encountered in dd_sort_state");
9222 if (fr->cutoff_scheme == ecutsGROUP)
9225 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9228 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9230 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9231 srenew(sort->ibuf, sort->ibuf_nalloc);
9234 /* Reorder the global cg index */
9235 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9236 /* Reorder the cginfo */
9237 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9238 /* Rebuild the local cg index */
9242 for (i = 0; i < dd->ncg_home; i++)
9244 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9245 ibuf[i+1] = ibuf[i] + cgsize;
9247 for (i = 0; i < dd->ncg_home+1; i++)
9249 dd->cgindex[i] = ibuf[i];
9254 for (i = 0; i < dd->ncg_home+1; i++)
9259 /* Set the home atom number */
9260 dd->nat_home = dd->cgindex[dd->ncg_home];
9262 if (fr->cutoff_scheme == ecutsVERLET)
9264 /* The atoms are now exactly in grid order, update the grid order */
9265 nbnxn_set_atomorder(fr->nbv->nbs);
9269 /* Copy the sorted ns cell indices back to the ns grid struct */
9270 for (i = 0; i < dd->ncg_home; i++)
9272 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9274 fr->ns.grid->nr = dd->ncg_home;
9278 static void add_dd_statistics(gmx_domdec_t *dd)
9280 gmx_domdec_comm_t *comm;
9285 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9287 comm->sum_nat[ddnat-ddnatZONE] +=
9288 comm->nat[ddnat] - comm->nat[ddnat-1];
9293 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9295 gmx_domdec_comm_t *comm;
9300 /* Reset all the statistics and counters for total run counting */
9301 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9303 comm->sum_nat[ddnat-ddnatZONE] = 0;
9307 comm->load_step = 0;
9310 clear_ivec(comm->load_lim);
9315 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9317 gmx_domdec_comm_t *comm;
9321 comm = cr->dd->comm;
9323 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9330 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");
9332 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9334 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9339 " av. #atoms communicated per step for force: %d x %.1f\n",
9343 if (cr->dd->vsite_comm)
9346 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9347 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9352 if (cr->dd->constraint_comm)
9355 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9356 1 + ir->nLincsIter, av);
9360 gmx_incons(" Unknown type for DD statistics");
9363 fprintf(fplog, "\n");
9365 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9367 print_dd_load_av(fplog, cr->dd);
9371 void dd_partition_system(FILE *fplog,
9374 gmx_bool bMasterState,
9376 t_state *state_global,
9377 gmx_mtop_t *top_global,
9379 t_state *state_local,
9382 gmx_localtop_t *top_local,
9385 gmx_shellfc_t shellfc,
9386 gmx_constr_t constr,
9388 gmx_wallcycle_t wcycle,
9392 gmx_domdec_comm_t *comm;
9393 gmx_ddbox_t ddbox = {0};
9395 gmx_int64_t step_pcoupl;
9396 rvec cell_ns_x0, cell_ns_x1;
9397 int i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9398 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bTurnOnDLB, bLogLoad;
9399 gmx_bool bRedist, bSortCG, bResortAll;
9400 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9404 wallcycle_start(wcycle, ewcDOMDEC);
9409 bBoxChanged = (bMasterState || DEFORM(*ir));
9410 if (ir->epc != epcNO)
9412 /* With nstpcouple > 1 pressure coupling happens.
9413 * one step after calculating the pressure.
9414 * Box scaling happens at the end of the MD step,
9415 * after the DD partitioning.
9416 * We therefore have to do DLB in the first partitioning
9417 * after an MD step where P-coupling occured.
9418 * We need to determine the last step in which p-coupling occurred.
9419 * MRS -- need to validate this for vv?
9424 step_pcoupl = step - 1;
9428 step_pcoupl = ((step - 1)/n)*n + 1;
9430 if (step_pcoupl >= comm->partition_step)
9436 bNStGlobalComm = (step % nstglobalcomm == 0);
9438 if (!comm->bDynLoadBal)
9444 /* Should we do dynamic load balacing this step?
9445 * Since it requires (possibly expensive) global communication,
9446 * we might want to do DLB less frequently.
9448 if (bBoxChanged || ir->epc != epcNO)
9450 bDoDLB = bBoxChanged;
9454 bDoDLB = bNStGlobalComm;
9458 /* Check if we have recorded loads on the nodes */
9459 if (comm->bRecordLoad && dd_load_count(comm) > 0)
9461 bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
9463 /* Print load every nstlog, first and last step to the log file */
9464 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9465 comm->n_load_collect == 0 ||
9467 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9469 /* Avoid extra communication due to verbose screen output
9470 * when nstglobalcomm is set.
9472 if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
9473 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9475 get_load_distribution(dd, wcycle);
9480 dd_print_load(fplog, dd, step-1);
9484 dd_print_load_verbose(dd);
9487 comm->n_load_collect++;
9489 if (bCheckWhetherToTurnDlbOn)
9491 /* Since the timings are node dependent, the master decides */
9494 /* Here we check if the max PME rank load is more than 0.98
9495 * the max PP force load. If so, PP DLB will not help,
9496 * since we are (almost) limited by PME. Furthermore,
9497 * DLB will cause a significant extra x/f redistribution
9498 * cost on the PME ranks, which will then surely result
9499 * in lower total performance.
9500 * This check might be fragile, since one measurement
9501 * below 0.98 (although only done once every 100 DD part.)
9502 * could turn on DLB for the rest of the run.
9504 if (cr->npmenodes > 0 &&
9505 dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
9512 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
9516 fprintf(debug, "step %s, imb loss %f\n",
9517 gmx_step_str(step, sbuf),
9518 dd_force_imb_perf_loss(dd));
9521 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9524 turn_on_dlb(fplog, cr, step);
9529 comm->n_load_have++;
9532 cgs_gl = &comm->cgs_gl;
9537 /* Clear the old state */
9538 clear_dd_indices(dd, 0, 0);
9541 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9542 TRUE, cgs_gl, state_global->x, &ddbox);
9544 get_cg_distribution(fplog, dd, cgs_gl,
9545 state_global->box, &ddbox, state_global->x);
9547 dd_distribute_state(dd, cgs_gl,
9548 state_global, state_local, f);
9550 dd_make_local_cgs(dd, &top_local->cgs);
9552 /* Ensure that we have space for the new distribution */
9553 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9555 if (fr->cutoff_scheme == ecutsGROUP)
9557 calc_cgcm(fplog, 0, dd->ncg_home,
9558 &top_local->cgs, state_local->x, fr->cg_cm);
9561 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9563 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9565 else if (state_local->ddp_count != dd->ddp_count)
9567 if (state_local->ddp_count > dd->ddp_count)
9569 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9572 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9574 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);
9577 /* Clear the old state */
9578 clear_dd_indices(dd, 0, 0);
9580 /* Build the new indices */
9581 rebuild_cgindex(dd, cgs_gl->index, state_local);
9582 make_dd_indices(dd, cgs_gl->index, 0);
9583 ncgindex_set = dd->ncg_home;
9585 if (fr->cutoff_scheme == ecutsGROUP)
9587 /* Redetermine the cg COMs */
9588 calc_cgcm(fplog, 0, dd->ncg_home,
9589 &top_local->cgs, state_local->x, fr->cg_cm);
9592 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9594 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9596 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9597 TRUE, &top_local->cgs, state_local->x, &ddbox);
9599 bRedist = comm->bDynLoadBal;
9603 /* We have the full state, only redistribute the cgs */
9605 /* Clear the non-home indices */
9606 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9609 /* Avoid global communication for dim's without pbc and -gcom */
9610 if (!bNStGlobalComm)
9612 copy_rvec(comm->box0, ddbox.box0 );
9613 copy_rvec(comm->box_size, ddbox.box_size);
9615 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9616 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9621 /* For dim's without pbc and -gcom */
9622 copy_rvec(ddbox.box0, comm->box0 );
9623 copy_rvec(ddbox.box_size, comm->box_size);
9625 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9628 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9630 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9633 /* Check if we should sort the charge groups */
9634 if (comm->nstSortCG > 0)
9636 bSortCG = (bMasterState ||
9637 (bRedist && (step % comm->nstSortCG == 0)));
9644 ncg_home_old = dd->ncg_home;
9649 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9651 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9653 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9655 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9658 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9660 &comm->cell_x0, &comm->cell_x1,
9661 dd->ncg_home, fr->cg_cm,
9662 cell_ns_x0, cell_ns_x1, &grid_density);
9666 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9669 switch (fr->cutoff_scheme)
9672 copy_ivec(fr->ns.grid->n, ncells_old);
9673 grid_first(fplog, fr->ns.grid, dd, &ddbox,
9674 state_local->box, cell_ns_x0, cell_ns_x1,
9675 fr->rlistlong, grid_density);
9678 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9681 gmx_incons("unimplemented");
9683 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9684 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9688 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9690 /* Sort the state on charge group position.
9691 * This enables exact restarts from this step.
9692 * It also improves performance by about 15% with larger numbers
9693 * of atoms per node.
9696 /* Fill the ns grid with the home cell,
9697 * so we can sort with the indices.
9699 set_zones_ncg_home(dd);
9701 switch (fr->cutoff_scheme)
9704 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9706 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9708 comm->zones.size[0].bb_x0,
9709 comm->zones.size[0].bb_x1,
9711 comm->zones.dens_zone0,
9714 ncg_moved, bRedist ? comm->moved : NULL,
9715 fr->nbv->grp[eintLocal].kernel_type,
9716 fr->nbv->grp[eintLocal].nbat);
9718 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9721 fill_grid(&comm->zones, fr->ns.grid, dd->ncg_home,
9722 0, dd->ncg_home, fr->cg_cm);
9724 copy_ivec(fr->ns.grid->n, ncells_new);
9727 gmx_incons("unimplemented");
9730 bResortAll = bMasterState;
9732 /* Check if we can user the old order and ns grid cell indices
9733 * of the charge groups to sort the charge groups efficiently.
9735 if (ncells_new[XX] != ncells_old[XX] ||
9736 ncells_new[YY] != ncells_old[YY] ||
9737 ncells_new[ZZ] != ncells_old[ZZ])
9744 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9745 gmx_step_str(step, sbuf), dd->ncg_home);
9747 dd_sort_state(dd, fr->cg_cm, fr, state_local,
9748 bResortAll ? -1 : ncg_home_old);
9749 /* Rebuild all the indices */
9750 ga2la_clear(dd->ga2la);
9753 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9756 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9758 /* Setup up the communication and communicate the coordinates */
9759 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9761 /* Set the indices */
9762 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9764 /* Set the charge group boundaries for neighbor searching */
9765 set_cg_boundaries(&comm->zones);
9767 if (fr->cutoff_scheme == ecutsVERLET)
9769 set_zones_size(dd, state_local->box, &ddbox,
9770 bSortCG ? 1 : 0, comm->zones.n);
9773 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9776 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9777 -1,state_local->x,state_local->box);
9780 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9782 /* Extract a local topology from the global topology */
9783 for (i = 0; i < dd->ndim; i++)
9785 np[dd->dim[i]] = comm->cd[i].np;
9787 dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
9788 comm->cellsize_min, np,
9790 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9791 vsite, top_global, top_local);
9793 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9795 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9797 /* Set up the special atom communication */
9798 n = comm->nat[ddnatZONE];
9799 for (i = ddnatZONE+1; i < ddnatNR; i++)
9804 if (vsite && vsite->n_intercg_vsite)
9806 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9810 if (dd->bInterCGcons || dd->bInterCGsettles)
9812 /* Only for inter-cg constraints we need special code */
9813 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9814 constr, ir->nProjOrder,
9815 top_local->idef.il);
9819 gmx_incons("Unknown special atom type setup");
9824 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9826 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9828 /* Make space for the extra coordinates for virtual site
9829 * or constraint communication.
9831 state_local->natoms = comm->nat[ddnatNR-1];
9832 if (state_local->natoms > state_local->nalloc)
9834 dd_realloc_state(state_local, f, state_local->natoms);
9837 if (fr->bF_NoVirSum)
9839 if (vsite && vsite->n_intercg_vsite)
9841 nat_f_novirsum = comm->nat[ddnatVSITE];
9845 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9847 nat_f_novirsum = dd->nat_tot;
9851 nat_f_novirsum = dd->nat_home;
9860 /* Set the number of atoms required for the force calculation.
9861 * Forces need to be constrained when using a twin-range setup
9862 * or with energy minimization. For simple simulations we could
9863 * avoid some allocation, zeroing and copying, but this is
9864 * probably not worth the complications ande checking.
9866 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9867 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9869 /* We make the all mdatoms up to nat_tot_con.
9870 * We could save some work by only setting invmass
9871 * between nat_tot and nat_tot_con.
9873 /* This call also sets the new number of home particles to dd->nat_home */
9874 atoms2md(top_global, ir,
9875 comm->nat[ddnatCON], dd->gatindex, dd->nat_home, mdatoms);
9877 /* Now we have the charges we can sort the FE interactions */
9878 dd_sort_local_top(dd, mdatoms, top_local);
9882 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9883 split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
9884 mdatoms, FALSE, vsite);
9889 /* Make the local shell stuff, currently no communication is done */
9890 make_local_shells(cr, mdatoms, shellfc);
9893 if (ir->implicit_solvent)
9895 make_local_gb(cr, fr->born, ir->gb_algorithm);
9898 setup_bonded_threading(fr, &top_local->idef);
9900 if (!(cr->duty & DUTY_PME))
9902 /* Send the charges and/or c6/sigmas to our PME only node */
9903 gmx_pme_send_parameters(cr,
9905 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
9906 mdatoms->chargeA, mdatoms->chargeB,
9907 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
9908 mdatoms->sigmaA, mdatoms->sigmaB,
9909 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9914 set_constraints(constr, top_local, ir, mdatoms, cr);
9919 /* Update the local pull groups */
9920 dd_make_local_pull_groups(cr, ir->pull_work, mdatoms);
9925 /* Update the local rotation groups */
9926 dd_make_local_rotation_groups(dd, ir->rot);
9929 if (ir->eSwapCoords != eswapNO)
9931 /* Update the local groups needed for ion swapping */
9932 dd_make_local_swap_groups(dd, ir->swap);
9935 /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
9936 dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
9938 add_dd_statistics(dd);
9940 /* Make sure we only count the cycles for this DD partitioning */
9941 clear_dd_cycle_counts(dd);
9943 /* Because the order of the atoms might have changed since
9944 * the last vsite construction, we need to communicate the constructing
9945 * atom coordinates again (for spreading the forces this MD step).
9947 dd_move_x_vsites(dd, state_local->box, state_local->x);
9949 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9951 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9953 dd_move_x(dd, state_local->box, state_local->x);
9954 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9955 -1, state_local->x, state_local->box);
9958 /* Store the partitioning step */
9959 comm->partition_step = step;
9961 /* Increase the DD partitioning counter */
9963 /* The state currently matches this DD partitioning count, store it */
9964 state_local->ddp_count = dd->ddp_count;
9967 /* The DD master node knows the complete cg distribution,
9968 * store the count so we can possibly skip the cg info communication.
9970 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9973 if (comm->DD_debug > 0)
9975 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9976 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9977 "after partitioning");
9980 wallcycle_stop(wcycle, ewcDOMDEC);