2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2008
5 * Copyright (c) 2012,2013, by the GROMACS development team, led by
6 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
7 * others, as listed in the AUTHORS file in the top-level source
8 * directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gmx_fatal.h"
49 #include "gmx_fatal_collective.h"
52 #include "domdec_network.h"
55 #include "chargegroup.h"
64 #include "pull_rotation.h"
65 #include "gmx_wallcycle.h"
69 #include "mtop_util.h"
71 #include "gmx_ga2la.h"
73 #include "nbnxn_search.h"
75 #include "gmx_omp_nthreads.h"
84 #define DDRANK(dd, rank) (rank)
85 #define DDMASTERRANK(dd) (dd->masterrank)
87 typedef struct gmx_domdec_master
89 /* The cell boundaries */
91 /* The global charge group division */
92 int *ncg; /* Number of home charge groups for each node */
93 int *index; /* Index of nnodes+1 into cg */
94 int *cg; /* Global charge group index */
95 int *nat; /* Number of home atoms for each node. */
96 int *ibuf; /* Buffer for communication */
97 rvec *vbuf; /* Buffer for state scattering and gathering */
98 } gmx_domdec_master_t;
102 /* The numbers of charge groups to send and receive for each cell
103 * that requires communication, the last entry contains the total
104 * number of atoms that needs to be communicated.
106 int nsend[DD_MAXIZONE+2];
107 int nrecv[DD_MAXIZONE+2];
108 /* The charge groups to send */
111 /* The atom range for non-in-place communication */
112 int cell2at0[DD_MAXIZONE];
113 int cell2at1[DD_MAXIZONE];
118 int np; /* Number of grid pulses in this dimension */
119 int np_dlb; /* For dlb, for use with edlbAUTO */
120 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
122 gmx_bool bInPlace; /* Can we communicate in place? */
123 } gmx_domdec_comm_dim_t;
127 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
128 real *cell_f; /* State var.: cell boundaries, box relative */
129 real *old_cell_f; /* Temp. var.: old cell size */
130 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
131 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
132 real *bound_min; /* Temp. var.: lower limit for cell boundary */
133 real *bound_max; /* Temp. var.: upper limit for cell boundary */
134 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
135 real *buf_ncd; /* Temp. var. */
138 #define DD_NLOAD_MAX 9
140 /* Here floats are accurate enough, since these variables
141 * only influence the load balancing, not the actual MD results.
168 gmx_cgsort_t *sort_new;
180 /* This enum determines the order of the coordinates.
181 * ddnatHOME and ddnatZONE should be first and second,
182 * the others can be ordered as wanted.
185 ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR
189 edlbAUTO, edlbNO, edlbYES, edlbNR
191 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
195 int dim; /* The dimension */
196 gmx_bool dim_match; /* Tells if DD and PME dims match */
197 int nslab; /* The number of PME slabs in this dimension */
198 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
199 int *pp_min; /* The minimum pp node location, size nslab */
200 int *pp_max; /* The maximum pp node location,size nslab */
201 int maxshift; /* The maximum shift for coordinate redistribution in PME */
206 real min0; /* The minimum bottom of this zone */
207 real max1; /* The maximum top of this zone */
208 real min1; /* The minimum top of this zone */
209 real mch0; /* The maximum bottom communicaton height for this zone */
210 real mch1; /* The maximum top communicaton height for this zone */
211 real p1_0; /* The bottom value of the first cell in this zone */
212 real p1_1; /* The top value of the first cell in this zone */
217 gmx_domdec_ind_t ind;
224 } dd_comm_setup_work_t;
226 typedef struct gmx_domdec_comm
228 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
229 * unless stated otherwise.
232 /* The number of decomposition dimensions for PME, 0: no PME */
234 /* The number of nodes doing PME (PP/PME or only PME) */
238 /* The communication setup including the PME only nodes */
239 gmx_bool bCartesianPP_PME;
242 int *pmenodes; /* size npmenodes */
243 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
244 * but with bCartesianPP_PME */
245 gmx_ddpme_t ddpme[2];
247 /* The DD particle-particle nodes only */
248 gmx_bool bCartesianPP;
249 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
251 /* The global charge groups */
254 /* Should we sort the cgs */
256 gmx_domdec_sort_t *sort;
258 /* Are there charge groups? */
261 /* Are there bonded and multi-body interactions between charge groups? */
262 gmx_bool bInterCGBondeds;
263 gmx_bool bInterCGMultiBody;
265 /* Data for the optional bonded interaction atom communication range */
272 /* Are we actually using DLB? */
273 gmx_bool bDynLoadBal;
275 /* Cell sizes for static load balancing, first index cartesian */
278 /* The width of the communicated boundaries */
281 /* The minimum cell size (including triclinic correction) */
283 /* For dlb, for use with edlbAUTO */
284 rvec cellsize_min_dlb;
285 /* The lower limit for the DD cell size with DLB */
287 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
288 gmx_bool bVacDLBNoLimit;
290 /* With PME load balancing we set limits on DLB */
291 gmx_bool bPMELoadBalDLBLimits;
292 /* DLB needs to take into account that we want to allow this maximum
293 * cut-off (for PME load balancing), this could limit cell boundaries.
295 real PMELoadBal_max_cutoff;
297 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
299 /* box0 and box_size are required with dim's without pbc and -gcom */
303 /* The cell boundaries */
307 /* The old location of the cell boundaries, to check cg displacements */
311 /* The communication setup and charge group boundaries for the zones */
312 gmx_domdec_zones_t zones;
314 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
315 * cell boundaries of neighboring cells for dynamic load balancing.
317 gmx_ddzone_t zone_d1[2];
318 gmx_ddzone_t zone_d2[2][2];
320 /* The coordinate/force communication setup and indices */
321 gmx_domdec_comm_dim_t cd[DIM];
322 /* The maximum number of cells to communicate with in one dimension */
325 /* Which cg distribution is stored on the master node */
326 int master_cg_ddp_count;
328 /* The number of cg's received from the direct neighbors */
329 int zone_ncg1[DD_MAXZONE];
331 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
334 /* Array for signalling if atoms have moved to another domain */
338 /* Communication buffer for general use */
342 /* Communication buffer for general use */
345 /* Temporary storage for thread parallel communication setup */
347 dd_comm_setup_work_t *dth;
349 /* Communication buffers only used with multiple grid pulses */
354 /* Communication buffers for local redistribution */
356 int cggl_flag_nalloc[DIM*2];
358 int cgcm_state_nalloc[DIM*2];
360 /* Cell sizes for dynamic load balancing */
361 gmx_domdec_root_t **root;
365 real cell_f_max0[DIM];
366 real cell_f_min1[DIM];
368 /* Stuff for load communication */
369 gmx_bool bRecordLoad;
370 gmx_domdec_load_t *load;
372 MPI_Comm *mpi_comm_load;
375 /* Maximum DLB scaling per load balancing step in percent */
379 float cycl[ddCyclNr];
380 int cycl_n[ddCyclNr];
381 float cycl_max[ddCyclNr];
382 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
386 /* Have often have did we have load measurements */
388 /* Have often have we collected the load measurements */
392 double sum_nat[ddnatNR-ddnatZONE];
402 /* The last partition step */
403 gmx_large_int_t partition_step;
411 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
414 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
415 #define DD_FLAG_NRCG 65535
416 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
417 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
419 /* Zone permutation required to obtain consecutive charge groups
420 * for neighbor searching.
422 static const int zone_perm[3][4] = { {0, 0, 0, 0}, {1, 0, 0, 0}, {3, 0, 1, 2} };
424 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
425 * components see only j zones with that component 0.
428 /* The DD zone order */
429 static const ivec dd_zo[DD_MAXZONE] =
430 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
435 static const ivec dd_zp3[dd_zp3n] = {{0, 0, 8}, {1, 3, 6}, {2, 5, 6}, {3, 5, 7}};
440 static const ivec dd_zp2[dd_zp2n] = {{0, 0, 4}, {1, 3, 4}};
445 static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
447 /* Factors used to avoid problems due to rounding issues */
448 #define DD_CELL_MARGIN 1.0001
449 #define DD_CELL_MARGIN2 1.00005
450 /* Factor to account for pressure scaling during nstlist steps */
451 #define DD_PRES_SCALE_MARGIN 1.02
453 /* Allowed performance loss before we DLB or warn */
454 #define DD_PERF_LOSS 0.05
456 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
458 /* Use separate MPI send and receive commands
459 * when nnodes <= GMX_DD_NNODES_SENDRECV.
460 * This saves memory (and some copying for small nnodes).
461 * For high parallelization scatter and gather calls are used.
463 #define GMX_DD_NNODES_SENDRECV 4
467 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
469 static void index2xyz(ivec nc,int ind,ivec xyz)
471 xyz[XX] = ind % nc[XX];
472 xyz[YY] = (ind / nc[XX]) % nc[YY];
473 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
477 /* This order is required to minimize the coordinate communication in PME
478 * which uses decomposition in the x direction.
480 #define dd_index(n, i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
482 static void ddindex2xyz(ivec nc, int ind, ivec xyz)
484 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
485 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
486 xyz[ZZ] = ind % nc[ZZ];
489 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
494 ddindex = dd_index(dd->nc, c);
495 if (dd->comm->bCartesianPP_PME)
497 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
499 else if (dd->comm->bCartesianPP)
502 MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
513 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox, t_inputrec *ir)
515 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
518 int ddglatnr(gmx_domdec_t *dd, int i)
528 if (i >= dd->comm->nat[ddnatNR-1])
530 gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
532 atnr = dd->gatindex[i] + 1;
538 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
540 return &dd->comm->cgs_gl;
543 static void vec_rvec_init(vec_rvec_t *v)
549 static void vec_rvec_check_alloc(vec_rvec_t *v, int n)
553 v->nalloc = over_alloc_dd(n);
554 srenew(v->v, v->nalloc);
558 void dd_store_state(gmx_domdec_t *dd, t_state *state)
562 if (state->ddp_count != dd->ddp_count)
564 gmx_incons("The state does not the domain decomposition state");
567 state->ncg_gl = dd->ncg_home;
568 if (state->ncg_gl > state->cg_gl_nalloc)
570 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
571 srenew(state->cg_gl, state->cg_gl_nalloc);
573 for (i = 0; i < state->ncg_gl; i++)
575 state->cg_gl[i] = dd->index_gl[i];
578 state->ddp_count_cg_gl = dd->ddp_count;
581 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
583 return &dd->comm->zones;
586 void dd_get_ns_ranges(gmx_domdec_t *dd, int icg,
587 int *jcg0, int *jcg1, ivec shift0, ivec shift1)
589 gmx_domdec_zones_t *zones;
592 zones = &dd->comm->zones;
595 while (icg >= zones->izone[izone].cg1)
604 else if (izone < zones->nizone)
606 *jcg0 = zones->izone[izone].jcg0;
610 gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
611 icg, izone, zones->nizone);
614 *jcg1 = zones->izone[izone].jcg1;
616 for (d = 0; d < dd->ndim; d++)
619 shift0[dim] = zones->izone[izone].shift0[dim];
620 shift1[dim] = zones->izone[izone].shift1[dim];
621 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
623 /* A conservative approach, this can be optimized */
630 int dd_natoms_vsite(gmx_domdec_t *dd)
632 return dd->comm->nat[ddnatVSITE];
635 void dd_get_constraint_range(gmx_domdec_t *dd, int *at_start, int *at_end)
637 *at_start = dd->comm->nat[ddnatCON-1];
638 *at_end = dd->comm->nat[ddnatCON];
641 void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[])
643 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
644 int *index, *cgindex;
645 gmx_domdec_comm_t *comm;
646 gmx_domdec_comm_dim_t *cd;
647 gmx_domdec_ind_t *ind;
648 rvec shift = {0, 0, 0}, *buf, *rbuf;
649 gmx_bool bPBC, bScrew;
653 cgindex = dd->cgindex;
658 nat_tot = dd->nat_home;
659 for (d = 0; d < dd->ndim; d++)
661 bPBC = (dd->ci[dd->dim[d]] == 0);
662 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
665 copy_rvec(box[dd->dim[d]], shift);
668 for (p = 0; p < cd->np; p++)
675 for (i = 0; i < ind->nsend[nzone]; i++)
677 at0 = cgindex[index[i]];
678 at1 = cgindex[index[i]+1];
679 for (j = at0; j < at1; j++)
681 copy_rvec(x[j], buf[n]);
688 for (i = 0; i < ind->nsend[nzone]; i++)
690 at0 = cgindex[index[i]];
691 at1 = cgindex[index[i]+1];
692 for (j = at0; j < at1; j++)
694 /* We need to shift the coordinates */
695 rvec_add(x[j], shift, buf[n]);
702 for (i = 0; i < ind->nsend[nzone]; i++)
704 at0 = cgindex[index[i]];
705 at1 = cgindex[index[i]+1];
706 for (j = at0; j < at1; j++)
709 buf[n][XX] = x[j][XX] + shift[XX];
711 * This operation requires a special shift force
712 * treatment, which is performed in calc_vir.
714 buf[n][YY] = box[YY][YY] - x[j][YY];
715 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
727 rbuf = comm->vbuf2.v;
729 /* Send and receive the coordinates */
730 dd_sendrecv_rvec(dd, d, dddirBackward,
731 buf, ind->nsend[nzone+1],
732 rbuf, ind->nrecv[nzone+1]);
736 for (zone = 0; zone < nzone; zone++)
738 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
740 copy_rvec(rbuf[j], x[i]);
745 nat_tot += ind->nrecv[nzone+1];
751 void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
753 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
754 int *index, *cgindex;
755 gmx_domdec_comm_t *comm;
756 gmx_domdec_comm_dim_t *cd;
757 gmx_domdec_ind_t *ind;
761 gmx_bool bPBC, bScrew;
765 cgindex = dd->cgindex;
770 nzone = comm->zones.n/2;
771 nat_tot = dd->nat_tot;
772 for (d = dd->ndim-1; d >= 0; d--)
774 bPBC = (dd->ci[dd->dim[d]] == 0);
775 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
776 if (fshift == NULL && !bScrew)
780 /* Determine which shift vector we need */
786 for (p = cd->np-1; p >= 0; p--)
789 nat_tot -= ind->nrecv[nzone+1];
796 sbuf = comm->vbuf2.v;
798 for (zone = 0; zone < nzone; zone++)
800 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
802 copy_rvec(f[i], sbuf[j]);
807 /* Communicate the forces */
808 dd_sendrecv_rvec(dd, d, dddirForward,
809 sbuf, ind->nrecv[nzone+1],
810 buf, ind->nsend[nzone+1]);
812 /* Add the received forces */
816 for (i = 0; i < ind->nsend[nzone]; i++)
818 at0 = cgindex[index[i]];
819 at1 = cgindex[index[i]+1];
820 for (j = at0; j < at1; j++)
822 rvec_inc(f[j], buf[n]);
829 for (i = 0; i < ind->nsend[nzone]; i++)
831 at0 = cgindex[index[i]];
832 at1 = cgindex[index[i]+1];
833 for (j = at0; j < at1; j++)
835 rvec_inc(f[j], buf[n]);
836 /* Add this force to the shift force */
837 rvec_inc(fshift[is], buf[n]);
844 for (i = 0; i < ind->nsend[nzone]; i++)
846 at0 = cgindex[index[i]];
847 at1 = cgindex[index[i]+1];
848 for (j = at0; j < at1; j++)
850 /* Rotate the force */
851 f[j][XX] += buf[n][XX];
852 f[j][YY] -= buf[n][YY];
853 f[j][ZZ] -= buf[n][ZZ];
856 /* Add this force to the shift force */
857 rvec_inc(fshift[is], buf[n]);
868 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
870 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
871 int *index, *cgindex;
872 gmx_domdec_comm_t *comm;
873 gmx_domdec_comm_dim_t *cd;
874 gmx_domdec_ind_t *ind;
879 cgindex = dd->cgindex;
881 buf = &comm->vbuf.v[0][0];
884 nat_tot = dd->nat_home;
885 for (d = 0; d < dd->ndim; d++)
888 for (p = 0; p < cd->np; p++)
893 for (i = 0; i < ind->nsend[nzone]; i++)
895 at0 = cgindex[index[i]];
896 at1 = cgindex[index[i]+1];
897 for (j = at0; j < at1; j++)
910 rbuf = &comm->vbuf2.v[0][0];
912 /* Send and receive the coordinates */
913 dd_sendrecv_real(dd, d, dddirBackward,
914 buf, ind->nsend[nzone+1],
915 rbuf, ind->nrecv[nzone+1]);
919 for (zone = 0; zone < nzone; zone++)
921 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
928 nat_tot += ind->nrecv[nzone+1];
934 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
936 int nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
937 int *index, *cgindex;
938 gmx_domdec_comm_t *comm;
939 gmx_domdec_comm_dim_t *cd;
940 gmx_domdec_ind_t *ind;
945 cgindex = dd->cgindex;
947 buf = &comm->vbuf.v[0][0];
950 nzone = comm->zones.n/2;
951 nat_tot = dd->nat_tot;
952 for (d = dd->ndim-1; d >= 0; d--)
955 for (p = cd->np-1; p >= 0; p--)
958 nat_tot -= ind->nrecv[nzone+1];
965 sbuf = &comm->vbuf2.v[0][0];
967 for (zone = 0; zone < nzone; zone++)
969 for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
976 /* Communicate the forces */
977 dd_sendrecv_real(dd, d, dddirForward,
978 sbuf, ind->nrecv[nzone+1],
979 buf, ind->nsend[nzone+1]);
981 /* Add the received forces */
983 for (i = 0; i < ind->nsend[nzone]; i++)
985 at0 = cgindex[index[i]];
986 at1 = cgindex[index[i]+1];
987 for (j = at0; j < at1; j++)
998 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
1000 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",
1002 zone->min0, zone->max1,
1003 zone->mch0, zone->mch0,
1004 zone->p1_0, zone->p1_1);
1008 #define DDZONECOMM_MAXZONE 5
1009 #define DDZONECOMM_BUFSIZE 3
1011 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
1012 int ddimind, int direction,
1013 gmx_ddzone_t *buf_s, int n_s,
1014 gmx_ddzone_t *buf_r, int n_r)
1016 #define ZBS DDZONECOMM_BUFSIZE
1017 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
1018 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1021 for (i = 0; i < n_s; i++)
1023 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1024 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1025 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1026 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1027 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1028 vbuf_s[i*ZBS+1][2] = 0;
1029 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1030 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1031 vbuf_s[i*ZBS+2][2] = 0;
1034 dd_sendrecv_rvec(dd, ddimind, direction,
1038 for (i = 0; i < n_r; i++)
1040 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1041 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1042 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1043 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1044 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1045 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1046 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1052 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
1053 rvec cell_ns_x0, rvec cell_ns_x1)
1055 int d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
1057 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1058 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1059 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1060 rvec extr_s[2], extr_r[2];
1062 real dist_d, c = 0, det;
1063 gmx_domdec_comm_t *comm;
1064 gmx_bool bPBC, bUse;
1068 for (d = 1; d < dd->ndim; d++)
1071 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1072 zp->min0 = cell_ns_x0[dim];
1073 zp->max1 = cell_ns_x1[dim];
1074 zp->min1 = cell_ns_x1[dim];
1075 zp->mch0 = cell_ns_x0[dim];
1076 zp->mch1 = cell_ns_x1[dim];
1077 zp->p1_0 = cell_ns_x0[dim];
1078 zp->p1_1 = cell_ns_x1[dim];
1081 for (d = dd->ndim-2; d >= 0; d--)
1084 bPBC = (dim < ddbox->npbcdim);
1086 /* Use an rvec to store two reals */
1087 extr_s[d][0] = comm->cell_f0[d+1];
1088 extr_s[d][1] = comm->cell_f1[d+1];
1089 extr_s[d][2] = comm->cell_f1[d+1];
1092 /* Store the extremes in the backward sending buffer,
1093 * so the get updated separately from the forward communication.
1095 for (d1 = d; d1 < dd->ndim-1; d1++)
1097 /* We invert the order to be able to use the same loop for buf_e */
1098 buf_s[pos].min0 = extr_s[d1][1];
1099 buf_s[pos].max1 = extr_s[d1][0];
1100 buf_s[pos].min1 = extr_s[d1][2];
1101 buf_s[pos].mch0 = 0;
1102 buf_s[pos].mch1 = 0;
1103 /* Store the cell corner of the dimension we communicate along */
1104 buf_s[pos].p1_0 = comm->cell_x0[dim];
1105 buf_s[pos].p1_1 = 0;
1109 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1112 if (dd->ndim == 3 && d == 0)
1114 buf_s[pos] = comm->zone_d2[0][1];
1116 buf_s[pos] = comm->zone_d1[0];
1120 /* We only need to communicate the extremes
1121 * in the forward direction
1123 npulse = comm->cd[d].np;
1126 /* Take the minimum to avoid double communication */
1127 npulse_min = min(npulse, dd->nc[dim]-1-npulse);
1131 /* Without PBC we should really not communicate over
1132 * the boundaries, but implementing that complicates
1133 * the communication setup and therefore we simply
1134 * do all communication, but ignore some data.
1136 npulse_min = npulse;
1138 for (p = 0; p < npulse_min; p++)
1140 /* Communicate the extremes forward */
1141 bUse = (bPBC || dd->ci[dim] > 0);
1143 dd_sendrecv_rvec(dd, d, dddirForward,
1144 extr_s+d, dd->ndim-d-1,
1145 extr_r+d, dd->ndim-d-1);
1149 for (d1 = d; d1 < dd->ndim-1; d1++)
1151 extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
1152 extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
1153 extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
1159 for (p = 0; p < npulse; p++)
1161 /* Communicate all the zone information backward */
1162 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1164 dd_sendrecv_ddzone(dd, d, dddirBackward,
1171 for (d1 = d+1; d1 < dd->ndim; d1++)
1173 /* Determine the decrease of maximum required
1174 * communication height along d1 due to the distance along d,
1175 * this avoids a lot of useless atom communication.
1177 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1179 if (ddbox->tric_dir[dim])
1181 /* c is the off-diagonal coupling between the cell planes
1182 * along directions d and d1.
1184 c = ddbox->v[dim][dd->dim[d1]][dim];
1190 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1193 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1197 /* A negative value signals out of range */
1203 /* Accumulate the extremes over all pulses */
1204 for (i = 0; i < buf_size; i++)
1208 buf_e[i] = buf_r[i];
1214 buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
1215 buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
1216 buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
1219 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1227 if (bUse && dh[d1] >= 0)
1229 buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
1230 buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
1233 /* Copy the received buffer to the send buffer,
1234 * to pass the data through with the next pulse.
1236 buf_s[i] = buf_r[i];
1238 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1239 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1241 /* Store the extremes */
1244 for (d1 = d; d1 < dd->ndim-1; d1++)
1246 extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
1247 extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
1248 extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
1252 if (d == 1 || (d == 0 && dd->ndim == 3))
1254 for (i = d; i < 2; i++)
1256 comm->zone_d2[1-d][i] = buf_e[pos];
1262 comm->zone_d1[1] = buf_e[pos];
1272 for (i = 0; i < 2; i++)
1276 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
1278 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
1279 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
1285 for (i = 0; i < 2; i++)
1287 for (j = 0; j < 2; j++)
1291 print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
1293 cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
1294 cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
1298 for (d = 1; d < dd->ndim; d++)
1300 comm->cell_f_max0[d] = extr_s[d-1][0];
1301 comm->cell_f_min1[d] = extr_s[d-1][1];
1304 fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1305 d, comm->cell_f_max0[d], comm->cell_f_min1[d]);
1310 static void dd_collect_cg(gmx_domdec_t *dd,
1311 t_state *state_local)
1313 gmx_domdec_master_t *ma = NULL;
1314 int buf2[2], *ibuf, i, ncg_home = 0, *cg = NULL, nat_home = 0;
1317 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1319 /* The master has the correct distribution */
1323 if (state_local->ddp_count == dd->ddp_count)
1325 ncg_home = dd->ncg_home;
1327 nat_home = dd->nat_home;
1329 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1331 cgs_gl = &dd->comm->cgs_gl;
1333 ncg_home = state_local->ncg_gl;
1334 cg = state_local->cg_gl;
1336 for (i = 0; i < ncg_home; i++)
1338 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1343 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1346 buf2[0] = dd->ncg_home;
1347 buf2[1] = dd->nat_home;
1357 /* Collect the charge group and atom counts on the master */
1358 dd_gather(dd, 2*sizeof(int), buf2, ibuf);
1363 for (i = 0; i < dd->nnodes; i++)
1365 ma->ncg[i] = ma->ibuf[2*i];
1366 ma->nat[i] = ma->ibuf[2*i+1];
1367 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1370 /* Make byte counts and indices */
1371 for (i = 0; i < dd->nnodes; i++)
1373 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1374 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1378 fprintf(debug, "Initial charge group distribution: ");
1379 for (i = 0; i < dd->nnodes; i++)
1381 fprintf(debug, " %d", ma->ncg[i]);
1383 fprintf(debug, "\n");
1387 /* Collect the charge group indices on the master */
1389 dd->ncg_home*sizeof(int), dd->index_gl,
1390 DDMASTER(dd) ? ma->ibuf : NULL,
1391 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1392 DDMASTER(dd) ? ma->cg : NULL);
1394 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1397 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1400 gmx_domdec_master_t *ma;
1401 int n, i, c, a, nalloc = 0;
1410 MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1411 dd->rank, dd->mpi_comm_all);
1416 /* Copy the master coordinates to the global array */
1417 cgs_gl = &dd->comm->cgs_gl;
1419 n = DDMASTERRANK(dd);
1421 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1423 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1425 copy_rvec(lv[a++], v[c]);
1429 for (n = 0; n < dd->nnodes; n++)
1433 if (ma->nat[n] > nalloc)
1435 nalloc = over_alloc_dd(ma->nat[n]);
1436 srenew(buf, nalloc);
1439 MPI_Recv(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE, DDRANK(dd, n),
1440 n, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1443 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1445 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1447 copy_rvec(buf[a++], v[c]);
1456 static void get_commbuffer_counts(gmx_domdec_t *dd,
1457 int **counts, int **disps)
1459 gmx_domdec_master_t *ma;
1464 /* Make the rvec count and displacment arrays */
1466 *disps = ma->ibuf + dd->nnodes;
1467 for (n = 0; n < dd->nnodes; n++)
1469 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1470 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1474 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1477 gmx_domdec_master_t *ma;
1478 int *rcounts = NULL, *disps = NULL;
1487 get_commbuffer_counts(dd, &rcounts, &disps);
1492 dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv, rcounts, disps, buf);
1496 cgs_gl = &dd->comm->cgs_gl;
1499 for (n = 0; n < dd->nnodes; n++)
1501 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1503 for (c = cgs_gl->index[ma->cg[i]]; c < cgs_gl->index[ma->cg[i]+1]; c++)
1505 copy_rvec(buf[a++], v[c]);
1512 void dd_collect_vec(gmx_domdec_t *dd,
1513 t_state *state_local, rvec *lv, rvec *v)
1515 gmx_domdec_master_t *ma;
1516 int n, i, c, a, nalloc = 0;
1519 dd_collect_cg(dd, state_local);
1521 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1523 dd_collect_vec_sendrecv(dd, lv, v);
1527 dd_collect_vec_gatherv(dd, lv, v);
1532 void dd_collect_state(gmx_domdec_t *dd,
1533 t_state *state_local, t_state *state)
1537 nh = state->nhchainlength;
1541 for (i = 0; i < efptNR; i++)
1543 state->lambda[i] = state_local->lambda[i];
1545 state->fep_state = state_local->fep_state;
1546 state->veta = state_local->veta;
1547 state->vol0 = state_local->vol0;
1548 copy_mat(state_local->box, state->box);
1549 copy_mat(state_local->boxv, state->boxv);
1550 copy_mat(state_local->svir_prev, state->svir_prev);
1551 copy_mat(state_local->fvir_prev, state->fvir_prev);
1552 copy_mat(state_local->pres_prev, state->pres_prev);
1554 for (i = 0; i < state_local->ngtc; i++)
1556 for (j = 0; j < nh; j++)
1558 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1559 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1561 state->therm_integral[i] = state_local->therm_integral[i];
1563 for (i = 0; i < state_local->nnhpres; i++)
1565 for (j = 0; j < nh; j++)
1567 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1568 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1572 for (est = 0; est < estNR; est++)
1574 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1579 dd_collect_vec(dd, state_local, state_local->x, state->x);
1582 dd_collect_vec(dd, state_local, state_local->v, state->v);
1585 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1588 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1591 if (state->nrngi == 1)
1595 for (i = 0; i < state_local->nrng; i++)
1597 state->ld_rng[i] = state_local->ld_rng[i];
1603 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1604 state_local->ld_rng, state->ld_rng);
1608 if (state->nrngi == 1)
1612 state->ld_rngi[0] = state_local->ld_rngi[0];
1617 dd_gather(dd, sizeof(state->ld_rngi[0]),
1618 state_local->ld_rngi, state->ld_rngi);
1621 case estDISRE_INITF:
1622 case estDISRE_RM3TAV:
1623 case estORIRE_INITF:
1627 gmx_incons("Unknown state entry encountered in dd_collect_state");
1633 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1639 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1642 state->nalloc = over_alloc_dd(nalloc);
1644 for (est = 0; est < estNR; est++)
1646 if (EST_DISTR(est) && (state->flags & (1<<est)))
1651 srenew(state->x, state->nalloc);
1654 srenew(state->v, state->nalloc);
1657 srenew(state->sd_X, state->nalloc);
1660 srenew(state->cg_p, state->nalloc);
1664 case estDISRE_INITF:
1665 case estDISRE_RM3TAV:
1666 case estORIRE_INITF:
1668 /* No reallocation required */
1671 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1678 srenew(*f, state->nalloc);
1682 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1685 if (nalloc > fr->cg_nalloc)
1689 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1691 fr->cg_nalloc = over_alloc_dd(nalloc);
1692 srenew(fr->cginfo, fr->cg_nalloc);
1693 if (fr->cutoff_scheme == ecutsGROUP)
1695 srenew(fr->cg_cm, fr->cg_nalloc);
1698 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1700 /* We don't use charge groups, we use x in state to set up
1701 * the atom communication.
1703 dd_realloc_state(state, f, nalloc);
1707 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1710 gmx_domdec_master_t *ma;
1711 int n, i, c, a, nalloc = 0;
1718 for (n = 0; n < dd->nnodes; n++)
1722 if (ma->nat[n] > nalloc)
1724 nalloc = over_alloc_dd(ma->nat[n]);
1725 srenew(buf, nalloc);
1727 /* Use lv as a temporary buffer */
1729 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1731 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1733 copy_rvec(v[c], buf[a++]);
1736 if (a != ma->nat[n])
1738 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1743 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1744 DDRANK(dd, n), n, dd->mpi_comm_all);
1749 n = DDMASTERRANK(dd);
1751 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1753 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1755 copy_rvec(v[c], lv[a++]);
1762 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1763 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1768 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1771 gmx_domdec_master_t *ma;
1772 int *scounts = NULL, *disps = NULL;
1773 int n, i, c, a, nalloc = 0;
1780 get_commbuffer_counts(dd, &scounts, &disps);
1784 for (n = 0; n < dd->nnodes; n++)
1786 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1788 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1790 copy_rvec(v[c], buf[a++]);
1796 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1799 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1801 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1803 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1807 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1811 static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
1814 dd_bcast(dd, sizeof(int), &dfhist->bEquil);
1815 dd_bcast(dd, sizeof(int), &dfhist->nlambda);
1816 dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
1818 if (dfhist->nlambda > 0)
1820 int nlam = dfhist->nlambda;
1821 dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
1822 dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
1823 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
1824 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
1825 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
1826 dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
1828 for (i = 0; i<nlam; i++) {
1829 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
1830 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
1831 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
1832 dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
1833 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
1834 dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
1839 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1840 t_state *state, t_state *state_local,
1845 nh = state->nhchainlength;
1849 for (i = 0; i < efptNR; i++)
1851 state_local->lambda[i] = state->lambda[i];
1853 state_local->fep_state = state->fep_state;
1854 state_local->veta = state->veta;
1855 state_local->vol0 = state->vol0;
1856 copy_mat(state->box, state_local->box);
1857 copy_mat(state->box_rel, state_local->box_rel);
1858 copy_mat(state->boxv, state_local->boxv);
1859 copy_mat(state->svir_prev, state_local->svir_prev);
1860 copy_mat(state->fvir_prev, state_local->fvir_prev);
1861 copy_df_history(&state_local->dfhist,&state->dfhist);
1862 for (i = 0; i < state_local->ngtc; i++)
1864 for (j = 0; j < nh; j++)
1866 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1867 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1869 state_local->therm_integral[i] = state->therm_integral[i];
1871 for (i = 0; i < state_local->nnhpres; i++)
1873 for (j = 0; j < nh; j++)
1875 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1876 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1880 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1881 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1882 dd_bcast(dd, sizeof(real), &state_local->veta);
1883 dd_bcast(dd, sizeof(real), &state_local->vol0);
1884 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1885 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1886 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1887 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1888 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1889 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1890 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1891 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1892 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1893 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1895 /* communicate df_history -- required for restarting from checkpoint */
1896 dd_distribute_dfhist(dd,&state_local->dfhist);
1898 if (dd->nat_home > state_local->nalloc)
1900 dd_realloc_state(state_local, f, dd->nat_home);
1902 for (i = 0; i < estNR; i++)
1904 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1909 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1912 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1915 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1918 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1921 if (state->nrngi == 1)
1924 state_local->nrng*sizeof(state_local->ld_rng[0]),
1925 state->ld_rng, state_local->ld_rng);
1930 state_local->nrng*sizeof(state_local->ld_rng[0]),
1931 state->ld_rng, state_local->ld_rng);
1935 if (state->nrngi == 1)
1937 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1938 state->ld_rngi, state_local->ld_rngi);
1942 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1943 state->ld_rngi, state_local->ld_rngi);
1946 case estDISRE_INITF:
1947 case estDISRE_RM3TAV:
1948 case estORIRE_INITF:
1950 /* Not implemented yet */
1953 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1959 static char dim2char(int dim)
1965 case XX: c = 'X'; break;
1966 case YY: c = 'Y'; break;
1967 case ZZ: c = 'Z'; break;
1968 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1974 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1975 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1977 rvec grid_s[2], *grid_r = NULL, cx, r;
1978 char fname[STRLEN], format[STRLEN], buf[22];
1980 int a, i, d, z, y, x;
1984 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1985 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1989 snew(grid_r, 2*dd->nnodes);
1992 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1996 for (d = 0; d < DIM; d++)
1998 for (i = 0; i < DIM; i++)
2006 if (d < ddbox->npbcdim && dd->nc[d] > 1)
2008 tric[d][i] = box[i][d]/box[i][i];
2017 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
2018 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2019 out = gmx_fio_fopen(fname, "w");
2020 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2022 for (i = 0; i < dd->nnodes; i++)
2024 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
2025 for (d = 0; d < DIM; d++)
2027 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
2029 for (z = 0; z < 2; z++)
2031 for (y = 0; y < 2; y++)
2033 for (x = 0; x < 2; x++)
2035 cx[XX] = grid_r[i*2+x][XX];
2036 cx[YY] = grid_r[i*2+y][YY];
2037 cx[ZZ] = grid_r[i*2+z][ZZ];
2039 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
2040 ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
2044 for (d = 0; d < DIM; d++)
2046 for (x = 0; x < 4; x++)
2050 case 0: y = 1 + i*8 + 2*x; break;
2051 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2052 case 2: y = 1 + i*8 + x; break;
2054 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2058 gmx_fio_fclose(out);
2063 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2064 gmx_mtop_t *mtop, t_commrec *cr,
2065 int natoms, rvec x[], matrix box)
2067 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2069 int i, ii, resnr, c;
2070 char *atomname, *resname;
2077 natoms = dd->comm->nat[ddnatVSITE];
2080 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2082 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2083 sprintf(format4, "%s%s\n", pdbformat4, "%6.2f%6.2f");
2085 out = gmx_fio_fopen(fname, "w");
2087 fprintf(out, "TITLE %s\n", title);
2088 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2089 for (i = 0; i < natoms; i++)
2091 ii = dd->gatindex[i];
2092 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2093 if (i < dd->comm->nat[ddnatZONE])
2096 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2102 else if (i < dd->comm->nat[ddnatVSITE])
2104 b = dd->comm->zones.n;
2108 b = dd->comm->zones.n + 1;
2110 fprintf(out, strlen(atomname) < 4 ? format : format4,
2111 "ATOM", (ii+1)%100000,
2112 atomname, resname, ' ', resnr%10000, ' ',
2113 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2115 fprintf(out, "TER\n");
2117 gmx_fio_fclose(out);
2120 real dd_cutoff_mbody(gmx_domdec_t *dd)
2122 gmx_domdec_comm_t *comm;
2129 if (comm->bInterCGBondeds)
2131 if (comm->cutoff_mbody > 0)
2133 r = comm->cutoff_mbody;
2137 /* cutoff_mbody=0 means we do not have DLB */
2138 r = comm->cellsize_min[dd->dim[0]];
2139 for (di = 1; di < dd->ndim; di++)
2141 r = min(r, comm->cellsize_min[dd->dim[di]]);
2143 if (comm->bBondComm)
2145 r = max(r, comm->cutoff_mbody);
2149 r = min(r, comm->cutoff);
2157 real dd_cutoff_twobody(gmx_domdec_t *dd)
2161 r_mb = dd_cutoff_mbody(dd);
2163 return max(dd->comm->cutoff, r_mb);
2167 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2171 nc = dd->nc[dd->comm->cartpmedim];
2172 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2173 copy_ivec(coord, coord_pme);
2174 coord_pme[dd->comm->cartpmedim] =
2175 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2178 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2180 /* Here we assign a PME node to communicate with this DD node
2181 * by assuming that the major index of both is x.
2182 * We add cr->npmenodes/2 to obtain an even distribution.
2184 return (ddindex*npme + npme/2)/ndd;
2187 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2189 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2192 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2194 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2197 static int *dd_pmenodes(t_commrec *cr)
2202 snew(pmenodes, cr->npmenodes);
2204 for (i = 0; i < cr->dd->nnodes; i++)
2206 p0 = cr_ddindex2pmeindex(cr, i);
2207 p1 = cr_ddindex2pmeindex(cr, i+1);
2208 if (i+1 == cr->dd->nnodes || p1 > p0)
2212 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2214 pmenodes[n] = i + 1 + n;
2222 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2225 ivec coords, coords_pme, nc;
2230 if (dd->comm->bCartesian) {
2231 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2232 dd_coords2pmecoords(dd,coords,coords_pme);
2233 copy_ivec(dd->ntot,nc);
2234 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2235 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2237 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2239 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2245 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2250 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2252 gmx_domdec_comm_t *comm;
2254 int ddindex, nodeid = -1;
2256 comm = cr->dd->comm;
2261 if (comm->bCartesianPP_PME)
2264 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2269 ddindex = dd_index(cr->dd->nc, coords);
2270 if (comm->bCartesianPP)
2272 nodeid = comm->ddindex2simnodeid[ddindex];
2278 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2290 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2293 gmx_domdec_comm_t *comm;
2294 ivec coord, coord_pme;
2301 /* This assumes a uniform x domain decomposition grid cell size */
2302 if (comm->bCartesianPP_PME)
2305 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2306 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2308 /* This is a PP node */
2309 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2310 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2314 else if (comm->bCartesianPP)
2316 if (sim_nodeid < dd->nnodes)
2318 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2323 /* This assumes DD cells with identical x coordinates
2324 * are numbered sequentially.
2326 if (dd->comm->pmenodes == NULL)
2328 if (sim_nodeid < dd->nnodes)
2330 /* The DD index equals the nodeid */
2331 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2337 while (sim_nodeid > dd->comm->pmenodes[i])
2341 if (sim_nodeid < dd->comm->pmenodes[i])
2343 pmenode = dd->comm->pmenodes[i];
2351 void get_pme_nnodes(const gmx_domdec_t *dd,
2352 int *npmenodes_x, int *npmenodes_y)
2356 *npmenodes_x = dd->comm->npmenodes_x;
2357 *npmenodes_y = dd->comm->npmenodes_y;
2366 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2368 gmx_bool bPMEOnlyNode;
2370 if (DOMAINDECOMP(cr))
2372 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2376 bPMEOnlyNode = FALSE;
2379 return bPMEOnlyNode;
2382 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2383 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2387 ivec coord, coord_pme;
2391 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2394 for (x = 0; x < dd->nc[XX]; x++)
2396 for (y = 0; y < dd->nc[YY]; y++)
2398 for (z = 0; z < dd->nc[ZZ]; z++)
2400 if (dd->comm->bCartesianPP_PME)
2405 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2406 if (dd->ci[XX] == coord_pme[XX] &&
2407 dd->ci[YY] == coord_pme[YY] &&
2408 dd->ci[ZZ] == coord_pme[ZZ])
2410 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2415 /* The slab corresponds to the nodeid in the PME group */
2416 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2418 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2425 /* The last PP-only node is the peer node */
2426 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2430 fprintf(debug, "Receive coordinates from PP nodes:");
2431 for (x = 0; x < *nmy_ddnodes; x++)
2433 fprintf(debug, " %d", (*my_ddnodes)[x]);
2435 fprintf(debug, "\n");
2439 static gmx_bool receive_vir_ener(t_commrec *cr)
2441 gmx_domdec_comm_t *comm;
2442 int pmenode, coords[DIM], rank;
2446 if (cr->npmenodes < cr->dd->nnodes)
2448 comm = cr->dd->comm;
2449 if (comm->bCartesianPP_PME)
2451 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2453 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2454 coords[comm->cartpmedim]++;
2455 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2457 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2458 if (dd_simnode2pmenode(cr, rank) == pmenode)
2460 /* This is not the last PP node for pmenode */
2468 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2469 if (cr->sim_nodeid+1 < cr->nnodes &&
2470 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2472 /* This is not the last PP node for pmenode */
2481 static void set_zones_ncg_home(gmx_domdec_t *dd)
2483 gmx_domdec_zones_t *zones;
2486 zones = &dd->comm->zones;
2488 zones->cg_range[0] = 0;
2489 for (i = 1; i < zones->n+1; i++)
2491 zones->cg_range[i] = dd->ncg_home;
2493 /* zone_ncg1[0] should always be equal to ncg_home */
2494 dd->comm->zone_ncg1[0] = dd->ncg_home;
2497 static void rebuild_cgindex(gmx_domdec_t *dd,
2498 const int *gcgs_index, t_state *state)
2500 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2503 dd_cg_gl = dd->index_gl;
2504 cgindex = dd->cgindex;
2507 for (i = 0; i < state->ncg_gl; i++)
2511 dd_cg_gl[i] = cg_gl;
2512 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2516 dd->ncg_home = state->ncg_gl;
2519 set_zones_ncg_home(dd);
2522 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2524 while (cg >= cginfo_mb->cg_end)
2529 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2532 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2533 t_forcerec *fr, char *bLocalCG)
2535 cginfo_mb_t *cginfo_mb;
2541 cginfo_mb = fr->cginfo_mb;
2542 cginfo = fr->cginfo;
2544 for (cg = cg0; cg < cg1; cg++)
2546 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2550 if (bLocalCG != NULL)
2552 for (cg = cg0; cg < cg1; cg++)
2554 bLocalCG[index_gl[cg]] = TRUE;
2559 static void make_dd_indices(gmx_domdec_t *dd,
2560 const int *gcgs_index, int cg_start)
2562 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2563 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2568 bLocalCG = dd->comm->bLocalCG;
2570 if (dd->nat_tot > dd->gatindex_nalloc)
2572 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2573 srenew(dd->gatindex, dd->gatindex_nalloc);
2576 nzone = dd->comm->zones.n;
2577 zone2cg = dd->comm->zones.cg_range;
2578 zone_ncg1 = dd->comm->zone_ncg1;
2579 index_gl = dd->index_gl;
2580 gatindex = dd->gatindex;
2581 bCGs = dd->comm->bCGs;
2583 if (zone2cg[1] != dd->ncg_home)
2585 gmx_incons("dd->ncg_zone is not up to date");
2588 /* Make the local to global and global to local atom index */
2589 a = dd->cgindex[cg_start];
2590 for (zone = 0; zone < nzone; zone++)
2598 cg0 = zone2cg[zone];
2600 cg1 = zone2cg[zone+1];
2601 cg1_p1 = cg0 + zone_ncg1[zone];
2603 for (cg = cg0; cg < cg1; cg++)
2608 /* Signal that this cg is from more than one pulse away */
2611 cg_gl = index_gl[cg];
2614 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2617 ga2la_set(dd->ga2la, a_gl, a, zone1);
2623 gatindex[a] = cg_gl;
2624 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2631 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2634 int ncg, i, ngl, nerr;
2637 if (bLocalCG == NULL)
2641 for (i = 0; i < dd->ncg_tot; i++)
2643 if (!bLocalCG[dd->index_gl[i]])
2646 "DD node %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);
2651 for (i = 0; i < ncg_sys; i++)
2658 if (ngl != dd->ncg_tot)
2660 fprintf(stderr, "DD node %d, %s: In bLocalCG %d cgs are marked as local, whereas there are %d\n", dd->rank, where, ngl, dd->ncg_tot);
2667 static void check_index_consistency(gmx_domdec_t *dd,
2668 int natoms_sys, int ncg_sys,
2671 int nerr, ngl, i, a, cell;
2676 if (dd->comm->DD_debug > 1)
2678 snew(have, natoms_sys);
2679 for (a = 0; a < dd->nat_tot; a++)
2681 if (have[dd->gatindex[a]] > 0)
2683 fprintf(stderr, "DD node %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
2687 have[dd->gatindex[a]] = a + 1;
2693 snew(have, dd->nat_tot);
2696 for (i = 0; i < natoms_sys; i++)
2698 if (ga2la_get(dd->ga2la, i, &a, &cell))
2700 if (a >= dd->nat_tot)
2702 fprintf(stderr, "DD node %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);
2708 if (dd->gatindex[a] != i)
2710 fprintf(stderr, "DD node %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);
2717 if (ngl != dd->nat_tot)
2720 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2721 dd->rank, where, ngl, dd->nat_tot);
2723 for (a = 0; a < dd->nat_tot; a++)
2728 "DD node %d, %s: local atom %d, global %d has no global index\n",
2729 dd->rank, where, a+1, dd->gatindex[a]+1);
2734 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2738 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2739 dd->rank, where, nerr);
2743 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2750 /* Clear the whole list without searching */
2751 ga2la_clear(dd->ga2la);
2755 for (i = a_start; i < dd->nat_tot; i++)
2757 ga2la_del(dd->ga2la, dd->gatindex[i]);
2761 bLocalCG = dd->comm->bLocalCG;
2764 for (i = cg_start; i < dd->ncg_tot; i++)
2766 bLocalCG[dd->index_gl[i]] = FALSE;
2770 dd_clear_local_vsite_indices(dd);
2772 if (dd->constraints)
2774 dd_clear_local_constraint_indices(dd);
2778 /* This function should be used for moving the domain boudaries during DLB,
2779 * for obtaining the minimum cell size. It checks the initially set limit
2780 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2781 * and, possibly, a longer cut-off limit set for PME load balancing.
2783 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2787 cellsize_min = comm->cellsize_min[dim];
2789 if (!comm->bVacDLBNoLimit)
2791 /* The cut-off might have changed, e.g. by PME load balacning,
2792 * from the value used to set comm->cellsize_min, so check it.
2794 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2796 if (comm->bPMELoadBalDLBLimits)
2798 /* Check for the cut-off limit set by the PME load balancing */
2799 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2803 return cellsize_min;
2806 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2809 real grid_jump_limit;
2811 /* The distance between the boundaries of cells at distance
2812 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2813 * and by the fact that cells should not be shifted by more than
2814 * half their size, such that cg's only shift by one cell
2815 * at redecomposition.
2817 grid_jump_limit = comm->cellsize_limit;
2818 if (!comm->bVacDLBNoLimit)
2820 if (comm->bPMELoadBalDLBLimits)
2822 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2824 grid_jump_limit = max(grid_jump_limit,
2825 cutoff/comm->cd[dim_ind].np);
2828 return grid_jump_limit;
2831 static gmx_bool check_grid_jump(gmx_large_int_t step,
2837 gmx_domdec_comm_t *comm;
2846 for (d = 1; d < dd->ndim; d++)
2849 limit = grid_jump_limit(comm, cutoff, d);
2850 bfac = ddbox->box_size[dim];
2851 if (ddbox->tric_dir[dim])
2853 bfac *= ddbox->skew_fac[dim];
2855 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2856 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2864 /* This error should never be triggered under normal
2865 * circumstances, but you never know ...
2867 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 less nodes might avoid this issue.",
2868 gmx_step_str(step, buf),
2869 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2877 static int dd_load_count(gmx_domdec_comm_t *comm)
2879 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2882 static float dd_force_load(gmx_domdec_comm_t *comm)
2889 if (comm->eFlop > 1)
2891 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2896 load = comm->cycl[ddCyclF];
2897 if (comm->cycl_n[ddCyclF] > 1)
2899 /* Subtract the maximum of the last n cycle counts
2900 * to get rid of possible high counts due to other soures,
2901 * for instance system activity, that would otherwise
2902 * affect the dynamic load balancing.
2904 load -= comm->cycl_max[ddCyclF];
2911 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2913 gmx_domdec_comm_t *comm;
2918 snew(*dim_f, dd->nc[dim]+1);
2920 for (i = 1; i < dd->nc[dim]; i++)
2922 if (comm->slb_frac[dim])
2924 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2928 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2931 (*dim_f)[dd->nc[dim]] = 1;
2934 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2936 int pmeindex, slab, nso, i;
2939 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2945 ddpme->dim = dimind;
2947 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2949 ddpme->nslab = (ddpme->dim == 0 ?
2950 dd->comm->npmenodes_x :
2951 dd->comm->npmenodes_y);
2953 if (ddpme->nslab <= 1)
2958 nso = dd->comm->npmenodes/ddpme->nslab;
2959 /* Determine for each PME slab the PP location range for dimension dim */
2960 snew(ddpme->pp_min, ddpme->nslab);
2961 snew(ddpme->pp_max, ddpme->nslab);
2962 for (slab = 0; slab < ddpme->nslab; slab++)
2964 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2965 ddpme->pp_max[slab] = 0;
2967 for (i = 0; i < dd->nnodes; i++)
2969 ddindex2xyz(dd->nc, i, xyz);
2970 /* For y only use our y/z slab.
2971 * This assumes that the PME x grid size matches the DD grid size.
2973 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2975 pmeindex = ddindex2pmeindex(dd, i);
2978 slab = pmeindex/nso;
2982 slab = pmeindex % ddpme->nslab;
2984 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2985 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2989 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2992 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2994 if (dd->comm->ddpme[0].dim == XX)
2996 return dd->comm->ddpme[0].maxshift;
3004 int dd_pme_maxshift_y(gmx_domdec_t *dd)
3006 if (dd->comm->ddpme[0].dim == YY)
3008 return dd->comm->ddpme[0].maxshift;
3010 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
3012 return dd->comm->ddpme[1].maxshift;
3020 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
3021 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
3023 gmx_domdec_comm_t *comm;
3026 real range, pme_boundary;
3030 nc = dd->nc[ddpme->dim];
3033 if (!ddpme->dim_match)
3035 /* PP decomposition is not along dim: the worst situation */
3038 else if (ns <= 3 || (bUniform && ns == nc))
3040 /* The optimal situation */
3045 /* We need to check for all pme nodes which nodes they
3046 * could possibly need to communicate with.
3048 xmin = ddpme->pp_min;
3049 xmax = ddpme->pp_max;
3050 /* Allow for atoms to be maximally 2/3 times the cut-off
3051 * out of their DD cell. This is a reasonable balance between
3052 * between performance and support for most charge-group/cut-off
3055 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3056 /* Avoid extra communication when we are exactly at a boundary */
3060 for (s = 0; s < ns; s++)
3062 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3063 pme_boundary = (real)s/ns;
3066 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3068 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3072 pme_boundary = (real)(s+1)/ns;
3075 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3077 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3084 ddpme->maxshift = sh;
3088 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3089 ddpme->dim, ddpme->maxshift);
3093 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3097 for (d = 0; d < dd->ndim; d++)
3100 if (dim < ddbox->nboundeddim &&
3101 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3102 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3104 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",
3105 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3106 dd->nc[dim], dd->comm->cellsize_limit);
3111 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3112 gmx_bool bMaster, ivec npulse)
3114 gmx_domdec_comm_t *comm;
3117 real *cell_x, cell_dx, cellsize;
3121 for (d = 0; d < DIM; d++)
3123 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3125 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3128 cell_dx = ddbox->box_size[d]/dd->nc[d];
3131 for (j = 0; j < dd->nc[d]+1; j++)
3133 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3138 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3139 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3141 cellsize = cell_dx*ddbox->skew_fac[d];
3142 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3146 cellsize_min[d] = cellsize;
3150 /* Statically load balanced grid */
3151 /* Also when we are not doing a master distribution we determine
3152 * all cell borders in a loop to obtain identical values
3153 * to the master distribution case and to determine npulse.
3157 cell_x = dd->ma->cell_x[d];
3161 snew(cell_x, dd->nc[d]+1);
3163 cell_x[0] = ddbox->box0[d];
3164 for (j = 0; j < dd->nc[d]; j++)
3166 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3167 cell_x[j+1] = cell_x[j] + cell_dx;
3168 cellsize = cell_dx*ddbox->skew_fac[d];
3169 while (cellsize*npulse[d] < comm->cutoff &&
3170 npulse[d] < dd->nc[d]-1)
3174 cellsize_min[d] = min(cellsize_min[d], cellsize);
3178 comm->cell_x0[d] = cell_x[dd->ci[d]];
3179 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3183 /* The following limitation is to avoid that a cell would receive
3184 * some of its own home charge groups back over the periodic boundary.
3185 * Double charge groups cause trouble with the global indices.
3187 if (d < ddbox->npbcdim &&
3188 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3190 gmx_fatal_collective(FARGS, NULL, dd,
3191 "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",
3192 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3194 dd->nc[d], dd->nc[d],
3195 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3199 if (!comm->bDynLoadBal)
3201 copy_rvec(cellsize_min, comm->cellsize_min);
3204 for (d = 0; d < comm->npmedecompdim; d++)
3206 set_pme_maxshift(dd, &comm->ddpme[d],
3207 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3208 comm->ddpme[d].slb_dim_f);
3213 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3214 int d, int dim, gmx_domdec_root_t *root,
3216 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3218 gmx_domdec_comm_t *comm;
3219 int ncd, i, j, nmin, nmin_old;
3220 gmx_bool bLimLo, bLimHi;
3222 real fac, halfway, cellsize_limit_f_i, region_size;
3223 gmx_bool bPBC, bLastHi = FALSE;
3224 int nrange[] = {range[0], range[1]};
3226 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3232 bPBC = (dim < ddbox->npbcdim);
3234 cell_size = root->buf_ncd;
3238 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3241 /* First we need to check if the scaling does not make cells
3242 * smaller than the smallest allowed size.
3243 * We need to do this iteratively, since if a cell is too small,
3244 * it needs to be enlarged, which makes all the other cells smaller,
3245 * which could in turn make another cell smaller than allowed.
3247 for (i = range[0]; i < range[1]; i++)
3249 root->bCellMin[i] = FALSE;
3255 /* We need the total for normalization */
3257 for (i = range[0]; i < range[1]; i++)
3259 if (root->bCellMin[i] == FALSE)
3261 fac += cell_size[i];
3264 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3265 /* Determine the cell boundaries */
3266 for (i = range[0]; i < range[1]; i++)
3268 if (root->bCellMin[i] == FALSE)
3270 cell_size[i] *= fac;
3271 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3273 cellsize_limit_f_i = 0;
3277 cellsize_limit_f_i = cellsize_limit_f;
3279 if (cell_size[i] < cellsize_limit_f_i)
3281 root->bCellMin[i] = TRUE;
3282 cell_size[i] = cellsize_limit_f_i;
3286 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3289 while (nmin > nmin_old);
3292 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3293 /* For this check we should not use DD_CELL_MARGIN,
3294 * but a slightly smaller factor,
3295 * since rounding could get use below the limit.
3297 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3300 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",
3301 gmx_step_str(step, buf),
3302 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3303 ncd, comm->cellsize_min[dim]);
3306 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3310 /* Check if the boundary did not displace more than halfway
3311 * each of the cells it bounds, as this could cause problems,
3312 * especially when the differences between cell sizes are large.
3313 * If changes are applied, they will not make cells smaller
3314 * than the cut-off, as we check all the boundaries which
3315 * might be affected by a change and if the old state was ok,
3316 * the cells will at most be shrunk back to their old size.
3318 for (i = range[0]+1; i < range[1]; i++)
3320 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3321 if (root->cell_f[i] < halfway)
3323 root->cell_f[i] = halfway;
3324 /* Check if the change also causes shifts of the next boundaries */
3325 for (j = i+1; j < range[1]; j++)
3327 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3329 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3333 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3334 if (root->cell_f[i] > halfway)
3336 root->cell_f[i] = halfway;
3337 /* Check if the change also causes shifts of the next boundaries */
3338 for (j = i-1; j >= range[0]+1; j--)
3340 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3342 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3349 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3350 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3351 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3352 * for a and b nrange is used */
3355 /* Take care of the staggering of the cell boundaries */
3358 for (i = range[0]; i < range[1]; i++)
3360 root->cell_f_max0[i] = root->cell_f[i];
3361 root->cell_f_min1[i] = root->cell_f[i+1];
3366 for (i = range[0]+1; i < range[1]; i++)
3368 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3369 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3370 if (bLimLo && bLimHi)
3372 /* Both limits violated, try the best we can */
3373 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3374 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3375 nrange[0] = range[0];
3377 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3380 nrange[1] = range[1];
3381 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3387 /* root->cell_f[i] = root->bound_min[i]; */
3388 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3391 else if (bLimHi && !bLastHi)
3394 if (nrange[1] < range[1]) /* found a LimLo before */
3396 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3397 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3398 nrange[0] = nrange[1];
3400 root->cell_f[i] = root->bound_max[i];
3402 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3404 nrange[1] = range[1];
3407 if (nrange[1] < range[1]) /* found last a LimLo */
3409 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3410 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3411 nrange[0] = nrange[1];
3412 nrange[1] = range[1];
3413 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3415 else if (nrange[0] > range[0]) /* found at least one LimHi */
3417 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3424 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3425 int d, int dim, gmx_domdec_root_t *root,
3426 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3427 gmx_bool bUniform, gmx_large_int_t step)
3429 gmx_domdec_comm_t *comm;
3430 int ncd, d1, i, j, pos;
3432 real load_aver, load_i, imbalance, change, change_max, sc;
3433 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3437 int range[] = { 0, 0 };
3441 /* Convert the maximum change from the input percentage to a fraction */
3442 change_limit = comm->dlb_scale_lim*0.01;
3446 bPBC = (dim < ddbox->npbcdim);
3448 cell_size = root->buf_ncd;
3450 /* Store the original boundaries */
3451 for (i = 0; i < ncd+1; i++)
3453 root->old_cell_f[i] = root->cell_f[i];
3457 for (i = 0; i < ncd; i++)
3459 cell_size[i] = 1.0/ncd;
3462 else if (dd_load_count(comm))
3464 load_aver = comm->load[d].sum_m/ncd;
3466 for (i = 0; i < ncd; i++)
3468 /* Determine the relative imbalance of cell i */
3469 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3470 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3471 /* Determine the change of the cell size using underrelaxation */
3472 change = -relax*imbalance;
3473 change_max = max(change_max, max(change, -change));
3475 /* Limit the amount of scaling.
3476 * We need to use the same rescaling for all cells in one row,
3477 * otherwise the load balancing might not converge.
3480 if (change_max > change_limit)
3482 sc *= change_limit/change_max;
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 = -sc*imbalance;
3491 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3495 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3496 cellsize_limit_f *= DD_CELL_MARGIN;
3497 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3498 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3499 if (ddbox->tric_dir[dim])
3501 cellsize_limit_f /= ddbox->skew_fac[dim];
3502 dist_min_f /= ddbox->skew_fac[dim];
3504 if (bDynamicBox && d > 0)
3506 dist_min_f *= DD_PRES_SCALE_MARGIN;
3508 if (d > 0 && !bUniform)
3510 /* Make sure that the grid is not shifted too much */
3511 for (i = 1; i < ncd; i++)
3513 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3515 gmx_incons("Inconsistent DD boundary staggering limits!");
3517 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3518 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3521 root->bound_min[i] += 0.5*space;
3523 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3524 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3527 root->bound_max[i] += 0.5*space;
3532 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3534 root->cell_f_max0[i-1] + dist_min_f,
3535 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3536 root->cell_f_min1[i] - dist_min_f);
3541 root->cell_f[0] = 0;
3542 root->cell_f[ncd] = 1;
3543 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3546 /* After the checks above, the cells should obey the cut-off
3547 * restrictions, but it does not hurt to check.
3549 for (i = 0; i < ncd; i++)
3553 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3554 dim, i, root->cell_f[i], root->cell_f[i+1]);
3557 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3558 root->cell_f[i+1] - root->cell_f[i] <
3559 cellsize_limit_f/DD_CELL_MARGIN)
3563 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3564 gmx_step_str(step, buf), dim2char(dim), i,
3565 (root->cell_f[i+1] - root->cell_f[i])
3566 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3571 /* Store the cell boundaries of the lower dimensions at the end */
3572 for (d1 = 0; d1 < d; d1++)
3574 root->cell_f[pos++] = comm->cell_f0[d1];
3575 root->cell_f[pos++] = comm->cell_f1[d1];
3578 if (d < comm->npmedecompdim)
3580 /* The master determines the maximum shift for
3581 * the coordinate communication between separate PME nodes.
3583 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3585 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3588 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3592 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3593 gmx_ddbox_t *ddbox, int dimind)
3595 gmx_domdec_comm_t *comm;
3600 /* Set the cell dimensions */
3601 dim = dd->dim[dimind];
3602 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3603 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3604 if (dim >= ddbox->nboundeddim)
3606 comm->cell_x0[dim] += ddbox->box0[dim];
3607 comm->cell_x1[dim] += ddbox->box0[dim];
3611 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3612 int d, int dim, real *cell_f_row,
3615 gmx_domdec_comm_t *comm;
3621 /* Each node would only need to know two fractions,
3622 * but it is probably cheaper to broadcast the whole array.
3624 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3625 0, comm->mpi_comm_load[d]);
3627 /* Copy the fractions for this dimension from the buffer */
3628 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3629 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3630 /* The whole array was communicated, so set the buffer position */
3631 pos = dd->nc[dim] + 1;
3632 for (d1 = 0; d1 <= d; d1++)
3636 /* Copy the cell fractions of the lower dimensions */
3637 comm->cell_f0[d1] = cell_f_row[pos++];
3638 comm->cell_f1[d1] = cell_f_row[pos++];
3640 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3642 /* Convert the communicated shift from float to int */
3643 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3646 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3650 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3651 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3652 gmx_bool bUniform, gmx_large_int_t step)
3654 gmx_domdec_comm_t *comm;
3656 gmx_bool bRowMember, bRowRoot;
3661 for (d = 0; d < dd->ndim; d++)
3666 for (d1 = d; d1 < dd->ndim; d1++)
3668 if (dd->ci[dd->dim[d1]] > 0)
3681 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3682 ddbox, bDynamicBox, bUniform, step);
3683 cell_f_row = comm->root[d]->cell_f;
3687 cell_f_row = comm->cell_f_row;
3689 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3694 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3698 /* This function assumes the box is static and should therefore
3699 * not be called when the box has changed since the last
3700 * call to dd_partition_system.
3702 for (d = 0; d < dd->ndim; d++)
3704 relative_to_absolute_cell_bounds(dd, ddbox, d);
3710 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3711 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3712 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3713 gmx_wallcycle_t wcycle)
3715 gmx_domdec_comm_t *comm;
3722 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3723 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3724 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3726 else if (bDynamicBox)
3728 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3731 /* Set the dimensions for which no DD is used */
3732 for (dim = 0; dim < DIM; dim++)
3734 if (dd->nc[dim] == 1)
3736 comm->cell_x0[dim] = 0;
3737 comm->cell_x1[dim] = ddbox->box_size[dim];
3738 if (dim >= ddbox->nboundeddim)
3740 comm->cell_x0[dim] += ddbox->box0[dim];
3741 comm->cell_x1[dim] += ddbox->box0[dim];
3747 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3750 gmx_domdec_comm_dim_t *cd;
3752 for (d = 0; d < dd->ndim; d++)
3754 cd = &dd->comm->cd[d];
3755 np = npulse[dd->dim[d]];
3756 if (np > cd->np_nalloc)
3760 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3761 dim2char(dd->dim[d]), np);
3763 if (DDMASTER(dd) && cd->np_nalloc > 0)
3765 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3767 srenew(cd->ind, np);
3768 for (i = cd->np_nalloc; i < np; i++)
3770 cd->ind[i].index = NULL;
3771 cd->ind[i].nalloc = 0;
3780 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3781 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3782 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3783 gmx_wallcycle_t wcycle)
3785 gmx_domdec_comm_t *comm;
3791 /* Copy the old cell boundaries for the cg displacement check */
3792 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3793 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3795 if (comm->bDynLoadBal)
3799 check_box_size(dd, ddbox);
3801 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3805 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3806 realloc_comm_ind(dd, npulse);
3811 for (d = 0; d < DIM; d++)
3813 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3814 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3819 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3821 rvec cell_ns_x0, rvec cell_ns_x1,
3822 gmx_large_int_t step)
3824 gmx_domdec_comm_t *comm;
3829 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3831 dim = dd->dim[dim_ind];
3833 /* Without PBC we don't have restrictions on the outer cells */
3834 if (!(dim >= ddbox->npbcdim &&
3835 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3836 comm->bDynLoadBal &&
3837 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3838 comm->cellsize_min[dim])
3841 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",
3842 gmx_step_str(step, buf), dim2char(dim),
3843 comm->cell_x1[dim] - comm->cell_x0[dim],
3844 ddbox->skew_fac[dim],
3845 dd->comm->cellsize_min[dim],
3846 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3850 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3852 /* Communicate the boundaries and update cell_ns_x0/1 */
3853 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3854 if (dd->bGridJump && dd->ndim > 1)
3856 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3861 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3865 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3873 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3874 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3883 static void check_screw_box(matrix box)
3885 /* Mathematical limitation */
3886 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3888 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3891 /* Limitation due to the asymmetry of the eighth shell method */
3892 if (box[ZZ][YY] != 0)
3894 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3898 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3899 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3902 gmx_domdec_master_t *ma;
3903 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3904 int i, icg, j, k, k0, k1, d, npbcdim;
3906 rvec box_size, cg_cm;
3908 real nrcg, inv_ncg, pos_d;
3910 gmx_bool bUnbounded, bScrew;
3914 if (tmp_ind == NULL)
3916 snew(tmp_nalloc, dd->nnodes);
3917 snew(tmp_ind, dd->nnodes);
3918 for (i = 0; i < dd->nnodes; i++)
3920 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3921 snew(tmp_ind[i], tmp_nalloc[i]);
3925 /* Clear the count */
3926 for (i = 0; i < dd->nnodes; i++)
3932 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3934 cgindex = cgs->index;
3936 /* Compute the center of geometry for all charge groups */
3937 for (icg = 0; icg < cgs->nr; icg++)
3940 k1 = cgindex[icg+1];
3944 copy_rvec(pos[k0], cg_cm);
3951 for (k = k0; (k < k1); k++)
3953 rvec_inc(cg_cm, pos[k]);
3955 for (d = 0; (d < DIM); d++)
3957 cg_cm[d] *= inv_ncg;
3960 /* Put the charge group in the box and determine the cell index */
3961 for (d = DIM-1; d >= 0; d--)
3964 if (d < dd->npbcdim)
3966 bScrew = (dd->bScrewPBC && d == XX);
3967 if (tric_dir[d] && dd->nc[d] > 1)
3969 /* Use triclinic coordintates for this dimension */
3970 for (j = d+1; j < DIM; j++)
3972 pos_d += cg_cm[j]*tcm[j][d];
3975 while (pos_d >= box[d][d])
3978 rvec_dec(cg_cm, box[d]);
3981 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3982 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3984 for (k = k0; (k < k1); k++)
3986 rvec_dec(pos[k], box[d]);
3989 pos[k][YY] = box[YY][YY] - pos[k][YY];
3990 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3997 rvec_inc(cg_cm, box[d]);
4000 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
4001 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
4003 for (k = k0; (k < k1); k++)
4005 rvec_inc(pos[k], box[d]);
4008 pos[k][YY] = box[YY][YY] - pos[k][YY];
4009 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
4014 /* This could be done more efficiently */
4016 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
4021 i = dd_index(dd->nc, ind);
4022 if (ma->ncg[i] == tmp_nalloc[i])
4024 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
4025 srenew(tmp_ind[i], tmp_nalloc[i]);
4027 tmp_ind[i][ma->ncg[i]] = icg;
4029 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4033 for (i = 0; i < dd->nnodes; i++)
4036 for (k = 0; k < ma->ncg[i]; k++)
4038 ma->cg[k1++] = tmp_ind[i][k];
4041 ma->index[dd->nnodes] = k1;
4043 for (i = 0; i < dd->nnodes; i++)
4053 fprintf(fplog, "Charge group distribution at step %s:",
4054 gmx_step_str(step, buf));
4055 for (i = 0; i < dd->nnodes; i++)
4057 fprintf(fplog, " %d", ma->ncg[i]);
4059 fprintf(fplog, "\n");
4063 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4064 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4067 gmx_domdec_master_t *ma = NULL;
4070 int *ibuf, buf2[2] = { 0, 0 };
4071 gmx_bool bMaster = DDMASTER(dd);
4078 check_screw_box(box);
4081 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4083 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4084 for (i = 0; i < dd->nnodes; i++)
4086 ma->ibuf[2*i] = ma->ncg[i];
4087 ma->ibuf[2*i+1] = ma->nat[i];
4095 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4097 dd->ncg_home = buf2[0];
4098 dd->nat_home = buf2[1];
4099 dd->ncg_tot = dd->ncg_home;
4100 dd->nat_tot = dd->nat_home;
4101 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4103 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4104 srenew(dd->index_gl, dd->cg_nalloc);
4105 srenew(dd->cgindex, dd->cg_nalloc+1);
4109 for (i = 0; i < dd->nnodes; i++)
4111 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4112 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4117 DDMASTER(dd) ? ma->ibuf : NULL,
4118 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4119 DDMASTER(dd) ? ma->cg : NULL,
4120 dd->ncg_home*sizeof(int), dd->index_gl);
4122 /* Determine the home charge group sizes */
4124 for (i = 0; i < dd->ncg_home; i++)
4126 cg_gl = dd->index_gl[i];
4128 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4133 fprintf(debug, "Home charge groups:\n");
4134 for (i = 0; i < dd->ncg_home; i++)
4136 fprintf(debug, " %d", dd->index_gl[i]);
4139 fprintf(debug, "\n");
4142 fprintf(debug, "\n");
4146 static int compact_and_copy_vec_at(int ncg, int *move,
4149 rvec *src, gmx_domdec_comm_t *comm,
4152 int m, icg, i, i0, i1, nrcg;
4158 for (m = 0; m < DIM*2; m++)
4164 for (icg = 0; icg < ncg; icg++)
4166 i1 = cgindex[icg+1];
4172 /* Compact the home array in place */
4173 for (i = i0; i < i1; i++)
4175 copy_rvec(src[i], src[home_pos++]);
4181 /* Copy to the communication buffer */
4183 pos_vec[m] += 1 + vec*nrcg;
4184 for (i = i0; i < i1; i++)
4186 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4188 pos_vec[m] += (nvec - vec - 1)*nrcg;
4192 home_pos += i1 - i0;
4200 static int compact_and_copy_vec_cg(int ncg, int *move,
4202 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4205 int m, icg, i0, i1, nrcg;
4211 for (m = 0; m < DIM*2; m++)
4217 for (icg = 0; icg < ncg; icg++)
4219 i1 = cgindex[icg+1];
4225 /* Compact the home array in place */
4226 copy_rvec(src[icg], src[home_pos++]);
4232 /* Copy to the communication buffer */
4233 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4234 pos_vec[m] += 1 + nrcg*nvec;
4246 static int compact_ind(int ncg, int *move,
4247 int *index_gl, int *cgindex,
4249 gmx_ga2la_t ga2la, char *bLocalCG,
4252 int cg, nat, a0, a1, a, a_gl;
4257 for (cg = 0; cg < ncg; cg++)
4263 /* Compact the home arrays in place.
4264 * Anything that can be done here avoids access to global arrays.
4266 cgindex[home_pos] = nat;
4267 for (a = a0; a < a1; a++)
4270 gatindex[nat] = a_gl;
4271 /* The cell number stays 0, so we don't need to set it */
4272 ga2la_change_la(ga2la, a_gl, nat);
4275 index_gl[home_pos] = index_gl[cg];
4276 cginfo[home_pos] = cginfo[cg];
4277 /* The charge group remains local, so bLocalCG does not change */
4282 /* Clear the global indices */
4283 for (a = a0; a < a1; a++)
4285 ga2la_del(ga2la, gatindex[a]);
4289 bLocalCG[index_gl[cg]] = FALSE;
4293 cgindex[home_pos] = nat;
4298 static void clear_and_mark_ind(int ncg, int *move,
4299 int *index_gl, int *cgindex, int *gatindex,
4300 gmx_ga2la_t ga2la, char *bLocalCG,
4305 for (cg = 0; cg < ncg; cg++)
4311 /* Clear the global indices */
4312 for (a = a0; a < a1; a++)
4314 ga2la_del(ga2la, gatindex[a]);
4318 bLocalCG[index_gl[cg]] = FALSE;
4320 /* Signal that this cg has moved using the ns cell index.
4321 * Here we set it to -1. fill_grid will change it
4322 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4324 cell_index[cg] = -1;
4329 static void print_cg_move(FILE *fplog,
4331 gmx_large_int_t step, int cg, int dim, int dir,
4332 gmx_bool bHaveLimitdAndCMOld, real limitd,
4333 rvec cm_old, rvec cm_new, real pos_d)
4335 gmx_domdec_comm_t *comm;
4340 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4341 if (bHaveLimitdAndCMOld)
4343 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4344 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4348 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4349 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4351 fprintf(fplog, "distance out of cell %f\n",
4352 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4353 if (bHaveLimitdAndCMOld)
4355 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4356 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4358 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4359 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4360 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4362 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4363 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4365 comm->cell_x0[dim], comm->cell_x1[dim]);
4368 static void cg_move_error(FILE *fplog,
4370 gmx_large_int_t step, int cg, int dim, int dir,
4371 gmx_bool bHaveLimitdAndCMOld, real limitd,
4372 rvec cm_old, rvec cm_new, real pos_d)
4376 print_cg_move(fplog, dd, step, cg, dim, dir,
4377 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4379 print_cg_move(stderr, dd, step, cg, dim, dir,
4380 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4382 "A charge group moved too far between two domain decomposition steps\n"
4383 "This usually means that your system is not well equilibrated");
4386 static void rotate_state_atom(t_state *state, int a)
4390 for (est = 0; est < estNR; est++)
4392 if (EST_DISTR(est) && (state->flags & (1<<est)))
4397 /* Rotate the complete state; for a rectangular box only */
4398 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4399 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4402 state->v[a][YY] = -state->v[a][YY];
4403 state->v[a][ZZ] = -state->v[a][ZZ];
4406 state->sd_X[a][YY] = -state->sd_X[a][YY];
4407 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4410 state->cg_p[a][YY] = -state->cg_p[a][YY];
4411 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4413 case estDISRE_INITF:
4414 case estDISRE_RM3TAV:
4415 case estORIRE_INITF:
4417 /* These are distances, so not affected by rotation */
4420 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4426 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4428 if (natoms > comm->moved_nalloc)
4430 /* Contents should be preserved here */
4431 comm->moved_nalloc = over_alloc_dd(natoms);
4432 srenew(comm->moved, comm->moved_nalloc);
4438 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4441 ivec tric_dir, matrix tcm,
4442 rvec cell_x0, rvec cell_x1,
4443 rvec limitd, rvec limit0, rvec limit1,
4445 int cg_start, int cg_end,
4450 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4451 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4455 real inv_ncg, pos_d;
4458 npbcdim = dd->npbcdim;
4460 for (cg = cg_start; cg < cg_end; cg++)
4467 copy_rvec(state->x[k0], cm_new);
4474 for (k = k0; (k < k1); k++)
4476 rvec_inc(cm_new, state->x[k]);
4478 for (d = 0; (d < DIM); d++)
4480 cm_new[d] = inv_ncg*cm_new[d];
4485 /* Do pbc and check DD cell boundary crossings */
4486 for (d = DIM-1; d >= 0; d--)
4490 bScrew = (dd->bScrewPBC && d == XX);
4491 /* Determine the location of this cg in lattice coordinates */
4495 for (d2 = d+1; d2 < DIM; d2++)
4497 pos_d += cm_new[d2]*tcm[d2][d];
4500 /* Put the charge group in the triclinic unit-cell */
4501 if (pos_d >= cell_x1[d])
4503 if (pos_d >= limit1[d])
4505 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4506 cg_cm[cg], cm_new, pos_d);
4509 if (dd->ci[d] == dd->nc[d] - 1)
4511 rvec_dec(cm_new, state->box[d]);
4514 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4515 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4517 for (k = k0; (k < k1); k++)
4519 rvec_dec(state->x[k], state->box[d]);
4522 rotate_state_atom(state, k);
4527 else if (pos_d < cell_x0[d])
4529 if (pos_d < limit0[d])
4531 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4532 cg_cm[cg], cm_new, pos_d);
4537 rvec_inc(cm_new, state->box[d]);
4540 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4541 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4543 for (k = k0; (k < k1); k++)
4545 rvec_inc(state->x[k], state->box[d]);
4548 rotate_state_atom(state, k);
4554 else if (d < npbcdim)
4556 /* Put the charge group in the rectangular unit-cell */
4557 while (cm_new[d] >= state->box[d][d])
4559 rvec_dec(cm_new, state->box[d]);
4560 for (k = k0; (k < k1); k++)
4562 rvec_dec(state->x[k], state->box[d]);
4565 while (cm_new[d] < 0)
4567 rvec_inc(cm_new, state->box[d]);
4568 for (k = k0; (k < k1); k++)
4570 rvec_inc(state->x[k], state->box[d]);
4576 copy_rvec(cm_new, cg_cm[cg]);
4578 /* Determine where this cg should go */
4581 for (d = 0; d < dd->ndim; d++)
4586 flag |= DD_FLAG_FW(d);
4592 else if (dev[dim] == -1)
4594 flag |= DD_FLAG_BW(d);
4597 if (dd->nc[dim] > 2)
4608 /* Temporarily store the flag in move */
4609 move[cg] = mc + flag;
4613 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4614 gmx_domdec_t *dd, ivec tric_dir,
4615 t_state *state, rvec **f,
4616 t_forcerec *fr, t_mdatoms *md,
4624 int ncg[DIM*2], nat[DIM*2];
4625 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4626 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4627 int sbuf[2], rbuf[2];
4628 int home_pos_cg, home_pos_at, buf_pos;
4630 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4633 real inv_ncg, pos_d;
4635 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4637 cginfo_mb_t *cginfo_mb;
4638 gmx_domdec_comm_t *comm;
4640 int nthread, thread;
4644 check_screw_box(state->box);
4648 if (fr->cutoff_scheme == ecutsGROUP)
4653 for (i = 0; i < estNR; i++)
4659 case estX: /* Always present */ break;
4660 case estV: bV = (state->flags & (1<<i)); break;
4661 case estSDX: bSDX = (state->flags & (1<<i)); break;
4662 case estCGP: bCGP = (state->flags & (1<<i)); break;
4665 case estDISRE_INITF:
4666 case estDISRE_RM3TAV:
4667 case estORIRE_INITF:
4669 /* No processing required */
4672 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4677 if (dd->ncg_tot > comm->nalloc_int)
4679 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4680 srenew(comm->buf_int, comm->nalloc_int);
4682 move = comm->buf_int;
4684 /* Clear the count */
4685 for (c = 0; c < dd->ndim*2; c++)
4691 npbcdim = dd->npbcdim;
4693 for (d = 0; (d < DIM); d++)
4695 limitd[d] = dd->comm->cellsize_min[d];
4696 if (d >= npbcdim && dd->ci[d] == 0)
4698 cell_x0[d] = -GMX_FLOAT_MAX;
4702 cell_x0[d] = comm->cell_x0[d];
4704 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4706 cell_x1[d] = GMX_FLOAT_MAX;
4710 cell_x1[d] = comm->cell_x1[d];
4714 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4715 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4719 /* We check after communication if a charge group moved
4720 * more than one cell. Set the pre-comm check limit to float_max.
4722 limit0[d] = -GMX_FLOAT_MAX;
4723 limit1[d] = GMX_FLOAT_MAX;
4727 make_tric_corr_matrix(npbcdim, state->box, tcm);
4729 cgindex = dd->cgindex;
4731 nthread = gmx_omp_nthreads_get(emntDomdec);
4733 /* Compute the center of geometry for all home charge groups
4734 * and put them in the box and determine where they should go.
4736 #pragma omp parallel for num_threads(nthread) schedule(static)
4737 for (thread = 0; thread < nthread; thread++)
4739 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4740 cell_x0, cell_x1, limitd, limit0, limit1,
4742 ( thread *dd->ncg_home)/nthread,
4743 ((thread+1)*dd->ncg_home)/nthread,
4744 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4748 for (cg = 0; cg < dd->ncg_home; cg++)
4753 flag = mc & ~DD_FLAG_NRCG;
4754 mc = mc & DD_FLAG_NRCG;
4757 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4759 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4760 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4762 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4763 /* We store the cg size in the lower 16 bits
4764 * and the place where the charge group should go
4765 * in the next 6 bits. This saves some communication volume.
4767 nrcg = cgindex[cg+1] - cgindex[cg];
4768 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4774 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4775 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4778 for (i = 0; i < dd->ndim*2; i++)
4780 *ncg_moved += ncg[i];
4797 /* Make sure the communication buffers are large enough */
4798 for (mc = 0; mc < dd->ndim*2; mc++)
4800 nvr = ncg[mc] + nat[mc]*nvec;
4801 if (nvr > comm->cgcm_state_nalloc[mc])
4803 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4804 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4808 switch (fr->cutoff_scheme)
4811 /* Recalculating cg_cm might be cheaper than communicating,
4812 * but that could give rise to rounding issues.
4815 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4816 nvec, cg_cm, comm, bCompact);
4819 /* Without charge groups we send the moved atom coordinates
4820 * over twice. This is so the code below can be used without
4821 * many conditionals for both for with and without charge groups.
4824 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4825 nvec, state->x, comm, FALSE);
4828 home_pos_cg -= *ncg_moved;
4832 gmx_incons("unimplemented");
4838 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4839 nvec, vec++, state->x, comm, bCompact);
4842 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4843 nvec, vec++, state->v, comm, bCompact);
4847 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4848 nvec, vec++, state->sd_X, comm, bCompact);
4852 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4853 nvec, vec++, state->cg_p, comm, bCompact);
4858 compact_ind(dd->ncg_home, move,
4859 dd->index_gl, dd->cgindex, dd->gatindex,
4860 dd->ga2la, comm->bLocalCG,
4865 if (fr->cutoff_scheme == ecutsVERLET)
4867 moved = get_moved(comm, dd->ncg_home);
4869 for (k = 0; k < dd->ncg_home; k++)
4876 moved = fr->ns.grid->cell_index;
4879 clear_and_mark_ind(dd->ncg_home, move,
4880 dd->index_gl, dd->cgindex, dd->gatindex,
4881 dd->ga2la, comm->bLocalCG,
4885 cginfo_mb = fr->cginfo_mb;
4887 *ncg_stay_home = home_pos_cg;
4888 for (d = 0; d < dd->ndim; d++)
4894 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4897 /* Communicate the cg and atom counts */
4902 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4903 d, dir, sbuf[0], sbuf[1]);
4905 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4907 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4909 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4910 srenew(comm->buf_int, comm->nalloc_int);
4913 /* Communicate the charge group indices, sizes and flags */
4914 dd_sendrecv_int(dd, d, dir,
4915 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4916 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4918 nvs = ncg[cdd] + nat[cdd]*nvec;
4919 i = rbuf[0] + rbuf[1] *nvec;
4920 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4922 /* Communicate cgcm and state */
4923 dd_sendrecv_rvec(dd, d, dir,
4924 comm->cgcm_state[cdd], nvs,
4925 comm->vbuf.v+nvr, i);
4926 ncg_recv += rbuf[0];
4927 nat_recv += rbuf[1];
4931 /* Process the received charge groups */
4933 for (cg = 0; cg < ncg_recv; cg++)
4935 flag = comm->buf_int[cg*DD_CGIBS+1];
4937 if (dim >= npbcdim && dd->nc[dim] > 2)
4939 /* No pbc in this dim and more than one domain boundary.
4940 * We do a separate check if a charge group didn't move too far.
4942 if (((flag & DD_FLAG_FW(d)) &&
4943 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4944 ((flag & DD_FLAG_BW(d)) &&
4945 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4947 cg_move_error(fplog, dd, step, cg, dim,
4948 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4950 comm->vbuf.v[buf_pos],
4951 comm->vbuf.v[buf_pos],
4952 comm->vbuf.v[buf_pos][dim]);
4959 /* Check which direction this cg should go */
4960 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4964 /* The cell boundaries for dimension d2 are not equal
4965 * for each cell row of the lower dimension(s),
4966 * therefore we might need to redetermine where
4967 * this cg should go.
4970 /* If this cg crosses the box boundary in dimension d2
4971 * we can use the communicated flag, so we do not
4972 * have to worry about pbc.
4974 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4975 (flag & DD_FLAG_FW(d2))) ||
4976 (dd->ci[dim2] == 0 &&
4977 (flag & DD_FLAG_BW(d2)))))
4979 /* Clear the two flags for this dimension */
4980 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4981 /* Determine the location of this cg
4982 * in lattice coordinates
4984 pos_d = comm->vbuf.v[buf_pos][dim2];
4987 for (d3 = dim2+1; d3 < DIM; d3++)
4990 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4993 /* Check of we are not at the box edge.
4994 * pbc is only handled in the first step above,
4995 * but this check could move over pbc while
4996 * the first step did not due to different rounding.
4998 if (pos_d >= cell_x1[dim2] &&
4999 dd->ci[dim2] != dd->nc[dim2]-1)
5001 flag |= DD_FLAG_FW(d2);
5003 else if (pos_d < cell_x0[dim2] &&
5006 flag |= DD_FLAG_BW(d2);
5008 comm->buf_int[cg*DD_CGIBS+1] = flag;
5011 /* Set to which neighboring cell this cg should go */
5012 if (flag & DD_FLAG_FW(d2))
5016 else if (flag & DD_FLAG_BW(d2))
5018 if (dd->nc[dd->dim[d2]] > 2)
5030 nrcg = flag & DD_FLAG_NRCG;
5033 if (home_pos_cg+1 > dd->cg_nalloc)
5035 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5036 srenew(dd->index_gl, dd->cg_nalloc);
5037 srenew(dd->cgindex, dd->cg_nalloc+1);
5039 /* Set the global charge group index and size */
5040 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5041 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5042 /* Copy the state from the buffer */
5043 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5044 if (fr->cutoff_scheme == ecutsGROUP)
5047 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5051 /* Set the cginfo */
5052 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5053 dd->index_gl[home_pos_cg]);
5056 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5059 if (home_pos_at+nrcg > state->nalloc)
5061 dd_realloc_state(state, f, home_pos_at+nrcg);
5063 for (i = 0; i < nrcg; i++)
5065 copy_rvec(comm->vbuf.v[buf_pos++],
5066 state->x[home_pos_at+i]);
5070 for (i = 0; i < nrcg; i++)
5072 copy_rvec(comm->vbuf.v[buf_pos++],
5073 state->v[home_pos_at+i]);
5078 for (i = 0; i < nrcg; i++)
5080 copy_rvec(comm->vbuf.v[buf_pos++],
5081 state->sd_X[home_pos_at+i]);
5086 for (i = 0; i < nrcg; i++)
5088 copy_rvec(comm->vbuf.v[buf_pos++],
5089 state->cg_p[home_pos_at+i]);
5093 home_pos_at += nrcg;
5097 /* Reallocate the buffers if necessary */
5098 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5100 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5101 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5103 nvr = ncg[mc] + nat[mc]*nvec;
5104 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5106 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5107 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5109 /* Copy from the receive to the send buffers */
5110 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5111 comm->buf_int + cg*DD_CGIBS,
5112 DD_CGIBS*sizeof(int));
5113 memcpy(comm->cgcm_state[mc][nvr],
5114 comm->vbuf.v[buf_pos],
5115 (1+nrcg*nvec)*sizeof(rvec));
5116 buf_pos += 1 + nrcg*nvec;
5123 /* With sorting (!bCompact) the indices are now only partially up to date
5124 * and ncg_home and nat_home are not the real count, since there are
5125 * "holes" in the arrays for the charge groups that moved to neighbors.
5127 if (fr->cutoff_scheme == ecutsVERLET)
5129 moved = get_moved(comm, home_pos_cg);
5131 for (i = dd->ncg_home; i < home_pos_cg; i++)
5136 dd->ncg_home = home_pos_cg;
5137 dd->nat_home = home_pos_at;
5142 "Finished repartitioning: cgs moved out %d, new home %d\n",
5143 *ncg_moved, dd->ncg_home-*ncg_moved);
5148 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5150 dd->comm->cycl[ddCycl] += cycles;
5151 dd->comm->cycl_n[ddCycl]++;
5152 if (cycles > dd->comm->cycl_max[ddCycl])
5154 dd->comm->cycl_max[ddCycl] = cycles;
5158 static double force_flop_count(t_nrnb *nrnb)
5165 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5167 /* To get closer to the real timings, we half the count
5168 * for the normal loops and again half it for water loops.
5171 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5173 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5177 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5180 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5183 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5185 sum += nrnb->n[i]*cost_nrnb(i);
5188 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5190 sum += nrnb->n[i]*cost_nrnb(i);
5196 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5198 if (dd->comm->eFlop)
5200 dd->comm->flop -= force_flop_count(nrnb);
5203 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5205 if (dd->comm->eFlop)
5207 dd->comm->flop += force_flop_count(nrnb);
5212 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5216 for (i = 0; i < ddCyclNr; i++)
5218 dd->comm->cycl[i] = 0;
5219 dd->comm->cycl_n[i] = 0;
5220 dd->comm->cycl_max[i] = 0;
5223 dd->comm->flop_n = 0;
5226 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5228 gmx_domdec_comm_t *comm;
5229 gmx_domdec_load_t *load;
5230 gmx_domdec_root_t *root = NULL;
5231 int d, dim, cid, i, pos;
5232 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5237 fprintf(debug, "get_load_distribution start\n");
5240 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5244 bSepPME = (dd->pme_nodeid >= 0);
5246 for (d = dd->ndim-1; d >= 0; d--)
5249 /* Check if we participate in the communication in this dimension */
5250 if (d == dd->ndim-1 ||
5251 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5253 load = &comm->load[d];
5256 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5259 if (d == dd->ndim-1)
5261 sbuf[pos++] = dd_force_load(comm);
5262 sbuf[pos++] = sbuf[0];
5265 sbuf[pos++] = sbuf[0];
5266 sbuf[pos++] = cell_frac;
5269 sbuf[pos++] = comm->cell_f_max0[d];
5270 sbuf[pos++] = comm->cell_f_min1[d];
5275 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5276 sbuf[pos++] = comm->cycl[ddCyclPME];
5281 sbuf[pos++] = comm->load[d+1].sum;
5282 sbuf[pos++] = comm->load[d+1].max;
5285 sbuf[pos++] = comm->load[d+1].sum_m;
5286 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5287 sbuf[pos++] = comm->load[d+1].flags;
5290 sbuf[pos++] = comm->cell_f_max0[d];
5291 sbuf[pos++] = comm->cell_f_min1[d];
5296 sbuf[pos++] = comm->load[d+1].mdf;
5297 sbuf[pos++] = comm->load[d+1].pme;
5301 /* Communicate a row in DD direction d.
5302 * The communicators are setup such that the root always has rank 0.
5305 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5306 load->load, load->nload*sizeof(float), MPI_BYTE,
5307 0, comm->mpi_comm_load[d]);
5309 if (dd->ci[dim] == dd->master_ci[dim])
5311 /* We are the root, process this row */
5312 if (comm->bDynLoadBal)
5314 root = comm->root[d];
5324 for (i = 0; i < dd->nc[dim]; i++)
5326 load->sum += load->load[pos++];
5327 load->max = max(load->max, load->load[pos]);
5333 /* This direction could not be load balanced properly,
5334 * therefore we need to use the maximum iso the average load.
5336 load->sum_m = max(load->sum_m, load->load[pos]);
5340 load->sum_m += load->load[pos];
5343 load->cvol_min = min(load->cvol_min, load->load[pos]);
5347 load->flags = (int)(load->load[pos++] + 0.5);
5351 root->cell_f_max0[i] = load->load[pos++];
5352 root->cell_f_min1[i] = load->load[pos++];
5357 load->mdf = max(load->mdf, load->load[pos]);
5359 load->pme = max(load->pme, load->load[pos]);
5363 if (comm->bDynLoadBal && root->bLimited)
5365 load->sum_m *= dd->nc[dim];
5366 load->flags |= (1<<d);
5374 comm->nload += dd_load_count(comm);
5375 comm->load_step += comm->cycl[ddCyclStep];
5376 comm->load_sum += comm->load[0].sum;
5377 comm->load_max += comm->load[0].max;
5378 if (comm->bDynLoadBal)
5380 for (d = 0; d < dd->ndim; d++)
5382 if (comm->load[0].flags & (1<<d))
5384 comm->load_lim[d]++;
5390 comm->load_mdf += comm->load[0].mdf;
5391 comm->load_pme += comm->load[0].pme;
5395 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5399 fprintf(debug, "get_load_distribution finished\n");
5403 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5405 /* Return the relative performance loss on the total run time
5406 * due to the force calculation load imbalance.
5408 if (dd->comm->nload > 0)
5411 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5412 (dd->comm->load_step*dd->nnodes);
5420 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5423 int npp, npme, nnodes, d, limp;
5424 float imbal, pme_f_ratio, lossf, lossp = 0;
5426 gmx_domdec_comm_t *comm;
5429 if (DDMASTER(dd) && comm->nload > 0)
5432 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5433 nnodes = npp + npme;
5434 imbal = comm->load_max*npp/comm->load_sum - 1;
5435 lossf = dd_force_imb_perf_loss(dd);
5436 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5437 fprintf(fplog, "%s", buf);
5438 fprintf(stderr, "\n");
5439 fprintf(stderr, "%s", buf);
5440 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5441 fprintf(fplog, "%s", buf);
5442 fprintf(stderr, "%s", buf);
5444 if (comm->bDynLoadBal)
5446 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5447 for (d = 0; d < dd->ndim; d++)
5449 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5450 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5456 sprintf(buf+strlen(buf), "\n");
5457 fprintf(fplog, "%s", buf);
5458 fprintf(stderr, "%s", buf);
5462 pme_f_ratio = comm->load_pme/comm->load_mdf;
5463 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5466 lossp *= (float)npme/(float)nnodes;
5470 lossp *= (float)npp/(float)nnodes;
5472 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5473 fprintf(fplog, "%s", buf);
5474 fprintf(stderr, "%s", buf);
5475 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5476 fprintf(fplog, "%s", buf);
5477 fprintf(stderr, "%s", buf);
5479 fprintf(fplog, "\n");
5480 fprintf(stderr, "\n");
5482 if (lossf >= DD_PERF_LOSS)
5485 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5486 " in the domain decomposition.\n", lossf*100);
5487 if (!comm->bDynLoadBal)
5489 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5493 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5495 fprintf(fplog, "%s\n", buf);
5496 fprintf(stderr, "%s\n", buf);
5498 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5501 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5502 " had %s work to do than the PP nodes.\n"
5503 " You might want to %s the number of PME nodes\n"
5504 " or %s the cut-off and the grid spacing.\n",
5506 (lossp < 0) ? "less" : "more",
5507 (lossp < 0) ? "decrease" : "increase",
5508 (lossp < 0) ? "decrease" : "increase");
5509 fprintf(fplog, "%s\n", buf);
5510 fprintf(stderr, "%s\n", buf);
5515 static float dd_vol_min(gmx_domdec_t *dd)
5517 return dd->comm->load[0].cvol_min*dd->nnodes;
5520 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5522 return dd->comm->load[0].flags;
5525 static float dd_f_imbal(gmx_domdec_t *dd)
5527 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5530 float dd_pme_f_ratio(gmx_domdec_t *dd)
5532 if (dd->comm->cycl_n[ddCyclPME] > 0)
5534 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5542 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5547 flags = dd_load_flags(dd);
5551 "DD load balancing is limited by minimum cell size in dimension");
5552 for (d = 0; d < dd->ndim; d++)
5556 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5559 fprintf(fplog, "\n");
5561 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5562 if (dd->comm->bDynLoadBal)
5564 fprintf(fplog, " vol min/aver %5.3f%c",
5565 dd_vol_min(dd), flags ? '!' : ' ');
5567 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5568 if (dd->comm->cycl_n[ddCyclPME])
5570 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5572 fprintf(fplog, "\n\n");
5575 static void dd_print_load_verbose(gmx_domdec_t *dd)
5577 if (dd->comm->bDynLoadBal)
5579 fprintf(stderr, "vol %4.2f%c ",
5580 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5582 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5583 if (dd->comm->cycl_n[ddCyclPME])
5585 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5590 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5595 gmx_domdec_root_t *root;
5596 gmx_bool bPartOfGroup = FALSE;
5598 dim = dd->dim[dim_ind];
5599 copy_ivec(loc, loc_c);
5600 for (i = 0; i < dd->nc[dim]; i++)
5603 rank = dd_index(dd->nc, loc_c);
5604 if (rank == dd->rank)
5606 /* This process is part of the group */
5607 bPartOfGroup = TRUE;
5610 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5614 dd->comm->mpi_comm_load[dim_ind] = c_row;
5615 if (dd->comm->eDLB != edlbNO)
5617 if (dd->ci[dim] == dd->master_ci[dim])
5619 /* This is the root process of this row */
5620 snew(dd->comm->root[dim_ind], 1);
5621 root = dd->comm->root[dim_ind];
5622 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5623 snew(root->old_cell_f, dd->nc[dim]+1);
5624 snew(root->bCellMin, dd->nc[dim]);
5627 snew(root->cell_f_max0, dd->nc[dim]);
5628 snew(root->cell_f_min1, dd->nc[dim]);
5629 snew(root->bound_min, dd->nc[dim]);
5630 snew(root->bound_max, dd->nc[dim]);
5632 snew(root->buf_ncd, dd->nc[dim]);
5636 /* This is not a root process, we only need to receive cell_f */
5637 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5640 if (dd->ci[dim] == dd->master_ci[dim])
5642 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5648 static void make_load_communicators(gmx_domdec_t *dd)
5651 int dim0, dim1, i, j;
5656 fprintf(debug, "Making load communicators\n");
5659 snew(dd->comm->load, dd->ndim);
5660 snew(dd->comm->mpi_comm_load, dd->ndim);
5663 make_load_communicator(dd, 0, loc);
5667 for (i = 0; i < dd->nc[dim0]; i++)
5670 make_load_communicator(dd, 1, loc);
5676 for (i = 0; i < dd->nc[dim0]; i++)
5680 for (j = 0; j < dd->nc[dim1]; j++)
5683 make_load_communicator(dd, 2, loc);
5690 fprintf(debug, "Finished making load communicators\n");
5695 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5698 int d, dim, i, j, m;
5701 ivec dd_zp[DD_MAXIZONE];
5702 gmx_domdec_zones_t *zones;
5703 gmx_domdec_ns_ranges_t *izone;
5705 for (d = 0; d < dd->ndim; d++)
5708 copy_ivec(dd->ci, tmp);
5709 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5710 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5711 copy_ivec(dd->ci, tmp);
5712 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5713 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5716 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5719 dd->neighbor[d][1]);
5725 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5727 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5728 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5735 for (i = 0; i < nzonep; i++)
5737 copy_ivec(dd_zp3[i], dd_zp[i]);
5743 for (i = 0; i < nzonep; i++)
5745 copy_ivec(dd_zp2[i], dd_zp[i]);
5751 for (i = 0; i < nzonep; i++)
5753 copy_ivec(dd_zp1[i], dd_zp[i]);
5757 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5762 zones = &dd->comm->zones;
5764 for (i = 0; i < nzone; i++)
5767 clear_ivec(zones->shift[i]);
5768 for (d = 0; d < dd->ndim; d++)
5770 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5775 for (i = 0; i < nzone; i++)
5777 for (d = 0; d < DIM; d++)
5779 s[d] = dd->ci[d] - zones->shift[i][d];
5784 else if (s[d] >= dd->nc[d])
5790 zones->nizone = nzonep;
5791 for (i = 0; i < zones->nizone; i++)
5793 if (dd_zp[i][0] != i)
5795 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5797 izone = &zones->izone[i];
5798 izone->j0 = dd_zp[i][1];
5799 izone->j1 = dd_zp[i][2];
5800 for (dim = 0; dim < DIM; dim++)
5802 if (dd->nc[dim] == 1)
5804 /* All shifts should be allowed */
5805 izone->shift0[dim] = -1;
5806 izone->shift1[dim] = 1;
5811 izone->shift0[d] = 0;
5812 izone->shift1[d] = 0;
5813 for(j=izone->j0; j<izone->j1; j++) {
5814 if (dd->shift[j][d] > dd->shift[i][d])
5815 izone->shift0[d] = -1;
5816 if (dd->shift[j][d] < dd->shift[i][d])
5817 izone->shift1[d] = 1;
5823 /* Assume the shift are not more than 1 cell */
5824 izone->shift0[dim] = 1;
5825 izone->shift1[dim] = -1;
5826 for (j = izone->j0; j < izone->j1; j++)
5828 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5829 if (shift_diff < izone->shift0[dim])
5831 izone->shift0[dim] = shift_diff;
5833 if (shift_diff > izone->shift1[dim])
5835 izone->shift1[dim] = shift_diff;
5842 if (dd->comm->eDLB != edlbNO)
5844 snew(dd->comm->root, dd->ndim);
5847 if (dd->comm->bRecordLoad)
5849 make_load_communicators(dd);
5853 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5856 gmx_domdec_comm_t *comm;
5867 if (comm->bCartesianPP)
5869 /* Set up cartesian communication for the particle-particle part */
5872 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5873 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5876 for (i = 0; i < DIM; i++)
5880 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5882 /* We overwrite the old communicator with the new cartesian one */
5883 cr->mpi_comm_mygroup = comm_cart;
5886 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5887 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5889 if (comm->bCartesianPP_PME)
5891 /* Since we want to use the original cartesian setup for sim,
5892 * and not the one after split, we need to make an index.
5894 snew(comm->ddindex2ddnodeid, dd->nnodes);
5895 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5896 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5897 /* Get the rank of the DD master,
5898 * above we made sure that the master node is a PP node.
5908 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5910 else if (comm->bCartesianPP)
5912 if (cr->npmenodes == 0)
5914 /* The PP communicator is also
5915 * the communicator for this simulation
5917 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5919 cr->nodeid = dd->rank;
5921 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5923 /* We need to make an index to go from the coordinates
5924 * to the nodeid of this simulation.
5926 snew(comm->ddindex2simnodeid, dd->nnodes);
5927 snew(buf, dd->nnodes);
5928 if (cr->duty & DUTY_PP)
5930 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5932 /* Communicate the ddindex to simulation nodeid index */
5933 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5934 cr->mpi_comm_mysim);
5937 /* Determine the master coordinates and rank.
5938 * The DD master should be the same node as the master of this sim.
5940 for (i = 0; i < dd->nnodes; i++)
5942 if (comm->ddindex2simnodeid[i] == 0)
5944 ddindex2xyz(dd->nc, i, dd->master_ci);
5945 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5950 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5955 /* No Cartesian communicators */
5956 /* We use the rank in dd->comm->all as DD index */
5957 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5958 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5960 clear_ivec(dd->master_ci);
5967 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5968 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5973 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5974 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5978 static void receive_ddindex2simnodeid(t_commrec *cr)
5982 gmx_domdec_comm_t *comm;
5989 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5991 snew(comm->ddindex2simnodeid, dd->nnodes);
5992 snew(buf, dd->nnodes);
5993 if (cr->duty & DUTY_PP)
5995 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5998 /* Communicate the ddindex to simulation nodeid index */
5999 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
6000 cr->mpi_comm_mysim);
6007 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
6008 int ncg, int natoms)
6010 gmx_domdec_master_t *ma;
6015 snew(ma->ncg, dd->nnodes);
6016 snew(ma->index, dd->nnodes+1);
6018 snew(ma->nat, dd->nnodes);
6019 snew(ma->ibuf, dd->nnodes*2);
6020 snew(ma->cell_x, DIM);
6021 for (i = 0; i < DIM; i++)
6023 snew(ma->cell_x[i], dd->nc[i]+1);
6026 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6032 snew(ma->vbuf, natoms);
6038 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
6042 gmx_domdec_comm_t *comm;
6053 if (comm->bCartesianPP)
6055 for (i = 1; i < DIM; i++)
6057 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6059 if (bDiv[YY] || bDiv[ZZ])
6061 comm->bCartesianPP_PME = TRUE;
6062 /* If we have 2D PME decomposition, which is always in x+y,
6063 * we stack the PME only nodes in z.
6064 * Otherwise we choose the direction that provides the thinnest slab
6065 * of PME only nodes as this will have the least effect
6066 * on the PP communication.
6067 * But for the PME communication the opposite might be better.
6069 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6071 dd->nc[YY] > dd->nc[ZZ]))
6073 comm->cartpmedim = ZZ;
6077 comm->cartpmedim = YY;
6079 comm->ntot[comm->cartpmedim]
6080 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6084 fprintf(fplog, "#pmenodes (%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]);
6086 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6091 if (comm->bCartesianPP_PME)
6095 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]);
6098 for (i = 0; i < DIM; i++)
6102 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6105 MPI_Comm_rank(comm_cart, &rank);
6106 if (MASTERNODE(cr) && rank != 0)
6108 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6111 /* With this assigment we loose the link to the original communicator
6112 * which will usually be MPI_COMM_WORLD, unless have multisim.
6114 cr->mpi_comm_mysim = comm_cart;
6115 cr->sim_nodeid = rank;
6117 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6121 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6122 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6125 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6129 if (cr->npmenodes == 0 ||
6130 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6132 cr->duty = DUTY_PME;
6135 /* Split the sim communicator into PP and PME only nodes */
6136 MPI_Comm_split(cr->mpi_comm_mysim,
6138 dd_index(comm->ntot, dd->ci),
6139 &cr->mpi_comm_mygroup);
6143 switch (dd_node_order)
6148 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6151 case ddnoINTERLEAVE:
6152 /* Interleave the PP-only and PME-only nodes,
6153 * as on clusters with dual-core machines this will double
6154 * the communication bandwidth of the PME processes
6155 * and thus speed up the PP <-> PME and inter PME communication.
6159 fprintf(fplog, "Interleaving PP and PME nodes\n");
6161 comm->pmenodes = dd_pmenodes(cr);
6166 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6169 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6171 cr->duty = DUTY_PME;
6178 /* Split the sim communicator into PP and PME only nodes */
6179 MPI_Comm_split(cr->mpi_comm_mysim,
6182 &cr->mpi_comm_mygroup);
6183 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6189 fprintf(fplog, "This is a %s only node\n\n",
6190 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6194 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6197 gmx_domdec_comm_t *comm;
6203 copy_ivec(dd->nc, comm->ntot);
6205 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6206 comm->bCartesianPP_PME = FALSE;
6208 /* Reorder the nodes by default. This might change the MPI ranks.
6209 * Real reordering is only supported on very few architectures,
6210 * Blue Gene is one of them.
6212 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6214 if (cr->npmenodes > 0)
6216 /* Split the communicator into a PP and PME part */
6217 split_communicator(fplog, cr, dd_node_order, CartReorder);
6218 if (comm->bCartesianPP_PME)
6220 /* We (possibly) reordered the nodes in split_communicator,
6221 * so it is no longer required in make_pp_communicator.
6223 CartReorder = FALSE;
6228 /* All nodes do PP and PME */
6230 /* We do not require separate communicators */
6231 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6235 if (cr->duty & DUTY_PP)
6237 /* Copy or make a new PP communicator */
6238 make_pp_communicator(fplog, cr, CartReorder);
6242 receive_ddindex2simnodeid(cr);
6245 if (!(cr->duty & DUTY_PME))
6247 /* Set up the commnuication to our PME node */
6248 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6249 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6252 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6253 dd->pme_nodeid, dd->pme_receive_vir_ener);
6258 dd->pme_nodeid = -1;
6263 dd->ma = init_gmx_domdec_master_t(dd,
6265 comm->cgs_gl.index[comm->cgs_gl.nr]);
6269 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6271 real *slb_frac, tot;
6276 if (nc > 1 && size_string != NULL)
6280 fprintf(fplog, "Using static load balancing for the %s direction\n",
6285 for (i = 0; i < nc; i++)
6288 sscanf(size_string, "%lf%n", &dbl, &n);
6291 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6300 fprintf(fplog, "Relative cell sizes:");
6302 for (i = 0; i < nc; i++)
6307 fprintf(fplog, " %5.3f", slb_frac[i]);
6312 fprintf(fplog, "\n");
6319 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6322 gmx_mtop_ilistloop_t iloop;
6326 iloop = gmx_mtop_ilistloop_init(mtop);
6327 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6329 for (ftype = 0; ftype < F_NRE; ftype++)
6331 if ((interaction_function[ftype].flags & IF_BOND) &&
6334 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6342 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6348 val = getenv(env_var);
6351 if (sscanf(val, "%d", &nst) <= 0)
6357 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6365 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6369 fprintf(stderr, "\n%s\n", warn_string);
6373 fprintf(fplog, "\n%s\n", warn_string);
6377 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6378 t_inputrec *ir, FILE *fplog)
6380 if (ir->ePBC == epbcSCREW &&
6381 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6383 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6386 if (ir->ns_type == ensSIMPLE)
6388 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6391 if (ir->nstlist == 0)
6393 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6396 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6398 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6402 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6407 r = ddbox->box_size[XX];
6408 for (di = 0; di < dd->ndim; di++)
6411 /* Check using the initial average cell size */
6412 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6418 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6419 const char *dlb_opt, gmx_bool bRecordLoad,
6420 unsigned long Flags, t_inputrec *ir)
6428 case 'a': eDLB = edlbAUTO; break;
6429 case 'n': eDLB = edlbNO; break;
6430 case 'y': eDLB = edlbYES; break;
6431 default: gmx_incons("Unknown dlb_opt");
6434 if (Flags & MD_RERUN)
6439 if (!EI_DYNAMICS(ir->eI))
6441 if (eDLB == edlbYES)
6443 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6444 dd_warning(cr, fplog, buf);
6452 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6457 if (Flags & MD_REPRODUCIBLE)
6464 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6468 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6471 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6479 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6484 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6486 /* Decomposition order z,y,x */
6489 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6491 for (dim = DIM-1; dim >= 0; dim--)
6493 if (dd->nc[dim] > 1)
6495 dd->dim[dd->ndim++] = dim;
6501 /* Decomposition order x,y,z */
6502 for (dim = 0; dim < DIM; dim++)
6504 if (dd->nc[dim] > 1)
6506 dd->dim[dd->ndim++] = dim;
6512 static gmx_domdec_comm_t *init_dd_comm()
6514 gmx_domdec_comm_t *comm;
6518 snew(comm->cggl_flag, DIM*2);
6519 snew(comm->cgcm_state, DIM*2);
6520 for (i = 0; i < DIM*2; i++)
6522 comm->cggl_flag_nalloc[i] = 0;
6523 comm->cgcm_state_nalloc[i] = 0;
6526 comm->nalloc_int = 0;
6527 comm->buf_int = NULL;
6529 vec_rvec_init(&comm->vbuf);
6531 comm->n_load_have = 0;
6532 comm->n_load_collect = 0;
6534 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6536 comm->sum_nat[i] = 0;
6540 comm->load_step = 0;
6543 clear_ivec(comm->load_lim);
6550 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6551 unsigned long Flags,
6553 real comm_distance_min, real rconstr,
6554 const char *dlb_opt, real dlb_scale,
6555 const char *sizex, const char *sizey, const char *sizez,
6556 gmx_mtop_t *mtop, t_inputrec *ir,
6557 matrix box, rvec *x,
6559 int *npme_x, int *npme_y)
6562 gmx_domdec_comm_t *comm;
6565 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6572 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6577 dd->comm = init_dd_comm();
6579 snew(comm->cggl_flag, DIM*2);
6580 snew(comm->cgcm_state, DIM*2);
6582 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6583 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6585 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6586 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6587 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6588 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6589 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6590 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6591 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6592 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6594 dd->pme_recv_f_alloc = 0;
6595 dd->pme_recv_f_buf = NULL;
6597 if (dd->bSendRecv2 && fplog)
6599 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");
6605 fprintf(fplog, "Will load balance based on FLOP count\n");
6607 if (comm->eFlop > 1)
6609 srand(1+cr->nodeid);
6611 comm->bRecordLoad = TRUE;
6615 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6619 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6621 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6624 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6626 dd->bGridJump = comm->bDynLoadBal;
6627 comm->bPMELoadBalDLBLimits = FALSE;
6629 if (comm->nstSortCG)
6633 if (comm->nstSortCG == 1)
6635 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6639 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6643 snew(comm->sort, 1);
6649 fprintf(fplog, "Will not sort the charge groups\n");
6653 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6655 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6656 if (comm->bInterCGBondeds)
6658 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6662 comm->bInterCGMultiBody = FALSE;
6665 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6666 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6668 if (ir->rlistlong == 0)
6670 /* Set the cut-off to some very large value,
6671 * so we don't need if statements everywhere in the code.
6672 * We use sqrt, since the cut-off is squared in some places.
6674 comm->cutoff = GMX_CUTOFF_INF;
6678 comm->cutoff = ir->rlistlong;
6680 comm->cutoff_mbody = 0;
6682 comm->cellsize_limit = 0;
6683 comm->bBondComm = FALSE;
6685 if (comm->bInterCGBondeds)
6687 if (comm_distance_min > 0)
6689 comm->cutoff_mbody = comm_distance_min;
6690 if (Flags & MD_DDBONDCOMM)
6692 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6696 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6698 r_bonded_limit = comm->cutoff_mbody;
6700 else if (ir->bPeriodicMols)
6702 /* Can not easily determine the required cut-off */
6703 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");
6704 comm->cutoff_mbody = comm->cutoff/2;
6705 r_bonded_limit = comm->cutoff_mbody;
6711 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6712 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6714 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6715 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6717 /* We use an initial margin of 10% for the minimum cell size,
6718 * except when we are just below the non-bonded cut-off.
6720 if (Flags & MD_DDBONDCOMM)
6722 if (max(r_2b, r_mb) > comm->cutoff)
6724 r_bonded = max(r_2b, r_mb);
6725 r_bonded_limit = 1.1*r_bonded;
6726 comm->bBondComm = TRUE;
6731 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6733 /* We determine cutoff_mbody later */
6737 /* No special bonded communication,
6738 * simply increase the DD cut-off.
6740 r_bonded_limit = 1.1*max(r_2b, r_mb);
6741 comm->cutoff_mbody = r_bonded_limit;
6742 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6745 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6749 "Minimum cell size due to bonded interactions: %.3f nm\n",
6750 comm->cellsize_limit);
6754 if (dd->bInterCGcons && rconstr <= 0)
6756 /* There is a cell size limit due to the constraints (P-LINCS) */
6757 rconstr = constr_r_max(fplog, mtop, ir);
6761 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6763 if (rconstr > comm->cellsize_limit)
6765 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6769 else if (rconstr > 0 && fplog)
6771 /* Here we do not check for dd->bInterCGcons,
6772 * because one can also set a cell size limit for virtual sites only
6773 * and at this point we don't know yet if there are intercg v-sites.
6776 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6779 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6781 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6785 copy_ivec(nc, dd->nc);
6786 set_dd_dim(fplog, dd);
6787 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6789 if (cr->npmenodes == -1)
6793 acs = average_cellsize_min(dd, ddbox);
6794 if (acs < comm->cellsize_limit)
6798 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6800 gmx_fatal_collective(FARGS, cr, NULL,
6801 "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",
6802 acs, comm->cellsize_limit);
6807 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6809 /* We need to choose the optimal DD grid and possibly PME nodes */
6810 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6811 comm->eDLB != edlbNO, dlb_scale,
6812 comm->cellsize_limit, comm->cutoff,
6813 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6815 if (dd->nc[XX] == 0)
6817 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6818 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6819 !bC ? "-rdd" : "-rcon",
6820 comm->eDLB != edlbNO ? " or -dds" : "",
6821 bC ? " or your LINCS settings" : "");
6823 gmx_fatal_collective(FARGS, cr, NULL,
6824 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6826 "Look in the log file for details on the domain decomposition",
6827 cr->nnodes-cr->npmenodes, limit, buf);
6829 set_dd_dim(fplog, dd);
6835 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6836 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6839 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6840 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6842 gmx_fatal_collective(FARGS, cr, NULL,
6843 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6844 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6846 if (cr->npmenodes > dd->nnodes)
6848 gmx_fatal_collective(FARGS, cr, NULL,
6849 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6851 if (cr->npmenodes > 0)
6853 comm->npmenodes = cr->npmenodes;
6857 comm->npmenodes = dd->nnodes;
6860 if (EEL_PME(ir->coulombtype))
6862 /* The following choices should match those
6863 * in comm_cost_est in domdec_setup.c.
6864 * Note that here the checks have to take into account
6865 * that the decomposition might occur in a different order than xyz
6866 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6867 * in which case they will not match those in comm_cost_est,
6868 * but since that is mainly for testing purposes that's fine.
6870 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6871 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6872 getenv("GMX_PMEONEDD") == NULL)
6874 comm->npmedecompdim = 2;
6875 comm->npmenodes_x = dd->nc[XX];
6876 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6880 /* In case nc is 1 in both x and y we could still choose to
6881 * decompose pme in y instead of x, but we use x for simplicity.
6883 comm->npmedecompdim = 1;
6884 if (dd->dim[0] == YY)
6886 comm->npmenodes_x = 1;
6887 comm->npmenodes_y = comm->npmenodes;
6891 comm->npmenodes_x = comm->npmenodes;
6892 comm->npmenodes_y = 1;
6897 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6898 comm->npmenodes_x, comm->npmenodes_y, 1);
6903 comm->npmedecompdim = 0;
6904 comm->npmenodes_x = 0;
6905 comm->npmenodes_y = 0;
6908 /* Technically we don't need both of these,
6909 * but it simplifies code not having to recalculate it.
6911 *npme_x = comm->npmenodes_x;
6912 *npme_y = comm->npmenodes_y;
6914 snew(comm->slb_frac, DIM);
6915 if (comm->eDLB == edlbNO)
6917 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6918 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6919 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6922 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6924 if (comm->bBondComm || comm->eDLB != edlbNO)
6926 /* Set the bonded communication distance to halfway
6927 * the minimum and the maximum,
6928 * since the extra communication cost is nearly zero.
6930 acs = average_cellsize_min(dd, ddbox);
6931 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6932 if (comm->eDLB != edlbNO)
6934 /* Check if this does not limit the scaling */
6935 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6937 if (!comm->bBondComm)
6939 /* Without bBondComm do not go beyond the n.b. cut-off */
6940 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6941 if (comm->cellsize_limit >= comm->cutoff)
6943 /* We don't loose a lot of efficieny
6944 * when increasing it to the n.b. cut-off.
6945 * It can even be slightly faster, because we need
6946 * less checks for the communication setup.
6948 comm->cutoff_mbody = comm->cutoff;
6951 /* Check if we did not end up below our original limit */
6952 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6954 if (comm->cutoff_mbody > comm->cellsize_limit)
6956 comm->cellsize_limit = comm->cutoff_mbody;
6959 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6964 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6965 "cellsize limit %f\n",
6966 comm->bBondComm, comm->cellsize_limit);
6971 check_dd_restrictions(cr, dd, ir, fplog);
6974 comm->partition_step = INT_MIN;
6977 clear_dd_cycle_counts(dd);
6982 static void set_dlb_limits(gmx_domdec_t *dd)
6987 for (d = 0; d < dd->ndim; d++)
6989 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6990 dd->comm->cellsize_min[dd->dim[d]] =
6991 dd->comm->cellsize_min_dlb[dd->dim[d]];
6996 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6999 gmx_domdec_comm_t *comm;
7009 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);
7012 cellsize_min = comm->cellsize_min[dd->dim[0]];
7013 for (d = 1; d < dd->ndim; d++)
7015 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
7018 if (cellsize_min < comm->cellsize_limit*1.05)
7020 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");
7022 /* Change DLB from "auto" to "no". */
7023 comm->eDLB = edlbNO;
7028 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
7029 comm->bDynLoadBal = TRUE;
7030 dd->bGridJump = TRUE;
7034 /* We can set the required cell size info here,
7035 * so we do not need to communicate this.
7036 * The grid is completely uniform.
7038 for (d = 0; d < dd->ndim; d++)
7042 comm->load[d].sum_m = comm->load[d].sum;
7044 nc = dd->nc[dd->dim[d]];
7045 for (i = 0; i < nc; i++)
7047 comm->root[d]->cell_f[i] = i/(real)nc;
7050 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7051 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7054 comm->root[d]->cell_f[nc] = 1.0;
7059 static char *init_bLocalCG(gmx_mtop_t *mtop)
7064 ncg = ncg_mtop(mtop);
7065 snew(bLocalCG, ncg);
7066 for (cg = 0; cg < ncg; cg++)
7068 bLocalCG[cg] = FALSE;
7074 void dd_init_bondeds(FILE *fplog,
7075 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7076 gmx_vsite_t *vsite, gmx_constr_t constr,
7077 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7079 gmx_domdec_comm_t *comm;
7083 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7087 if (comm->bBondComm)
7089 /* Communicate atoms beyond the cut-off for bonded interactions */
7092 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7094 comm->bLocalCG = init_bLocalCG(mtop);
7098 /* Only communicate atoms based on cut-off */
7099 comm->cglink = NULL;
7100 comm->bLocalCG = NULL;
7104 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7106 gmx_bool bDynLoadBal, real dlb_scale,
7109 gmx_domdec_comm_t *comm;
7124 fprintf(fplog, "The maximum number of communication pulses is:");
7125 for (d = 0; d < dd->ndim; d++)
7127 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7129 fprintf(fplog, "\n");
7130 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7131 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7132 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7133 for (d = 0; d < DIM; d++)
7137 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7144 comm->cellsize_min_dlb[d]/
7145 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7147 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7150 fprintf(fplog, "\n");
7154 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7155 fprintf(fplog, "The initial number of communication pulses is:");
7156 for (d = 0; d < dd->ndim; d++)
7158 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7160 fprintf(fplog, "\n");
7161 fprintf(fplog, "The initial domain decomposition cell size is:");
7162 for (d = 0; d < DIM; d++)
7166 fprintf(fplog, " %c %.2f nm",
7167 dim2char(d), dd->comm->cellsize_min[d]);
7170 fprintf(fplog, "\n\n");
7173 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7175 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7176 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7177 "non-bonded interactions", "", comm->cutoff);
7181 limit = dd->comm->cellsize_limit;
7185 if (dynamic_dd_box(ddbox, ir))
7187 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7189 limit = dd->comm->cellsize_min[XX];
7190 for (d = 1; d < DIM; d++)
7192 limit = min(limit, dd->comm->cellsize_min[d]);
7196 if (comm->bInterCGBondeds)
7198 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7199 "two-body bonded interactions", "(-rdd)",
7200 max(comm->cutoff, comm->cutoff_mbody));
7201 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7202 "multi-body bonded interactions", "(-rdd)",
7203 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7207 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7208 "virtual site constructions", "(-rcon)", limit);
7210 if (dd->constraint_comm)
7212 sprintf(buf, "atoms separated by up to %d constraints",
7214 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7215 buf, "(-rcon)", limit);
7217 fprintf(fplog, "\n");
7223 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7225 const t_inputrec *ir,
7226 const gmx_ddbox_t *ddbox)
7228 gmx_domdec_comm_t *comm;
7229 int d, dim, npulse, npulse_d_max, npulse_d;
7234 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7236 /* Determine the maximum number of comm. pulses in one dimension */
7238 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7240 /* Determine the maximum required number of grid pulses */
7241 if (comm->cellsize_limit >= comm->cutoff)
7243 /* Only a single pulse is required */
7246 else if (!bNoCutOff && comm->cellsize_limit > 0)
7248 /* We round down slightly here to avoid overhead due to the latency
7249 * of extra communication calls when the cut-off
7250 * would be only slightly longer than the cell size.
7251 * Later cellsize_limit is redetermined,
7252 * so we can not miss interactions due to this rounding.
7254 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7258 /* There is no cell size limit */
7259 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7262 if (!bNoCutOff && npulse > 1)
7264 /* See if we can do with less pulses, based on dlb_scale */
7266 for (d = 0; d < dd->ndim; d++)
7269 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7270 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7271 npulse_d_max = max(npulse_d_max, npulse_d);
7273 npulse = min(npulse, npulse_d_max);
7276 /* This env var can override npulse */
7277 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7284 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7285 for (d = 0; d < dd->ndim; d++)
7287 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7288 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7289 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7290 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7291 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7293 comm->bVacDLBNoLimit = FALSE;
7297 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7298 if (!comm->bVacDLBNoLimit)
7300 comm->cellsize_limit = max(comm->cellsize_limit,
7301 comm->cutoff/comm->maxpulse);
7303 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7304 /* Set the minimum cell size for each DD dimension */
7305 for (d = 0; d < dd->ndim; d++)
7307 if (comm->bVacDLBNoLimit ||
7308 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7310 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7314 comm->cellsize_min_dlb[dd->dim[d]] =
7315 comm->cutoff/comm->cd[d].np_dlb;
7318 if (comm->cutoff_mbody <= 0)
7320 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7322 if (comm->bDynLoadBal)
7328 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7330 /* If each molecule is a single charge group
7331 * or we use domain decomposition for each periodic dimension,
7332 * we do not need to take pbc into account for the bonded interactions.
7334 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7337 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7340 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7341 t_inputrec *ir, t_forcerec *fr,
7344 gmx_domdec_comm_t *comm;
7350 /* Initialize the thread data.
7351 * This can not be done in init_domain_decomposition,
7352 * as the numbers of threads is determined later.
7354 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7357 snew(comm->dth, comm->nth);
7360 if (EEL_PME(ir->coulombtype))
7362 init_ddpme(dd, &comm->ddpme[0], 0);
7363 if (comm->npmedecompdim >= 2)
7365 init_ddpme(dd, &comm->ddpme[1], 1);
7370 comm->npmenodes = 0;
7371 if (dd->pme_nodeid >= 0)
7373 gmx_fatal_collective(FARGS, NULL, dd,
7374 "Can not have separate PME nodes without PME electrostatics");
7380 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7382 if (comm->eDLB != edlbNO)
7384 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7387 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7388 if (comm->eDLB == edlbAUTO)
7392 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7394 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7397 if (ir->ePBC == epbcNONE)
7399 vol_frac = 1 - 1/(double)dd->nnodes;
7404 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7408 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7410 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7412 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7415 static gmx_bool test_dd_cutoff(t_commrec *cr,
7416 t_state *state, t_inputrec *ir,
7427 set_ddbox(dd, FALSE, cr, ir, state->box,
7428 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7432 for (d = 0; d < dd->ndim; d++)
7436 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7437 if (dynamic_dd_box(&ddbox, ir))
7439 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7442 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7444 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7445 dd->comm->cd[d].np_dlb > 0)
7447 if (np > dd->comm->cd[d].np_dlb)
7452 /* If a current local cell size is smaller than the requested
7453 * cut-off, we could still fix it, but this gets very complicated.
7454 * Without fixing here, we might actually need more checks.
7456 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7463 if (dd->comm->eDLB != edlbNO)
7465 /* If DLB is not active yet, we don't need to check the grid jumps.
7466 * Actually we shouldn't, because then the grid jump data is not set.
7468 if (dd->comm->bDynLoadBal &&
7469 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7474 gmx_sumi(1, &LocallyLimited, cr);
7476 if (LocallyLimited > 0)
7485 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7488 gmx_bool bCutoffAllowed;
7490 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7494 cr->dd->comm->cutoff = cutoff_req;
7497 return bCutoffAllowed;
7500 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7502 gmx_domdec_comm_t *comm;
7504 comm = cr->dd->comm;
7506 /* Turn on the DLB limiting (might have been on already) */
7507 comm->bPMELoadBalDLBLimits = TRUE;
7509 /* Change the cut-off limit */
7510 comm->PMELoadBal_max_cutoff = comm->cutoff;
7513 static void merge_cg_buffers(int ncell,
7514 gmx_domdec_comm_dim_t *cd, int pulse,
7516 int *index_gl, int *recv_i,
7517 rvec *cg_cm, rvec *recv_vr,
7519 cginfo_mb_t *cginfo_mb, int *cginfo)
7521 gmx_domdec_ind_t *ind, *ind_p;
7522 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7523 int shift, shift_at;
7525 ind = &cd->ind[pulse];
7527 /* First correct the already stored data */
7528 shift = ind->nrecv[ncell];
7529 for (cell = ncell-1; cell >= 0; cell--)
7531 shift -= ind->nrecv[cell];
7534 /* Move the cg's present from previous grid pulses */
7535 cg0 = ncg_cell[ncell+cell];
7536 cg1 = ncg_cell[ncell+cell+1];
7537 cgindex[cg1+shift] = cgindex[cg1];
7538 for (cg = cg1-1; cg >= cg0; cg--)
7540 index_gl[cg+shift] = index_gl[cg];
7541 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7542 cgindex[cg+shift] = cgindex[cg];
7543 cginfo[cg+shift] = cginfo[cg];
7545 /* Correct the already stored send indices for the shift */
7546 for (p = 1; p <= pulse; p++)
7548 ind_p = &cd->ind[p];
7550 for (c = 0; c < cell; c++)
7552 cg0 += ind_p->nsend[c];
7554 cg1 = cg0 + ind_p->nsend[cell];
7555 for (cg = cg0; cg < cg1; cg++)
7557 ind_p->index[cg] += shift;
7563 /* Merge in the communicated buffers */
7567 for (cell = 0; cell < ncell; cell++)
7569 cg1 = ncg_cell[ncell+cell+1] + shift;
7572 /* Correct the old cg indices */
7573 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7575 cgindex[cg+1] += shift_at;
7578 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7580 /* Copy this charge group from the buffer */
7581 index_gl[cg1] = recv_i[cg0];
7582 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7583 /* Add it to the cgindex */
7584 cg_gl = index_gl[cg1];
7585 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7586 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7587 cgindex[cg1+1] = cgindex[cg1] + nat;
7592 shift += ind->nrecv[cell];
7593 ncg_cell[ncell+cell+1] = cg1;
7597 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7598 int nzone, int cg0, const int *cgindex)
7602 /* Store the atom block boundaries for easy copying of communication buffers
7605 for (zone = 0; zone < nzone; zone++)
7607 for (p = 0; p < cd->np; p++)
7609 cd->ind[p].cell2at0[zone] = cgindex[cg];
7610 cg += cd->ind[p].nrecv[zone];
7611 cd->ind[p].cell2at1[zone] = cgindex[cg];
7616 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7622 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7624 if (!bLocalCG[link->a[i]])
7633 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7635 real c[DIM][4]; /* the corners for the non-bonded communication */
7636 real cr0; /* corner for rounding */
7637 real cr1[4]; /* corners for rounding */
7638 real bc[DIM]; /* corners for bounded communication */
7639 real bcr1; /* corner for rounding for bonded communication */
7642 /* Determine the corners of the domain(s) we are communicating with */
7644 set_dd_corners(const gmx_domdec_t *dd,
7645 int dim0, int dim1, int dim2,
7649 const gmx_domdec_comm_t *comm;
7650 const gmx_domdec_zones_t *zones;
7655 zones = &comm->zones;
7657 /* Keep the compiler happy */
7661 /* The first dimension is equal for all cells */
7662 c->c[0][0] = comm->cell_x0[dim0];
7665 c->bc[0] = c->c[0][0];
7670 /* This cell row is only seen from the first row */
7671 c->c[1][0] = comm->cell_x0[dim1];
7672 /* All rows can see this row */
7673 c->c[1][1] = comm->cell_x0[dim1];
7676 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7679 /* For the multi-body distance we need the maximum */
7680 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7683 /* Set the upper-right corner for rounding */
7684 c->cr0 = comm->cell_x1[dim0];
7689 for (j = 0; j < 4; j++)
7691 c->c[2][j] = comm->cell_x0[dim2];
7695 /* Use the maximum of the i-cells that see a j-cell */
7696 for (i = 0; i < zones->nizone; i++)
7698 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7704 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7710 /* For the multi-body distance we need the maximum */
7711 c->bc[2] = comm->cell_x0[dim2];
7712 for (i = 0; i < 2; i++)
7714 for (j = 0; j < 2; j++)
7716 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7722 /* Set the upper-right corner for rounding */
7723 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7724 * Only cell (0,0,0) can see cell 7 (1,1,1)
7726 c->cr1[0] = comm->cell_x1[dim1];
7727 c->cr1[3] = comm->cell_x1[dim1];
7730 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7733 /* For the multi-body distance we need the maximum */
7734 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7741 /* Determine which cg's we need to send in this pulse from this zone */
7743 get_zone_pulse_cgs(gmx_domdec_t *dd,
7744 int zonei, int zone,
7746 const int *index_gl,
7748 int dim, int dim_ind,
7749 int dim0, int dim1, int dim2,
7750 real r_comm2, real r_bcomm2,
7754 real skew_fac2_d, real skew_fac_01,
7755 rvec *v_d, rvec *v_0, rvec *v_1,
7756 const dd_corners_t *c,
7758 gmx_bool bDistBonded,
7764 gmx_domdec_ind_t *ind,
7765 int **ibuf, int *ibuf_nalloc,
7771 gmx_domdec_comm_t *comm;
7773 gmx_bool bDistMB_pulse;
7775 real r2, rb2, r, tric_sh;
7778 int nsend_z, nsend, nat;
7782 bScrew = (dd->bScrewPBC && dim == XX);
7784 bDistMB_pulse = (bDistMB && bDistBonded);
7790 for (cg = cg0; cg < cg1; cg++)
7794 if (tric_dist[dim_ind] == 0)
7796 /* Rectangular direction, easy */
7797 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7804 r = cg_cm[cg][dim] - c->bc[dim_ind];
7810 /* Rounding gives at most a 16% reduction
7811 * in communicated atoms
7813 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7815 r = cg_cm[cg][dim0] - c->cr0;
7816 /* This is the first dimension, so always r >= 0 */
7823 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7825 r = cg_cm[cg][dim1] - c->cr1[zone];
7832 r = cg_cm[cg][dim1] - c->bcr1;
7842 /* Triclinic direction, more complicated */
7845 /* Rounding, conservative as the skew_fac multiplication
7846 * will slightly underestimate the distance.
7848 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7850 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7851 for (i = dim0+1; i < DIM; i++)
7853 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7855 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7858 rb[dim0] = rn[dim0];
7861 /* Take care that the cell planes along dim0 might not
7862 * be orthogonal to those along dim1 and dim2.
7864 for (i = 1; i <= dim_ind; i++)
7867 if (normal[dim0][dimd] > 0)
7869 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7872 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7877 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7879 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7881 for (i = dim1+1; i < DIM; i++)
7883 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7885 rn[dim1] += tric_sh;
7888 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7889 /* Take care of coupling of the distances
7890 * to the planes along dim0 and dim1 through dim2.
7892 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7893 /* Take care that the cell planes along dim1
7894 * might not be orthogonal to that along dim2.
7896 if (normal[dim1][dim2] > 0)
7898 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7904 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7907 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7908 /* Take care of coupling of the distances
7909 * to the planes along dim0 and dim1 through dim2.
7911 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7912 /* Take care that the cell planes along dim1
7913 * might not be orthogonal to that along dim2.
7915 if (normal[dim1][dim2] > 0)
7917 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7922 /* The distance along the communication direction */
7923 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7925 for (i = dim+1; i < DIM; i++)
7927 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7932 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7933 /* Take care of coupling of the distances
7934 * to the planes along dim0 and dim1 through dim2.
7936 if (dim_ind == 1 && zonei == 1)
7938 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7944 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7947 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7948 /* Take care of coupling of the distances
7949 * to the planes along dim0 and dim1 through dim2.
7951 if (dim_ind == 1 && zonei == 1)
7953 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7961 ((bDistMB && rb2 < r_bcomm2) ||
7962 (bDist2B && r2 < r_bcomm2)) &&
7964 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7965 missing_link(comm->cglink, index_gl[cg],
7968 /* Make an index to the local charge groups */
7969 if (nsend+1 > ind->nalloc)
7971 ind->nalloc = over_alloc_large(nsend+1);
7972 srenew(ind->index, ind->nalloc);
7974 if (nsend+1 > *ibuf_nalloc)
7976 *ibuf_nalloc = over_alloc_large(nsend+1);
7977 srenew(*ibuf, *ibuf_nalloc);
7979 ind->index[nsend] = cg;
7980 (*ibuf)[nsend] = index_gl[cg];
7982 vec_rvec_check_alloc(vbuf, nsend+1);
7984 if (dd->ci[dim] == 0)
7986 /* Correct cg_cm for pbc */
7987 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7990 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7991 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7996 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7999 nat += cgindex[cg+1] - cgindex[cg];
8005 *nsend_z_ptr = nsend_z;
8008 static void setup_dd_communication(gmx_domdec_t *dd,
8009 matrix box, gmx_ddbox_t *ddbox,
8010 t_forcerec *fr, t_state *state, rvec **f)
8012 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
8013 int nzone, nzone_send, zone, zonei, cg0, cg1;
8014 int c, i, j, cg, cg_gl, nrcg;
8015 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
8016 gmx_domdec_comm_t *comm;
8017 gmx_domdec_zones_t *zones;
8018 gmx_domdec_comm_dim_t *cd;
8019 gmx_domdec_ind_t *ind;
8020 cginfo_mb_t *cginfo_mb;
8021 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
8022 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
8023 dd_corners_t corners;
8025 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
8026 real skew_fac2_d, skew_fac_01;
8033 fprintf(debug, "Setting up DD communication\n");
8038 switch (fr->cutoff_scheme)
8047 gmx_incons("unimplemented");
8051 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8053 dim = dd->dim[dim_ind];
8055 /* Check if we need to use triclinic distances */
8056 tric_dist[dim_ind] = 0;
8057 for (i = 0; i <= dim_ind; i++)
8059 if (ddbox->tric_dir[dd->dim[i]])
8061 tric_dist[dim_ind] = 1;
8066 bBondComm = comm->bBondComm;
8068 /* Do we need to determine extra distances for multi-body bondeds? */
8069 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8071 /* Do we need to determine extra distances for only two-body bondeds? */
8072 bDist2B = (bBondComm && !bDistMB);
8074 r_comm2 = sqr(comm->cutoff);
8075 r_bcomm2 = sqr(comm->cutoff_mbody);
8079 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8082 zones = &comm->zones;
8085 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8086 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8088 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8090 /* Triclinic stuff */
8091 normal = ddbox->normal;
8095 v_0 = ddbox->v[dim0];
8096 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8098 /* Determine the coupling coefficient for the distances
8099 * to the cell planes along dim0 and dim1 through dim2.
8100 * This is required for correct rounding.
8103 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8106 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8112 v_1 = ddbox->v[dim1];
8115 zone_cg_range = zones->cg_range;
8116 index_gl = dd->index_gl;
8117 cgindex = dd->cgindex;
8118 cginfo_mb = fr->cginfo_mb;
8120 zone_cg_range[0] = 0;
8121 zone_cg_range[1] = dd->ncg_home;
8122 comm->zone_ncg1[0] = dd->ncg_home;
8123 pos_cg = dd->ncg_home;
8125 nat_tot = dd->nat_home;
8127 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8129 dim = dd->dim[dim_ind];
8130 cd = &comm->cd[dim_ind];
8132 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8134 /* No pbc in this dimension, the first node should not comm. */
8142 v_d = ddbox->v[dim];
8143 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8145 cd->bInPlace = TRUE;
8146 for (p = 0; p < cd->np; p++)
8148 /* Only atoms communicated in the first pulse are used
8149 * for multi-body bonded interactions or for bBondComm.
8151 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8156 for (zone = 0; zone < nzone_send; zone++)
8158 if (tric_dist[dim_ind] && dim_ind > 0)
8160 /* Determine slightly more optimized skew_fac's
8162 * This reduces the number of communicated atoms
8163 * by about 10% for 3D DD of rhombic dodecahedra.
8165 for (dimd = 0; dimd < dim; dimd++)
8167 sf2_round[dimd] = 1;
8168 if (ddbox->tric_dir[dimd])
8170 for (i = dd->dim[dimd]+1; i < DIM; i++)
8172 /* If we are shifted in dimension i
8173 * and the cell plane is tilted forward
8174 * in dimension i, skip this coupling.
8176 if (!(zones->shift[nzone+zone][i] &&
8177 ddbox->v[dimd][i][dimd] >= 0))
8180 sqr(ddbox->v[dimd][i][dimd]);
8183 sf2_round[dimd] = 1/sf2_round[dimd];
8188 zonei = zone_perm[dim_ind][zone];
8191 /* Here we permutate the zones to obtain a convenient order
8192 * for neighbor searching
8194 cg0 = zone_cg_range[zonei];
8195 cg1 = zone_cg_range[zonei+1];
8199 /* Look only at the cg's received in the previous grid pulse
8201 cg1 = zone_cg_range[nzone+zone+1];
8202 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8205 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8206 for (th = 0; th < comm->nth; th++)
8208 gmx_domdec_ind_t *ind_p;
8209 int **ibuf_p, *ibuf_nalloc_p;
8211 int *nsend_p, *nat_p;
8217 /* Thread 0 writes in the comm buffers */
8219 ibuf_p = &comm->buf_int;
8220 ibuf_nalloc_p = &comm->nalloc_int;
8221 vbuf_p = &comm->vbuf;
8224 nsend_zone_p = &ind->nsend[zone];
8228 /* Other threads write into temp buffers */
8229 ind_p = &comm->dth[th].ind;
8230 ibuf_p = &comm->dth[th].ibuf;
8231 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8232 vbuf_p = &comm->dth[th].vbuf;
8233 nsend_p = &comm->dth[th].nsend;
8234 nat_p = &comm->dth[th].nat;
8235 nsend_zone_p = &comm->dth[th].nsend_zone;
8237 comm->dth[th].nsend = 0;
8238 comm->dth[th].nat = 0;
8239 comm->dth[th].nsend_zone = 0;
8249 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8250 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8253 /* Get the cg's for this pulse in this zone */
8254 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8256 dim, dim_ind, dim0, dim1, dim2,
8259 normal, skew_fac2_d, skew_fac_01,
8260 v_d, v_0, v_1, &corners, sf2_round,
8261 bDistBonded, bBondComm,
8265 ibuf_p, ibuf_nalloc_p,
8271 /* Append data of threads>=1 to the communication buffers */
8272 for (th = 1; th < comm->nth; th++)
8274 dd_comm_setup_work_t *dth;
8277 dth = &comm->dth[th];
8279 ns1 = nsend + dth->nsend_zone;
8280 if (ns1 > ind->nalloc)
8282 ind->nalloc = over_alloc_dd(ns1);
8283 srenew(ind->index, ind->nalloc);
8285 if (ns1 > comm->nalloc_int)
8287 comm->nalloc_int = over_alloc_dd(ns1);
8288 srenew(comm->buf_int, comm->nalloc_int);
8290 if (ns1 > comm->vbuf.nalloc)
8292 comm->vbuf.nalloc = over_alloc_dd(ns1);
8293 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8296 for (i = 0; i < dth->nsend_zone; i++)
8298 ind->index[nsend] = dth->ind.index[i];
8299 comm->buf_int[nsend] = dth->ibuf[i];
8300 copy_rvec(dth->vbuf.v[i],
8301 comm->vbuf.v[nsend]);
8305 ind->nsend[zone] += dth->nsend_zone;
8308 /* Clear the counts in case we do not have pbc */
8309 for (zone = nzone_send; zone < nzone; zone++)
8311 ind->nsend[zone] = 0;
8313 ind->nsend[nzone] = nsend;
8314 ind->nsend[nzone+1] = nat;
8315 /* Communicate the number of cg's and atoms to receive */
8316 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8317 ind->nsend, nzone+2,
8318 ind->nrecv, nzone+2);
8320 /* The rvec buffer is also required for atom buffers of size nsend
8321 * in dd_move_x and dd_move_f.
8323 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8327 /* We can receive in place if only the last zone is not empty */
8328 for (zone = 0; zone < nzone-1; zone++)
8330 if (ind->nrecv[zone] > 0)
8332 cd->bInPlace = FALSE;
8337 /* The int buffer is only required here for the cg indices */
8338 if (ind->nrecv[nzone] > comm->nalloc_int2)
8340 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8341 srenew(comm->buf_int2, comm->nalloc_int2);
8343 /* The rvec buffer is also required for atom buffers
8344 * of size nrecv in dd_move_x and dd_move_f.
8346 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8347 vec_rvec_check_alloc(&comm->vbuf2, i);
8351 /* Make space for the global cg indices */
8352 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8353 || dd->cg_nalloc == 0)
8355 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8356 srenew(index_gl, dd->cg_nalloc);
8357 srenew(cgindex, dd->cg_nalloc+1);
8359 /* Communicate the global cg indices */
8362 recv_i = index_gl + pos_cg;
8366 recv_i = comm->buf_int2;
8368 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8369 comm->buf_int, nsend,
8370 recv_i, ind->nrecv[nzone]);
8372 /* Make space for cg_cm */
8373 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8374 if (fr->cutoff_scheme == ecutsGROUP)
8382 /* Communicate cg_cm */
8385 recv_vr = cg_cm + pos_cg;
8389 recv_vr = comm->vbuf2.v;
8391 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8392 comm->vbuf.v, nsend,
8393 recv_vr, ind->nrecv[nzone]);
8395 /* Make the charge group index */
8398 zone = (p == 0 ? 0 : nzone - 1);
8399 while (zone < nzone)
8401 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8403 cg_gl = index_gl[pos_cg];
8404 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8405 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8406 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8409 /* Update the charge group presence,
8410 * so we can use it in the next pass of the loop.
8412 comm->bLocalCG[cg_gl] = TRUE;
8418 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8421 zone_cg_range[nzone+zone] = pos_cg;
8426 /* This part of the code is never executed with bBondComm. */
8427 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8428 index_gl, recv_i, cg_cm, recv_vr,
8429 cgindex, fr->cginfo_mb, fr->cginfo);
8430 pos_cg += ind->nrecv[nzone];
8432 nat_tot += ind->nrecv[nzone+1];
8436 /* Store the atom block for easy copying of communication buffers */
8437 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8441 dd->index_gl = index_gl;
8442 dd->cgindex = cgindex;
8444 dd->ncg_tot = zone_cg_range[zones->n];
8445 dd->nat_tot = nat_tot;
8446 comm->nat[ddnatHOME] = dd->nat_home;
8447 for (i = ddnatZONE; i < ddnatNR; i++)
8449 comm->nat[i] = dd->nat_tot;
8454 /* We don't need to update cginfo, since that was alrady done above.
8455 * So we pass NULL for the forcerec.
8457 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8458 NULL, comm->bLocalCG);
8463 fprintf(debug, "Finished setting up DD communication, zones:");
8464 for (c = 0; c < zones->n; c++)
8466 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8468 fprintf(debug, "\n");
8472 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8476 for (c = 0; c < zones->nizone; c++)
8478 zones->izone[c].cg1 = zones->cg_range[c+1];
8479 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8480 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8484 static void set_zones_size(gmx_domdec_t *dd,
8485 matrix box, const gmx_ddbox_t *ddbox,
8486 int zone_start, int zone_end)
8488 gmx_domdec_comm_t *comm;
8489 gmx_domdec_zones_t *zones;
8491 int z, zi, zj0, zj1, d, dim;
8494 real size_j, add_tric;
8499 zones = &comm->zones;
8501 /* Do we need to determine extra distances for multi-body bondeds? */
8502 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8504 for (z = zone_start; z < zone_end; z++)
8506 /* Copy cell limits to zone limits.
8507 * Valid for non-DD dims and non-shifted dims.
8509 copy_rvec(comm->cell_x0, zones->size[z].x0);
8510 copy_rvec(comm->cell_x1, zones->size[z].x1);
8513 for (d = 0; d < dd->ndim; d++)
8517 for (z = 0; z < zones->n; z++)
8519 /* With a staggered grid we have different sizes
8520 * for non-shifted dimensions.
8522 if (dd->bGridJump && zones->shift[z][dim] == 0)
8526 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8527 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8531 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8532 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8538 rcmbs = comm->cutoff_mbody;
8539 if (ddbox->tric_dir[dim])
8541 rcs /= ddbox->skew_fac[dim];
8542 rcmbs /= ddbox->skew_fac[dim];
8545 /* Set the lower limit for the shifted zone dimensions */
8546 for (z = zone_start; z < zone_end; z++)
8548 if (zones->shift[z][dim] > 0)
8551 if (!dd->bGridJump || d == 0)
8553 zones->size[z].x0[dim] = comm->cell_x1[dim];
8554 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8558 /* Here we take the lower limit of the zone from
8559 * the lowest domain of the zone below.
8563 zones->size[z].x0[dim] =
8564 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8570 zones->size[z].x0[dim] =
8571 zones->size[zone_perm[2][z-4]].x0[dim];
8575 zones->size[z].x0[dim] =
8576 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8579 /* A temporary limit, is updated below */
8580 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8584 for (zi = 0; zi < zones->nizone; zi++)
8586 if (zones->shift[zi][dim] == 0)
8588 /* This takes the whole zone into account.
8589 * With multiple pulses this will lead
8590 * to a larger zone then strictly necessary.
8592 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8593 zones->size[zi].x1[dim]+rcmbs);
8601 /* Loop over the i-zones to set the upper limit of each
8604 for (zi = 0; zi < zones->nizone; zi++)
8606 if (zones->shift[zi][dim] == 0)
8608 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8610 if (zones->shift[z][dim] > 0)
8612 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8613 zones->size[zi].x1[dim]+rcs);
8620 for (z = zone_start; z < zone_end; z++)
8622 /* Initialization only required to keep the compiler happy */
8623 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8626 /* To determine the bounding box for a zone we need to find
8627 * the extreme corners of 4, 2 or 1 corners.
8629 nc = 1 << (ddbox->npbcdim - 1);
8631 for (c = 0; c < nc; c++)
8633 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8637 corner[YY] = zones->size[z].x0[YY];
8641 corner[YY] = zones->size[z].x1[YY];
8645 corner[ZZ] = zones->size[z].x0[ZZ];
8649 corner[ZZ] = zones->size[z].x1[ZZ];
8651 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8653 /* With 1D domain decomposition the cg's are not in
8654 * the triclinic box, but triclinic x-y and rectangular y-z.
8655 * Shift y back, so it will later end up at 0.
8657 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8659 /* Apply the triclinic couplings */
8660 for (i = YY; i < ddbox->npbcdim; i++)
8662 for (j = XX; j < i; j++)
8664 corner[j] += corner[i]*box[i][j]/box[i][i];
8669 copy_rvec(corner, corner_min);
8670 copy_rvec(corner, corner_max);
8674 for (i = 0; i < DIM; i++)
8676 corner_min[i] = min(corner_min[i], corner[i]);
8677 corner_max[i] = max(corner_max[i], corner[i]);
8681 /* Copy the extreme cornes without offset along x */
8682 for (i = 0; i < DIM; i++)
8684 zones->size[z].bb_x0[i] = corner_min[i];
8685 zones->size[z].bb_x1[i] = corner_max[i];
8687 /* Add the offset along x */
8688 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8689 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8692 if (zone_start == 0)
8695 for (dim = 0; dim < DIM; dim++)
8697 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8699 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8704 for (z = zone_start; z < zone_end; z++)
8706 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8708 zones->size[z].x0[XX], zones->size[z].x1[XX],
8709 zones->size[z].x0[YY], zones->size[z].x1[YY],
8710 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8711 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8713 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8714 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8715 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8720 static int comp_cgsort(const void *a, const void *b)
8724 gmx_cgsort_t *cga, *cgb;
8725 cga = (gmx_cgsort_t *)a;
8726 cgb = (gmx_cgsort_t *)b;
8728 comp = cga->nsc - cgb->nsc;
8731 comp = cga->ind_gl - cgb->ind_gl;
8737 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8742 /* Order the data */
8743 for (i = 0; i < n; i++)
8745 buf[i] = a[sort[i].ind];
8748 /* Copy back to the original array */
8749 for (i = 0; i < n; i++)
8755 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8760 /* Order the data */
8761 for (i = 0; i < n; i++)
8763 copy_rvec(v[sort[i].ind], buf[i]);
8766 /* Copy back to the original array */
8767 for (i = 0; i < n; i++)
8769 copy_rvec(buf[i], v[i]);
8773 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8776 int a, atot, cg, cg0, cg1, i;
8778 if (cgindex == NULL)
8780 /* Avoid the useless loop of the atoms within a cg */
8781 order_vec_cg(ncg, sort, v, buf);
8786 /* Order the data */
8788 for (cg = 0; cg < ncg; cg++)
8790 cg0 = cgindex[sort[cg].ind];
8791 cg1 = cgindex[sort[cg].ind+1];
8792 for (i = cg0; i < cg1; i++)
8794 copy_rvec(v[i], buf[a]);
8800 /* Copy back to the original array */
8801 for (a = 0; a < atot; a++)
8803 copy_rvec(buf[a], v[a]);
8807 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8808 int nsort_new, gmx_cgsort_t *sort_new,
8809 gmx_cgsort_t *sort1)
8813 /* The new indices are not very ordered, so we qsort them */
8814 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8816 /* sort2 is already ordered, so now we can merge the two arrays */
8820 while (i2 < nsort2 || i_new < nsort_new)
8824 sort1[i1++] = sort_new[i_new++];
8826 else if (i_new == nsort_new)
8828 sort1[i1++] = sort2[i2++];
8830 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8831 (sort2[i2].nsc == sort_new[i_new].nsc &&
8832 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8834 sort1[i1++] = sort2[i2++];
8838 sort1[i1++] = sort_new[i_new++];
8843 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8845 gmx_domdec_sort_t *sort;
8846 gmx_cgsort_t *cgsort, *sort_i;
8847 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8848 int sort_last, sort_skip;
8850 sort = dd->comm->sort;
8852 a = fr->ns.grid->cell_index;
8854 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8856 if (ncg_home_old >= 0)
8858 /* The charge groups that remained in the same ns grid cell
8859 * are completely ordered. So we can sort efficiently by sorting
8860 * the charge groups that did move into the stationary list.
8865 for (i = 0; i < dd->ncg_home; i++)
8867 /* Check if this cg did not move to another node */
8870 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8872 /* This cg is new on this node or moved ns grid cell */
8873 if (nsort_new >= sort->sort_new_nalloc)
8875 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8876 srenew(sort->sort_new, sort->sort_new_nalloc);
8878 sort_i = &(sort->sort_new[nsort_new++]);
8882 /* This cg did not move */
8883 sort_i = &(sort->sort2[nsort2++]);
8885 /* Sort on the ns grid cell indices
8886 * and the global topology index.
8887 * index_gl is irrelevant with cell ns,
8888 * but we set it here anyhow to avoid a conditional.
8891 sort_i->ind_gl = dd->index_gl[i];
8898 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8901 /* Sort efficiently */
8902 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8907 cgsort = sort->sort;
8909 for (i = 0; i < dd->ncg_home; i++)
8911 /* Sort on the ns grid cell indices
8912 * and the global topology index
8914 cgsort[i].nsc = a[i];
8915 cgsort[i].ind_gl = dd->index_gl[i];
8917 if (cgsort[i].nsc < moved)
8924 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8926 /* Determine the order of the charge groups using qsort */
8927 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8933 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8936 int ncg_new, i, *a, na;
8938 sort = dd->comm->sort->sort;
8940 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8943 for (i = 0; i < na; i++)
8947 sort[ncg_new].ind = a[i];
8955 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8956 rvec *cgcm, t_forcerec *fr, t_state *state,
8959 gmx_domdec_sort_t *sort;
8960 gmx_cgsort_t *cgsort, *sort_i;
8962 int ncg_new, i, *ibuf, cgsize;
8965 sort = dd->comm->sort;
8967 if (dd->ncg_home > sort->sort_nalloc)
8969 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8970 srenew(sort->sort, sort->sort_nalloc);
8971 srenew(sort->sort2, sort->sort_nalloc);
8973 cgsort = sort->sort;
8975 switch (fr->cutoff_scheme)
8978 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8981 ncg_new = dd_sort_order_nbnxn(dd, fr);
8984 gmx_incons("unimplemented");
8988 /* We alloc with the old size, since cgindex is still old */
8989 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8990 vbuf = dd->comm->vbuf.v;
8994 cgindex = dd->cgindex;
9001 /* Remove the charge groups which are no longer at home here */
9002 dd->ncg_home = ncg_new;
9005 fprintf(debug, "Set the new home charge group count to %d\n",
9009 /* Reorder the state */
9010 for (i = 0; i < estNR; i++)
9012 if (EST_DISTR(i) && (state->flags & (1<<i)))
9017 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
9020 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
9023 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
9026 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
9030 case estDISRE_INITF:
9031 case estDISRE_RM3TAV:
9032 case estORIRE_INITF:
9034 /* No ordering required */
9037 gmx_incons("Unknown state entry encountered in dd_sort_state");
9042 if (fr->cutoff_scheme == ecutsGROUP)
9045 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9048 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9050 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9051 srenew(sort->ibuf, sort->ibuf_nalloc);
9054 /* Reorder the global cg index */
9055 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9056 /* Reorder the cginfo */
9057 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9058 /* Rebuild the local cg index */
9062 for (i = 0; i < dd->ncg_home; i++)
9064 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9065 ibuf[i+1] = ibuf[i] + cgsize;
9067 for (i = 0; i < dd->ncg_home+1; i++)
9069 dd->cgindex[i] = ibuf[i];
9074 for (i = 0; i < dd->ncg_home+1; i++)
9079 /* Set the home atom number */
9080 dd->nat_home = dd->cgindex[dd->ncg_home];
9082 if (fr->cutoff_scheme == ecutsVERLET)
9084 /* The atoms are now exactly in grid order, update the grid order */
9085 nbnxn_set_atomorder(fr->nbv->nbs);
9089 /* Copy the sorted ns cell indices back to the ns grid struct */
9090 for (i = 0; i < dd->ncg_home; i++)
9092 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9094 fr->ns.grid->nr = dd->ncg_home;
9098 static void add_dd_statistics(gmx_domdec_t *dd)
9100 gmx_domdec_comm_t *comm;
9105 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9107 comm->sum_nat[ddnat-ddnatZONE] +=
9108 comm->nat[ddnat] - comm->nat[ddnat-1];
9113 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9115 gmx_domdec_comm_t *comm;
9120 /* Reset all the statistics and counters for total run counting */
9121 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9123 comm->sum_nat[ddnat-ddnatZONE] = 0;
9127 comm->load_step = 0;
9130 clear_ivec(comm->load_lim);
9135 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9137 gmx_domdec_comm_t *comm;
9141 comm = cr->dd->comm;
9143 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9150 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");
9152 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9154 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9159 " av. #atoms communicated per step for force: %d x %.1f\n",
9163 if (cr->dd->vsite_comm)
9166 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9167 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9172 if (cr->dd->constraint_comm)
9175 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9176 1 + ir->nLincsIter, av);
9180 gmx_incons(" Unknown type for DD statistics");
9183 fprintf(fplog, "\n");
9185 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9187 print_dd_load_av(fplog, cr->dd);
9191 void dd_partition_system(FILE *fplog,
9192 gmx_large_int_t step,
9194 gmx_bool bMasterState,
9196 t_state *state_global,
9197 gmx_mtop_t *top_global,
9199 t_state *state_local,
9202 gmx_localtop_t *top_local,
9205 gmx_shellfc_t shellfc,
9206 gmx_constr_t constr,
9208 gmx_wallcycle_t wcycle,
9212 gmx_domdec_comm_t *comm;
9213 gmx_ddbox_t ddbox = {0};
9215 gmx_large_int_t step_pcoupl;
9216 rvec cell_ns_x0, cell_ns_x1;
9217 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9218 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9219 gmx_bool bRedist, bSortCG, bResortAll;
9220 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9227 bBoxChanged = (bMasterState || DEFORM(*ir));
9228 if (ir->epc != epcNO)
9230 /* With nstpcouple > 1 pressure coupling happens.
9231 * one step after calculating the pressure.
9232 * Box scaling happens at the end of the MD step,
9233 * after the DD partitioning.
9234 * We therefore have to do DLB in the first partitioning
9235 * after an MD step where P-coupling occured.
9236 * We need to determine the last step in which p-coupling occurred.
9237 * MRS -- need to validate this for vv?
9242 step_pcoupl = step - 1;
9246 step_pcoupl = ((step - 1)/n)*n + 1;
9248 if (step_pcoupl >= comm->partition_step)
9254 bNStGlobalComm = (step % nstglobalcomm == 0);
9256 if (!comm->bDynLoadBal)
9262 /* Should we do dynamic load balacing this step?
9263 * Since it requires (possibly expensive) global communication,
9264 * we might want to do DLB less frequently.
9266 if (bBoxChanged || ir->epc != epcNO)
9268 bDoDLB = bBoxChanged;
9272 bDoDLB = bNStGlobalComm;
9276 /* Check if we have recorded loads on the nodes */
9277 if (comm->bRecordLoad && dd_load_count(comm))
9279 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9281 /* Check if we should use DLB at the second partitioning
9282 * and every 100 partitionings,
9283 * so the extra communication cost is negligible.
9285 n = max(100, nstglobalcomm);
9286 bCheckDLB = (comm->n_load_collect == 0 ||
9287 comm->n_load_have % n == n-1);
9294 /* Print load every nstlog, first and last step to the log file */
9295 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9296 comm->n_load_collect == 0 ||
9298 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9300 /* Avoid extra communication due to verbose screen output
9301 * when nstglobalcomm is set.
9303 if (bDoDLB || bLogLoad || bCheckDLB ||
9304 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9306 get_load_distribution(dd, wcycle);
9311 dd_print_load(fplog, dd, step-1);
9315 dd_print_load_verbose(dd);
9318 comm->n_load_collect++;
9322 /* Since the timings are node dependent, the master decides */
9326 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9329 fprintf(debug, "step %s, imb loss %f\n",
9330 gmx_step_str(step, sbuf),
9331 dd_force_imb_perf_loss(dd));
9334 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9337 turn_on_dlb(fplog, cr, step);
9342 comm->n_load_have++;
9345 cgs_gl = &comm->cgs_gl;
9350 /* Clear the old state */
9351 clear_dd_indices(dd, 0, 0);
9354 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9355 TRUE, cgs_gl, state_global->x, &ddbox);
9357 get_cg_distribution(fplog, step, dd, cgs_gl,
9358 state_global->box, &ddbox, state_global->x);
9360 dd_distribute_state(dd, cgs_gl,
9361 state_global, state_local, f);
9363 dd_make_local_cgs(dd, &top_local->cgs);
9365 /* Ensure that we have space for the new distribution */
9366 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9368 if (fr->cutoff_scheme == ecutsGROUP)
9370 calc_cgcm(fplog, 0, dd->ncg_home,
9371 &top_local->cgs, state_local->x, fr->cg_cm);
9374 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9376 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9378 else if (state_local->ddp_count != dd->ddp_count)
9380 if (state_local->ddp_count > dd->ddp_count)
9382 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9385 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9387 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);
9390 /* Clear the old state */
9391 clear_dd_indices(dd, 0, 0);
9393 /* Build the new indices */
9394 rebuild_cgindex(dd, cgs_gl->index, state_local);
9395 make_dd_indices(dd, cgs_gl->index, 0);
9396 ncgindex_set = dd->ncg_home;
9398 if (fr->cutoff_scheme == ecutsGROUP)
9400 /* Redetermine the cg COMs */
9401 calc_cgcm(fplog, 0, dd->ncg_home,
9402 &top_local->cgs, state_local->x, fr->cg_cm);
9405 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9407 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9409 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9410 TRUE, &top_local->cgs, state_local->x, &ddbox);
9412 bRedist = comm->bDynLoadBal;
9416 /* We have the full state, only redistribute the cgs */
9418 /* Clear the non-home indices */
9419 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9422 /* Avoid global communication for dim's without pbc and -gcom */
9423 if (!bNStGlobalComm)
9425 copy_rvec(comm->box0, ddbox.box0 );
9426 copy_rvec(comm->box_size, ddbox.box_size);
9428 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9429 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9434 /* For dim's without pbc and -gcom */
9435 copy_rvec(ddbox.box0, comm->box0 );
9436 copy_rvec(ddbox.box_size, comm->box_size);
9438 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9441 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9443 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9446 /* Check if we should sort the charge groups */
9447 if (comm->nstSortCG > 0)
9449 bSortCG = (bMasterState ||
9450 (bRedist && (step % comm->nstSortCG == 0)));
9457 ncg_home_old = dd->ncg_home;
9462 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9464 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9465 state_local, f, fr, mdatoms,
9466 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9468 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9471 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9473 &comm->cell_x0, &comm->cell_x1,
9474 dd->ncg_home, fr->cg_cm,
9475 cell_ns_x0, cell_ns_x1, &grid_density);
9479 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9482 switch (fr->cutoff_scheme)
9485 copy_ivec(fr->ns.grid->n, ncells_old);
9486 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9487 state_local->box, cell_ns_x0, cell_ns_x1,
9488 fr->rlistlong, grid_density);
9491 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9494 gmx_incons("unimplemented");
9496 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9497 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9501 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9503 /* Sort the state on charge group position.
9504 * This enables exact restarts from this step.
9505 * It also improves performance by about 15% with larger numbers
9506 * of atoms per node.
9509 /* Fill the ns grid with the home cell,
9510 * so we can sort with the indices.
9512 set_zones_ncg_home(dd);
9514 switch (fr->cutoff_scheme)
9517 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9519 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9521 comm->zones.size[0].bb_x0,
9522 comm->zones.size[0].bb_x1,
9524 comm->zones.dens_zone0,
9527 ncg_moved, bRedist ? comm->moved : NULL,
9528 fr->nbv->grp[eintLocal].kernel_type,
9529 fr->nbv->grp[eintLocal].nbat);
9531 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9534 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9535 0, dd->ncg_home, fr->cg_cm);
9537 copy_ivec(fr->ns.grid->n, ncells_new);
9540 gmx_incons("unimplemented");
9543 bResortAll = bMasterState;
9545 /* Check if we can user the old order and ns grid cell indices
9546 * of the charge groups to sort the charge groups efficiently.
9548 if (ncells_new[XX] != ncells_old[XX] ||
9549 ncells_new[YY] != ncells_old[YY] ||
9550 ncells_new[ZZ] != ncells_old[ZZ])
9557 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9558 gmx_step_str(step, sbuf), dd->ncg_home);
9560 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9561 bResortAll ? -1 : ncg_home_old);
9562 /* Rebuild all the indices */
9563 ga2la_clear(dd->ga2la);
9566 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9569 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9571 /* Setup up the communication and communicate the coordinates */
9572 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9574 /* Set the indices */
9575 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9577 /* Set the charge group boundaries for neighbor searching */
9578 set_cg_boundaries(&comm->zones);
9580 if (fr->cutoff_scheme == ecutsVERLET)
9582 set_zones_size(dd, state_local->box, &ddbox,
9583 bSortCG ? 1 : 0, comm->zones.n);
9586 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9589 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9590 -1,state_local->x,state_local->box);
9593 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9595 /* Extract a local topology from the global topology */
9596 for (i = 0; i < dd->ndim; i++)
9598 np[dd->dim[i]] = comm->cd[i].np;
9600 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9601 comm->cellsize_min, np,
9603 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9604 vsite, top_global, top_local);
9606 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9608 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9610 /* Set up the special atom communication */
9611 n = comm->nat[ddnatZONE];
9612 for (i = ddnatZONE+1; i < ddnatNR; i++)
9617 if (vsite && vsite->n_intercg_vsite)
9619 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9623 if (dd->bInterCGcons || dd->bInterCGsettles)
9625 /* Only for inter-cg constraints we need special code */
9626 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9627 constr, ir->nProjOrder,
9628 top_local->idef.il);
9632 gmx_incons("Unknown special atom type setup");
9637 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9639 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9641 /* Make space for the extra coordinates for virtual site
9642 * or constraint communication.
9644 state_local->natoms = comm->nat[ddnatNR-1];
9645 if (state_local->natoms > state_local->nalloc)
9647 dd_realloc_state(state_local, f, state_local->natoms);
9650 if (fr->bF_NoVirSum)
9652 if (vsite && vsite->n_intercg_vsite)
9654 nat_f_novirsum = comm->nat[ddnatVSITE];
9658 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9660 nat_f_novirsum = dd->nat_tot;
9664 nat_f_novirsum = dd->nat_home;
9673 /* Set the number of atoms required for the force calculation.
9674 * Forces need to be constrained when using a twin-range setup
9675 * or with energy minimization. For simple simulations we could
9676 * avoid some allocation, zeroing and copying, but this is
9677 * probably not worth the complications ande checking.
9679 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9680 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9682 /* We make the all mdatoms up to nat_tot_con.
9683 * We could save some work by only setting invmass
9684 * between nat_tot and nat_tot_con.
9686 /* This call also sets the new number of home particles to dd->nat_home */
9687 atoms2md(top_global, ir,
9688 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9690 /* Now we have the charges we can sort the FE interactions */
9691 dd_sort_local_top(dd, mdatoms, top_local);
9695 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9696 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9701 /* Make the local shell stuff, currently no communication is done */
9702 make_local_shells(cr, mdatoms, shellfc);
9705 if (ir->implicit_solvent)
9707 make_local_gb(cr, fr->born, ir->gb_algorithm);
9710 setup_bonded_threading(fr, &top_local->idef);
9712 if (!(cr->duty & DUTY_PME))
9714 /* Send the charges to our PME only node */
9715 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9716 mdatoms->chargeA, mdatoms->chargeB,
9717 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9722 set_constraints(constr, top_local, ir, mdatoms, cr);
9725 if (ir->ePull != epullNO)
9727 /* Update the local pull groups */
9728 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9733 /* Update the local rotation groups */
9734 dd_make_local_rotation_groups(dd, ir->rot);
9738 add_dd_statistics(dd);
9740 /* Make sure we only count the cycles for this DD partitioning */
9741 clear_dd_cycle_counts(dd);
9743 /* Because the order of the atoms might have changed since
9744 * the last vsite construction, we need to communicate the constructing
9745 * atom coordinates again (for spreading the forces this MD step).
9747 dd_move_x_vsites(dd, state_local->box, state_local->x);
9749 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9751 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9753 dd_move_x(dd, state_local->box, state_local->x);
9754 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9755 -1, state_local->x, state_local->box);
9758 /* Store the partitioning step */
9759 comm->partition_step = step;
9761 /* Increase the DD partitioning counter */
9763 /* The state currently matches this DD partitioning count, store it */
9764 state_local->ddp_count = dd->ddp_count;
9767 /* The DD master node knows the complete cg distribution,
9768 * store the count so we can possibly skip the cg info communication.
9770 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9773 if (comm->DD_debug > 0)
9775 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9776 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9777 "after partitioning");