1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
30 #include "gmx_fatal.h"
31 #include "gmx_fatal_collective.h"
34 #include "domdec_network.h"
37 #include "chargegroup.h"
46 #include "pull_rotation.h"
47 #include "gmx_wallcycle.h"
51 #include "mtop_util.h"
53 #include "gmx_ga2la.h"
56 #include "nbnxn_search.h"
58 #include "gmx_omp_nthreads.h"
67 #define DDRANK(dd,rank) (rank)
68 #define DDMASTERRANK(dd) (dd->masterrank)
70 typedef struct gmx_domdec_master
72 /* The cell boundaries */
74 /* The global charge group division */
75 int *ncg; /* Number of home charge groups for each node */
76 int *index; /* Index of nnodes+1 into cg */
77 int *cg; /* Global charge group index */
78 int *nat; /* Number of home atoms for each node. */
79 int *ibuf; /* Buffer for communication */
80 rvec *vbuf; /* Buffer for state scattering and gathering */
81 } gmx_domdec_master_t;
85 /* The numbers of charge groups to send and receive for each cell
86 * that requires communication, the last entry contains the total
87 * number of atoms that needs to be communicated.
89 int nsend[DD_MAXIZONE+2];
90 int nrecv[DD_MAXIZONE+2];
91 /* The charge groups to send */
94 /* The atom range for non-in-place communication */
95 int cell2at0[DD_MAXIZONE];
96 int cell2at1[DD_MAXIZONE];
101 int np; /* Number of grid pulses in this dimension */
102 int np_dlb; /* For dlb, for use with edlbAUTO */
103 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
105 gmx_bool bInPlace; /* Can we communicate in place? */
106 } gmx_domdec_comm_dim_t;
110 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
111 real *cell_f; /* State var.: cell boundaries, box relative */
112 real *old_cell_f; /* Temp. var.: old cell size */
113 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
114 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
115 real *bound_min; /* Temp. var.: lower limit for cell boundary */
116 real *bound_max; /* Temp. var.: upper limit for cell boundary */
117 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
118 real *buf_ncd; /* Temp. var. */
121 #define DD_NLOAD_MAX 9
123 /* Here floats are accurate enough, since these variables
124 * only influence the load balancing, not the actual MD results.
151 gmx_cgsort_t *sort_new;
163 /* This enum determines the order of the coordinates.
164 * ddnatHOME and ddnatZONE should be first and second,
165 * the others can be ordered as wanted.
167 enum { ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR };
169 enum { edlbAUTO, edlbNO, edlbYES, edlbNR };
170 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
174 int dim; /* The dimension */
175 gmx_bool dim_match;/* Tells if DD and PME dims match */
176 int nslab; /* The number of PME slabs in this dimension */
177 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
178 int *pp_min; /* The minimum pp node location, size nslab */
179 int *pp_max; /* The maximum pp node location,size nslab */
180 int maxshift; /* The maximum shift for coordinate redistribution in PME */
185 real min0; /* The minimum bottom of this zone */
186 real max1; /* The maximum top of this zone */
187 real min1; /* The minimum top of this zone */
188 real mch0; /* The maximum bottom communicaton height for this zone */
189 real mch1; /* The maximum top communicaton height for this zone */
190 real p1_0; /* The bottom value of the first cell in this zone */
191 real p1_1; /* The top value of the first cell in this zone */
196 gmx_domdec_ind_t ind;
203 } dd_comm_setup_work_t;
205 typedef struct gmx_domdec_comm
207 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
208 * unless stated otherwise.
211 /* The number of decomposition dimensions for PME, 0: no PME */
213 /* The number of nodes doing PME (PP/PME or only PME) */
217 /* The communication setup including the PME only nodes */
218 gmx_bool bCartesianPP_PME;
221 int *pmenodes; /* size npmenodes */
222 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
223 * but with bCartesianPP_PME */
224 gmx_ddpme_t ddpme[2];
226 /* The DD particle-particle nodes only */
227 gmx_bool bCartesianPP;
228 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
230 /* The global charge groups */
233 /* Should we sort the cgs */
235 gmx_domdec_sort_t *sort;
237 /* Are there charge groups? */
240 /* Are there bonded and multi-body interactions between charge groups? */
241 gmx_bool bInterCGBondeds;
242 gmx_bool bInterCGMultiBody;
244 /* Data for the optional bonded interaction atom communication range */
251 /* Are we actually using DLB? */
252 gmx_bool bDynLoadBal;
254 /* Cell sizes for static load balancing, first index cartesian */
257 /* The width of the communicated boundaries */
260 /* The minimum cell size (including triclinic correction) */
262 /* For dlb, for use with edlbAUTO */
263 rvec cellsize_min_dlb;
264 /* The lower limit for the DD cell size with DLB */
266 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
267 gmx_bool bVacDLBNoLimit;
269 /* With PME load balancing we set limits on DLB */
270 gmx_bool bPMELoadBalDLBLimits;
271 /* DLB needs to take into account that we want to allow this maximum
272 * cut-off (for PME load balancing), this could limit cell boundaries.
274 real PMELoadBal_max_cutoff;
276 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
278 /* box0 and box_size are required with dim's without pbc and -gcom */
282 /* The cell boundaries */
286 /* The old location of the cell boundaries, to check cg displacements */
290 /* The communication setup and charge group boundaries for the zones */
291 gmx_domdec_zones_t zones;
293 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
294 * cell boundaries of neighboring cells for dynamic load balancing.
296 gmx_ddzone_t zone_d1[2];
297 gmx_ddzone_t zone_d2[2][2];
299 /* The coordinate/force communication setup and indices */
300 gmx_domdec_comm_dim_t cd[DIM];
301 /* The maximum number of cells to communicate with in one dimension */
304 /* Which cg distribution is stored on the master node */
305 int master_cg_ddp_count;
307 /* The number of cg's received from the direct neighbors */
308 int zone_ncg1[DD_MAXZONE];
310 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
313 /* Array for signalling if atoms have moved to another domain */
317 /* Communication buffer for general use */
321 /* Communication buffer for general use */
324 /* Temporary storage for thread parallel communication setup */
326 dd_comm_setup_work_t *dth;
328 /* Communication buffers only used with multiple grid pulses */
333 /* Communication buffers for local redistribution */
335 int cggl_flag_nalloc[DIM*2];
337 int cgcm_state_nalloc[DIM*2];
339 /* Cell sizes for dynamic load balancing */
340 gmx_domdec_root_t **root;
344 real cell_f_max0[DIM];
345 real cell_f_min1[DIM];
347 /* Stuff for load communication */
348 gmx_bool bRecordLoad;
349 gmx_domdec_load_t *load;
351 MPI_Comm *mpi_comm_load;
354 /* Maximum DLB scaling per load balancing step in percent */
358 float cycl[ddCyclNr];
359 int cycl_n[ddCyclNr];
360 float cycl_max[ddCyclNr];
361 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
365 /* Have often have did we have load measurements */
367 /* Have often have we collected the load measurements */
371 double sum_nat[ddnatNR-ddnatZONE];
381 /* The last partition step */
382 gmx_large_int_t partition_step;
390 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
393 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
394 #define DD_FLAG_NRCG 65535
395 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
396 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
398 /* Zone permutation required to obtain consecutive charge groups
399 * for neighbor searching.
401 static const int zone_perm[3][4] = { {0,0,0,0},{1,0,0,0},{3,0,1,2} };
403 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
404 * components see only j zones with that component 0.
407 /* The DD zone order */
408 static const ivec dd_zo[DD_MAXZONE] =
409 {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,1,1},{0,0,1},{1,0,1},{1,1,1}};
414 static const ivec dd_zp3[dd_zp3n] = {{0,0,8},{1,3,6},{2,5,6},{3,5,7}};
419 static const ivec dd_zp2[dd_zp2n] = {{0,0,4},{1,3,4}};
424 static const ivec dd_zp1[dd_zp1n] = {{0,0,2}};
426 /* Factors used to avoid problems due to rounding issues */
427 #define DD_CELL_MARGIN 1.0001
428 #define DD_CELL_MARGIN2 1.00005
429 /* Factor to account for pressure scaling during nstlist steps */
430 #define DD_PRES_SCALE_MARGIN 1.02
432 /* Allowed performance loss before we DLB or warn */
433 #define DD_PERF_LOSS 0.05
435 #define DD_CELL_F_SIZE(dd,di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
437 /* Use separate MPI send and receive commands
438 * when nnodes <= GMX_DD_NNODES_SENDRECV.
439 * This saves memory (and some copying for small nnodes).
440 * For high parallelization scatter and gather calls are used.
442 #define GMX_DD_NNODES_SENDRECV 4
446 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
448 static void index2xyz(ivec nc,int ind,ivec xyz)
450 xyz[XX] = ind % nc[XX];
451 xyz[YY] = (ind / nc[XX]) % nc[YY];
452 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
456 /* This order is required to minimize the coordinate communication in PME
457 * which uses decomposition in the x direction.
459 #define dd_index(n,i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
461 static void ddindex2xyz(ivec nc,int ind,ivec xyz)
463 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
464 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
465 xyz[ZZ] = ind % nc[ZZ];
468 static int ddcoord2ddnodeid(gmx_domdec_t *dd,ivec c)
473 ddindex = dd_index(dd->nc,c);
474 if (dd->comm->bCartesianPP_PME)
476 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
478 else if (dd->comm->bCartesianPP)
481 MPI_Cart_rank(dd->mpi_comm_all,c,&ddnodeid);
492 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox,t_inputrec *ir)
494 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
497 int ddglatnr(gmx_domdec_t *dd,int i)
507 if (i >= dd->comm->nat[ddnatNR-1])
509 gmx_fatal(FARGS,"glatnr called with %d, which is larger than the local number of atoms (%d)",i,dd->comm->nat[ddnatNR-1]);
511 atnr = dd->gatindex[i] + 1;
517 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
519 return &dd->comm->cgs_gl;
522 static void vec_rvec_init(vec_rvec_t *v)
528 static void vec_rvec_check_alloc(vec_rvec_t *v,int n)
532 v->nalloc = over_alloc_dd(n);
533 srenew(v->v,v->nalloc);
537 void dd_store_state(gmx_domdec_t *dd,t_state *state)
541 if (state->ddp_count != dd->ddp_count)
543 gmx_incons("The state does not the domain decomposition state");
546 state->ncg_gl = dd->ncg_home;
547 if (state->ncg_gl > state->cg_gl_nalloc)
549 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
550 srenew(state->cg_gl,state->cg_gl_nalloc);
552 for(i=0; i<state->ncg_gl; i++)
554 state->cg_gl[i] = dd->index_gl[i];
557 state->ddp_count_cg_gl = dd->ddp_count;
560 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
562 return &dd->comm->zones;
565 void dd_get_ns_ranges(gmx_domdec_t *dd,int icg,
566 int *jcg0,int *jcg1,ivec shift0,ivec shift1)
568 gmx_domdec_zones_t *zones;
571 zones = &dd->comm->zones;
574 while (icg >= zones->izone[izone].cg1)
583 else if (izone < zones->nizone)
585 *jcg0 = zones->izone[izone].jcg0;
589 gmx_fatal(FARGS,"DD icg %d out of range: izone (%d) >= nizone (%d)",
590 icg,izone,zones->nizone);
593 *jcg1 = zones->izone[izone].jcg1;
595 for(d=0; d<dd->ndim; d++)
598 shift0[dim] = zones->izone[izone].shift0[dim];
599 shift1[dim] = zones->izone[izone].shift1[dim];
600 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
602 /* A conservative approach, this can be optimized */
609 int dd_natoms_vsite(gmx_domdec_t *dd)
611 return dd->comm->nat[ddnatVSITE];
614 void dd_get_constraint_range(gmx_domdec_t *dd,int *at_start,int *at_end)
616 *at_start = dd->comm->nat[ddnatCON-1];
617 *at_end = dd->comm->nat[ddnatCON];
620 void dd_move_x(gmx_domdec_t *dd,matrix box,rvec x[])
622 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
624 gmx_domdec_comm_t *comm;
625 gmx_domdec_comm_dim_t *cd;
626 gmx_domdec_ind_t *ind;
627 rvec shift={0,0,0},*buf,*rbuf;
628 gmx_bool bPBC,bScrew;
632 cgindex = dd->cgindex;
637 nat_tot = dd->nat_home;
638 for(d=0; d<dd->ndim; d++)
640 bPBC = (dd->ci[dd->dim[d]] == 0);
641 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
644 copy_rvec(box[dd->dim[d]],shift);
647 for(p=0; p<cd->np; p++)
654 for(i=0; i<ind->nsend[nzone]; i++)
656 at0 = cgindex[index[i]];
657 at1 = cgindex[index[i]+1];
658 for(j=at0; j<at1; j++)
660 copy_rvec(x[j],buf[n]);
667 for(i=0; i<ind->nsend[nzone]; i++)
669 at0 = cgindex[index[i]];
670 at1 = cgindex[index[i]+1];
671 for(j=at0; j<at1; j++)
673 /* We need to shift the coordinates */
674 rvec_add(x[j],shift,buf[n]);
681 for(i=0; i<ind->nsend[nzone]; i++)
683 at0 = cgindex[index[i]];
684 at1 = cgindex[index[i]+1];
685 for(j=at0; j<at1; j++)
688 buf[n][XX] = x[j][XX] + shift[XX];
690 * This operation requires a special shift force
691 * treatment, which is performed in calc_vir.
693 buf[n][YY] = box[YY][YY] - x[j][YY];
694 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
706 rbuf = comm->vbuf2.v;
708 /* Send and receive the coordinates */
709 dd_sendrecv_rvec(dd, d, dddirBackward,
710 buf, ind->nsend[nzone+1],
711 rbuf, ind->nrecv[nzone+1]);
715 for(zone=0; zone<nzone; zone++)
717 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
719 copy_rvec(rbuf[j],x[i]);
724 nat_tot += ind->nrecv[nzone+1];
730 void dd_move_f(gmx_domdec_t *dd,rvec f[],rvec *fshift)
732 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
734 gmx_domdec_comm_t *comm;
735 gmx_domdec_comm_dim_t *cd;
736 gmx_domdec_ind_t *ind;
740 gmx_bool bPBC,bScrew;
744 cgindex = dd->cgindex;
749 nzone = comm->zones.n/2;
750 nat_tot = dd->nat_tot;
751 for(d=dd->ndim-1; d>=0; d--)
753 bPBC = (dd->ci[dd->dim[d]] == 0);
754 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
755 if (fshift == NULL && !bScrew)
759 /* Determine which shift vector we need */
765 for(p=cd->np-1; p>=0; p--) {
767 nat_tot -= ind->nrecv[nzone+1];
774 sbuf = comm->vbuf2.v;
776 for(zone=0; zone<nzone; zone++)
778 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
780 copy_rvec(f[i],sbuf[j]);
785 /* Communicate the forces */
786 dd_sendrecv_rvec(dd, d, dddirForward,
787 sbuf, ind->nrecv[nzone+1],
788 buf, ind->nsend[nzone+1]);
790 /* Add the received forces */
794 for(i=0; i<ind->nsend[nzone]; i++)
796 at0 = cgindex[index[i]];
797 at1 = cgindex[index[i]+1];
798 for(j=at0; j<at1; j++)
800 rvec_inc(f[j],buf[n]);
807 for(i=0; i<ind->nsend[nzone]; i++)
809 at0 = cgindex[index[i]];
810 at1 = cgindex[index[i]+1];
811 for(j=at0; j<at1; j++)
813 rvec_inc(f[j],buf[n]);
814 /* Add this force to the shift force */
815 rvec_inc(fshift[is],buf[n]);
822 for(i=0; i<ind->nsend[nzone]; i++)
824 at0 = cgindex[index[i]];
825 at1 = cgindex[index[i]+1];
826 for(j=at0; j<at1; j++)
828 /* Rotate the force */
829 f[j][XX] += buf[n][XX];
830 f[j][YY] -= buf[n][YY];
831 f[j][ZZ] -= buf[n][ZZ];
834 /* Add this force to the shift force */
835 rvec_inc(fshift[is],buf[n]);
846 void dd_atom_spread_real(gmx_domdec_t *dd,real v[])
848 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
850 gmx_domdec_comm_t *comm;
851 gmx_domdec_comm_dim_t *cd;
852 gmx_domdec_ind_t *ind;
857 cgindex = dd->cgindex;
859 buf = &comm->vbuf.v[0][0];
862 nat_tot = dd->nat_home;
863 for(d=0; d<dd->ndim; d++)
866 for(p=0; p<cd->np; p++)
871 for(i=0; i<ind->nsend[nzone]; i++)
873 at0 = cgindex[index[i]];
874 at1 = cgindex[index[i]+1];
875 for(j=at0; j<at1; j++)
888 rbuf = &comm->vbuf2.v[0][0];
890 /* Send and receive the coordinates */
891 dd_sendrecv_real(dd, d, dddirBackward,
892 buf, ind->nsend[nzone+1],
893 rbuf, ind->nrecv[nzone+1]);
897 for(zone=0; zone<nzone; zone++)
899 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
906 nat_tot += ind->nrecv[nzone+1];
912 void dd_atom_sum_real(gmx_domdec_t *dd,real v[])
914 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
916 gmx_domdec_comm_t *comm;
917 gmx_domdec_comm_dim_t *cd;
918 gmx_domdec_ind_t *ind;
923 cgindex = dd->cgindex;
925 buf = &comm->vbuf.v[0][0];
928 nzone = comm->zones.n/2;
929 nat_tot = dd->nat_tot;
930 for(d=dd->ndim-1; d>=0; d--)
933 for(p=cd->np-1; p>=0; p--) {
935 nat_tot -= ind->nrecv[nzone+1];
942 sbuf = &comm->vbuf2.v[0][0];
944 for(zone=0; zone<nzone; zone++)
946 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
953 /* Communicate the forces */
954 dd_sendrecv_real(dd, d, dddirForward,
955 sbuf, ind->nrecv[nzone+1],
956 buf, ind->nsend[nzone+1]);
958 /* Add the received forces */
960 for(i=0; i<ind->nsend[nzone]; i++)
962 at0 = cgindex[index[i]];
963 at1 = cgindex[index[i]+1];
964 for(j=at0; j<at1; j++)
975 static void print_ddzone(FILE *fp,int d,int i,int j,gmx_ddzone_t *zone)
977 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",
979 zone->min0,zone->max1,
980 zone->mch0,zone->mch0,
981 zone->p1_0,zone->p1_1);
985 #define DDZONECOMM_MAXZONE 5
986 #define DDZONECOMM_BUFSIZE 3
988 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
989 int ddimind,int direction,
990 gmx_ddzone_t *buf_s,int n_s,
991 gmx_ddzone_t *buf_r,int n_r)
993 #define ZBS DDZONECOMM_BUFSIZE
994 rvec vbuf_s[DDZONECOMM_MAXZONE*ZBS];
995 rvec vbuf_r[DDZONECOMM_MAXZONE*ZBS];
1000 vbuf_s[i*ZBS ][0] = buf_s[i].min0;
1001 vbuf_s[i*ZBS ][1] = buf_s[i].max1;
1002 vbuf_s[i*ZBS ][2] = buf_s[i].min1;
1003 vbuf_s[i*ZBS+1][0] = buf_s[i].mch0;
1004 vbuf_s[i*ZBS+1][1] = buf_s[i].mch1;
1005 vbuf_s[i*ZBS+1][2] = 0;
1006 vbuf_s[i*ZBS+2][0] = buf_s[i].p1_0;
1007 vbuf_s[i*ZBS+2][1] = buf_s[i].p1_1;
1008 vbuf_s[i*ZBS+2][2] = 0;
1011 dd_sendrecv_rvec(dd, ddimind, direction,
1015 for(i=0; i<n_r; i++)
1017 buf_r[i].min0 = vbuf_r[i*ZBS ][0];
1018 buf_r[i].max1 = vbuf_r[i*ZBS ][1];
1019 buf_r[i].min1 = vbuf_r[i*ZBS ][2];
1020 buf_r[i].mch0 = vbuf_r[i*ZBS+1][0];
1021 buf_r[i].mch1 = vbuf_r[i*ZBS+1][1];
1022 buf_r[i].p1_0 = vbuf_r[i*ZBS+2][0];
1023 buf_r[i].p1_1 = vbuf_r[i*ZBS+2][1];
1029 static void dd_move_cellx(gmx_domdec_t *dd,gmx_ddbox_t *ddbox,
1030 rvec cell_ns_x0,rvec cell_ns_x1)
1032 int d,d1,dim,dim1,pos,buf_size,i,j,k,p,npulse,npulse_min;
1034 gmx_ddzone_t buf_s[DDZONECOMM_MAXZONE];
1035 gmx_ddzone_t buf_r[DDZONECOMM_MAXZONE];
1036 gmx_ddzone_t buf_e[DDZONECOMM_MAXZONE];
1037 rvec extr_s[2],extr_r[2];
1039 real dist_d,c=0,det;
1040 gmx_domdec_comm_t *comm;
1045 for(d=1; d<dd->ndim; d++)
1048 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
1049 zp->min0 = cell_ns_x0[dim];
1050 zp->max1 = cell_ns_x1[dim];
1051 zp->min1 = cell_ns_x1[dim];
1052 zp->mch0 = cell_ns_x0[dim];
1053 zp->mch1 = cell_ns_x1[dim];
1054 zp->p1_0 = cell_ns_x0[dim];
1055 zp->p1_1 = cell_ns_x1[dim];
1058 for(d=dd->ndim-2; d>=0; d--)
1061 bPBC = (dim < ddbox->npbcdim);
1063 /* Use an rvec to store two reals */
1064 extr_s[d][0] = comm->cell_f0[d+1];
1065 extr_s[d][1] = comm->cell_f1[d+1];
1066 extr_s[d][2] = comm->cell_f1[d+1];
1069 /* Store the extremes in the backward sending buffer,
1070 * so the get updated separately from the forward communication.
1072 for(d1=d; d1<dd->ndim-1; d1++)
1074 /* We invert the order to be able to use the same loop for buf_e */
1075 buf_s[pos].min0 = extr_s[d1][1];
1076 buf_s[pos].max1 = extr_s[d1][0];
1077 buf_s[pos].min1 = extr_s[d1][2];
1078 buf_s[pos].mch0 = 0;
1079 buf_s[pos].mch1 = 0;
1080 /* Store the cell corner of the dimension we communicate along */
1081 buf_s[pos].p1_0 = comm->cell_x0[dim];
1082 buf_s[pos].p1_1 = 0;
1086 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1089 if (dd->ndim == 3 && d == 0)
1091 buf_s[pos] = comm->zone_d2[0][1];
1093 buf_s[pos] = comm->zone_d1[0];
1097 /* We only need to communicate the extremes
1098 * in the forward direction
1100 npulse = comm->cd[d].np;
1103 /* Take the minimum to avoid double communication */
1104 npulse_min = min(npulse,dd->nc[dim]-1-npulse);
1108 /* Without PBC we should really not communicate over
1109 * the boundaries, but implementing that complicates
1110 * the communication setup and therefore we simply
1111 * do all communication, but ignore some data.
1113 npulse_min = npulse;
1115 for(p=0; p<npulse_min; p++)
1117 /* Communicate the extremes forward */
1118 bUse = (bPBC || dd->ci[dim] > 0);
1120 dd_sendrecv_rvec(dd, d, dddirForward,
1121 extr_s+d, dd->ndim-d-1,
1122 extr_r+d, dd->ndim-d-1);
1126 for(d1=d; d1<dd->ndim-1; d1++)
1128 extr_s[d1][0] = max(extr_s[d1][0],extr_r[d1][0]);
1129 extr_s[d1][1] = min(extr_s[d1][1],extr_r[d1][1]);
1130 extr_s[d1][2] = min(extr_s[d1][2],extr_r[d1][2]);
1136 for(p=0; p<npulse; p++)
1138 /* Communicate all the zone information backward */
1139 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1141 dd_sendrecv_ddzone(dd, d, dddirBackward,
1148 for(d1=d+1; d1<dd->ndim; d1++)
1150 /* Determine the decrease of maximum required
1151 * communication height along d1 due to the distance along d,
1152 * this avoids a lot of useless atom communication.
1154 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1156 if (ddbox->tric_dir[dim])
1158 /* c is the off-diagonal coupling between the cell planes
1159 * along directions d and d1.
1161 c = ddbox->v[dim][dd->dim[d1]][dim];
1167 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1170 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1174 /* A negative value signals out of range */
1180 /* Accumulate the extremes over all pulses */
1181 for(i=0; i<buf_size; i++)
1185 buf_e[i] = buf_r[i];
1191 buf_e[i].min0 = min(buf_e[i].min0,buf_r[i].min0);
1192 buf_e[i].max1 = max(buf_e[i].max1,buf_r[i].max1);
1193 buf_e[i].min1 = min(buf_e[i].min1,buf_r[i].min1);
1196 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1204 if (bUse && dh[d1] >= 0)
1206 buf_e[i].mch0 = max(buf_e[i].mch0,buf_r[i].mch0-dh[d1]);
1207 buf_e[i].mch1 = max(buf_e[i].mch1,buf_r[i].mch1-dh[d1]);
1210 /* Copy the received buffer to the send buffer,
1211 * to pass the data through with the next pulse.
1213 buf_s[i] = buf_r[i];
1215 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1216 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1218 /* Store the extremes */
1221 for(d1=d; d1<dd->ndim-1; d1++)
1223 extr_s[d1][1] = min(extr_s[d1][1],buf_e[pos].min0);
1224 extr_s[d1][0] = max(extr_s[d1][0],buf_e[pos].max1);
1225 extr_s[d1][2] = min(extr_s[d1][2],buf_e[pos].min1);
1229 if (d == 1 || (d == 0 && dd->ndim == 3))
1233 comm->zone_d2[1-d][i] = buf_e[pos];
1239 comm->zone_d1[1] = buf_e[pos];
1253 print_ddzone(debug,1,i,0,&comm->zone_d1[i]);
1255 cell_ns_x0[dim] = min(cell_ns_x0[dim],comm->zone_d1[i].min0);
1256 cell_ns_x1[dim] = max(cell_ns_x1[dim],comm->zone_d1[i].max1);
1268 print_ddzone(debug,2,i,j,&comm->zone_d2[i][j]);
1270 cell_ns_x0[dim] = min(cell_ns_x0[dim],comm->zone_d2[i][j].min0);
1271 cell_ns_x1[dim] = max(cell_ns_x1[dim],comm->zone_d2[i][j].max1);
1275 for(d=1; d<dd->ndim; d++)
1277 comm->cell_f_max0[d] = extr_s[d-1][0];
1278 comm->cell_f_min1[d] = extr_s[d-1][1];
1281 fprintf(debug,"Cell fraction d %d, max0 %f, min1 %f\n",
1282 d,comm->cell_f_max0[d],comm->cell_f_min1[d]);
1287 static void dd_collect_cg(gmx_domdec_t *dd,
1288 t_state *state_local)
1290 gmx_domdec_master_t *ma=NULL;
1291 int buf2[2],*ibuf,i,ncg_home=0,*cg=NULL,nat_home=0;
1294 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1296 /* The master has the correct distribution */
1300 if (state_local->ddp_count == dd->ddp_count)
1302 ncg_home = dd->ncg_home;
1304 nat_home = dd->nat_home;
1306 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1308 cgs_gl = &dd->comm->cgs_gl;
1310 ncg_home = state_local->ncg_gl;
1311 cg = state_local->cg_gl;
1313 for(i=0; i<ncg_home; i++)
1315 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1320 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1323 buf2[0] = dd->ncg_home;
1324 buf2[1] = dd->nat_home;
1334 /* Collect the charge group and atom counts on the master */
1335 dd_gather(dd,2*sizeof(int),buf2,ibuf);
1340 for(i=0; i<dd->nnodes; i++)
1342 ma->ncg[i] = ma->ibuf[2*i];
1343 ma->nat[i] = ma->ibuf[2*i+1];
1344 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1347 /* Make byte counts and indices */
1348 for(i=0; i<dd->nnodes; i++)
1350 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1351 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1355 fprintf(debug,"Initial charge group distribution: ");
1356 for(i=0; i<dd->nnodes; i++)
1357 fprintf(debug," %d",ma->ncg[i]);
1358 fprintf(debug,"\n");
1362 /* Collect the charge group indices on the master */
1364 dd->ncg_home*sizeof(int),dd->index_gl,
1365 DDMASTER(dd) ? ma->ibuf : NULL,
1366 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1367 DDMASTER(dd) ? ma->cg : NULL);
1369 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1372 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1375 gmx_domdec_master_t *ma;
1376 int n,i,c,a,nalloc=0;
1385 MPI_Send(lv,dd->nat_home*sizeof(rvec),MPI_BYTE,DDMASTERRANK(dd),
1386 dd->rank,dd->mpi_comm_all);
1389 /* Copy the master coordinates to the global array */
1390 cgs_gl = &dd->comm->cgs_gl;
1392 n = DDMASTERRANK(dd);
1394 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1396 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1398 copy_rvec(lv[a++],v[c]);
1402 for(n=0; n<dd->nnodes; n++)
1406 if (ma->nat[n] > nalloc)
1408 nalloc = over_alloc_dd(ma->nat[n]);
1412 MPI_Recv(buf,ma->nat[n]*sizeof(rvec),MPI_BYTE,DDRANK(dd,n),
1413 n,dd->mpi_comm_all,MPI_STATUS_IGNORE);
1416 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1418 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1420 copy_rvec(buf[a++],v[c]);
1429 static void get_commbuffer_counts(gmx_domdec_t *dd,
1430 int **counts,int **disps)
1432 gmx_domdec_master_t *ma;
1437 /* Make the rvec count and displacment arrays */
1439 *disps = ma->ibuf + dd->nnodes;
1440 for(n=0; n<dd->nnodes; n++)
1442 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1443 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1447 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1450 gmx_domdec_master_t *ma;
1451 int *rcounts=NULL,*disps=NULL;
1460 get_commbuffer_counts(dd,&rcounts,&disps);
1465 dd_gatherv(dd,dd->nat_home*sizeof(rvec),lv,rcounts,disps,buf);
1469 cgs_gl = &dd->comm->cgs_gl;
1472 for(n=0; n<dd->nnodes; n++)
1474 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1476 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1478 copy_rvec(buf[a++],v[c]);
1485 void dd_collect_vec(gmx_domdec_t *dd,
1486 t_state *state_local,rvec *lv,rvec *v)
1488 gmx_domdec_master_t *ma;
1489 int n,i,c,a,nalloc=0;
1492 dd_collect_cg(dd,state_local);
1494 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1496 dd_collect_vec_sendrecv(dd,lv,v);
1500 dd_collect_vec_gatherv(dd,lv,v);
1505 void dd_collect_state(gmx_domdec_t *dd,
1506 t_state *state_local,t_state *state)
1510 nh = state->nhchainlength;
1514 for (i=0;i<efptNR;i++) {
1515 state->lambda[i] = state_local->lambda[i];
1517 state->fep_state = state_local->fep_state;
1518 state->veta = state_local->veta;
1519 state->vol0 = state_local->vol0;
1520 copy_mat(state_local->box,state->box);
1521 copy_mat(state_local->boxv,state->boxv);
1522 copy_mat(state_local->svir_prev,state->svir_prev);
1523 copy_mat(state_local->fvir_prev,state->fvir_prev);
1524 copy_mat(state_local->pres_prev,state->pres_prev);
1527 for(i=0; i<state_local->ngtc; i++)
1529 for(j=0; j<nh; j++) {
1530 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1531 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1533 state->therm_integral[i] = state_local->therm_integral[i];
1535 for(i=0; i<state_local->nnhpres; i++)
1537 for(j=0; j<nh; j++) {
1538 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1539 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1543 for(est=0; est<estNR; est++)
1545 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1549 dd_collect_vec(dd,state_local,state_local->x,state->x);
1552 dd_collect_vec(dd,state_local,state_local->v,state->v);
1555 dd_collect_vec(dd,state_local,state_local->sd_X,state->sd_X);
1558 dd_collect_vec(dd,state_local,state_local->cg_p,state->cg_p);
1561 if (state->nrngi == 1)
1565 for(i=0; i<state_local->nrng; i++)
1567 state->ld_rng[i] = state_local->ld_rng[i];
1573 dd_gather(dd,state_local->nrng*sizeof(state->ld_rng[0]),
1574 state_local->ld_rng,state->ld_rng);
1578 if (state->nrngi == 1)
1582 state->ld_rngi[0] = state_local->ld_rngi[0];
1587 dd_gather(dd,sizeof(state->ld_rngi[0]),
1588 state_local->ld_rngi,state->ld_rngi);
1591 case estDISRE_INITF:
1592 case estDISRE_RM3TAV:
1593 case estORIRE_INITF:
1597 gmx_incons("Unknown state entry encountered in dd_collect_state");
1603 static void dd_realloc_state(t_state *state,rvec **f,int nalloc)
1609 fprintf(debug,"Reallocating state: currently %d, required %d, allocating %d\n",state->nalloc,nalloc,over_alloc_dd(nalloc));
1612 state->nalloc = over_alloc_dd(nalloc);
1614 for(est=0; est<estNR; est++)
1616 if (EST_DISTR(est) && (state->flags & (1<<est)))
1620 srenew(state->x,state->nalloc);
1623 srenew(state->v,state->nalloc);
1626 srenew(state->sd_X,state->nalloc);
1629 srenew(state->cg_p,state->nalloc);
1633 case estDISRE_INITF:
1634 case estDISRE_RM3TAV:
1635 case estORIRE_INITF:
1637 /* No reallocation required */
1640 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1647 srenew(*f,state->nalloc);
1651 static void dd_check_alloc_ncg(t_forcerec *fr,t_state *state,rvec **f,
1654 if (nalloc > fr->cg_nalloc)
1658 fprintf(debug,"Reallocating forcerec: currently %d, required %d, allocating %d\n",fr->cg_nalloc,nalloc,over_alloc_dd(nalloc));
1660 fr->cg_nalloc = over_alloc_dd(nalloc);
1661 srenew(fr->cginfo,fr->cg_nalloc);
1662 if (fr->cutoff_scheme == ecutsGROUP)
1664 srenew(fr->cg_cm,fr->cg_nalloc);
1667 if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
1669 /* We don't use charge groups, we use x in state to set up
1670 * the atom communication.
1672 dd_realloc_state(state,f,nalloc);
1676 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd,t_block *cgs,
1679 gmx_domdec_master_t *ma;
1680 int n,i,c,a,nalloc=0;
1687 for(n=0; n<dd->nnodes; n++)
1691 if (ma->nat[n] > nalloc)
1693 nalloc = over_alloc_dd(ma->nat[n]);
1696 /* Use lv as a temporary buffer */
1698 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1700 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1702 copy_rvec(v[c],buf[a++]);
1705 if (a != ma->nat[n])
1707 gmx_fatal(FARGS,"Internal error a (%d) != nat (%d)",
1712 MPI_Send(buf,ma->nat[n]*sizeof(rvec),MPI_BYTE,
1713 DDRANK(dd,n),n,dd->mpi_comm_all);
1718 n = DDMASTERRANK(dd);
1720 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1722 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1724 copy_rvec(v[c],lv[a++]);
1731 MPI_Recv(lv,dd->nat_home*sizeof(rvec),MPI_BYTE,DDMASTERRANK(dd),
1732 MPI_ANY_TAG,dd->mpi_comm_all,MPI_STATUS_IGNORE);
1737 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd,t_block *cgs,
1740 gmx_domdec_master_t *ma;
1741 int *scounts=NULL,*disps=NULL;
1742 int n,i,c,a,nalloc=0;
1749 get_commbuffer_counts(dd,&scounts,&disps);
1753 for(n=0; n<dd->nnodes; n++)
1755 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1757 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1759 copy_rvec(v[c],buf[a++]);
1765 dd_scatterv(dd,scounts,disps,buf,dd->nat_home*sizeof(rvec),lv);
1768 static void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv)
1770 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1772 dd_distribute_vec_sendrecv(dd,cgs,v,lv);
1776 dd_distribute_vec_scatterv(dd,cgs,v,lv);
1780 static void dd_distribute_state(gmx_domdec_t *dd,t_block *cgs,
1781 t_state *state,t_state *state_local,
1786 nh = state->nhchainlength;
1790 for(i=0;i<efptNR;i++)
1792 state_local->lambda[i] = state->lambda[i];
1794 state_local->fep_state = state->fep_state;
1795 state_local->veta = state->veta;
1796 state_local->vol0 = state->vol0;
1797 copy_mat(state->box,state_local->box);
1798 copy_mat(state->box_rel,state_local->box_rel);
1799 copy_mat(state->boxv,state_local->boxv);
1800 copy_mat(state->svir_prev,state_local->svir_prev);
1801 copy_mat(state->fvir_prev,state_local->fvir_prev);
1802 for(i=0; i<state_local->ngtc; i++)
1804 for(j=0; j<nh; j++) {
1805 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1806 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1808 state_local->therm_integral[i] = state->therm_integral[i];
1810 for(i=0; i<state_local->nnhpres; i++)
1812 for(j=0; j<nh; j++) {
1813 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1814 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1818 dd_bcast(dd,((efptNR)*sizeof(real)),state_local->lambda);
1819 dd_bcast(dd,sizeof(int),&state_local->fep_state);
1820 dd_bcast(dd,sizeof(real),&state_local->veta);
1821 dd_bcast(dd,sizeof(real),&state_local->vol0);
1822 dd_bcast(dd,sizeof(state_local->box),state_local->box);
1823 dd_bcast(dd,sizeof(state_local->box_rel),state_local->box_rel);
1824 dd_bcast(dd,sizeof(state_local->boxv),state_local->boxv);
1825 dd_bcast(dd,sizeof(state_local->svir_prev),state_local->svir_prev);
1826 dd_bcast(dd,sizeof(state_local->fvir_prev),state_local->fvir_prev);
1827 dd_bcast(dd,((state_local->ngtc*nh)*sizeof(double)),state_local->nosehoover_xi);
1828 dd_bcast(dd,((state_local->ngtc*nh)*sizeof(double)),state_local->nosehoover_vxi);
1829 dd_bcast(dd,state_local->ngtc*sizeof(double),state_local->therm_integral);
1830 dd_bcast(dd,((state_local->nnhpres*nh)*sizeof(double)),state_local->nhpres_xi);
1831 dd_bcast(dd,((state_local->nnhpres*nh)*sizeof(double)),state_local->nhpres_vxi);
1833 if (dd->nat_home > state_local->nalloc)
1835 dd_realloc_state(state_local,f,dd->nat_home);
1837 for(i=0; i<estNR; i++)
1839 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1843 dd_distribute_vec(dd,cgs,state->x,state_local->x);
1846 dd_distribute_vec(dd,cgs,state->v,state_local->v);
1849 dd_distribute_vec(dd,cgs,state->sd_X,state_local->sd_X);
1852 dd_distribute_vec(dd,cgs,state->cg_p,state_local->cg_p);
1855 if (state->nrngi == 1)
1858 state_local->nrng*sizeof(state_local->ld_rng[0]),
1859 state->ld_rng,state_local->ld_rng);
1864 state_local->nrng*sizeof(state_local->ld_rng[0]),
1865 state->ld_rng,state_local->ld_rng);
1869 if (state->nrngi == 1)
1871 dd_bcastc(dd,sizeof(state_local->ld_rngi[0]),
1872 state->ld_rngi,state_local->ld_rngi);
1876 dd_scatter(dd,sizeof(state_local->ld_rngi[0]),
1877 state->ld_rngi,state_local->ld_rngi);
1880 case estDISRE_INITF:
1881 case estDISRE_RM3TAV:
1882 case estORIRE_INITF:
1884 /* Not implemented yet */
1887 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1893 static char dim2char(int dim)
1899 case XX: c = 'X'; break;
1900 case YY: c = 'Y'; break;
1901 case ZZ: c = 'Z'; break;
1902 default: gmx_fatal(FARGS,"Unknown dim %d",dim);
1908 static void write_dd_grid_pdb(const char *fn,gmx_large_int_t step,
1909 gmx_domdec_t *dd,matrix box,gmx_ddbox_t *ddbox)
1911 rvec grid_s[2],*grid_r=NULL,cx,r;
1912 char fname[STRLEN],format[STRLEN],buf[22];
1918 copy_rvec(dd->comm->cell_x0,grid_s[0]);
1919 copy_rvec(dd->comm->cell_x1,grid_s[1]);
1923 snew(grid_r,2*dd->nnodes);
1926 dd_gather(dd,2*sizeof(rvec),grid_s[0],DDMASTER(dd) ? grid_r[0] : NULL);
1930 for(d=0; d<DIM; d++)
1932 for(i=0; i<DIM; i++)
1940 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1942 tric[d][i] = box[i][d]/box[i][i];
1951 sprintf(fname,"%s_%s.pdb",fn,gmx_step_str(step,buf));
1952 sprintf(format,"%s%s\n",get_pdbformat(),"%6.2f%6.2f");
1953 out = gmx_fio_fopen(fname,"w");
1954 gmx_write_pdb_box(out,dd->bScrewPBC ? epbcSCREW : epbcXYZ,box);
1956 for(i=0; i<dd->nnodes; i++)
1958 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1959 for(d=0; d<DIM; d++)
1961 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1969 cx[XX] = grid_r[i*2+x][XX];
1970 cx[YY] = grid_r[i*2+y][YY];
1971 cx[ZZ] = grid_r[i*2+z][ZZ];
1973 fprintf(out,format,"ATOM",a++,"CA","GLY",' ',1+i,
1974 10*r[XX],10*r[YY],10*r[ZZ],1.0,vol);
1978 for(d=0; d<DIM; d++)
1984 case 0: y = 1 + i*8 + 2*x; break;
1985 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
1986 case 2: y = 1 + i*8 + x; break;
1988 fprintf(out,"%6s%5d%5d\n","CONECT",y,y+(1<<d));
1992 gmx_fio_fclose(out);
1997 void write_dd_pdb(const char *fn,gmx_large_int_t step,const char *title,
1998 gmx_mtop_t *mtop,t_commrec *cr,
1999 int natoms,rvec x[],matrix box)
2001 char fname[STRLEN],format[STRLEN],format4[STRLEN],buf[22];
2004 char *atomname,*resname;
2011 natoms = dd->comm->nat[ddnatVSITE];
2014 sprintf(fname,"%s_%s_n%d.pdb",fn,gmx_step_str(step,buf),cr->sim_nodeid);
2016 sprintf(format,"%s%s\n",get_pdbformat(),"%6.2f%6.2f");
2017 sprintf(format4,"%s%s\n",get_pdbformat4(),"%6.2f%6.2f");
2019 out = gmx_fio_fopen(fname,"w");
2021 fprintf(out,"TITLE %s\n",title);
2022 gmx_write_pdb_box(out,dd->bScrewPBC ? epbcSCREW : epbcXYZ,box);
2023 for(i=0; i<natoms; i++)
2025 ii = dd->gatindex[i];
2026 gmx_mtop_atominfo_global(mtop,ii,&atomname,&resnr,&resname);
2027 if (i < dd->comm->nat[ddnatZONE])
2030 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
2036 else if (i < dd->comm->nat[ddnatVSITE])
2038 b = dd->comm->zones.n;
2042 b = dd->comm->zones.n + 1;
2044 fprintf(out,strlen(atomname)<4 ? format : format4,
2045 "ATOM",(ii+1)%100000,
2046 atomname,resname,' ',resnr%10000,' ',
2047 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],1.0,b);
2049 fprintf(out,"TER\n");
2051 gmx_fio_fclose(out);
2054 real dd_cutoff_mbody(gmx_domdec_t *dd)
2056 gmx_domdec_comm_t *comm;
2063 if (comm->bInterCGBondeds)
2065 if (comm->cutoff_mbody > 0)
2067 r = comm->cutoff_mbody;
2071 /* cutoff_mbody=0 means we do not have DLB */
2072 r = comm->cellsize_min[dd->dim[0]];
2073 for(di=1; di<dd->ndim; di++)
2075 r = min(r,comm->cellsize_min[dd->dim[di]]);
2077 if (comm->bBondComm)
2079 r = max(r,comm->cutoff_mbody);
2083 r = min(r,comm->cutoff);
2091 real dd_cutoff_twobody(gmx_domdec_t *dd)
2095 r_mb = dd_cutoff_mbody(dd);
2097 return max(dd->comm->cutoff,r_mb);
2101 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd,ivec coord,ivec coord_pme)
2105 nc = dd->nc[dd->comm->cartpmedim];
2106 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2107 copy_ivec(coord,coord_pme);
2108 coord_pme[dd->comm->cartpmedim] =
2109 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2112 static int low_ddindex2pmeindex(int ndd,int npme,int ddindex)
2114 /* Here we assign a PME node to communicate with this DD node
2115 * by assuming that the major index of both is x.
2116 * We add cr->npmenodes/2 to obtain an even distribution.
2118 return (ddindex*npme + npme/2)/ndd;
2121 static int ddindex2pmeindex(const gmx_domdec_t *dd,int ddindex)
2123 return low_ddindex2pmeindex(dd->nnodes,dd->comm->npmenodes,ddindex);
2126 static int cr_ddindex2pmeindex(const t_commrec *cr,int ddindex)
2128 return low_ddindex2pmeindex(cr->dd->nnodes,cr->npmenodes,ddindex);
2131 static int *dd_pmenodes(t_commrec *cr)
2136 snew(pmenodes,cr->npmenodes);
2138 for(i=0; i<cr->dd->nnodes; i++) {
2139 p0 = cr_ddindex2pmeindex(cr,i);
2140 p1 = cr_ddindex2pmeindex(cr,i+1);
2141 if (i+1 == cr->dd->nnodes || p1 > p0) {
2143 fprintf(debug,"pmenode[%d] = %d\n",n,i+1+n);
2144 pmenodes[n] = i + 1 + n;
2152 static int gmx_ddcoord2pmeindex(t_commrec *cr,int x,int y,int z)
2155 ivec coords,coords_pme,nc;
2160 if (dd->comm->bCartesian) {
2161 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2162 dd_coords2pmecoords(dd,coords,coords_pme);
2163 copy_ivec(dd->ntot,nc);
2164 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2165 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2167 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2169 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2175 slab = ddindex2pmeindex(dd,dd_index(dd->nc,coords));
2180 static int ddcoord2simnodeid(t_commrec *cr,int x,int y,int z)
2182 gmx_domdec_comm_t *comm;
2184 int ddindex,nodeid=-1;
2186 comm = cr->dd->comm;
2191 if (comm->bCartesianPP_PME)
2194 MPI_Cart_rank(cr->mpi_comm_mysim,coords,&nodeid);
2199 ddindex = dd_index(cr->dd->nc,coords);
2200 if (comm->bCartesianPP)
2202 nodeid = comm->ddindex2simnodeid[ddindex];
2208 nodeid = ddindex + gmx_ddcoord2pmeindex(cr,x,y,z);
2220 static int dd_simnode2pmenode(t_commrec *cr,int sim_nodeid)
2223 gmx_domdec_comm_t *comm;
2224 ivec coord,coord_pme;
2231 /* This assumes a uniform x domain decomposition grid cell size */
2232 if (comm->bCartesianPP_PME)
2235 MPI_Cart_coords(cr->mpi_comm_mysim,sim_nodeid,DIM,coord);
2236 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2238 /* This is a PP node */
2239 dd_cart_coord2pmecoord(dd,coord,coord_pme);
2240 MPI_Cart_rank(cr->mpi_comm_mysim,coord_pme,&pmenode);
2244 else if (comm->bCartesianPP)
2246 if (sim_nodeid < dd->nnodes)
2248 pmenode = dd->nnodes + ddindex2pmeindex(dd,sim_nodeid);
2253 /* This assumes DD cells with identical x coordinates
2254 * are numbered sequentially.
2256 if (dd->comm->pmenodes == NULL)
2258 if (sim_nodeid < dd->nnodes)
2260 /* The DD index equals the nodeid */
2261 pmenode = dd->nnodes + ddindex2pmeindex(dd,sim_nodeid);
2267 while (sim_nodeid > dd->comm->pmenodes[i])
2271 if (sim_nodeid < dd->comm->pmenodes[i])
2273 pmenode = dd->comm->pmenodes[i];
2281 gmx_bool gmx_pmeonlynode(t_commrec *cr,int sim_nodeid)
2283 gmx_bool bPMEOnlyNode;
2285 if (DOMAINDECOMP(cr))
2287 bPMEOnlyNode = (dd_simnode2pmenode(cr,sim_nodeid) == -1);
2291 bPMEOnlyNode = FALSE;
2294 return bPMEOnlyNode;
2297 void get_pme_ddnodes(t_commrec *cr,int pmenodeid,
2298 int *nmy_ddnodes,int **my_ddnodes,int *node_peer)
2302 ivec coord,coord_pme;
2306 snew(*my_ddnodes,(dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2309 for(x=0; x<dd->nc[XX]; x++)
2311 for(y=0; y<dd->nc[YY]; y++)
2313 for(z=0; z<dd->nc[ZZ]; z++)
2315 if (dd->comm->bCartesianPP_PME)
2320 dd_cart_coord2pmecoord(dd,coord,coord_pme);
2321 if (dd->ci[XX] == coord_pme[XX] &&
2322 dd->ci[YY] == coord_pme[YY] &&
2323 dd->ci[ZZ] == coord_pme[ZZ])
2324 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr,x,y,z);
2328 /* The slab corresponds to the nodeid in the PME group */
2329 if (gmx_ddcoord2pmeindex(cr,x,y,z) == pmenodeid)
2331 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr,x,y,z);
2338 /* The last PP-only node is the peer node */
2339 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2343 fprintf(debug,"Receive coordinates from PP nodes:");
2344 for(x=0; x<*nmy_ddnodes; x++)
2346 fprintf(debug," %d",(*my_ddnodes)[x]);
2348 fprintf(debug,"\n");
2352 static gmx_bool receive_vir_ener(t_commrec *cr)
2354 gmx_domdec_comm_t *comm;
2355 int pmenode,coords[DIM],rank;
2359 if (cr->npmenodes < cr->dd->nnodes)
2361 comm = cr->dd->comm;
2362 if (comm->bCartesianPP_PME)
2364 pmenode = dd_simnode2pmenode(cr,cr->sim_nodeid);
2366 MPI_Cart_coords(cr->mpi_comm_mysim,cr->sim_nodeid,DIM,coords);
2367 coords[comm->cartpmedim]++;
2368 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2370 MPI_Cart_rank(cr->mpi_comm_mysim,coords,&rank);
2371 if (dd_simnode2pmenode(cr,rank) == pmenode)
2373 /* This is not the last PP node for pmenode */
2381 pmenode = dd_simnode2pmenode(cr,cr->sim_nodeid);
2382 if (cr->sim_nodeid+1 < cr->nnodes &&
2383 dd_simnode2pmenode(cr,cr->sim_nodeid+1) == pmenode)
2385 /* This is not the last PP node for pmenode */
2394 static void set_zones_ncg_home(gmx_domdec_t *dd)
2396 gmx_domdec_zones_t *zones;
2399 zones = &dd->comm->zones;
2401 zones->cg_range[0] = 0;
2402 for(i=1; i<zones->n+1; i++)
2404 zones->cg_range[i] = dd->ncg_home;
2408 static void rebuild_cgindex(gmx_domdec_t *dd,
2409 const int *gcgs_index,t_state *state)
2411 int nat,i,*ind,*dd_cg_gl,*cgindex,cg_gl;
2414 dd_cg_gl = dd->index_gl;
2415 cgindex = dd->cgindex;
2418 for(i=0; i<state->ncg_gl; i++)
2422 dd_cg_gl[i] = cg_gl;
2423 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2427 dd->ncg_home = state->ncg_gl;
2430 set_zones_ncg_home(dd);
2433 static int ddcginfo(const cginfo_mb_t *cginfo_mb,int cg)
2435 while (cg >= cginfo_mb->cg_end)
2440 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2443 static void dd_set_cginfo(int *index_gl,int cg0,int cg1,
2444 t_forcerec *fr,char *bLocalCG)
2446 cginfo_mb_t *cginfo_mb;
2452 cginfo_mb = fr->cginfo_mb;
2453 cginfo = fr->cginfo;
2455 for(cg=cg0; cg<cg1; cg++)
2457 cginfo[cg] = ddcginfo(cginfo_mb,index_gl[cg]);
2461 if (bLocalCG != NULL)
2463 for(cg=cg0; cg<cg1; cg++)
2465 bLocalCG[index_gl[cg]] = TRUE;
2470 static void make_dd_indices(gmx_domdec_t *dd,
2471 const int *gcgs_index,int cg_start)
2473 int nzone,zone,zone1,cg0,cg1,cg1_p1,cg,cg_gl,a,a_gl;
2474 int *zone2cg,*zone_ncg1,*index_gl,*gatindex;
2479 bLocalCG = dd->comm->bLocalCG;
2481 if (dd->nat_tot > dd->gatindex_nalloc)
2483 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2484 srenew(dd->gatindex,dd->gatindex_nalloc);
2487 nzone = dd->comm->zones.n;
2488 zone2cg = dd->comm->zones.cg_range;
2489 zone_ncg1 = dd->comm->zone_ncg1;
2490 index_gl = dd->index_gl;
2491 gatindex = dd->gatindex;
2492 bCGs = dd->comm->bCGs;
2494 if (zone2cg[1] != dd->ncg_home)
2496 gmx_incons("dd->ncg_zone is not up to date");
2499 /* Make the local to global and global to local atom index */
2500 a = dd->cgindex[cg_start];
2501 for(zone=0; zone<nzone; zone++)
2509 cg0 = zone2cg[zone];
2511 cg1 = zone2cg[zone+1];
2512 cg1_p1 = cg0 + zone_ncg1[zone];
2514 for(cg=cg0; cg<cg1; cg++)
2519 /* Signal that this cg is from more than one pulse away */
2522 cg_gl = index_gl[cg];
2525 for(a_gl=gcgs_index[cg_gl]; a_gl<gcgs_index[cg_gl+1]; a_gl++)
2528 ga2la_set(dd->ga2la,a_gl,a,zone1);
2534 gatindex[a] = cg_gl;
2535 ga2la_set(dd->ga2la,cg_gl,a,zone1);
2542 static int check_bLocalCG(gmx_domdec_t *dd,int ncg_sys,const char *bLocalCG,
2548 if (bLocalCG == NULL)
2552 for(i=0; i<dd->ncg_tot; i++)
2554 if (!bLocalCG[dd->index_gl[i]])
2557 "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);
2562 for(i=0; i<ncg_sys; i++)
2569 if (ngl != dd->ncg_tot)
2571 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);
2578 static void check_index_consistency(gmx_domdec_t *dd,
2579 int natoms_sys,int ncg_sys,
2582 int nerr,ngl,i,a,cell;
2587 if (dd->comm->DD_debug > 1)
2589 snew(have,natoms_sys);
2590 for(a=0; a<dd->nat_tot; a++)
2592 if (have[dd->gatindex[a]] > 0)
2594 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);
2598 have[dd->gatindex[a]] = a + 1;
2604 snew(have,dd->nat_tot);
2607 for(i=0; i<natoms_sys; i++)
2609 if (ga2la_get(dd->ga2la,i,&a,&cell))
2611 if (a >= dd->nat_tot)
2613 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);
2619 if (dd->gatindex[a] != i)
2621 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);
2628 if (ngl != dd->nat_tot)
2631 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2632 dd->rank,where,ngl,dd->nat_tot);
2634 for(a=0; a<dd->nat_tot; a++)
2639 "DD node %d, %s: local atom %d, global %d has no global index\n",
2640 dd->rank,where,a+1,dd->gatindex[a]+1);
2645 nerr += check_bLocalCG(dd,ncg_sys,dd->comm->bLocalCG,where);
2648 gmx_fatal(FARGS,"DD node %d, %s: %d atom/cg index inconsistencies",
2649 dd->rank,where,nerr);
2653 static void clear_dd_indices(gmx_domdec_t *dd,int cg_start,int a_start)
2660 /* Clear the whole list without searching */
2661 ga2la_clear(dd->ga2la);
2665 for(i=a_start; i<dd->nat_tot; i++)
2667 ga2la_del(dd->ga2la,dd->gatindex[i]);
2671 bLocalCG = dd->comm->bLocalCG;
2674 for(i=cg_start; i<dd->ncg_tot; i++)
2676 bLocalCG[dd->index_gl[i]] = FALSE;
2680 dd_clear_local_vsite_indices(dd);
2682 if (dd->constraints)
2684 dd_clear_local_constraint_indices(dd);
2688 /* This function should be used for moving the domain boudaries during DLB,
2689 * for obtaining the minimum cell size. It checks the initially set limit
2690 * comm->cellsize_min, for bonded and initial non-bonded cut-offs,
2691 * and, possibly, a longer cut-off limit set for PME load balancing.
2693 static real cellsize_min_dlb(gmx_domdec_comm_t *comm,int dim_ind,int dim)
2697 cellsize_min = comm->cellsize_min[dim];
2699 if (!comm->bVacDLBNoLimit && comm->bPMELoadBalDLBLimits)
2701 cellsize_min = max(cellsize_min,
2702 comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
2705 return cellsize_min;
2708 static real grid_jump_limit(gmx_domdec_comm_t *comm,real cutoff,
2711 real grid_jump_limit;
2713 /* The distance between the boundaries of cells at distance
2714 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2715 * and by the fact that cells should not be shifted by more than
2716 * half their size, such that cg's only shift by one cell
2717 * at redecomposition.
2719 grid_jump_limit = comm->cellsize_limit;
2720 if (!comm->bVacDLBNoLimit)
2722 if (comm->bPMELoadBalDLBLimits)
2724 cutoff = max(cutoff,comm->PMELoadBal_max_cutoff);
2726 grid_jump_limit = max(grid_jump_limit,
2727 cutoff/comm->cd[dim_ind].np);
2730 return grid_jump_limit;
2733 static gmx_bool check_grid_jump(gmx_large_int_t step,
2739 gmx_domdec_comm_t *comm;
2748 for(d=1; d<dd->ndim; d++)
2751 limit = grid_jump_limit(comm,cutoff,d);
2752 bfac = ddbox->box_size[dim];
2753 if (ddbox->tric_dir[dim])
2755 bfac *= ddbox->skew_fac[dim];
2757 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2758 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2766 /* This error should never be triggered under normal
2767 * circumstances, but you never know ...
2769 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.",
2770 gmx_step_str(step,buf),
2771 dim2char(dim),dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
2779 static int dd_load_count(gmx_domdec_comm_t *comm)
2781 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2784 static float dd_force_load(gmx_domdec_comm_t *comm)
2791 if (comm->eFlop > 1)
2793 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2798 load = comm->cycl[ddCyclF];
2799 if (comm->cycl_n[ddCyclF] > 1)
2801 /* Subtract the maximum of the last n cycle counts
2802 * to get rid of possible high counts due to other soures,
2803 * for instance system activity, that would otherwise
2804 * affect the dynamic load balancing.
2806 load -= comm->cycl_max[ddCyclF];
2813 static void set_slb_pme_dim_f(gmx_domdec_t *dd,int dim,real **dim_f)
2815 gmx_domdec_comm_t *comm;
2820 snew(*dim_f,dd->nc[dim]+1);
2822 for(i=1; i<dd->nc[dim]; i++)
2824 if (comm->slb_frac[dim])
2826 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2830 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2833 (*dim_f)[dd->nc[dim]] = 1;
2836 static void init_ddpme(gmx_domdec_t *dd,gmx_ddpme_t *ddpme,int dimind)
2838 int pmeindex,slab,nso,i;
2841 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2847 ddpme->dim = dimind;
2849 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2851 ddpme->nslab = (ddpme->dim == 0 ?
2852 dd->comm->npmenodes_x :
2853 dd->comm->npmenodes_y);
2855 if (ddpme->nslab <= 1)
2860 nso = dd->comm->npmenodes/ddpme->nslab;
2861 /* Determine for each PME slab the PP location range for dimension dim */
2862 snew(ddpme->pp_min,ddpme->nslab);
2863 snew(ddpme->pp_max,ddpme->nslab);
2864 for(slab=0; slab<ddpme->nslab; slab++) {
2865 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2866 ddpme->pp_max[slab] = 0;
2868 for(i=0; i<dd->nnodes; i++) {
2869 ddindex2xyz(dd->nc,i,xyz);
2870 /* For y only use our y/z slab.
2871 * This assumes that the PME x grid size matches the DD grid size.
2873 if (dimind == 0 || xyz[XX] == dd->ci[XX]) {
2874 pmeindex = ddindex2pmeindex(dd,i);
2876 slab = pmeindex/nso;
2878 slab = pmeindex % ddpme->nslab;
2880 ddpme->pp_min[slab] = min(ddpme->pp_min[slab],xyz[dimind]);
2881 ddpme->pp_max[slab] = max(ddpme->pp_max[slab],xyz[dimind]);
2885 set_slb_pme_dim_f(dd,ddpme->dim,&ddpme->slb_dim_f);
2888 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2890 if (dd->comm->ddpme[0].dim == XX)
2892 return dd->comm->ddpme[0].maxshift;
2900 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2902 if (dd->comm->ddpme[0].dim == YY)
2904 return dd->comm->ddpme[0].maxshift;
2906 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2908 return dd->comm->ddpme[1].maxshift;
2916 static void set_pme_maxshift(gmx_domdec_t *dd,gmx_ddpme_t *ddpme,
2917 gmx_bool bUniform,gmx_ddbox_t *ddbox,real *cell_f)
2919 gmx_domdec_comm_t *comm;
2922 real range,pme_boundary;
2926 nc = dd->nc[ddpme->dim];
2929 if (!ddpme->dim_match)
2931 /* PP decomposition is not along dim: the worst situation */
2934 else if (ns <= 3 || (bUniform && ns == nc))
2936 /* The optimal situation */
2941 /* We need to check for all pme nodes which nodes they
2942 * could possibly need to communicate with.
2944 xmin = ddpme->pp_min;
2945 xmax = ddpme->pp_max;
2946 /* Allow for atoms to be maximally 2/3 times the cut-off
2947 * out of their DD cell. This is a reasonable balance between
2948 * between performance and support for most charge-group/cut-off
2951 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2952 /* Avoid extra communication when we are exactly at a boundary */
2958 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2959 pme_boundary = (real)s/ns;
2962 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
2964 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
2968 pme_boundary = (real)(s+1)/ns;
2971 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
2973 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
2980 ddpme->maxshift = sh;
2984 fprintf(debug,"PME slab communication range for dim %d is %d\n",
2985 ddpme->dim,ddpme->maxshift);
2989 static void check_box_size(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
2993 for(d=0; d<dd->ndim; d++)
2996 if (dim < ddbox->nboundeddim &&
2997 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
2998 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
3000 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",
3001 dim2char(dim),ddbox->box_size[dim],ddbox->skew_fac[dim],
3002 dd->nc[dim],dd->comm->cellsize_limit);
3007 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd,gmx_ddbox_t *ddbox,
3008 gmx_bool bMaster,ivec npulse)
3010 gmx_domdec_comm_t *comm;
3013 real *cell_x,cell_dx,cellsize;
3017 for(d=0; d<DIM; d++)
3019 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
3021 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
3024 cell_dx = ddbox->box_size[d]/dd->nc[d];
3027 for(j=0; j<dd->nc[d]+1; j++)
3029 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
3034 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
3035 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
3037 cellsize = cell_dx*ddbox->skew_fac[d];
3038 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
3042 cellsize_min[d] = cellsize;
3046 /* Statically load balanced grid */
3047 /* Also when we are not doing a master distribution we determine
3048 * all cell borders in a loop to obtain identical values
3049 * to the master distribution case and to determine npulse.
3053 cell_x = dd->ma->cell_x[d];
3057 snew(cell_x,dd->nc[d]+1);
3059 cell_x[0] = ddbox->box0[d];
3060 for(j=0; j<dd->nc[d]; j++)
3062 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
3063 cell_x[j+1] = cell_x[j] + cell_dx;
3064 cellsize = cell_dx*ddbox->skew_fac[d];
3065 while (cellsize*npulse[d] < comm->cutoff &&
3066 npulse[d] < dd->nc[d]-1)
3070 cellsize_min[d] = min(cellsize_min[d],cellsize);
3074 comm->cell_x0[d] = cell_x[dd->ci[d]];
3075 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
3079 /* The following limitation is to avoid that a cell would receive
3080 * some of its own home charge groups back over the periodic boundary.
3081 * Double charge groups cause trouble with the global indices.
3083 if (d < ddbox->npbcdim &&
3084 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
3086 gmx_fatal_collective(FARGS,NULL,dd,
3087 "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",
3088 dim2char(d),ddbox->box_size[d],ddbox->skew_fac[d],
3090 dd->nc[d],dd->nc[d],
3091 dd->nnodes > dd->nc[d] ? "cells" : "processors");
3095 if (!comm->bDynLoadBal)
3097 copy_rvec(cellsize_min,comm->cellsize_min);
3100 for(d=0; d<comm->npmedecompdim; d++)
3102 set_pme_maxshift(dd,&comm->ddpme[d],
3103 comm->slb_frac[dd->dim[d]]==NULL,ddbox,
3104 comm->ddpme[d].slb_dim_f);
3109 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
3110 int d,int dim,gmx_domdec_root_t *root,
3112 gmx_bool bUniform,gmx_large_int_t step, real cellsize_limit_f, int range[])
3114 gmx_domdec_comm_t *comm;
3115 int ncd,i,j,nmin,nmin_old;
3116 gmx_bool bLimLo,bLimHi;
3118 real fac,halfway,cellsize_limit_f_i,region_size;
3119 gmx_bool bPBC,bLastHi=FALSE;
3120 int nrange[]={range[0],range[1]};
3122 region_size= root->cell_f[range[1]]-root->cell_f[range[0]];
3128 bPBC = (dim < ddbox->npbcdim);
3130 cell_size = root->buf_ncd;
3134 fprintf(debug,"enforce_limits: %d %d\n",range[0],range[1]);
3137 /* First we need to check if the scaling does not make cells
3138 * smaller than the smallest allowed size.
3139 * We need to do this iteratively, since if a cell is too small,
3140 * it needs to be enlarged, which makes all the other cells smaller,
3141 * which could in turn make another cell smaller than allowed.
3143 for(i=range[0]; i<range[1]; i++)
3145 root->bCellMin[i] = FALSE;
3151 /* We need the total for normalization */
3153 for(i=range[0]; i<range[1]; i++)
3155 if (root->bCellMin[i] == FALSE)
3157 fac += cell_size[i];
3160 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3161 /* Determine the cell boundaries */
3162 for(i=range[0]; i<range[1]; i++)
3164 if (root->bCellMin[i] == FALSE)
3166 cell_size[i] *= fac;
3167 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3169 cellsize_limit_f_i = 0;
3173 cellsize_limit_f_i = cellsize_limit_f;
3175 if (cell_size[i] < cellsize_limit_f_i)
3177 root->bCellMin[i] = TRUE;
3178 cell_size[i] = cellsize_limit_f_i;
3182 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3185 while (nmin > nmin_old);
3188 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3189 /* For this check we should not use DD_CELL_MARGIN,
3190 * but a slightly smaller factor,
3191 * since rounding could get use below the limit.
3193 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3196 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",
3197 gmx_step_str(step,buf),
3198 dim2char(dim),ddbox->box_size[dim],ddbox->skew_fac[dim],
3199 ncd,comm->cellsize_min[dim]);
3202 root->bLimited = (nmin > 0) || (range[0]>0) || (range[1]<ncd);
3206 /* Check if the boundary did not displace more than halfway
3207 * each of the cells it bounds, as this could cause problems,
3208 * especially when the differences between cell sizes are large.
3209 * If changes are applied, they will not make cells smaller
3210 * than the cut-off, as we check all the boundaries which
3211 * might be affected by a change and if the old state was ok,
3212 * the cells will at most be shrunk back to their old size.
3214 for(i=range[0]+1; i<range[1]; i++)
3216 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3217 if (root->cell_f[i] < halfway)
3219 root->cell_f[i] = halfway;
3220 /* Check if the change also causes shifts of the next boundaries */
3221 for(j=i+1; j<range[1]; j++)
3223 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3224 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3227 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3228 if (root->cell_f[i] > halfway)
3230 root->cell_f[i] = halfway;
3231 /* Check if the change also causes shifts of the next boundaries */
3232 for(j=i-1; j>=range[0]+1; j--)
3234 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3235 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3241 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3242 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3243 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3244 * for a and b nrange is used */
3247 /* Take care of the staggering of the cell boundaries */
3250 for(i=range[0]; i<range[1]; i++)
3252 root->cell_f_max0[i] = root->cell_f[i];
3253 root->cell_f_min1[i] = root->cell_f[i+1];
3258 for(i=range[0]+1; i<range[1]; i++)
3260 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3261 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3262 if (bLimLo && bLimHi)
3264 /* Both limits violated, try the best we can */
3265 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3266 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3269 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3273 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3279 /* root->cell_f[i] = root->bound_min[i]; */
3280 nrange[1]=i; /* only store violation location. There could be a LimLo violation following with an higher index */
3283 else if (bLimHi && !bLastHi)
3286 if (nrange[1] < range[1]) /* found a LimLo before */
3288 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3289 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3290 nrange[0]=nrange[1];
3292 root->cell_f[i] = root->bound_max[i];
3294 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3299 if (nrange[1] < range[1]) /* found last a LimLo */
3301 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3302 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3303 nrange[0]=nrange[1];
3305 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3307 else if (nrange[0] > range[0]) /* found at least one LimHi */
3309 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3316 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3317 int d,int dim,gmx_domdec_root_t *root,
3318 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3319 gmx_bool bUniform,gmx_large_int_t step)
3321 gmx_domdec_comm_t *comm;
3324 real load_aver,load_i,imbalance,change,change_max,sc;
3325 real cellsize_limit_f,dist_min_f,dist_min_f_hard,space;
3329 int range[] = { 0, 0 };
3333 /* Convert the maximum change from the input percentage to a fraction */
3334 change_limit = comm->dlb_scale_lim*0.01;
3338 bPBC = (dim < ddbox->npbcdim);
3340 cell_size = root->buf_ncd;
3342 /* Store the original boundaries */
3343 for(i=0; i<ncd+1; i++)
3345 root->old_cell_f[i] = root->cell_f[i];
3348 for(i=0; i<ncd; i++)
3350 cell_size[i] = 1.0/ncd;
3353 else if (dd_load_count(comm))
3355 load_aver = comm->load[d].sum_m/ncd;
3357 for(i=0; i<ncd; i++)
3359 /* Determine the relative imbalance of cell i */
3360 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3361 imbalance = (load_i - load_aver)/(load_aver>0 ? load_aver : 1);
3362 /* Determine the change of the cell size using underrelaxation */
3363 change = -relax*imbalance;
3364 change_max = max(change_max,max(change,-change));
3366 /* Limit the amount of scaling.
3367 * We need to use the same rescaling for all cells in one row,
3368 * otherwise the load balancing might not converge.
3371 if (change_max > change_limit)
3373 sc *= change_limit/change_max;
3375 for(i=0; i<ncd; i++)
3377 /* Determine the relative imbalance of cell i */
3378 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3379 imbalance = (load_i - load_aver)/(load_aver>0 ? load_aver : 1);
3380 /* Determine the change of the cell size using underrelaxation */
3381 change = -sc*imbalance;
3382 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3386 cellsize_limit_f = cellsize_min_dlb(comm,d,dim)/ddbox->box_size[dim];
3387 cellsize_limit_f *= DD_CELL_MARGIN;
3388 dist_min_f_hard = grid_jump_limit(comm,comm->cutoff,d)/ddbox->box_size[dim];
3389 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3390 if (ddbox->tric_dir[dim])
3392 cellsize_limit_f /= ddbox->skew_fac[dim];
3393 dist_min_f /= ddbox->skew_fac[dim];
3395 if (bDynamicBox && d > 0)
3397 dist_min_f *= DD_PRES_SCALE_MARGIN;
3399 if (d > 0 && !bUniform)
3401 /* Make sure that the grid is not shifted too much */
3402 for(i=1; i<ncd; i++) {
3403 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3405 gmx_incons("Inconsistent DD boundary staggering limits!");
3407 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3408 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3410 root->bound_min[i] += 0.5*space;
3412 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3413 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3415 root->bound_max[i] += 0.5*space;
3420 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3422 root->cell_f_max0[i-1] + dist_min_f,
3423 root->bound_min[i],root->cell_f[i],root->bound_max[i],
3424 root->cell_f_min1[i] - dist_min_f);
3429 root->cell_f[0] = 0;
3430 root->cell_f[ncd] = 1;
3431 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3434 /* After the checks above, the cells should obey the cut-off
3435 * restrictions, but it does not hurt to check.
3437 for(i=0; i<ncd; i++)
3441 fprintf(debug,"Relative bounds dim %d cell %d: %f %f\n",
3442 dim,i,root->cell_f[i],root->cell_f[i+1]);
3445 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3446 root->cell_f[i+1] - root->cell_f[i] <
3447 cellsize_limit_f/DD_CELL_MARGIN)
3451 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3452 gmx_step_str(step,buf),dim2char(dim),i,
3453 (root->cell_f[i+1] - root->cell_f[i])
3454 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3459 /* Store the cell boundaries of the lower dimensions at the end */
3460 for(d1=0; d1<d; d1++)
3462 root->cell_f[pos++] = comm->cell_f0[d1];
3463 root->cell_f[pos++] = comm->cell_f1[d1];
3466 if (d < comm->npmedecompdim)
3468 /* The master determines the maximum shift for
3469 * the coordinate communication between separate PME nodes.
3471 set_pme_maxshift(dd,&comm->ddpme[d],bUniform,ddbox,root->cell_f);
3473 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3476 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3480 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3481 gmx_ddbox_t *ddbox,int dimind)
3483 gmx_domdec_comm_t *comm;
3488 /* Set the cell dimensions */
3489 dim = dd->dim[dimind];
3490 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3491 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3492 if (dim >= ddbox->nboundeddim)
3494 comm->cell_x0[dim] += ddbox->box0[dim];
3495 comm->cell_x1[dim] += ddbox->box0[dim];
3499 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3500 int d,int dim,real *cell_f_row,
3503 gmx_domdec_comm_t *comm;
3509 /* Each node would only need to know two fractions,
3510 * but it is probably cheaper to broadcast the whole array.
3512 MPI_Bcast(cell_f_row,DD_CELL_F_SIZE(dd,d)*sizeof(real),MPI_BYTE,
3513 0,comm->mpi_comm_load[d]);
3515 /* Copy the fractions for this dimension from the buffer */
3516 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3517 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3518 /* The whole array was communicated, so set the buffer position */
3519 pos = dd->nc[dim] + 1;
3520 for(d1=0; d1<=d; d1++)
3524 /* Copy the cell fractions of the lower dimensions */
3525 comm->cell_f0[d1] = cell_f_row[pos++];
3526 comm->cell_f1[d1] = cell_f_row[pos++];
3528 relative_to_absolute_cell_bounds(dd,ddbox,d1);
3530 /* Convert the communicated shift from float to int */
3531 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3534 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3538 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3539 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3540 gmx_bool bUniform,gmx_large_int_t step)
3542 gmx_domdec_comm_t *comm;
3544 gmx_bool bRowMember,bRowRoot;
3549 for(d=0; d<dd->ndim; d++)
3554 for(d1=d; d1<dd->ndim; d1++)
3556 if (dd->ci[dd->dim[d1]] > 0)
3569 set_dd_cell_sizes_dlb_root(dd,d,dim,comm->root[d],
3570 ddbox,bDynamicBox,bUniform,step);
3571 cell_f_row = comm->root[d]->cell_f;
3575 cell_f_row = comm->cell_f_row;
3577 distribute_dd_cell_sizes_dlb(dd,d,dim,cell_f_row,ddbox);
3582 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
3586 /* This function assumes the box is static and should therefore
3587 * not be called when the box has changed since the last
3588 * call to dd_partition_system.
3590 for(d=0; d<dd->ndim; d++)
3592 relative_to_absolute_cell_bounds(dd,ddbox,d);
3598 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3599 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3600 gmx_bool bUniform,gmx_bool bDoDLB,gmx_large_int_t step,
3601 gmx_wallcycle_t wcycle)
3603 gmx_domdec_comm_t *comm;
3610 wallcycle_start(wcycle,ewcDDCOMMBOUND);
3611 set_dd_cell_sizes_dlb_change(dd,ddbox,bDynamicBox,bUniform,step);
3612 wallcycle_stop(wcycle,ewcDDCOMMBOUND);
3614 else if (bDynamicBox)
3616 set_dd_cell_sizes_dlb_nochange(dd,ddbox);
3619 /* Set the dimensions for which no DD is used */
3620 for(dim=0; dim<DIM; dim++) {
3621 if (dd->nc[dim] == 1) {
3622 comm->cell_x0[dim] = 0;
3623 comm->cell_x1[dim] = ddbox->box_size[dim];
3624 if (dim >= ddbox->nboundeddim)
3626 comm->cell_x0[dim] += ddbox->box0[dim];
3627 comm->cell_x1[dim] += ddbox->box0[dim];
3633 static void realloc_comm_ind(gmx_domdec_t *dd,ivec npulse)
3636 gmx_domdec_comm_dim_t *cd;
3638 for(d=0; d<dd->ndim; d++)
3640 cd = &dd->comm->cd[d];
3641 np = npulse[dd->dim[d]];
3642 if (np > cd->np_nalloc)
3646 fprintf(debug,"(Re)allocing cd for %c to %d pulses\n",
3647 dim2char(dd->dim[d]),np);
3649 if (DDMASTER(dd) && cd->np_nalloc > 0)
3651 fprintf(stderr,"\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n",dim2char(dd->dim[d]),np);
3654 for(i=cd->np_nalloc; i<np; i++)
3656 cd->ind[i].index = NULL;
3657 cd->ind[i].nalloc = 0;
3666 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3667 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3668 gmx_bool bUniform,gmx_bool bDoDLB,gmx_large_int_t step,
3669 gmx_wallcycle_t wcycle)
3671 gmx_domdec_comm_t *comm;
3677 /* Copy the old cell boundaries for the cg displacement check */
3678 copy_rvec(comm->cell_x0,comm->old_cell_x0);
3679 copy_rvec(comm->cell_x1,comm->old_cell_x1);
3681 if (comm->bDynLoadBal)
3685 check_box_size(dd,ddbox);
3687 set_dd_cell_sizes_dlb(dd,ddbox,bDynamicBox,bUniform,bDoDLB,step,wcycle);
3691 set_dd_cell_sizes_slb(dd,ddbox,FALSE,npulse);
3692 realloc_comm_ind(dd,npulse);
3697 for(d=0; d<DIM; d++)
3699 fprintf(debug,"cell_x[%d] %f - %f skew_fac %f\n",
3700 d,comm->cell_x0[d],comm->cell_x1[d],ddbox->skew_fac[d]);
3705 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3707 rvec cell_ns_x0,rvec cell_ns_x1,
3708 gmx_large_int_t step)
3710 gmx_domdec_comm_t *comm;
3715 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
3717 dim = dd->dim[dim_ind];
3719 /* Without PBC we don't have restrictions on the outer cells */
3720 if (!(dim >= ddbox->npbcdim &&
3721 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3722 comm->bDynLoadBal &&
3723 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3724 comm->cellsize_min[dim])
3727 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",
3728 gmx_step_str(step,buf),dim2char(dim),
3729 comm->cell_x1[dim] - comm->cell_x0[dim],
3730 ddbox->skew_fac[dim],
3731 dd->comm->cellsize_min[dim],
3732 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
3736 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3738 /* Communicate the boundaries and update cell_ns_x0/1 */
3739 dd_move_cellx(dd,ddbox,cell_ns_x0,cell_ns_x1);
3740 if (dd->bGridJump && dd->ndim > 1)
3742 check_grid_jump(step,dd,dd->comm->cutoff,ddbox,TRUE);
3747 static void make_tric_corr_matrix(int npbcdim,matrix box,matrix tcm)
3751 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3759 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3760 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3769 static void check_screw_box(matrix box)
3771 /* Mathematical limitation */
3772 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3774 gmx_fatal(FARGS,"With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3777 /* Limitation due to the asymmetry of the eighth shell method */
3778 if (box[ZZ][YY] != 0)
3780 gmx_fatal(FARGS,"pbc=screw with non-zero box_zy is not supported");
3784 static void distribute_cg(FILE *fplog,gmx_large_int_t step,
3785 matrix box,ivec tric_dir,t_block *cgs,rvec pos[],
3788 gmx_domdec_master_t *ma;
3789 int **tmp_ind=NULL,*tmp_nalloc=NULL;
3790 int i,icg,j,k,k0,k1,d,npbcdim;
3792 rvec box_size,cg_cm;
3794 real nrcg,inv_ncg,pos_d;
3796 gmx_bool bUnbounded,bScrew;
3800 if (tmp_ind == NULL)
3802 snew(tmp_nalloc,dd->nnodes);
3803 snew(tmp_ind,dd->nnodes);
3804 for(i=0; i<dd->nnodes; i++)
3806 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3807 snew(tmp_ind[i],tmp_nalloc[i]);
3811 /* Clear the count */
3812 for(i=0; i<dd->nnodes; i++)
3818 make_tric_corr_matrix(dd->npbcdim,box,tcm);
3820 cgindex = cgs->index;
3822 /* Compute the center of geometry for all charge groups */
3823 for(icg=0; icg<cgs->nr; icg++)
3826 k1 = cgindex[icg+1];
3830 copy_rvec(pos[k0],cg_cm);
3837 for(k=k0; (k<k1); k++)
3839 rvec_inc(cg_cm,pos[k]);
3841 for(d=0; (d<DIM); d++)
3843 cg_cm[d] *= inv_ncg;
3846 /* Put the charge group in the box and determine the cell index */
3847 for(d=DIM-1; d>=0; d--) {
3849 if (d < dd->npbcdim)
3851 bScrew = (dd->bScrewPBC && d == XX);
3852 if (tric_dir[d] && dd->nc[d] > 1)
3854 /* Use triclinic coordintates for this dimension */
3855 for(j=d+1; j<DIM; j++)
3857 pos_d += cg_cm[j]*tcm[j][d];
3860 while(pos_d >= box[d][d])
3863 rvec_dec(cg_cm,box[d]);
3866 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3867 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3869 for(k=k0; (k<k1); k++)
3871 rvec_dec(pos[k],box[d]);
3874 pos[k][YY] = box[YY][YY] - pos[k][YY];
3875 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3882 rvec_inc(cg_cm,box[d]);
3885 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3886 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3888 for(k=k0; (k<k1); k++)
3890 rvec_inc(pos[k],box[d]);
3892 pos[k][YY] = box[YY][YY] - pos[k][YY];
3893 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3898 /* This could be done more efficiently */
3900 while(ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3905 i = dd_index(dd->nc,ind);
3906 if (ma->ncg[i] == tmp_nalloc[i])
3908 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3909 srenew(tmp_ind[i],tmp_nalloc[i]);
3911 tmp_ind[i][ma->ncg[i]] = icg;
3913 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3917 for(i=0; i<dd->nnodes; i++)
3920 for(k=0; k<ma->ncg[i]; k++)
3922 ma->cg[k1++] = tmp_ind[i][k];
3925 ma->index[dd->nnodes] = k1;
3927 for(i=0; i<dd->nnodes; i++)
3937 fprintf(fplog,"Charge group distribution at step %s:",
3938 gmx_step_str(step,buf));
3939 for(i=0; i<dd->nnodes; i++)
3941 fprintf(fplog," %d",ma->ncg[i]);
3943 fprintf(fplog,"\n");
3947 static void get_cg_distribution(FILE *fplog,gmx_large_int_t step,gmx_domdec_t *dd,
3948 t_block *cgs,matrix box,gmx_ddbox_t *ddbox,
3951 gmx_domdec_master_t *ma=NULL;
3954 int *ibuf,buf2[2] = { 0, 0 };
3955 gmx_bool bMaster = DDMASTER(dd);
3962 check_screw_box(box);
3965 set_dd_cell_sizes_slb(dd,ddbox,TRUE,npulse);
3967 distribute_cg(fplog,step,box,ddbox->tric_dir,cgs,pos,dd);
3968 for(i=0; i<dd->nnodes; i++)
3970 ma->ibuf[2*i] = ma->ncg[i];
3971 ma->ibuf[2*i+1] = ma->nat[i];
3979 dd_scatter(dd,2*sizeof(int),ibuf,buf2);
3981 dd->ncg_home = buf2[0];
3982 dd->nat_home = buf2[1];
3983 dd->ncg_tot = dd->ncg_home;
3984 dd->nat_tot = dd->nat_home;
3985 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
3987 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
3988 srenew(dd->index_gl,dd->cg_nalloc);
3989 srenew(dd->cgindex,dd->cg_nalloc+1);
3993 for(i=0; i<dd->nnodes; i++)
3995 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
3996 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
4001 DDMASTER(dd) ? ma->ibuf : NULL,
4002 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
4003 DDMASTER(dd) ? ma->cg : NULL,
4004 dd->ncg_home*sizeof(int),dd->index_gl);
4006 /* Determine the home charge group sizes */
4008 for(i=0; i<dd->ncg_home; i++)
4010 cg_gl = dd->index_gl[i];
4012 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
4017 fprintf(debug,"Home charge groups:\n");
4018 for(i=0; i<dd->ncg_home; i++)
4020 fprintf(debug," %d",dd->index_gl[i]);
4022 fprintf(debug,"\n");
4024 fprintf(debug,"\n");
4028 static int compact_and_copy_vec_at(int ncg,int *move,
4031 rvec *src,gmx_domdec_comm_t *comm,
4034 int m,icg,i,i0,i1,nrcg;
4040 for(m=0; m<DIM*2; m++)
4046 for(icg=0; icg<ncg; icg++)
4048 i1 = cgindex[icg+1];
4054 /* Compact the home array in place */
4055 for(i=i0; i<i1; i++)
4057 copy_rvec(src[i],src[home_pos++]);
4063 /* Copy to the communication buffer */
4065 pos_vec[m] += 1 + vec*nrcg;
4066 for(i=i0; i<i1; i++)
4068 copy_rvec(src[i],comm->cgcm_state[m][pos_vec[m]++]);
4070 pos_vec[m] += (nvec - vec - 1)*nrcg;
4074 home_pos += i1 - i0;
4082 static int compact_and_copy_vec_cg(int ncg,int *move,
4084 int nvec,rvec *src,gmx_domdec_comm_t *comm,
4087 int m,icg,i0,i1,nrcg;
4093 for(m=0; m<DIM*2; m++)
4099 for(icg=0; icg<ncg; icg++)
4101 i1 = cgindex[icg+1];
4107 /* Compact the home array in place */
4108 copy_rvec(src[icg],src[home_pos++]);
4114 /* Copy to the communication buffer */
4115 copy_rvec(src[icg],comm->cgcm_state[m][pos_vec[m]]);
4116 pos_vec[m] += 1 + nrcg*nvec;
4128 static int compact_ind(int ncg,int *move,
4129 int *index_gl,int *cgindex,
4131 gmx_ga2la_t ga2la,char *bLocalCG,
4134 int cg,nat,a0,a1,a,a_gl;
4139 for(cg=0; cg<ncg; cg++)
4145 /* Compact the home arrays in place.
4146 * Anything that can be done here avoids access to global arrays.
4148 cgindex[home_pos] = nat;
4149 for(a=a0; a<a1; a++)
4152 gatindex[nat] = a_gl;
4153 /* The cell number stays 0, so we don't need to set it */
4154 ga2la_change_la(ga2la,a_gl,nat);
4157 index_gl[home_pos] = index_gl[cg];
4158 cginfo[home_pos] = cginfo[cg];
4159 /* The charge group remains local, so bLocalCG does not change */
4164 /* Clear the global indices */
4165 for(a=a0; a<a1; a++)
4167 ga2la_del(ga2la,gatindex[a]);
4171 bLocalCG[index_gl[cg]] = FALSE;
4175 cgindex[home_pos] = nat;
4180 static void clear_and_mark_ind(int ncg,int *move,
4181 int *index_gl,int *cgindex,int *gatindex,
4182 gmx_ga2la_t ga2la,char *bLocalCG,
4187 for(cg=0; cg<ncg; cg++)
4193 /* Clear the global indices */
4194 for(a=a0; a<a1; a++)
4196 ga2la_del(ga2la,gatindex[a]);
4200 bLocalCG[index_gl[cg]] = FALSE;
4202 /* Signal that this cg has moved using the ns cell index.
4203 * Here we set it to -1. fill_grid will change it
4204 * from -1 to NSGRID_SIGNAL_MOVED_FAC*grid->ncells.
4206 cell_index[cg] = -1;
4211 static void print_cg_move(FILE *fplog,
4213 gmx_large_int_t step,int cg,int dim,int dir,
4214 gmx_bool bHaveLimitdAndCMOld,real limitd,
4215 rvec cm_old,rvec cm_new,real pos_d)
4217 gmx_domdec_comm_t *comm;
4222 fprintf(fplog,"\nStep %s:\n",gmx_step_str(step,buf));
4223 if (bHaveLimitdAndCMOld)
4225 fprintf(fplog,"The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
4226 ddglatnr(dd,dd->cgindex[cg]),limitd,dim2char(dim));
4230 fprintf(fplog,"The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4231 ddglatnr(dd,dd->cgindex[cg]),dim2char(dim));
4233 fprintf(fplog,"distance out of cell %f\n",
4234 dir==1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4235 if (bHaveLimitdAndCMOld)
4237 fprintf(fplog,"Old coordinates: %8.3f %8.3f %8.3f\n",
4238 cm_old[XX],cm_old[YY],cm_old[ZZ]);
4240 fprintf(fplog,"New coordinates: %8.3f %8.3f %8.3f\n",
4241 cm_new[XX],cm_new[YY],cm_new[ZZ]);
4242 fprintf(fplog,"Old cell boundaries in direction %c: %8.3f %8.3f\n",
4244 comm->old_cell_x0[dim],comm->old_cell_x1[dim]);
4245 fprintf(fplog,"New cell boundaries in direction %c: %8.3f %8.3f\n",
4247 comm->cell_x0[dim],comm->cell_x1[dim]);
4250 static void cg_move_error(FILE *fplog,
4252 gmx_large_int_t step,int cg,int dim,int dir,
4253 gmx_bool bHaveLimitdAndCMOld,real limitd,
4254 rvec cm_old,rvec cm_new,real pos_d)
4258 print_cg_move(fplog, dd,step,cg,dim,dir,
4259 bHaveLimitdAndCMOld,limitd,cm_old,cm_new,pos_d);
4261 print_cg_move(stderr,dd,step,cg,dim,dir,
4262 bHaveLimitdAndCMOld,limitd,cm_old,cm_new,pos_d);
4264 "A charge group moved too far between two domain decomposition steps\n"
4265 "This usually means that your system is not well equilibrated");
4268 static void rotate_state_atom(t_state *state,int a)
4272 for(est=0; est<estNR; est++)
4274 if (EST_DISTR(est) && (state->flags & (1<<est))) {
4277 /* Rotate the complete state; for a rectangular box only */
4278 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4279 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4282 state->v[a][YY] = -state->v[a][YY];
4283 state->v[a][ZZ] = -state->v[a][ZZ];
4286 state->sd_X[a][YY] = -state->sd_X[a][YY];
4287 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4290 state->cg_p[a][YY] = -state->cg_p[a][YY];
4291 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4293 case estDISRE_INITF:
4294 case estDISRE_RM3TAV:
4295 case estORIRE_INITF:
4297 /* These are distances, so not affected by rotation */
4300 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4306 static int *get_moved(gmx_domdec_comm_t *comm,int natoms)
4308 if (natoms > comm->moved_nalloc)
4310 /* Contents should be preserved here */
4311 comm->moved_nalloc = over_alloc_dd(natoms);
4312 srenew(comm->moved,comm->moved_nalloc);
4318 static void calc_cg_move(FILE *fplog,gmx_large_int_t step,
4321 ivec tric_dir,matrix tcm,
4322 rvec cell_x0,rvec cell_x1,
4323 rvec limitd,rvec limit0,rvec limit1,
4325 int cg_start,int cg_end,
4330 int c,i,cg,k,k0,k1,d,dim,dim2,dir,d2,d3,d4,cell_d;
4331 int mc,cdd,nrcg,ncg_recv,nat_recv,nvs,nvr,nvec,vec;
4338 npbcdim = dd->npbcdim;
4340 for(cg=cg_start; cg<cg_end; cg++)
4347 copy_rvec(state->x[k0],cm_new);
4354 for(k=k0; (k<k1); k++)
4356 rvec_inc(cm_new,state->x[k]);
4358 for(d=0; (d<DIM); d++)
4360 cm_new[d] = inv_ncg*cm_new[d];
4365 /* Do pbc and check DD cell boundary crossings */
4366 for(d=DIM-1; d>=0; d--)
4370 bScrew = (dd->bScrewPBC && d == XX);
4371 /* Determine the location of this cg in lattice coordinates */
4375 for(d2=d+1; d2<DIM; d2++)
4377 pos_d += cm_new[d2]*tcm[d2][d];
4380 /* Put the charge group in the triclinic unit-cell */
4381 if (pos_d >= cell_x1[d])
4383 if (pos_d >= limit1[d])
4385 cg_move_error(fplog,dd,step,cg,d,1,TRUE,limitd[d],
4386 cg_cm[cg],cm_new,pos_d);
4389 if (dd->ci[d] == dd->nc[d] - 1)
4391 rvec_dec(cm_new,state->box[d]);
4394 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4395 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4397 for(k=k0; (k<k1); k++)
4399 rvec_dec(state->x[k],state->box[d]);
4402 rotate_state_atom(state,k);
4407 else if (pos_d < cell_x0[d])
4409 if (pos_d < limit0[d])
4411 cg_move_error(fplog,dd,step,cg,d,-1,TRUE,limitd[d],
4412 cg_cm[cg],cm_new,pos_d);
4417 rvec_inc(cm_new,state->box[d]);
4420 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4421 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4423 for(k=k0; (k<k1); k++)
4425 rvec_inc(state->x[k],state->box[d]);
4428 rotate_state_atom(state,k);
4434 else if (d < npbcdim)
4436 /* Put the charge group in the rectangular unit-cell */
4437 while (cm_new[d] >= state->box[d][d])
4439 rvec_dec(cm_new,state->box[d]);
4440 for(k=k0; (k<k1); k++)
4442 rvec_dec(state->x[k],state->box[d]);
4445 while (cm_new[d] < 0)
4447 rvec_inc(cm_new,state->box[d]);
4448 for(k=k0; (k<k1); k++)
4450 rvec_inc(state->x[k],state->box[d]);
4456 copy_rvec(cm_new,cg_cm[cg]);
4458 /* Determine where this cg should go */
4461 for(d=0; d<dd->ndim; d++)
4466 flag |= DD_FLAG_FW(d);
4472 else if (dev[dim] == -1)
4474 flag |= DD_FLAG_BW(d);
4476 if (dd->nc[dim] > 2)
4487 /* Temporarily store the flag in move */
4488 move[cg] = mc + flag;
4492 static void dd_redistribute_cg(FILE *fplog,gmx_large_int_t step,
4493 gmx_domdec_t *dd,ivec tric_dir,
4494 t_state *state,rvec **f,
4495 t_forcerec *fr,t_mdatoms *md,
4503 int ncg[DIM*2],nat[DIM*2];
4504 int c,i,cg,k,k0,k1,d,dim,dim2,dir,d2,d3,d4,cell_d;
4505 int mc,cdd,nrcg,ncg_recv,nat_recv,nvs,nvr,nvec,vec;
4506 int sbuf[2],rbuf[2];
4507 int home_pos_cg,home_pos_at,buf_pos;
4509 gmx_bool bV=FALSE,bSDX=FALSE,bCGP=FALSE;
4514 rvec *cg_cm=NULL,cell_x0,cell_x1,limitd,limit0,limit1,cm_new;
4516 cginfo_mb_t *cginfo_mb;
4517 gmx_domdec_comm_t *comm;
4523 check_screw_box(state->box);
4527 if (fr->cutoff_scheme == ecutsGROUP)
4532 for(i=0; i<estNR; i++)
4538 case estX: /* Always present */ break;
4539 case estV: bV = (state->flags & (1<<i)); break;
4540 case estSDX: bSDX = (state->flags & (1<<i)); break;
4541 case estCGP: bCGP = (state->flags & (1<<i)); break;
4544 case estDISRE_INITF:
4545 case estDISRE_RM3TAV:
4546 case estORIRE_INITF:
4548 /* No processing required */
4551 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4556 if (dd->ncg_tot > comm->nalloc_int)
4558 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4559 srenew(comm->buf_int,comm->nalloc_int);
4561 move = comm->buf_int;
4563 /* Clear the count */
4564 for(c=0; c<dd->ndim*2; c++)
4570 npbcdim = dd->npbcdim;
4572 for(d=0; (d<DIM); d++)
4574 limitd[d] = dd->comm->cellsize_min[d];
4575 if (d >= npbcdim && dd->ci[d] == 0)
4577 cell_x0[d] = -GMX_FLOAT_MAX;
4581 cell_x0[d] = comm->cell_x0[d];
4583 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4585 cell_x1[d] = GMX_FLOAT_MAX;
4589 cell_x1[d] = comm->cell_x1[d];
4593 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4594 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4598 /* We check after communication if a charge group moved
4599 * more than one cell. Set the pre-comm check limit to float_max.
4601 limit0[d] = -GMX_FLOAT_MAX;
4602 limit1[d] = GMX_FLOAT_MAX;
4606 make_tric_corr_matrix(npbcdim,state->box,tcm);
4608 cgindex = dd->cgindex;
4610 nthread = gmx_omp_nthreads_get(emntDomdec);
4612 /* Compute the center of geometry for all home charge groups
4613 * and put them in the box and determine where they should go.
4615 #pragma omp parallel for num_threads(nthread) schedule(static)
4616 for(thread=0; thread<nthread; thread++)
4618 calc_cg_move(fplog,step,dd,state,tric_dir,tcm,
4619 cell_x0,cell_x1,limitd,limit0,limit1,
4621 ( thread *dd->ncg_home)/nthread,
4622 ((thread+1)*dd->ncg_home)/nthread,
4623 fr->cutoff_scheme==ecutsGROUP ? cg_cm : state->x,
4627 for(cg=0; cg<dd->ncg_home; cg++)
4632 flag = mc & ~DD_FLAG_NRCG;
4633 mc = mc & DD_FLAG_NRCG;
4636 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4638 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4639 srenew(comm->cggl_flag[mc],comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4641 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4642 /* We store the cg size in the lower 16 bits
4643 * and the place where the charge group should go
4644 * in the next 6 bits. This saves some communication volume.
4646 nrcg = cgindex[cg+1] - cgindex[cg];
4647 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4653 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
4654 inc_nrnb(nrnb,eNR_RESETX,dd->ncg_home);
4657 for(i=0; i<dd->ndim*2; i++)
4659 *ncg_moved += ncg[i];
4676 /* Make sure the communication buffers are large enough */
4677 for(mc=0; mc<dd->ndim*2; mc++)
4679 nvr = ncg[mc] + nat[mc]*nvec;
4680 if (nvr > comm->cgcm_state_nalloc[mc])
4682 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4683 srenew(comm->cgcm_state[mc],comm->cgcm_state_nalloc[mc]);
4687 switch (fr->cutoff_scheme)
4690 /* Recalculating cg_cm might be cheaper than communicating,
4691 * but that could give rise to rounding issues.
4694 compact_and_copy_vec_cg(dd->ncg_home,move,cgindex,
4695 nvec,cg_cm,comm,bCompact);
4698 /* Without charge groups we send the moved atom coordinates
4699 * over twice. This is so the code below can be used without
4700 * many conditionals for both for with and without charge groups.
4703 compact_and_copy_vec_cg(dd->ncg_home,move,cgindex,
4704 nvec,state->x,comm,FALSE);
4707 home_pos_cg -= *ncg_moved;
4711 gmx_incons("unimplemented");
4717 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4718 nvec,vec++,state->x,comm,bCompact);
4721 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4722 nvec,vec++,state->v,comm,bCompact);
4726 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4727 nvec,vec++,state->sd_X,comm,bCompact);
4731 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4732 nvec,vec++,state->cg_p,comm,bCompact);
4737 compact_ind(dd->ncg_home,move,
4738 dd->index_gl,dd->cgindex,dd->gatindex,
4739 dd->ga2la,comm->bLocalCG,
4744 if (fr->cutoff_scheme == ecutsVERLET)
4746 moved = get_moved(comm,dd->ncg_home);
4748 for(k=0; k<dd->ncg_home; k++)
4755 moved = fr->ns.grid->cell_index;
4758 clear_and_mark_ind(dd->ncg_home,move,
4759 dd->index_gl,dd->cgindex,dd->gatindex,
4760 dd->ga2la,comm->bLocalCG,
4764 cginfo_mb = fr->cginfo_mb;
4766 *ncg_stay_home = home_pos_cg;
4767 for(d=0; d<dd->ndim; d++)
4773 for(dir=0; dir<(dd->nc[dim]==2 ? 1 : 2); dir++)
4776 /* Communicate the cg and atom counts */
4781 fprintf(debug,"Sending ddim %d dir %d: ncg %d nat %d\n",
4782 d,dir,sbuf[0],sbuf[1]);
4784 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4786 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4788 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4789 srenew(comm->buf_int,comm->nalloc_int);
4792 /* Communicate the charge group indices, sizes and flags */
4793 dd_sendrecv_int(dd, d, dir,
4794 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4795 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4797 nvs = ncg[cdd] + nat[cdd]*nvec;
4798 i = rbuf[0] + rbuf[1] *nvec;
4799 vec_rvec_check_alloc(&comm->vbuf,nvr+i);
4801 /* Communicate cgcm and state */
4802 dd_sendrecv_rvec(dd, d, dir,
4803 comm->cgcm_state[cdd], nvs,
4804 comm->vbuf.v+nvr, i);
4805 ncg_recv += rbuf[0];
4806 nat_recv += rbuf[1];
4810 /* Process the received charge groups */
4812 for(cg=0; cg<ncg_recv; cg++)
4814 flag = comm->buf_int[cg*DD_CGIBS+1];
4816 if (dim >= npbcdim && dd->nc[dim] > 2)
4818 /* No pbc in this dim and more than one domain boundary.
4819 * We do a separate check if a charge group didn't move too far.
4821 if (((flag & DD_FLAG_FW(d)) &&
4822 comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
4823 ((flag & DD_FLAG_BW(d)) &&
4824 comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
4826 cg_move_error(fplog,dd,step,cg,dim,
4827 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4829 comm->vbuf.v[buf_pos],
4830 comm->vbuf.v[buf_pos],
4831 comm->vbuf.v[buf_pos][dim]);
4838 /* Check which direction this cg should go */
4839 for(d2=d+1; (d2<dd->ndim && mc==-1); d2++)
4843 /* The cell boundaries for dimension d2 are not equal
4844 * for each cell row of the lower dimension(s),
4845 * therefore we might need to redetermine where
4846 * this cg should go.
4849 /* If this cg crosses the box boundary in dimension d2
4850 * we can use the communicated flag, so we do not
4851 * have to worry about pbc.
4853 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4854 (flag & DD_FLAG_FW(d2))) ||
4855 (dd->ci[dim2] == 0 &&
4856 (flag & DD_FLAG_BW(d2)))))
4858 /* Clear the two flags for this dimension */
4859 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4860 /* Determine the location of this cg
4861 * in lattice coordinates
4863 pos_d = comm->vbuf.v[buf_pos][dim2];
4866 for(d3=dim2+1; d3<DIM; d3++)
4869 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4872 /* Check of we are not at the box edge.
4873 * pbc is only handled in the first step above,
4874 * but this check could move over pbc while
4875 * the first step did not due to different rounding.
4877 if (pos_d >= cell_x1[dim2] &&
4878 dd->ci[dim2] != dd->nc[dim2]-1)
4880 flag |= DD_FLAG_FW(d2);
4882 else if (pos_d < cell_x0[dim2] &&
4885 flag |= DD_FLAG_BW(d2);
4887 comm->buf_int[cg*DD_CGIBS+1] = flag;
4890 /* Set to which neighboring cell this cg should go */
4891 if (flag & DD_FLAG_FW(d2))
4895 else if (flag & DD_FLAG_BW(d2))
4897 if (dd->nc[dd->dim[d2]] > 2)
4909 nrcg = flag & DD_FLAG_NRCG;
4912 if (home_pos_cg+1 > dd->cg_nalloc)
4914 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4915 srenew(dd->index_gl,dd->cg_nalloc);
4916 srenew(dd->cgindex,dd->cg_nalloc+1);
4918 /* Set the global charge group index and size */
4919 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4920 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4921 /* Copy the state from the buffer */
4922 dd_check_alloc_ncg(fr,state,f,home_pos_cg+1);
4923 if (fr->cutoff_scheme == ecutsGROUP)
4926 copy_rvec(comm->vbuf.v[buf_pos],cg_cm[home_pos_cg]);
4930 /* Set the cginfo */
4931 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4932 dd->index_gl[home_pos_cg]);
4935 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4938 if (home_pos_at+nrcg > state->nalloc)
4940 dd_realloc_state(state,f,home_pos_at+nrcg);
4942 for(i=0; i<nrcg; i++)
4944 copy_rvec(comm->vbuf.v[buf_pos++],
4945 state->x[home_pos_at+i]);
4949 for(i=0; i<nrcg; i++)
4951 copy_rvec(comm->vbuf.v[buf_pos++],
4952 state->v[home_pos_at+i]);
4957 for(i=0; i<nrcg; i++)
4959 copy_rvec(comm->vbuf.v[buf_pos++],
4960 state->sd_X[home_pos_at+i]);
4965 for(i=0; i<nrcg; i++)
4967 copy_rvec(comm->vbuf.v[buf_pos++],
4968 state->cg_p[home_pos_at+i]);
4972 home_pos_at += nrcg;
4976 /* Reallocate the buffers if necessary */
4977 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4979 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4980 srenew(comm->cggl_flag[mc],comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4982 nvr = ncg[mc] + nat[mc]*nvec;
4983 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
4985 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
4986 srenew(comm->cgcm_state[mc],comm->cgcm_state_nalloc[mc]);
4988 /* Copy from the receive to the send buffers */
4989 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
4990 comm->buf_int + cg*DD_CGIBS,
4991 DD_CGIBS*sizeof(int));
4992 memcpy(comm->cgcm_state[mc][nvr],
4993 comm->vbuf.v[buf_pos],
4994 (1+nrcg*nvec)*sizeof(rvec));
4995 buf_pos += 1 + nrcg*nvec;
5002 /* With sorting (!bCompact) the indices are now only partially up to date
5003 * and ncg_home and nat_home are not the real count, since there are
5004 * "holes" in the arrays for the charge groups that moved to neighbors.
5006 if (fr->cutoff_scheme == ecutsVERLET)
5008 moved = get_moved(comm,home_pos_cg);
5010 for(i=dd->ncg_home; i<home_pos_cg; i++)
5015 dd->ncg_home = home_pos_cg;
5016 dd->nat_home = home_pos_at;
5021 "Finished repartitioning: cgs moved out %d, new home %d\n",
5022 *ncg_moved,dd->ncg_home-*ncg_moved);
5027 void dd_cycles_add(gmx_domdec_t *dd,float cycles,int ddCycl)
5029 dd->comm->cycl[ddCycl] += cycles;
5030 dd->comm->cycl_n[ddCycl]++;
5031 if (cycles > dd->comm->cycl_max[ddCycl])
5033 dd->comm->cycl_max[ddCycl] = cycles;
5037 static double force_flop_count(t_nrnb *nrnb)
5044 for(i=0; i<eNR_NBKERNEL_FREE_ENERGY; i++)
5046 /* To get closer to the real timings, we half the count
5047 * for the normal loops and again half it for water loops.
5050 if (strstr(name,"W3") != NULL || strstr(name,"W4") != NULL)
5052 sum += nrnb->n[i]*0.25*cost_nrnb(i);
5056 sum += nrnb->n[i]*0.50*cost_nrnb(i);
5059 for(i=eNR_NBKERNEL_FREE_ENERGY; i<=eNR_NB14; i++)
5062 if (strstr(name,"W3") != NULL || strstr(name,"W4") != NULL)
5063 sum += nrnb->n[i]*cost_nrnb(i);
5065 for(i=eNR_BONDS; i<=eNR_WALLS; i++)
5067 sum += nrnb->n[i]*cost_nrnb(i);
5073 void dd_force_flop_start(gmx_domdec_t *dd,t_nrnb *nrnb)
5075 if (dd->comm->eFlop)
5077 dd->comm->flop -= force_flop_count(nrnb);
5080 void dd_force_flop_stop(gmx_domdec_t *dd,t_nrnb *nrnb)
5082 if (dd->comm->eFlop)
5084 dd->comm->flop += force_flop_count(nrnb);
5089 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
5093 for(i=0; i<ddCyclNr; i++)
5095 dd->comm->cycl[i] = 0;
5096 dd->comm->cycl_n[i] = 0;
5097 dd->comm->cycl_max[i] = 0;
5100 dd->comm->flop_n = 0;
5103 static void get_load_distribution(gmx_domdec_t *dd,gmx_wallcycle_t wcycle)
5105 gmx_domdec_comm_t *comm;
5106 gmx_domdec_load_t *load;
5107 gmx_domdec_root_t *root=NULL;
5108 int d,dim,cid,i,pos;
5109 float cell_frac=0,sbuf[DD_NLOAD_MAX];
5114 fprintf(debug,"get_load_distribution start\n");
5117 wallcycle_start(wcycle,ewcDDCOMMLOAD);
5121 bSepPME = (dd->pme_nodeid >= 0);
5123 for(d=dd->ndim-1; d>=0; d--)
5126 /* Check if we participate in the communication in this dimension */
5127 if (d == dd->ndim-1 ||
5128 (dd->ci[dd->dim[d+1]]==0 && dd->ci[dd->dim[dd->ndim-1]]==0))
5130 load = &comm->load[d];
5133 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
5136 if (d == dd->ndim-1)
5138 sbuf[pos++] = dd_force_load(comm);
5139 sbuf[pos++] = sbuf[0];
5142 sbuf[pos++] = sbuf[0];
5143 sbuf[pos++] = cell_frac;
5146 sbuf[pos++] = comm->cell_f_max0[d];
5147 sbuf[pos++] = comm->cell_f_min1[d];
5152 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
5153 sbuf[pos++] = comm->cycl[ddCyclPME];
5158 sbuf[pos++] = comm->load[d+1].sum;
5159 sbuf[pos++] = comm->load[d+1].max;
5162 sbuf[pos++] = comm->load[d+1].sum_m;
5163 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
5164 sbuf[pos++] = comm->load[d+1].flags;
5167 sbuf[pos++] = comm->cell_f_max0[d];
5168 sbuf[pos++] = comm->cell_f_min1[d];
5173 sbuf[pos++] = comm->load[d+1].mdf;
5174 sbuf[pos++] = comm->load[d+1].pme;
5178 /* Communicate a row in DD direction d.
5179 * The communicators are setup such that the root always has rank 0.
5182 MPI_Gather(sbuf ,load->nload*sizeof(float),MPI_BYTE,
5183 load->load,load->nload*sizeof(float),MPI_BYTE,
5184 0,comm->mpi_comm_load[d]);
5186 if (dd->ci[dim] == dd->master_ci[dim])
5188 /* We are the root, process this row */
5189 if (comm->bDynLoadBal)
5191 root = comm->root[d];
5201 for(i=0; i<dd->nc[dim]; i++)
5203 load->sum += load->load[pos++];
5204 load->max = max(load->max,load->load[pos]);
5210 /* This direction could not be load balanced properly,
5211 * therefore we need to use the maximum iso the average load.
5213 load->sum_m = max(load->sum_m,load->load[pos]);
5217 load->sum_m += load->load[pos];
5220 load->cvol_min = min(load->cvol_min,load->load[pos]);
5224 load->flags = (int)(load->load[pos++] + 0.5);
5228 root->cell_f_max0[i] = load->load[pos++];
5229 root->cell_f_min1[i] = load->load[pos++];
5234 load->mdf = max(load->mdf,load->load[pos]);
5236 load->pme = max(load->pme,load->load[pos]);
5240 if (comm->bDynLoadBal && root->bLimited)
5242 load->sum_m *= dd->nc[dim];
5243 load->flags |= (1<<d);
5251 comm->nload += dd_load_count(comm);
5252 comm->load_step += comm->cycl[ddCyclStep];
5253 comm->load_sum += comm->load[0].sum;
5254 comm->load_max += comm->load[0].max;
5255 if (comm->bDynLoadBal)
5257 for(d=0; d<dd->ndim; d++)
5259 if (comm->load[0].flags & (1<<d))
5261 comm->load_lim[d]++;
5267 comm->load_mdf += comm->load[0].mdf;
5268 comm->load_pme += comm->load[0].pme;
5272 wallcycle_stop(wcycle,ewcDDCOMMLOAD);
5276 fprintf(debug,"get_load_distribution finished\n");
5280 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5282 /* Return the relative performance loss on the total run time
5283 * due to the force calculation load imbalance.
5285 if (dd->comm->nload > 0)
5288 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5289 (dd->comm->load_step*dd->nnodes);
5297 static void print_dd_load_av(FILE *fplog,gmx_domdec_t *dd)
5300 int npp,npme,nnodes,d,limp;
5301 float imbal,pme_f_ratio,lossf,lossp=0;
5303 gmx_domdec_comm_t *comm;
5306 if (DDMASTER(dd) && comm->nload > 0)
5309 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5310 nnodes = npp + npme;
5311 imbal = comm->load_max*npp/comm->load_sum - 1;
5312 lossf = dd_force_imb_perf_loss(dd);
5313 sprintf(buf," Average load imbalance: %.1f %%\n",imbal*100);
5314 fprintf(fplog,"%s",buf);
5315 fprintf(stderr,"\n");
5316 fprintf(stderr,"%s",buf);
5317 sprintf(buf," Part of the total run time spent waiting due to load imbalance: %.1f %%\n",lossf*100);
5318 fprintf(fplog,"%s",buf);
5319 fprintf(stderr,"%s",buf);
5321 if (comm->bDynLoadBal)
5323 sprintf(buf," Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5324 for(d=0; d<dd->ndim; d++)
5326 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5327 sprintf(buf+strlen(buf)," %c %d %%",dim2char(dd->dim[d]),limp);
5333 sprintf(buf+strlen(buf),"\n");
5334 fprintf(fplog,"%s",buf);
5335 fprintf(stderr,"%s",buf);
5339 pme_f_ratio = comm->load_pme/comm->load_mdf;
5340 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5343 lossp *= (float)npme/(float)nnodes;
5347 lossp *= (float)npp/(float)nnodes;
5349 sprintf(buf," Average PME mesh/force load: %5.3f\n",pme_f_ratio);
5350 fprintf(fplog,"%s",buf);
5351 fprintf(stderr,"%s",buf);
5352 sprintf(buf," Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n",fabs(lossp)*100);
5353 fprintf(fplog,"%s",buf);
5354 fprintf(stderr,"%s",buf);
5356 fprintf(fplog,"\n");
5357 fprintf(stderr,"\n");
5359 if (lossf >= DD_PERF_LOSS)
5362 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
5363 " in the domain decomposition.\n",lossf*100);
5364 if (!comm->bDynLoadBal)
5366 sprintf(buf+strlen(buf)," You might want to use dynamic load balancing (option -dlb.)\n");
5370 sprintf(buf+strlen(buf)," You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5372 fprintf(fplog,"%s\n",buf);
5373 fprintf(stderr,"%s\n",buf);
5375 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5378 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5379 " had %s work to do than the PP nodes.\n"
5380 " You might want to %s the number of PME nodes\n"
5381 " or %s the cut-off and the grid spacing.\n",
5383 (lossp < 0) ? "less" : "more",
5384 (lossp < 0) ? "decrease" : "increase",
5385 (lossp < 0) ? "decrease" : "increase");
5386 fprintf(fplog,"%s\n",buf);
5387 fprintf(stderr,"%s\n",buf);
5392 static float dd_vol_min(gmx_domdec_t *dd)
5394 return dd->comm->load[0].cvol_min*dd->nnodes;
5397 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5399 return dd->comm->load[0].flags;
5402 static float dd_f_imbal(gmx_domdec_t *dd)
5404 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5407 float dd_pme_f_ratio(gmx_domdec_t *dd)
5409 if (dd->comm->cycl_n[ddCyclPME] > 0)
5411 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5419 static void dd_print_load(FILE *fplog,gmx_domdec_t *dd,gmx_large_int_t step)
5424 flags = dd_load_flags(dd);
5428 "DD load balancing is limited by minimum cell size in dimension");
5429 for(d=0; d<dd->ndim; d++)
5433 fprintf(fplog," %c",dim2char(dd->dim[d]));
5436 fprintf(fplog,"\n");
5438 fprintf(fplog,"DD step %s",gmx_step_str(step,buf));
5439 if (dd->comm->bDynLoadBal)
5441 fprintf(fplog," vol min/aver %5.3f%c",
5442 dd_vol_min(dd),flags ? '!' : ' ');
5444 fprintf(fplog," load imb.: force %4.1f%%",dd_f_imbal(dd)*100);
5445 if (dd->comm->cycl_n[ddCyclPME])
5447 fprintf(fplog," pme mesh/force %5.3f",dd_pme_f_ratio(dd));
5449 fprintf(fplog,"\n\n");
5452 static void dd_print_load_verbose(gmx_domdec_t *dd)
5454 if (dd->comm->bDynLoadBal)
5456 fprintf(stderr,"vol %4.2f%c ",
5457 dd_vol_min(dd),dd_load_flags(dd) ? '!' : ' ');
5459 fprintf(stderr,"imb F %2d%% ",(int)(dd_f_imbal(dd)*100+0.5));
5460 if (dd->comm->cycl_n[ddCyclPME])
5462 fprintf(stderr,"pme/F %4.2f ",dd_pme_f_ratio(dd));
5467 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind,ivec loc)
5472 gmx_domdec_root_t *root;
5473 gmx_bool bPartOfGroup = FALSE;
5475 dim = dd->dim[dim_ind];
5476 copy_ivec(loc,loc_c);
5477 for(i=0; i<dd->nc[dim]; i++)
5480 rank = dd_index(dd->nc,loc_c);
5481 if (rank == dd->rank)
5483 /* This process is part of the group */
5484 bPartOfGroup = TRUE;
5487 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup?0:MPI_UNDEFINED, dd->rank,
5491 dd->comm->mpi_comm_load[dim_ind] = c_row;
5492 if (dd->comm->eDLB != edlbNO)
5494 if (dd->ci[dim] == dd->master_ci[dim])
5496 /* This is the root process of this row */
5497 snew(dd->comm->root[dim_ind],1);
5498 root = dd->comm->root[dim_ind];
5499 snew(root->cell_f,DD_CELL_F_SIZE(dd,dim_ind));
5500 snew(root->old_cell_f,dd->nc[dim]+1);
5501 snew(root->bCellMin,dd->nc[dim]);
5504 snew(root->cell_f_max0,dd->nc[dim]);
5505 snew(root->cell_f_min1,dd->nc[dim]);
5506 snew(root->bound_min,dd->nc[dim]);
5507 snew(root->bound_max,dd->nc[dim]);
5509 snew(root->buf_ncd,dd->nc[dim]);
5513 /* This is not a root process, we only need to receive cell_f */
5514 snew(dd->comm->cell_f_row,DD_CELL_F_SIZE(dd,dim_ind));
5517 if (dd->ci[dim] == dd->master_ci[dim])
5519 snew(dd->comm->load[dim_ind].load,dd->nc[dim]*DD_NLOAD_MAX);
5525 static void make_load_communicators(gmx_domdec_t *dd)
5532 fprintf(debug,"Making load communicators\n");
5534 snew(dd->comm->load,dd->ndim);
5535 snew(dd->comm->mpi_comm_load,dd->ndim);
5538 make_load_communicator(dd,0,loc);
5541 for(i=0; i<dd->nc[dim0]; i++) {
5543 make_load_communicator(dd,1,loc);
5548 for(i=0; i<dd->nc[dim0]; i++) {
5551 for(j=0; j<dd->nc[dim1]; j++) {
5553 make_load_communicator(dd,2,loc);
5559 fprintf(debug,"Finished making load communicators\n");
5563 void setup_dd_grid(FILE *fplog,gmx_domdec_t *dd)
5569 ivec dd_zp[DD_MAXIZONE];
5570 gmx_domdec_zones_t *zones;
5571 gmx_domdec_ns_ranges_t *izone;
5573 for(d=0; d<dd->ndim; d++)
5576 copy_ivec(dd->ci,tmp);
5577 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5578 dd->neighbor[d][0] = ddcoord2ddnodeid(dd,tmp);
5579 copy_ivec(dd->ci,tmp);
5580 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5581 dd->neighbor[d][1] = ddcoord2ddnodeid(dd,tmp);
5584 fprintf(debug,"DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5587 dd->neighbor[d][1]);
5593 fprintf(fplog,"\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5595 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],
5596 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5603 for(i=0; i<nzonep; i++)
5605 copy_ivec(dd_zp3[i],dd_zp[i]);
5611 for(i=0; i<nzonep; i++)
5613 copy_ivec(dd_zp2[i],dd_zp[i]);
5619 for(i=0; i<nzonep; i++)
5621 copy_ivec(dd_zp1[i],dd_zp[i]);
5625 gmx_fatal(FARGS,"Can only do 1, 2 or 3D domain decomposition");
5630 zones = &dd->comm->zones;
5632 for(i=0; i<nzone; i++)
5635 clear_ivec(zones->shift[i]);
5636 for(d=0; d<dd->ndim; d++)
5638 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5643 for(i=0; i<nzone; i++)
5645 for(d=0; d<DIM; d++)
5647 s[d] = dd->ci[d] - zones->shift[i][d];
5652 else if (s[d] >= dd->nc[d])
5658 zones->nizone = nzonep;
5659 for(i=0; i<zones->nizone; i++)
5661 if (dd_zp[i][0] != i)
5663 gmx_fatal(FARGS,"Internal inconsistency in the dd grid setup");
5665 izone = &zones->izone[i];
5666 izone->j0 = dd_zp[i][1];
5667 izone->j1 = dd_zp[i][2];
5668 for(dim=0; dim<DIM; dim++)
5670 if (dd->nc[dim] == 1)
5672 /* All shifts should be allowed */
5673 izone->shift0[dim] = -1;
5674 izone->shift1[dim] = 1;
5679 izone->shift0[d] = 0;
5680 izone->shift1[d] = 0;
5681 for(j=izone->j0; j<izone->j1; j++) {
5682 if (dd->shift[j][d] > dd->shift[i][d])
5683 izone->shift0[d] = -1;
5684 if (dd->shift[j][d] < dd->shift[i][d])
5685 izone->shift1[d] = 1;
5691 /* Assume the shift are not more than 1 cell */
5692 izone->shift0[dim] = 1;
5693 izone->shift1[dim] = -1;
5694 for(j=izone->j0; j<izone->j1; j++)
5696 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5697 if (shift_diff < izone->shift0[dim])
5699 izone->shift0[dim] = shift_diff;
5701 if (shift_diff > izone->shift1[dim])
5703 izone->shift1[dim] = shift_diff;
5710 if (dd->comm->eDLB != edlbNO)
5712 snew(dd->comm->root,dd->ndim);
5715 if (dd->comm->bRecordLoad)
5717 make_load_communicators(dd);
5721 static void make_pp_communicator(FILE *fplog,t_commrec *cr,int reorder)
5724 gmx_domdec_comm_t *comm;
5735 if (comm->bCartesianPP)
5737 /* Set up cartesian communication for the particle-particle part */
5740 fprintf(fplog,"Will use a Cartesian communicator: %d x %d x %d\n",
5741 dd->nc[XX],dd->nc[YY],dd->nc[ZZ]);
5744 for(i=0; i<DIM; i++)
5748 MPI_Cart_create(cr->mpi_comm_mygroup,DIM,dd->nc,periods,reorder,
5750 /* We overwrite the old communicator with the new cartesian one */
5751 cr->mpi_comm_mygroup = comm_cart;
5754 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5755 MPI_Comm_rank(dd->mpi_comm_all,&dd->rank);
5757 if (comm->bCartesianPP_PME)
5759 /* Since we want to use the original cartesian setup for sim,
5760 * and not the one after split, we need to make an index.
5762 snew(comm->ddindex2ddnodeid,dd->nnodes);
5763 comm->ddindex2ddnodeid[dd_index(dd->nc,dd->ci)] = dd->rank;
5764 gmx_sumi(dd->nnodes,comm->ddindex2ddnodeid,cr);
5765 /* Get the rank of the DD master,
5766 * above we made sure that the master node is a PP node.
5776 MPI_Allreduce(&rank,&dd->masterrank,1,MPI_INT,MPI_SUM,dd->mpi_comm_all);
5778 else if (comm->bCartesianPP)
5780 if (cr->npmenodes == 0)
5782 /* The PP communicator is also
5783 * the communicator for this simulation
5785 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5787 cr->nodeid = dd->rank;
5789 MPI_Cart_coords(dd->mpi_comm_all,dd->rank,DIM,dd->ci);
5791 /* We need to make an index to go from the coordinates
5792 * to the nodeid of this simulation.
5794 snew(comm->ddindex2simnodeid,dd->nnodes);
5795 snew(buf,dd->nnodes);
5796 if (cr->duty & DUTY_PP)
5798 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5800 /* Communicate the ddindex to simulation nodeid index */
5801 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5802 cr->mpi_comm_mysim);
5805 /* Determine the master coordinates and rank.
5806 * The DD master should be the same node as the master of this sim.
5808 for(i=0; i<dd->nnodes; i++)
5810 if (comm->ddindex2simnodeid[i] == 0)
5812 ddindex2xyz(dd->nc,i,dd->master_ci);
5813 MPI_Cart_rank(dd->mpi_comm_all,dd->master_ci,&dd->masterrank);
5818 fprintf(debug,"The master rank is %d\n",dd->masterrank);
5823 /* No Cartesian communicators */
5824 /* We use the rank in dd->comm->all as DD index */
5825 ddindex2xyz(dd->nc,dd->rank,dd->ci);
5826 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5828 clear_ivec(dd->master_ci);
5835 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5836 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5841 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5842 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5846 static void receive_ddindex2simnodeid(t_commrec *cr)
5850 gmx_domdec_comm_t *comm;
5857 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5859 snew(comm->ddindex2simnodeid,dd->nnodes);
5860 snew(buf,dd->nnodes);
5861 if (cr->duty & DUTY_PP)
5863 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5866 /* Communicate the ddindex to simulation nodeid index */
5867 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5868 cr->mpi_comm_mysim);
5875 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5878 gmx_domdec_master_t *ma;
5883 snew(ma->ncg,dd->nnodes);
5884 snew(ma->index,dd->nnodes+1);
5886 snew(ma->nat,dd->nnodes);
5887 snew(ma->ibuf,dd->nnodes*2);
5888 snew(ma->cell_x,DIM);
5889 for(i=0; i<DIM; i++)
5891 snew(ma->cell_x[i],dd->nc[i]+1);
5894 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5900 snew(ma->vbuf,natoms);
5906 static void split_communicator(FILE *fplog,t_commrec *cr,int dd_node_order,
5910 gmx_domdec_comm_t *comm;
5921 if (comm->bCartesianPP)
5923 for(i=1; i<DIM; i++)
5925 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5927 if (bDiv[YY] || bDiv[ZZ])
5929 comm->bCartesianPP_PME = TRUE;
5930 /* If we have 2D PME decomposition, which is always in x+y,
5931 * we stack the PME only nodes in z.
5932 * Otherwise we choose the direction that provides the thinnest slab
5933 * of PME only nodes as this will have the least effect
5934 * on the PP communication.
5935 * But for the PME communication the opposite might be better.
5937 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
5939 dd->nc[YY] > dd->nc[ZZ]))
5941 comm->cartpmedim = ZZ;
5945 comm->cartpmedim = YY;
5947 comm->ntot[comm->cartpmedim]
5948 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
5952 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]);
5954 "Will not use a Cartesian communicator for PP <-> PME\n\n");
5959 if (comm->bCartesianPP_PME)
5963 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]);
5966 for(i=0; i<DIM; i++)
5970 MPI_Cart_create(cr->mpi_comm_mysim,DIM,comm->ntot,periods,reorder,
5973 MPI_Comm_rank(comm_cart,&rank);
5974 if (MASTERNODE(cr) && rank != 0)
5976 gmx_fatal(FARGS,"MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
5979 /* With this assigment we loose the link to the original communicator
5980 * which will usually be MPI_COMM_WORLD, unless have multisim.
5982 cr->mpi_comm_mysim = comm_cart;
5983 cr->sim_nodeid = rank;
5985 MPI_Cart_coords(cr->mpi_comm_mysim,cr->sim_nodeid,DIM,dd->ci);
5989 fprintf(fplog,"Cartesian nodeid %d, coordinates %d %d %d\n\n",
5990 cr->sim_nodeid,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5993 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
5997 if (cr->npmenodes == 0 ||
5998 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
6000 cr->duty = DUTY_PME;
6003 /* Split the sim communicator into PP and PME only nodes */
6004 MPI_Comm_split(cr->mpi_comm_mysim,
6006 dd_index(comm->ntot,dd->ci),
6007 &cr->mpi_comm_mygroup);
6011 switch (dd_node_order)
6016 fprintf(fplog,"Order of the nodes: PP first, PME last\n");
6019 case ddnoINTERLEAVE:
6020 /* Interleave the PP-only and PME-only nodes,
6021 * as on clusters with dual-core machines this will double
6022 * the communication bandwidth of the PME processes
6023 * and thus speed up the PP <-> PME and inter PME communication.
6027 fprintf(fplog,"Interleaving PP and PME nodes\n");
6029 comm->pmenodes = dd_pmenodes(cr);
6034 gmx_fatal(FARGS,"Unknown dd_node_order=%d",dd_node_order);
6037 if (dd_simnode2pmenode(cr,cr->sim_nodeid) == -1)
6039 cr->duty = DUTY_PME;
6046 /* Split the sim communicator into PP and PME only nodes */
6047 MPI_Comm_split(cr->mpi_comm_mysim,
6050 &cr->mpi_comm_mygroup);
6051 MPI_Comm_rank(cr->mpi_comm_mygroup,&cr->nodeid);
6057 fprintf(fplog,"This is a %s only node\n\n",
6058 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
6062 void make_dd_communicators(FILE *fplog,t_commrec *cr,int dd_node_order)
6065 gmx_domdec_comm_t *comm;
6071 copy_ivec(dd->nc,comm->ntot);
6073 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
6074 comm->bCartesianPP_PME = FALSE;
6076 /* Reorder the nodes by default. This might change the MPI ranks.
6077 * Real reordering is only supported on very few architectures,
6078 * Blue Gene is one of them.
6080 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
6082 if (cr->npmenodes > 0)
6084 /* Split the communicator into a PP and PME part */
6085 split_communicator(fplog,cr,dd_node_order,CartReorder);
6086 if (comm->bCartesianPP_PME)
6088 /* We (possibly) reordered the nodes in split_communicator,
6089 * so it is no longer required in make_pp_communicator.
6091 CartReorder = FALSE;
6096 /* All nodes do PP and PME */
6098 /* We do not require separate communicators */
6099 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
6103 if (cr->duty & DUTY_PP)
6105 /* Copy or make a new PP communicator */
6106 make_pp_communicator(fplog,cr,CartReorder);
6110 receive_ddindex2simnodeid(cr);
6113 if (!(cr->duty & DUTY_PME))
6115 /* Set up the commnuication to our PME node */
6116 dd->pme_nodeid = dd_simnode2pmenode(cr,cr->sim_nodeid);
6117 dd->pme_receive_vir_ener = receive_vir_ener(cr);
6120 fprintf(debug,"My pme_nodeid %d receive ener %d\n",
6121 dd->pme_nodeid,dd->pme_receive_vir_ener);
6126 dd->pme_nodeid = -1;
6131 dd->ma = init_gmx_domdec_master_t(dd,
6133 comm->cgs_gl.index[comm->cgs_gl.nr]);
6137 static real *get_slb_frac(FILE *fplog,const char *dir,int nc,const char *size_string)
6144 if (nc > 1 && size_string != NULL)
6148 fprintf(fplog,"Using static load balancing for the %s direction\n",
6153 for (i=0; i<nc; i++)
6156 sscanf(size_string,"%lf%n",&dbl,&n);
6159 gmx_fatal(FARGS,"Incorrect or not enough DD cell size entries for direction %s: '%s'",dir,size_string);
6168 fprintf(fplog,"Relative cell sizes:");
6170 for (i=0; i<nc; i++)
6175 fprintf(fplog," %5.3f",slb_frac[i]);
6180 fprintf(fplog,"\n");
6187 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
6190 gmx_mtop_ilistloop_t iloop;
6194 iloop = gmx_mtop_ilistloop_init(mtop);
6195 while (gmx_mtop_ilistloop_next(iloop,&il,&nmol))
6197 for(ftype=0; ftype<F_NRE; ftype++)
6199 if ((interaction_function[ftype].flags & IF_BOND) &&
6202 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
6210 static int dd_nst_env(FILE *fplog,const char *env_var,int def)
6216 val = getenv(env_var);
6219 if (sscanf(val,"%d",&nst) <= 0)
6225 fprintf(fplog,"Found env.var. %s = %s, using value %d\n",
6233 static void dd_warning(t_commrec *cr,FILE *fplog,const char *warn_string)
6237 fprintf(stderr,"\n%s\n",warn_string);
6241 fprintf(fplog,"\n%s\n",warn_string);
6245 static void check_dd_restrictions(t_commrec *cr,gmx_domdec_t *dd,
6246 t_inputrec *ir,FILE *fplog)
6248 if (ir->ePBC == epbcSCREW &&
6249 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6251 gmx_fatal(FARGS,"With pbc=%s can only do domain decomposition in the x-direction",epbc_names[ir->ePBC]);
6254 if (ir->ns_type == ensSIMPLE)
6256 gmx_fatal(FARGS,"Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6259 if (ir->nstlist == 0)
6261 gmx_fatal(FARGS,"Domain decomposition does not work with nstlist=0");
6264 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6266 dd_warning(cr,fplog,"comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6270 static real average_cellsize_min(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
6275 r = ddbox->box_size[XX];
6276 for(di=0; di<dd->ndim; di++)
6279 /* Check using the initial average cell size */
6280 r = min(r,ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6286 static int check_dlb_support(FILE *fplog,t_commrec *cr,
6287 const char *dlb_opt,gmx_bool bRecordLoad,
6288 unsigned long Flags,t_inputrec *ir)
6296 case 'a': eDLB = edlbAUTO; break;
6297 case 'n': eDLB = edlbNO; break;
6298 case 'y': eDLB = edlbYES; break;
6299 default: gmx_incons("Unknown dlb_opt");
6302 if (Flags & MD_RERUN)
6307 if (!EI_DYNAMICS(ir->eI))
6309 if (eDLB == edlbYES)
6311 sprintf(buf,"NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n",EI(ir->eI));
6312 dd_warning(cr,fplog,buf);
6320 dd_warning(cr,fplog,"NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6325 if (Flags & MD_REPRODUCIBLE)
6332 dd_warning(cr,fplog,"NOTE: reproducibility requested, will not use dynamic load balancing\n");
6336 dd_warning(cr,fplog,"WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6339 gmx_fatal(FARGS,"Death horror: undefined case (%d) for load balancing choice",eDLB);
6347 static void set_dd_dim(FILE *fplog,gmx_domdec_t *dd)
6352 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6354 /* Decomposition order z,y,x */
6357 fprintf(fplog,"Using domain decomposition order z, y, x\n");
6359 for(dim=DIM-1; dim>=0; dim--)
6361 if (dd->nc[dim] > 1)
6363 dd->dim[dd->ndim++] = dim;
6369 /* Decomposition order x,y,z */
6370 for(dim=0; dim<DIM; dim++)
6372 if (dd->nc[dim] > 1)
6374 dd->dim[dd->ndim++] = dim;
6380 static gmx_domdec_comm_t *init_dd_comm()
6382 gmx_domdec_comm_t *comm;
6386 snew(comm->cggl_flag,DIM*2);
6387 snew(comm->cgcm_state,DIM*2);
6388 for(i=0; i<DIM*2; i++)
6390 comm->cggl_flag_nalloc[i] = 0;
6391 comm->cgcm_state_nalloc[i] = 0;
6394 comm->nalloc_int = 0;
6395 comm->buf_int = NULL;
6397 vec_rvec_init(&comm->vbuf);
6399 comm->n_load_have = 0;
6400 comm->n_load_collect = 0;
6402 for(i=0; i<ddnatNR-ddnatZONE; i++)
6404 comm->sum_nat[i] = 0;
6408 comm->load_step = 0;
6411 clear_ivec(comm->load_lim);
6418 gmx_domdec_t *init_domain_decomposition(FILE *fplog,t_commrec *cr,
6419 unsigned long Flags,
6421 real comm_distance_min,real rconstr,
6422 const char *dlb_opt,real dlb_scale,
6423 const char *sizex,const char *sizey,const char *sizez,
6424 gmx_mtop_t *mtop,t_inputrec *ir,
6427 int *npme_x,int *npme_y)
6430 gmx_domdec_comm_t *comm;
6433 real r_2b,r_mb,r_bonded=-1,r_bonded_limit=-1,limit,acs;
6440 "\nInitializing Domain Decomposition on %d nodes\n",cr->nnodes);
6445 dd->comm = init_dd_comm();
6447 snew(comm->cggl_flag,DIM*2);
6448 snew(comm->cgcm_state,DIM*2);
6450 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6451 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6453 dd->bSendRecv2 = dd_nst_env(fplog,"GMX_DD_SENDRECV2",0);
6454 comm->dlb_scale_lim = dd_nst_env(fplog,"GMX_DLB_MAX",10);
6455 comm->eFlop = dd_nst_env(fplog,"GMX_DLB_FLOP",0);
6456 recload = dd_nst_env(fplog,"GMX_DD_LOAD",1);
6457 comm->nstSortCG = dd_nst_env(fplog,"GMX_DD_SORT",1);
6458 comm->nstDDDump = dd_nst_env(fplog,"GMX_DD_DUMP",0);
6459 comm->nstDDDumpGrid = dd_nst_env(fplog,"GMX_DD_DUMP_GRID",0);
6460 comm->DD_debug = dd_nst_env(fplog,"GMX_DD_DEBUG",0);
6462 dd->pme_recv_f_alloc = 0;
6463 dd->pme_recv_f_buf = NULL;
6465 if (dd->bSendRecv2 && fplog)
6467 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");
6473 fprintf(fplog,"Will load balance based on FLOP count\n");
6475 if (comm->eFlop > 1)
6477 srand(1+cr->nodeid);
6479 comm->bRecordLoad = TRUE;
6483 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6487 comm->eDLB = check_dlb_support(fplog,cr,dlb_opt,comm->bRecordLoad,Flags,ir);
6489 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6492 fprintf(fplog,"Dynamic load balancing: %s\n",edlb_names[comm->eDLB]);
6494 dd->bGridJump = comm->bDynLoadBal;
6495 comm->bPMELoadBalDLBLimits = FALSE;
6497 if (comm->nstSortCG)
6501 if (comm->nstSortCG == 1)
6503 fprintf(fplog,"Will sort the charge groups at every domain (re)decomposition\n");
6507 fprintf(fplog,"Will sort the charge groups every %d steps\n",
6517 fprintf(fplog,"Will not sort the charge groups\n");
6521 comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
6523 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6524 if (comm->bInterCGBondeds)
6526 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6530 comm->bInterCGMultiBody = FALSE;
6533 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6534 dd->bInterCGsettles = inter_charge_group_settles(mtop);
6536 if (ir->rlistlong == 0)
6538 /* Set the cut-off to some very large value,
6539 * so we don't need if statements everywhere in the code.
6540 * We use sqrt, since the cut-off is squared in some places.
6542 comm->cutoff = GMX_CUTOFF_INF;
6546 comm->cutoff = ir->rlistlong;
6548 comm->cutoff_mbody = 0;
6550 comm->cellsize_limit = 0;
6551 comm->bBondComm = FALSE;
6553 if (comm->bInterCGBondeds)
6555 if (comm_distance_min > 0)
6557 comm->cutoff_mbody = comm_distance_min;
6558 if (Flags & MD_DDBONDCOMM)
6560 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6564 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6566 r_bonded_limit = comm->cutoff_mbody;
6568 else if (ir->bPeriodicMols)
6570 /* Can not easily determine the required cut-off */
6571 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");
6572 comm->cutoff_mbody = comm->cutoff/2;
6573 r_bonded_limit = comm->cutoff_mbody;
6579 dd_bonded_cg_distance(fplog,dd,mtop,ir,x,box,
6580 Flags & MD_DDBONDCHECK,&r_2b,&r_mb);
6582 gmx_bcast(sizeof(r_2b),&r_2b,cr);
6583 gmx_bcast(sizeof(r_mb),&r_mb,cr);
6585 /* We use an initial margin of 10% for the minimum cell size,
6586 * except when we are just below the non-bonded cut-off.
6588 if (Flags & MD_DDBONDCOMM)
6590 if (max(r_2b,r_mb) > comm->cutoff)
6592 r_bonded = max(r_2b,r_mb);
6593 r_bonded_limit = 1.1*r_bonded;
6594 comm->bBondComm = TRUE;
6599 r_bonded_limit = min(1.1*r_bonded,comm->cutoff);
6601 /* We determine cutoff_mbody later */
6605 /* No special bonded communication,
6606 * simply increase the DD cut-off.
6608 r_bonded_limit = 1.1*max(r_2b,r_mb);
6609 comm->cutoff_mbody = r_bonded_limit;
6610 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6613 comm->cellsize_limit = max(comm->cellsize_limit,r_bonded_limit);
6617 "Minimum cell size due to bonded interactions: %.3f nm\n",
6618 comm->cellsize_limit);
6622 if (dd->bInterCGcons && rconstr <= 0)
6624 /* There is a cell size limit due to the constraints (P-LINCS) */
6625 rconstr = constr_r_max(fplog,mtop,ir);
6629 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6631 if (rconstr > comm->cellsize_limit)
6633 fprintf(fplog,"This distance will limit the DD cell size, you can override this with -rcon\n");
6637 else if (rconstr > 0 && fplog)
6639 /* Here we do not check for dd->bInterCGcons,
6640 * because one can also set a cell size limit for virtual sites only
6641 * and at this point we don't know yet if there are intercg v-sites.
6644 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6647 comm->cellsize_limit = max(comm->cellsize_limit,rconstr);
6649 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6653 copy_ivec(nc,dd->nc);
6654 set_dd_dim(fplog,dd);
6655 set_ddbox_cr(cr,&dd->nc,ir,box,&comm->cgs_gl,x,ddbox);
6657 if (cr->npmenodes == -1)
6661 acs = average_cellsize_min(dd,ddbox);
6662 if (acs < comm->cellsize_limit)
6666 fprintf(fplog,"ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n",acs,comm->cellsize_limit);
6668 gmx_fatal_collective(FARGS,cr,NULL,
6669 "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",
6670 acs,comm->cellsize_limit);
6675 set_ddbox_cr(cr,NULL,ir,box,&comm->cgs_gl,x,ddbox);
6677 /* We need to choose the optimal DD grid and possibly PME nodes */
6678 limit = dd_choose_grid(fplog,cr,dd,ir,mtop,box,ddbox,
6679 comm->eDLB!=edlbNO,dlb_scale,
6680 comm->cellsize_limit,comm->cutoff,
6681 comm->bInterCGBondeds,comm->bInterCGMultiBody);
6683 if (dd->nc[XX] == 0)
6685 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6686 sprintf(buf,"Change the number of nodes or mdrun option %s%s%s",
6687 !bC ? "-rdd" : "-rcon",
6688 comm->eDLB!=edlbNO ? " or -dds" : "",
6689 bC ? " or your LINCS settings" : "");
6691 gmx_fatal_collective(FARGS,cr,NULL,
6692 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6694 "Look in the log file for details on the domain decomposition",
6695 cr->nnodes-cr->npmenodes,limit,buf);
6697 set_dd_dim(fplog,dd);
6703 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6704 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],cr->npmenodes);
6707 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6708 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6710 gmx_fatal_collective(FARGS,cr,NULL,
6711 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6712 dd->nnodes,cr->nnodes - cr->npmenodes,cr->nnodes);
6714 if (cr->npmenodes > dd->nnodes)
6716 gmx_fatal_collective(FARGS,cr,NULL,
6717 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.",cr->npmenodes,dd->nnodes);
6719 if (cr->npmenodes > 0)
6721 comm->npmenodes = cr->npmenodes;
6725 comm->npmenodes = dd->nnodes;
6728 if (EEL_PME(ir->coulombtype))
6730 /* The following choices should match those
6731 * in comm_cost_est in domdec_setup.c.
6732 * Note that here the checks have to take into account
6733 * that the decomposition might occur in a different order than xyz
6734 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6735 * in which case they will not match those in comm_cost_est,
6736 * but since that is mainly for testing purposes that's fine.
6738 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6739 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6740 getenv("GMX_PMEONEDD") == NULL)
6742 comm->npmedecompdim = 2;
6743 comm->npmenodes_x = dd->nc[XX];
6744 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6748 /* In case nc is 1 in both x and y we could still choose to
6749 * decompose pme in y instead of x, but we use x for simplicity.
6751 comm->npmedecompdim = 1;
6752 if (dd->dim[0] == YY)
6754 comm->npmenodes_x = 1;
6755 comm->npmenodes_y = comm->npmenodes;
6759 comm->npmenodes_x = comm->npmenodes;
6760 comm->npmenodes_y = 1;
6765 fprintf(fplog,"PME domain decomposition: %d x %d x %d\n",
6766 comm->npmenodes_x,comm->npmenodes_y,1);
6771 comm->npmedecompdim = 0;
6772 comm->npmenodes_x = 0;
6773 comm->npmenodes_y = 0;
6776 /* Technically we don't need both of these,
6777 * but it simplifies code not having to recalculate it.
6779 *npme_x = comm->npmenodes_x;
6780 *npme_y = comm->npmenodes_y;
6782 snew(comm->slb_frac,DIM);
6783 if (comm->eDLB == edlbNO)
6785 comm->slb_frac[XX] = get_slb_frac(fplog,"x",dd->nc[XX],sizex);
6786 comm->slb_frac[YY] = get_slb_frac(fplog,"y",dd->nc[YY],sizey);
6787 comm->slb_frac[ZZ] = get_slb_frac(fplog,"z",dd->nc[ZZ],sizez);
6790 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6792 if (comm->bBondComm || comm->eDLB != edlbNO)
6794 /* Set the bonded communication distance to halfway
6795 * the minimum and the maximum,
6796 * since the extra communication cost is nearly zero.
6798 acs = average_cellsize_min(dd,ddbox);
6799 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6800 if (comm->eDLB != edlbNO)
6802 /* Check if this does not limit the scaling */
6803 comm->cutoff_mbody = min(comm->cutoff_mbody,dlb_scale*acs);
6805 if (!comm->bBondComm)
6807 /* Without bBondComm do not go beyond the n.b. cut-off */
6808 comm->cutoff_mbody = min(comm->cutoff_mbody,comm->cutoff);
6809 if (comm->cellsize_limit >= comm->cutoff)
6811 /* We don't loose a lot of efficieny
6812 * when increasing it to the n.b. cut-off.
6813 * It can even be slightly faster, because we need
6814 * less checks for the communication setup.
6816 comm->cutoff_mbody = comm->cutoff;
6819 /* Check if we did not end up below our original limit */
6820 comm->cutoff_mbody = max(comm->cutoff_mbody,r_bonded_limit);
6822 if (comm->cutoff_mbody > comm->cellsize_limit)
6824 comm->cellsize_limit = comm->cutoff_mbody;
6827 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6832 fprintf(debug,"Bonded atom communication beyond the cut-off: %d\n"
6833 "cellsize limit %f\n",
6834 comm->bBondComm,comm->cellsize_limit);
6839 check_dd_restrictions(cr,dd,ir,fplog);
6842 comm->partition_step = INT_MIN;
6845 clear_dd_cycle_counts(dd);
6850 static void set_dlb_limits(gmx_domdec_t *dd)
6855 for(d=0; d<dd->ndim; d++)
6857 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6858 dd->comm->cellsize_min[dd->dim[d]] =
6859 dd->comm->cellsize_min_dlb[dd->dim[d]];
6864 static void turn_on_dlb(FILE *fplog,t_commrec *cr,gmx_large_int_t step)
6867 gmx_domdec_comm_t *comm;
6877 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);
6880 cellsize_min = comm->cellsize_min[dd->dim[0]];
6881 for(d=1; d<dd->ndim; d++)
6883 cellsize_min = min(cellsize_min,comm->cellsize_min[dd->dim[d]]);
6886 if (cellsize_min < comm->cellsize_limit*1.05)
6888 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");
6890 /* Change DLB from "auto" to "no". */
6891 comm->eDLB = edlbNO;
6896 dd_warning(cr,fplog,"NOTE: Turning on dynamic load balancing\n");
6897 comm->bDynLoadBal = TRUE;
6898 dd->bGridJump = TRUE;
6902 /* We can set the required cell size info here,
6903 * so we do not need to communicate this.
6904 * The grid is completely uniform.
6906 for(d=0; d<dd->ndim; d++)
6910 comm->load[d].sum_m = comm->load[d].sum;
6912 nc = dd->nc[dd->dim[d]];
6915 comm->root[d]->cell_f[i] = i/(real)nc;
6918 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6919 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6922 comm->root[d]->cell_f[nc] = 1.0;
6927 static char *init_bLocalCG(gmx_mtop_t *mtop)
6932 ncg = ncg_mtop(mtop);
6934 for(cg=0; cg<ncg; cg++)
6936 bLocalCG[cg] = FALSE;
6942 void dd_init_bondeds(FILE *fplog,
6943 gmx_domdec_t *dd,gmx_mtop_t *mtop,
6944 gmx_vsite_t *vsite,gmx_constr_t constr,
6945 t_inputrec *ir,gmx_bool bBCheck,cginfo_mb_t *cginfo_mb)
6947 gmx_domdec_comm_t *comm;
6951 dd_make_reverse_top(fplog,dd,mtop,vsite,constr,ir,bBCheck);
6955 if (comm->bBondComm)
6957 /* Communicate atoms beyond the cut-off for bonded interactions */
6960 comm->cglink = make_charge_group_links(mtop,dd,cginfo_mb);
6962 comm->bLocalCG = init_bLocalCG(mtop);
6966 /* Only communicate atoms based on cut-off */
6967 comm->cglink = NULL;
6968 comm->bLocalCG = NULL;
6972 static void print_dd_settings(FILE *fplog,gmx_domdec_t *dd,
6974 gmx_bool bDynLoadBal,real dlb_scale,
6977 gmx_domdec_comm_t *comm;
6992 fprintf(fplog,"The maximum number of communication pulses is:");
6993 for(d=0; d<dd->ndim; d++)
6995 fprintf(fplog," %c %d",dim2char(dd->dim[d]),comm->cd[d].np_dlb);
6997 fprintf(fplog,"\n");
6998 fprintf(fplog,"The minimum size for domain decomposition cells is %.3f nm\n",comm->cellsize_limit);
6999 fprintf(fplog,"The requested allowed shrink of DD cells (option -dds) is: %.2f\n",dlb_scale);
7000 fprintf(fplog,"The allowed shrink of domain decomposition cells is:");
7001 for(d=0; d<DIM; d++)
7005 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
7012 comm->cellsize_min_dlb[d]/
7013 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
7015 fprintf(fplog," %c %.2f",dim2char(d),shrink);
7018 fprintf(fplog,"\n");
7022 set_dd_cell_sizes_slb(dd,ddbox,FALSE,np);
7023 fprintf(fplog,"The initial number of communication pulses is:");
7024 for(d=0; d<dd->ndim; d++)
7026 fprintf(fplog," %c %d",dim2char(dd->dim[d]),np[dd->dim[d]]);
7028 fprintf(fplog,"\n");
7029 fprintf(fplog,"The initial domain decomposition cell size is:");
7030 for(d=0; d<DIM; d++) {
7033 fprintf(fplog," %c %.2f nm",
7034 dim2char(d),dd->comm->cellsize_min[d]);
7037 fprintf(fplog,"\n\n");
7040 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
7042 fprintf(fplog,"The maximum allowed distance for charge groups involved in interactions is:\n");
7043 fprintf(fplog,"%40s %-7s %6.3f nm\n",
7044 "non-bonded interactions","",comm->cutoff);
7048 limit = dd->comm->cellsize_limit;
7052 if (dynamic_dd_box(ddbox,ir))
7054 fprintf(fplog,"(the following are initial values, they could change due to box deformation)\n");
7056 limit = dd->comm->cellsize_min[XX];
7057 for(d=1; d<DIM; d++)
7059 limit = min(limit,dd->comm->cellsize_min[d]);
7063 if (comm->bInterCGBondeds)
7065 fprintf(fplog,"%40s %-7s %6.3f nm\n",
7066 "two-body bonded interactions","(-rdd)",
7067 max(comm->cutoff,comm->cutoff_mbody));
7068 fprintf(fplog,"%40s %-7s %6.3f nm\n",
7069 "multi-body bonded interactions","(-rdd)",
7070 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff,limit));
7074 fprintf(fplog,"%40s %-7s %6.3f nm\n",
7075 "virtual site constructions","(-rcon)",limit);
7077 if (dd->constraint_comm)
7079 sprintf(buf,"atoms separated by up to %d constraints",
7081 fprintf(fplog,"%40s %-7s %6.3f nm\n",
7082 buf,"(-rcon)",limit);
7084 fprintf(fplog,"\n");
7090 static void set_cell_limits_dlb(gmx_domdec_t *dd,
7092 const t_inputrec *ir,
7093 const gmx_ddbox_t *ddbox)
7095 gmx_domdec_comm_t *comm;
7096 int d,dim,npulse,npulse_d_max,npulse_d;
7101 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
7103 /* Determine the maximum number of comm. pulses in one dimension */
7105 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
7107 /* Determine the maximum required number of grid pulses */
7108 if (comm->cellsize_limit >= comm->cutoff)
7110 /* Only a single pulse is required */
7113 else if (!bNoCutOff && comm->cellsize_limit > 0)
7115 /* We round down slightly here to avoid overhead due to the latency
7116 * of extra communication calls when the cut-off
7117 * would be only slightly longer than the cell size.
7118 * Later cellsize_limit is redetermined,
7119 * so we can not miss interactions due to this rounding.
7121 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
7125 /* There is no cell size limit */
7126 npulse = max(dd->nc[XX]-1,max(dd->nc[YY]-1,dd->nc[ZZ]-1));
7129 if (!bNoCutOff && npulse > 1)
7131 /* See if we can do with less pulses, based on dlb_scale */
7133 for(d=0; d<dd->ndim; d++)
7136 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
7137 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
7138 npulse_d_max = max(npulse_d_max,npulse_d);
7140 npulse = min(npulse,npulse_d_max);
7143 /* This env var can override npulse */
7144 d = dd_nst_env(debug,"GMX_DD_NPULSE",0);
7151 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
7152 for(d=0; d<dd->ndim; d++)
7154 comm->cd[d].np_dlb = min(npulse,dd->nc[dd->dim[d]]-1);
7155 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
7156 snew(comm->cd[d].ind,comm->cd[d].np_nalloc);
7157 comm->maxpulse = max(comm->maxpulse,comm->cd[d].np_dlb);
7158 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
7160 comm->bVacDLBNoLimit = FALSE;
7164 /* cellsize_limit is set for LINCS in init_domain_decomposition */
7165 if (!comm->bVacDLBNoLimit)
7167 comm->cellsize_limit = max(comm->cellsize_limit,
7168 comm->cutoff/comm->maxpulse);
7170 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
7171 /* Set the minimum cell size for each DD dimension */
7172 for(d=0; d<dd->ndim; d++)
7174 if (comm->bVacDLBNoLimit ||
7175 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
7177 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
7181 comm->cellsize_min_dlb[dd->dim[d]] =
7182 comm->cutoff/comm->cd[d].np_dlb;
7185 if (comm->cutoff_mbody <= 0)
7187 comm->cutoff_mbody = min(comm->cutoff,comm->cellsize_limit);
7189 if (comm->bDynLoadBal)
7195 gmx_bool dd_bonded_molpbc(gmx_domdec_t *dd,int ePBC)
7197 /* If each molecule is a single charge group
7198 * or we use domain decomposition for each periodic dimension,
7199 * we do not need to take pbc into account for the bonded interactions.
7201 return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
7204 (dd->nc[ZZ]>1 || ePBC==epbcXY)));
7207 void set_dd_parameters(FILE *fplog,gmx_domdec_t *dd,real dlb_scale,
7208 t_inputrec *ir,t_forcerec *fr,
7211 gmx_domdec_comm_t *comm;
7217 /* Initialize the thread data.
7218 * This can not be done in init_domain_decomposition,
7219 * as the numbers of threads is determined later.
7221 comm->nth = gmx_omp_nthreads_get(emntDomdec);
7224 snew(comm->dth,comm->nth);
7227 if (EEL_PME(ir->coulombtype))
7229 init_ddpme(dd,&comm->ddpme[0],0);
7230 if (comm->npmedecompdim >= 2)
7232 init_ddpme(dd,&comm->ddpme[1],1);
7237 comm->npmenodes = 0;
7238 if (dd->pme_nodeid >= 0)
7240 gmx_fatal_collective(FARGS,NULL,dd,
7241 "Can not have separate PME nodes without PME electrostatics");
7247 fprintf(debug,"The DD cut-off is %f\n",comm->cutoff);
7249 if (comm->eDLB != edlbNO)
7251 set_cell_limits_dlb(dd,dlb_scale,ir,ddbox);
7254 print_dd_settings(fplog,dd,ir,comm->bDynLoadBal,dlb_scale,ddbox);
7255 if (comm->eDLB == edlbAUTO)
7259 fprintf(fplog,"When dynamic load balancing gets turned on, these settings will change to:\n");
7261 print_dd_settings(fplog,dd,ir,TRUE,dlb_scale,ddbox);
7264 if (ir->ePBC == epbcNONE)
7266 vol_frac = 1 - 1/(double)dd->nnodes;
7271 (1 + comm_box_frac(dd->nc,comm->cutoff,ddbox))/(double)dd->nnodes;
7275 fprintf(debug,"Volume fraction for all DD zones: %f\n",vol_frac);
7277 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7279 dd->ga2la = ga2la_init(natoms_tot,vol_frac*natoms_tot);
7282 static gmx_bool test_dd_cutoff(t_commrec *cr,
7283 t_state *state,t_inputrec *ir,
7294 set_ddbox(dd,FALSE,cr,ir,state->box,
7295 TRUE,&dd->comm->cgs_gl,state->x,&ddbox);
7299 for(d=0; d<dd->ndim; d++)
7303 inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
7304 if (dynamic_dd_box(&ddbox,ir))
7306 inv_cell_size *= DD_PRES_SCALE_MARGIN;
7309 np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
7311 if (dd->comm->eDLB != edlbNO && dim < ddbox.npbcdim &&
7312 dd->comm->cd[d].np_dlb > 0)
7314 if (np > dd->comm->cd[d].np_dlb)
7319 /* If a current local cell size is smaller than the requested
7320 * cut-off, we could still fix it, but this gets very complicated.
7321 * Without fixing here, we might actually need more checks.
7323 if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
7330 if (dd->comm->eDLB != edlbNO)
7332 /* If DLB is not active yet, we don't need to check the grid jumps.
7333 * Actually we shouldn't, because then the grid jump data is not set.
7335 if (dd->comm->bDynLoadBal &&
7336 check_grid_jump(0,dd,cutoff_req,&ddbox,FALSE))
7341 gmx_sumi(1,&LocallyLimited,cr);
7343 if (LocallyLimited > 0)
7352 gmx_bool change_dd_cutoff(t_commrec *cr,t_state *state,t_inputrec *ir,
7355 gmx_bool bCutoffAllowed;
7357 bCutoffAllowed = test_dd_cutoff(cr,state,ir,cutoff_req);
7361 cr->dd->comm->cutoff = cutoff_req;
7364 return bCutoffAllowed;
7367 void change_dd_dlb_cutoff_limit(t_commrec *cr)
7369 gmx_domdec_comm_t *comm;
7371 comm = cr->dd->comm;
7373 /* Turn on the DLB limiting (might have been on already) */
7374 comm->bPMELoadBalDLBLimits = TRUE;
7376 /* Change the cut-off limit */
7377 comm->PMELoadBal_max_cutoff = comm->cutoff;
7380 static void merge_cg_buffers(int ncell,
7381 gmx_domdec_comm_dim_t *cd, int pulse,
7383 int *index_gl, int *recv_i,
7384 rvec *cg_cm, rvec *recv_vr,
7386 cginfo_mb_t *cginfo_mb,int *cginfo)
7388 gmx_domdec_ind_t *ind,*ind_p;
7389 int p,cell,c,cg,cg0,cg1,cg_gl,nat;
7392 ind = &cd->ind[pulse];
7394 /* First correct the already stored data */
7395 shift = ind->nrecv[ncell];
7396 for(cell=ncell-1; cell>=0; cell--)
7398 shift -= ind->nrecv[cell];
7401 /* Move the cg's present from previous grid pulses */
7402 cg0 = ncg_cell[ncell+cell];
7403 cg1 = ncg_cell[ncell+cell+1];
7404 cgindex[cg1+shift] = cgindex[cg1];
7405 for(cg=cg1-1; cg>=cg0; cg--)
7407 index_gl[cg+shift] = index_gl[cg];
7408 copy_rvec(cg_cm[cg],cg_cm[cg+shift]);
7409 cgindex[cg+shift] = cgindex[cg];
7410 cginfo[cg+shift] = cginfo[cg];
7412 /* Correct the already stored send indices for the shift */
7413 for(p=1; p<=pulse; p++)
7415 ind_p = &cd->ind[p];
7417 for(c=0; c<cell; c++)
7419 cg0 += ind_p->nsend[c];
7421 cg1 = cg0 + ind_p->nsend[cell];
7422 for(cg=cg0; cg<cg1; cg++)
7424 ind_p->index[cg] += shift;
7430 /* Merge in the communicated buffers */
7434 for(cell=0; cell<ncell; cell++)
7436 cg1 = ncg_cell[ncell+cell+1] + shift;
7439 /* Correct the old cg indices */
7440 for(cg=ncg_cell[ncell+cell]; cg<cg1; cg++)
7442 cgindex[cg+1] += shift_at;
7445 for(cg=0; cg<ind->nrecv[cell]; cg++)
7447 /* Copy this charge group from the buffer */
7448 index_gl[cg1] = recv_i[cg0];
7449 copy_rvec(recv_vr[cg0],cg_cm[cg1]);
7450 /* Add it to the cgindex */
7451 cg_gl = index_gl[cg1];
7452 cginfo[cg1] = ddcginfo(cginfo_mb,cg_gl);
7453 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7454 cgindex[cg1+1] = cgindex[cg1] + nat;
7459 shift += ind->nrecv[cell];
7460 ncg_cell[ncell+cell+1] = cg1;
7464 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7465 int nzone,int cg0,const int *cgindex)
7469 /* Store the atom block boundaries for easy copying of communication buffers
7472 for(zone=0; zone<nzone; zone++)
7474 for(p=0; p<cd->np; p++) {
7475 cd->ind[p].cell2at0[zone] = cgindex[cg];
7476 cg += cd->ind[p].nrecv[zone];
7477 cd->ind[p].cell2at1[zone] = cgindex[cg];
7482 static gmx_bool missing_link(t_blocka *link,int cg_gl,char *bLocalCG)
7488 for(i=link->index[cg_gl]; i<link->index[cg_gl+1]; i++)
7490 if (!bLocalCG[link->a[i]])
7499 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
7501 real c[DIM][4]; /* the corners for the non-bonded communication */
7502 real cr0; /* corner for rounding */
7503 real cr1[4]; /* corners for rounding */
7504 real bc[DIM]; /* corners for bounded communication */
7505 real bcr1; /* corner for rounding for bonded communication */
7508 /* Determine the corners of the domain(s) we are communicating with */
7510 set_dd_corners(const gmx_domdec_t *dd,
7511 int dim0, int dim1, int dim2,
7515 const gmx_domdec_comm_t *comm;
7516 const gmx_domdec_zones_t *zones;
7521 zones = &comm->zones;
7523 /* Keep the compiler happy */
7527 /* The first dimension is equal for all cells */
7528 c->c[0][0] = comm->cell_x0[dim0];
7531 c->bc[0] = c->c[0][0];
7536 /* This cell row is only seen from the first row */
7537 c->c[1][0] = comm->cell_x0[dim1];
7538 /* All rows can see this row */
7539 c->c[1][1] = comm->cell_x0[dim1];
7542 c->c[1][1] = max(comm->cell_x0[dim1],comm->zone_d1[1].mch0);
7545 /* For the multi-body distance we need the maximum */
7546 c->bc[1] = max(comm->cell_x0[dim1],comm->zone_d1[1].p1_0);
7549 /* Set the upper-right corner for rounding */
7550 c->cr0 = comm->cell_x1[dim0];
7557 c->c[2][j] = comm->cell_x0[dim2];
7561 /* Use the maximum of the i-cells that see a j-cell */
7562 for(i=0; i<zones->nizone; i++)
7564 for(j=zones->izone[i].j0; j<zones->izone[i].j1; j++)
7570 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7576 /* For the multi-body distance we need the maximum */
7577 c->bc[2] = comm->cell_x0[dim2];
7582 c->bc[2] = max(c->bc[2],comm->zone_d2[i][j].p1_0);
7588 /* Set the upper-right corner for rounding */
7589 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7590 * Only cell (0,0,0) can see cell 7 (1,1,1)
7592 c->cr1[0] = comm->cell_x1[dim1];
7593 c->cr1[3] = comm->cell_x1[dim1];
7596 c->cr1[0] = max(comm->cell_x1[dim1],comm->zone_d1[1].mch1);
7599 /* For the multi-body distance we need the maximum */
7600 c->bcr1 = max(comm->cell_x1[dim1],comm->zone_d1[1].p1_1);
7607 /* Determine which cg's we need to send in this pulse from this zone */
7609 get_zone_pulse_cgs(gmx_domdec_t *dd,
7610 int zonei, int zone,
7612 const int *index_gl,
7614 int dim, int dim_ind,
7615 int dim0, int dim1, int dim2,
7616 real r_comm2, real r_bcomm2,
7620 real skew_fac2_d, real skew_fac_01,
7621 rvec *v_d, rvec *v_0, rvec *v_1,
7622 const dd_corners_t *c,
7624 gmx_bool bDistBonded,
7630 gmx_domdec_ind_t *ind,
7631 int **ibuf, int *ibuf_nalloc,
7637 gmx_domdec_comm_t *comm;
7639 gmx_bool bDistMB_pulse;
7641 real r2,rb2,r,tric_sh;
7644 int nsend_z,nsend,nat;
7648 bScrew = (dd->bScrewPBC && dim == XX);
7650 bDistMB_pulse = (bDistMB && bDistBonded);
7656 for(cg=cg0; cg<cg1; cg++)
7660 if (tric_dist[dim_ind] == 0)
7662 /* Rectangular direction, easy */
7663 r = cg_cm[cg][dim] - c->c[dim_ind][zone];
7670 r = cg_cm[cg][dim] - c->bc[dim_ind];
7676 /* Rounding gives at most a 16% reduction
7677 * in communicated atoms
7679 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7681 r = cg_cm[cg][dim0] - c->cr0;
7682 /* This is the first dimension, so always r >= 0 */
7689 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7691 r = cg_cm[cg][dim1] - c->cr1[zone];
7698 r = cg_cm[cg][dim1] - c->bcr1;
7708 /* Triclinic direction, more complicated */
7711 /* Rounding, conservative as the skew_fac multiplication
7712 * will slightly underestimate the distance.
7714 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7716 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
7717 for(i=dim0+1; i<DIM; i++)
7719 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7721 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7724 rb[dim0] = rn[dim0];
7727 /* Take care that the cell planes along dim0 might not
7728 * be orthogonal to those along dim1 and dim2.
7730 for(i=1; i<=dim_ind; i++)
7733 if (normal[dim0][dimd] > 0)
7735 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7738 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7743 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7745 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
7747 for(i=dim1+1; i<DIM; i++)
7749 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7751 rn[dim1] += tric_sh;
7754 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7755 /* Take care of coupling of the distances
7756 * to the planes along dim0 and dim1 through dim2.
7758 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7759 /* Take care that the cell planes along dim1
7760 * might not be orthogonal to that along dim2.
7762 if (normal[dim1][dim2] > 0)
7764 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7770 cg_cm[cg][dim1] - c->bcr1 + tric_sh;
7773 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7774 /* Take care of coupling of the distances
7775 * to the planes along dim0 and dim1 through dim2.
7777 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7778 /* Take care that the cell planes along dim1
7779 * might not be orthogonal to that along dim2.
7781 if (normal[dim1][dim2] > 0)
7783 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7788 /* The distance along the communication direction */
7789 rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
7791 for(i=dim+1; i<DIM; i++)
7793 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7798 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7799 /* Take care of coupling of the distances
7800 * to the planes along dim0 and dim1 through dim2.
7802 if (dim_ind == 1 && zonei == 1)
7804 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7810 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
7813 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7814 /* Take care of coupling of the distances
7815 * to the planes along dim0 and dim1 through dim2.
7817 if (dim_ind == 1 && zonei == 1)
7819 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7827 ((bDistMB && rb2 < r_bcomm2) ||
7828 (bDist2B && r2 < r_bcomm2)) &&
7830 (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
7831 missing_link(comm->cglink,index_gl[cg],
7834 /* Make an index to the local charge groups */
7835 if (nsend+1 > ind->nalloc)
7837 ind->nalloc = over_alloc_large(nsend+1);
7838 srenew(ind->index,ind->nalloc);
7840 if (nsend+1 > *ibuf_nalloc)
7842 *ibuf_nalloc = over_alloc_large(nsend+1);
7843 srenew(*ibuf,*ibuf_nalloc);
7845 ind->index[nsend] = cg;
7846 (*ibuf)[nsend] = index_gl[cg];
7848 vec_rvec_check_alloc(vbuf,nsend+1);
7850 if (dd->ci[dim] == 0)
7852 /* Correct cg_cm for pbc */
7853 rvec_add(cg_cm[cg],box[dim],vbuf->v[nsend]);
7856 vbuf->v[nsend][YY] = box[YY][YY] - vbuf->v[nsend][YY];
7857 vbuf->v[nsend][ZZ] = box[ZZ][ZZ] - vbuf->v[nsend][ZZ];
7862 copy_rvec(cg_cm[cg],vbuf->v[nsend]);
7865 nat += cgindex[cg+1] - cgindex[cg];
7871 *nsend_z_ptr = nsend_z;
7874 static void setup_dd_communication(gmx_domdec_t *dd,
7875 matrix box,gmx_ddbox_t *ddbox,
7876 t_forcerec *fr,t_state *state,rvec **f)
7878 int dim_ind,dim,dim0,dim1,dim2,dimd,p,nat_tot;
7879 int nzone,nzone_send,zone,zonei,cg0,cg1;
7880 int c,i,j,cg,cg_gl,nrcg;
7881 int *zone_cg_range,pos_cg,*index_gl,*cgindex,*recv_i;
7882 gmx_domdec_comm_t *comm;
7883 gmx_domdec_zones_t *zones;
7884 gmx_domdec_comm_dim_t *cd;
7885 gmx_domdec_ind_t *ind;
7886 cginfo_mb_t *cginfo_mb;
7887 gmx_bool bBondComm,bDist2B,bDistMB,bDistBonded;
7888 real r_mb,r_comm2,r_scomm2,r_bcomm2,r_0,r_1,r2inc,inv_ncg;
7889 dd_corners_t corners;
7891 rvec *cg_cm,*normal,*v_d,*v_0=NULL,*v_1=NULL,*recv_vr;
7892 real skew_fac2_d,skew_fac_01;
7899 fprintf(debug,"Setting up DD communication\n");
7904 switch (fr->cutoff_scheme)
7913 gmx_incons("unimplemented");
7917 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7919 dim = dd->dim[dim_ind];
7921 /* Check if we need to use triclinic distances */
7922 tric_dist[dim_ind] = 0;
7923 for(i=0; i<=dim_ind; i++)
7925 if (ddbox->tric_dir[dd->dim[i]])
7927 tric_dist[dim_ind] = 1;
7932 bBondComm = comm->bBondComm;
7934 /* Do we need to determine extra distances for multi-body bondeds? */
7935 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
7937 /* Do we need to determine extra distances for only two-body bondeds? */
7938 bDist2B = (bBondComm && !bDistMB);
7940 r_comm2 = sqr(comm->cutoff);
7941 r_bcomm2 = sqr(comm->cutoff_mbody);
7945 fprintf(debug,"bBondComm %d, r_bc %f\n",bBondComm,sqrt(r_bcomm2));
7948 zones = &comm->zones;
7951 dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
7952 dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
7954 set_dd_corners(dd,dim0,dim1,dim2,bDistMB,&corners);
7956 /* Triclinic stuff */
7957 normal = ddbox->normal;
7961 v_0 = ddbox->v[dim0];
7962 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
7964 /* Determine the coupling coefficient for the distances
7965 * to the cell planes along dim0 and dim1 through dim2.
7966 * This is required for correct rounding.
7969 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
7972 fprintf(debug,"\nskew_fac_01 %f\n",skew_fac_01);
7978 v_1 = ddbox->v[dim1];
7981 zone_cg_range = zones->cg_range;
7982 index_gl = dd->index_gl;
7983 cgindex = dd->cgindex;
7984 cginfo_mb = fr->cginfo_mb;
7986 zone_cg_range[0] = 0;
7987 zone_cg_range[1] = dd->ncg_home;
7988 comm->zone_ncg1[0] = dd->ncg_home;
7989 pos_cg = dd->ncg_home;
7991 nat_tot = dd->nat_home;
7993 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7995 dim = dd->dim[dim_ind];
7996 cd = &comm->cd[dim_ind];
7998 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
8000 /* No pbc in this dimension, the first node should not comm. */
8008 v_d = ddbox->v[dim];
8009 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
8011 cd->bInPlace = TRUE;
8012 for(p=0; p<cd->np; p++)
8014 /* Only atoms communicated in the first pulse are used
8015 * for multi-body bonded interactions or for bBondComm.
8017 bDistBonded = ((bDistMB || bDist2B) && p == 0);
8022 for(zone=0; zone<nzone_send; zone++)
8024 if (tric_dist[dim_ind] && dim_ind > 0)
8026 /* Determine slightly more optimized skew_fac's
8028 * This reduces the number of communicated atoms
8029 * by about 10% for 3D DD of rhombic dodecahedra.
8031 for(dimd=0; dimd<dim; dimd++)
8033 sf2_round[dimd] = 1;
8034 if (ddbox->tric_dir[dimd])
8036 for(i=dd->dim[dimd]+1; i<DIM; i++)
8038 /* If we are shifted in dimension i
8039 * and the cell plane is tilted forward
8040 * in dimension i, skip this coupling.
8042 if (!(zones->shift[nzone+zone][i] &&
8043 ddbox->v[dimd][i][dimd] >= 0))
8046 sqr(ddbox->v[dimd][i][dimd]);
8049 sf2_round[dimd] = 1/sf2_round[dimd];
8054 zonei = zone_perm[dim_ind][zone];
8057 /* Here we permutate the zones to obtain a convenient order
8058 * for neighbor searching
8060 cg0 = zone_cg_range[zonei];
8061 cg1 = zone_cg_range[zonei+1];
8065 /* Look only at the cg's received in the previous grid pulse
8067 cg1 = zone_cg_range[nzone+zone+1];
8068 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
8071 #pragma omp parallel for num_threads(comm->nth) schedule(static)
8072 for(th=0; th<comm->nth; th++)
8074 gmx_domdec_ind_t *ind_p;
8075 int **ibuf_p,*ibuf_nalloc_p;
8077 int *nsend_p,*nat_p;
8083 /* Thread 0 writes in the comm buffers */
8085 ibuf_p = &comm->buf_int;
8086 ibuf_nalloc_p = &comm->nalloc_int;
8087 vbuf_p = &comm->vbuf;
8090 nsend_zone_p = &ind->nsend[zone];
8094 /* Other threads write into temp buffers */
8095 ind_p = &comm->dth[th].ind;
8096 ibuf_p = &comm->dth[th].ibuf;
8097 ibuf_nalloc_p = &comm->dth[th].ibuf_nalloc;
8098 vbuf_p = &comm->dth[th].vbuf;
8099 nsend_p = &comm->dth[th].nsend;
8100 nat_p = &comm->dth[th].nat;
8101 nsend_zone_p = &comm->dth[th].nsend_zone;
8103 comm->dth[th].nsend = 0;
8104 comm->dth[th].nat = 0;
8105 comm->dth[th].nsend_zone = 0;
8115 cg0_th = cg0 + ((cg1 - cg0)* th )/comm->nth;
8116 cg1_th = cg0 + ((cg1 - cg0)*(th+1))/comm->nth;
8119 /* Get the cg's for this pulse in this zone */
8120 get_zone_pulse_cgs(dd,zonei,zone,cg0_th,cg1_th,
8122 dim,dim_ind,dim0,dim1,dim2,
8125 normal,skew_fac2_d,skew_fac_01,
8126 v_d,v_0,v_1,&corners,sf2_round,
8127 bDistBonded,bBondComm,
8131 ibuf_p,ibuf_nalloc_p,
8137 /* Append data of threads>=1 to the communication buffers */
8138 for(th=1; th<comm->nth; th++)
8140 dd_comm_setup_work_t *dth;
8143 dth = &comm->dth[th];
8145 ns1 = nsend + dth->nsend_zone;
8146 if (ns1 > ind->nalloc)
8148 ind->nalloc = over_alloc_dd(ns1);
8149 srenew(ind->index,ind->nalloc);
8151 if (ns1 > comm->nalloc_int)
8153 comm->nalloc_int = over_alloc_dd(ns1);
8154 srenew(comm->buf_int,comm->nalloc_int);
8156 if (ns1 > comm->vbuf.nalloc)
8158 comm->vbuf.nalloc = over_alloc_dd(ns1);
8159 srenew(comm->vbuf.v,comm->vbuf.nalloc);
8162 for(i=0; i<dth->nsend_zone; i++)
8164 ind->index[nsend] = dth->ind.index[i];
8165 comm->buf_int[nsend] = dth->ibuf[i];
8166 copy_rvec(dth->vbuf.v[i],
8167 comm->vbuf.v[nsend]);
8171 ind->nsend[zone] += dth->nsend_zone;
8174 /* Clear the counts in case we do not have pbc */
8175 for(zone=nzone_send; zone<nzone; zone++)
8177 ind->nsend[zone] = 0;
8179 ind->nsend[nzone] = nsend;
8180 ind->nsend[nzone+1] = nat;
8181 /* Communicate the number of cg's and atoms to receive */
8182 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8183 ind->nsend, nzone+2,
8184 ind->nrecv, nzone+2);
8186 /* The rvec buffer is also required for atom buffers of size nsend
8187 * in dd_move_x and dd_move_f.
8189 vec_rvec_check_alloc(&comm->vbuf,ind->nsend[nzone+1]);
8193 /* We can receive in place if only the last zone is not empty */
8194 for(zone=0; zone<nzone-1; zone++)
8196 if (ind->nrecv[zone] > 0)
8198 cd->bInPlace = FALSE;
8203 /* The int buffer is only required here for the cg indices */
8204 if (ind->nrecv[nzone] > comm->nalloc_int2)
8206 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
8207 srenew(comm->buf_int2,comm->nalloc_int2);
8209 /* The rvec buffer is also required for atom buffers
8210 * of size nrecv in dd_move_x and dd_move_f.
8212 i = max(cd->ind[0].nrecv[nzone+1],ind->nrecv[nzone+1]);
8213 vec_rvec_check_alloc(&comm->vbuf2,i);
8217 /* Make space for the global cg indices */
8218 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
8219 || dd->cg_nalloc == 0)
8221 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
8222 srenew(index_gl,dd->cg_nalloc);
8223 srenew(cgindex,dd->cg_nalloc+1);
8225 /* Communicate the global cg indices */
8228 recv_i = index_gl + pos_cg;
8232 recv_i = comm->buf_int2;
8234 dd_sendrecv_int(dd, dim_ind, dddirBackward,
8235 comm->buf_int, nsend,
8236 recv_i, ind->nrecv[nzone]);
8238 /* Make space for cg_cm */
8239 dd_check_alloc_ncg(fr,state,f,pos_cg + ind->nrecv[nzone]);
8240 if (fr->cutoff_scheme == ecutsGROUP)
8248 /* Communicate cg_cm */
8251 recv_vr = cg_cm + pos_cg;
8255 recv_vr = comm->vbuf2.v;
8257 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
8258 comm->vbuf.v, nsend,
8259 recv_vr, ind->nrecv[nzone]);
8261 /* Make the charge group index */
8264 zone = (p == 0 ? 0 : nzone - 1);
8265 while (zone < nzone)
8267 for(cg=0; cg<ind->nrecv[zone]; cg++)
8269 cg_gl = index_gl[pos_cg];
8270 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb,cg_gl);
8271 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
8272 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
8275 /* Update the charge group presence,
8276 * so we can use it in the next pass of the loop.
8278 comm->bLocalCG[cg_gl] = TRUE;
8284 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
8287 zone_cg_range[nzone+zone] = pos_cg;
8292 /* This part of the code is never executed with bBondComm. */
8293 merge_cg_buffers(nzone,cd,p,zone_cg_range,
8294 index_gl,recv_i,cg_cm,recv_vr,
8295 cgindex,fr->cginfo_mb,fr->cginfo);
8296 pos_cg += ind->nrecv[nzone];
8298 nat_tot += ind->nrecv[nzone+1];
8302 /* Store the atom block for easy copying of communication buffers */
8303 make_cell2at_index(cd,nzone,zone_cg_range[nzone],cgindex);
8307 dd->index_gl = index_gl;
8308 dd->cgindex = cgindex;
8310 dd->ncg_tot = zone_cg_range[zones->n];
8311 dd->nat_tot = nat_tot;
8312 comm->nat[ddnatHOME] = dd->nat_home;
8313 for(i=ddnatZONE; i<ddnatNR; i++)
8315 comm->nat[i] = dd->nat_tot;
8320 /* We don't need to update cginfo, since that was alrady done above.
8321 * So we pass NULL for the forcerec.
8323 dd_set_cginfo(dd->index_gl,dd->ncg_home,dd->ncg_tot,
8324 NULL,comm->bLocalCG);
8329 fprintf(debug,"Finished setting up DD communication, zones:");
8330 for(c=0; c<zones->n; c++)
8332 fprintf(debug," %d",zones->cg_range[c+1]-zones->cg_range[c]);
8334 fprintf(debug,"\n");
8338 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
8342 for(c=0; c<zones->nizone; c++)
8344 zones->izone[c].cg1 = zones->cg_range[c+1];
8345 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
8346 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
8350 static void set_zones_size(gmx_domdec_t *dd,
8351 matrix box,const gmx_ddbox_t *ddbox,
8352 int zone_start,int zone_end)
8354 gmx_domdec_comm_t *comm;
8355 gmx_domdec_zones_t *zones;
8357 int z,zi,zj0,zj1,d,dim;
8360 real size_j,add_tric;
8365 zones = &comm->zones;
8367 /* Do we need to determine extra distances for multi-body bondeds? */
8368 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
8370 for(z=zone_start; z<zone_end; z++)
8372 /* Copy cell limits to zone limits.
8373 * Valid for non-DD dims and non-shifted dims.
8375 copy_rvec(comm->cell_x0,zones->size[z].x0);
8376 copy_rvec(comm->cell_x1,zones->size[z].x1);
8379 for(d=0; d<dd->ndim; d++)
8383 for(z=0; z<zones->n; z++)
8385 /* With a staggered grid we have different sizes
8386 * for non-shifted dimensions.
8388 if (dd->bGridJump && zones->shift[z][dim] == 0)
8392 zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
8393 zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
8397 zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
8398 zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
8404 rcmbs = comm->cutoff_mbody;
8405 if (ddbox->tric_dir[dim])
8407 rcs /= ddbox->skew_fac[dim];
8408 rcmbs /= ddbox->skew_fac[dim];
8411 /* Set the lower limit for the shifted zone dimensions */
8412 for(z=zone_start; z<zone_end; z++)
8414 if (zones->shift[z][dim] > 0)
8417 if (!dd->bGridJump || d == 0)
8419 zones->size[z].x0[dim] = comm->cell_x1[dim];
8420 zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
8424 /* Here we take the lower limit of the zone from
8425 * the lowest domain of the zone below.
8429 zones->size[z].x0[dim] =
8430 comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
8436 zones->size[z].x0[dim] =
8437 zones->size[zone_perm[2][z-4]].x0[dim];
8441 zones->size[z].x0[dim] =
8442 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
8445 /* A temporary limit, is updated below */
8446 zones->size[z].x1[dim] = zones->size[z].x0[dim];
8450 for(zi=0; zi<zones->nizone; zi++)
8452 if (zones->shift[zi][dim] == 0)
8454 /* This takes the whole zone into account.
8455 * With multiple pulses this will lead
8456 * to a larger zone then strictly necessary.
8458 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8459 zones->size[zi].x1[dim]+rcmbs);
8467 /* Loop over the i-zones to set the upper limit of each
8470 for(zi=0; zi<zones->nizone; zi++)
8472 if (zones->shift[zi][dim] == 0)
8474 for(z=zones->izone[zi].j0; z<zones->izone[zi].j1; z++)
8476 if (zones->shift[z][dim] > 0)
8478 zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
8479 zones->size[zi].x1[dim]+rcs);
8486 for(z=zone_start; z<zone_end; z++)
8488 /* Initialization only required to keep the compiler happy */
8489 rvec corner_min={0,0,0},corner_max={0,0,0},corner;
8492 /* To determine the bounding box for a zone we need to find
8493 * the extreme corners of 4, 2 or 1 corners.
8495 nc = 1 << (ddbox->npbcdim - 1);
8499 /* Set up a zone corner at x=0, ignoring trilinic couplings */
8503 corner[YY] = zones->size[z].x0[YY];
8507 corner[YY] = zones->size[z].x1[YY];
8511 corner[ZZ] = zones->size[z].x0[ZZ];
8515 corner[ZZ] = zones->size[z].x1[ZZ];
8517 if (dd->ndim == 1 && box[ZZ][YY] != 0)
8519 /* With 1D domain decomposition the cg's are not in
8520 * the triclinic box, but triclinic x-y and rectangular y-z.
8521 * Shift y back, so it will later end up at 0.
8523 corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
8525 /* Apply the triclinic couplings */
8526 for(i=YY; i<ddbox->npbcdim; i++)
8530 corner[j] += corner[i]*box[i][j]/box[i][i];
8535 copy_rvec(corner,corner_min);
8536 copy_rvec(corner,corner_max);
8540 for(i=0; i<DIM; i++)
8542 corner_min[i] = min(corner_min[i],corner[i]);
8543 corner_max[i] = max(corner_max[i],corner[i]);
8547 /* Copy the extreme cornes without offset along x */
8548 for(i=0; i<DIM; i++)
8550 zones->size[z].bb_x0[i] = corner_min[i];
8551 zones->size[z].bb_x1[i] = corner_max[i];
8553 /* Add the offset along x */
8554 zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
8555 zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
8558 if (zone_start == 0)
8561 for(dim=0; dim<DIM; dim++)
8563 vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
8565 zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0])/vol;
8570 for(z=zone_start; z<zone_end; z++)
8572 fprintf(debug,"zone %d %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8574 zones->size[z].x0[XX],zones->size[z].x1[XX],
8575 zones->size[z].x0[YY],zones->size[z].x1[YY],
8576 zones->size[z].x0[ZZ],zones->size[z].x1[ZZ]);
8577 fprintf(debug,"zone %d bb %6.3f - %6.3f %6.3f - %6.3f %6.3f - %6.3f\n",
8579 zones->size[z].bb_x0[XX],zones->size[z].bb_x1[XX],
8580 zones->size[z].bb_x0[YY],zones->size[z].bb_x1[YY],
8581 zones->size[z].bb_x0[ZZ],zones->size[z].bb_x1[ZZ]);
8586 static int comp_cgsort(const void *a,const void *b)
8590 gmx_cgsort_t *cga,*cgb;
8591 cga = (gmx_cgsort_t *)a;
8592 cgb = (gmx_cgsort_t *)b;
8594 comp = cga->nsc - cgb->nsc;
8597 comp = cga->ind_gl - cgb->ind_gl;
8603 static void order_int_cg(int n,const gmx_cgsort_t *sort,
8608 /* Order the data */
8611 buf[i] = a[sort[i].ind];
8614 /* Copy back to the original array */
8621 static void order_vec_cg(int n,const gmx_cgsort_t *sort,
8626 /* Order the data */
8629 copy_rvec(v[sort[i].ind],buf[i]);
8632 /* Copy back to the original array */
8635 copy_rvec(buf[i],v[i]);
8639 static void order_vec_atom(int ncg,const int *cgindex,const gmx_cgsort_t *sort,
8642 int a,atot,cg,cg0,cg1,i;
8644 if (cgindex == NULL)
8646 /* Avoid the useless loop of the atoms within a cg */
8647 order_vec_cg(ncg,sort,v,buf);
8652 /* Order the data */
8654 for(cg=0; cg<ncg; cg++)
8656 cg0 = cgindex[sort[cg].ind];
8657 cg1 = cgindex[sort[cg].ind+1];
8658 for(i=cg0; i<cg1; i++)
8660 copy_rvec(v[i],buf[a]);
8666 /* Copy back to the original array */
8667 for(a=0; a<atot; a++)
8669 copy_rvec(buf[a],v[a]);
8673 static void ordered_sort(int nsort2,gmx_cgsort_t *sort2,
8674 int nsort_new,gmx_cgsort_t *sort_new,
8675 gmx_cgsort_t *sort1)
8679 /* The new indices are not very ordered, so we qsort them */
8680 qsort_threadsafe(sort_new,nsort_new,sizeof(sort_new[0]),comp_cgsort);
8682 /* sort2 is already ordered, so now we can merge the two arrays */
8686 while(i2 < nsort2 || i_new < nsort_new)
8690 sort1[i1++] = sort_new[i_new++];
8692 else if (i_new == nsort_new)
8694 sort1[i1++] = sort2[i2++];
8696 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
8697 (sort2[i2].nsc == sort_new[i_new].nsc &&
8698 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
8700 sort1[i1++] = sort2[i2++];
8704 sort1[i1++] = sort_new[i_new++];
8709 static int dd_sort_order(gmx_domdec_t *dd,t_forcerec *fr,int ncg_home_old)
8711 gmx_domdec_sort_t *sort;
8712 gmx_cgsort_t *cgsort,*sort_i;
8713 int ncg_new,nsort2,nsort_new,i,*a,moved,*ibuf;
8714 int sort_last,sort_skip;
8716 sort = dd->comm->sort;
8718 a = fr->ns.grid->cell_index;
8720 moved = NSGRID_SIGNAL_MOVED_FAC*fr->ns.grid->ncells;
8722 if (ncg_home_old >= 0)
8724 /* The charge groups that remained in the same ns grid cell
8725 * are completely ordered. So we can sort efficiently by sorting
8726 * the charge groups that did move into the stationary list.
8731 for(i=0; i<dd->ncg_home; i++)
8733 /* Check if this cg did not move to another node */
8736 if (i >= ncg_home_old || a[i] != sort->sort[i].nsc)
8738 /* This cg is new on this node or moved ns grid cell */
8739 if (nsort_new >= sort->sort_new_nalloc)
8741 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
8742 srenew(sort->sort_new,sort->sort_new_nalloc);
8744 sort_i = &(sort->sort_new[nsort_new++]);
8748 /* This cg did not move */
8749 sort_i = &(sort->sort2[nsort2++]);
8751 /* Sort on the ns grid cell indices
8752 * and the global topology index.
8753 * index_gl is irrelevant with cell ns,
8754 * but we set it here anyhow to avoid a conditional.
8757 sort_i->ind_gl = dd->index_gl[i];
8764 fprintf(debug,"ordered sort cgs: stationary %d moved %d\n",
8767 /* Sort efficiently */
8768 ordered_sort(nsort2,sort->sort2,nsort_new,sort->sort_new,
8773 cgsort = sort->sort;
8775 for(i=0; i<dd->ncg_home; i++)
8777 /* Sort on the ns grid cell indices
8778 * and the global topology index
8780 cgsort[i].nsc = a[i];
8781 cgsort[i].ind_gl = dd->index_gl[i];
8783 if (cgsort[i].nsc < moved)
8790 fprintf(debug,"qsort cgs: %d new home %d\n",dd->ncg_home,ncg_new);
8792 /* Determine the order of the charge groups using qsort */
8793 qsort_threadsafe(cgsort,dd->ncg_home,sizeof(cgsort[0]),comp_cgsort);
8799 static int dd_sort_order_nbnxn(gmx_domdec_t *dd,t_forcerec *fr)
8802 int ncg_new,i,*a,na;
8804 sort = dd->comm->sort->sort;
8806 nbnxn_get_atomorder(fr->nbv->nbs,&a,&na);
8813 sort[ncg_new].ind = a[i];
8821 static void dd_sort_state(gmx_domdec_t *dd,int ePBC,
8822 rvec *cgcm,t_forcerec *fr,t_state *state,
8825 gmx_domdec_sort_t *sort;
8826 gmx_cgsort_t *cgsort,*sort_i;
8828 int ncg_new,i,*ibuf,cgsize;
8831 sort = dd->comm->sort;
8833 if (dd->ncg_home > sort->sort_nalloc)
8835 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
8836 srenew(sort->sort,sort->sort_nalloc);
8837 srenew(sort->sort2,sort->sort_nalloc);
8839 cgsort = sort->sort;
8841 switch (fr->cutoff_scheme)
8844 ncg_new = dd_sort_order(dd,fr,ncg_home_old);
8847 ncg_new = dd_sort_order_nbnxn(dd,fr);
8850 gmx_incons("unimplemented");
8854 /* We alloc with the old size, since cgindex is still old */
8855 vec_rvec_check_alloc(&dd->comm->vbuf,dd->cgindex[dd->ncg_home]);
8856 vbuf = dd->comm->vbuf.v;
8860 cgindex = dd->cgindex;
8867 /* Remove the charge groups which are no longer at home here */
8868 dd->ncg_home = ncg_new;
8871 fprintf(debug,"Set the new home charge group count to %d\n",
8875 /* Reorder the state */
8876 for(i=0; i<estNR; i++)
8878 if (EST_DISTR(i) && (state->flags & (1<<i)))
8883 order_vec_atom(dd->ncg_home,cgindex,cgsort,state->x,vbuf);
8886 order_vec_atom(dd->ncg_home,cgindex,cgsort,state->v,vbuf);
8889 order_vec_atom(dd->ncg_home,cgindex,cgsort,state->sd_X,vbuf);
8892 order_vec_atom(dd->ncg_home,cgindex,cgsort,state->cg_p,vbuf);
8896 case estDISRE_INITF:
8897 case estDISRE_RM3TAV:
8898 case estORIRE_INITF:
8900 /* No ordering required */
8903 gmx_incons("Unknown state entry encountered in dd_sort_state");
8908 if (fr->cutoff_scheme == ecutsGROUP)
8911 order_vec_cg(dd->ncg_home,cgsort,cgcm,vbuf);
8914 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8916 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8917 srenew(sort->ibuf,sort->ibuf_nalloc);
8920 /* Reorder the global cg index */
8921 order_int_cg(dd->ncg_home,cgsort,dd->index_gl,ibuf);
8922 /* Reorder the cginfo */
8923 order_int_cg(dd->ncg_home,cgsort,fr->cginfo,ibuf);
8924 /* Rebuild the local cg index */
8928 for(i=0; i<dd->ncg_home; i++)
8930 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8931 ibuf[i+1] = ibuf[i] + cgsize;
8933 for(i=0; i<dd->ncg_home+1; i++)
8935 dd->cgindex[i] = ibuf[i];
8940 for(i=0; i<dd->ncg_home+1; i++)
8945 /* Set the home atom number */
8946 dd->nat_home = dd->cgindex[dd->ncg_home];
8948 if (fr->cutoff_scheme == ecutsVERLET)
8950 /* The atoms are now exactly in grid order, update the grid order */
8951 nbnxn_set_atomorder(fr->nbv->nbs);
8955 /* Copy the sorted ns cell indices back to the ns grid struct */
8956 for(i=0; i<dd->ncg_home; i++)
8958 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
8960 fr->ns.grid->nr = dd->ncg_home;
8964 static void add_dd_statistics(gmx_domdec_t *dd)
8966 gmx_domdec_comm_t *comm;
8971 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8973 comm->sum_nat[ddnat-ddnatZONE] +=
8974 comm->nat[ddnat] - comm->nat[ddnat-1];
8979 void reset_dd_statistics_counters(gmx_domdec_t *dd)
8981 gmx_domdec_comm_t *comm;
8986 /* Reset all the statistics and counters for total run counting */
8987 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8989 comm->sum_nat[ddnat-ddnatZONE] = 0;
8993 comm->load_step = 0;
8996 clear_ivec(comm->load_lim);
9001 void print_dd_statistics(t_commrec *cr,t_inputrec *ir,FILE *fplog)
9003 gmx_domdec_comm_t *comm;
9007 comm = cr->dd->comm;
9009 gmx_sumd(ddnatNR-ddnatZONE,comm->sum_nat,cr);
9016 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");
9018 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
9020 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
9025 " av. #atoms communicated per step for force: %d x %.1f\n",
9029 if (cr->dd->vsite_comm)
9032 " av. #atoms communicated per step for vsites: %d x %.1f\n",
9033 (EEL_PME(ir->coulombtype) || ir->coulombtype==eelEWALD) ? 3 : 2,
9038 if (cr->dd->constraint_comm)
9041 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
9042 1 + ir->nLincsIter,av);
9046 gmx_incons(" Unknown type for DD statistics");
9049 fprintf(fplog,"\n");
9051 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
9053 print_dd_load_av(fplog,cr->dd);
9057 void dd_partition_system(FILE *fplog,
9058 gmx_large_int_t step,
9060 gmx_bool bMasterState,
9062 t_state *state_global,
9063 gmx_mtop_t *top_global,
9065 t_state *state_local,
9068 gmx_localtop_t *top_local,
9071 gmx_shellfc_t shellfc,
9072 gmx_constr_t constr,
9074 gmx_wallcycle_t wcycle,
9078 gmx_domdec_comm_t *comm;
9079 gmx_ddbox_t ddbox={0};
9081 gmx_large_int_t step_pcoupl;
9082 rvec cell_ns_x0,cell_ns_x1;
9083 int i,j,n,cg0=0,ncg_home_old=-1,ncg_moved,nat_f_novirsum;
9084 gmx_bool bBoxChanged,bNStGlobalComm,bDoDLB,bCheckDLB,bTurnOnDLB,bLogLoad;
9085 gmx_bool bRedist,bSortCG,bResortAll;
9086 ivec ncells_old={0,0,0},ncells_new={0,0,0},np;
9093 bBoxChanged = (bMasterState || DEFORM(*ir));
9094 if (ir->epc != epcNO)
9096 /* With nstpcouple > 1 pressure coupling happens.
9097 * one step after calculating the pressure.
9098 * Box scaling happens at the end of the MD step,
9099 * after the DD partitioning.
9100 * We therefore have to do DLB in the first partitioning
9101 * after an MD step where P-coupling occured.
9102 * We need to determine the last step in which p-coupling occurred.
9103 * MRS -- need to validate this for vv?
9108 step_pcoupl = step - 1;
9112 step_pcoupl = ((step - 1)/n)*n + 1;
9114 if (step_pcoupl >= comm->partition_step)
9120 bNStGlobalComm = (step % nstglobalcomm == 0);
9122 if (!comm->bDynLoadBal)
9128 /* Should we do dynamic load balacing this step?
9129 * Since it requires (possibly expensive) global communication,
9130 * we might want to do DLB less frequently.
9132 if (bBoxChanged || ir->epc != epcNO)
9134 bDoDLB = bBoxChanged;
9138 bDoDLB = bNStGlobalComm;
9142 /* Check if we have recorded loads on the nodes */
9143 if (comm->bRecordLoad && dd_load_count(comm))
9145 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
9147 /* Check if we should use DLB at the second partitioning
9148 * and every 100 partitionings,
9149 * so the extra communication cost is negligible.
9151 n = max(100,nstglobalcomm);
9152 bCheckDLB = (comm->n_load_collect == 0 ||
9153 comm->n_load_have % n == n-1);
9160 /* Print load every nstlog, first and last step to the log file */
9161 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
9162 comm->n_load_collect == 0 ||
9164 (step + ir->nstlist > ir->init_step + ir->nsteps)));
9166 /* Avoid extra communication due to verbose screen output
9167 * when nstglobalcomm is set.
9169 if (bDoDLB || bLogLoad || bCheckDLB ||
9170 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
9172 get_load_distribution(dd,wcycle);
9177 dd_print_load(fplog,dd,step-1);
9181 dd_print_load_verbose(dd);
9184 comm->n_load_collect++;
9187 /* Since the timings are node dependent, the master decides */
9191 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
9194 fprintf(debug,"step %s, imb loss %f\n",
9195 gmx_step_str(step,sbuf),
9196 dd_force_imb_perf_loss(dd));
9199 dd_bcast(dd,sizeof(bTurnOnDLB),&bTurnOnDLB);
9202 turn_on_dlb(fplog,cr,step);
9207 comm->n_load_have++;
9210 cgs_gl = &comm->cgs_gl;
9215 /* Clear the old state */
9216 clear_dd_indices(dd,0,0);
9218 set_ddbox(dd,bMasterState,cr,ir,state_global->box,
9219 TRUE,cgs_gl,state_global->x,&ddbox);
9221 get_cg_distribution(fplog,step,dd,cgs_gl,
9222 state_global->box,&ddbox,state_global->x);
9224 dd_distribute_state(dd,cgs_gl,
9225 state_global,state_local,f);
9227 dd_make_local_cgs(dd,&top_local->cgs);
9229 /* Ensure that we have space for the new distribution */
9230 dd_check_alloc_ncg(fr,state_local,f,dd->ncg_home);
9232 if (fr->cutoff_scheme == ecutsGROUP)
9234 calc_cgcm(fplog,0,dd->ncg_home,
9235 &top_local->cgs,state_local->x,fr->cg_cm);
9238 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
9240 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
9244 else if (state_local->ddp_count != dd->ddp_count)
9246 if (state_local->ddp_count > dd->ddp_count)
9248 gmx_fatal(FARGS,"Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)",state_local->ddp_count,dd->ddp_count);
9251 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
9253 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);
9256 /* Clear the old state */
9257 clear_dd_indices(dd,0,0);
9259 /* Build the new indices */
9260 rebuild_cgindex(dd,cgs_gl->index,state_local);
9261 make_dd_indices(dd,cgs_gl->index,0);
9263 if (fr->cutoff_scheme == ecutsGROUP)
9265 /* Redetermine the cg COMs */
9266 calc_cgcm(fplog,0,dd->ncg_home,
9267 &top_local->cgs,state_local->x,fr->cg_cm);
9270 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
9272 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
9274 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
9275 TRUE,&top_local->cgs,state_local->x,&ddbox);
9277 bRedist = comm->bDynLoadBal;
9281 /* We have the full state, only redistribute the cgs */
9283 /* Clear the non-home indices */
9284 clear_dd_indices(dd,dd->ncg_home,dd->nat_home);
9286 /* Avoid global communication for dim's without pbc and -gcom */
9287 if (!bNStGlobalComm)
9289 copy_rvec(comm->box0 ,ddbox.box0 );
9290 copy_rvec(comm->box_size,ddbox.box_size);
9292 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
9293 bNStGlobalComm,&top_local->cgs,state_local->x,&ddbox);
9298 /* For dim's without pbc and -gcom */
9299 copy_rvec(ddbox.box0 ,comm->box0 );
9300 copy_rvec(ddbox.box_size,comm->box_size);
9302 set_dd_cell_sizes(dd,&ddbox,dynamic_dd_box(&ddbox,ir),bMasterState,bDoDLB,
9305 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
9307 write_dd_grid_pdb("dd_grid",step,dd,state_local->box,&ddbox);
9310 /* Check if we should sort the charge groups */
9311 if (comm->nstSortCG > 0)
9313 bSortCG = (bMasterState ||
9314 (bRedist && (step % comm->nstSortCG == 0)));
9321 ncg_home_old = dd->ncg_home;
9326 wallcycle_sub_start(wcycle,ewcsDD_REDIST);
9328 dd_redistribute_cg(fplog,step,dd,ddbox.tric_dir,
9329 state_local,f,fr,mdatoms,
9330 !bSortCG,nrnb,&cg0,&ncg_moved);
9332 wallcycle_sub_stop(wcycle,ewcsDD_REDIST);
9335 get_nsgrid_boundaries(ddbox.nboundeddim,state_local->box,
9337 &comm->cell_x0,&comm->cell_x1,
9338 dd->ncg_home,fr->cg_cm,
9339 cell_ns_x0,cell_ns_x1,&grid_density);
9343 comm_dd_ns_cell_sizes(dd,&ddbox,cell_ns_x0,cell_ns_x1,step);
9346 switch (fr->cutoff_scheme)
9349 copy_ivec(fr->ns.grid->n,ncells_old);
9350 grid_first(fplog,fr->ns.grid,dd,&ddbox,fr->ePBC,
9351 state_local->box,cell_ns_x0,cell_ns_x1,
9352 fr->rlistlong,grid_density);
9355 nbnxn_get_ncells(fr->nbv->nbs,&ncells_old[XX],&ncells_old[YY]);
9358 gmx_incons("unimplemented");
9360 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
9361 copy_ivec(ddbox.tric_dir,comm->tric_dir);
9365 wallcycle_sub_start(wcycle,ewcsDD_GRID);
9367 /* Sort the state on charge group position.
9368 * This enables exact restarts from this step.
9369 * It also improves performance by about 15% with larger numbers
9370 * of atoms per node.
9373 /* Fill the ns grid with the home cell,
9374 * so we can sort with the indices.
9376 set_zones_ncg_home(dd);
9378 switch (fr->cutoff_scheme)
9381 set_zones_size(dd,state_local->box,&ddbox,0,1);
9383 nbnxn_put_on_grid(fr->nbv->nbs,fr->ePBC,state_local->box,
9385 comm->zones.size[0].bb_x0,
9386 comm->zones.size[0].bb_x1,
9388 comm->zones.dens_zone0,
9391 ncg_moved,bRedist ? comm->moved : NULL,
9392 fr->nbv->grp[eintLocal].kernel_type,
9393 fr->nbv->grp[eintLocal].nbat);
9395 nbnxn_get_ncells(fr->nbv->nbs,&ncells_new[XX],&ncells_new[YY]);
9398 fill_grid(fplog,&comm->zones,fr->ns.grid,dd->ncg_home,
9399 0,dd->ncg_home,fr->cg_cm);
9401 copy_ivec(fr->ns.grid->n,ncells_new);
9404 gmx_incons("unimplemented");
9407 bResortAll = bMasterState;
9409 /* Check if we can user the old order and ns grid cell indices
9410 * of the charge groups to sort the charge groups efficiently.
9412 if (ncells_new[XX] != ncells_old[XX] ||
9413 ncells_new[YY] != ncells_old[YY] ||
9414 ncells_new[ZZ] != ncells_old[ZZ])
9421 fprintf(debug,"Step %s, sorting the %d home charge groups\n",
9422 gmx_step_str(step,sbuf),dd->ncg_home);
9424 dd_sort_state(dd,ir->ePBC,fr->cg_cm,fr,state_local,
9425 bResortAll ? -1 : ncg_home_old);
9426 /* Rebuild all the indices */
9428 ga2la_clear(dd->ga2la);
9430 wallcycle_sub_stop(wcycle,ewcsDD_GRID);
9433 wallcycle_sub_start(wcycle,ewcsDD_SETUPCOMM);
9435 /* Setup up the communication and communicate the coordinates */
9436 setup_dd_communication(dd,state_local->box,&ddbox,fr,state_local,f);
9438 /* Set the indices */
9439 make_dd_indices(dd,cgs_gl->index,cg0);
9441 /* Set the charge group boundaries for neighbor searching */
9442 set_cg_boundaries(&comm->zones);
9444 if (fr->cutoff_scheme == ecutsVERLET)
9446 set_zones_size(dd,state_local->box,&ddbox,
9447 bSortCG ? 1 : 0,comm->zones.n);
9450 wallcycle_sub_stop(wcycle,ewcsDD_SETUPCOMM);
9453 write_dd_pdb("dd_home",step,"dump",top_global,cr,
9454 -1,state_local->x,state_local->box);
9457 wallcycle_sub_start(wcycle,ewcsDD_MAKETOP);
9459 /* Extract a local topology from the global topology */
9460 for(i=0; i<dd->ndim; i++)
9462 np[dd->dim[i]] = comm->cd[i].np;
9464 dd_make_local_top(fplog,dd,&comm->zones,dd->npbcdim,state_local->box,
9465 comm->cellsize_min,np,
9467 fr->cutoff_scheme==ecutsGROUP ? fr->cg_cm : state_local->x,
9468 vsite,top_global,top_local);
9470 wallcycle_sub_stop(wcycle,ewcsDD_MAKETOP);
9472 wallcycle_sub_start(wcycle,ewcsDD_MAKECONSTR);
9474 /* Set up the special atom communication */
9475 n = comm->nat[ddnatZONE];
9476 for(i=ddnatZONE+1; i<ddnatNR; i++)
9481 if (vsite && vsite->n_intercg_vsite)
9483 n = dd_make_local_vsites(dd,n,top_local->idef.il);
9487 if (dd->bInterCGcons || dd->bInterCGsettles)
9489 /* Only for inter-cg constraints we need special code */
9490 n = dd_make_local_constraints(dd,n,top_global,fr->cginfo,
9491 constr,ir->nProjOrder,
9492 top_local->idef.il);
9496 gmx_incons("Unknown special atom type setup");
9501 wallcycle_sub_stop(wcycle,ewcsDD_MAKECONSTR);
9503 wallcycle_sub_start(wcycle,ewcsDD_TOPOTHER);
9505 /* Make space for the extra coordinates for virtual site
9506 * or constraint communication.
9508 state_local->natoms = comm->nat[ddnatNR-1];
9509 if (state_local->natoms > state_local->nalloc)
9511 dd_realloc_state(state_local,f,state_local->natoms);
9514 if (fr->bF_NoVirSum)
9516 if (vsite && vsite->n_intercg_vsite)
9518 nat_f_novirsum = comm->nat[ddnatVSITE];
9522 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
9524 nat_f_novirsum = dd->nat_tot;
9528 nat_f_novirsum = dd->nat_home;
9537 /* Set the number of atoms required for the force calculation.
9538 * Forces need to be constrained when using a twin-range setup
9539 * or with energy minimization. For simple simulations we could
9540 * avoid some allocation, zeroing and copying, but this is
9541 * probably not worth the complications ande checking.
9543 forcerec_set_ranges(fr,dd->ncg_home,dd->ncg_tot,
9544 dd->nat_tot,comm->nat[ddnatCON],nat_f_novirsum);
9546 /* We make the all mdatoms up to nat_tot_con.
9547 * We could save some work by only setting invmass
9548 * between nat_tot and nat_tot_con.
9550 /* This call also sets the new number of home particles to dd->nat_home */
9551 atoms2md(top_global,ir,
9552 comm->nat[ddnatCON],dd->gatindex,0,dd->nat_home,mdatoms);
9554 /* Now we have the charges we can sort the FE interactions */
9555 dd_sort_local_top(dd,mdatoms,top_local);
9559 /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
9560 split_vsites_over_threads(top_local->idef.il,mdatoms,FALSE,vsite);
9565 /* Make the local shell stuff, currently no communication is done */
9566 make_local_shells(cr,mdatoms,shellfc);
9569 if (ir->implicit_solvent)
9571 make_local_gb(cr,fr->born,ir->gb_algorithm);
9574 init_bonded_thread_force_reduction(fr,&top_local->idef);
9576 if (!(cr->duty & DUTY_PME))
9578 /* Send the charges to our PME only node */
9579 gmx_pme_send_q(cr,mdatoms->nChargePerturbed,
9580 mdatoms->chargeA,mdatoms->chargeB,
9581 dd_pme_maxshift_x(dd),dd_pme_maxshift_y(dd));
9586 set_constraints(constr,top_local,ir,mdatoms,cr);
9589 if (ir->ePull != epullNO)
9591 /* Update the local pull groups */
9592 dd_make_local_pull_groups(dd,ir->pull,mdatoms);
9597 /* Update the local rotation groups */
9598 dd_make_local_rotation_groups(dd,ir->rot);
9602 add_dd_statistics(dd);
9604 /* Make sure we only count the cycles for this DD partitioning */
9605 clear_dd_cycle_counts(dd);
9607 /* Because the order of the atoms might have changed since
9608 * the last vsite construction, we need to communicate the constructing
9609 * atom coordinates again (for spreading the forces this MD step).
9611 dd_move_x_vsites(dd,state_local->box,state_local->x);
9613 wallcycle_sub_stop(wcycle,ewcsDD_TOPOTHER);
9615 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
9617 dd_move_x(dd,state_local->box,state_local->x);
9618 write_dd_pdb("dd_dump",step,"dump",top_global,cr,
9619 -1,state_local->x,state_local->box);
9622 /* Store the partitioning step */
9623 comm->partition_step = step;
9625 /* Increase the DD partitioning counter */
9627 /* The state currently matches this DD partitioning count, store it */
9628 state_local->ddp_count = dd->ddp_count;
9631 /* The DD master node knows the complete cg distribution,
9632 * store the count so we can possibly skip the cg info communication.
9634 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
9637 if (comm->DD_debug > 0)
9639 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
9640 check_index_consistency(dd,top_global->natoms,ncg_mtop(top_global),
9641 "after partitioning");