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);
1555 for (i = 0; i < state_local->ngtc; i++)
1557 for (j = 0; j < nh; j++)
1559 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1560 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1562 state->therm_integral[i] = state_local->therm_integral[i];
1564 for (i = 0; i < state_local->nnhpres; i++)
1566 for (j = 0; j < nh; j++)
1568 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1569 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1573 for (est = 0; est < estNR; est++)
1575 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1580 dd_collect_vec(dd, state_local, state_local->x, state->x);
1583 dd_collect_vec(dd, state_local, state_local->v, state->v);
1586 dd_collect_vec(dd, state_local, state_local->sd_X, state->sd_X);
1589 dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
1592 if (state->nrngi == 1)
1596 for (i = 0; i < state_local->nrng; i++)
1598 state->ld_rng[i] = state_local->ld_rng[i];
1604 dd_gather(dd, state_local->nrng*sizeof(state->ld_rng[0]),
1605 state_local->ld_rng, state->ld_rng);
1609 if (state->nrngi == 1)
1613 state->ld_rngi[0] = state_local->ld_rngi[0];
1618 dd_gather(dd, sizeof(state->ld_rngi[0]),
1619 state_local->ld_rngi, state->ld_rngi);
1622 case estDISRE_INITF:
1623 case estDISRE_RM3TAV:
1624 case estORIRE_INITF:
1628 gmx_incons("Unknown state entry encountered in dd_collect_state");
1634 static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
1640 fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
1643 state->nalloc = over_alloc_dd(nalloc);
1645 for (est = 0; est < estNR; est++)
1647 if (EST_DISTR(est) && (state->flags & (1<<est)))
1652 srenew(state->x, state->nalloc);
1655 srenew(state->v, state->nalloc);
1658 srenew(state->sd_X, state->nalloc);
1661 srenew(state->cg_p, state->nalloc);
1665 case estDISRE_INITF:
1666 case estDISRE_RM3TAV:
1667 case estORIRE_INITF:
1669 /* No reallocation required */
1672 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1679 srenew(*f, state->nalloc);
1683 static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
1686 if (nalloc > fr->cg_nalloc)
1690 fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
1692 fr->cg_nalloc = over_alloc_dd(nalloc);
1693 srenew(fr->cginfo, fr->cg_nalloc);
1694 if (fr->cutoff_scheme == ecutsGROUP)
1696 srenew(fr->cg_cm, fr->cg_nalloc);
1699 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1701 /* We don't use charge groups, we use x in state to set up
1702 * the atom communication.
1704 dd_realloc_state(state, f, nalloc);
1708 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, t_block *cgs,
1711 gmx_domdec_master_t *ma;
1712 int n, i, c, a, nalloc = 0;
1719 for (n = 0; n < dd->nnodes; n++)
1723 if (ma->nat[n] > nalloc)
1725 nalloc = over_alloc_dd(ma->nat[n]);
1726 srenew(buf, nalloc);
1728 /* Use lv as a temporary buffer */
1730 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1732 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1734 copy_rvec(v[c], buf[a++]);
1737 if (a != ma->nat[n])
1739 gmx_fatal(FARGS, "Internal error a (%d) != nat (%d)",
1744 MPI_Send(buf, ma->nat[n]*sizeof(rvec), MPI_BYTE,
1745 DDRANK(dd, n), n, dd->mpi_comm_all);
1750 n = DDMASTERRANK(dd);
1752 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1754 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1756 copy_rvec(v[c], lv[a++]);
1763 MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
1764 MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
1769 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
1772 gmx_domdec_master_t *ma;
1773 int *scounts = NULL, *disps = NULL;
1774 int n, i, c, a, nalloc = 0;
1781 get_commbuffer_counts(dd, &scounts, &disps);
1785 for (n = 0; n < dd->nnodes; n++)
1787 for (i = ma->index[n]; i < ma->index[n+1]; i++)
1789 for (c = cgs->index[ma->cg[i]]; c < cgs->index[ma->cg[i]+1]; c++)
1791 copy_rvec(v[c], buf[a++]);
1797 dd_scatterv(dd, scounts, disps, buf, dd->nat_home*sizeof(rvec), lv);
1800 static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
1802 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1804 dd_distribute_vec_sendrecv(dd, cgs, v, lv);
1808 dd_distribute_vec_scatterv(dd, cgs, v, lv);
1812 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
1813 t_state *state, t_state *state_local,
1818 nh = state->nhchainlength;
1822 for (i = 0; i < efptNR; i++)
1824 state_local->lambda[i] = state->lambda[i];
1826 state_local->fep_state = state->fep_state;
1827 state_local->veta = state->veta;
1828 state_local->vol0 = state->vol0;
1829 copy_mat(state->box, state_local->box);
1830 copy_mat(state->box_rel, state_local->box_rel);
1831 copy_mat(state->boxv, state_local->boxv);
1832 copy_mat(state->svir_prev, state_local->svir_prev);
1833 copy_mat(state->fvir_prev, state_local->fvir_prev);
1834 for (i = 0; i < state_local->ngtc; i++)
1836 for (j = 0; j < nh; j++)
1838 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1839 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1841 state_local->therm_integral[i] = state->therm_integral[i];
1843 for (i = 0; i < state_local->nnhpres; i++)
1845 for (j = 0; j < nh; j++)
1847 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1848 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1852 dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
1853 dd_bcast(dd, sizeof(int), &state_local->fep_state);
1854 dd_bcast(dd, sizeof(real), &state_local->veta);
1855 dd_bcast(dd, sizeof(real), &state_local->vol0);
1856 dd_bcast(dd, sizeof(state_local->box), state_local->box);
1857 dd_bcast(dd, sizeof(state_local->box_rel), state_local->box_rel);
1858 dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
1859 dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
1860 dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
1861 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
1862 dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
1863 dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
1864 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
1865 dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
1867 if (dd->nat_home > state_local->nalloc)
1869 dd_realloc_state(state_local, f, dd->nat_home);
1871 for (i = 0; i < estNR; i++)
1873 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1878 dd_distribute_vec(dd, cgs, state->x, state_local->x);
1881 dd_distribute_vec(dd, cgs, state->v, state_local->v);
1884 dd_distribute_vec(dd, cgs, state->sd_X, state_local->sd_X);
1887 dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
1890 if (state->nrngi == 1)
1893 state_local->nrng*sizeof(state_local->ld_rng[0]),
1894 state->ld_rng, state_local->ld_rng);
1899 state_local->nrng*sizeof(state_local->ld_rng[0]),
1900 state->ld_rng, state_local->ld_rng);
1904 if (state->nrngi == 1)
1906 dd_bcastc(dd, sizeof(state_local->ld_rngi[0]),
1907 state->ld_rngi, state_local->ld_rngi);
1911 dd_scatter(dd, sizeof(state_local->ld_rngi[0]),
1912 state->ld_rngi, state_local->ld_rngi);
1915 case estDISRE_INITF:
1916 case estDISRE_RM3TAV:
1917 case estORIRE_INITF:
1919 /* Not implemented yet */
1922 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1928 static char dim2char(int dim)
1934 case XX: c = 'X'; break;
1935 case YY: c = 'Y'; break;
1936 case ZZ: c = 'Z'; break;
1937 default: gmx_fatal(FARGS, "Unknown dim %d", dim);
1943 static void write_dd_grid_pdb(const char *fn, gmx_large_int_t step,
1944 gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1946 rvec grid_s[2], *grid_r = NULL, cx, r;
1947 char fname[STRLEN], format[STRLEN], buf[22];
1949 int a, i, d, z, y, x;
1953 copy_rvec(dd->comm->cell_x0, grid_s[0]);
1954 copy_rvec(dd->comm->cell_x1, grid_s[1]);
1958 snew(grid_r, 2*dd->nnodes);
1961 dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
1965 for (d = 0; d < DIM; d++)
1967 for (i = 0; i < DIM; i++)
1975 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1977 tric[d][i] = box[i][d]/box[i][i];
1986 sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1987 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
1988 out = gmx_fio_fopen(fname, "w");
1989 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1991 for (i = 0; i < dd->nnodes; i++)
1993 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1994 for (d = 0; d < DIM; d++)
1996 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1998 for (z = 0; z < 2; z++)
2000 for (y = 0; y < 2; y++)
2002 for (x = 0; x < 2; x++)
2004 cx[XX] = grid_r[i*2+x][XX];
2005 cx[YY] = grid_r[i*2+y][YY];
2006 cx[ZZ] = grid_r[i*2+z][ZZ];
2008 fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
2009 ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
2013 for (d = 0; d < DIM; d++)
2015 for (x = 0; x < 4; x++)
2019 case 0: y = 1 + i*8 + 2*x; break;
2020 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
2021 case 2: y = 1 + i*8 + x; break;
2023 fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
2027 gmx_fio_fclose(out);
2032 void write_dd_pdb(const char *fn, gmx_large_int_t step, const char *title,
2033 gmx_mtop_t *mtop, t_commrec *cr,
2034 int natoms, rvec x[], matrix box)
2036 char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
2038 int i, ii, resnr, c;
2039 char *atomname, *resname;
2046 natoms = dd->comm->nat[ddnatVSITE];
2049 sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
2051 sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
2052 sprintf(format4, "%s%s\n", pdbformat4, "%6.2f%6.2f");
2054 out = gmx_fio_fopen(fname, "w");
2056 fprintf(out, "TITLE %s\n", title);
2057 gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
2058 for (i = 0; i < natoms; i++)
2060 ii = dd->gatindex[i];
2061 gmx_mtop_atominfo_global(mtop, ii, &atomname, &resnr, &resname);
2062 if (i < dd->comm->nat[ddnatZONE])
2065 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2071 else if (i < dd->comm->nat[ddnatVSITE])
2073 b = dd->comm->zones.n;
2077 b = dd->comm->zones.n + 1;
2079 fprintf(out, strlen(atomname) < 4 ? format : format4,
2080 "ATOM", (ii+1)%100000,
2081 atomname, resname, ' ', resnr%10000, ' ',
2082 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
2084 fprintf(out, "TER\n");
2086 gmx_fio_fclose(out);
2089 real dd_cutoff_mbody(gmx_domdec_t *dd)
2091 gmx_domdec_comm_t *comm;
2098 if (comm->bInterCGBondeds)
2100 if (comm->cutoff_mbody > 0)
2102 r = comm->cutoff_mbody;
2106 /* cutoff_mbody=0 means we do not have DLB */
2107 r = comm->cellsize_min[dd->dim[0]];
2108 for (di = 1; di < dd->ndim; di++)
2110 r = min(r, comm->cellsize_min[dd->dim[di]]);
2112 if (comm->bBondComm)
2114 r = max(r, comm->cutoff_mbody);
2118 r = min(r, comm->cutoff);
2126 real dd_cutoff_twobody(gmx_domdec_t *dd)
2130 r_mb = dd_cutoff_mbody(dd);
2132 return max(dd->comm->cutoff, r_mb);
2136 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd, ivec coord, ivec coord_pme)
2140 nc = dd->nc[dd->comm->cartpmedim];
2141 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2142 copy_ivec(coord, coord_pme);
2143 coord_pme[dd->comm->cartpmedim] =
2144 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2147 static int low_ddindex2pmeindex(int ndd, int npme, int ddindex)
2149 /* Here we assign a PME node to communicate with this DD node
2150 * by assuming that the major index of both is x.
2151 * We add cr->npmenodes/2 to obtain an even distribution.
2153 return (ddindex*npme + npme/2)/ndd;
2156 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
2158 return low_ddindex2pmeindex(dd->nnodes, dd->comm->npmenodes, ddindex);
2161 static int cr_ddindex2pmeindex(const t_commrec *cr, int ddindex)
2163 return low_ddindex2pmeindex(cr->dd->nnodes, cr->npmenodes, ddindex);
2166 static int *dd_pmenodes(t_commrec *cr)
2171 snew(pmenodes, cr->npmenodes);
2173 for (i = 0; i < cr->dd->nnodes; i++)
2175 p0 = cr_ddindex2pmeindex(cr, i);
2176 p1 = cr_ddindex2pmeindex(cr, i+1);
2177 if (i+1 == cr->dd->nnodes || p1 > p0)
2181 fprintf(debug, "pmenode[%d] = %d\n", n, i+1+n);
2183 pmenodes[n] = i + 1 + n;
2191 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
2194 ivec coords, coords_pme, nc;
2199 if (dd->comm->bCartesian) {
2200 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2201 dd_coords2pmecoords(dd,coords,coords_pme);
2202 copy_ivec(dd->ntot,nc);
2203 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2204 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2206 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2208 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2214 slab = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
2219 static int ddcoord2simnodeid(t_commrec *cr, int x, int y, int z)
2221 gmx_domdec_comm_t *comm;
2223 int ddindex, nodeid = -1;
2225 comm = cr->dd->comm;
2230 if (comm->bCartesianPP_PME)
2233 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
2238 ddindex = dd_index(cr->dd->nc, coords);
2239 if (comm->bCartesianPP)
2241 nodeid = comm->ddindex2simnodeid[ddindex];
2247 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
2259 static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
2262 gmx_domdec_comm_t *comm;
2263 ivec coord, coord_pme;
2270 /* This assumes a uniform x domain decomposition grid cell size */
2271 if (comm->bCartesianPP_PME)
2274 MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
2275 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2277 /* This is a PP node */
2278 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2279 MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
2283 else if (comm->bCartesianPP)
2285 if (sim_nodeid < dd->nnodes)
2287 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2292 /* This assumes DD cells with identical x coordinates
2293 * are numbered sequentially.
2295 if (dd->comm->pmenodes == NULL)
2297 if (sim_nodeid < dd->nnodes)
2299 /* The DD index equals the nodeid */
2300 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
2306 while (sim_nodeid > dd->comm->pmenodes[i])
2310 if (sim_nodeid < dd->comm->pmenodes[i])
2312 pmenode = dd->comm->pmenodes[i];
2320 void get_pme_nnodes(const gmx_domdec_t *dd,
2321 int *npmenodes_x, int *npmenodes_y)
2325 *npmenodes_x = dd->comm->npmenodes_x;
2326 *npmenodes_y = dd->comm->npmenodes_y;
2335 gmx_bool gmx_pmeonlynode(t_commrec *cr, int sim_nodeid)
2337 gmx_bool bPMEOnlyNode;
2339 if (DOMAINDECOMP(cr))
2341 bPMEOnlyNode = (dd_simnode2pmenode(cr, sim_nodeid) == -1);
2345 bPMEOnlyNode = FALSE;
2348 return bPMEOnlyNode;
2351 void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
2352 int *nmy_ddnodes, int **my_ddnodes, int *node_peer)
2356 ivec coord, coord_pme;
2360 snew(*my_ddnodes, (dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2363 for (x = 0; x < dd->nc[XX]; x++)
2365 for (y = 0; y < dd->nc[YY]; y++)
2367 for (z = 0; z < dd->nc[ZZ]; z++)
2369 if (dd->comm->bCartesianPP_PME)
2374 dd_cart_coord2pmecoord(dd, coord, coord_pme);
2375 if (dd->ci[XX] == coord_pme[XX] &&
2376 dd->ci[YY] == coord_pme[YY] &&
2377 dd->ci[ZZ] == coord_pme[ZZ])
2379 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2384 /* The slab corresponds to the nodeid in the PME group */
2385 if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
2387 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr, x, y, z);
2394 /* The last PP-only node is the peer node */
2395 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2399 fprintf(debug, "Receive coordinates from PP nodes:");
2400 for (x = 0; x < *nmy_ddnodes; x++)
2402 fprintf(debug, " %d", (*my_ddnodes)[x]);
2404 fprintf(debug, "\n");
2408 static gmx_bool receive_vir_ener(t_commrec *cr)
2410 gmx_domdec_comm_t *comm;
2411 int pmenode, coords[DIM], rank;
2415 if (cr->npmenodes < cr->dd->nnodes)
2417 comm = cr->dd->comm;
2418 if (comm->bCartesianPP_PME)
2420 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2422 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
2423 coords[comm->cartpmedim]++;
2424 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2426 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
2427 if (dd_simnode2pmenode(cr, rank) == pmenode)
2429 /* This is not the last PP node for pmenode */
2437 pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
2438 if (cr->sim_nodeid+1 < cr->nnodes &&
2439 dd_simnode2pmenode(cr, cr->sim_nodeid+1) == pmenode)
2441 /* This is not the last PP node for pmenode */
2450 static void set_zones_ncg_home(gmx_domdec_t *dd)
2452 gmx_domdec_zones_t *zones;
2455 zones = &dd->comm->zones;
2457 zones->cg_range[0] = 0;
2458 for (i = 1; i < zones->n+1; i++)
2460 zones->cg_range[i] = dd->ncg_home;
2462 /* zone_ncg1[0] should always be equal to ncg_home */
2463 dd->comm->zone_ncg1[0] = dd->ncg_home;
2466 static void rebuild_cgindex(gmx_domdec_t *dd,
2467 const int *gcgs_index, t_state *state)
2469 int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
2472 dd_cg_gl = dd->index_gl;
2473 cgindex = dd->cgindex;
2476 for (i = 0; i < state->ncg_gl; i++)
2480 dd_cg_gl[i] = cg_gl;
2481 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2485 dd->ncg_home = state->ncg_gl;
2488 set_zones_ncg_home(dd);
2491 static int ddcginfo(const cginfo_mb_t *cginfo_mb, int cg)
2493 while (cg >= cginfo_mb->cg_end)
2498 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2501 static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
2502 t_forcerec *fr, char *bLocalCG)
2504 cginfo_mb_t *cginfo_mb;
2510 cginfo_mb = fr->cginfo_mb;
2511 cginfo = fr->cginfo;
2513 for (cg = cg0; cg < cg1; cg++)
2515 cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
2519 if (bLocalCG != NULL)
2521 for (cg = cg0; cg < cg1; cg++)
2523 bLocalCG[index_gl[cg]] = TRUE;
2528 static void make_dd_indices(gmx_domdec_t *dd,
2529 const int *gcgs_index, int cg_start)
2531 int nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
2532 int *zone2cg, *zone_ncg1, *index_gl, *gatindex;
2537 bLocalCG = dd->comm->bLocalCG;
2539 if (dd->nat_tot > dd->gatindex_nalloc)
2541 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2542 srenew(dd->gatindex, dd->gatindex_nalloc);
2545 nzone = dd->comm->zones.n;
2546 zone2cg = dd->comm->zones.cg_range;
2547 zone_ncg1 = dd->comm->zone_ncg1;
2548 index_gl = dd->index_gl;
2549 gatindex = dd->gatindex;
2550 bCGs = dd->comm->bCGs;
2552 if (zone2cg[1] != dd->ncg_home)
2554 gmx_incons("dd->ncg_zone is not up to date");
2557 /* Make the local to global and global to local atom index */
2558 a = dd->cgindex[cg_start];
2559 for (zone = 0; zone < nzone; zone++)
2567 cg0 = zone2cg[zone];
2569 cg1 = zone2cg[zone+1];
2570 cg1_p1 = cg0 + zone_ncg1[zone];
2572 for (cg = cg0; cg < cg1; cg++)
2577 /* Signal that this cg is from more than one pulse away */
2580 cg_gl = index_gl[cg];
2583 for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
2586 ga2la_set(dd->ga2la, a_gl, a, zone1);
2592 gatindex[a] = cg_gl;
2593 ga2la_set(dd->ga2la, cg_gl, a, zone1);
2600 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
2603 int ncg, i, ngl, nerr;
2606 if (bLocalCG == NULL)
2610 for (i = 0; i < dd->ncg_tot; i++)
2612 if (!bLocalCG[dd->index_gl[i]])
2615 "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);
2620 for (i = 0; i < ncg_sys; i++)
2627 if (ngl != dd->ncg_tot)
2629 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);
2636 static void check_index_consistency(gmx_domdec_t *dd,
2637 int natoms_sys, int ncg_sys,
2640 int nerr, ngl, i, a, cell;
2645 if (dd->comm->DD_debug > 1)
2647 snew(have, natoms_sys);
2648 for (a = 0; a < dd->nat_tot; a++)
2650 if (have[dd->gatindex[a]] > 0)
2652 fprintf(stderr, "DD node %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
2656 have[dd->gatindex[a]] = a + 1;
2662 snew(have, dd->nat_tot);
2665 for (i = 0; i < natoms_sys; i++)
2667 if (ga2la_get(dd->ga2la, i, &a, &cell))
2669 if (a >= dd->nat_tot)
2671 fprintf(stderr, "DD 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);
2677 if (dd->gatindex[a] != i)
2679 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);
2686 if (ngl != dd->nat_tot)
2689 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2690 dd->rank, where, ngl, dd->nat_tot);
2692 for (a = 0; a < dd->nat_tot; a++)
2697 "DD node %d, %s: local atom %d, global %d has no global index\n",
2698 dd->rank, where, a+1, dd->gatindex[a]+1);
2703 nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
2707 gmx_fatal(FARGS, "DD node %d, %s: %d atom/cg index inconsistencies",
2708 dd->rank, where, nerr);
2712 static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
2719 /* Clear the whole list without searching */
2720 ga2la_clear(dd->ga2la);
2724 for (i = a_start; i < dd->nat_tot; i++)
2726 ga2la_del(dd->ga2la, dd->gatindex[i]);
2730 bLocalCG = dd->comm->bLocalCG;
2733 for (i = cg_start; i < dd->ncg_tot; i++)
2735 bLocalCG[dd->index_gl[i]] = FALSE;
2739 dd_clear_local_vsite_indices(dd);
2741 if (dd->constraints)
2743 dd_clear_local_constraint_indices(dd);
2747 /* This function should be used for moving the domain boudaries during DLB,
2748 * for obtaining the minimum cell size. It checks the initially set limit
2749 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2750 * and, possibly, a longer cut-off limit set for PME load balancing.
2752 static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
2756 cellsize_min = comm->cellsize_min[dim];
2758 if (!comm->bVacDLBNoLimit)
2760 /* The cut-off might have changed, e.g. by PME load balacning,
2761 * from the value used to set comm->cellsize_min, so check it.
2763 cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
2765 if (comm->bPMELoadBalDLBLimits)
2767 /* Check for the cut-off limit set by the PME load balancing */
2768 cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2772 return cellsize_min;
2775 static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
2778 real grid_jump_limit;
2780 /* The distance between the boundaries of cells at distance
2781 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2782 * and by the fact that cells should not be shifted by more than
2783 * half their size, such that cg's only shift by one cell
2784 * at redecomposition.
2786 grid_jump_limit = comm->cellsize_limit;
2787 if (!comm->bVacDLBNoLimit)
2789 if (comm->bPMELoadBalDLBLimits)
2791 cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
2793 grid_jump_limit = max(grid_jump_limit,
2794 cutoff/comm->cd[dim_ind].np);
2797 return grid_jump_limit;
2800 static gmx_bool check_grid_jump(gmx_large_int_t step,
2806 gmx_domdec_comm_t *comm;
2815 for (d = 1; d < dd->ndim; d++)
2818 limit = grid_jump_limit(comm, cutoff, d);
2819 bfac = ddbox->box_size[dim];
2820 if (ddbox->tric_dir[dim])
2822 bfac *= ddbox->skew_fac[dim];
2824 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2825 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2833 /* This error should never be triggered under normal
2834 * circumstances, but you never know ...
2836 gmx_fatal(FARGS, "Step %s: The domain decomposition grid has shifted too much in the %c-direction around cell %d %d %d. This should not have happened. Running with less nodes might avoid this issue.",
2837 gmx_step_str(step, buf),
2838 dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2846 static int dd_load_count(gmx_domdec_comm_t *comm)
2848 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2851 static float dd_force_load(gmx_domdec_comm_t *comm)
2858 if (comm->eFlop > 1)
2860 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2865 load = comm->cycl[ddCyclF];
2866 if (comm->cycl_n[ddCyclF] > 1)
2868 /* Subtract the maximum of the last n cycle counts
2869 * to get rid of possible high counts due to other soures,
2870 * for instance system activity, that would otherwise
2871 * affect the dynamic load balancing.
2873 load -= comm->cycl_max[ddCyclF];
2880 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
2882 gmx_domdec_comm_t *comm;
2887 snew(*dim_f, dd->nc[dim]+1);
2889 for (i = 1; i < dd->nc[dim]; i++)
2891 if (comm->slb_frac[dim])
2893 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2897 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2900 (*dim_f)[dd->nc[dim]] = 1;
2903 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
2905 int pmeindex, slab, nso, i;
2908 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2914 ddpme->dim = dimind;
2916 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2918 ddpme->nslab = (ddpme->dim == 0 ?
2919 dd->comm->npmenodes_x :
2920 dd->comm->npmenodes_y);
2922 if (ddpme->nslab <= 1)
2927 nso = dd->comm->npmenodes/ddpme->nslab;
2928 /* Determine for each PME slab the PP location range for dimension dim */
2929 snew(ddpme->pp_min, ddpme->nslab);
2930 snew(ddpme->pp_max, ddpme->nslab);
2931 for (slab = 0; slab < ddpme->nslab; slab++)
2933 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2934 ddpme->pp_max[slab] = 0;
2936 for (i = 0; i < dd->nnodes; i++)
2938 ddindex2xyz(dd->nc, i, xyz);
2939 /* For y only use our y/z slab.
2940 * This assumes that the PME x grid size matches the DD grid size.
2942 if (dimind == 0 || xyz[XX] == dd->ci[XX])
2944 pmeindex = ddindex2pmeindex(dd, i);
2947 slab = pmeindex/nso;
2951 slab = pmeindex % ddpme->nslab;
2953 ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
2954 ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
2958 set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
2961 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2963 if (dd->comm->ddpme[0].dim == XX)
2965 return dd->comm->ddpme[0].maxshift;
2973 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2975 if (dd->comm->ddpme[0].dim == YY)
2977 return dd->comm->ddpme[0].maxshift;
2979 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2981 return dd->comm->ddpme[1].maxshift;
2989 static void set_pme_maxshift(gmx_domdec_t *dd, gmx_ddpme_t *ddpme,
2990 gmx_bool bUniform, gmx_ddbox_t *ddbox, real *cell_f)
2992 gmx_domdec_comm_t *comm;
2995 real range, pme_boundary;
2999 nc = dd->nc[ddpme->dim];
3002 if (!ddpme->dim_match)
3004 /* PP decomposition is not along dim: the worst situation */
3007 else if (ns <= 3 || (bUniform && ns == nc))
3009 /* The optimal situation */
3014 /* We need to check for all pme nodes which nodes they
3015 * could possibly need to communicate with.
3017 xmin = ddpme->pp_min;
3018 xmax = ddpme->pp_max;
3019 /* Allow for atoms to be maximally 2/3 times the cut-off
3020 * out of their DD cell. This is a reasonable balance between
3021 * between performance and support for most charge-group/cut-off
3024 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
3025 /* Avoid extra communication when we are exactly at a boundary */
3029 for (s = 0; s < ns; s++)
3031 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
3032 pme_boundary = (real)s/ns;
3035 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
3037 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
3041 pme_boundary = (real)(s+1)/ns;
3044 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
3046 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
3053 ddpme->maxshift = sh;
3057 fprintf(debug, "PME slab communication range for dim %d is %d\n",
3058 ddpme->dim, ddpme->maxshift);
3062 static void check_box_size(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3066 for (d = 0; d < dd->ndim; d++)
3069 if (dim < ddbox->nboundeddim &&
3070 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
3071 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3073 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",
3074 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3075 dd->nc[dim], dd->comm->cellsize_limit);
3080 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
3081 gmx_bool bMaster, ivec npulse)
3083 gmx_domdec_comm_t *comm;
3086 real *cell_x, cell_dx, cellsize;
3090 for (d = 0; d < DIM; d++)
3092 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3094 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3097 cell_dx = ddbox->box_size[d]/dd->nc[d];
3100 for (j = 0; j < dd->nc[d]+1; j++)
3102 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3107 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3108 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3110 cellsize = cell_dx*ddbox->skew_fac[d];
3111 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3115 cellsize_min[d] = cellsize;
3119 /* Statically load balanced grid */
3120 /* Also when we are not doing a master distribution we determine
3121 * all cell borders in a loop to obtain identical values
3122 * to the master distribution case and to determine npulse.
3126 cell_x = dd->ma->cell_x[d];
3130 snew(cell_x, dd->nc[d]+1);
3132 cell_x[0] = ddbox->box0[d];
3133 for (j = 0; j < dd->nc[d]; j++)
3135 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3136 cell_x[j+1] = cell_x[j] + cell_dx;
3137 cellsize = cell_dx*ddbox->skew_fac[d];
3138 while (cellsize*npulse[d] < comm->cutoff &&
3139 npulse[d] < dd->nc[d]-1)
3143 cellsize_min[d] = min(cellsize_min[d], cellsize);
3147 comm->cell_x0[d] = cell_x[dd->ci[d]];
3148 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3152 /* The following limitation is to avoid that a cell would receive
3153 * some of its own home charge groups back over the periodic boundary.
3154 * Double charge groups cause trouble with the global indices.
3156 if (d < ddbox->npbcdim &&
3157 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3159 gmx_fatal_collective(FARGS, NULL, dd,
3160 "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",
3161 dim2char(d), ddbox->box_size[d], ddbox->skew_fac[d],
3163 dd->nc[d], dd->nc[d],
3164 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3168 if (!comm->bDynLoadBal)
3170 copy_rvec(cellsize_min, comm->cellsize_min);
3173 for (d = 0; d < comm->npmedecompdim; d++)
3175 set_pme_maxshift(dd, &comm->ddpme[d],
3176 comm->slb_frac[dd->dim[d]] == NULL, ddbox,
3177 comm->ddpme[d].slb_dim_f);
3182 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3183 int d, int dim, gmx_domdec_root_t *root,
3185 gmx_bool bUniform, gmx_large_int_t step, real cellsize_limit_f, int range[])
3187 gmx_domdec_comm_t *comm;
3188 int ncd, i, j, nmin, nmin_old;
3189 gmx_bool bLimLo, bLimHi;
3191 real fac, halfway, cellsize_limit_f_i, region_size;
3192 gmx_bool bPBC, bLastHi = FALSE;
3193 int nrange[] = {range[0], range[1]};
3195 region_size = root->cell_f[range[1]]-root->cell_f[range[0]];
3201 bPBC = (dim < ddbox->npbcdim);
3203 cell_size = root->buf_ncd;
3207 fprintf(debug, "enforce_limits: %d %d\n", range[0], range[1]);
3210 /* First we need to check if the scaling does not make cells
3211 * smaller than the smallest allowed size.
3212 * We need to do this iteratively, since if a cell is too small,
3213 * it needs to be enlarged, which makes all the other cells smaller,
3214 * which could in turn make another cell smaller than allowed.
3216 for (i = range[0]; i < range[1]; i++)
3218 root->bCellMin[i] = FALSE;
3224 /* We need the total for normalization */
3226 for (i = range[0]; i < range[1]; i++)
3228 if (root->bCellMin[i] == FALSE)
3230 fac += cell_size[i];
3233 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3234 /* Determine the cell boundaries */
3235 for (i = range[0]; i < range[1]; i++)
3237 if (root->bCellMin[i] == FALSE)
3239 cell_size[i] *= fac;
3240 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3242 cellsize_limit_f_i = 0;
3246 cellsize_limit_f_i = cellsize_limit_f;
3248 if (cell_size[i] < cellsize_limit_f_i)
3250 root->bCellMin[i] = TRUE;
3251 cell_size[i] = cellsize_limit_f_i;
3255 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3258 while (nmin > nmin_old);
3261 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3262 /* For this check we should not use DD_CELL_MARGIN,
3263 * but a slightly smaller factor,
3264 * since rounding could get use below the limit.
3266 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3269 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",
3270 gmx_step_str(step, buf),
3271 dim2char(dim), ddbox->box_size[dim], ddbox->skew_fac[dim],
3272 ncd, comm->cellsize_min[dim]);
3275 root->bLimited = (nmin > 0) || (range[0] > 0) || (range[1] < ncd);
3279 /* Check if the boundary did not displace more than halfway
3280 * each of the cells it bounds, as this could cause problems,
3281 * especially when the differences between cell sizes are large.
3282 * If changes are applied, they will not make cells smaller
3283 * than the cut-off, as we check all the boundaries which
3284 * might be affected by a change and if the old state was ok,
3285 * the cells will at most be shrunk back to their old size.
3287 for (i = range[0]+1; i < range[1]; i++)
3289 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3290 if (root->cell_f[i] < halfway)
3292 root->cell_f[i] = halfway;
3293 /* Check if the change also causes shifts of the next boundaries */
3294 for (j = i+1; j < range[1]; j++)
3296 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3298 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3302 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3303 if (root->cell_f[i] > halfway)
3305 root->cell_f[i] = halfway;
3306 /* Check if the change also causes shifts of the next boundaries */
3307 for (j = i-1; j >= range[0]+1; j--)
3309 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3311 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3318 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3319 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3320 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3321 * for a and b nrange is used */
3324 /* Take care of the staggering of the cell boundaries */
3327 for (i = range[0]; i < range[1]; i++)
3329 root->cell_f_max0[i] = root->cell_f[i];
3330 root->cell_f_min1[i] = root->cell_f[i+1];
3335 for (i = range[0]+1; i < range[1]; i++)
3337 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3338 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3339 if (bLimLo && bLimHi)
3341 /* Both limits violated, try the best we can */
3342 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3343 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3344 nrange[0] = range[0];
3346 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3349 nrange[1] = range[1];
3350 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3356 /* root->cell_f[i] = root->bound_min[i]; */
3357 nrange[1] = i; /* only store violation location. There could be a LimLo violation following with an higher index */
3360 else if (bLimHi && !bLastHi)
3363 if (nrange[1] < range[1]) /* found a LimLo before */
3365 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3366 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3367 nrange[0] = nrange[1];
3369 root->cell_f[i] = root->bound_max[i];
3371 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3373 nrange[1] = range[1];
3376 if (nrange[1] < range[1]) /* found last a LimLo */
3378 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3379 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3380 nrange[0] = nrange[1];
3381 nrange[1] = range[1];
3382 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3384 else if (nrange[0] > range[0]) /* found at least one LimHi */
3386 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3393 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3394 int d, int dim, gmx_domdec_root_t *root,
3395 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3396 gmx_bool bUniform, gmx_large_int_t step)
3398 gmx_domdec_comm_t *comm;
3399 int ncd, d1, i, j, pos;
3401 real load_aver, load_i, imbalance, change, change_max, sc;
3402 real cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
3406 int range[] = { 0, 0 };
3410 /* Convert the maximum change from the input percentage to a fraction */
3411 change_limit = comm->dlb_scale_lim*0.01;
3415 bPBC = (dim < ddbox->npbcdim);
3417 cell_size = root->buf_ncd;
3419 /* Store the original boundaries */
3420 for (i = 0; i < ncd+1; i++)
3422 root->old_cell_f[i] = root->cell_f[i];
3426 for (i = 0; i < ncd; i++)
3428 cell_size[i] = 1.0/ncd;
3431 else if (dd_load_count(comm))
3433 load_aver = comm->load[d].sum_m/ncd;
3435 for (i = 0; i < ncd; i++)
3437 /* Determine the relative imbalance of cell i */
3438 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3439 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3440 /* Determine the change of the cell size using underrelaxation */
3441 change = -relax*imbalance;
3442 change_max = max(change_max, max(change, -change));
3444 /* Limit the amount of scaling.
3445 * We need to use the same rescaling for all cells in one row,
3446 * otherwise the load balancing might not converge.
3449 if (change_max > change_limit)
3451 sc *= change_limit/change_max;
3453 for (i = 0; i < ncd; i++)
3455 /* Determine the relative imbalance of cell i */
3456 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3457 imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
3458 /* Determine the change of the cell size using underrelaxation */
3459 change = -sc*imbalance;
3460 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3464 cellsize_limit_f = cellsize_min_dlb(comm, d, dim)/ddbox->box_size[dim];
3465 cellsize_limit_f *= DD_CELL_MARGIN;
3466 dist_min_f_hard = grid_jump_limit(comm, comm->cutoff, d)/ddbox->box_size[dim];
3467 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3468 if (ddbox->tric_dir[dim])
3470 cellsize_limit_f /= ddbox->skew_fac[dim];
3471 dist_min_f /= ddbox->skew_fac[dim];
3473 if (bDynamicBox && d > 0)
3475 dist_min_f *= DD_PRES_SCALE_MARGIN;
3477 if (d > 0 && !bUniform)
3479 /* Make sure that the grid is not shifted too much */
3480 for (i = 1; i < ncd; i++)
3482 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3484 gmx_incons("Inconsistent DD boundary staggering limits!");
3486 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3487 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3490 root->bound_min[i] += 0.5*space;
3492 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3493 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3496 root->bound_max[i] += 0.5*space;
3501 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3503 root->cell_f_max0[i-1] + dist_min_f,
3504 root->bound_min[i], root->cell_f[i], root->bound_max[i],
3505 root->cell_f_min1[i] - dist_min_f);
3510 root->cell_f[0] = 0;
3511 root->cell_f[ncd] = 1;
3512 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3515 /* After the checks above, the cells should obey the cut-off
3516 * restrictions, but it does not hurt to check.
3518 for (i = 0; i < ncd; i++)
3522 fprintf(debug, "Relative bounds dim %d cell %d: %f %f\n",
3523 dim, i, root->cell_f[i], root->cell_f[i+1]);
3526 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3527 root->cell_f[i+1] - root->cell_f[i] <
3528 cellsize_limit_f/DD_CELL_MARGIN)
3532 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3533 gmx_step_str(step, buf), dim2char(dim), i,
3534 (root->cell_f[i+1] - root->cell_f[i])
3535 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3540 /* Store the cell boundaries of the lower dimensions at the end */
3541 for (d1 = 0; d1 < d; d1++)
3543 root->cell_f[pos++] = comm->cell_f0[d1];
3544 root->cell_f[pos++] = comm->cell_f1[d1];
3547 if (d < comm->npmedecompdim)
3549 /* The master determines the maximum shift for
3550 * the coordinate communication between separate PME nodes.
3552 set_pme_maxshift(dd, &comm->ddpme[d], bUniform, ddbox, root->cell_f);
3554 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3557 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3561 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3562 gmx_ddbox_t *ddbox, int dimind)
3564 gmx_domdec_comm_t *comm;
3569 /* Set the cell dimensions */
3570 dim = dd->dim[dimind];
3571 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3572 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3573 if (dim >= ddbox->nboundeddim)
3575 comm->cell_x0[dim] += ddbox->box0[dim];
3576 comm->cell_x1[dim] += ddbox->box0[dim];
3580 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3581 int d, int dim, real *cell_f_row,
3584 gmx_domdec_comm_t *comm;
3590 /* Each node would only need to know two fractions,
3591 * but it is probably cheaper to broadcast the whole array.
3593 MPI_Bcast(cell_f_row, DD_CELL_F_SIZE(dd, d)*sizeof(real), MPI_BYTE,
3594 0, comm->mpi_comm_load[d]);
3596 /* Copy the fractions for this dimension from the buffer */
3597 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3598 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3599 /* The whole array was communicated, so set the buffer position */
3600 pos = dd->nc[dim] + 1;
3601 for (d1 = 0; d1 <= d; d1++)
3605 /* Copy the cell fractions of the lower dimensions */
3606 comm->cell_f0[d1] = cell_f_row[pos++];
3607 comm->cell_f1[d1] = cell_f_row[pos++];
3609 relative_to_absolute_cell_bounds(dd, ddbox, d1);
3611 /* Convert the communicated shift from float to int */
3612 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3615 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3619 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3620 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3621 gmx_bool bUniform, gmx_large_int_t step)
3623 gmx_domdec_comm_t *comm;
3625 gmx_bool bRowMember, bRowRoot;
3630 for (d = 0; d < dd->ndim; d++)
3635 for (d1 = d; d1 < dd->ndim; d1++)
3637 if (dd->ci[dd->dim[d1]] > 0)
3650 set_dd_cell_sizes_dlb_root(dd, d, dim, comm->root[d],
3651 ddbox, bDynamicBox, bUniform, step);
3652 cell_f_row = comm->root[d]->cell_f;
3656 cell_f_row = comm->cell_f_row;
3658 distribute_dd_cell_sizes_dlb(dd, d, dim, cell_f_row, ddbox);
3663 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3667 /* This function assumes the box is static and should therefore
3668 * not be called when the box has changed since the last
3669 * call to dd_partition_system.
3671 for (d = 0; d < dd->ndim; d++)
3673 relative_to_absolute_cell_bounds(dd, ddbox, d);
3679 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3680 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3681 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3682 gmx_wallcycle_t wcycle)
3684 gmx_domdec_comm_t *comm;
3691 wallcycle_start(wcycle, ewcDDCOMMBOUND);
3692 set_dd_cell_sizes_dlb_change(dd, ddbox, bDynamicBox, bUniform, step);
3693 wallcycle_stop(wcycle, ewcDDCOMMBOUND);
3695 else if (bDynamicBox)
3697 set_dd_cell_sizes_dlb_nochange(dd, ddbox);
3700 /* Set the dimensions for which no DD is used */
3701 for (dim = 0; dim < DIM; dim++)
3703 if (dd->nc[dim] == 1)
3705 comm->cell_x0[dim] = 0;
3706 comm->cell_x1[dim] = ddbox->box_size[dim];
3707 if (dim >= ddbox->nboundeddim)
3709 comm->cell_x0[dim] += ddbox->box0[dim];
3710 comm->cell_x1[dim] += ddbox->box0[dim];
3716 static void realloc_comm_ind(gmx_domdec_t *dd, ivec npulse)
3719 gmx_domdec_comm_dim_t *cd;
3721 for (d = 0; d < dd->ndim; d++)
3723 cd = &dd->comm->cd[d];
3724 np = npulse[dd->dim[d]];
3725 if (np > cd->np_nalloc)
3729 fprintf(debug, "(Re)allocing cd for %c to %d pulses\n",
3730 dim2char(dd->dim[d]), np);
3732 if (DDMASTER(dd) && cd->np_nalloc > 0)
3734 fprintf(stderr, "\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n", dim2char(dd->dim[d]), np);
3736 srenew(cd->ind, np);
3737 for (i = cd->np_nalloc; i < np; i++)
3739 cd->ind[i].index = NULL;
3740 cd->ind[i].nalloc = 0;
3749 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3750 gmx_ddbox_t *ddbox, gmx_bool bDynamicBox,
3751 gmx_bool bUniform, gmx_bool bDoDLB, gmx_large_int_t step,
3752 gmx_wallcycle_t wcycle)
3754 gmx_domdec_comm_t *comm;
3760 /* Copy the old cell boundaries for the cg displacement check */
3761 copy_rvec(comm->cell_x0, comm->old_cell_x0);
3762 copy_rvec(comm->cell_x1, comm->old_cell_x1);
3764 if (comm->bDynLoadBal)
3768 check_box_size(dd, ddbox);
3770 set_dd_cell_sizes_dlb(dd, ddbox, bDynamicBox, bUniform, bDoDLB, step, wcycle);
3774 set_dd_cell_sizes_slb(dd, ddbox, FALSE, npulse);
3775 realloc_comm_ind(dd, npulse);
3780 for (d = 0; d < DIM; d++)
3782 fprintf(debug, "cell_x[%d] %f - %f skew_fac %f\n",
3783 d, comm->cell_x0[d], comm->cell_x1[d], ddbox->skew_fac[d]);
3788 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3790 rvec cell_ns_x0, rvec cell_ns_x1,
3791 gmx_large_int_t step)
3793 gmx_domdec_comm_t *comm;
3798 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
3800 dim = dd->dim[dim_ind];
3802 /* Without PBC we don't have restrictions on the outer cells */
3803 if (!(dim >= ddbox->npbcdim &&
3804 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3805 comm->bDynLoadBal &&
3806 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3807 comm->cellsize_min[dim])
3810 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",
3811 gmx_step_str(step, buf), dim2char(dim),
3812 comm->cell_x1[dim] - comm->cell_x0[dim],
3813 ddbox->skew_fac[dim],
3814 dd->comm->cellsize_min[dim],
3815 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3819 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3821 /* Communicate the boundaries and update cell_ns_x0/1 */
3822 dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
3823 if (dd->bGridJump && dd->ndim > 1)
3825 check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
3830 static void make_tric_corr_matrix(int npbcdim, matrix box, matrix tcm)
3834 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3842 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3843 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3852 static void check_screw_box(matrix box)
3854 /* Mathematical limitation */
3855 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3857 gmx_fatal(FARGS, "With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3860 /* Limitation due to the asymmetry of the eighth shell method */
3861 if (box[ZZ][YY] != 0)
3863 gmx_fatal(FARGS, "pbc=screw with non-zero box_zy is not supported");
3867 static void distribute_cg(FILE *fplog, gmx_large_int_t step,
3868 matrix box, ivec tric_dir, t_block *cgs, rvec pos[],
3871 gmx_domdec_master_t *ma;
3872 int **tmp_ind = NULL, *tmp_nalloc = NULL;
3873 int i, icg, j, k, k0, k1, d, npbcdim;
3875 rvec box_size, cg_cm;
3877 real nrcg, inv_ncg, pos_d;
3879 gmx_bool bUnbounded, bScrew;
3883 if (tmp_ind == NULL)
3885 snew(tmp_nalloc, dd->nnodes);
3886 snew(tmp_ind, dd->nnodes);
3887 for (i = 0; i < dd->nnodes; i++)
3889 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3890 snew(tmp_ind[i], tmp_nalloc[i]);
3894 /* Clear the count */
3895 for (i = 0; i < dd->nnodes; i++)
3901 make_tric_corr_matrix(dd->npbcdim, box, tcm);
3903 cgindex = cgs->index;
3905 /* Compute the center of geometry for all charge groups */
3906 for (icg = 0; icg < cgs->nr; icg++)
3909 k1 = cgindex[icg+1];
3913 copy_rvec(pos[k0], cg_cm);
3920 for (k = k0; (k < k1); k++)
3922 rvec_inc(cg_cm, pos[k]);
3924 for (d = 0; (d < DIM); d++)
3926 cg_cm[d] *= inv_ncg;
3929 /* Put the charge group in the box and determine the cell index */
3930 for (d = DIM-1; d >= 0; d--)
3933 if (d < dd->npbcdim)
3935 bScrew = (dd->bScrewPBC && d == XX);
3936 if (tric_dir[d] && dd->nc[d] > 1)
3938 /* Use triclinic coordintates for this dimension */
3939 for (j = d+1; j < DIM; j++)
3941 pos_d += cg_cm[j]*tcm[j][d];
3944 while (pos_d >= box[d][d])
3947 rvec_dec(cg_cm, box[d]);
3950 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3951 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3953 for (k = k0; (k < k1); k++)
3955 rvec_dec(pos[k], box[d]);
3958 pos[k][YY] = box[YY][YY] - pos[k][YY];
3959 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3966 rvec_inc(cg_cm, box[d]);
3969 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3970 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3972 for (k = k0; (k < k1); k++)
3974 rvec_inc(pos[k], box[d]);
3977 pos[k][YY] = box[YY][YY] - pos[k][YY];
3978 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3983 /* This could be done more efficiently */
3985 while (ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3990 i = dd_index(dd->nc, ind);
3991 if (ma->ncg[i] == tmp_nalloc[i])
3993 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3994 srenew(tmp_ind[i], tmp_nalloc[i]);
3996 tmp_ind[i][ma->ncg[i]] = icg;
3998 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
4002 for (i = 0; i < dd->nnodes; i++)
4005 for (k = 0; k < ma->ncg[i]; k++)
4007 ma->cg[k1++] = tmp_ind[i][k];
4010 ma->index[dd->nnodes] = k1;
4012 for (i = 0; i < dd->nnodes; i++)
4022 fprintf(fplog, "Charge group distribution at step %s:",
4023 gmx_step_str(step, buf));
4024 for (i = 0; i < dd->nnodes; i++)
4026 fprintf(fplog, " %d", ma->ncg[i]);
4028 fprintf(fplog, "\n");
4032 static void get_cg_distribution(FILE *fplog, gmx_large_int_t step, gmx_domdec_t *dd,
4033 t_block *cgs, matrix box, gmx_ddbox_t *ddbox,
4036 gmx_domdec_master_t *ma = NULL;
4039 int *ibuf, buf2[2] = { 0, 0 };
4040 gmx_bool bMaster = DDMASTER(dd);
4047 check_screw_box(box);
4050 set_dd_cell_sizes_slb(dd, ddbox, TRUE, npulse);
4052 distribute_cg(fplog, step, box, ddbox->tric_dir, cgs, pos, dd);
4053 for (i = 0; i < dd->nnodes; i++)
4055 ma->ibuf[2*i] = ma->ncg[i];
4056 ma->ibuf[2*i+1] = ma->nat[i];
4064 dd_scatter(dd, 2*sizeof(int), ibuf, buf2);
4066 dd->ncg_home = buf2[0];
4067 dd->nat_home = buf2[1];
4068 dd->ncg_tot = dd->ncg_home;
4069 dd->nat_tot = dd->nat_home;
4070 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
4072 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
4073 srenew(dd->index_gl, dd->cg_nalloc);
4074 srenew(dd->cgindex, dd->cg_nalloc+1);
4078 for (i = 0; i < dd->nnodes; i++)
4080 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
4081 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4086 DDMASTER(dd) ? ma->ibuf : NULL,
4087 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4088 DDMASTER(dd) ? ma->cg : NULL,
4089 dd->ncg_home*sizeof(int), dd->index_gl);
4091 /* Determine the home charge group sizes */
4093 for (i = 0; i < dd->ncg_home; i++)
4095 cg_gl = dd->index_gl[i];
4097 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4102 fprintf(debug, "Home charge groups:\n");
4103 for (i = 0; i < dd->ncg_home; i++)
4105 fprintf(debug, " %d", dd->index_gl[i]);
4108 fprintf(debug, "\n");
4111 fprintf(debug, "\n");
4115 static int compact_and_copy_vec_at(int ncg, int *move,
4118 rvec *src, gmx_domdec_comm_t *comm,
4121 int m, icg, i, i0, i1, nrcg;
4127 for (m = 0; m < DIM*2; m++)
4133 for (icg = 0; icg < ncg; icg++)
4135 i1 = cgindex[icg+1];
4141 /* Compact the home array in place */
4142 for (i = i0; i < i1; i++)
4144 copy_rvec(src[i], src[home_pos++]);
4150 /* Copy to the communication buffer */
4152 pos_vec[m] += 1 + vec*nrcg;
4153 for (i = i0; i < i1; i++)
4155 copy_rvec(src[i], comm->cgcm_state[m][pos_vec[m]++]);
4157 pos_vec[m] += (nvec - vec - 1)*nrcg;
4161 home_pos += i1 - i0;
4169 static int compact_and_copy_vec_cg(int ncg, int *move,
4171 int nvec, rvec *src, gmx_domdec_comm_t *comm,
4174 int m, icg, i0, i1, nrcg;
4180 for (m = 0; m < DIM*2; m++)
4186 for (icg = 0; icg < ncg; icg++)
4188 i1 = cgindex[icg+1];
4194 /* Compact the home array in place */
4195 copy_rvec(src[icg], src[home_pos++]);
4201 /* Copy to the communication buffer */
4202 copy_rvec(src[icg], comm->cgcm_state[m][pos_vec[m]]);
4203 pos_vec[m] += 1 + nrcg*nvec;
4215 static int compact_ind(int ncg, int *move,
4216 int *index_gl, int *cgindex,
4218 gmx_ga2la_t ga2la, char *bLocalCG,
4221 int cg, nat, a0, a1, a, a_gl;
4226 for (cg = 0; cg < ncg; cg++)
4232 /* Compact the home arrays in place.
4233 * Anything that can be done here avoids access to global arrays.
4235 cgindex[home_pos] = nat;
4236 for (a = a0; a < a1; a++)
4239 gatindex[nat] = a_gl;
4240 /* The cell number stays 0, so we don't need to set it */
4241 ga2la_change_la(ga2la, a_gl, nat);
4244 index_gl[home_pos] = index_gl[cg];
4245 cginfo[home_pos] = cginfo[cg];
4246 /* The charge group remains local, so bLocalCG does not change */
4251 /* Clear the global indices */
4252 for (a = a0; a < a1; a++)
4254 ga2la_del(ga2la, gatindex[a]);
4258 bLocalCG[index_gl[cg]] = FALSE;
4262 cgindex[home_pos] = nat;
4267 static void clear_and_mark_ind(int ncg, int *move,
4268 int *index_gl, int *cgindex, int *gatindex,
4269 gmx_ga2la_t ga2la, char *bLocalCG,
4274 for (cg = 0; cg < ncg; cg++)
4280 /* Clear the global indices */
4281 for (a = a0; a < a1; a++)
4283 ga2la_del(ga2la, gatindex[a]);
4287 bLocalCG[index_gl[cg]] = FALSE;
4289 /* Signal that this cg has moved using the ns cell index.
4290 * Here we set it to -1. fill_grid will change it
4291 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4293 cell_index[cg] = -1;
4298 static void print_cg_move(FILE *fplog,
4300 gmx_large_int_t step, int cg, int dim, int dir,
4301 gmx_bool bHaveLimitdAndCMOld, real limitd,
4302 rvec cm_old, rvec cm_new, real pos_d)
4304 gmx_domdec_comm_t *comm;
4309 fprintf(fplog, "\nStep %s:\n", gmx_step_str(step, buf));
4310 if (bHaveLimitdAndCMOld)
4312 fprintf(fplog, "The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4313 ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
4317 fprintf(fplog, "The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4318 ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
4320 fprintf(fplog, "distance out of cell %f\n",
4321 dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4322 if (bHaveLimitdAndCMOld)
4324 fprintf(fplog, "Old coordinates: %8.3f %8.3f %8.3f\n",
4325 cm_old[XX], cm_old[YY], cm_old[ZZ]);
4327 fprintf(fplog, "New coordinates: %8.3f %8.3f %8.3f\n",
4328 cm_new[XX], cm_new[YY], cm_new[ZZ]);
4329 fprintf(fplog, "Old cell boundaries in direction %c: %8.3f %8.3f\n",
4331 comm->old_cell_x0[dim], comm->old_cell_x1[dim]);
4332 fprintf(fplog, "New cell boundaries in direction %c: %8.3f %8.3f\n",
4334 comm->cell_x0[dim], comm->cell_x1[dim]);
4337 static void cg_move_error(FILE *fplog,
4339 gmx_large_int_t step, int cg, int dim, int dir,
4340 gmx_bool bHaveLimitdAndCMOld, real limitd,
4341 rvec cm_old, rvec cm_new, real pos_d)
4345 print_cg_move(fplog, dd, step, cg, dim, dir,
4346 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4348 print_cg_move(stderr, dd, step, cg, dim, dir,
4349 bHaveLimitdAndCMOld, limitd, cm_old, cm_new, pos_d);
4351 "A charge group moved too far between two domain decomposition steps\n"
4352 "This usually means that your system is not well equilibrated");
4355 static void rotate_state_atom(t_state *state, int a)
4359 for (est = 0; est < estNR; est++)
4361 if (EST_DISTR(est) && (state->flags & (1<<est)))
4366 /* Rotate the complete state; for a rectangular box only */
4367 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4368 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4371 state->v[a][YY] = -state->v[a][YY];
4372 state->v[a][ZZ] = -state->v[a][ZZ];
4375 state->sd_X[a][YY] = -state->sd_X[a][YY];
4376 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4379 state->cg_p[a][YY] = -state->cg_p[a][YY];
4380 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4382 case estDISRE_INITF:
4383 case estDISRE_RM3TAV:
4384 case estORIRE_INITF:
4386 /* These are distances, so not affected by rotation */
4389 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4395 static int *get_moved(gmx_domdec_comm_t *comm, int natoms)
4397 if (natoms > comm->moved_nalloc)
4399 /* Contents should be preserved here */
4400 comm->moved_nalloc = over_alloc_dd(natoms);
4401 srenew(comm->moved, comm->moved_nalloc);
4407 static void calc_cg_move(FILE *fplog, gmx_large_int_t step,
4410 ivec tric_dir, matrix tcm,
4411 rvec cell_x0, rvec cell_x1,
4412 rvec limitd, rvec limit0, rvec limit1,
4414 int cg_start, int cg_end,
4419 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4420 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4424 real inv_ncg, pos_d;
4427 npbcdim = dd->npbcdim;
4429 for (cg = cg_start; cg < cg_end; cg++)
4436 copy_rvec(state->x[k0], cm_new);
4443 for (k = k0; (k < k1); k++)
4445 rvec_inc(cm_new, state->x[k]);
4447 for (d = 0; (d < DIM); d++)
4449 cm_new[d] = inv_ncg*cm_new[d];
4454 /* Do pbc and check DD cell boundary crossings */
4455 for (d = DIM-1; d >= 0; d--)
4459 bScrew = (dd->bScrewPBC && d == XX);
4460 /* Determine the location of this cg in lattice coordinates */
4464 for (d2 = d+1; d2 < DIM; d2++)
4466 pos_d += cm_new[d2]*tcm[d2][d];
4469 /* Put the charge group in the triclinic unit-cell */
4470 if (pos_d >= cell_x1[d])
4472 if (pos_d >= limit1[d])
4474 cg_move_error(fplog, dd, step, cg, d, 1, TRUE, limitd[d],
4475 cg_cm[cg], cm_new, pos_d);
4478 if (dd->ci[d] == dd->nc[d] - 1)
4480 rvec_dec(cm_new, state->box[d]);
4483 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4484 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4486 for (k = k0; (k < k1); k++)
4488 rvec_dec(state->x[k], state->box[d]);
4491 rotate_state_atom(state, k);
4496 else if (pos_d < cell_x0[d])
4498 if (pos_d < limit0[d])
4500 cg_move_error(fplog, dd, step, cg, d, -1, TRUE, limitd[d],
4501 cg_cm[cg], cm_new, pos_d);
4506 rvec_inc(cm_new, state->box[d]);
4509 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4510 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4512 for (k = k0; (k < k1); k++)
4514 rvec_inc(state->x[k], state->box[d]);
4517 rotate_state_atom(state, k);
4523 else if (d < npbcdim)
4525 /* Put the charge group in the rectangular unit-cell */
4526 while (cm_new[d] >= state->box[d][d])
4528 rvec_dec(cm_new, state->box[d]);
4529 for (k = k0; (k < k1); k++)
4531 rvec_dec(state->x[k], state->box[d]);
4534 while (cm_new[d] < 0)
4536 rvec_inc(cm_new, state->box[d]);
4537 for (k = k0; (k < k1); k++)
4539 rvec_inc(state->x[k], state->box[d]);
4545 copy_rvec(cm_new, cg_cm[cg]);
4547 /* Determine where this cg should go */
4550 for (d = 0; d < dd->ndim; d++)
4555 flag |= DD_FLAG_FW(d);
4561 else if (dev[dim] == -1)
4563 flag |= DD_FLAG_BW(d);
4566 if (dd->nc[dim] > 2)
4577 /* Temporarily store the flag in move */
4578 move[cg] = mc + flag;
4582 static void dd_redistribute_cg(FILE *fplog, gmx_large_int_t step,
4583 gmx_domdec_t *dd, ivec tric_dir,
4584 t_state *state, rvec **f,
4585 t_forcerec *fr, t_mdatoms *md,
4593 int ncg[DIM*2], nat[DIM*2];
4594 int c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
4595 int mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
4596 int sbuf[2], rbuf[2];
4597 int home_pos_cg, home_pos_at, buf_pos;
4599 gmx_bool bV = FALSE, bSDX = FALSE, bCGP = FALSE;
4602 real inv_ncg, pos_d;
4604 rvec *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
4606 cginfo_mb_t *cginfo_mb;
4607 gmx_domdec_comm_t *comm;
4609 int nthread, thread;
4613 check_screw_box(state->box);
4617 if (fr->cutoff_scheme == ecutsGROUP)
4622 for (i = 0; i < estNR; i++)
4628 case estX: /* Always present */ break;
4629 case estV: bV = (state->flags & (1<<i)); break;
4630 case estSDX: bSDX = (state->flags & (1<<i)); break;
4631 case estCGP: bCGP = (state->flags & (1<<i)); break;
4634 case estDISRE_INITF:
4635 case estDISRE_RM3TAV:
4636 case estORIRE_INITF:
4638 /* No processing required */
4641 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4646 if (dd->ncg_tot > comm->nalloc_int)
4648 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4649 srenew(comm->buf_int, comm->nalloc_int);
4651 move = comm->buf_int;
4653 /* Clear the count */
4654 for (c = 0; c < dd->ndim*2; c++)
4660 npbcdim = dd->npbcdim;
4662 for (d = 0; (d < DIM); d++)
4664 limitd[d] = dd->comm->cellsize_min[d];
4665 if (d >= npbcdim && dd->ci[d] == 0)
4667 cell_x0[d] = -GMX_FLOAT_MAX;
4671 cell_x0[d] = comm->cell_x0[d];
4673 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4675 cell_x1[d] = GMX_FLOAT_MAX;
4679 cell_x1[d] = comm->cell_x1[d];
4683 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4684 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4688 /* We check after communication if a charge group moved
4689 * more than one cell. Set the pre-comm check limit to float_max.
4691 limit0[d] = -GMX_FLOAT_MAX;
4692 limit1[d] = GMX_FLOAT_MAX;
4696 make_tric_corr_matrix(npbcdim, state->box, tcm);
4698 cgindex = dd->cgindex;
4700 nthread = gmx_omp_nthreads_get(emntDomdec);
4702 /* Compute the center of geometry for all home charge groups
4703 * and put them in the box and determine where they should go.
4705 #pragma omp parallel for num_threads(nthread) schedule(static)
4706 for (thread = 0; thread < nthread; thread++)
4708 calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
4709 cell_x0, cell_x1, limitd, limit0, limit1,
4711 ( thread *dd->ncg_home)/nthread,
4712 ((thread+1)*dd->ncg_home)/nthread,
4713 fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
4717 for (cg = 0; cg < dd->ncg_home; cg++)
4722 flag = mc & ~DD_FLAG_NRCG;
4723 mc = mc & DD_FLAG_NRCG;
4726 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4728 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4729 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4731 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4732 /* We store the cg size in the lower 16 bits
4733 * and the place where the charge group should go
4734 * in the next 6 bits. This saves some communication volume.
4736 nrcg = cgindex[cg+1] - cgindex[cg];
4737 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4743 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
4744 inc_nrnb(nrnb, eNR_RESETX, dd->ncg_home);
4747 for (i = 0; i < dd->ndim*2; i++)
4749 *ncg_moved += ncg[i];
4766 /* Make sure the communication buffers are large enough */
4767 for (mc = 0; mc < dd->ndim*2; mc++)
4769 nvr = ncg[mc] + nat[mc]*nvec;
4770 if (nvr > comm->cgcm_state_nalloc[mc])
4772 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4773 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
4777 switch (fr->cutoff_scheme)
4780 /* Recalculating cg_cm might be cheaper than communicating,
4781 * but that could give rise to rounding issues.
4784 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4785 nvec, cg_cm, comm, bCompact);
4788 /* Without charge groups we send the moved atom coordinates
4789 * over twice. This is so the code below can be used without
4790 * many conditionals for both for with and without charge groups.
4793 compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
4794 nvec, state->x, comm, FALSE);
4797 home_pos_cg -= *ncg_moved;
4801 gmx_incons("unimplemented");
4807 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4808 nvec, vec++, state->x, comm, bCompact);
4811 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4812 nvec, vec++, state->v, comm, bCompact);
4816 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4817 nvec, vec++, state->sd_X, comm, bCompact);
4821 compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
4822 nvec, vec++, state->cg_p, comm, bCompact);
4827 compact_ind(dd->ncg_home, move,
4828 dd->index_gl, dd->cgindex, dd->gatindex,
4829 dd->ga2la, comm->bLocalCG,
4834 if (fr->cutoff_scheme == ecutsVERLET)
4836 moved = get_moved(comm, dd->ncg_home);
4838 for (k = 0; k < dd->ncg_home; k++)
4845 moved = fr->ns.grid->cell_index;
4848 clear_and_mark_ind(dd->ncg_home, move,
4849 dd->index_gl, dd->cgindex, dd->gatindex,
4850 dd->ga2la, comm->bLocalCG,
4854 cginfo_mb = fr->cginfo_mb;
4856 *ncg_stay_home = home_pos_cg;
4857 for (d = 0; d < dd->ndim; d++)
4863 for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
4866 /* Communicate the cg and atom counts */
4871 fprintf(debug, "Sending ddim %d dir %d: ncg %d nat %d\n",
4872 d, dir, sbuf[0], sbuf[1]);
4874 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4876 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4878 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4879 srenew(comm->buf_int, comm->nalloc_int);
4882 /* Communicate the charge group indices, sizes and flags */
4883 dd_sendrecv_int(dd, d, dir,
4884 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4885 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4887 nvs = ncg[cdd] + nat[cdd]*nvec;
4888 i = rbuf[0] + rbuf[1] *nvec;
4889 vec_rvec_check_alloc(&comm->vbuf, nvr+i);
4891 /* Communicate cgcm and state */
4892 dd_sendrecv_rvec(dd, d, dir,
4893 comm->cgcm_state[cdd], nvs,
4894 comm->vbuf.v+nvr, i);
4895 ncg_recv += rbuf[0];
4896 nat_recv += rbuf[1];
4900 /* Process the received charge groups */
4902 for (cg = 0; cg < ncg_recv; cg++)
4904 flag = comm->buf_int[cg*DD_CGIBS+1];
4906 if (dim >= npbcdim && dd->nc[dim] > 2)
4908 /* No pbc in this dim and more than one domain boundary.
4909 * We do a separate check if a charge group didn't move too far.
4911 if (((flag & DD_FLAG_FW(d)) &&
4912 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4913 ((flag & DD_FLAG_BW(d)) &&
4914 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4916 cg_move_error(fplog, dd, step, cg, dim,
4917 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4919 comm->vbuf.v[buf_pos],
4920 comm->vbuf.v[buf_pos],
4921 comm->vbuf.v[buf_pos][dim]);
4928 /* Check which direction this cg should go */
4929 for (d2 = d+1; (d2 < dd->ndim && mc == -1); d2++)
4933 /* The cell boundaries for dimension d2 are not equal
4934 * for each cell row of the lower dimension(s),
4935 * therefore we might need to redetermine where
4936 * this cg should go.
4939 /* If this cg crosses the box boundary in dimension d2
4940 * we can use the communicated flag, so we do not
4941 * have to worry about pbc.
4943 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4944 (flag & DD_FLAG_FW(d2))) ||
4945 (dd->ci[dim2] == 0 &&
4946 (flag & DD_FLAG_BW(d2)))))
4948 /* Clear the two flags for this dimension */
4949 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4950 /* Determine the location of this cg
4951 * in lattice coordinates
4953 pos_d = comm->vbuf.v[buf_pos][dim2];
4956 for (d3 = dim2+1; d3 < DIM; d3++)
4959 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4962 /* Check of we are not at the box edge.
4963 * pbc is only handled in the first step above,
4964 * but this check could move over pbc while
4965 * the first step did not due to different rounding.
4967 if (pos_d >= cell_x1[dim2] &&
4968 dd->ci[dim2] != dd->nc[dim2]-1)
4970 flag |= DD_FLAG_FW(d2);
4972 else if (pos_d < cell_x0[dim2] &&
4975 flag |= DD_FLAG_BW(d2);
4977 comm->buf_int[cg*DD_CGIBS+1] = flag;
4980 /* Set to which neighboring cell this cg should go */
4981 if (flag & DD_FLAG_FW(d2))
4985 else if (flag & DD_FLAG_BW(d2))
4987 if (dd->nc[dd->dim[d2]] > 2)
4999 nrcg = flag & DD_FLAG_NRCG;
5002 if (home_pos_cg+1 > dd->cg_nalloc)
5004 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
5005 srenew(dd->index_gl, dd->cg_nalloc);
5006 srenew(dd->cgindex, dd->cg_nalloc+1);
5008 /* Set the global charge group index and size */
5009 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
5010 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
5011 /* Copy the state from the buffer */
5012 dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
5013 if (fr->cutoff_scheme == ecutsGROUP)
5016 copy_rvec(comm->vbuf.v[buf_pos], cg_cm[home_pos_cg]);
5020 /* Set the cginfo */
5021 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
5022 dd->index_gl[home_pos_cg]);
5025 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
5028 if (home_pos_at+nrcg > state->nalloc)
5030 dd_realloc_state(state, f, home_pos_at+nrcg);
5032 for (i = 0; i < nrcg; i++)
5034 copy_rvec(comm->vbuf.v[buf_pos++],
5035 state->x[home_pos_at+i]);
5039 for (i = 0; i < nrcg; i++)
5041 copy_rvec(comm->vbuf.v[buf_pos++],
5042 state->v[home_pos_at+i]);
5047 for (i = 0; i < nrcg; i++)
5049 copy_rvec(comm->vbuf.v[buf_pos++],
5050 state->sd_X[home_pos_at+i]);
5055 for (i = 0; i < nrcg; i++)
5057 copy_rvec(comm->vbuf.v[buf_pos++],
5058 state->cg_p[home_pos_at+i]);
5062 home_pos_at += nrcg;
5066 /* Reallocate the buffers if necessary */
5067 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
5069 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
5070 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
5072 nvr = ncg[mc] + nat[mc]*nvec;
5073 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
5075 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
5076 srenew(comm->cgcm_state[mc], comm->cgcm_state_nalloc[mc]);
5078 /* Copy from the receive to the send buffers */
5079 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
5080 comm->buf_int + cg*DD_CGIBS,
5081 DD_CGIBS*sizeof(int));
5082 memcpy(comm->cgcm_state[mc][nvr],
5083 comm->vbuf.v[buf_pos],
5084 (1+nrcg*nvec)*sizeof(rvec));
5085 buf_pos += 1 + nrcg*nvec;
5092 /* With sorting (!bCompact) the indices are now only partially up to date
5093 * and ncg_home and nat_home are not the real count, since there are
5094 * "holes" in the arrays for the charge groups that moved to neighbors.
5096 if (fr->cutoff_scheme == ecutsVERLET)
5098 moved = get_moved(comm, home_pos_cg);
5100 for (i = dd->ncg_home; i < home_pos_cg; i++)
5105 dd->ncg_home = home_pos_cg;
5106 dd->nat_home = home_pos_at;
5111 "Finished repartitioning: cgs moved out %d, new home %d\n",
5112 *ncg_moved, dd->ncg_home-*ncg_moved);
5117 void dd_cycles_add(gmx_domdec_t *dd, float cycles, int ddCycl)
5119 dd->comm->cycl[ddCycl] += cycles;
5120 dd->comm->cycl_n[ddCycl]++;
5121 if (cycles > dd->comm->cycl_max[ddCycl])
5123 dd->comm->cycl_max[ddCycl] = cycles;
5127 static double force_flop_count(t_nrnb *nrnb)
5134 for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
5136 /* To get closer to the real timings, we half the count
5137 * for the normal loops and again half it for water loops.
5140 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5142 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5146 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5149 for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
5152 if (strstr(name, "W3") != NULL || strstr(name, "W4") != NULL)
5154 sum += nrnb->n[i]*cost_nrnb(i);
5157 for (i = eNR_BONDS; i <= eNR_WALLS; i++)
5159 sum += nrnb->n[i]*cost_nrnb(i);
5165 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
5167 if (dd->comm->eFlop)
5169 dd->comm->flop -= force_flop_count(nrnb);
5172 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
5174 if (dd->comm->eFlop)
5176 dd->comm->flop += force_flop_count(nrnb);
5181 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5185 for (i = 0; i < ddCyclNr; i++)
5187 dd->comm->cycl[i] = 0;
5188 dd->comm->cycl_n[i] = 0;
5189 dd->comm->cycl_max[i] = 0;
5192 dd->comm->flop_n = 0;
5195 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
5197 gmx_domdec_comm_t *comm;
5198 gmx_domdec_load_t *load;
5199 gmx_domdec_root_t *root = NULL;
5200 int d, dim, cid, i, pos;
5201 float cell_frac = 0, sbuf[DD_NLOAD_MAX];
5206 fprintf(debug, "get_load_distribution start\n");
5209 wallcycle_start(wcycle, ewcDDCOMMLOAD);
5213 bSepPME = (dd->pme_nodeid >= 0);
5215 for (d = dd->ndim-1; d >= 0; d--)
5218 /* Check if we participate in the communication in this dimension */
5219 if (d == dd->ndim-1 ||
5220 (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
5222 load = &comm->load[d];
5225 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5228 if (d == dd->ndim-1)
5230 sbuf[pos++] = dd_force_load(comm);
5231 sbuf[pos++] = sbuf[0];
5234 sbuf[pos++] = sbuf[0];
5235 sbuf[pos++] = cell_frac;
5238 sbuf[pos++] = comm->cell_f_max0[d];
5239 sbuf[pos++] = comm->cell_f_min1[d];
5244 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5245 sbuf[pos++] = comm->cycl[ddCyclPME];
5250 sbuf[pos++] = comm->load[d+1].sum;
5251 sbuf[pos++] = comm->load[d+1].max;
5254 sbuf[pos++] = comm->load[d+1].sum_m;
5255 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5256 sbuf[pos++] = comm->load[d+1].flags;
5259 sbuf[pos++] = comm->cell_f_max0[d];
5260 sbuf[pos++] = comm->cell_f_min1[d];
5265 sbuf[pos++] = comm->load[d+1].mdf;
5266 sbuf[pos++] = comm->load[d+1].pme;
5270 /* Communicate a row in DD direction d.
5271 * The communicators are setup such that the root always has rank 0.
5274 MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
5275 load->load, load->nload*sizeof(float), MPI_BYTE,
5276 0, comm->mpi_comm_load[d]);
5278 if (dd->ci[dim] == dd->master_ci[dim])
5280 /* We are the root, process this row */
5281 if (comm->bDynLoadBal)
5283 root = comm->root[d];
5293 for (i = 0; i < dd->nc[dim]; i++)
5295 load->sum += load->load[pos++];
5296 load->max = max(load->max, load->load[pos]);
5302 /* This direction could not be load balanced properly,
5303 * therefore we need to use the maximum iso the average load.
5305 load->sum_m = max(load->sum_m, load->load[pos]);
5309 load->sum_m += load->load[pos];
5312 load->cvol_min = min(load->cvol_min, load->load[pos]);
5316 load->flags = (int)(load->load[pos++] + 0.5);
5320 root->cell_f_max0[i] = load->load[pos++];
5321 root->cell_f_min1[i] = load->load[pos++];
5326 load->mdf = max(load->mdf, load->load[pos]);
5328 load->pme = max(load->pme, load->load[pos]);
5332 if (comm->bDynLoadBal && root->bLimited)
5334 load->sum_m *= dd->nc[dim];
5335 load->flags |= (1<<d);
5343 comm->nload += dd_load_count(comm);
5344 comm->load_step += comm->cycl[ddCyclStep];
5345 comm->load_sum += comm->load[0].sum;
5346 comm->load_max += comm->load[0].max;
5347 if (comm->bDynLoadBal)
5349 for (d = 0; d < dd->ndim; d++)
5351 if (comm->load[0].flags & (1<<d))
5353 comm->load_lim[d]++;
5359 comm->load_mdf += comm->load[0].mdf;
5360 comm->load_pme += comm->load[0].pme;
5364 wallcycle_stop(wcycle, ewcDDCOMMLOAD);
5368 fprintf(debug, "get_load_distribution finished\n");
5372 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5374 /* Return the relative performance loss on the total run time
5375 * due to the force calculation load imbalance.
5377 if (dd->comm->nload > 0)
5380 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5381 (dd->comm->load_step*dd->nnodes);
5389 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
5392 int npp, npme, nnodes, d, limp;
5393 float imbal, pme_f_ratio, lossf, lossp = 0;
5395 gmx_domdec_comm_t *comm;
5398 if (DDMASTER(dd) && comm->nload > 0)
5401 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5402 nnodes = npp + npme;
5403 imbal = comm->load_max*npp/comm->load_sum - 1;
5404 lossf = dd_force_imb_perf_loss(dd);
5405 sprintf(buf, " Average load imbalance: %.1f %%\n", imbal*100);
5406 fprintf(fplog, "%s", buf);
5407 fprintf(stderr, "\n");
5408 fprintf(stderr, "%s", buf);
5409 sprintf(buf, " Part of the total run time spent waiting due to load imbalance: %.1f %%\n", lossf*100);
5410 fprintf(fplog, "%s", buf);
5411 fprintf(stderr, "%s", buf);
5413 if (comm->bDynLoadBal)
5415 sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5416 for (d = 0; d < dd->ndim; d++)
5418 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5419 sprintf(buf+strlen(buf), " %c %d %%", dim2char(dd->dim[d]), limp);
5425 sprintf(buf+strlen(buf), "\n");
5426 fprintf(fplog, "%s", buf);
5427 fprintf(stderr, "%s", buf);
5431 pme_f_ratio = comm->load_pme/comm->load_mdf;
5432 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5435 lossp *= (float)npme/(float)nnodes;
5439 lossp *= (float)npp/(float)nnodes;
5441 sprintf(buf, " Average PME mesh/force load: %5.3f\n", pme_f_ratio);
5442 fprintf(fplog, "%s", buf);
5443 fprintf(stderr, "%s", buf);
5444 sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", fabs(lossp)*100);
5445 fprintf(fplog, "%s", buf);
5446 fprintf(stderr, "%s", buf);
5448 fprintf(fplog, "\n");
5449 fprintf(stderr, "\n");
5451 if (lossf >= DD_PERF_LOSS)
5454 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5455 " in the domain decomposition.\n", lossf*100);
5456 if (!comm->bDynLoadBal)
5458 sprintf(buf+strlen(buf), " You might want to use dynamic load balancing (option -dlb.)\n");
5462 sprintf(buf+strlen(buf), " You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5464 fprintf(fplog, "%s\n", buf);
5465 fprintf(stderr, "%s\n", buf);
5467 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5470 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5471 " had %s work to do than the PP nodes.\n"
5472 " You might want to %s the number of PME nodes\n"
5473 " or %s the cut-off and the grid spacing.\n",
5475 (lossp < 0) ? "less" : "more",
5476 (lossp < 0) ? "decrease" : "increase",
5477 (lossp < 0) ? "decrease" : "increase");
5478 fprintf(fplog, "%s\n", buf);
5479 fprintf(stderr, "%s\n", buf);
5484 static float dd_vol_min(gmx_domdec_t *dd)
5486 return dd->comm->load[0].cvol_min*dd->nnodes;
5489 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5491 return dd->comm->load[0].flags;
5494 static float dd_f_imbal(gmx_domdec_t *dd)
5496 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5499 float dd_pme_f_ratio(gmx_domdec_t *dd)
5501 if (dd->comm->cycl_n[ddCyclPME] > 0)
5503 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5511 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, gmx_large_int_t step)
5516 flags = dd_load_flags(dd);
5520 "DD load balancing is limited by minimum cell size in dimension");
5521 for (d = 0; d < dd->ndim; d++)
5525 fprintf(fplog, " %c", dim2char(dd->dim[d]));
5528 fprintf(fplog, "\n");
5530 fprintf(fplog, "DD step %s", gmx_step_str(step, buf));
5531 if (dd->comm->bDynLoadBal)
5533 fprintf(fplog, " vol min/aver %5.3f%c",
5534 dd_vol_min(dd), flags ? '!' : ' ');
5536 fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
5537 if (dd->comm->cycl_n[ddCyclPME])
5539 fprintf(fplog, " pme mesh/force %5.3f", dd_pme_f_ratio(dd));
5541 fprintf(fplog, "\n\n");
5544 static void dd_print_load_verbose(gmx_domdec_t *dd)
5546 if (dd->comm->bDynLoadBal)
5548 fprintf(stderr, "vol %4.2f%c ",
5549 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
5551 fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
5552 if (dd->comm->cycl_n[ddCyclPME])
5554 fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
5559 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
5564 gmx_domdec_root_t *root;
5565 gmx_bool bPartOfGroup = FALSE;
5567 dim = dd->dim[dim_ind];
5568 copy_ivec(loc, loc_c);
5569 for (i = 0; i < dd->nc[dim]; i++)
5572 rank = dd_index(dd->nc, loc_c);
5573 if (rank == dd->rank)
5575 /* This process is part of the group */
5576 bPartOfGroup = TRUE;
5579 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
5583 dd->comm->mpi_comm_load[dim_ind] = c_row;
5584 if (dd->comm->eDLB != edlbNO)
5586 if (dd->ci[dim] == dd->master_ci[dim])
5588 /* This is the root process of this row */
5589 snew(dd->comm->root[dim_ind], 1);
5590 root = dd->comm->root[dim_ind];
5591 snew(root->cell_f, DD_CELL_F_SIZE(dd, dim_ind));
5592 snew(root->old_cell_f, dd->nc[dim]+1);
5593 snew(root->bCellMin, dd->nc[dim]);
5596 snew(root->cell_f_max0, dd->nc[dim]);
5597 snew(root->cell_f_min1, dd->nc[dim]);
5598 snew(root->bound_min, dd->nc[dim]);
5599 snew(root->bound_max, dd->nc[dim]);
5601 snew(root->buf_ncd, dd->nc[dim]);
5605 /* This is not a root process, we only need to receive cell_f */
5606 snew(dd->comm->cell_f_row, DD_CELL_F_SIZE(dd, dim_ind));
5609 if (dd->ci[dim] == dd->master_ci[dim])
5611 snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
5617 static void make_load_communicators(gmx_domdec_t *dd)
5620 int dim0, dim1, i, j;
5625 fprintf(debug, "Making load communicators\n");
5628 snew(dd->comm->load, dd->ndim);
5629 snew(dd->comm->mpi_comm_load, dd->ndim);
5632 make_load_communicator(dd, 0, loc);
5636 for (i = 0; i < dd->nc[dim0]; i++)
5639 make_load_communicator(dd, 1, loc);
5645 for (i = 0; i < dd->nc[dim0]; i++)
5649 for (j = 0; j < dd->nc[dim1]; j++)
5652 make_load_communicator(dd, 2, loc);
5659 fprintf(debug, "Finished making load communicators\n");
5664 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
5667 int d, dim, i, j, m;
5670 ivec dd_zp[DD_MAXIZONE];
5671 gmx_domdec_zones_t *zones;
5672 gmx_domdec_ns_ranges_t *izone;
5674 for (d = 0; d < dd->ndim; d++)
5677 copy_ivec(dd->ci, tmp);
5678 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5679 dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
5680 copy_ivec(dd->ci, tmp);
5681 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5682 dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
5685 fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5688 dd->neighbor[d][1]);
5694 fprintf(fplog, "\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5696 dd->nc[XX], dd->nc[YY], dd->nc[ZZ],
5697 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5704 for (i = 0; i < nzonep; i++)
5706 copy_ivec(dd_zp3[i], dd_zp[i]);
5712 for (i = 0; i < nzonep; i++)
5714 copy_ivec(dd_zp2[i], dd_zp[i]);
5720 for (i = 0; i < nzonep; i++)
5722 copy_ivec(dd_zp1[i], dd_zp[i]);
5726 gmx_fatal(FARGS, "Can only do 1, 2 or 3D domain decomposition");
5731 zones = &dd->comm->zones;
5733 for (i = 0; i < nzone; i++)
5736 clear_ivec(zones->shift[i]);
5737 for (d = 0; d < dd->ndim; d++)
5739 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5744 for (i = 0; i < nzone; i++)
5746 for (d = 0; d < DIM; d++)
5748 s[d] = dd->ci[d] - zones->shift[i][d];
5753 else if (s[d] >= dd->nc[d])
5759 zones->nizone = nzonep;
5760 for (i = 0; i < zones->nizone; i++)
5762 if (dd_zp[i][0] != i)
5764 gmx_fatal(FARGS, "Internal inconsistency in the dd grid setup");
5766 izone = &zones->izone[i];
5767 izone->j0 = dd_zp[i][1];
5768 izone->j1 = dd_zp[i][2];
5769 for (dim = 0; dim < DIM; dim++)
5771 if (dd->nc[dim] == 1)
5773 /* All shifts should be allowed */
5774 izone->shift0[dim] = -1;
5775 izone->shift1[dim] = 1;
5780 izone->shift0[d] = 0;
5781 izone->shift1[d] = 0;
5782 for(j=izone->j0; j<izone->j1; j++) {
5783 if (dd->shift[j][d] > dd->shift[i][d])
5784 izone->shift0[d] = -1;
5785 if (dd->shift[j][d] < dd->shift[i][d])
5786 izone->shift1[d] = 1;
5792 /* Assume the shift are not more than 1 cell */
5793 izone->shift0[dim] = 1;
5794 izone->shift1[dim] = -1;
5795 for (j = izone->j0; j < izone->j1; j++)
5797 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5798 if (shift_diff < izone->shift0[dim])
5800 izone->shift0[dim] = shift_diff;
5802 if (shift_diff > izone->shift1[dim])
5804 izone->shift1[dim] = shift_diff;
5811 if (dd->comm->eDLB != edlbNO)
5813 snew(dd->comm->root, dd->ndim);
5816 if (dd->comm->bRecordLoad)
5818 make_load_communicators(dd);
5822 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int reorder)
5825 gmx_domdec_comm_t *comm;
5836 if (comm->bCartesianPP)
5838 /* Set up cartesian communication for the particle-particle part */
5841 fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
5842 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
5845 for (i = 0; i < DIM; i++)
5849 MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
5851 /* We overwrite the old communicator with the new cartesian one */
5852 cr->mpi_comm_mygroup = comm_cart;
5855 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5856 MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
5858 if (comm->bCartesianPP_PME)
5860 /* Since we want to use the original cartesian setup for sim,
5861 * and not the one after split, we need to make an index.
5863 snew(comm->ddindex2ddnodeid, dd->nnodes);
5864 comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
5865 gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
5866 /* Get the rank of the DD master,
5867 * above we made sure that the master node is a PP node.
5877 MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
5879 else if (comm->bCartesianPP)
5881 if (cr->npmenodes == 0)
5883 /* The PP communicator is also
5884 * the communicator for this simulation
5886 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5888 cr->nodeid = dd->rank;
5890 MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
5892 /* We need to make an index to go from the coordinates
5893 * to the nodeid of this simulation.
5895 snew(comm->ddindex2simnodeid, dd->nnodes);
5896 snew(buf, dd->nnodes);
5897 if (cr->duty & DUTY_PP)
5899 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5901 /* Communicate the ddindex to simulation nodeid index */
5902 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5903 cr->mpi_comm_mysim);
5906 /* Determine the master coordinates and rank.
5907 * The DD master should be the same node as the master of this sim.
5909 for (i = 0; i < dd->nnodes; i++)
5911 if (comm->ddindex2simnodeid[i] == 0)
5913 ddindex2xyz(dd->nc, i, dd->master_ci);
5914 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
5919 fprintf(debug, "The master rank is %d\n", dd->masterrank);
5924 /* No Cartesian communicators */
5925 /* We use the rank in dd->comm->all as DD index */
5926 ddindex2xyz(dd->nc, dd->rank, dd->ci);
5927 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5929 clear_ivec(dd->master_ci);
5936 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5937 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5942 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5943 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
5947 static void receive_ddindex2simnodeid(t_commrec *cr)
5951 gmx_domdec_comm_t *comm;
5958 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5960 snew(comm->ddindex2simnodeid, dd->nnodes);
5961 snew(buf, dd->nnodes);
5962 if (cr->duty & DUTY_PP)
5964 buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
5967 /* Communicate the ddindex to simulation nodeid index */
5968 MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
5969 cr->mpi_comm_mysim);
5976 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5977 int ncg, int natoms)
5979 gmx_domdec_master_t *ma;
5984 snew(ma->ncg, dd->nnodes);
5985 snew(ma->index, dd->nnodes+1);
5987 snew(ma->nat, dd->nnodes);
5988 snew(ma->ibuf, dd->nnodes*2);
5989 snew(ma->cell_x, DIM);
5990 for (i = 0; i < DIM; i++)
5992 snew(ma->cell_x[i], dd->nc[i]+1);
5995 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
6001 snew(ma->vbuf, natoms);
6007 static void split_communicator(FILE *fplog, t_commrec *cr, int dd_node_order,
6011 gmx_domdec_comm_t *comm;
6022 if (comm->bCartesianPP)
6024 for (i = 1; i < DIM; i++)
6026 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
6028 if (bDiv[YY] || bDiv[ZZ])
6030 comm->bCartesianPP_PME = TRUE;
6031 /* If we have 2D PME decomposition, which is always in x+y,
6032 * we stack the PME only nodes in z.
6033 * Otherwise we choose the direction that provides the thinnest slab
6034 * of PME only nodes as this will have the least effect
6035 * on the PP communication.
6036 * But for the PME communication the opposite might be better.
6038 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
6040 dd->nc[YY] > dd->nc[ZZ]))
6042 comm->cartpmedim = ZZ;
6046 comm->cartpmedim = YY;
6048 comm->ntot[comm->cartpmedim]
6049 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
6053 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]);
6055 "Will not use a Cartesian communicator for PP <-> PME\n\n");
6060 if (comm->bCartesianPP_PME)
6064 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]);
6067 for (i = 0; i < DIM; i++)
6071 MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
6074 MPI_Comm_rank(comm_cart, &rank);
6075 if (MASTERNODE(cr) && rank != 0)
6077 gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
6080 /* With this assigment we loose the link to the original communicator
6081 * which will usually be MPI_COMM_WORLD, unless have multisim.
6083 cr->mpi_comm_mysim = comm_cart;
6084 cr->sim_nodeid = rank;
6086 MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
6090 fprintf(fplog, "Cartesian nodeid %d, coordinates %d %d %d\n\n",
6091 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
6094 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
6098 if (cr->npmenodes == 0 ||
6099 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6101 cr->duty = DUTY_PME;
6104 /* Split the sim communicator into PP and PME only nodes */
6105 MPI_Comm_split(cr->mpi_comm_mysim,
6107 dd_index(comm->ntot, dd->ci),
6108 &cr->mpi_comm_mygroup);
6112 switch (dd_node_order)
6117 fprintf(fplog, "Order of the nodes: PP first, PME last\n");
6120 case ddnoINTERLEAVE:
6121 /* Interleave the PP-only and PME-only nodes,
6122 * as on clusters with dual-core machines this will double
6123 * the communication bandwidth of the PME processes
6124 * and thus speed up the PP <-> PME and inter PME communication.
6128 fprintf(fplog, "Interleaving PP and PME nodes\n");
6130 comm->pmenodes = dd_pmenodes(cr);
6135 gmx_fatal(FARGS, "Unknown dd_node_order=%d", dd_node_order);
6138 if (dd_simnode2pmenode(cr, cr->sim_nodeid) == -1)
6140 cr->duty = DUTY_PME;
6147 /* Split the sim communicator into PP and PME only nodes */
6148 MPI_Comm_split(cr->mpi_comm_mysim,
6151 &cr->mpi_comm_mygroup);
6152 MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
6158 fprintf(fplog, "This is a %s only node\n\n",
6159 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6163 void make_dd_communicators(FILE *fplog, t_commrec *cr, int dd_node_order)
6166 gmx_domdec_comm_t *comm;
6172 copy_ivec(dd->nc, comm->ntot);
6174 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6175 comm->bCartesianPP_PME = FALSE;
6177 /* Reorder the nodes by default. This might change the MPI ranks.
6178 * Real reordering is only supported on very few architectures,
6179 * Blue Gene is one of them.
6181 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6183 if (cr->npmenodes > 0)
6185 /* Split the communicator into a PP and PME part */
6186 split_communicator(fplog, cr, dd_node_order, CartReorder);
6187 if (comm->bCartesianPP_PME)
6189 /* We (possibly) reordered the nodes in split_communicator,
6190 * so it is no longer required in make_pp_communicator.
6192 CartReorder = FALSE;
6197 /* All nodes do PP and PME */
6199 /* We do not require separate communicators */
6200 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6204 if (cr->duty & DUTY_PP)
6206 /* Copy or make a new PP communicator */
6207 make_pp_communicator(fplog, cr, CartReorder);
6211 receive_ddindex2simnodeid(cr);
6214 if (!(cr->duty & DUTY_PME))
6216 /* Set up the commnuication to our PME node */
6217 dd->pme_nodeid = dd_simnode2pmenode(cr, cr->sim_nodeid);
6218 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6221 fprintf(debug, "My pme_nodeid %d receive ener %d\n",
6222 dd->pme_nodeid, dd->pme_receive_vir_ener);
6227 dd->pme_nodeid = -1;
6232 dd->ma = init_gmx_domdec_master_t(dd,
6234 comm->cgs_gl.index[comm->cgs_gl.nr]);
6238 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
6240 real *slb_frac, tot;
6245 if (nc > 1 && size_string != NULL)
6249 fprintf(fplog, "Using static load balancing for the %s direction\n",
6254 for (i = 0; i < nc; i++)
6257 sscanf(size_string, "%lf%n", &dbl, &n);
6260 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
6269 fprintf(fplog, "Relative cell sizes:");
6271 for (i = 0; i < nc; i++)
6276 fprintf(fplog, " %5.3f", slb_frac[i]);
6281 fprintf(fplog, "\n");
6288 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6291 gmx_mtop_ilistloop_t iloop;
6295 iloop = gmx_mtop_ilistloop_init(mtop);
6296 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
6298 for (ftype = 0; ftype < F_NRE; ftype++)
6300 if ((interaction_function[ftype].flags & IF_BOND) &&
6303 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6311 static int dd_nst_env(FILE *fplog, const char *env_var, int def)
6317 val = getenv(env_var);
6320 if (sscanf(val, "%d", &nst) <= 0)
6326 fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
6334 static void dd_warning(t_commrec *cr, FILE *fplog, const char *warn_string)
6338 fprintf(stderr, "\n%s\n", warn_string);
6342 fprintf(fplog, "\n%s\n", warn_string);
6346 static void check_dd_restrictions(t_commrec *cr, gmx_domdec_t *dd,
6347 t_inputrec *ir, FILE *fplog)
6349 if (ir->ePBC == epbcSCREW &&
6350 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6352 gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
6355 if (ir->ns_type == ensSIMPLE)
6357 gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6360 if (ir->nstlist == 0)
6362 gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
6365 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6367 dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6371 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
6376 r = ddbox->box_size[XX];
6377 for (di = 0; di < dd->ndim; di++)
6380 /* Check using the initial average cell size */
6381 r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6387 static int check_dlb_support(FILE *fplog, t_commrec *cr,
6388 const char *dlb_opt, gmx_bool bRecordLoad,
6389 unsigned long Flags, t_inputrec *ir)
6397 case 'a': eDLB = edlbAUTO; break;
6398 case 'n': eDLB = edlbNO; break;
6399 case 'y': eDLB = edlbYES; break;
6400 default: gmx_incons("Unknown dlb_opt");
6403 if (Flags & MD_RERUN)
6408 if (!EI_DYNAMICS(ir->eI))
6410 if (eDLB == edlbYES)
6412 sprintf(buf, "NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n", EI(ir->eI));
6413 dd_warning(cr, fplog, buf);
6421 dd_warning(cr, fplog, "NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6426 if (Flags & MD_REPRODUCIBLE)
6433 dd_warning(cr, fplog, "NOTE: reproducibility requested, will not use dynamic load balancing\n");
6437 dd_warning(cr, fplog, "WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6440 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", eDLB);
6448 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
6453 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6455 /* Decomposition order z,y,x */
6458 fprintf(fplog, "Using domain decomposition order z, y, x\n");
6460 for (dim = DIM-1; dim >= 0; dim--)
6462 if (dd->nc[dim] > 1)
6464 dd->dim[dd->ndim++] = dim;
6470 /* Decomposition order x,y,z */
6471 for (dim = 0; dim < DIM; dim++)
6473 if (dd->nc[dim] > 1)
6475 dd->dim[dd->ndim++] = dim;
6481 static gmx_domdec_comm_t *init_dd_comm()
6483 gmx_domdec_comm_t *comm;
6487 snew(comm->cggl_flag, DIM*2);
6488 snew(comm->cgcm_state, DIM*2);
6489 for (i = 0; i < DIM*2; i++)
6491 comm->cggl_flag_nalloc[i] = 0;
6492 comm->cgcm_state_nalloc[i] = 0;
6495 comm->nalloc_int = 0;
6496 comm->buf_int = NULL;
6498 vec_rvec_init(&comm->vbuf);
6500 comm->n_load_have = 0;
6501 comm->n_load_collect = 0;
6503 for (i = 0; i < ddnatNR-ddnatZONE; i++)
6505 comm->sum_nat[i] = 0;
6509 comm->load_step = 0;
6512 clear_ivec(comm->load_lim);
6519 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
6520 unsigned long Flags,
6522 real comm_distance_min, real rconstr,
6523 const char *dlb_opt, real dlb_scale,
6524 const char *sizex, const char *sizey, const char *sizez,
6525 gmx_mtop_t *mtop, t_inputrec *ir,
6526 matrix box, rvec *x,
6528 int *npme_x, int *npme_y)
6531 gmx_domdec_comm_t *comm;
6534 real r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
6541 "\nInitializing Domain Decomposition on %d nodes\n", cr->nnodes);
6546 dd->comm = init_dd_comm();
6548 snew(comm->cggl_flag, DIM*2);
6549 snew(comm->cgcm_state, DIM*2);
6551 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6552 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6554 dd->bSendRecv2 = dd_nst_env(fplog, "GMX_DD_SENDRECV2", 0);
6555 comm->dlb_scale_lim = dd_nst_env(fplog, "GMX_DLB_MAX", 10);
6556 comm->eFlop = dd_nst_env(fplog, "GMX_DLB_FLOP", 0);
6557 recload = dd_nst_env(fplog, "GMX_DD_LOAD", 1);
6558 comm->nstSortCG = dd_nst_env(fplog, "GMX_DD_SORT", 1);
6559 comm->nstDDDump = dd_nst_env(fplog, "GMX_DD_DUMP", 0);
6560 comm->nstDDDumpGrid = dd_nst_env(fplog, "GMX_DD_DUMP_GRID", 0);
6561 comm->DD_debug = dd_nst_env(fplog, "GMX_DD_DEBUG", 0);
6563 dd->pme_recv_f_alloc = 0;
6564 dd->pme_recv_f_buf = NULL;
6566 if (dd->bSendRecv2 && fplog)
6568 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");
6574 fprintf(fplog, "Will load balance based on FLOP count\n");
6576 if (comm->eFlop > 1)
6578 srand(1+cr->nodeid);
6580 comm->bRecordLoad = TRUE;
6584 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6588 comm->eDLB = check_dlb_support(fplog, cr, dlb_opt, comm->bRecordLoad, Flags, ir);
6590 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6593 fprintf(fplog, "Dynamic load balancing: %s\n", edlb_names[comm->eDLB]);
6595 dd->bGridJump = comm->bDynLoadBal;
6596 comm->bPMELoadBalDLBLimits = FALSE;
6598 if (comm->nstSortCG)
6602 if (comm->nstSortCG == 1)
6604 fprintf(fplog, "Will sort the charge groups at every domain (re)decomposition\n");
6608 fprintf(fplog, "Will sort the charge groups every %d steps\n",
6612 snew(comm->sort, 1);
6618 fprintf(fplog, "Will not sort the charge groups\n");
6622 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6624 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6625 if (comm->bInterCGBondeds)
6627 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6631 comm->bInterCGMultiBody = FALSE;
6634 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6635 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6637 if (ir->rlistlong == 0)
6639 /* Set the cut-off to some very large value,
6640 * so we don't need if statements everywhere in the code.
6641 * We use sqrt, since the cut-off is squared in some places.
6643 comm->cutoff = GMX_CUTOFF_INF;
6647 comm->cutoff = ir->rlistlong;
6649 comm->cutoff_mbody = 0;
6651 comm->cellsize_limit = 0;
6652 comm->bBondComm = FALSE;
6654 if (comm->bInterCGBondeds)
6656 if (comm_distance_min > 0)
6658 comm->cutoff_mbody = comm_distance_min;
6659 if (Flags & MD_DDBONDCOMM)
6661 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6665 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6667 r_bonded_limit = comm->cutoff_mbody;
6669 else if (ir->bPeriodicMols)
6671 /* Can not easily determine the required cut-off */
6672 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");
6673 comm->cutoff_mbody = comm->cutoff/2;
6674 r_bonded_limit = comm->cutoff_mbody;
6680 dd_bonded_cg_distance(fplog, dd, mtop, ir, x, box,
6681 Flags & MD_DDBONDCHECK, &r_2b, &r_mb);
6683 gmx_bcast(sizeof(r_2b), &r_2b, cr);
6684 gmx_bcast(sizeof(r_mb), &r_mb, cr);
6686 /* We use an initial margin of 10% for the minimum cell size,
6687 * except when we are just below the non-bonded cut-off.
6689 if (Flags & MD_DDBONDCOMM)
6691 if (max(r_2b, r_mb) > comm->cutoff)
6693 r_bonded = max(r_2b, r_mb);
6694 r_bonded_limit = 1.1*r_bonded;
6695 comm->bBondComm = TRUE;
6700 r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
6702 /* We determine cutoff_mbody later */
6706 /* No special bonded communication,
6707 * simply increase the DD cut-off.
6709 r_bonded_limit = 1.1*max(r_2b, r_mb);
6710 comm->cutoff_mbody = r_bonded_limit;
6711 comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
6714 comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
6718 "Minimum cell size due to bonded interactions: %.3f nm\n",
6719 comm->cellsize_limit);
6723 if (dd->bInterCGcons && rconstr <= 0)
6725 /* There is a cell size limit due to the constraints (P-LINCS) */
6726 rconstr = constr_r_max(fplog, mtop, ir);
6730 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6732 if (rconstr > comm->cellsize_limit)
6734 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
6738 else if (rconstr > 0 && fplog)
6740 /* Here we do not check for dd->bInterCGcons,
6741 * because one can also set a cell size limit for virtual sites only
6742 * and at this point we don't know yet if there are intercg v-sites.
6745 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6748 comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
6750 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6754 copy_ivec(nc, dd->nc);
6755 set_dd_dim(fplog, dd);
6756 set_ddbox_cr(cr, &dd->nc, ir, box, &comm->cgs_gl, x, ddbox);
6758 if (cr->npmenodes == -1)
6762 acs = average_cellsize_min(dd, ddbox);
6763 if (acs < comm->cellsize_limit)
6767 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
6769 gmx_fatal_collective(FARGS, cr, NULL,
6770 "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",
6771 acs, comm->cellsize_limit);
6776 set_ddbox_cr(cr, NULL, ir, box, &comm->cgs_gl, x, ddbox);
6778 /* We need to choose the optimal DD grid and possibly PME nodes */
6779 limit = dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
6780 comm->eDLB != edlbNO, dlb_scale,
6781 comm->cellsize_limit, comm->cutoff,
6782 comm->bInterCGBondeds, comm->bInterCGMultiBody);
6784 if (dd->nc[XX] == 0)
6786 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6787 sprintf(buf, "Change the number of nodes or mdrun option %s%s%s",
6788 !bC ? "-rdd" : "-rcon",
6789 comm->eDLB != edlbNO ? " or -dds" : "",
6790 bC ? " or your LINCS settings" : "");
6792 gmx_fatal_collective(FARGS, cr, NULL,
6793 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6795 "Look in the log file for details on the domain decomposition",
6796 cr->nnodes-cr->npmenodes, limit, buf);
6798 set_dd_dim(fplog, dd);
6804 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6805 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
6808 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6809 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6811 gmx_fatal_collective(FARGS, cr, NULL,
6812 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6813 dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
6815 if (cr->npmenodes > dd->nnodes)
6817 gmx_fatal_collective(FARGS, cr, NULL,
6818 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.", cr->npmenodes, dd->nnodes);
6820 if (cr->npmenodes > 0)
6822 comm->npmenodes = cr->npmenodes;
6826 comm->npmenodes = dd->nnodes;
6829 if (EEL_PME(ir->coulombtype))
6831 /* The following choices should match those
6832 * in comm_cost_est in domdec_setup.c.
6833 * Note that here the checks have to take into account
6834 * that the decomposition might occur in a different order than xyz
6835 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6836 * in which case they will not match those in comm_cost_est,
6837 * but since that is mainly for testing purposes that's fine.
6839 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6840 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6841 getenv("GMX_PMEONEDD") == NULL)
6843 comm->npmedecompdim = 2;
6844 comm->npmenodes_x = dd->nc[XX];
6845 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6849 /* In case nc is 1 in both x and y we could still choose to
6850 * decompose pme in y instead of x, but we use x for simplicity.
6852 comm->npmedecompdim = 1;
6853 if (dd->dim[0] == YY)
6855 comm->npmenodes_x = 1;
6856 comm->npmenodes_y = comm->npmenodes;
6860 comm->npmenodes_x = comm->npmenodes;
6861 comm->npmenodes_y = 1;
6866 fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
6867 comm->npmenodes_x, comm->npmenodes_y, 1);
6872 comm->npmedecompdim = 0;
6873 comm->npmenodes_x = 0;
6874 comm->npmenodes_y = 0;
6877 /* Technically we don't need both of these,
6878 * but it simplifies code not having to recalculate it.
6880 *npme_x = comm->npmenodes_x;
6881 *npme_y = comm->npmenodes_y;
6883 snew(comm->slb_frac, DIM);
6884 if (comm->eDLB == edlbNO)
6886 comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], sizex);
6887 comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], sizey);
6888 comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], sizez);
6891 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6893 if (comm->bBondComm || comm->eDLB != edlbNO)
6895 /* Set the bonded communication distance to halfway
6896 * the minimum and the maximum,
6897 * since the extra communication cost is nearly zero.
6899 acs = average_cellsize_min(dd, ddbox);
6900 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6901 if (comm->eDLB != edlbNO)
6903 /* Check if this does not limit the scaling */
6904 comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
6906 if (!comm->bBondComm)
6908 /* Without bBondComm do not go beyond the n.b. cut-off */
6909 comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
6910 if (comm->cellsize_limit >= comm->cutoff)
6912 /* We don't loose a lot of efficieny
6913 * when increasing it to the n.b. cut-off.
6914 * It can even be slightly faster, because we need
6915 * less checks for the communication setup.
6917 comm->cutoff_mbody = comm->cutoff;
6920 /* Check if we did not end up below our original limit */
6921 comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
6923 if (comm->cutoff_mbody > comm->cellsize_limit)
6925 comm->cellsize_limit = comm->cutoff_mbody;
6928 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6933 fprintf(debug, "Bonded atom communication beyond the cut-off: %d\n"
6934 "cellsize limit %f\n",
6935 comm->bBondComm, comm->cellsize_limit);
6940 check_dd_restrictions(cr, dd, ir, fplog);
6943 comm->partition_step = INT_MIN;
6946 clear_dd_cycle_counts(dd);
6951 static void set_dlb_limits(gmx_domdec_t *dd)
6956 for (d = 0; d < dd->ndim; d++)
6958 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6959 dd->comm->cellsize_min[dd->dim[d]] =
6960 dd->comm->cellsize_min_dlb[dd->dim[d]];
6965 static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_large_int_t step)
6968 gmx_domdec_comm_t *comm;
6978 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);
6981 cellsize_min = comm->cellsize_min[dd->dim[0]];
6982 for (d = 1; d < dd->ndim; d++)
6984 cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
6987 if (cellsize_min < comm->cellsize_limit*1.05)
6989 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");
6991 /* Change DLB from "auto" to "no". */
6992 comm->eDLB = edlbNO;
6997 dd_warning(cr, fplog, "NOTE: Turning on dynamic load balancing\n");
6998 comm->bDynLoadBal = TRUE;
6999 dd->bGridJump = TRUE;
7003 /* We can set the required cell size info here,
7004 * so we do not need to communicate this.
7005 * The grid is completely uniform.
7007 for (d = 0; d < dd->ndim; d++)
7011 comm->load[d].sum_m = comm->load[d].sum;
7013 nc = dd->nc[dd->dim[d]];
7014 for (i = 0; i < nc; i++)
7016 comm->root[d]->cell_f[i] = i/(real)nc;
7019 comm->root[d]->cell_f_max0[i] = i /(real)nc;
7020 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
7023 comm->root[d]->cell_f[nc] = 1.0;
7028 static char *init_bLocalCG(gmx_mtop_t *mtop)
7033 ncg = ncg_mtop(mtop);
7034 snew(bLocalCG, ncg);
7035 for (cg = 0; cg < ncg; cg++)
7037 bLocalCG[cg] = FALSE;
7043 void dd_init_bondeds(FILE *fplog,
7044 gmx_domdec_t *dd, gmx_mtop_t *mtop,
7045 gmx_vsite_t *vsite, gmx_constr_t constr,
7046 t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
7048 gmx_domdec_comm_t *comm;
7052 dd_make_reverse_top(fplog, dd, mtop, vsite, constr, ir, bBCheck);
7056 if (comm->bBondComm)
7058 /* Communicate atoms beyond the cut-off for bonded interactions */
7061 comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
7063 comm->bLocalCG = init_bLocalCG(mtop);
7067 /* Only communicate atoms based on cut-off */
7068 comm->cglink = NULL;
7069 comm->bLocalCG = NULL;
7073 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
7075 gmx_bool bDynLoadBal, real dlb_scale,
7078 gmx_domdec_comm_t *comm;
7093 fprintf(fplog, "The maximum number of communication pulses is:");
7094 for (d = 0; d < dd->ndim; d++)
7096 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
7098 fprintf(fplog, "\n");
7099 fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
7100 fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
7101 fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
7102 for (d = 0; d < DIM; d++)
7106 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7113 comm->cellsize_min_dlb[d]/
7114 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7116 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
7119 fprintf(fplog, "\n");
7123 set_dd_cell_sizes_slb(dd, ddbox, FALSE, np);
7124 fprintf(fplog, "The initial number of communication pulses is:");
7125 for (d = 0; d < dd->ndim; d++)
7127 fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
7129 fprintf(fplog, "\n");
7130 fprintf(fplog, "The initial domain decomposition cell size is:");
7131 for (d = 0; d < DIM; d++)
7135 fprintf(fplog, " %c %.2f nm",
7136 dim2char(d), dd->comm->cellsize_min[d]);
7139 fprintf(fplog, "\n\n");
7142 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7144 fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
7145 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7146 "non-bonded interactions", "", comm->cutoff);
7150 limit = dd->comm->cellsize_limit;
7154 if (dynamic_dd_box(ddbox, ir))
7156 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
7158 limit = dd->comm->cellsize_min[XX];
7159 for (d = 1; d < DIM; d++)
7161 limit = min(limit, dd->comm->cellsize_min[d]);
7165 if (comm->bInterCGBondeds)
7167 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7168 "two-body bonded interactions", "(-rdd)",
7169 max(comm->cutoff, comm->cutoff_mbody));
7170 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7171 "multi-body bonded interactions", "(-rdd)",
7172 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
7176 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7177 "virtual site constructions", "(-rcon)", limit);
7179 if (dd->constraint_comm)
7181 sprintf(buf, "atoms separated by up to %d constraints",
7183 fprintf(fplog, "%40s %-7s %6.3f nm\n",
7184 buf, "(-rcon)", limit);
7186 fprintf(fplog, "\n");
7192 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7194 const t_inputrec *ir,
7195 const gmx_ddbox_t *ddbox)
7197 gmx_domdec_comm_t *comm;
7198 int d, dim, npulse, npulse_d_max, npulse_d;
7203 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7205 /* Determine the maximum number of comm. pulses in one dimension */
7207 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7209 /* Determine the maximum required number of grid pulses */
7210 if (comm->cellsize_limit >= comm->cutoff)
7212 /* Only a single pulse is required */
7215 else if (!bNoCutOff && comm->cellsize_limit > 0)
7217 /* We round down slightly here to avoid overhead due to the latency
7218 * of extra communication calls when the cut-off
7219 * would be only slightly longer than the cell size.
7220 * Later cellsize_limit is redetermined,
7221 * so we can not miss interactions due to this rounding.
7223 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7227 /* There is no cell size limit */
7228 npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
7231 if (!bNoCutOff && npulse > 1)
7233 /* See if we can do with less pulses, based on dlb_scale */
7235 for (d = 0; d < dd->ndim; d++)
7238 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7239 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7240 npulse_d_max = max(npulse_d_max, npulse_d);
7242 npulse = min(npulse, npulse_d_max);
7245 /* This env var can override npulse */
7246 d = dd_nst_env(debug, "GMX_DD_NPULSE", 0);
7253 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7254 for (d = 0; d < dd->ndim; d++)
7256 comm->cd[d].np_dlb = min(npulse, dd->nc[dd->dim[d]]-1);
7257 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7258 snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
7259 comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
7260 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7262 comm->bVacDLBNoLimit = FALSE;
7266 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7267 if (!comm->bVacDLBNoLimit)
7269 comm->cellsize_limit = max(comm->cellsize_limit,
7270 comm->cutoff/comm->maxpulse);
7272 comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
7273 /* Set the minimum cell size for each DD dimension */
7274 for (d = 0; d < dd->ndim; d++)
7276 if (comm->bVacDLBNoLimit ||
7277 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7279 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7283 comm->cellsize_min_dlb[dd->dim[d]] =
7284 comm->cutoff/comm->cd[d].np_dlb;
7287 if (comm->cutoff_mbody <= 0)
7289 comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
7291 if (comm->bDynLoadBal)
7297 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd, int ePBC)
7299 /* If each molecule is a single charge group
7300 * or we use domain decomposition for each periodic dimension,
7301 * we do not need to take pbc into account for the bonded interactions.
7303 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7306 (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
7309 void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
7310 t_inputrec *ir, t_forcerec *fr,
7313 gmx_domdec_comm_t *comm;
7319 /* Initialize the thread data.
7320 * This can not be done in init_domain_decomposition,
7321 * as the numbers of threads is determined later.
7323 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7326 snew(comm->dth, comm->nth);
7329 if (EEL_PME(ir->coulombtype))
7331 init_ddpme(dd, &comm->ddpme[0], 0);
7332 if (comm->npmedecompdim >= 2)
7334 init_ddpme(dd, &comm->ddpme[1], 1);
7339 comm->npmenodes = 0;
7340 if (dd->pme_nodeid >= 0)
7342 gmx_fatal_collective(FARGS, NULL, dd,
7343 "Can not have separate PME nodes without PME electrostatics");
7349 fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
7351 if (comm->eDLB != edlbNO)
7353 set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
7356 print_dd_settings(fplog, dd, ir, comm->bDynLoadBal, dlb_scale, ddbox);
7357 if (comm->eDLB == edlbAUTO)
7361 fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
7363 print_dd_settings(fplog, dd, ir, TRUE, dlb_scale, ddbox);
7366 if (ir->ePBC == epbcNONE)
7368 vol_frac = 1 - 1/(double)dd->nnodes;
7373 (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
7377 fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
7379 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7381 dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
7384 static gmx_bool test_dd_cutoff(t_commrec *cr,
7385 t_state *state, t_inputrec *ir,
7396 set_ddbox(dd, FALSE, cr, ir, state->box,
7397 TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
7401 for (d = 0; d < dd->ndim; d++)
7405 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7406 if (dynamic_dd_box(&ddbox, ir))
7408 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7411 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7413 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7414 dd->comm->cd[d].np_dlb > 0)
7416 if (np > dd->comm->cd[d].np_dlb)
7421 /* If a current local cell size is smaller than the requested
7422 * cut-off, we could still fix it, but this gets very complicated.
7423 * Without fixing here, we might actually need more checks.
7425 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7432 if (dd->comm->eDLB != edlbNO)
7434 /* If DLB is not active yet, we don't need to check the grid jumps.
7435 * Actually we shouldn't, because then the grid jump data is not set.
7437 if (dd->comm->bDynLoadBal &&
7438 check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
7443 gmx_sumi(1, &LocallyLimited, cr);
7445 if (LocallyLimited > 0)
7454 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, t_inputrec *ir,
7457 gmx_bool bCutoffAllowed;
7459 bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
7463 cr->dd->comm->cutoff = cutoff_req;
7466 return bCutoffAllowed;
7469 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7471 gmx_domdec_comm_t *comm;
7473 comm = cr->dd->comm;
7475 /* Turn on the DLB limiting (might have been on already) */
7476 comm->bPMELoadBalDLBLimits = TRUE;
7478 /* Change the cut-off limit */
7479 comm->PMELoadBal_max_cutoff = comm->cutoff;
7482 static void merge_cg_buffers(int ncell,
7483 gmx_domdec_comm_dim_t *cd, int pulse,
7485 int *index_gl, int *recv_i,
7486 rvec *cg_cm, rvec *recv_vr,
7488 cginfo_mb_t *cginfo_mb, int *cginfo)
7490 gmx_domdec_ind_t *ind, *ind_p;
7491 int p, cell, c, cg, cg0, cg1, cg_gl, nat;
7492 int shift, shift_at;
7494 ind = &cd->ind[pulse];
7496 /* First correct the already stored data */
7497 shift = ind->nrecv[ncell];
7498 for (cell = ncell-1; cell >= 0; cell--)
7500 shift -= ind->nrecv[cell];
7503 /* Move the cg's present from previous grid pulses */
7504 cg0 = ncg_cell[ncell+cell];
7505 cg1 = ncg_cell[ncell+cell+1];
7506 cgindex[cg1+shift] = cgindex[cg1];
7507 for (cg = cg1-1; cg >= cg0; cg--)
7509 index_gl[cg+shift] = index_gl[cg];
7510 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
7511 cgindex[cg+shift] = cgindex[cg];
7512 cginfo[cg+shift] = cginfo[cg];
7514 /* Correct the already stored send indices for the shift */
7515 for (p = 1; p <= pulse; p++)
7517 ind_p = &cd->ind[p];
7519 for (c = 0; c < cell; c++)
7521 cg0 += ind_p->nsend[c];
7523 cg1 = cg0 + ind_p->nsend[cell];
7524 for (cg = cg0; cg < cg1; cg++)
7526 ind_p->index[cg] += shift;
7532 /* Merge in the communicated buffers */
7536 for (cell = 0; cell < ncell; cell++)
7538 cg1 = ncg_cell[ncell+cell+1] + shift;
7541 /* Correct the old cg indices */
7542 for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
7544 cgindex[cg+1] += shift_at;
7547 for (cg = 0; cg < ind->nrecv[cell]; cg++)
7549 /* Copy this charge group from the buffer */
7550 index_gl[cg1] = recv_i[cg0];
7551 copy_rvec(recv_vr[cg0], cg_cm[cg1]);
7552 /* Add it to the cgindex */
7553 cg_gl = index_gl[cg1];
7554 cginfo[cg1] = ddcginfo(cginfo_mb, cg_gl);
7555 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7556 cgindex[cg1+1] = cgindex[cg1] + nat;
7561 shift += ind->nrecv[cell];
7562 ncg_cell[ncell+cell+1] = cg1;
7566 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7567 int nzone, int cg0, const int *cgindex)
7571 /* Store the atom block boundaries for easy copying of communication buffers
7574 for (zone = 0; zone < nzone; zone++)
7576 for (p = 0; p < cd->np; p++)
7578 cd->ind[p].cell2at0[zone] = cgindex[cg];
7579 cg += cd->ind[p].nrecv[zone];
7580 cd->ind[p].cell2at1[zone] = cgindex[cg];
7585 static gmx_bool missing_link(t_blocka *link, int cg_gl, char *bLocalCG)
7591 for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
7593 if (!bLocalCG[link->a[i]])
7602 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7604 real c[DIM][4]; /* the corners for the non-bonded communication */
7605 real cr0; /* corner for rounding */
7606 real cr1[4]; /* corners for rounding */
7607 real bc[DIM]; /* corners for bounded communication */
7608 real bcr1; /* corner for rounding for bonded communication */
7611 /* Determine the corners of the domain(s) we are communicating with */
7613 set_dd_corners(const gmx_domdec_t *dd,
7614 int dim0, int dim1, int dim2,
7618 const gmx_domdec_comm_t *comm;
7619 const gmx_domdec_zones_t *zones;
7624 zones = &comm->zones;
7626 /* Keep the compiler happy */
7630 /* The first dimension is equal for all cells */
7631 c->c[0][0] = comm->cell_x0[dim0];
7634 c->bc[0] = c->c[0][0];
7639 /* This cell row is only seen from the first row */
7640 c->c[1][0] = comm->cell_x0[dim1];
7641 /* All rows can see this row */
7642 c->c[1][1] = comm->cell_x0[dim1];
7645 c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
7648 /* For the multi-body distance we need the maximum */
7649 c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
7652 /* Set the upper-right corner for rounding */
7653 c->cr0 = comm->cell_x1[dim0];
7658 for (j = 0; j < 4; j++)
7660 c->c[2][j] = comm->cell_x0[dim2];
7664 /* Use the maximum of the i-cells that see a j-cell */
7665 for (i = 0; i < zones->nizone; i++)
7667 for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
7673 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7679 /* For the multi-body distance we need the maximum */
7680 c->bc[2] = comm->cell_x0[dim2];
7681 for (i = 0; i < 2; i++)
7683 for (j = 0; j < 2; j++)
7685 c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
7691 /* Set the upper-right corner for rounding */
7692 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7693 * Only cell (0,0,0) can see cell 7 (1,1,1)
7695 c->cr1[0] = comm->cell_x1[dim1];
7696 c->cr1[3] = comm->cell_x1[dim1];
7699 c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
7702 /* For the multi-body distance we need the maximum */
7703 c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
7710 /* Determine which cg's we need to send in this pulse from this zone */
7712 get_zone_pulse_cgs(gmx_domdec_t *dd,
7713 int zonei, int zone,
7715 const int *index_gl,
7717 int dim, int dim_ind,
7718 int dim0, int dim1, int dim2,
7719 real r_comm2, real r_bcomm2,
7723 real skew_fac2_d, real skew_fac_01,
7724 rvec *v_d, rvec *v_0, rvec *v_1,
7725 const dd_corners_t *c,
7727 gmx_bool bDistBonded,
7733 gmx_domdec_ind_t *ind,
7734 int **ibuf, int *ibuf_nalloc,
7740 gmx_domdec_comm_t *comm;
7742 gmx_bool bDistMB_pulse;
7744 real r2, rb2, r, tric_sh;
7747 int nsend_z, nsend, nat;
7751 bScrew = (dd->bScrewPBC && dim == XX);
7753 bDistMB_pulse = (bDistMB && bDistBonded);
7759 for (cg = cg0; cg < cg1; cg++)
7763 if (tric_dist[dim_ind] == 0)
7765 /* Rectangular direction, easy */
7766 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7773 r = cg_cm[cg][dim] - c->bc[dim_ind];
7779 /* Rounding gives at most a 16% reduction
7780 * in communicated atoms
7782 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7784 r = cg_cm[cg][dim0] - c->cr0;
7785 /* This is the first dimension, so always r >= 0 */
7792 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7794 r = cg_cm[cg][dim1] - c->cr1[zone];
7801 r = cg_cm[cg][dim1] - c->bcr1;
7811 /* Triclinic direction, more complicated */
7814 /* Rounding, conservative as the skew_fac multiplication
7815 * will slightly underestimate the distance.
7817 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7819 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7820 for (i = dim0+1; i < DIM; i++)
7822 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7824 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7827 rb[dim0] = rn[dim0];
7830 /* Take care that the cell planes along dim0 might not
7831 * be orthogonal to those along dim1 and dim2.
7833 for (i = 1; i <= dim_ind; i++)
7836 if (normal[dim0][dimd] > 0)
7838 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7841 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7846 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7848 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7850 for (i = dim1+1; i < DIM; i++)
7852 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7854 rn[dim1] += tric_sh;
7857 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7858 /* Take care of coupling of the distances
7859 * to the planes along dim0 and dim1 through dim2.
7861 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7862 /* Take care that the cell planes along dim1
7863 * might not be orthogonal to that along dim2.
7865 if (normal[dim1][dim2] > 0)
7867 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7873 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7876 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7877 /* Take care of coupling of the distances
7878 * to the planes along dim0 and dim1 through dim2.
7880 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7881 /* Take care that the cell planes along dim1
7882 * might not be orthogonal to that along dim2.
7884 if (normal[dim1][dim2] > 0)
7886 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7891 /* The distance along the communication direction */
7892 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7894 for (i = dim+1; i < DIM; i++)
7896 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7901 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7902 /* Take care of coupling of the distances
7903 * to the planes along dim0 and dim1 through dim2.
7905 if (dim_ind == 1 && zonei == 1)
7907 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7913 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7916 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7917 /* Take care of coupling of the distances
7918 * to the planes along dim0 and dim1 through dim2.
7920 if (dim_ind == 1 && zonei == 1)
7922 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7930 ((bDistMB && rb2 < r_bcomm2) ||
7931 (bDist2B && r2 < r_bcomm2)) &&
7933 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7934 missing_link(comm->cglink, index_gl[cg],
7937 /* Make an index to the local charge groups */
7938 if (nsend+1 > ind->nalloc)
7940 ind->nalloc = over_alloc_large(nsend+1);
7941 srenew(ind->index, ind->nalloc);
7943 if (nsend+1 > *ibuf_nalloc)
7945 *ibuf_nalloc = over_alloc_large(nsend+1);
7946 srenew(*ibuf, *ibuf_nalloc);
7948 ind->index[nsend] = cg;
7949 (*ibuf)[nsend] = index_gl[cg];
7951 vec_rvec_check_alloc(vbuf, nsend+1);
7953 if (dd->ci[dim] == 0)
7955 /* Correct cg_cm for pbc */
7956 rvec_add(cg_cm[cg], box[dim], vbuf->v[nsend]);
7959 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7960 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7965 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
7968 nat += cgindex[cg+1] - cgindex[cg];
7974 *nsend_z_ptr = nsend_z;
7977 static void setup_dd_communication(gmx_domdec_t *dd,
7978 matrix box, gmx_ddbox_t *ddbox,
7979 t_forcerec *fr, t_state *state, rvec **f)
7981 int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
7982 int nzone, nzone_send, zone, zonei, cg0, cg1;
7983 int c, i, j, cg, cg_gl, nrcg;
7984 int *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
7985 gmx_domdec_comm_t *comm;
7986 gmx_domdec_zones_t *zones;
7987 gmx_domdec_comm_dim_t *cd;
7988 gmx_domdec_ind_t *ind;
7989 cginfo_mb_t *cginfo_mb;
7990 gmx_bool bBondComm, bDist2B, bDistMB, bDistBonded;
7991 real r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
7992 dd_corners_t corners;
7994 rvec *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
7995 real skew_fac2_d, skew_fac_01;
8002 fprintf(debug, "Setting up DD communication\n");
8007 switch (fr->cutoff_scheme)
8016 gmx_incons("unimplemented");
8020 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8022 dim = dd->dim[dim_ind];
8024 /* Check if we need to use triclinic distances */
8025 tric_dist[dim_ind] = 0;
8026 for (i = 0; i <= dim_ind; i++)
8028 if (ddbox->tric_dir[dd->dim[i]])
8030 tric_dist[dim_ind] = 1;
8035 bBondComm = comm->bBondComm;
8037 /* Do we need to determine extra distances for multi-body bondeds? */
8038 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8040 /* Do we need to determine extra distances for only two-body bondeds? */
8041 bDist2B = (bBondComm && !bDistMB);
8043 r_comm2 = sqr(comm->cutoff);
8044 r_bcomm2 = sqr(comm->cutoff_mbody);
8048 fprintf(debug, "bBondComm %d, r_bc %f\n", bBondComm, sqrt(r_bcomm2));
8051 zones = &comm->zones;
8054 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
8055 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
8057 set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
8059 /* Triclinic stuff */
8060 normal = ddbox->normal;
8064 v_0 = ddbox->v[dim0];
8065 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
8067 /* Determine the coupling coefficient for the distances
8068 * to the cell planes along dim0 and dim1 through dim2.
8069 * This is required for correct rounding.
8072 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
8075 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
8081 v_1 = ddbox->v[dim1];
8084 zone_cg_range = zones->cg_range;
8085 index_gl = dd->index_gl;
8086 cgindex = dd->cgindex;
8087 cginfo_mb = fr->cginfo_mb;
8089 zone_cg_range[0] = 0;
8090 zone_cg_range[1] = dd->ncg_home;
8091 comm->zone_ncg1[0] = dd->ncg_home;
8092 pos_cg = dd->ncg_home;
8094 nat_tot = dd->nat_home;
8096 for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
8098 dim = dd->dim[dim_ind];
8099 cd = &comm->cd[dim_ind];
8101 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8103 /* No pbc in this dimension, the first node should not comm. */
8111 v_d = ddbox->v[dim];
8112 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8114 cd->bInPlace = TRUE;
8115 for (p = 0; p < cd->np; p++)
8117 /* Only atoms communicated in the first pulse are used
8118 * for multi-body bonded interactions or for bBondComm.
8120 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8125 for (zone = 0; zone < nzone_send; zone++)
8127 if (tric_dist[dim_ind] && dim_ind > 0)
8129 /* Determine slightly more optimized skew_fac's
8131 * This reduces the number of communicated atoms
8132 * by about 10% for 3D DD of rhombic dodecahedra.
8134 for (dimd = 0; dimd < dim; dimd++)
8136 sf2_round[dimd] = 1;
8137 if (ddbox->tric_dir[dimd])
8139 for (i = dd->dim[dimd]+1; i < DIM; i++)
8141 /* If we are shifted in dimension i
8142 * and the cell plane is tilted forward
8143 * in dimension i, skip this coupling.
8145 if (!(zones->shift[nzone+zone][i] &&
8146 ddbox->v[dimd][i][dimd] >= 0))
8149 sqr(ddbox->v[dimd][i][dimd]);
8152 sf2_round[dimd] = 1/sf2_round[dimd];
8157 zonei = zone_perm[dim_ind][zone];
8160 /* Here we permutate the zones to obtain a convenient order
8161 * for neighbor searching
8163 cg0 = zone_cg_range[zonei];
8164 cg1 = zone_cg_range[zonei+1];
8168 /* Look only at the cg's received in the previous grid pulse
8170 cg1 = zone_cg_range[nzone+zone+1];
8171 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8174 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8175 for (th = 0; th < comm->nth; th++)
8177 gmx_domdec_ind_t *ind_p;
8178 int **ibuf_p, *ibuf_nalloc_p;
8180 int *nsend_p, *nat_p;
8186 /* Thread 0 writes in the comm buffers */
8188 ibuf_p = &comm->buf_int;
8189 ibuf_nalloc_p = &comm->nalloc_int;
8190 vbuf_p = &comm->vbuf;
8193 nsend_zone_p = &ind->nsend[zone];
8197 /* Other threads write into temp buffers */
8198 ind_p = &comm->dth[th].ind;
8199 ibuf_p = &comm->dth[th].ibuf;
8200 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8201 vbuf_p = &comm->dth[th].vbuf;
8202 nsend_p = &comm->dth[th].nsend;
8203 nat_p = &comm->dth[th].nat;
8204 nsend_zone_p = &comm->dth[th].nsend_zone;
8206 comm->dth[th].nsend = 0;
8207 comm->dth[th].nat = 0;
8208 comm->dth[th].nsend_zone = 0;
8218 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8219 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8222 /* Get the cg's for this pulse in this zone */
8223 get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
8225 dim, dim_ind, dim0, dim1, dim2,
8228 normal, skew_fac2_d, skew_fac_01,
8229 v_d, v_0, v_1, &corners, sf2_round,
8230 bDistBonded, bBondComm,
8234 ibuf_p, ibuf_nalloc_p,
8240 /* Append data of threads>=1 to the communication buffers */
8241 for (th = 1; th < comm->nth; th++)
8243 dd_comm_setup_work_t *dth;
8246 dth = &comm->dth[th];
8248 ns1 = nsend + dth->nsend_zone;
8249 if (ns1 > ind->nalloc)
8251 ind->nalloc = over_alloc_dd(ns1);
8252 srenew(ind->index, ind->nalloc);
8254 if (ns1 > comm->nalloc_int)
8256 comm->nalloc_int = over_alloc_dd(ns1);
8257 srenew(comm->buf_int, comm->nalloc_int);
8259 if (ns1 > comm->vbuf.nalloc)
8261 comm->vbuf.nalloc = over_alloc_dd(ns1);
8262 srenew(comm->vbuf.v, comm->vbuf.nalloc);
8265 for (i = 0; i < dth->nsend_zone; i++)
8267 ind->index[nsend] = dth->ind.index[i];
8268 comm->buf_int[nsend] = dth->ibuf[i];
8269 copy_rvec(dth->vbuf.v[i],
8270 comm->vbuf.v[nsend]);
8274 ind->nsend[zone] += dth->nsend_zone;
8277 /* Clear the counts in case we do not have pbc */
8278 for (zone = nzone_send; zone < nzone; zone++)
8280 ind->nsend[zone] = 0;
8282 ind->nsend[nzone] = nsend;
8283 ind->nsend[nzone+1] = nat;
8284 /* Communicate the number of cg's and atoms to receive */
8285 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8286 ind->nsend, nzone+2,
8287 ind->nrecv, nzone+2);
8289 /* The rvec buffer is also required for atom buffers of size nsend
8290 * in dd_move_x and dd_move_f.
8292 vec_rvec_check_alloc(&comm->vbuf, ind->nsend[nzone+1]);
8296 /* We can receive in place if only the last zone is not empty */
8297 for (zone = 0; zone < nzone-1; zone++)
8299 if (ind->nrecv[zone] > 0)
8301 cd->bInPlace = FALSE;
8306 /* The int buffer is only required here for the cg indices */
8307 if (ind->nrecv[nzone] > comm->nalloc_int2)
8309 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8310 srenew(comm->buf_int2, comm->nalloc_int2);
8312 /* The rvec buffer is also required for atom buffers
8313 * of size nrecv in dd_move_x and dd_move_f.
8315 i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
8316 vec_rvec_check_alloc(&comm->vbuf2, i);
8320 /* Make space for the global cg indices */
8321 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8322 || dd->cg_nalloc == 0)
8324 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8325 srenew(index_gl, dd->cg_nalloc);
8326 srenew(cgindex, dd->cg_nalloc+1);
8328 /* Communicate the global cg indices */
8331 recv_i = index_gl + pos_cg;
8335 recv_i = comm->buf_int2;
8337 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8338 comm->buf_int, nsend,
8339 recv_i, ind->nrecv[nzone]);
8341 /* Make space for cg_cm */
8342 dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
8343 if (fr->cutoff_scheme == ecutsGROUP)
8351 /* Communicate cg_cm */
8354 recv_vr = cg_cm + pos_cg;
8358 recv_vr = comm->vbuf2.v;
8360 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8361 comm->vbuf.v, nsend,
8362 recv_vr, ind->nrecv[nzone]);
8364 /* Make the charge group index */
8367 zone = (p == 0 ? 0 : nzone - 1);
8368 while (zone < nzone)
8370 for (cg = 0; cg < ind->nrecv[zone]; cg++)
8372 cg_gl = index_gl[pos_cg];
8373 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
8374 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8375 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8378 /* Update the charge group presence,
8379 * so we can use it in the next pass of the loop.
8381 comm->bLocalCG[cg_gl] = TRUE;
8387 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8390 zone_cg_range[nzone+zone] = pos_cg;
8395 /* This part of the code is never executed with bBondComm. */
8396 merge_cg_buffers(nzone, cd, p, zone_cg_range,
8397 index_gl, recv_i, cg_cm, recv_vr,
8398 cgindex, fr->cginfo_mb, fr->cginfo);
8399 pos_cg += ind->nrecv[nzone];
8401 nat_tot += ind->nrecv[nzone+1];
8405 /* Store the atom block for easy copying of communication buffers */
8406 make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
8410 dd->index_gl = index_gl;
8411 dd->cgindex = cgindex;
8413 dd->ncg_tot = zone_cg_range[zones->n];
8414 dd->nat_tot = nat_tot;
8415 comm->nat[ddnatHOME] = dd->nat_home;
8416 for (i = ddnatZONE; i < ddnatNR; i++)
8418 comm->nat[i] = dd->nat_tot;
8423 /* We don't need to update cginfo, since that was alrady done above.
8424 * So we pass NULL for the forcerec.
8426 dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
8427 NULL, comm->bLocalCG);
8432 fprintf(debug, "Finished setting up DD communication, zones:");
8433 for (c = 0; c < zones->n; c++)
8435 fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
8437 fprintf(debug, "\n");
8441 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8445 for (c = 0; c < zones->nizone; c++)
8447 zones->izone[c].cg1 = zones->cg_range[c+1];
8448 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8449 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8453 static void set_zones_size(gmx_domdec_t *dd,
8454 matrix box, const gmx_ddbox_t *ddbox,
8455 int zone_start, int zone_end)
8457 gmx_domdec_comm_t *comm;
8458 gmx_domdec_zones_t *zones;
8460 int z, zi, zj0, zj1, d, dim;
8463 real size_j, add_tric;
8468 zones = &comm->zones;
8470 /* Do we need to determine extra distances for multi-body bondeds? */
8471 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8473 for (z = zone_start; z < zone_end; z++)
8475 /* Copy cell limits to zone limits.
8476 * Valid for non-DD dims and non-shifted dims.
8478 copy_rvec(comm->cell_x0, zones->size[z].x0);
8479 copy_rvec(comm->cell_x1, zones->size[z].x1);
8482 for (d = 0; d < dd->ndim; d++)
8486 for (z = 0; z < zones->n; z++)
8488 /* With a staggered grid we have different sizes
8489 * for non-shifted dimensions.
8491 if (dd->bGridJump && zones->shift[z][dim] == 0)
8495 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8496 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8500 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8501 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8507 rcmbs = comm->cutoff_mbody;
8508 if (ddbox->tric_dir[dim])
8510 rcs /= ddbox->skew_fac[dim];
8511 rcmbs /= ddbox->skew_fac[dim];
8514 /* Set the lower limit for the shifted zone dimensions */
8515 for (z = zone_start; z < zone_end; z++)
8517 if (zones->shift[z][dim] > 0)
8520 if (!dd->bGridJump || d == 0)
8522 zones->size[z].x0[dim] = comm->cell_x1[dim];
8523 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8527 /* Here we take the lower limit of the zone from
8528 * the lowest domain of the zone below.
8532 zones->size[z].x0[dim] =
8533 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8539 zones->size[z].x0[dim] =
8540 zones->size[zone_perm[2][z-4]].x0[dim];
8544 zones->size[z].x0[dim] =
8545 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8548 /* A temporary limit, is updated below */
8549 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8553 for (zi = 0; zi < zones->nizone; zi++)
8555 if (zones->shift[zi][dim] == 0)
8557 /* This takes the whole zone into account.
8558 * With multiple pulses this will lead
8559 * to a larger zone then strictly necessary.
8561 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8562 zones->size[zi].x1[dim]+rcmbs);
8570 /* Loop over the i-zones to set the upper limit of each
8573 for (zi = 0; zi < zones->nizone; zi++)
8575 if (zones->shift[zi][dim] == 0)
8577 for (z = zones->izone[zi].j0; z < zones->izone[zi].j1; z++)
8579 if (zones->shift[z][dim] > 0)
8581 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8582 zones->size[zi].x1[dim]+rcs);
8589 for (z = zone_start; z < zone_end; z++)
8591 /* Initialization only required to keep the compiler happy */
8592 rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
8595 /* To determine the bounding box for a zone we need to find
8596 * the extreme corners of 4, 2 or 1 corners.
8598 nc = 1 << (ddbox->npbcdim - 1);
8600 for (c = 0; c < nc; c++)
8602 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8606 corner[YY] = zones->size[z].x0[YY];
8610 corner[YY] = zones->size[z].x1[YY];
8614 corner[ZZ] = zones->size[z].x0[ZZ];
8618 corner[ZZ] = zones->size[z].x1[ZZ];
8620 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8622 /* With 1D domain decomposition the cg's are not in
8623 * the triclinic box, but triclinic x-y and rectangular y-z.
8624 * Shift y back, so it will later end up at 0.
8626 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8628 /* Apply the triclinic couplings */
8629 for (i = YY; i < ddbox->npbcdim; i++)
8631 for (j = XX; j < i; j++)
8633 corner[j] += corner[i]*box[i][j]/box[i][i];
8638 copy_rvec(corner, corner_min);
8639 copy_rvec(corner, corner_max);
8643 for (i = 0; i < DIM; i++)
8645 corner_min[i] = min(corner_min[i], corner[i]);
8646 corner_max[i] = max(corner_max[i], corner[i]);
8650 /* Copy the extreme cornes without offset along x */
8651 for (i = 0; i < DIM; i++)
8653 zones->size[z].bb_x0[i] = corner_min[i];
8654 zones->size[z].bb_x1[i] = corner_max[i];
8656 /* Add the offset along x */
8657 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8658 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8661 if (zone_start == 0)
8664 for (dim = 0; dim < DIM; dim++)
8666 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8668 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8673 for (z = zone_start; z < zone_end; z++)
8675 fprintf(debug, "zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8677 zones->size[z].x0[XX], zones->size[z].x1[XX],
8678 zones->size[z].x0[YY], zones->size[z].x1[YY],
8679 zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
8680 fprintf(debug, "zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8682 zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
8683 zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
8684 zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
8689 static int comp_cgsort(const void *a, const void *b)
8693 gmx_cgsort_t *cga, *cgb;
8694 cga = (gmx_cgsort_t *)a;
8695 cgb = (gmx_cgsort_t *)b;
8697 comp = cga->nsc - cgb->nsc;
8700 comp = cga->ind_gl - cgb->ind_gl;
8706 static void order_int_cg(int n, const gmx_cgsort_t *sort,
8711 /* Order the data */
8712 for (i = 0; i < n; i++)
8714 buf[i] = a[sort[i].ind];
8717 /* Copy back to the original array */
8718 for (i = 0; i < n; i++)
8724 static void order_vec_cg(int n, const gmx_cgsort_t *sort,
8729 /* Order the data */
8730 for (i = 0; i < n; i++)
8732 copy_rvec(v[sort[i].ind], buf[i]);
8735 /* Copy back to the original array */
8736 for (i = 0; i < n; i++)
8738 copy_rvec(buf[i], v[i]);
8742 static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
8745 int a, atot, cg, cg0, cg1, i;
8747 if (cgindex == NULL)
8749 /* Avoid the useless loop of the atoms within a cg */
8750 order_vec_cg(ncg, sort, v, buf);
8755 /* Order the data */
8757 for (cg = 0; cg < ncg; cg++)
8759 cg0 = cgindex[sort[cg].ind];
8760 cg1 = cgindex[sort[cg].ind+1];
8761 for (i = cg0; i < cg1; i++)
8763 copy_rvec(v[i], buf[a]);
8769 /* Copy back to the original array */
8770 for (a = 0; a < atot; a++)
8772 copy_rvec(buf[a], v[a]);
8776 static void ordered_sort(int nsort2, gmx_cgsort_t *sort2,
8777 int nsort_new, gmx_cgsort_t *sort_new,
8778 gmx_cgsort_t *sort1)
8782 /* The new indices are not very ordered, so we qsort them */
8783 qsort_threadsafe(sort_new, nsort_new, sizeof(sort_new[0]), comp_cgsort);
8785 /* sort2 is already ordered, so now we can merge the two arrays */
8789 while (i2 < nsort2 || i_new < nsort_new)
8793 sort1[i1++] = sort_new[i_new++];
8795 else if (i_new == nsort_new)
8797 sort1[i1++] = sort2[i2++];
8799 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8800 (sort2[i2].nsc == sort_new[i_new].nsc &&
8801 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8803 sort1[i1++] = sort2[i2++];
8807 sort1[i1++] = sort_new[i_new++];
8812 static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
8814 gmx_domdec_sort_t *sort;
8815 gmx_cgsort_t *cgsort, *sort_i;
8816 int ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
8817 int sort_last, sort_skip;
8819 sort = dd->comm->sort;
8821 a = fr->ns.grid->cell_index;
8823 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8825 if (ncg_home_old >= 0)
8827 /* The charge groups that remained in the same ns grid cell
8828 * are completely ordered. So we can sort efficiently by sorting
8829 * the charge groups that did move into the stationary list.
8834 for (i = 0; i < dd->ncg_home; i++)
8836 /* Check if this cg did not move to another node */
8839 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8841 /* This cg is new on this node or moved ns grid cell */
8842 if (nsort_new >= sort->sort_new_nalloc)
8844 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8845 srenew(sort->sort_new, sort->sort_new_nalloc);
8847 sort_i = &(sort->sort_new[nsort_new++]);
8851 /* This cg did not move */
8852 sort_i = &(sort->sort2[nsort2++]);
8854 /* Sort on the ns grid cell indices
8855 * and the global topology index.
8856 * index_gl is irrelevant with cell ns,
8857 * but we set it here anyhow to avoid a conditional.
8860 sort_i->ind_gl = dd->index_gl[i];
8867 fprintf(debug, "ordered sort cgs: stationary %d moved %d\n",
8870 /* Sort efficiently */
8871 ordered_sort(nsort2, sort->sort2, nsort_new, sort->sort_new,
8876 cgsort = sort->sort;
8878 for (i = 0; i < dd->ncg_home; i++)
8880 /* Sort on the ns grid cell indices
8881 * and the global topology index
8883 cgsort[i].nsc = a[i];
8884 cgsort[i].ind_gl = dd->index_gl[i];
8886 if (cgsort[i].nsc < moved)
8893 fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, ncg_new);
8895 /* Determine the order of the charge groups using qsort */
8896 qsort_threadsafe(cgsort, dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
8902 static int dd_sort_order_nbnxn(gmx_domdec_t *dd, t_forcerec *fr)
8905 int ncg_new, i, *a, na;
8907 sort = dd->comm->sort->sort;
8909 nbnxn_get_atomorder(fr->nbv->nbs, &a, &na);
8912 for (i = 0; i < na; i++)
8916 sort[ncg_new].ind = a[i];
8924 static void dd_sort_state(gmx_domdec_t *dd, int ePBC,
8925 rvec *cgcm, t_forcerec *fr, t_state *state,
8928 gmx_domdec_sort_t *sort;
8929 gmx_cgsort_t *cgsort, *sort_i;
8931 int ncg_new, i, *ibuf, cgsize;
8934 sort = dd->comm->sort;
8936 if (dd->ncg_home > sort->sort_nalloc)
8938 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8939 srenew(sort->sort, sort->sort_nalloc);
8940 srenew(sort->sort2, sort->sort_nalloc);
8942 cgsort = sort->sort;
8944 switch (fr->cutoff_scheme)
8947 ncg_new = dd_sort_order(dd, fr, ncg_home_old);
8950 ncg_new = dd_sort_order_nbnxn(dd, fr);
8953 gmx_incons("unimplemented");
8957 /* We alloc with the old size, since cgindex is still old */
8958 vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
8959 vbuf = dd->comm->vbuf.v;
8963 cgindex = dd->cgindex;
8970 /* Remove the charge groups which are no longer at home here */
8971 dd->ncg_home = ncg_new;
8974 fprintf(debug, "Set the new home charge group count to %d\n",
8978 /* Reorder the state */
8979 for (i = 0; i < estNR; i++)
8981 if (EST_DISTR(i) && (state->flags & (1<<i)))
8986 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
8989 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
8992 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->sd_X, vbuf);
8995 order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
8999 case estDISRE_INITF:
9000 case estDISRE_RM3TAV:
9001 case estORIRE_INITF:
9003 /* No ordering required */
9006 gmx_incons("Unknown state entry encountered in dd_sort_state");
9011 if (fr->cutoff_scheme == ecutsGROUP)
9014 order_vec_cg(dd->ncg_home, cgsort, cgcm, vbuf);
9017 if (dd->ncg_home+1 > sort->ibuf_nalloc)
9019 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
9020 srenew(sort->ibuf, sort->ibuf_nalloc);
9023 /* Reorder the global cg index */
9024 order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
9025 /* Reorder the cginfo */
9026 order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
9027 /* Rebuild the local cg index */
9031 for (i = 0; i < dd->ncg_home; i++)
9033 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
9034 ibuf[i+1] = ibuf[i] + cgsize;
9036 for (i = 0; i < dd->ncg_home+1; i++)
9038 dd->cgindex[i] = ibuf[i];
9043 for (i = 0; i < dd->ncg_home+1; i++)
9048 /* Set the home atom number */
9049 dd->nat_home = dd->cgindex[dd->ncg_home];
9051 if (fr->cutoff_scheme == ecutsVERLET)
9053 /* The atoms are now exactly in grid order, update the grid order */
9054 nbnxn_set_atomorder(fr->nbv->nbs);
9058 /* Copy the sorted ns cell indices back to the ns grid struct */
9059 for (i = 0; i < dd->ncg_home; i++)
9061 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
9063 fr->ns.grid->nr = dd->ncg_home;
9067 static void add_dd_statistics(gmx_domdec_t *dd)
9069 gmx_domdec_comm_t *comm;
9074 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9076 comm->sum_nat[ddnat-ddnatZONE] +=
9077 comm->nat[ddnat] - comm->nat[ddnat-1];
9082 void reset_dd_statistics_counters(gmx_domdec_t *dd)
9084 gmx_domdec_comm_t *comm;
9089 /* Reset all the statistics and counters for total run counting */
9090 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9092 comm->sum_nat[ddnat-ddnatZONE] = 0;
9096 comm->load_step = 0;
9099 clear_ivec(comm->load_lim);
9104 void print_dd_statistics(t_commrec *cr, t_inputrec *ir, FILE *fplog)
9106 gmx_domdec_comm_t *comm;
9110 comm = cr->dd->comm;
9112 gmx_sumd(ddnatNR-ddnatZONE, comm->sum_nat, cr);
9119 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");
9121 for (ddnat = ddnatZONE; ddnat < ddnatNR; ddnat++)
9123 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9128 " av. #atoms communicated per step for force: %d x %.1f\n",
9132 if (cr->dd->vsite_comm)
9135 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9136 (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
9141 if (cr->dd->constraint_comm)
9144 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9145 1 + ir->nLincsIter, av);
9149 gmx_incons(" Unknown type for DD statistics");
9152 fprintf(fplog, "\n");
9154 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9156 print_dd_load_av(fplog, cr->dd);
9160 void dd_partition_system(FILE *fplog,
9161 gmx_large_int_t step,
9163 gmx_bool bMasterState,
9165 t_state *state_global,
9166 gmx_mtop_t *top_global,
9168 t_state *state_local,
9171 gmx_localtop_t *top_local,
9174 gmx_shellfc_t shellfc,
9175 gmx_constr_t constr,
9177 gmx_wallcycle_t wcycle,
9181 gmx_domdec_comm_t *comm;
9182 gmx_ddbox_t ddbox = {0};
9184 gmx_large_int_t step_pcoupl;
9185 rvec cell_ns_x0, cell_ns_x1;
9186 int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
9187 gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
9188 gmx_bool bRedist, bSortCG, bResortAll;
9189 ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
9196 bBoxChanged = (bMasterState || DEFORM(*ir));
9197 if (ir->epc != epcNO)
9199 /* With nstpcouple > 1 pressure coupling happens.
9200 * one step after calculating the pressure.
9201 * Box scaling happens at the end of the MD step,
9202 * after the DD partitioning.
9203 * We therefore have to do DLB in the first partitioning
9204 * after an MD step where P-coupling occured.
9205 * We need to determine the last step in which p-coupling occurred.
9206 * MRS -- need to validate this for vv?
9211 step_pcoupl = step - 1;
9215 step_pcoupl = ((step - 1)/n)*n + 1;
9217 if (step_pcoupl >= comm->partition_step)
9223 bNStGlobalComm = (step % nstglobalcomm == 0);
9225 if (!comm->bDynLoadBal)
9231 /* Should we do dynamic load balacing this step?
9232 * Since it requires (possibly expensive) global communication,
9233 * we might want to do DLB less frequently.
9235 if (bBoxChanged || ir->epc != epcNO)
9237 bDoDLB = bBoxChanged;
9241 bDoDLB = bNStGlobalComm;
9245 /* Check if we have recorded loads on the nodes */
9246 if (comm->bRecordLoad && dd_load_count(comm))
9248 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9250 /* Check if we should use DLB at the second partitioning
9251 * and every 100 partitionings,
9252 * so the extra communication cost is negligible.
9254 n = max(100, nstglobalcomm);
9255 bCheckDLB = (comm->n_load_collect == 0 ||
9256 comm->n_load_have % n == n-1);
9263 /* Print load every nstlog, first and last step to the log file */
9264 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9265 comm->n_load_collect == 0 ||
9267 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9269 /* Avoid extra communication due to verbose screen output
9270 * when nstglobalcomm is set.
9272 if (bDoDLB || bLogLoad || bCheckDLB ||
9273 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9275 get_load_distribution(dd, wcycle);
9280 dd_print_load(fplog, dd, step-1);
9284 dd_print_load_verbose(dd);
9287 comm->n_load_collect++;
9291 /* Since the timings are node dependent, the master decides */
9295 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9298 fprintf(debug, "step %s, imb loss %f\n",
9299 gmx_step_str(step, sbuf),
9300 dd_force_imb_perf_loss(dd));
9303 dd_bcast(dd, sizeof(bTurnOnDLB), &bTurnOnDLB);
9306 turn_on_dlb(fplog, cr, step);
9311 comm->n_load_have++;
9314 cgs_gl = &comm->cgs_gl;
9319 /* Clear the old state */
9320 clear_dd_indices(dd, 0, 0);
9323 set_ddbox(dd, bMasterState, cr, ir, state_global->box,
9324 TRUE, cgs_gl, state_global->x, &ddbox);
9326 get_cg_distribution(fplog, step, dd, cgs_gl,
9327 state_global->box, &ddbox, state_global->x);
9329 dd_distribute_state(dd, cgs_gl,
9330 state_global, state_local, f);
9332 dd_make_local_cgs(dd, &top_local->cgs);
9334 /* Ensure that we have space for the new distribution */
9335 dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
9337 if (fr->cutoff_scheme == ecutsGROUP)
9339 calc_cgcm(fplog, 0, dd->ncg_home,
9340 &top_local->cgs, state_local->x, fr->cg_cm);
9343 inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
9345 dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
9347 else if (state_local->ddp_count != dd->ddp_count)
9349 if (state_local->ddp_count > dd->ddp_count)
9351 gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)", state_local->ddp_count, dd->ddp_count);
9354 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9356 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);
9359 /* Clear the old state */
9360 clear_dd_indices(dd, 0, 0);
9362 /* Build the new indices */
9363 rebuild_cgindex(dd, cgs_gl->index, state_local);
9364 make_dd_indices(dd, cgs_gl->index, 0);
9365 ncgindex_set = dd->ncg_home;
9367 if (fr->cutoff_scheme == ecutsGROUP)
9369 /* Redetermine the cg COMs */
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 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9379 TRUE, &top_local->cgs, state_local->x, &ddbox);
9381 bRedist = comm->bDynLoadBal;
9385 /* We have the full state, only redistribute the cgs */
9387 /* Clear the non-home indices */
9388 clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
9391 /* Avoid global communication for dim's without pbc and -gcom */
9392 if (!bNStGlobalComm)
9394 copy_rvec(comm->box0, ddbox.box0 );
9395 copy_rvec(comm->box_size, ddbox.box_size);
9397 set_ddbox(dd, bMasterState, cr, ir, state_local->box,
9398 bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
9403 /* For dim's without pbc and -gcom */
9404 copy_rvec(ddbox.box0, comm->box0 );
9405 copy_rvec(ddbox.box_size, comm->box_size);
9407 set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
9410 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9412 write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
9415 /* Check if we should sort the charge groups */
9416 if (comm->nstSortCG > 0)
9418 bSortCG = (bMasterState ||
9419 (bRedist && (step % comm->nstSortCG == 0)));
9426 ncg_home_old = dd->ncg_home;
9431 wallcycle_sub_start(wcycle, ewcsDD_REDIST);
9433 dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
9434 state_local, f, fr, mdatoms,
9435 !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
9437 wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
9440 get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
9442 &comm->cell_x0, &comm->cell_x1,
9443 dd->ncg_home, fr->cg_cm,
9444 cell_ns_x0, cell_ns_x1, &grid_density);
9448 comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
9451 switch (fr->cutoff_scheme)
9454 copy_ivec(fr->ns.grid->n, ncells_old);
9455 grid_first(fplog, fr->ns.grid, dd, &ddbox, fr->ePBC,
9456 state_local->box, cell_ns_x0, cell_ns_x1,
9457 fr->rlistlong, grid_density);
9460 nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
9463 gmx_incons("unimplemented");
9465 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9466 copy_ivec(ddbox.tric_dir, comm->tric_dir);
9470 wallcycle_sub_start(wcycle, ewcsDD_GRID);
9472 /* Sort the state on charge group position.
9473 * This enables exact restarts from this step.
9474 * It also improves performance by about 15% with larger numbers
9475 * of atoms per node.
9478 /* Fill the ns grid with the home cell,
9479 * so we can sort with the indices.
9481 set_zones_ncg_home(dd);
9483 switch (fr->cutoff_scheme)
9486 set_zones_size(dd, state_local->box, &ddbox, 0, 1);
9488 nbnxn_put_on_grid(fr->nbv->nbs, fr->ePBC, state_local->box,
9490 comm->zones.size[0].bb_x0,
9491 comm->zones.size[0].bb_x1,
9493 comm->zones.dens_zone0,
9496 ncg_moved, bRedist ? comm->moved : NULL,
9497 fr->nbv->grp[eintLocal].kernel_type,
9498 fr->nbv->grp[eintLocal].nbat);
9500 nbnxn_get_ncells(fr->nbv->nbs, &ncells_new[XX], &ncells_new[YY]);
9503 fill_grid(fplog, &comm->zones, fr->ns.grid, dd->ncg_home,
9504 0, dd->ncg_home, fr->cg_cm);
9506 copy_ivec(fr->ns.grid->n, ncells_new);
9509 gmx_incons("unimplemented");
9512 bResortAll = bMasterState;
9514 /* Check if we can user the old order and ns grid cell indices
9515 * of the charge groups to sort the charge groups efficiently.
9517 if (ncells_new[XX] != ncells_old[XX] ||
9518 ncells_new[YY] != ncells_old[YY] ||
9519 ncells_new[ZZ] != ncells_old[ZZ])
9526 fprintf(debug, "Step %s, sorting the %d home charge groups\n",
9527 gmx_step_str(step, sbuf), dd->ncg_home);
9529 dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
9530 bResortAll ? -1 : ncg_home_old);
9531 /* Rebuild all the indices */
9532 ga2la_clear(dd->ga2la);
9535 wallcycle_sub_stop(wcycle, ewcsDD_GRID);
9538 wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
9540 /* Setup up the communication and communicate the coordinates */
9541 setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
9543 /* Set the indices */
9544 make_dd_indices(dd, cgs_gl->index, ncgindex_set);
9546 /* Set the charge group boundaries for neighbor searching */
9547 set_cg_boundaries(&comm->zones);
9549 if (fr->cutoff_scheme == ecutsVERLET)
9551 set_zones_size(dd, state_local->box, &ddbox,
9552 bSortCG ? 1 : 0, comm->zones.n);
9555 wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
9558 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9559 -1,state_local->x,state_local->box);
9562 wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
9564 /* Extract a local topology from the global topology */
9565 for (i = 0; i < dd->ndim; i++)
9567 np[dd->dim[i]] = comm->cd[i].np;
9569 dd_make_local_top(fplog, dd, &comm->zones, dd->npbcdim, state_local->box,
9570 comm->cellsize_min, np,
9572 fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
9573 vsite, top_global, top_local);
9575 wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
9577 wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
9579 /* Set up the special atom communication */
9580 n = comm->nat[ddnatZONE];
9581 for (i = ddnatZONE+1; i < ddnatNR; i++)
9586 if (vsite && vsite->n_intercg_vsite)
9588 n = dd_make_local_vsites(dd, n, top_local->idef.il);
9592 if (dd->bInterCGcons || dd->bInterCGsettles)
9594 /* Only for inter-cg constraints we need special code */
9595 n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
9596 constr, ir->nProjOrder,
9597 top_local->idef.il);
9601 gmx_incons("Unknown special atom type setup");
9606 wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
9608 wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
9610 /* Make space for the extra coordinates for virtual site
9611 * or constraint communication.
9613 state_local->natoms = comm->nat[ddnatNR-1];
9614 if (state_local->natoms > state_local->nalloc)
9616 dd_realloc_state(state_local, f, state_local->natoms);
9619 if (fr->bF_NoVirSum)
9621 if (vsite && vsite->n_intercg_vsite)
9623 nat_f_novirsum = comm->nat[ddnatVSITE];
9627 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9629 nat_f_novirsum = dd->nat_tot;
9633 nat_f_novirsum = dd->nat_home;
9642 /* Set the number of atoms required for the force calculation.
9643 * Forces need to be constrained when using a twin-range setup
9644 * or with energy minimization. For simple simulations we could
9645 * avoid some allocation, zeroing and copying, but this is
9646 * probably not worth the complications ande checking.
9648 forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
9649 dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
9651 /* We make the all mdatoms up to nat_tot_con.
9652 * We could save some work by only setting invmass
9653 * between nat_tot and nat_tot_con.
9655 /* This call also sets the new number of home particles to dd->nat_home */
9656 atoms2md(top_global, ir,
9657 comm->nat[ddnatCON], dd->gatindex, 0, dd->nat_home, mdatoms);
9659 /* Now we have the charges we can sort the FE interactions */
9660 dd_sort_local_top(dd, mdatoms, top_local);
9664 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9665 split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
9670 /* Make the local shell stuff, currently no communication is done */
9671 make_local_shells(cr, mdatoms, shellfc);
9674 if (ir->implicit_solvent)
9676 make_local_gb(cr, fr->born, ir->gb_algorithm);
9679 init_bonded_thread_force_reduction(fr, &top_local->idef);
9681 if (!(cr->duty & DUTY_PME))
9683 /* Send the charges to our PME only node */
9684 gmx_pme_send_q(cr, mdatoms->nChargePerturbed,
9685 mdatoms->chargeA, mdatoms->chargeB,
9686 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
9691 set_constraints(constr, top_local, ir, mdatoms, cr);
9694 if (ir->ePull != epullNO)
9696 /* Update the local pull groups */
9697 dd_make_local_pull_groups(dd, ir->pull, mdatoms);
9702 /* Update the local rotation groups */
9703 dd_make_local_rotation_groups(dd, ir->rot);
9707 add_dd_statistics(dd);
9709 /* Make sure we only count the cycles for this DD partitioning */
9710 clear_dd_cycle_counts(dd);
9712 /* Because the order of the atoms might have changed since
9713 * the last vsite construction, we need to communicate the constructing
9714 * atom coordinates again (for spreading the forces this MD step).
9716 dd_move_x_vsites(dd, state_local->box, state_local->x);
9718 wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
9720 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9722 dd_move_x(dd, state_local->box, state_local->x);
9723 write_dd_pdb("dd_dump", step, "dump", top_global, cr,
9724 -1, state_local->x, state_local->box);
9727 /* Store the partitioning step */
9728 comm->partition_step = step;
9730 /* Increase the DD partitioning counter */
9732 /* The state currently matches this DD partitioning count, store it */
9733 state_local->ddp_count = dd->ddp_count;
9736 /* The DD master node knows the complete cg distribution,
9737 * store the count so we can possibly skip the cg info communication.
9739 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9742 if (comm->DD_debug > 0)
9744 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9745 check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
9746 "after partitioning");