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
32 #include "domdec_network.h"
35 #include "chargegroup.h"
44 #include "pull_rotation.h"
45 #include "gmx_wallcycle.h"
49 #include "mtop_util.h"
51 #include "gmx_ga2la.h"
62 #define DDRANK(dd,rank) (rank)
63 #define DDMASTERRANK(dd) (dd->masterrank)
65 typedef struct gmx_domdec_master
67 /* The cell boundaries */
69 /* The global charge group division */
70 int *ncg; /* Number of home charge groups for each node */
71 int *index; /* Index of nnodes+1 into cg */
72 int *cg; /* Global charge group index */
73 int *nat; /* Number of home atoms for each node. */
74 int *ibuf; /* Buffer for communication */
75 rvec *vbuf; /* Buffer for state scattering and gathering */
76 } gmx_domdec_master_t;
80 /* The numbers of charge groups to send and receive for each cell
81 * that requires communication, the last entry contains the total
82 * number of atoms that needs to be communicated.
84 int nsend[DD_MAXIZONE+2];
85 int nrecv[DD_MAXIZONE+2];
86 /* The charge groups to send */
89 /* The atom range for non-in-place communication */
90 int cell2at0[DD_MAXIZONE];
91 int cell2at1[DD_MAXIZONE];
96 int np; /* Number of grid pulses in this dimension */
97 int np_dlb; /* For dlb, for use with edlbAUTO */
98 gmx_domdec_ind_t *ind; /* The indices to communicate, size np */
100 gmx_bool bInPlace; /* Can we communicate in place? */
101 } gmx_domdec_comm_dim_t;
105 gmx_bool *bCellMin; /* Temp. var.: is this cell size at the limit */
106 real *cell_f; /* State var.: cell boundaries, box relative */
107 real *old_cell_f; /* Temp. var.: old cell size */
108 real *cell_f_max0; /* State var.: max lower boundary, incl neighbors */
109 real *cell_f_min1; /* State var.: min upper boundary, incl neighbors */
110 real *bound_min; /* Temp. var.: lower limit for cell boundary */
111 real *bound_max; /* Temp. var.: upper limit for cell boundary */
112 gmx_bool bLimited; /* State var.: is DLB limited in this dim and row */
113 real *buf_ncd; /* Temp. var. */
116 #define DD_NLOAD_MAX 9
118 /* Here floats are accurate enough, since these variables
119 * only influence the load balancing, not the actual MD results.
143 gmx_cgsort_t *sort1,*sort2;
145 gmx_cgsort_t *sort_new;
157 /* This enum determines the order of the coordinates.
158 * ddnatHOME and ddnatZONE should be first and second,
159 * the others can be ordered as wanted.
161 enum { ddnatHOME, ddnatZONE, ddnatVSITE, ddnatCON, ddnatNR };
163 enum { edlbAUTO, edlbNO, edlbYES, edlbNR };
164 const char *edlb_names[edlbNR] = { "auto", "no", "yes" };
168 int dim; /* The dimension */
169 gmx_bool dim_match;/* Tells if DD and PME dims match */
170 int nslab; /* The number of PME slabs in this dimension */
171 real *slb_dim_f; /* Cell sizes for determining the PME comm. with SLB */
172 int *pp_min; /* The minimum pp node location, size nslab */
173 int *pp_max; /* The maximum pp node location,size nslab */
174 int maxshift; /* The maximum shift for coordinate redistribution in PME */
179 real min0; /* The minimum bottom of this zone */
180 real max1; /* The maximum top of this zone */
181 real mch0; /* The maximum bottom communicaton height for this zone */
182 real mch1; /* The maximum top communicaton height for this zone */
183 real p1_0; /* The bottom value of the first cell in this zone */
184 real p1_1; /* The top value of the first cell in this zone */
187 typedef struct gmx_domdec_comm
189 /* All arrays are indexed with 0 to dd->ndim (not Cartesian indexing),
190 * unless stated otherwise.
193 /* The number of decomposition dimensions for PME, 0: no PME */
195 /* The number of nodes doing PME (PP/PME or only PME) */
199 /* The communication setup including the PME only nodes */
200 gmx_bool bCartesianPP_PME;
203 int *pmenodes; /* size npmenodes */
204 int *ddindex2simnodeid; /* size npmenodes, only with bCartesianPP
205 * but with bCartesianPP_PME */
206 gmx_ddpme_t ddpme[2];
208 /* The DD particle-particle nodes only */
209 gmx_bool bCartesianPP;
210 int *ddindex2ddnodeid; /* size npmenode, only with bCartesianPP_PME */
212 /* The global charge groups */
215 /* Should we sort the cgs */
217 gmx_domdec_sort_t *sort;
219 /* Are there bonded and multi-body interactions between charge groups? */
220 gmx_bool bInterCGBondeds;
221 gmx_bool bInterCGMultiBody;
223 /* Data for the optional bonded interaction atom communication range */
230 /* Are we actually using DLB? */
231 gmx_bool bDynLoadBal;
233 /* Cell sizes for static load balancing, first index cartesian */
236 /* The width of the communicated boundaries */
239 /* The minimum cell size (including triclinic correction) */
241 /* For dlb, for use with edlbAUTO */
242 rvec cellsize_min_dlb;
243 /* The lower limit for the DD cell size with DLB */
245 /* Effectively no NB cut-off limit with DLB for systems without PBC? */
246 gmx_bool bVacDLBNoLimit;
248 /* tric_dir is only stored here because dd_get_ns_ranges needs it */
250 /* box0 and box_size are required with dim's without pbc and -gcom */
254 /* The cell boundaries */
258 /* The old location of the cell boundaries, to check cg displacements */
262 /* The communication setup and charge group boundaries for the zones */
263 gmx_domdec_zones_t zones;
265 /* The zone limits for DD dimensions 1 and 2 (not 0), determined from
266 * cell boundaries of neighboring cells for dynamic load balancing.
268 gmx_ddzone_t zone_d1[2];
269 gmx_ddzone_t zone_d2[2][2];
271 /* The coordinate/force communication setup and indices */
272 gmx_domdec_comm_dim_t cd[DIM];
273 /* The maximum number of cells to communicate with in one dimension */
276 /* Which cg distribution is stored on the master node */
277 int master_cg_ddp_count;
279 /* The number of cg's received from the direct neighbors */
280 int zone_ncg1[DD_MAXZONE];
282 /* The atom counts, the range for each type t is nat[t-1] <= at < nat[t] */
285 /* Communication buffer for general use */
289 /* Communication buffer for general use */
292 /* Communication buffers only used with multiple grid pulses */
297 /* Communication buffers for local redistribution */
299 int cggl_flag_nalloc[DIM*2];
301 int cgcm_state_nalloc[DIM*2];
303 /* Cell sizes for dynamic load balancing */
304 gmx_domdec_root_t **root;
308 real cell_f_max0[DIM];
309 real cell_f_min1[DIM];
311 /* Stuff for load communication */
312 gmx_bool bRecordLoad;
313 gmx_domdec_load_t *load;
315 MPI_Comm *mpi_comm_load;
318 /* Maximum DLB scaling per load balancing step in percent */
322 float cycl[ddCyclNr];
323 int cycl_n[ddCyclNr];
324 float cycl_max[ddCyclNr];
325 /* Flop counter (0=no,1=yes,2=with (eFlop-1)*5% noise */
329 /* Have often have did we have load measurements */
331 /* Have often have we collected the load measurements */
335 double sum_nat[ddnatNR-ddnatZONE];
345 /* The last partition step */
346 gmx_large_int_t globalcomm_step;
354 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
357 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
358 #define DD_FLAG_NRCG 65535
359 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
360 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
362 /* Zone permutation required to obtain consecutive charge groups
363 * for neighbor searching.
365 static const int zone_perm[3][4] = { {0,0,0,0},{1,0,0,0},{3,0,1,2} };
367 /* dd_zo and dd_zp3/dd_zp2 are set up such that i zones with non-zero
368 * components see only j zones with that component 0.
371 /* The DD zone order */
372 static const ivec dd_zo[DD_MAXZONE] =
373 {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,1,1},{0,0,1},{1,0,1},{1,1,1}};
378 static const ivec dd_zp3[dd_zp3n] = {{0,0,8},{1,3,6},{2,5,6},{3,5,7}};
383 static const ivec dd_zp2[dd_zp2n] = {{0,0,4},{1,3,4}};
388 static const ivec dd_zp1[dd_zp1n] = {{0,0,2}};
390 /* Factors used to avoid problems due to rounding issues */
391 #define DD_CELL_MARGIN 1.0001
392 #define DD_CELL_MARGIN2 1.00005
393 /* Factor to account for pressure scaling during nstlist steps */
394 #define DD_PRES_SCALE_MARGIN 1.02
396 /* Allowed performance loss before we DLB or warn */
397 #define DD_PERF_LOSS 0.05
399 #define DD_CELL_F_SIZE(dd,di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
401 /* Use separate MPI send and receive commands
402 * when nnodes <= GMX_DD_NNODES_SENDRECV.
403 * This saves memory (and some copying for small nnodes).
404 * For high parallelization scatter and gather calls are used.
406 #define GMX_DD_NNODES_SENDRECV 4
410 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
412 static void index2xyz(ivec nc,int ind,ivec xyz)
414 xyz[XX] = ind % nc[XX];
415 xyz[YY] = (ind / nc[XX]) % nc[YY];
416 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
420 /* This order is required to minimize the coordinate communication in PME
421 * which uses decomposition in the x direction.
423 #define dd_index(n,i) ((((i)[XX]*(n)[YY] + (i)[YY])*(n)[ZZ]) + (i)[ZZ])
425 static void ddindex2xyz(ivec nc,int ind,ivec xyz)
427 xyz[XX] = ind / (nc[YY]*nc[ZZ]);
428 xyz[YY] = (ind / nc[ZZ]) % nc[YY];
429 xyz[ZZ] = ind % nc[ZZ];
432 static int ddcoord2ddnodeid(gmx_domdec_t *dd,ivec c)
437 ddindex = dd_index(dd->nc,c);
438 if (dd->comm->bCartesianPP_PME)
440 ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
442 else if (dd->comm->bCartesianPP)
445 MPI_Cart_rank(dd->mpi_comm_all,c,&ddnodeid);
456 static gmx_bool dynamic_dd_box(gmx_ddbox_t *ddbox,t_inputrec *ir)
458 return (ddbox->nboundeddim < DIM || DYNAMIC_BOX(*ir));
461 int ddglatnr(gmx_domdec_t *dd,int i)
471 if (i >= dd->comm->nat[ddnatNR-1])
473 gmx_fatal(FARGS,"glatnr called with %d, which is larger than the local number of atoms (%d)",i,dd->comm->nat[ddnatNR-1]);
475 atnr = dd->gatindex[i] + 1;
481 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
483 return &dd->comm->cgs_gl;
486 static void vec_rvec_init(vec_rvec_t *v)
492 static void vec_rvec_check_alloc(vec_rvec_t *v,int n)
496 v->nalloc = over_alloc_dd(n);
497 srenew(v->v,v->nalloc);
501 void dd_store_state(gmx_domdec_t *dd,t_state *state)
505 if (state->ddp_count != dd->ddp_count)
507 gmx_incons("The state does not the domain decomposition state");
510 state->ncg_gl = dd->ncg_home;
511 if (state->ncg_gl > state->cg_gl_nalloc)
513 state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
514 srenew(state->cg_gl,state->cg_gl_nalloc);
516 for(i=0; i<state->ncg_gl; i++)
518 state->cg_gl[i] = dd->index_gl[i];
521 state->ddp_count_cg_gl = dd->ddp_count;
524 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
526 return &dd->comm->zones;
529 void dd_get_ns_ranges(gmx_domdec_t *dd,int icg,
530 int *jcg0,int *jcg1,ivec shift0,ivec shift1)
532 gmx_domdec_zones_t *zones;
535 zones = &dd->comm->zones;
538 while (icg >= zones->izone[izone].cg1)
547 else if (izone < zones->nizone)
549 *jcg0 = zones->izone[izone].jcg0;
553 gmx_fatal(FARGS,"DD icg %d out of range: izone (%d) >= nizone (%d)",
554 icg,izone,zones->nizone);
557 *jcg1 = zones->izone[izone].jcg1;
559 for(d=0; d<dd->ndim; d++)
562 shift0[dim] = zones->izone[izone].shift0[dim];
563 shift1[dim] = zones->izone[izone].shift1[dim];
564 if (dd->comm->tric_dir[dim] || (dd->bGridJump && d > 0))
566 /* A conservative approach, this can be optimized */
573 int dd_natoms_vsite(gmx_domdec_t *dd)
575 return dd->comm->nat[ddnatVSITE];
578 void dd_get_constraint_range(gmx_domdec_t *dd,int *at_start,int *at_end)
580 *at_start = dd->comm->nat[ddnatCON-1];
581 *at_end = dd->comm->nat[ddnatCON];
584 void dd_move_x(gmx_domdec_t *dd,matrix box,rvec x[])
586 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
588 gmx_domdec_comm_t *comm;
589 gmx_domdec_comm_dim_t *cd;
590 gmx_domdec_ind_t *ind;
591 rvec shift={0,0,0},*buf,*rbuf;
592 gmx_bool bPBC,bScrew;
596 cgindex = dd->cgindex;
601 nat_tot = dd->nat_home;
602 for(d=0; d<dd->ndim; d++)
604 bPBC = (dd->ci[dd->dim[d]] == 0);
605 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
608 copy_rvec(box[dd->dim[d]],shift);
611 for(p=0; p<cd->np; p++)
618 for(i=0; i<ind->nsend[nzone]; i++)
620 at0 = cgindex[index[i]];
621 at1 = cgindex[index[i]+1];
622 for(j=at0; j<at1; j++)
624 copy_rvec(x[j],buf[n]);
631 for(i=0; i<ind->nsend[nzone]; i++)
633 at0 = cgindex[index[i]];
634 at1 = cgindex[index[i]+1];
635 for(j=at0; j<at1; j++)
637 /* We need to shift the coordinates */
638 rvec_add(x[j],shift,buf[n]);
645 for(i=0; i<ind->nsend[nzone]; i++)
647 at0 = cgindex[index[i]];
648 at1 = cgindex[index[i]+1];
649 for(j=at0; j<at1; j++)
652 buf[n][XX] = x[j][XX] + shift[XX];
654 * This operation requires a special shift force
655 * treatment, which is performed in calc_vir.
657 buf[n][YY] = box[YY][YY] - x[j][YY];
658 buf[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
670 rbuf = comm->vbuf2.v;
672 /* Send and receive the coordinates */
673 dd_sendrecv_rvec(dd, d, dddirBackward,
674 buf, ind->nsend[nzone+1],
675 rbuf, ind->nrecv[nzone+1]);
679 for(zone=0; zone<nzone; zone++)
681 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
683 copy_rvec(rbuf[j],x[i]);
688 nat_tot += ind->nrecv[nzone+1];
694 void dd_move_f(gmx_domdec_t *dd,rvec f[],rvec *fshift)
696 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
698 gmx_domdec_comm_t *comm;
699 gmx_domdec_comm_dim_t *cd;
700 gmx_domdec_ind_t *ind;
704 gmx_bool bPBC,bScrew;
708 cgindex = dd->cgindex;
713 nzone = comm->zones.n/2;
714 nat_tot = dd->nat_tot;
715 for(d=dd->ndim-1; d>=0; d--)
717 bPBC = (dd->ci[dd->dim[d]] == 0);
718 bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
719 if (fshift == NULL && !bScrew)
723 /* Determine which shift vector we need */
729 for(p=cd->np-1; p>=0; p--) {
731 nat_tot -= ind->nrecv[nzone+1];
738 sbuf = comm->vbuf2.v;
740 for(zone=0; zone<nzone; zone++)
742 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
744 copy_rvec(f[i],sbuf[j]);
749 /* Communicate the forces */
750 dd_sendrecv_rvec(dd, d, dddirForward,
751 sbuf, ind->nrecv[nzone+1],
752 buf, ind->nsend[nzone+1]);
754 /* Add the received forces */
758 for(i=0; i<ind->nsend[nzone]; i++)
760 at0 = cgindex[index[i]];
761 at1 = cgindex[index[i]+1];
762 for(j=at0; j<at1; j++)
764 rvec_inc(f[j],buf[n]);
771 for(i=0; i<ind->nsend[nzone]; i++)
773 at0 = cgindex[index[i]];
774 at1 = cgindex[index[i]+1];
775 for(j=at0; j<at1; j++)
777 rvec_inc(f[j],buf[n]);
778 /* Add this force to the shift force */
779 rvec_inc(fshift[is],buf[n]);
786 for(i=0; i<ind->nsend[nzone]; i++)
788 at0 = cgindex[index[i]];
789 at1 = cgindex[index[i]+1];
790 for(j=at0; j<at1; j++)
792 /* Rotate the force */
793 f[j][XX] += buf[n][XX];
794 f[j][YY] -= buf[n][YY];
795 f[j][ZZ] -= buf[n][ZZ];
798 /* Add this force to the shift force */
799 rvec_inc(fshift[is],buf[n]);
810 void dd_atom_spread_real(gmx_domdec_t *dd,real v[])
812 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
814 gmx_domdec_comm_t *comm;
815 gmx_domdec_comm_dim_t *cd;
816 gmx_domdec_ind_t *ind;
821 cgindex = dd->cgindex;
823 buf = &comm->vbuf.v[0][0];
826 nat_tot = dd->nat_home;
827 for(d=0; d<dd->ndim; d++)
830 for(p=0; p<cd->np; p++)
835 for(i=0; i<ind->nsend[nzone]; i++)
837 at0 = cgindex[index[i]];
838 at1 = cgindex[index[i]+1];
839 for(j=at0; j<at1; j++)
852 rbuf = &comm->vbuf2.v[0][0];
854 /* Send and receive the coordinates */
855 dd_sendrecv_real(dd, d, dddirBackward,
856 buf, ind->nsend[nzone+1],
857 rbuf, ind->nrecv[nzone+1]);
861 for(zone=0; zone<nzone; zone++)
863 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
870 nat_tot += ind->nrecv[nzone+1];
876 void dd_atom_sum_real(gmx_domdec_t *dd,real v[])
878 int nzone,nat_tot,n,d,p,i,j,at0,at1,zone;
880 gmx_domdec_comm_t *comm;
881 gmx_domdec_comm_dim_t *cd;
882 gmx_domdec_ind_t *ind;
887 cgindex = dd->cgindex;
889 buf = &comm->vbuf.v[0][0];
892 nzone = comm->zones.n/2;
893 nat_tot = dd->nat_tot;
894 for(d=dd->ndim-1; d>=0; d--)
897 for(p=cd->np-1; p>=0; p--) {
899 nat_tot -= ind->nrecv[nzone+1];
906 sbuf = &comm->vbuf2.v[0][0];
908 for(zone=0; zone<nzone; zone++)
910 for(i=ind->cell2at0[zone]; i<ind->cell2at1[zone]; i++)
917 /* Communicate the forces */
918 dd_sendrecv_real(dd, d, dddirForward,
919 sbuf, ind->nrecv[nzone+1],
920 buf, ind->nsend[nzone+1]);
922 /* Add the received forces */
924 for(i=0; i<ind->nsend[nzone]; i++)
926 at0 = cgindex[index[i]];
927 at1 = cgindex[index[i]+1];
928 for(j=at0; j<at1; j++)
939 static void print_ddzone(FILE *fp,int d,int i,int j,gmx_ddzone_t *zone)
941 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",
943 zone->min0,zone->max1,
944 zone->mch0,zone->mch0,
945 zone->p1_0,zone->p1_1);
948 static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
949 int ddimind,int direction,
950 gmx_ddzone_t *buf_s,int n_s,
951 gmx_ddzone_t *buf_r,int n_r)
953 rvec vbuf_s[5*2],vbuf_r[5*2];
958 vbuf_s[i*2 ][0] = buf_s[i].min0;
959 vbuf_s[i*2 ][1] = buf_s[i].max1;
960 vbuf_s[i*2 ][2] = buf_s[i].mch0;
961 vbuf_s[i*2+1][0] = buf_s[i].mch1;
962 vbuf_s[i*2+1][1] = buf_s[i].p1_0;
963 vbuf_s[i*2+1][2] = buf_s[i].p1_1;
966 dd_sendrecv_rvec(dd, ddimind, direction,
972 buf_r[i].min0 = vbuf_r[i*2 ][0];
973 buf_r[i].max1 = vbuf_r[i*2 ][1];
974 buf_r[i].mch0 = vbuf_r[i*2 ][2];
975 buf_r[i].mch1 = vbuf_r[i*2+1][0];
976 buf_r[i].p1_0 = vbuf_r[i*2+1][1];
977 buf_r[i].p1_1 = vbuf_r[i*2+1][2];
981 static void dd_move_cellx(gmx_domdec_t *dd,gmx_ddbox_t *ddbox,
982 rvec cell_ns_x0,rvec cell_ns_x1)
984 int d,d1,dim,dim1,pos,buf_size,i,j,k,p,npulse,npulse_min;
985 gmx_ddzone_t *zp,buf_s[5],buf_r[5],buf_e[5];
986 rvec extr_s[2],extr_r[2];
989 gmx_domdec_comm_t *comm;
994 for(d=1; d<dd->ndim; d++)
997 zp = (d == 1) ? &comm->zone_d1[0] : &comm->zone_d2[0][0];
998 zp->min0 = cell_ns_x0[dim];
999 zp->max1 = cell_ns_x1[dim];
1000 zp->mch0 = cell_ns_x0[dim];
1001 zp->mch1 = cell_ns_x1[dim];
1002 zp->p1_0 = cell_ns_x0[dim];
1003 zp->p1_1 = cell_ns_x1[dim];
1006 for(d=dd->ndim-2; d>=0; d--)
1009 bPBC = (dim < ddbox->npbcdim);
1011 /* Use an rvec to store two reals */
1012 extr_s[d][0] = comm->cell_f0[d+1];
1013 extr_s[d][1] = comm->cell_f1[d+1];
1017 /* Store the extremes in the backward sending buffer,
1018 * so the get updated separately from the forward communication.
1020 for(d1=d; d1<dd->ndim-1; d1++)
1022 /* We invert the order to be able to use the same loop for buf_e */
1023 buf_s[pos].min0 = extr_s[d1][1];
1024 buf_s[pos].max1 = extr_s[d1][0];
1025 buf_s[pos].mch0 = 0;
1026 buf_s[pos].mch1 = 0;
1027 /* Store the cell corner of the dimension we communicate along */
1028 buf_s[pos].p1_0 = comm->cell_x0[dim];
1029 buf_s[pos].p1_1 = 0;
1033 buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
1036 if (dd->ndim == 3 && d == 0)
1038 buf_s[pos] = comm->zone_d2[0][1];
1040 buf_s[pos] = comm->zone_d1[0];
1044 /* We only need to communicate the extremes
1045 * in the forward direction
1047 npulse = comm->cd[d].np;
1050 /* Take the minimum to avoid double communication */
1051 npulse_min = min(npulse,dd->nc[dim]-1-npulse);
1055 /* Without PBC we should really not communicate over
1056 * the boundaries, but implementing that complicates
1057 * the communication setup and therefore we simply
1058 * do all communication, but ignore some data.
1060 npulse_min = npulse;
1062 for(p=0; p<npulse_min; p++)
1064 /* Communicate the extremes forward */
1065 bUse = (bPBC || dd->ci[dim] > 0);
1067 dd_sendrecv_rvec(dd, d, dddirForward,
1068 extr_s+d, dd->ndim-d-1,
1069 extr_r+d, dd->ndim-d-1);
1073 for(d1=d; d1<dd->ndim-1; d1++)
1075 extr_s[d1][0] = max(extr_s[d1][0],extr_r[d1][0]);
1076 extr_s[d1][1] = min(extr_s[d1][1],extr_r[d1][1]);
1082 for(p=0; p<npulse; p++)
1084 /* Communicate all the zone information backward */
1085 bUse = (bPBC || dd->ci[dim] < dd->nc[dim] - 1);
1087 dd_sendrecv_ddzone(dd, d, dddirBackward,
1094 for(d1=d+1; d1<dd->ndim; d1++)
1096 /* Determine the decrease of maximum required
1097 * communication height along d1 due to the distance along d,
1098 * this avoids a lot of useless atom communication.
1100 dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
1102 if (ddbox->tric_dir[dim])
1104 /* c is the off-diagonal coupling between the cell planes
1105 * along directions d and d1.
1107 c = ddbox->v[dim][dd->dim[d1]][dim];
1113 det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
1116 dh[d1] = comm->cutoff - (c*dist_d + sqrt(det))/(1 + c*c);
1120 /* A negative value signals out of range */
1126 /* Accumulate the extremes over all pulses */
1127 for(i=0; i<buf_size; i++)
1131 buf_e[i] = buf_r[i];
1137 buf_e[i].min0 = min(buf_e[i].min0,buf_r[i].min0);
1138 buf_e[i].max1 = max(buf_e[i].max1,buf_r[i].max1);
1141 if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
1149 if (bUse && dh[d1] >= 0)
1151 buf_e[i].mch0 = max(buf_e[i].mch0,buf_r[i].mch0-dh[d1]);
1152 buf_e[i].mch1 = max(buf_e[i].mch1,buf_r[i].mch1-dh[d1]);
1155 /* Copy the received buffer to the send buffer,
1156 * to pass the data through with the next pulse.
1158 buf_s[i] = buf_r[i];
1160 if (((bPBC || dd->ci[dim]+npulse < dd->nc[dim]) && p == npulse-1) ||
1161 (!bPBC && dd->ci[dim]+1+p == dd->nc[dim]-1))
1163 /* Store the extremes */
1166 for(d1=d; d1<dd->ndim-1; d1++)
1168 extr_s[d1][1] = min(extr_s[d1][1],buf_e[pos].min0);
1169 extr_s[d1][0] = max(extr_s[d1][0],buf_e[pos].max1);
1173 if (d == 1 || (d == 0 && dd->ndim == 3))
1177 comm->zone_d2[1-d][i] = buf_e[pos];
1183 comm->zone_d1[1] = buf_e[pos];
1197 print_ddzone(debug,1,i,0,&comm->zone_d1[i]);
1199 cell_ns_x0[dim] = min(cell_ns_x0[dim],comm->zone_d1[i].min0);
1200 cell_ns_x1[dim] = max(cell_ns_x1[dim],comm->zone_d1[i].max1);
1212 print_ddzone(debug,2,i,j,&comm->zone_d2[i][j]);
1214 cell_ns_x0[dim] = min(cell_ns_x0[dim],comm->zone_d2[i][j].min0);
1215 cell_ns_x1[dim] = max(cell_ns_x1[dim],comm->zone_d2[i][j].max1);
1219 for(d=1; d<dd->ndim; d++)
1221 comm->cell_f_max0[d] = extr_s[d-1][0];
1222 comm->cell_f_min1[d] = extr_s[d-1][1];
1225 fprintf(debug,"Cell fraction d %d, max0 %f, min1 %f\n",
1226 d,comm->cell_f_max0[d],comm->cell_f_min1[d]);
1231 static void dd_collect_cg(gmx_domdec_t *dd,
1232 t_state *state_local)
1234 gmx_domdec_master_t *ma=NULL;
1235 int buf2[2],*ibuf,i,ncg_home=0,*cg=NULL,nat_home=0;
1238 if (state_local->ddp_count == dd->comm->master_cg_ddp_count)
1240 /* The master has the correct distribution */
1244 if (state_local->ddp_count == dd->ddp_count)
1246 ncg_home = dd->ncg_home;
1248 nat_home = dd->nat_home;
1250 else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
1252 cgs_gl = &dd->comm->cgs_gl;
1254 ncg_home = state_local->ncg_gl;
1255 cg = state_local->cg_gl;
1257 for(i=0; i<ncg_home; i++)
1259 nat_home += cgs_gl->index[cg[i]+1] - cgs_gl->index[cg[i]];
1264 gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
1267 buf2[0] = dd->ncg_home;
1268 buf2[1] = dd->nat_home;
1278 /* Collect the charge group and atom counts on the master */
1279 dd_gather(dd,2*sizeof(int),buf2,ibuf);
1284 for(i=0; i<dd->nnodes; i++)
1286 ma->ncg[i] = ma->ibuf[2*i];
1287 ma->nat[i] = ma->ibuf[2*i+1];
1288 ma->index[i+1] = ma->index[i] + ma->ncg[i];
1291 /* Make byte counts and indices */
1292 for(i=0; i<dd->nnodes; i++)
1294 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
1295 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
1299 fprintf(debug,"Initial charge group distribution: ");
1300 for(i=0; i<dd->nnodes; i++)
1301 fprintf(debug," %d",ma->ncg[i]);
1302 fprintf(debug,"\n");
1306 /* Collect the charge group indices on the master */
1308 dd->ncg_home*sizeof(int),dd->index_gl,
1309 DDMASTER(dd) ? ma->ibuf : NULL,
1310 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
1311 DDMASTER(dd) ? ma->cg : NULL);
1313 dd->comm->master_cg_ddp_count = state_local->ddp_count;
1316 static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
1319 gmx_domdec_master_t *ma;
1320 int n,i,c,a,nalloc=0;
1329 MPI_Send(lv,dd->nat_home*sizeof(rvec),MPI_BYTE,DDMASTERRANK(dd),
1330 dd->rank,dd->mpi_comm_all);
1333 /* Copy the master coordinates to the global array */
1334 cgs_gl = &dd->comm->cgs_gl;
1336 n = DDMASTERRANK(dd);
1338 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1340 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1342 copy_rvec(lv[a++],v[c]);
1346 for(n=0; n<dd->nnodes; n++)
1350 if (ma->nat[n] > nalloc)
1352 nalloc = over_alloc_dd(ma->nat[n]);
1356 MPI_Recv(buf,ma->nat[n]*sizeof(rvec),MPI_BYTE,DDRANK(dd,n),
1357 n,dd->mpi_comm_all,MPI_STATUS_IGNORE);
1360 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1362 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1364 copy_rvec(buf[a++],v[c]);
1373 static void get_commbuffer_counts(gmx_domdec_t *dd,
1374 int **counts,int **disps)
1376 gmx_domdec_master_t *ma;
1381 /* Make the rvec count and displacment arrays */
1383 *disps = ma->ibuf + dd->nnodes;
1384 for(n=0; n<dd->nnodes; n++)
1386 (*counts)[n] = ma->nat[n]*sizeof(rvec);
1387 (*disps)[n] = (n == 0 ? 0 : (*disps)[n-1] + (*counts)[n-1]);
1391 static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
1394 gmx_domdec_master_t *ma;
1395 int *rcounts=NULL,*disps=NULL;
1404 get_commbuffer_counts(dd,&rcounts,&disps);
1409 dd_gatherv(dd,dd->nat_home*sizeof(rvec),lv,rcounts,disps,buf);
1413 cgs_gl = &dd->comm->cgs_gl;
1416 for(n=0; n<dd->nnodes; n++)
1418 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1420 for(c=cgs_gl->index[ma->cg[i]]; c<cgs_gl->index[ma->cg[i]+1]; c++)
1422 copy_rvec(buf[a++],v[c]);
1429 void dd_collect_vec(gmx_domdec_t *dd,
1430 t_state *state_local,rvec *lv,rvec *v)
1432 gmx_domdec_master_t *ma;
1433 int n,i,c,a,nalloc=0;
1436 dd_collect_cg(dd,state_local);
1438 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1440 dd_collect_vec_sendrecv(dd,lv,v);
1444 dd_collect_vec_gatherv(dd,lv,v);
1449 void dd_collect_state(gmx_domdec_t *dd,
1450 t_state *state_local,t_state *state)
1454 nh = state->nhchainlength;
1458 for (i=0;i<efptNR;i++) {
1459 state->lambda[i] = state_local->lambda[i];
1461 state->fep_state = state_local->fep_state;
1462 state->veta = state_local->veta;
1463 state->vol0 = state_local->vol0;
1464 copy_mat(state_local->box,state->box);
1465 copy_mat(state_local->boxv,state->boxv);
1466 copy_mat(state_local->svir_prev,state->svir_prev);
1467 copy_mat(state_local->fvir_prev,state->fvir_prev);
1468 copy_mat(state_local->pres_prev,state->pres_prev);
1471 for(i=0; i<state_local->ngtc; i++)
1473 for(j=0; j<nh; j++) {
1474 state->nosehoover_xi[i*nh+j] = state_local->nosehoover_xi[i*nh+j];
1475 state->nosehoover_vxi[i*nh+j] = state_local->nosehoover_vxi[i*nh+j];
1477 state->therm_integral[i] = state_local->therm_integral[i];
1479 for(i=0; i<state_local->nnhpres; i++)
1481 for(j=0; j<nh; j++) {
1482 state->nhpres_xi[i*nh+j] = state_local->nhpres_xi[i*nh+j];
1483 state->nhpres_vxi[i*nh+j] = state_local->nhpres_vxi[i*nh+j];
1487 for(est=0; est<estNR; est++)
1489 if (EST_DISTR(est) && (state_local->flags & (1<<est)))
1493 dd_collect_vec(dd,state_local,state_local->x,state->x);
1496 dd_collect_vec(dd,state_local,state_local->v,state->v);
1499 dd_collect_vec(dd,state_local,state_local->sd_X,state->sd_X);
1502 dd_collect_vec(dd,state_local,state_local->cg_p,state->cg_p);
1505 if (state->nrngi == 1)
1509 for(i=0; i<state_local->nrng; i++)
1511 state->ld_rng[i] = state_local->ld_rng[i];
1517 dd_gather(dd,state_local->nrng*sizeof(state->ld_rng[0]),
1518 state_local->ld_rng,state->ld_rng);
1522 if (state->nrngi == 1)
1526 state->ld_rngi[0] = state_local->ld_rngi[0];
1531 dd_gather(dd,sizeof(state->ld_rngi[0]),
1532 state_local->ld_rngi,state->ld_rngi);
1535 case estDISRE_INITF:
1536 case estDISRE_RM3TAV:
1537 case estORIRE_INITF:
1541 gmx_incons("Unknown state entry encountered in dd_collect_state");
1547 static void dd_realloc_fr_cg(t_forcerec *fr,int nalloc)
1551 fprintf(debug,"Reallocating forcerec: currently %d, required %d, allocating %d\n",fr->cg_nalloc,nalloc,over_alloc_dd(nalloc));
1553 fr->cg_nalloc = over_alloc_dd(nalloc);
1554 srenew(fr->cg_cm,fr->cg_nalloc);
1555 srenew(fr->cginfo,fr->cg_nalloc);
1558 static void dd_realloc_state(t_state *state,rvec **f,int nalloc)
1564 fprintf(debug,"Reallocating state: currently %d, required %d, allocating %d\n",state->nalloc,nalloc,over_alloc_dd(nalloc));
1567 state->nalloc = over_alloc_dd(nalloc);
1569 for(est=0; est<estNR; est++)
1571 if (EST_DISTR(est) && (state->flags & (1<<est)))
1575 srenew(state->x,state->nalloc);
1578 srenew(state->v,state->nalloc);
1581 srenew(state->sd_X,state->nalloc);
1584 srenew(state->cg_p,state->nalloc);
1588 case estDISRE_INITF:
1589 case estDISRE_RM3TAV:
1590 case estORIRE_INITF:
1592 /* No reallocation required */
1595 gmx_incons("Unknown state entry encountered in dd_realloc_state");
1602 srenew(*f,state->nalloc);
1606 static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd,t_block *cgs,
1609 gmx_domdec_master_t *ma;
1610 int n,i,c,a,nalloc=0;
1617 for(n=0; n<dd->nnodes; n++)
1621 if (ma->nat[n] > nalloc)
1623 nalloc = over_alloc_dd(ma->nat[n]);
1626 /* Use lv as a temporary buffer */
1628 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1630 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1632 copy_rvec(v[c],buf[a++]);
1635 if (a != ma->nat[n])
1637 gmx_fatal(FARGS,"Internal error a (%d) != nat (%d)",
1642 MPI_Send(buf,ma->nat[n]*sizeof(rvec),MPI_BYTE,
1643 DDRANK(dd,n),n,dd->mpi_comm_all);
1648 n = DDMASTERRANK(dd);
1650 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1652 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1654 copy_rvec(v[c],lv[a++]);
1661 MPI_Recv(lv,dd->nat_home*sizeof(rvec),MPI_BYTE,DDMASTERRANK(dd),
1662 MPI_ANY_TAG,dd->mpi_comm_all,MPI_STATUS_IGNORE);
1667 static void dd_distribute_vec_scatterv(gmx_domdec_t *dd,t_block *cgs,
1670 gmx_domdec_master_t *ma;
1671 int *scounts=NULL,*disps=NULL;
1672 int n,i,c,a,nalloc=0;
1679 get_commbuffer_counts(dd,&scounts,&disps);
1683 for(n=0; n<dd->nnodes; n++)
1685 for(i=ma->index[n]; i<ma->index[n+1]; i++)
1687 for(c=cgs->index[ma->cg[i]]; c<cgs->index[ma->cg[i]+1]; c++)
1689 copy_rvec(v[c],buf[a++]);
1695 dd_scatterv(dd,scounts,disps,buf,dd->nat_home*sizeof(rvec),lv);
1698 static void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv)
1700 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
1702 dd_distribute_vec_sendrecv(dd,cgs,v,lv);
1706 dd_distribute_vec_scatterv(dd,cgs,v,lv);
1710 static void dd_distribute_state(gmx_domdec_t *dd,t_block *cgs,
1711 t_state *state,t_state *state_local,
1716 nh = state->nhchainlength;
1720 for(i=0;i<efptNR;i++)
1722 state_local->lambda[i] = state->lambda[i];
1724 state_local->fep_state = state->fep_state;
1725 state_local->veta = state->veta;
1726 state_local->vol0 = state->vol0;
1727 copy_mat(state->box,state_local->box);
1728 copy_mat(state->box_rel,state_local->box_rel);
1729 copy_mat(state->boxv,state_local->boxv);
1730 copy_mat(state->svir_prev,state_local->svir_prev);
1731 copy_mat(state->fvir_prev,state_local->fvir_prev);
1732 for(i=0; i<state_local->ngtc; i++)
1734 for(j=0; j<nh; j++) {
1735 state_local->nosehoover_xi[i*nh+j] = state->nosehoover_xi[i*nh+j];
1736 state_local->nosehoover_vxi[i*nh+j] = state->nosehoover_vxi[i*nh+j];
1738 state_local->therm_integral[i] = state->therm_integral[i];
1740 for(i=0; i<state_local->nnhpres; i++)
1742 for(j=0; j<nh; j++) {
1743 state_local->nhpres_xi[i*nh+j] = state->nhpres_xi[i*nh+j];
1744 state_local->nhpres_vxi[i*nh+j] = state->nhpres_vxi[i*nh+j];
1748 dd_bcast(dd,((efptNR)*sizeof(real)),state_local->lambda);
1749 dd_bcast(dd,sizeof(int),&state_local->fep_state);
1750 dd_bcast(dd,sizeof(real),&state_local->veta);
1751 dd_bcast(dd,sizeof(real),&state_local->vol0);
1752 dd_bcast(dd,sizeof(state_local->box),state_local->box);
1753 dd_bcast(dd,sizeof(state_local->box_rel),state_local->box_rel);
1754 dd_bcast(dd,sizeof(state_local->boxv),state_local->boxv);
1755 dd_bcast(dd,sizeof(state_local->svir_prev),state_local->svir_prev);
1756 dd_bcast(dd,sizeof(state_local->fvir_prev),state_local->fvir_prev);
1757 dd_bcast(dd,((state_local->ngtc*nh)*sizeof(double)),state_local->nosehoover_xi);
1758 dd_bcast(dd,((state_local->ngtc*nh)*sizeof(double)),state_local->nosehoover_vxi);
1759 dd_bcast(dd,state_local->ngtc*sizeof(double),state_local->therm_integral);
1760 dd_bcast(dd,((state_local->nnhpres*nh)*sizeof(double)),state_local->nhpres_xi);
1761 dd_bcast(dd,((state_local->nnhpres*nh)*sizeof(double)),state_local->nhpres_vxi);
1763 if (dd->nat_home > state_local->nalloc)
1765 dd_realloc_state(state_local,f,dd->nat_home);
1767 for(i=0; i<estNR; i++)
1769 if (EST_DISTR(i) && (state_local->flags & (1<<i)))
1773 dd_distribute_vec(dd,cgs,state->x,state_local->x);
1776 dd_distribute_vec(dd,cgs,state->v,state_local->v);
1779 dd_distribute_vec(dd,cgs,state->sd_X,state_local->sd_X);
1782 dd_distribute_vec(dd,cgs,state->cg_p,state_local->cg_p);
1785 if (state->nrngi == 1)
1788 state_local->nrng*sizeof(state_local->ld_rng[0]),
1789 state->ld_rng,state_local->ld_rng);
1794 state_local->nrng*sizeof(state_local->ld_rng[0]),
1795 state->ld_rng,state_local->ld_rng);
1799 if (state->nrngi == 1)
1801 dd_bcastc(dd,sizeof(state_local->ld_rngi[0]),
1802 state->ld_rngi,state_local->ld_rngi);
1806 dd_scatter(dd,sizeof(state_local->ld_rngi[0]),
1807 state->ld_rngi,state_local->ld_rngi);
1810 case estDISRE_INITF:
1811 case estDISRE_RM3TAV:
1812 case estORIRE_INITF:
1814 /* Not implemented yet */
1817 gmx_incons("Unknown state entry encountered in dd_distribute_state");
1823 static char dim2char(int dim)
1829 case XX: c = 'X'; break;
1830 case YY: c = 'Y'; break;
1831 case ZZ: c = 'Z'; break;
1832 default: gmx_fatal(FARGS,"Unknown dim %d",dim);
1838 static void write_dd_grid_pdb(const char *fn,gmx_large_int_t step,
1839 gmx_domdec_t *dd,matrix box,gmx_ddbox_t *ddbox)
1841 rvec grid_s[2],*grid_r=NULL,cx,r;
1842 char fname[STRLEN],format[STRLEN],buf[22];
1848 copy_rvec(dd->comm->cell_x0,grid_s[0]);
1849 copy_rvec(dd->comm->cell_x1,grid_s[1]);
1853 snew(grid_r,2*dd->nnodes);
1856 dd_gather(dd,2*sizeof(rvec),grid_s[0],DDMASTER(dd) ? grid_r[0] : NULL);
1860 for(d=0; d<DIM; d++)
1862 for(i=0; i<DIM; i++)
1870 if (d < ddbox->npbcdim && dd->nc[d] > 1)
1872 tric[d][i] = box[i][d]/box[i][i];
1881 sprintf(fname,"%s_%s.pdb",fn,gmx_step_str(step,buf));
1882 sprintf(format,"%s%s\n",get_pdbformat(),"%6.2f%6.2f");
1883 out = gmx_fio_fopen(fname,"w");
1884 gmx_write_pdb_box(out,dd->bScrewPBC ? epbcSCREW : epbcXYZ,box);
1886 for(i=0; i<dd->nnodes; i++)
1888 vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1889 for(d=0; d<DIM; d++)
1891 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1899 cx[XX] = grid_r[i*2+x][XX];
1900 cx[YY] = grid_r[i*2+y][YY];
1901 cx[ZZ] = grid_r[i*2+z][ZZ];
1903 fprintf(out,format,"ATOM",a++,"CA","GLY",' ',1+i,
1904 10*r[XX],10*r[YY],10*r[ZZ],1.0,vol);
1908 for(d=0; d<DIM; d++)
1914 case 0: y = 1 + i*8 + 2*x; break;
1915 case 1: y = 1 + i*8 + 2*x - (x % 2); break;
1916 case 2: y = 1 + i*8 + x; break;
1918 fprintf(out,"%6s%5d%5d\n","CONECT",y,y+(1<<d));
1922 gmx_fio_fclose(out);
1927 void write_dd_pdb(const char *fn,gmx_large_int_t step,const char *title,
1928 gmx_mtop_t *mtop,t_commrec *cr,
1929 int natoms,rvec x[],matrix box)
1931 char fname[STRLEN],format[STRLEN],format4[STRLEN],buf[22];
1934 char *atomname,*resname;
1941 natoms = dd->comm->nat[ddnatVSITE];
1944 sprintf(fname,"%s_%s_n%d.pdb",fn,gmx_step_str(step,buf),cr->sim_nodeid);
1946 sprintf(format,"%s%s\n",get_pdbformat(),"%6.2f%6.2f");
1947 sprintf(format4,"%s%s\n",get_pdbformat4(),"%6.2f%6.2f");
1949 out = gmx_fio_fopen(fname,"w");
1951 fprintf(out,"TITLE %s\n",title);
1952 gmx_write_pdb_box(out,dd->bScrewPBC ? epbcSCREW : epbcXYZ,box);
1953 for(i=0; i<natoms; i++)
1955 ii = dd->gatindex[i];
1956 gmx_mtop_atominfo_global(mtop,ii,&atomname,&resnr,&resname);
1957 if (i < dd->comm->nat[ddnatZONE])
1960 while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
1966 else if (i < dd->comm->nat[ddnatVSITE])
1968 b = dd->comm->zones.n;
1972 b = dd->comm->zones.n + 1;
1974 fprintf(out,strlen(atomname)<4 ? format : format4,
1975 "ATOM",(ii+1)%100000,
1976 atomname,resname,' ',resnr%10000,' ',
1977 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],1.0,b);
1979 fprintf(out,"TER\n");
1981 gmx_fio_fclose(out);
1984 real dd_cutoff_mbody(gmx_domdec_t *dd)
1986 gmx_domdec_comm_t *comm;
1993 if (comm->bInterCGBondeds)
1995 if (comm->cutoff_mbody > 0)
1997 r = comm->cutoff_mbody;
2001 /* cutoff_mbody=0 means we do not have DLB */
2002 r = comm->cellsize_min[dd->dim[0]];
2003 for(di=1; di<dd->ndim; di++)
2005 r = min(r,comm->cellsize_min[dd->dim[di]]);
2007 if (comm->bBondComm)
2009 r = max(r,comm->cutoff_mbody);
2013 r = min(r,comm->cutoff);
2021 real dd_cutoff_twobody(gmx_domdec_t *dd)
2025 r_mb = dd_cutoff_mbody(dd);
2027 return max(dd->comm->cutoff,r_mb);
2031 static void dd_cart_coord2pmecoord(gmx_domdec_t *dd,ivec coord,ivec coord_pme)
2035 nc = dd->nc[dd->comm->cartpmedim];
2036 ntot = dd->comm->ntot[dd->comm->cartpmedim];
2037 copy_ivec(coord,coord_pme);
2038 coord_pme[dd->comm->cartpmedim] =
2039 nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
2042 static int low_ddindex2pmeindex(int ndd,int npme,int ddindex)
2044 /* Here we assign a PME node to communicate with this DD node
2045 * by assuming that the major index of both is x.
2046 * We add cr->npmenodes/2 to obtain an even distribution.
2048 return (ddindex*npme + npme/2)/ndd;
2051 static int ddindex2pmeindex(const gmx_domdec_t *dd,int ddindex)
2053 return low_ddindex2pmeindex(dd->nnodes,dd->comm->npmenodes,ddindex);
2056 static int cr_ddindex2pmeindex(const t_commrec *cr,int ddindex)
2058 return low_ddindex2pmeindex(cr->dd->nnodes,cr->npmenodes,ddindex);
2061 static int *dd_pmenodes(t_commrec *cr)
2066 snew(pmenodes,cr->npmenodes);
2068 for(i=0; i<cr->dd->nnodes; i++) {
2069 p0 = cr_ddindex2pmeindex(cr,i);
2070 p1 = cr_ddindex2pmeindex(cr,i+1);
2071 if (i+1 == cr->dd->nnodes || p1 > p0) {
2073 fprintf(debug,"pmenode[%d] = %d\n",n,i+1+n);
2074 pmenodes[n] = i + 1 + n;
2082 static int gmx_ddcoord2pmeindex(t_commrec *cr,int x,int y,int z)
2085 ivec coords,coords_pme,nc;
2090 if (dd->comm->bCartesian) {
2091 gmx_ddindex2xyz(dd->nc,ddindex,coords);
2092 dd_coords2pmecoords(dd,coords,coords_pme);
2093 copy_ivec(dd->ntot,nc);
2094 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2095 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
2097 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
2099 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
2105 slab = ddindex2pmeindex(dd,dd_index(dd->nc,coords));
2110 static int ddcoord2simnodeid(t_commrec *cr,int x,int y,int z)
2112 gmx_domdec_comm_t *comm;
2114 int ddindex,nodeid=-1;
2116 comm = cr->dd->comm;
2121 if (comm->bCartesianPP_PME)
2124 MPI_Cart_rank(cr->mpi_comm_mysim,coords,&nodeid);
2129 ddindex = dd_index(cr->dd->nc,coords);
2130 if (comm->bCartesianPP)
2132 nodeid = comm->ddindex2simnodeid[ddindex];
2138 nodeid = ddindex + gmx_ddcoord2pmeindex(cr,x,y,z);
2150 static int dd_simnode2pmenode(t_commrec *cr,int sim_nodeid)
2153 gmx_domdec_comm_t *comm;
2154 ivec coord,coord_pme;
2161 /* This assumes a uniform x domain decomposition grid cell size */
2162 if (comm->bCartesianPP_PME)
2165 MPI_Cart_coords(cr->mpi_comm_mysim,sim_nodeid,DIM,coord);
2166 if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
2168 /* This is a PP node */
2169 dd_cart_coord2pmecoord(dd,coord,coord_pme);
2170 MPI_Cart_rank(cr->mpi_comm_mysim,coord_pme,&pmenode);
2174 else if (comm->bCartesianPP)
2176 if (sim_nodeid < dd->nnodes)
2178 pmenode = dd->nnodes + ddindex2pmeindex(dd,sim_nodeid);
2183 /* This assumes DD cells with identical x coordinates
2184 * are numbered sequentially.
2186 if (dd->comm->pmenodes == NULL)
2188 if (sim_nodeid < dd->nnodes)
2190 /* The DD index equals the nodeid */
2191 pmenode = dd->nnodes + ddindex2pmeindex(dd,sim_nodeid);
2197 while (sim_nodeid > dd->comm->pmenodes[i])
2201 if (sim_nodeid < dd->comm->pmenodes[i])
2203 pmenode = dd->comm->pmenodes[i];
2211 gmx_bool gmx_pmeonlynode(t_commrec *cr,int sim_nodeid)
2213 gmx_bool bPMEOnlyNode;
2215 if (DOMAINDECOMP(cr))
2217 bPMEOnlyNode = (dd_simnode2pmenode(cr,sim_nodeid) == -1);
2221 bPMEOnlyNode = FALSE;
2224 return bPMEOnlyNode;
2227 void get_pme_ddnodes(t_commrec *cr,int pmenodeid,
2228 int *nmy_ddnodes,int **my_ddnodes,int *node_peer)
2232 ivec coord,coord_pme;
2236 snew(*my_ddnodes,(dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
2239 for(x=0; x<dd->nc[XX]; x++)
2241 for(y=0; y<dd->nc[YY]; y++)
2243 for(z=0; z<dd->nc[ZZ]; z++)
2245 if (dd->comm->bCartesianPP_PME)
2250 dd_cart_coord2pmecoord(dd,coord,coord_pme);
2251 if (dd->ci[XX] == coord_pme[XX] &&
2252 dd->ci[YY] == coord_pme[YY] &&
2253 dd->ci[ZZ] == coord_pme[ZZ])
2254 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr,x,y,z);
2258 /* The slab corresponds to the nodeid in the PME group */
2259 if (gmx_ddcoord2pmeindex(cr,x,y,z) == pmenodeid)
2261 (*my_ddnodes)[(*nmy_ddnodes)++] = ddcoord2simnodeid(cr,x,y,z);
2268 /* The last PP-only node is the peer node */
2269 *node_peer = (*my_ddnodes)[*nmy_ddnodes-1];
2273 fprintf(debug,"Receive coordinates from PP nodes:");
2274 for(x=0; x<*nmy_ddnodes; x++)
2276 fprintf(debug," %d",(*my_ddnodes)[x]);
2278 fprintf(debug,"\n");
2282 static gmx_bool receive_vir_ener(t_commrec *cr)
2284 gmx_domdec_comm_t *comm;
2285 int pmenode,coords[DIM],rank;
2289 if (cr->npmenodes < cr->dd->nnodes)
2291 comm = cr->dd->comm;
2292 if (comm->bCartesianPP_PME)
2294 pmenode = dd_simnode2pmenode(cr,cr->sim_nodeid);
2296 MPI_Cart_coords(cr->mpi_comm_mysim,cr->sim_nodeid,DIM,coords);
2297 coords[comm->cartpmedim]++;
2298 if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
2300 MPI_Cart_rank(cr->mpi_comm_mysim,coords,&rank);
2301 if (dd_simnode2pmenode(cr,rank) == pmenode)
2303 /* This is not the last PP node for pmenode */
2311 pmenode = dd_simnode2pmenode(cr,cr->sim_nodeid);
2312 if (cr->sim_nodeid+1 < cr->nnodes &&
2313 dd_simnode2pmenode(cr,cr->sim_nodeid+1) == pmenode)
2315 /* This is not the last PP node for pmenode */
2324 static void set_zones_ncg_home(gmx_domdec_t *dd)
2326 gmx_domdec_zones_t *zones;
2329 zones = &dd->comm->zones;
2331 zones->cg_range[0] = 0;
2332 for(i=1; i<zones->n+1; i++)
2334 zones->cg_range[i] = dd->ncg_home;
2338 static void rebuild_cgindex(gmx_domdec_t *dd,int *gcgs_index,t_state *state)
2340 int nat,i,*ind,*dd_cg_gl,*cgindex,cg_gl;
2343 dd_cg_gl = dd->index_gl;
2344 cgindex = dd->cgindex;
2347 for(i=0; i<state->ncg_gl; i++)
2351 dd_cg_gl[i] = cg_gl;
2352 nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
2356 dd->ncg_home = state->ncg_gl;
2359 set_zones_ncg_home(dd);
2362 static int ddcginfo(const cginfo_mb_t *cginfo_mb,int cg)
2364 while (cg >= cginfo_mb->cg_end)
2369 return cginfo_mb->cginfo[(cg - cginfo_mb->cg_start) % cginfo_mb->cg_mod];
2372 static void dd_set_cginfo(int *index_gl,int cg0,int cg1,
2373 t_forcerec *fr,char *bLocalCG)
2375 cginfo_mb_t *cginfo_mb;
2381 cginfo_mb = fr->cginfo_mb;
2382 cginfo = fr->cginfo;
2384 for(cg=cg0; cg<cg1; cg++)
2386 cginfo[cg] = ddcginfo(cginfo_mb,index_gl[cg]);
2390 if (bLocalCG != NULL)
2392 for(cg=cg0; cg<cg1; cg++)
2394 bLocalCG[index_gl[cg]] = TRUE;
2399 static void make_dd_indices(gmx_domdec_t *dd,int *gcgs_index,int cg_start)
2401 int nzone,zone,zone1,cg0,cg,cg_gl,a,a_gl;
2402 int *zone2cg,*zone_ncg1,*index_gl,*gatindex;
2406 bLocalCG = dd->comm->bLocalCG;
2408 if (dd->nat_tot > dd->gatindex_nalloc)
2410 dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
2411 srenew(dd->gatindex,dd->gatindex_nalloc);
2414 nzone = dd->comm->zones.n;
2415 zone2cg = dd->comm->zones.cg_range;
2416 zone_ncg1 = dd->comm->zone_ncg1;
2417 index_gl = dd->index_gl;
2418 gatindex = dd->gatindex;
2420 if (zone2cg[1] != dd->ncg_home)
2422 gmx_incons("dd->ncg_zone is not up to date");
2425 /* Make the local to global and global to local atom index */
2426 a = dd->cgindex[cg_start];
2427 for(zone=0; zone<nzone; zone++)
2435 cg0 = zone2cg[zone];
2437 for(cg=cg0; cg<zone2cg[zone+1]; cg++)
2440 if (cg - cg0 >= zone_ncg1[zone])
2442 /* Signal that this cg is from more than one zone away */
2445 cg_gl = index_gl[cg];
2446 for(a_gl=gcgs_index[cg_gl]; a_gl<gcgs_index[cg_gl+1]; a_gl++)
2449 ga2la_set(dd->ga2la,a_gl,a,zone1);
2456 static int check_bLocalCG(gmx_domdec_t *dd,int ncg_sys,const char *bLocalCG,
2462 if (bLocalCG == NULL)
2466 for(i=0; i<dd->ncg_tot; i++)
2468 if (!bLocalCG[dd->index_gl[i]])
2471 "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);
2476 for(i=0; i<ncg_sys; i++)
2483 if (ngl != dd->ncg_tot)
2485 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);
2492 static void check_index_consistency(gmx_domdec_t *dd,
2493 int natoms_sys,int ncg_sys,
2496 int nerr,ngl,i,a,cell;
2501 if (dd->comm->DD_debug > 1)
2503 snew(have,natoms_sys);
2504 for(a=0; a<dd->nat_tot; a++)
2506 if (have[dd->gatindex[a]] > 0)
2508 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);
2512 have[dd->gatindex[a]] = a + 1;
2518 snew(have,dd->nat_tot);
2521 for(i=0; i<natoms_sys; i++)
2523 if (ga2la_get(dd->ga2la,i,&a,&cell))
2525 if (a >= dd->nat_tot)
2527 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);
2533 if (dd->gatindex[a] != i)
2535 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);
2542 if (ngl != dd->nat_tot)
2545 "DD node %d, %s: %d global atom indices, %d local atoms\n",
2546 dd->rank,where,ngl,dd->nat_tot);
2548 for(a=0; a<dd->nat_tot; a++)
2553 "DD node %d, %s: local atom %d, global %d has no global index\n",
2554 dd->rank,where,a+1,dd->gatindex[a]+1);
2559 nerr += check_bLocalCG(dd,ncg_sys,dd->comm->bLocalCG,where);
2562 gmx_fatal(FARGS,"DD node %d, %s: %d atom/cg index inconsistencies",
2563 dd->rank,where,nerr);
2567 static void clear_dd_indices(gmx_domdec_t *dd,int cg_start,int a_start)
2574 /* Clear the whole list without searching */
2575 ga2la_clear(dd->ga2la);
2579 for(i=a_start; i<dd->nat_tot; i++)
2581 ga2la_del(dd->ga2la,dd->gatindex[i]);
2585 bLocalCG = dd->comm->bLocalCG;
2588 for(i=cg_start; i<dd->ncg_tot; i++)
2590 bLocalCG[dd->index_gl[i]] = FALSE;
2594 dd_clear_local_vsite_indices(dd);
2596 if (dd->constraints)
2598 dd_clear_local_constraint_indices(dd);
2602 static real grid_jump_limit(gmx_domdec_comm_t *comm,int dim_ind)
2604 real grid_jump_limit;
2606 /* The distance between the boundaries of cells at distance
2607 * x+-1,y+-1 or y+-1,z+-1 is limited by the cut-off restrictions
2608 * and by the fact that cells should not be shifted by more than
2609 * half their size, such that cg's only shift by one cell
2610 * at redecomposition.
2612 grid_jump_limit = comm->cellsize_limit;
2613 if (!comm->bVacDLBNoLimit)
2615 grid_jump_limit = max(grid_jump_limit,
2616 comm->cutoff/comm->cd[dim_ind].np);
2619 return grid_jump_limit;
2622 static void check_grid_jump(gmx_large_int_t step,gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
2624 gmx_domdec_comm_t *comm;
2630 for(d=1; d<dd->ndim; d++)
2633 limit = grid_jump_limit(comm,d);
2634 bfac = ddbox->box_size[dim];
2635 if (ddbox->tric_dir[dim])
2637 bfac *= ddbox->skew_fac[dim];
2639 if ((comm->cell_f1[d] - comm->cell_f_max0[d])*bfac < limit ||
2640 (comm->cell_f0[d] - comm->cell_f_min1[d])*bfac > -limit)
2643 gmx_fatal(FARGS,"Step %s: The domain decomposition grid has shifted too much in the %c-direction around cell %d %d %d\n",
2644 gmx_step_str(step,buf),
2645 dim2char(dim),dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
2650 static int dd_load_count(gmx_domdec_comm_t *comm)
2652 return (comm->eFlop ? comm->flop_n : comm->cycl_n[ddCyclF]);
2655 static float dd_force_load(gmx_domdec_comm_t *comm)
2662 if (comm->eFlop > 1)
2664 load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
2669 load = comm->cycl[ddCyclF];
2670 if (comm->cycl_n[ddCyclF] > 1)
2672 /* Subtract the maximum of the last n cycle counts
2673 * to get rid of possible high counts due to other soures,
2674 * for instance system activity, that would otherwise
2675 * affect the dynamic load balancing.
2677 load -= comm->cycl_max[ddCyclF];
2684 static void set_slb_pme_dim_f(gmx_domdec_t *dd,int dim,real **dim_f)
2686 gmx_domdec_comm_t *comm;
2691 snew(*dim_f,dd->nc[dim]+1);
2693 for(i=1; i<dd->nc[dim]; i++)
2695 if (comm->slb_frac[dim])
2697 (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
2701 (*dim_f)[i] = (real)i/(real)dd->nc[dim];
2704 (*dim_f)[dd->nc[dim]] = 1;
2707 static void init_ddpme(gmx_domdec_t *dd,gmx_ddpme_t *ddpme,int dimind)
2709 int pmeindex,slab,nso,i;
2712 if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
2718 ddpme->dim = dimind;
2720 ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
2722 ddpme->nslab = (ddpme->dim == 0 ?
2723 dd->comm->npmenodes_x :
2724 dd->comm->npmenodes_y);
2726 if (ddpme->nslab <= 1)
2731 nso = dd->comm->npmenodes/ddpme->nslab;
2732 /* Determine for each PME slab the PP location range for dimension dim */
2733 snew(ddpme->pp_min,ddpme->nslab);
2734 snew(ddpme->pp_max,ddpme->nslab);
2735 for(slab=0; slab<ddpme->nslab; slab++) {
2736 ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
2737 ddpme->pp_max[slab] = 0;
2739 for(i=0; i<dd->nnodes; i++) {
2740 ddindex2xyz(dd->nc,i,xyz);
2741 /* For y only use our y/z slab.
2742 * This assumes that the PME x grid size matches the DD grid size.
2744 if (dimind == 0 || xyz[XX] == dd->ci[XX]) {
2745 pmeindex = ddindex2pmeindex(dd,i);
2747 slab = pmeindex/nso;
2749 slab = pmeindex % ddpme->nslab;
2751 ddpme->pp_min[slab] = min(ddpme->pp_min[slab],xyz[dimind]);
2752 ddpme->pp_max[slab] = max(ddpme->pp_max[slab],xyz[dimind]);
2756 set_slb_pme_dim_f(dd,ddpme->dim,&ddpme->slb_dim_f);
2759 int dd_pme_maxshift_x(gmx_domdec_t *dd)
2761 if (dd->comm->ddpme[0].dim == XX)
2763 return dd->comm->ddpme[0].maxshift;
2771 int dd_pme_maxshift_y(gmx_domdec_t *dd)
2773 if (dd->comm->ddpme[0].dim == YY)
2775 return dd->comm->ddpme[0].maxshift;
2777 else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
2779 return dd->comm->ddpme[1].maxshift;
2787 static void set_pme_maxshift(gmx_domdec_t *dd,gmx_ddpme_t *ddpme,
2788 gmx_bool bUniform,gmx_ddbox_t *ddbox,real *cell_f)
2790 gmx_domdec_comm_t *comm;
2793 real range,pme_boundary;
2797 nc = dd->nc[ddpme->dim];
2800 if (!ddpme->dim_match)
2802 /* PP decomposition is not along dim: the worst situation */
2805 else if (ns <= 3 || (bUniform && ns == nc))
2807 /* The optimal situation */
2812 /* We need to check for all pme nodes which nodes they
2813 * could possibly need to communicate with.
2815 xmin = ddpme->pp_min;
2816 xmax = ddpme->pp_max;
2817 /* Allow for atoms to be maximally 2/3 times the cut-off
2818 * out of their DD cell. This is a reasonable balance between
2819 * between performance and support for most charge-group/cut-off
2822 range = 2.0/3.0*comm->cutoff/ddbox->box_size[ddpme->dim];
2823 /* Avoid extra communication when we are exactly at a boundary */
2829 /* PME slab s spreads atoms between box frac. s/ns and (s+1)/ns */
2830 pme_boundary = (real)s/ns;
2833 cell_f[xmax[s-(sh+1) ]+1] + range > pme_boundary) ||
2835 cell_f[xmax[s-(sh+1)+ns]+1] - 1 + range > pme_boundary)))
2839 pme_boundary = (real)(s+1)/ns;
2842 cell_f[xmin[s+(sh+1) ] ] - range < pme_boundary) ||
2844 cell_f[xmin[s+(sh+1)-ns] ] + 1 - range < pme_boundary)))
2851 ddpme->maxshift = sh;
2855 fprintf(debug,"PME slab communication range for dim %d is %d\n",
2856 ddpme->dim,ddpme->maxshift);
2860 static void check_box_size(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
2864 for(d=0; d<dd->ndim; d++)
2867 if (dim < ddbox->nboundeddim &&
2868 ddbox->box_size[dim]*ddbox->skew_fac[dim] <
2869 dd->nc[dim]*dd->comm->cellsize_limit*DD_CELL_MARGIN)
2871 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",
2872 dim2char(dim),ddbox->box_size[dim],ddbox->skew_fac[dim],
2873 dd->nc[dim],dd->comm->cellsize_limit);
2878 static void set_dd_cell_sizes_slb(gmx_domdec_t *dd,gmx_ddbox_t *ddbox,
2879 gmx_bool bMaster,ivec npulse)
2881 gmx_domdec_comm_t *comm;
2884 real *cell_x,cell_dx,cellsize;
2888 for(d=0; d<DIM; d++)
2890 cellsize_min[d] = ddbox->box_size[d]*ddbox->skew_fac[d];
2892 if (dd->nc[d] == 1 || comm->slb_frac[d] == NULL)
2895 cell_dx = ddbox->box_size[d]/dd->nc[d];
2898 for(j=0; j<dd->nc[d]+1; j++)
2900 dd->ma->cell_x[d][j] = ddbox->box0[d] + j*cell_dx;
2905 comm->cell_x0[d] = ddbox->box0[d] + (dd->ci[d] )*cell_dx;
2906 comm->cell_x1[d] = ddbox->box0[d] + (dd->ci[d]+1)*cell_dx;
2908 cellsize = cell_dx*ddbox->skew_fac[d];
2909 while (cellsize*npulse[d] < comm->cutoff && npulse[d] < dd->nc[d]-1)
2913 cellsize_min[d] = cellsize;
2917 /* Statically load balanced grid */
2918 /* Also when we are not doing a master distribution we determine
2919 * all cell borders in a loop to obtain identical values
2920 * to the master distribution case and to determine npulse.
2924 cell_x = dd->ma->cell_x[d];
2928 snew(cell_x,dd->nc[d]+1);
2930 cell_x[0] = ddbox->box0[d];
2931 for(j=0; j<dd->nc[d]; j++)
2933 cell_dx = ddbox->box_size[d]*comm->slb_frac[d][j];
2934 cell_x[j+1] = cell_x[j] + cell_dx;
2935 cellsize = cell_dx*ddbox->skew_fac[d];
2936 while (cellsize*npulse[d] < comm->cutoff &&
2937 npulse[d] < dd->nc[d]-1)
2941 cellsize_min[d] = min(cellsize_min[d],cellsize);
2945 comm->cell_x0[d] = cell_x[dd->ci[d]];
2946 comm->cell_x1[d] = cell_x[dd->ci[d]+1];
2950 /* The following limitation is to avoid that a cell would receive
2951 * some of its own home charge groups back over the periodic boundary.
2952 * Double charge groups cause trouble with the global indices.
2954 if (d < ddbox->npbcdim &&
2955 dd->nc[d] > 1 && npulse[d] >= dd->nc[d])
2957 gmx_fatal_collective(FARGS,NULL,dd,
2958 "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",
2959 dim2char(d),ddbox->box_size[d],ddbox->skew_fac[d],
2961 dd->nc[d],dd->nc[d],
2962 dd->nnodes > dd->nc[d] ? "cells" : "processors");
2966 if (!comm->bDynLoadBal)
2968 copy_rvec(cellsize_min,comm->cellsize_min);
2971 for(d=0; d<comm->npmedecompdim; d++)
2973 set_pme_maxshift(dd,&comm->ddpme[d],
2974 comm->slb_frac[dd->dim[d]]==NULL,ddbox,
2975 comm->ddpme[d].slb_dim_f);
2980 static void dd_cell_sizes_dlb_root_enforce_limits(gmx_domdec_t *dd,
2981 int d,int dim,gmx_domdec_root_t *root,
2983 gmx_bool bUniform,gmx_large_int_t step, real cellsize_limit_f, int range[])
2985 gmx_domdec_comm_t *comm;
2986 int ncd,i,j,nmin,nmin_old;
2987 gmx_bool bLimLo,bLimHi;
2989 real fac,halfway,cellsize_limit_f_i,region_size;
2990 gmx_bool bPBC,bLastHi=FALSE;
2991 int nrange[]={range[0],range[1]};
2993 region_size= root->cell_f[range[1]]-root->cell_f[range[0]];
2999 bPBC = (dim < ddbox->npbcdim);
3001 cell_size = root->buf_ncd;
3005 fprintf(debug,"enforce_limits: %d %d\n",range[0],range[1]);
3008 /* First we need to check if the scaling does not make cells
3009 * smaller than the smallest allowed size.
3010 * We need to do this iteratively, since if a cell is too small,
3011 * it needs to be enlarged, which makes all the other cells smaller,
3012 * which could in turn make another cell smaller than allowed.
3014 for(i=range[0]; i<range[1]; i++)
3016 root->bCellMin[i] = FALSE;
3022 /* We need the total for normalization */
3024 for(i=range[0]; i<range[1]; i++)
3026 if (root->bCellMin[i] == FALSE)
3028 fac += cell_size[i];
3031 fac = ( region_size - nmin*cellsize_limit_f)/fac; /* substracting cells already set to cellsize_limit_f */
3032 /* Determine the cell boundaries */
3033 for(i=range[0]; i<range[1]; i++)
3035 if (root->bCellMin[i] == FALSE)
3037 cell_size[i] *= fac;
3038 if (!bPBC && (i == 0 || i == dd->nc[dim] -1))
3040 cellsize_limit_f_i = 0;
3044 cellsize_limit_f_i = cellsize_limit_f;
3046 if (cell_size[i] < cellsize_limit_f_i)
3048 root->bCellMin[i] = TRUE;
3049 cell_size[i] = cellsize_limit_f_i;
3053 root->cell_f[i+1] = root->cell_f[i] + cell_size[i];
3056 while (nmin > nmin_old);
3059 cell_size[i] = root->cell_f[i+1] - root->cell_f[i];
3060 /* For this check we should not use DD_CELL_MARGIN,
3061 * but a slightly smaller factor,
3062 * since rounding could get use below the limit.
3064 if (bPBC && cell_size[i] < cellsize_limit_f*DD_CELL_MARGIN2/DD_CELL_MARGIN)
3067 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",
3068 gmx_step_str(step,buf),
3069 dim2char(dim),ddbox->box_size[dim],ddbox->skew_fac[dim],
3070 ncd,comm->cellsize_min[dim]);
3073 root->bLimited = (nmin > 0) || (range[0]>0) || (range[1]<ncd);
3077 /* Check if the boundary did not displace more than halfway
3078 * each of the cells it bounds, as this could cause problems,
3079 * especially when the differences between cell sizes are large.
3080 * If changes are applied, they will not make cells smaller
3081 * than the cut-off, as we check all the boundaries which
3082 * might be affected by a change and if the old state was ok,
3083 * the cells will at most be shrunk back to their old size.
3085 for(i=range[0]+1; i<range[1]; i++)
3087 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i-1]);
3088 if (root->cell_f[i] < halfway)
3090 root->cell_f[i] = halfway;
3091 /* Check if the change also causes shifts of the next boundaries */
3092 for(j=i+1; j<range[1]; j++)
3094 if (root->cell_f[j] < root->cell_f[j-1] + cellsize_limit_f)
3095 root->cell_f[j] = root->cell_f[j-1] + cellsize_limit_f;
3098 halfway = 0.5*(root->old_cell_f[i] + root->old_cell_f[i+1]);
3099 if (root->cell_f[i] > halfway)
3101 root->cell_f[i] = halfway;
3102 /* Check if the change also causes shifts of the next boundaries */
3103 for(j=i-1; j>=range[0]+1; j--)
3105 if (root->cell_f[j] > root->cell_f[j+1] - cellsize_limit_f)
3106 root->cell_f[j] = root->cell_f[j+1] - cellsize_limit_f;
3112 /* nrange is defined as [lower, upper) range for new call to enforce_limits */
3113 /* find highest violation of LimLo (a) and the following violation of LimHi (thus the lowest following) (b)
3114 * then call enforce_limits for (oldb,a), (a,b). In the next step: (b,nexta). oldb and nexta can be the boundaries.
3115 * for a and b nrange is used */
3118 /* Take care of the staggering of the cell boundaries */
3121 for(i=range[0]; i<range[1]; i++)
3123 root->cell_f_max0[i] = root->cell_f[i];
3124 root->cell_f_min1[i] = root->cell_f[i+1];
3129 for(i=range[0]+1; i<range[1]; i++)
3131 bLimLo = (root->cell_f[i] < root->bound_min[i]);
3132 bLimHi = (root->cell_f[i] > root->bound_max[i]);
3133 if (bLimLo && bLimHi)
3135 /* Both limits violated, try the best we can */
3136 /* For this case we split the original range (range) in two parts and care about the other limitiations in the next iteration. */
3137 root->cell_f[i] = 0.5*(root->bound_min[i] + root->bound_max[i]);
3140 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3144 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3150 /* root->cell_f[i] = root->bound_min[i]; */
3151 nrange[1]=i; /* only store violation location. There could be a LimLo violation following with an higher index */
3154 else if (bLimHi && !bLastHi)
3157 if (nrange[1] < range[1]) /* found a LimLo before */
3159 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3160 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3161 nrange[0]=nrange[1];
3163 root->cell_f[i] = root->bound_max[i];
3165 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3170 if (nrange[1] < range[1]) /* found last a LimLo */
3172 root->cell_f[nrange[1]] = root->bound_min[nrange[1]];
3173 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3174 nrange[0]=nrange[1];
3176 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3178 else if (nrange[0] > range[0]) /* found at least one LimHi */
3180 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, nrange);
3187 static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
3188 int d,int dim,gmx_domdec_root_t *root,
3189 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3190 gmx_bool bUniform,gmx_large_int_t step)
3192 gmx_domdec_comm_t *comm;
3195 real load_aver,load_i,imbalance,change,change_max,sc;
3196 real cellsize_limit_f,dist_min_f,dist_min_f_hard,space;
3200 int range[] = { 0, 0 };
3204 /* Convert the maximum change from the input percentage to a fraction */
3205 change_limit = comm->dlb_scale_lim*0.01;
3209 bPBC = (dim < ddbox->npbcdim);
3211 cell_size = root->buf_ncd;
3213 /* Store the original boundaries */
3214 for(i=0; i<ncd+1; i++)
3216 root->old_cell_f[i] = root->cell_f[i];
3219 for(i=0; i<ncd; i++)
3221 cell_size[i] = 1.0/ncd;
3224 else if (dd_load_count(comm))
3226 load_aver = comm->load[d].sum_m/ncd;
3228 for(i=0; i<ncd; i++)
3230 /* Determine the relative imbalance of cell i */
3231 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3232 imbalance = (load_i - load_aver)/(load_aver>0 ? load_aver : 1);
3233 /* Determine the change of the cell size using underrelaxation */
3234 change = -relax*imbalance;
3235 change_max = max(change_max,max(change,-change));
3237 /* Limit the amount of scaling.
3238 * We need to use the same rescaling for all cells in one row,
3239 * otherwise the load balancing might not converge.
3242 if (change_max > change_limit)
3244 sc *= change_limit/change_max;
3246 for(i=0; i<ncd; i++)
3248 /* Determine the relative imbalance of cell i */
3249 load_i = comm->load[d].load[i*comm->load[d].nload+2];
3250 imbalance = (load_i - load_aver)/(load_aver>0 ? load_aver : 1);
3251 /* Determine the change of the cell size using underrelaxation */
3252 change = -sc*imbalance;
3253 cell_size[i] = (root->cell_f[i+1]-root->cell_f[i])*(1 + change);
3257 cellsize_limit_f = comm->cellsize_min[dim]/ddbox->box_size[dim];
3258 cellsize_limit_f *= DD_CELL_MARGIN;
3259 dist_min_f_hard = grid_jump_limit(comm,d)/ddbox->box_size[dim];
3260 dist_min_f = dist_min_f_hard * DD_CELL_MARGIN;
3261 if (ddbox->tric_dir[dim])
3263 cellsize_limit_f /= ddbox->skew_fac[dim];
3264 dist_min_f /= ddbox->skew_fac[dim];
3266 if (bDynamicBox && d > 0)
3268 dist_min_f *= DD_PRES_SCALE_MARGIN;
3270 if (d > 0 && !bUniform)
3272 /* Make sure that the grid is not shifted too much */
3273 for(i=1; i<ncd; i++) {
3274 if (root->cell_f_min1[i] - root->cell_f_max0[i-1] < 2 * dist_min_f_hard)
3276 gmx_incons("Inconsistent DD boundary staggering limits!");
3278 root->bound_min[i] = root->cell_f_max0[i-1] + dist_min_f;
3279 space = root->cell_f[i] - (root->cell_f_max0[i-1] + dist_min_f);
3281 root->bound_min[i] += 0.5*space;
3283 root->bound_max[i] = root->cell_f_min1[i] - dist_min_f;
3284 space = root->cell_f[i] - (root->cell_f_min1[i] - dist_min_f);
3286 root->bound_max[i] += 0.5*space;
3291 "dim %d boundary %d %.3f < %.3f < %.3f < %.3f < %.3f\n",
3293 root->cell_f_max0[i-1] + dist_min_f,
3294 root->bound_min[i],root->cell_f[i],root->bound_max[i],
3295 root->cell_f_min1[i] - dist_min_f);
3300 root->cell_f[0] = 0;
3301 root->cell_f[ncd] = 1;
3302 dd_cell_sizes_dlb_root_enforce_limits(dd, d, dim, root, ddbox, bUniform, step, cellsize_limit_f, range);
3305 /* After the checks above, the cells should obey the cut-off
3306 * restrictions, but it does not hurt to check.
3308 for(i=0; i<ncd; i++)
3312 fprintf(debug,"Relative bounds dim %d cell %d: %f %f\n",
3313 dim,i,root->cell_f[i],root->cell_f[i+1]);
3316 if ((bPBC || (i != 0 && i != dd->nc[dim]-1)) &&
3317 root->cell_f[i+1] - root->cell_f[i] <
3318 cellsize_limit_f/DD_CELL_MARGIN)
3322 "\nWARNING step %s: direction %c, cell %d too small: %f\n",
3323 gmx_step_str(step,buf),dim2char(dim),i,
3324 (root->cell_f[i+1] - root->cell_f[i])
3325 *ddbox->box_size[dim]*ddbox->skew_fac[dim]);
3330 /* Store the cell boundaries of the lower dimensions at the end */
3331 for(d1=0; d1<d; d1++)
3333 root->cell_f[pos++] = comm->cell_f0[d1];
3334 root->cell_f[pos++] = comm->cell_f1[d1];
3337 if (d < comm->npmedecompdim)
3339 /* The master determines the maximum shift for
3340 * the coordinate communication between separate PME nodes.
3342 set_pme_maxshift(dd,&comm->ddpme[d],bUniform,ddbox,root->cell_f);
3344 root->cell_f[pos++] = comm->ddpme[0].maxshift;
3347 root->cell_f[pos++] = comm->ddpme[1].maxshift;
3351 static void relative_to_absolute_cell_bounds(gmx_domdec_t *dd,
3352 gmx_ddbox_t *ddbox,int dimind)
3354 gmx_domdec_comm_t *comm;
3359 /* Set the cell dimensions */
3360 dim = dd->dim[dimind];
3361 comm->cell_x0[dim] = comm->cell_f0[dimind]*ddbox->box_size[dim];
3362 comm->cell_x1[dim] = comm->cell_f1[dimind]*ddbox->box_size[dim];
3363 if (dim >= ddbox->nboundeddim)
3365 comm->cell_x0[dim] += ddbox->box0[dim];
3366 comm->cell_x1[dim] += ddbox->box0[dim];
3370 static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3371 int d,int dim,real *cell_f_row,
3374 gmx_domdec_comm_t *comm;
3380 /* Each node would only need to know two fractions,
3381 * but it is probably cheaper to broadcast the whole array.
3383 MPI_Bcast(cell_f_row,DD_CELL_F_SIZE(dd,d)*sizeof(real),MPI_BYTE,
3384 0,comm->mpi_comm_load[d]);
3386 /* Copy the fractions for this dimension from the buffer */
3387 comm->cell_f0[d] = cell_f_row[dd->ci[dim] ];
3388 comm->cell_f1[d] = cell_f_row[dd->ci[dim]+1];
3389 /* The whole array was communicated, so set the buffer position */
3390 pos = dd->nc[dim] + 1;
3391 for(d1=0; d1<=d; d1++)
3395 /* Copy the cell fractions of the lower dimensions */
3396 comm->cell_f0[d1] = cell_f_row[pos++];
3397 comm->cell_f1[d1] = cell_f_row[pos++];
3399 relative_to_absolute_cell_bounds(dd,ddbox,d1);
3401 /* Convert the communicated shift from float to int */
3402 comm->ddpme[0].maxshift = (int)(cell_f_row[pos++] + 0.5);
3405 comm->ddpme[1].maxshift = (int)(cell_f_row[pos++] + 0.5);
3409 static void set_dd_cell_sizes_dlb_change(gmx_domdec_t *dd,
3410 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3411 gmx_bool bUniform,gmx_large_int_t step)
3413 gmx_domdec_comm_t *comm;
3415 gmx_bool bRowMember,bRowRoot;
3420 for(d=0; d<dd->ndim; d++)
3425 for(d1=d; d1<dd->ndim; d1++)
3427 if (dd->ci[dd->dim[d1]] > 0)
3440 set_dd_cell_sizes_dlb_root(dd,d,dim,comm->root[d],
3441 ddbox,bDynamicBox,bUniform,step);
3442 cell_f_row = comm->root[d]->cell_f;
3446 cell_f_row = comm->cell_f_row;
3448 distribute_dd_cell_sizes_dlb(dd,d,dim,cell_f_row,ddbox);
3453 static void set_dd_cell_sizes_dlb_nochange(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
3457 /* This function assumes the box is static and should therefore
3458 * not be called when the box has changed since the last
3459 * call to dd_partition_system.
3461 for(d=0; d<dd->ndim; d++)
3463 relative_to_absolute_cell_bounds(dd,ddbox,d);
3469 static void set_dd_cell_sizes_dlb(gmx_domdec_t *dd,
3470 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3471 gmx_bool bUniform,gmx_bool bDoDLB,gmx_large_int_t step,
3472 gmx_wallcycle_t wcycle)
3474 gmx_domdec_comm_t *comm;
3481 wallcycle_start(wcycle,ewcDDCOMMBOUND);
3482 set_dd_cell_sizes_dlb_change(dd,ddbox,bDynamicBox,bUniform,step);
3483 wallcycle_stop(wcycle,ewcDDCOMMBOUND);
3485 else if (bDynamicBox)
3487 set_dd_cell_sizes_dlb_nochange(dd,ddbox);
3490 /* Set the dimensions for which no DD is used */
3491 for(dim=0; dim<DIM; dim++) {
3492 if (dd->nc[dim] == 1) {
3493 comm->cell_x0[dim] = 0;
3494 comm->cell_x1[dim] = ddbox->box_size[dim];
3495 if (dim >= ddbox->nboundeddim)
3497 comm->cell_x0[dim] += ddbox->box0[dim];
3498 comm->cell_x1[dim] += ddbox->box0[dim];
3504 static void realloc_comm_ind(gmx_domdec_t *dd,ivec npulse)
3507 gmx_domdec_comm_dim_t *cd;
3509 for(d=0; d<dd->ndim; d++)
3511 cd = &dd->comm->cd[d];
3512 np = npulse[dd->dim[d]];
3513 if (np > cd->np_nalloc)
3517 fprintf(debug,"(Re)allocing cd for %c to %d pulses\n",
3518 dim2char(dd->dim[d]),np);
3520 if (DDMASTER(dd) && cd->np_nalloc > 0)
3522 fprintf(stderr,"\nIncreasing the number of cell to communicate in dimension %c to %d for the first time\n",dim2char(dd->dim[d]),np);
3525 for(i=cd->np_nalloc; i<np; i++)
3527 cd->ind[i].index = NULL;
3528 cd->ind[i].nalloc = 0;
3537 static void set_dd_cell_sizes(gmx_domdec_t *dd,
3538 gmx_ddbox_t *ddbox,gmx_bool bDynamicBox,
3539 gmx_bool bUniform,gmx_bool bDoDLB,gmx_large_int_t step,
3540 gmx_wallcycle_t wcycle)
3542 gmx_domdec_comm_t *comm;
3548 /* Copy the old cell boundaries for the cg displacement check */
3549 copy_rvec(comm->cell_x0,comm->old_cell_x0);
3550 copy_rvec(comm->cell_x1,comm->old_cell_x1);
3552 if (comm->bDynLoadBal)
3556 check_box_size(dd,ddbox);
3558 set_dd_cell_sizes_dlb(dd,ddbox,bDynamicBox,bUniform,bDoDLB,step,wcycle);
3562 set_dd_cell_sizes_slb(dd,ddbox,FALSE,npulse);
3563 realloc_comm_ind(dd,npulse);
3568 for(d=0; d<DIM; d++)
3570 fprintf(debug,"cell_x[%d] %f - %f skew_fac %f\n",
3571 d,comm->cell_x0[d],comm->cell_x1[d],ddbox->skew_fac[d]);
3576 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
3578 rvec cell_ns_x0,rvec cell_ns_x1,
3579 gmx_large_int_t step)
3581 gmx_domdec_comm_t *comm;
3586 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
3588 dim = dd->dim[dim_ind];
3590 /* Without PBC we don't have restrictions on the outer cells */
3591 if (!(dim >= ddbox->npbcdim &&
3592 (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
3593 comm->bDynLoadBal &&
3594 (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
3595 comm->cellsize_min[dim])
3598 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",
3599 gmx_step_str(step,buf),dim2char(dim),
3600 comm->cell_x1[dim] - comm->cell_x0[dim],
3601 ddbox->skew_fac[dim],
3602 dd->comm->cellsize_min[dim],
3603 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
3607 if ((dd->bGridJump && dd->ndim > 1) || ddbox->nboundeddim < DIM)
3609 /* Communicate the boundaries and update cell_ns_x0/1 */
3610 dd_move_cellx(dd,ddbox,cell_ns_x0,cell_ns_x1);
3611 if (dd->bGridJump && dd->ndim > 1)
3613 check_grid_jump(step,dd,ddbox);
3618 static void make_tric_corr_matrix(int npbcdim,matrix box,matrix tcm)
3622 tcm[YY][XX] = -box[YY][XX]/box[YY][YY];
3630 tcm[ZZ][XX] = -(box[ZZ][YY]*tcm[YY][XX] + box[ZZ][XX])/box[ZZ][ZZ];
3631 tcm[ZZ][YY] = -box[ZZ][YY]/box[ZZ][ZZ];
3640 static void check_screw_box(matrix box)
3642 /* Mathematical limitation */
3643 if (box[YY][XX] != 0 || box[ZZ][XX] != 0)
3645 gmx_fatal(FARGS,"With screw pbc the unit cell can not have non-zero off-diagonal x-components");
3648 /* Limitation due to the asymmetry of the eighth shell method */
3649 if (box[ZZ][YY] != 0)
3651 gmx_fatal(FARGS,"pbc=screw with non-zero box_zy is not supported");
3655 static void distribute_cg(FILE *fplog,gmx_large_int_t step,
3656 matrix box,ivec tric_dir,t_block *cgs,rvec pos[],
3659 gmx_domdec_master_t *ma;
3660 int **tmp_ind=NULL,*tmp_nalloc=NULL;
3661 int i,icg,j,k,k0,k1,d,npbcdim;
3663 rvec box_size,cg_cm;
3665 real nrcg,inv_ncg,pos_d;
3667 gmx_bool bUnbounded,bScrew;
3671 if (tmp_ind == NULL)
3673 snew(tmp_nalloc,dd->nnodes);
3674 snew(tmp_ind,dd->nnodes);
3675 for(i=0; i<dd->nnodes; i++)
3677 tmp_nalloc[i] = over_alloc_large(cgs->nr/dd->nnodes+1);
3678 snew(tmp_ind[i],tmp_nalloc[i]);
3682 /* Clear the count */
3683 for(i=0; i<dd->nnodes; i++)
3689 make_tric_corr_matrix(dd->npbcdim,box,tcm);
3691 cgindex = cgs->index;
3693 /* Compute the center of geometry for all charge groups */
3694 for(icg=0; icg<cgs->nr; icg++)
3697 k1 = cgindex[icg+1];
3701 copy_rvec(pos[k0],cg_cm);
3708 for(k=k0; (k<k1); k++)
3710 rvec_inc(cg_cm,pos[k]);
3712 for(d=0; (d<DIM); d++)
3714 cg_cm[d] *= inv_ncg;
3717 /* Put the charge group in the box and determine the cell index */
3718 for(d=DIM-1; d>=0; d--) {
3720 if (d < dd->npbcdim)
3722 bScrew = (dd->bScrewPBC && d == XX);
3723 if (tric_dir[d] && dd->nc[d] > 1)
3725 /* Use triclinic coordintates for this dimension */
3726 for(j=d+1; j<DIM; j++)
3728 pos_d += cg_cm[j]*tcm[j][d];
3731 while(pos_d >= box[d][d])
3734 rvec_dec(cg_cm,box[d]);
3737 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3738 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3740 for(k=k0; (k<k1); k++)
3742 rvec_dec(pos[k],box[d]);
3745 pos[k][YY] = box[YY][YY] - pos[k][YY];
3746 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3753 rvec_inc(cg_cm,box[d]);
3756 cg_cm[YY] = box[YY][YY] - cg_cm[YY];
3757 cg_cm[ZZ] = box[ZZ][ZZ] - cg_cm[ZZ];
3759 for(k=k0; (k<k1); k++)
3761 rvec_inc(pos[k],box[d]);
3763 pos[k][YY] = box[YY][YY] - pos[k][YY];
3764 pos[k][ZZ] = box[ZZ][ZZ] - pos[k][ZZ];
3769 /* This could be done more efficiently */
3771 while(ind[d]+1 < dd->nc[d] && pos_d >= ma->cell_x[d][ind[d]+1])
3776 i = dd_index(dd->nc,ind);
3777 if (ma->ncg[i] == tmp_nalloc[i])
3779 tmp_nalloc[i] = over_alloc_large(ma->ncg[i]+1);
3780 srenew(tmp_ind[i],tmp_nalloc[i]);
3782 tmp_ind[i][ma->ncg[i]] = icg;
3784 ma->nat[i] += cgindex[icg+1] - cgindex[icg];
3788 for(i=0; i<dd->nnodes; i++)
3791 for(k=0; k<ma->ncg[i]; k++)
3793 ma->cg[k1++] = tmp_ind[i][k];
3796 ma->index[dd->nnodes] = k1;
3798 for(i=0; i<dd->nnodes; i++)
3808 fprintf(fplog,"Charge group distribution at step %s:",
3809 gmx_step_str(step,buf));
3810 for(i=0; i<dd->nnodes; i++)
3812 fprintf(fplog," %d",ma->ncg[i]);
3814 fprintf(fplog,"\n");
3818 static void get_cg_distribution(FILE *fplog,gmx_large_int_t step,gmx_domdec_t *dd,
3819 t_block *cgs,matrix box,gmx_ddbox_t *ddbox,
3822 gmx_domdec_master_t *ma=NULL;
3825 int *ibuf,buf2[2] = { 0, 0 };
3826 gmx_bool bMaster = DDMASTER(dd);
3833 check_screw_box(box);
3836 set_dd_cell_sizes_slb(dd,ddbox,TRUE,npulse);
3838 distribute_cg(fplog,step,box,ddbox->tric_dir,cgs,pos,dd);
3839 for(i=0; i<dd->nnodes; i++)
3841 ma->ibuf[2*i] = ma->ncg[i];
3842 ma->ibuf[2*i+1] = ma->nat[i];
3850 dd_scatter(dd,2*sizeof(int),ibuf,buf2);
3852 dd->ncg_home = buf2[0];
3853 dd->nat_home = buf2[1];
3854 dd->ncg_tot = dd->ncg_home;
3855 dd->nat_tot = dd->nat_home;
3856 if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
3858 dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
3859 srenew(dd->index_gl,dd->cg_nalloc);
3860 srenew(dd->cgindex,dd->cg_nalloc+1);
3864 for(i=0; i<dd->nnodes; i++)
3866 ma->ibuf[i] = ma->ncg[i]*sizeof(int);
3867 ma->ibuf[dd->nnodes+i] = ma->index[i]*sizeof(int);
3872 DDMASTER(dd) ? ma->ibuf : NULL,
3873 DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
3874 DDMASTER(dd) ? ma->cg : NULL,
3875 dd->ncg_home*sizeof(int),dd->index_gl);
3877 /* Determine the home charge group sizes */
3879 for(i=0; i<dd->ncg_home; i++)
3881 cg_gl = dd->index_gl[i];
3883 dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
3888 fprintf(debug,"Home charge groups:\n");
3889 for(i=0; i<dd->ncg_home; i++)
3891 fprintf(debug," %d",dd->index_gl[i]);
3893 fprintf(debug,"\n");
3895 fprintf(debug,"\n");
3899 static int compact_and_copy_vec_at(int ncg,int *move,
3902 rvec *src,gmx_domdec_comm_t *comm,
3905 int m,icg,i,i0,i1,nrcg;
3911 for(m=0; m<DIM*2; m++)
3917 for(icg=0; icg<ncg; icg++)
3919 i1 = cgindex[icg+1];
3925 /* Compact the home array in place */
3926 for(i=i0; i<i1; i++)
3928 copy_rvec(src[i],src[home_pos++]);
3934 /* Copy to the communication buffer */
3936 pos_vec[m] += 1 + vec*nrcg;
3937 for(i=i0; i<i1; i++)
3939 copy_rvec(src[i],comm->cgcm_state[m][pos_vec[m]++]);
3941 pos_vec[m] += (nvec - vec - 1)*nrcg;
3945 home_pos += i1 - i0;
3953 static int compact_and_copy_vec_cg(int ncg,int *move,
3955 int nvec,rvec *src,gmx_domdec_comm_t *comm,
3958 int m,icg,i0,i1,nrcg;
3964 for(m=0; m<DIM*2; m++)
3970 for(icg=0; icg<ncg; icg++)
3972 i1 = cgindex[icg+1];
3978 /* Compact the home array in place */
3979 copy_rvec(src[icg],src[home_pos++]);
3985 /* Copy to the communication buffer */
3986 copy_rvec(src[icg],comm->cgcm_state[m][pos_vec[m]]);
3987 pos_vec[m] += 1 + nrcg*nvec;
3999 static int compact_ind(int ncg,int *move,
4000 int *index_gl,int *cgindex,
4002 gmx_ga2la_t ga2la,char *bLocalCG,
4005 int cg,nat,a0,a1,a,a_gl;
4010 for(cg=0; cg<ncg; cg++)
4016 /* Compact the home arrays in place.
4017 * Anything that can be done here avoids access to global arrays.
4019 cgindex[home_pos] = nat;
4020 for(a=a0; a<a1; a++)
4023 gatindex[nat] = a_gl;
4024 /* The cell number stays 0, so we don't need to set it */
4025 ga2la_change_la(ga2la,a_gl,nat);
4028 index_gl[home_pos] = index_gl[cg];
4029 cginfo[home_pos] = cginfo[cg];
4030 /* The charge group remains local, so bLocalCG does not change */
4035 /* Clear the global indices */
4036 for(a=a0; a<a1; a++)
4038 ga2la_del(ga2la,gatindex[a]);
4042 bLocalCG[index_gl[cg]] = FALSE;
4046 cgindex[home_pos] = nat;
4051 static void clear_and_mark_ind(int ncg,int *move,
4052 int *index_gl,int *cgindex,int *gatindex,
4053 gmx_ga2la_t ga2la,char *bLocalCG,
4058 for(cg=0; cg<ncg; cg++)
4064 /* Clear the global indices */
4065 for(a=a0; a<a1; a++)
4067 ga2la_del(ga2la,gatindex[a]);
4071 bLocalCG[index_gl[cg]] = FALSE;
4073 /* Signal that this cg has moved using the ns cell index.
4074 * Here we set it to -1.
4075 * fill_grid will change it from -1 to 4*grid->ncells.
4077 cell_index[cg] = -1;
4082 static void print_cg_move(FILE *fplog,
4084 gmx_large_int_t step,int cg,int dim,int dir,
4085 gmx_bool bHaveLimitdAndCMOld,real limitd,
4086 rvec cm_old,rvec cm_new,real pos_d)
4088 gmx_domdec_comm_t *comm;
4093 fprintf(fplog,"\nStep %s:\n",gmx_step_str(step,buf));
4094 if (bHaveLimitdAndCMOld)
4096 fprintf(fplog,"The charge group starting at atom %d moved than the distance allowed by the domain decomposition (%f) in direction %c\n",
4097 ddglatnr(dd,dd->cgindex[cg]),limitd,dim2char(dim));
4101 fprintf(fplog,"The charge group starting at atom %d moved than the distance allowed by the domain decomposition in direction %c\n",
4102 ddglatnr(dd,dd->cgindex[cg]),dim2char(dim));
4104 fprintf(fplog,"distance out of cell %f\n",
4105 dir==1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
4106 if (bHaveLimitdAndCMOld)
4108 fprintf(fplog,"Old coordinates: %8.3f %8.3f %8.3f\n",
4109 cm_old[XX],cm_old[YY],cm_old[ZZ]);
4111 fprintf(fplog,"New coordinates: %8.3f %8.3f %8.3f\n",
4112 cm_new[XX],cm_new[YY],cm_new[ZZ]);
4113 fprintf(fplog,"Old cell boundaries in direction %c: %8.3f %8.3f\n",
4115 comm->old_cell_x0[dim],comm->old_cell_x1[dim]);
4116 fprintf(fplog,"New cell boundaries in direction %c: %8.3f %8.3f\n",
4118 comm->cell_x0[dim],comm->cell_x1[dim]);
4121 static void cg_move_error(FILE *fplog,
4123 gmx_large_int_t step,int cg,int dim,int dir,
4124 gmx_bool bHaveLimitdAndCMOld,real limitd,
4125 rvec cm_old,rvec cm_new,real pos_d)
4129 print_cg_move(fplog, dd,step,cg,dim,dir,
4130 bHaveLimitdAndCMOld,limitd,cm_old,cm_new,pos_d);
4132 print_cg_move(stderr,dd,step,cg,dim,dir,
4133 bHaveLimitdAndCMOld,limitd,cm_old,cm_new,pos_d);
4135 "A charge group moved too far between two domain decomposition steps\n"
4136 "This usually means that your system is not well equilibrated");
4139 static void rotate_state_atom(t_state *state,int a)
4143 for(est=0; est<estNR; est++)
4145 if (EST_DISTR(est) && (state->flags & (1<<est))) {
4148 /* Rotate the complete state; for a rectangular box only */
4149 state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
4150 state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
4153 state->v[a][YY] = -state->v[a][YY];
4154 state->v[a][ZZ] = -state->v[a][ZZ];
4157 state->sd_X[a][YY] = -state->sd_X[a][YY];
4158 state->sd_X[a][ZZ] = -state->sd_X[a][ZZ];
4161 state->cg_p[a][YY] = -state->cg_p[a][YY];
4162 state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
4164 case estDISRE_INITF:
4165 case estDISRE_RM3TAV:
4166 case estORIRE_INITF:
4168 /* These are distances, so not affected by rotation */
4171 gmx_incons("Unknown state entry encountered in rotate_state_atom");
4177 static int dd_redistribute_cg(FILE *fplog,gmx_large_int_t step,
4178 gmx_domdec_t *dd,ivec tric_dir,
4179 t_state *state,rvec **f,
4180 t_forcerec *fr,t_mdatoms *md,
4186 int ncg[DIM*2],nat[DIM*2];
4187 int c,i,cg,k,k0,k1,d,dim,dim2,dir,d2,d3,d4,cell_d;
4188 int mc,cdd,nrcg,ncg_recv,nat_recv,nvs,nvr,nvec,vec;
4189 int sbuf[2],rbuf[2];
4190 int home_pos_cg,home_pos_at,ncg_stay_home,buf_pos;
4192 gmx_bool bV=FALSE,bSDX=FALSE,bCGP=FALSE;
4197 rvec *cg_cm,cell_x0,cell_x1,limitd,limit0,limit1,cm_new;
4199 cginfo_mb_t *cginfo_mb;
4200 gmx_domdec_comm_t *comm;
4204 check_screw_box(state->box);
4210 for(i=0; i<estNR; i++)
4216 case estX: /* Always present */ break;
4217 case estV: bV = (state->flags & (1<<i)); break;
4218 case estSDX: bSDX = (state->flags & (1<<i)); break;
4219 case estCGP: bCGP = (state->flags & (1<<i)); break;
4222 case estDISRE_INITF:
4223 case estDISRE_RM3TAV:
4224 case estORIRE_INITF:
4226 /* No processing required */
4229 gmx_incons("Unknown state entry encountered in dd_redistribute_cg");
4234 if (dd->ncg_tot > comm->nalloc_int)
4236 comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
4237 srenew(comm->buf_int,comm->nalloc_int);
4239 move = comm->buf_int;
4241 /* Clear the count */
4242 for(c=0; c<dd->ndim*2; c++)
4248 npbcdim = dd->npbcdim;
4250 for(d=0; (d<DIM); d++)
4252 limitd[d] = dd->comm->cellsize_min[d];
4253 if (d >= npbcdim && dd->ci[d] == 0)
4255 cell_x0[d] = -GMX_FLOAT_MAX;
4259 cell_x0[d] = comm->cell_x0[d];
4261 if (d >= npbcdim && dd->ci[d] == dd->nc[d] - 1)
4263 cell_x1[d] = GMX_FLOAT_MAX;
4267 cell_x1[d] = comm->cell_x1[d];
4271 limit0[d] = comm->old_cell_x0[d] - limitd[d];
4272 limit1[d] = comm->old_cell_x1[d] + limitd[d];
4276 /* We check after communication if a charge group moved
4277 * more than one cell. Set the pre-comm check limit to float_max.
4279 limit0[d] = -GMX_FLOAT_MAX;
4280 limit1[d] = GMX_FLOAT_MAX;
4284 make_tric_corr_matrix(npbcdim,state->box,tcm);
4286 cgindex = dd->cgindex;
4288 /* Compute the center of geometry for all home charge groups
4289 * and put them in the box and determine where they should go.
4291 for(cg=0; cg<dd->ncg_home; cg++)
4298 copy_rvec(state->x[k0],cm_new);
4305 for(k=k0; (k<k1); k++)
4307 rvec_inc(cm_new,state->x[k]);
4309 for(d=0; (d<DIM); d++)
4311 cm_new[d] = inv_ncg*cm_new[d];
4316 /* Do pbc and check DD cell boundary crossings */
4317 for(d=DIM-1; d>=0; d--)
4321 bScrew = (dd->bScrewPBC && d == XX);
4322 /* Determine the location of this cg in lattice coordinates */
4326 for(d2=d+1; d2<DIM; d2++)
4328 pos_d += cm_new[d2]*tcm[d2][d];
4331 /* Put the charge group in the triclinic unit-cell */
4332 if (pos_d >= cell_x1[d])
4334 if (pos_d >= limit1[d])
4336 cg_move_error(fplog,dd,step,cg,d,1,TRUE,limitd[d],
4337 cg_cm[cg],cm_new,pos_d);
4340 if (dd->ci[d] == dd->nc[d] - 1)
4342 rvec_dec(cm_new,state->box[d]);
4345 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4346 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4348 for(k=k0; (k<k1); k++)
4350 rvec_dec(state->x[k],state->box[d]);
4353 rotate_state_atom(state,k);
4358 else if (pos_d < cell_x0[d])
4360 if (pos_d < limit0[d])
4362 cg_move_error(fplog,dd,step,cg,d,-1,TRUE,limitd[d],
4363 cg_cm[cg],cm_new,pos_d);
4368 rvec_inc(cm_new,state->box[d]);
4371 cm_new[YY] = state->box[YY][YY] - cm_new[YY];
4372 cm_new[ZZ] = state->box[ZZ][ZZ] - cm_new[ZZ];
4374 for(k=k0; (k<k1); k++)
4376 rvec_inc(state->x[k],state->box[d]);
4379 rotate_state_atom(state,k);
4385 else if (d < npbcdim)
4387 /* Put the charge group in the rectangular unit-cell */
4388 while (cm_new[d] >= state->box[d][d])
4390 rvec_dec(cm_new,state->box[d]);
4391 for(k=k0; (k<k1); k++)
4393 rvec_dec(state->x[k],state->box[d]);
4396 while (cm_new[d] < 0)
4398 rvec_inc(cm_new,state->box[d]);
4399 for(k=k0; (k<k1); k++)
4401 rvec_inc(state->x[k],state->box[d]);
4407 copy_rvec(cm_new,cg_cm[cg]);
4409 /* Determine where this cg should go */
4412 for(d=0; d<dd->ndim; d++)
4417 flag |= DD_FLAG_FW(d);
4423 else if (dev[dim] == -1)
4425 flag |= DD_FLAG_BW(d);
4427 if (dd->nc[dim] > 2)
4441 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4443 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4444 srenew(comm->cggl_flag[mc],comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4446 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS ] = dd->index_gl[cg];
4447 /* We store the cg size in the lower 16 bits
4448 * and the place where the charge group should go
4449 * in the next 6 bits. This saves some communication volume.
4451 comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
4457 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
4458 inc_nrnb(nrnb,eNR_RESETX,dd->ncg_home);
4474 /* Make sure the communication buffers are large enough */
4475 for(mc=0; mc<dd->ndim*2; mc++)
4477 nvr = ncg[mc] + nat[mc]*nvec;
4478 if (nvr > comm->cgcm_state_nalloc[mc])
4480 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr);
4481 srenew(comm->cgcm_state[mc],comm->cgcm_state_nalloc[mc]);
4485 /* Recalculating cg_cm might be cheaper than communicating,
4486 * but that could give rise to rounding issues.
4489 compact_and_copy_vec_cg(dd->ncg_home,move,cgindex,
4490 nvec,cg_cm,comm,bCompact);
4494 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4495 nvec,vec++,state->x,comm,bCompact);
4498 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4499 nvec,vec++,state->v,comm,bCompact);
4503 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4504 nvec,vec++,state->sd_X,comm,bCompact);
4508 compact_and_copy_vec_at(dd->ncg_home,move,cgindex,
4509 nvec,vec++,state->cg_p,comm,bCompact);
4514 compact_ind(dd->ncg_home,move,
4515 dd->index_gl,dd->cgindex,dd->gatindex,
4516 dd->ga2la,comm->bLocalCG,
4521 clear_and_mark_ind(dd->ncg_home,move,
4522 dd->index_gl,dd->cgindex,dd->gatindex,
4523 dd->ga2la,comm->bLocalCG,
4524 fr->ns.grid->cell_index);
4527 cginfo_mb = fr->cginfo_mb;
4529 ncg_stay_home = home_pos_cg;
4530 for(d=0; d<dd->ndim; d++)
4536 for(dir=0; dir<(dd->nc[dim]==2 ? 1 : 2); dir++)
4539 /* Communicate the cg and atom counts */
4544 fprintf(debug,"Sending ddim %d dir %d: ncg %d nat %d\n",
4545 d,dir,sbuf[0],sbuf[1]);
4547 dd_sendrecv_int(dd, d, dir, sbuf, 2, rbuf, 2);
4549 if ((ncg_recv+rbuf[0])*DD_CGIBS > comm->nalloc_int)
4551 comm->nalloc_int = over_alloc_dd((ncg_recv+rbuf[0])*DD_CGIBS);
4552 srenew(comm->buf_int,comm->nalloc_int);
4555 /* Communicate the charge group indices, sizes and flags */
4556 dd_sendrecv_int(dd, d, dir,
4557 comm->cggl_flag[cdd], sbuf[0]*DD_CGIBS,
4558 comm->buf_int+ncg_recv*DD_CGIBS, rbuf[0]*DD_CGIBS);
4560 nvs = ncg[cdd] + nat[cdd]*nvec;
4561 i = rbuf[0] + rbuf[1] *nvec;
4562 vec_rvec_check_alloc(&comm->vbuf,nvr+i);
4564 /* Communicate cgcm and state */
4565 dd_sendrecv_rvec(dd, d, dir,
4566 comm->cgcm_state[cdd], nvs,
4567 comm->vbuf.v+nvr, i);
4568 ncg_recv += rbuf[0];
4569 nat_recv += rbuf[1];
4573 /* Process the received charge groups */
4575 for(cg=0; cg<ncg_recv; cg++)
4577 flag = comm->buf_int[cg*DD_CGIBS+1];
4579 if (dim >= npbcdim && dd->nc[dim] > 2)
4581 /* No pbc in this dim and more than one domain boundary.
4582 * We to a separate check if a charge did not move too far.
4584 if (((flag & DD_FLAG_FW(d)) &&
4585 comm->vbuf.v[buf_pos][d] > cell_x1[dim]) ||
4586 ((flag & DD_FLAG_BW(d)) &&
4587 comm->vbuf.v[buf_pos][d] < cell_x0[dim]))
4589 cg_move_error(fplog,dd,step,cg,d,
4590 (flag & DD_FLAG_FW(d)) ? 1 : 0,
4592 comm->vbuf.v[buf_pos],
4593 comm->vbuf.v[buf_pos],
4594 comm->vbuf.v[buf_pos][d]);
4601 /* Check which direction this cg should go */
4602 for(d2=d+1; (d2<dd->ndim && mc==-1); d2++)
4606 /* The cell boundaries for dimension d2 are not equal
4607 * for each cell row of the lower dimension(s),
4608 * therefore we might need to redetermine where
4609 * this cg should go.
4612 /* If this cg crosses the box boundary in dimension d2
4613 * we can use the communicated flag, so we do not
4614 * have to worry about pbc.
4616 if (!((dd->ci[dim2] == dd->nc[dim2]-1 &&
4617 (flag & DD_FLAG_FW(d2))) ||
4618 (dd->ci[dim2] == 0 &&
4619 (flag & DD_FLAG_BW(d2)))))
4621 /* Clear the two flags for this dimension */
4622 flag &= ~(DD_FLAG_FW(d2) | DD_FLAG_BW(d2));
4623 /* Determine the location of this cg
4624 * in lattice coordinates
4626 pos_d = comm->vbuf.v[buf_pos][dim2];
4629 for(d3=dim2+1; d3<DIM; d3++)
4632 comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
4635 /* Check of we are not at the box edge.
4636 * pbc is only handled in the first step above,
4637 * but this check could move over pbc while
4638 * the first step did not due to different rounding.
4640 if (pos_d >= cell_x1[dim2] &&
4641 dd->ci[dim2] != dd->nc[dim2]-1)
4643 flag |= DD_FLAG_FW(d2);
4645 else if (pos_d < cell_x0[dim2] &&
4648 flag |= DD_FLAG_BW(d2);
4650 comm->buf_int[cg*DD_CGIBS+1] = flag;
4653 /* Set to which neighboring cell this cg should go */
4654 if (flag & DD_FLAG_FW(d2))
4658 else if (flag & DD_FLAG_BW(d2))
4660 if (dd->nc[dd->dim[d2]] > 2)
4672 nrcg = flag & DD_FLAG_NRCG;
4675 if (home_pos_cg+1 > dd->cg_nalloc)
4677 dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
4678 srenew(dd->index_gl,dd->cg_nalloc);
4679 srenew(dd->cgindex,dd->cg_nalloc+1);
4681 /* Set the global charge group index and size */
4682 dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
4683 dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
4684 /* Copy the state from the buffer */
4685 if (home_pos_cg >= fr->cg_nalloc)
4687 dd_realloc_fr_cg(fr,home_pos_cg+1);
4690 copy_rvec(comm->vbuf.v[buf_pos++],cg_cm[home_pos_cg]);
4691 /* Set the cginfo */
4692 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
4693 dd->index_gl[home_pos_cg]);
4696 comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
4699 if (home_pos_at+nrcg > state->nalloc)
4701 dd_realloc_state(state,f,home_pos_at+nrcg);
4703 for(i=0; i<nrcg; i++)
4705 copy_rvec(comm->vbuf.v[buf_pos++],
4706 state->x[home_pos_at+i]);
4710 for(i=0; i<nrcg; i++)
4712 copy_rvec(comm->vbuf.v[buf_pos++],
4713 state->v[home_pos_at+i]);
4718 for(i=0; i<nrcg; i++)
4720 copy_rvec(comm->vbuf.v[buf_pos++],
4721 state->sd_X[home_pos_at+i]);
4726 for(i=0; i<nrcg; i++)
4728 copy_rvec(comm->vbuf.v[buf_pos++],
4729 state->cg_p[home_pos_at+i]);
4733 home_pos_at += nrcg;
4737 /* Reallocate the buffers if necessary */
4738 if (ncg[mc]+1 > comm->cggl_flag_nalloc[mc])
4740 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
4741 srenew(comm->cggl_flag[mc],comm->cggl_flag_nalloc[mc]*DD_CGIBS);
4743 nvr = ncg[mc] + nat[mc]*nvec;
4744 if (nvr + 1 + nrcg*nvec > comm->cgcm_state_nalloc[mc])
4746 comm->cgcm_state_nalloc[mc] = over_alloc_dd(nvr + 1 + nrcg*nvec);
4747 srenew(comm->cgcm_state[mc],comm->cgcm_state_nalloc[mc]);
4749 /* Copy from the receive to the send buffers */
4750 memcpy(comm->cggl_flag[mc] + ncg[mc]*DD_CGIBS,
4751 comm->buf_int + cg*DD_CGIBS,
4752 DD_CGIBS*sizeof(int));
4753 memcpy(comm->cgcm_state[mc][nvr],
4754 comm->vbuf.v[buf_pos],
4755 (1+nrcg*nvec)*sizeof(rvec));
4756 buf_pos += 1 + nrcg*nvec;
4763 /* With sorting (!bCompact) the indices are now only partially up to date
4764 * and ncg_home and nat_home are not the real count, since there are
4765 * "holes" in the arrays for the charge groups that moved to neighbors.
4767 dd->ncg_home = home_pos_cg;
4768 dd->nat_home = home_pos_at;
4772 fprintf(debug,"Finished repartitioning\n");
4775 return ncg_stay_home;
4778 void dd_cycles_add(gmx_domdec_t *dd,float cycles,int ddCycl)
4780 dd->comm->cycl[ddCycl] += cycles;
4781 dd->comm->cycl_n[ddCycl]++;
4782 if (cycles > dd->comm->cycl_max[ddCycl])
4784 dd->comm->cycl_max[ddCycl] = cycles;
4788 static double force_flop_count(t_nrnb *nrnb)
4795 for(i=eNR_NBKERNEL010; i<eNR_NBKERNEL_FREE_ENERGY; i++)
4797 /* To get closer to the real timings, we half the count
4798 * for the normal loops and again half it for water loops.
4801 if (strstr(name,"W3") != NULL || strstr(name,"W4") != NULL)
4803 sum += nrnb->n[i]*0.25*cost_nrnb(i);
4807 sum += nrnb->n[i]*0.50*cost_nrnb(i);
4810 for(i=eNR_NBKERNEL_FREE_ENERGY; i<=eNR_NB14; i++)
4813 if (strstr(name,"W3") != NULL || strstr(name,"W4") != NULL)
4814 sum += nrnb->n[i]*cost_nrnb(i);
4816 for(i=eNR_BONDS; i<=eNR_WALLS; i++)
4818 sum += nrnb->n[i]*cost_nrnb(i);
4824 void dd_force_flop_start(gmx_domdec_t *dd,t_nrnb *nrnb)
4826 if (dd->comm->eFlop)
4828 dd->comm->flop -= force_flop_count(nrnb);
4831 void dd_force_flop_stop(gmx_domdec_t *dd,t_nrnb *nrnb)
4833 if (dd->comm->eFlop)
4835 dd->comm->flop += force_flop_count(nrnb);
4840 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
4844 for(i=0; i<ddCyclNr; i++)
4846 dd->comm->cycl[i] = 0;
4847 dd->comm->cycl_n[i] = 0;
4848 dd->comm->cycl_max[i] = 0;
4851 dd->comm->flop_n = 0;
4854 static void get_load_distribution(gmx_domdec_t *dd,gmx_wallcycle_t wcycle)
4856 gmx_domdec_comm_t *comm;
4857 gmx_domdec_load_t *load;
4858 gmx_domdec_root_t *root=NULL;
4859 int d,dim,cid,i,pos;
4860 float cell_frac=0,sbuf[DD_NLOAD_MAX];
4865 fprintf(debug,"get_load_distribution start\n");
4868 wallcycle_start(wcycle,ewcDDCOMMLOAD);
4872 bSepPME = (dd->pme_nodeid >= 0);
4874 for(d=dd->ndim-1; d>=0; d--)
4877 /* Check if we participate in the communication in this dimension */
4878 if (d == dd->ndim-1 ||
4879 (dd->ci[dd->dim[d+1]]==0 && dd->ci[dd->dim[dd->ndim-1]]==0))
4881 load = &comm->load[d];
4884 cell_frac = comm->cell_f1[d] - comm->cell_f0[d];
4887 if (d == dd->ndim-1)
4889 sbuf[pos++] = dd_force_load(comm);
4890 sbuf[pos++] = sbuf[0];
4893 sbuf[pos++] = sbuf[0];
4894 sbuf[pos++] = cell_frac;
4897 sbuf[pos++] = comm->cell_f_max0[d];
4898 sbuf[pos++] = comm->cell_f_min1[d];
4903 sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
4904 sbuf[pos++] = comm->cycl[ddCyclPME];
4909 sbuf[pos++] = comm->load[d+1].sum;
4910 sbuf[pos++] = comm->load[d+1].max;
4913 sbuf[pos++] = comm->load[d+1].sum_m;
4914 sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
4915 sbuf[pos++] = comm->load[d+1].flags;
4918 sbuf[pos++] = comm->cell_f_max0[d];
4919 sbuf[pos++] = comm->cell_f_min1[d];
4924 sbuf[pos++] = comm->load[d+1].mdf;
4925 sbuf[pos++] = comm->load[d+1].pme;
4929 /* Communicate a row in DD direction d.
4930 * The communicators are setup such that the root always has rank 0.
4933 MPI_Gather(sbuf ,load->nload*sizeof(float),MPI_BYTE,
4934 load->load,load->nload*sizeof(float),MPI_BYTE,
4935 0,comm->mpi_comm_load[d]);
4937 if (dd->ci[dim] == dd->master_ci[dim])
4939 /* We are the root, process this row */
4940 if (comm->bDynLoadBal)
4942 root = comm->root[d];
4952 for(i=0; i<dd->nc[dim]; i++)
4954 load->sum += load->load[pos++];
4955 load->max = max(load->max,load->load[pos]);
4961 /* This direction could not be load balanced properly,
4962 * therefore we need to use the maximum iso the average load.
4964 load->sum_m = max(load->sum_m,load->load[pos]);
4968 load->sum_m += load->load[pos];
4971 load->cvol_min = min(load->cvol_min,load->load[pos]);
4975 load->flags = (int)(load->load[pos++] + 0.5);
4979 root->cell_f_max0[i] = load->load[pos++];
4980 root->cell_f_min1[i] = load->load[pos++];
4985 load->mdf = max(load->mdf,load->load[pos]);
4987 load->pme = max(load->pme,load->load[pos]);
4991 if (comm->bDynLoadBal && root->bLimited)
4993 load->sum_m *= dd->nc[dim];
4994 load->flags |= (1<<d);
5002 comm->nload += dd_load_count(comm);
5003 comm->load_step += comm->cycl[ddCyclStep];
5004 comm->load_sum += comm->load[0].sum;
5005 comm->load_max += comm->load[0].max;
5006 if (comm->bDynLoadBal)
5008 for(d=0; d<dd->ndim; d++)
5010 if (comm->load[0].flags & (1<<d))
5012 comm->load_lim[d]++;
5018 comm->load_mdf += comm->load[0].mdf;
5019 comm->load_pme += comm->load[0].pme;
5023 wallcycle_stop(wcycle,ewcDDCOMMLOAD);
5027 fprintf(debug,"get_load_distribution finished\n");
5031 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
5033 /* Return the relative performance loss on the total run time
5034 * due to the force calculation load imbalance.
5036 if (dd->comm->nload > 0)
5039 (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
5040 (dd->comm->load_step*dd->nnodes);
5048 static void print_dd_load_av(FILE *fplog,gmx_domdec_t *dd)
5051 int npp,npme,nnodes,d,limp;
5052 float imbal,pme_f_ratio,lossf,lossp=0;
5054 gmx_domdec_comm_t *comm;
5057 if (DDMASTER(dd) && comm->nload > 0)
5060 npme = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
5061 nnodes = npp + npme;
5062 imbal = comm->load_max*npp/comm->load_sum - 1;
5063 lossf = dd_force_imb_perf_loss(dd);
5064 sprintf(buf," Average load imbalance: %.1f %%\n",imbal*100);
5065 fprintf(fplog,"%s",buf);
5066 fprintf(stderr,"\n");
5067 fprintf(stderr,"%s",buf);
5068 sprintf(buf," Part of the total run time spent waiting due to load imbalance: %.1f %%\n",lossf*100);
5069 fprintf(fplog,"%s",buf);
5070 fprintf(stderr,"%s",buf);
5072 if (comm->bDynLoadBal)
5074 sprintf(buf," Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
5075 for(d=0; d<dd->ndim; d++)
5077 limp = (200*comm->load_lim[d]+1)/(2*comm->nload);
5078 sprintf(buf+strlen(buf)," %c %d %%",dim2char(dd->dim[d]),limp);
5084 sprintf(buf+strlen(buf),"\n");
5085 fprintf(fplog,"%s",buf);
5086 fprintf(stderr,"%s",buf);
5090 pme_f_ratio = comm->load_pme/comm->load_mdf;
5091 lossp = (comm->load_pme -comm->load_mdf)/comm->load_step;
5094 lossp *= (float)npme/(float)nnodes;
5098 lossp *= (float)npp/(float)nnodes;
5100 sprintf(buf," Average PME mesh/force load: %5.3f\n",pme_f_ratio);
5101 fprintf(fplog,"%s",buf);
5102 fprintf(stderr,"%s",buf);
5103 sprintf(buf," Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n",fabs(lossp)*100);
5104 fprintf(fplog,"%s",buf);
5105 fprintf(stderr,"%s",buf);
5107 fprintf(fplog,"\n");
5108 fprintf(stderr,"\n");
5110 if (lossf >= DD_PERF_LOSS)
5113 "NOTE: %.1f %% performance was lost due to load imbalance\n"
5114 " in the domain decomposition.\n",lossf*100);
5115 if (!comm->bDynLoadBal)
5117 sprintf(buf+strlen(buf)," You might want to use dynamic load balancing (option -dlb.)\n");
5121 sprintf(buf+strlen(buf)," You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
5123 fprintf(fplog,"%s\n",buf);
5124 fprintf(stderr,"%s\n",buf);
5126 if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
5129 "NOTE: %.1f %% performance was lost because the PME nodes\n"
5130 " had %s work to do than the PP nodes.\n"
5131 " You might want to %s the number of PME nodes\n"
5132 " or %s the cut-off and the grid spacing.\n",
5134 (lossp < 0) ? "less" : "more",
5135 (lossp < 0) ? "decrease" : "increase",
5136 (lossp < 0) ? "decrease" : "increase");
5137 fprintf(fplog,"%s\n",buf);
5138 fprintf(stderr,"%s\n",buf);
5143 static float dd_vol_min(gmx_domdec_t *dd)
5145 return dd->comm->load[0].cvol_min*dd->nnodes;
5148 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
5150 return dd->comm->load[0].flags;
5153 static float dd_f_imbal(gmx_domdec_t *dd)
5155 return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1;
5158 static float dd_pme_f_ratio(gmx_domdec_t *dd)
5160 return dd->comm->load[0].pme/dd->comm->load[0].mdf;
5163 static void dd_print_load(FILE *fplog,gmx_domdec_t *dd,gmx_large_int_t step)
5168 flags = dd_load_flags(dd);
5172 "DD load balancing is limited by minimum cell size in dimension");
5173 for(d=0; d<dd->ndim; d++)
5177 fprintf(fplog," %c",dim2char(dd->dim[d]));
5180 fprintf(fplog,"\n");
5182 fprintf(fplog,"DD step %s",gmx_step_str(step,buf));
5183 if (dd->comm->bDynLoadBal)
5185 fprintf(fplog," vol min/aver %5.3f%c",
5186 dd_vol_min(dd),flags ? '!' : ' ');
5188 fprintf(fplog," load imb.: force %4.1f%%",dd_f_imbal(dd)*100);
5189 if (dd->comm->cycl_n[ddCyclPME])
5191 fprintf(fplog," pme mesh/force %5.3f",dd_pme_f_ratio(dd));
5193 fprintf(fplog,"\n\n");
5196 static void dd_print_load_verbose(gmx_domdec_t *dd)
5198 if (dd->comm->bDynLoadBal)
5200 fprintf(stderr,"vol %4.2f%c ",
5201 dd_vol_min(dd),dd_load_flags(dd) ? '!' : ' ');
5203 fprintf(stderr,"imb F %2d%% ",(int)(dd_f_imbal(dd)*100+0.5));
5204 if (dd->comm->cycl_n[ddCyclPME])
5206 fprintf(stderr,"pme/F %4.2f ",dd_pme_f_ratio(dd));
5211 static void make_load_communicator(gmx_domdec_t *dd,MPI_Group g_all,
5212 int dim_ind,ivec loc)
5214 MPI_Group g_row = MPI_GROUP_EMPTY;
5218 gmx_domdec_root_t *root;
5219 gmx_bool bPartOfGroup = FALSE;
5221 dim = dd->dim[dim_ind];
5222 copy_ivec(loc,loc_c);
5223 snew(rank,dd->nc[dim]);
5224 for(i=0; i<dd->nc[dim]; i++)
5227 rank[i] = dd_index(dd->nc,loc_c);
5228 if (rank[i] == dd->rank)
5230 /* This process is part of the group */
5231 bPartOfGroup = TRUE;
5236 MPI_Group_incl(g_all,dd->nc[dim],rank,&g_row);
5238 MPI_Comm_create(dd->mpi_comm_all,g_row,&c_row);
5241 dd->comm->mpi_comm_load[dim_ind] = c_row;
5242 if (dd->comm->eDLB != edlbNO)
5244 if (dd->ci[dim] == dd->master_ci[dim])
5246 /* This is the root process of this row */
5247 snew(dd->comm->root[dim_ind],1);
5248 root = dd->comm->root[dim_ind];
5249 snew(root->cell_f,DD_CELL_F_SIZE(dd,dim_ind));
5250 snew(root->old_cell_f,dd->nc[dim]+1);
5251 snew(root->bCellMin,dd->nc[dim]);
5254 snew(root->cell_f_max0,dd->nc[dim]);
5255 snew(root->cell_f_min1,dd->nc[dim]);
5256 snew(root->bound_min,dd->nc[dim]);
5257 snew(root->bound_max,dd->nc[dim]);
5259 snew(root->buf_ncd,dd->nc[dim]);
5263 /* This is not a root process, we only need to receive cell_f */
5264 snew(dd->comm->cell_f_row,DD_CELL_F_SIZE(dd,dim_ind));
5267 if (dd->ci[dim] == dd->master_ci[dim])
5269 snew(dd->comm->load[dim_ind].load,dd->nc[dim]*DD_NLOAD_MAX);
5276 static void make_load_communicators(gmx_domdec_t *dd)
5284 fprintf(debug,"Making load communicators\n");
5286 MPI_Comm_group(dd->mpi_comm_all,&g_all);
5288 snew(dd->comm->load,dd->ndim);
5289 snew(dd->comm->mpi_comm_load,dd->ndim);
5292 make_load_communicator(dd,g_all,0,loc);
5295 for(i=0; i<dd->nc[dim0]; i++) {
5297 make_load_communicator(dd,g_all,1,loc);
5302 for(i=0; i<dd->nc[dim0]; i++) {
5305 for(j=0; j<dd->nc[dim1]; j++) {
5307 make_load_communicator(dd,g_all,2,loc);
5312 MPI_Group_free(&g_all);
5315 fprintf(debug,"Finished making load communicators\n");
5319 void setup_dd_grid(FILE *fplog,gmx_domdec_t *dd)
5325 ivec dd_zp[DD_MAXIZONE];
5326 gmx_domdec_zones_t *zones;
5327 gmx_domdec_ns_ranges_t *izone;
5329 for(d=0; d<dd->ndim; d++)
5332 copy_ivec(dd->ci,tmp);
5333 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5334 dd->neighbor[d][0] = ddcoord2ddnodeid(dd,tmp);
5335 copy_ivec(dd->ci,tmp);
5336 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5337 dd->neighbor[d][1] = ddcoord2ddnodeid(dd,tmp);
5340 fprintf(debug,"DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5343 dd->neighbor[d][1]);
5349 fprintf(stderr,"Making %dD domain decomposition %d x %d x %d\n",
5350 dd->ndim,dd->nc[XX],dd->nc[YY],dd->nc[ZZ]);
5354 fprintf(fplog,"\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5356 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],
5357 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5364 for(i=0; i<nzonep; i++)
5366 copy_ivec(dd_zp3[i],dd_zp[i]);
5372 for(i=0; i<nzonep; i++)
5374 copy_ivec(dd_zp2[i],dd_zp[i]);
5380 for(i=0; i<nzonep; i++)
5382 copy_ivec(dd_zp1[i],dd_zp[i]);
5386 gmx_fatal(FARGS,"Can only do 1, 2 or 3D domain decomposition");
5391 zones = &dd->comm->zones;
5393 for(i=0; i<nzone; i++)
5396 clear_ivec(zones->shift[i]);
5397 for(d=0; d<dd->ndim; d++)
5399 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5404 for(i=0; i<nzone; i++)
5406 for(d=0; d<DIM; d++)
5408 s[d] = dd->ci[d] - zones->shift[i][d];
5413 else if (s[d] >= dd->nc[d])
5419 zones->nizone = nzonep;
5420 for(i=0; i<zones->nizone; i++)
5422 if (dd_zp[i][0] != i)
5424 gmx_fatal(FARGS,"Internal inconsistency in the dd grid setup");
5426 izone = &zones->izone[i];
5427 izone->j0 = dd_zp[i][1];
5428 izone->j1 = dd_zp[i][2];
5429 for(dim=0; dim<DIM; dim++)
5431 if (dd->nc[dim] == 1)
5433 /* All shifts should be allowed */
5434 izone->shift0[dim] = -1;
5435 izone->shift1[dim] = 1;
5440 izone->shift0[d] = 0;
5441 izone->shift1[d] = 0;
5442 for(j=izone->j0; j<izone->j1; j++) {
5443 if (dd->shift[j][d] > dd->shift[i][d])
5444 izone->shift0[d] = -1;
5445 if (dd->shift[j][d] < dd->shift[i][d])
5446 izone->shift1[d] = 1;
5452 /* Assume the shift are not more than 1 cell */
5453 izone->shift0[dim] = 1;
5454 izone->shift1[dim] = -1;
5455 for(j=izone->j0; j<izone->j1; j++)
5457 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5458 if (shift_diff < izone->shift0[dim])
5460 izone->shift0[dim] = shift_diff;
5462 if (shift_diff > izone->shift1[dim])
5464 izone->shift1[dim] = shift_diff;
5471 if (dd->comm->eDLB != edlbNO)
5473 snew(dd->comm->root,dd->ndim);
5476 if (dd->comm->bRecordLoad)
5478 make_load_communicators(dd);
5482 static void make_pp_communicator(FILE *fplog,t_commrec *cr,int reorder)
5485 gmx_domdec_comm_t *comm;
5496 if (comm->bCartesianPP)
5498 /* Set up cartesian communication for the particle-particle part */
5501 fprintf(fplog,"Will use a Cartesian communicator: %d x %d x %d\n",
5502 dd->nc[XX],dd->nc[YY],dd->nc[ZZ]);
5505 for(i=0; i<DIM; i++)
5509 MPI_Cart_create(cr->mpi_comm_mygroup,DIM,dd->nc,periods,reorder,
5511 /* We overwrite the old communicator with the new cartesian one */
5512 cr->mpi_comm_mygroup = comm_cart;
5515 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5516 MPI_Comm_rank(dd->mpi_comm_all,&dd->rank);
5518 if (comm->bCartesianPP_PME)
5520 /* Since we want to use the original cartesian setup for sim,
5521 * and not the one after split, we need to make an index.
5523 snew(comm->ddindex2ddnodeid,dd->nnodes);
5524 comm->ddindex2ddnodeid[dd_index(dd->nc,dd->ci)] = dd->rank;
5525 gmx_sumi(dd->nnodes,comm->ddindex2ddnodeid,cr);
5526 /* Get the rank of the DD master,
5527 * above we made sure that the master node is a PP node.
5537 MPI_Allreduce(&rank,&dd->masterrank,1,MPI_INT,MPI_SUM,dd->mpi_comm_all);
5539 else if (comm->bCartesianPP)
5541 if (cr->npmenodes == 0)
5543 /* The PP communicator is also
5544 * the communicator for this simulation
5546 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5548 cr->nodeid = dd->rank;
5550 MPI_Cart_coords(dd->mpi_comm_all,dd->rank,DIM,dd->ci);
5552 /* We need to make an index to go from the coordinates
5553 * to the nodeid of this simulation.
5555 snew(comm->ddindex2simnodeid,dd->nnodes);
5556 snew(buf,dd->nnodes);
5557 if (cr->duty & DUTY_PP)
5559 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5561 /* Communicate the ddindex to simulation nodeid index */
5562 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5563 cr->mpi_comm_mysim);
5566 /* Determine the master coordinates and rank.
5567 * The DD master should be the same node as the master of this sim.
5569 for(i=0; i<dd->nnodes; i++)
5571 if (comm->ddindex2simnodeid[i] == 0)
5573 ddindex2xyz(dd->nc,i,dd->master_ci);
5574 MPI_Cart_rank(dd->mpi_comm_all,dd->master_ci,&dd->masterrank);
5579 fprintf(debug,"The master rank is %d\n",dd->masterrank);
5584 /* No Cartesian communicators */
5585 /* We use the rank in dd->comm->all as DD index */
5586 ddindex2xyz(dd->nc,dd->rank,dd->ci);
5587 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5589 clear_ivec(dd->master_ci);
5596 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5597 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5602 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5603 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5607 static void receive_ddindex2simnodeid(t_commrec *cr)
5611 gmx_domdec_comm_t *comm;
5618 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5620 snew(comm->ddindex2simnodeid,dd->nnodes);
5621 snew(buf,dd->nnodes);
5622 if (cr->duty & DUTY_PP)
5624 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5627 /* Communicate the ddindex to simulation nodeid index */
5628 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5629 cr->mpi_comm_mysim);
5636 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5639 gmx_domdec_master_t *ma;
5644 snew(ma->ncg,dd->nnodes);
5645 snew(ma->index,dd->nnodes+1);
5647 snew(ma->nat,dd->nnodes);
5648 snew(ma->ibuf,dd->nnodes*2);
5649 snew(ma->cell_x,DIM);
5650 for(i=0; i<DIM; i++)
5652 snew(ma->cell_x[i],dd->nc[i]+1);
5655 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5661 snew(ma->vbuf,natoms);
5667 static void split_communicator(FILE *fplog,t_commrec *cr,int dd_node_order,
5671 gmx_domdec_comm_t *comm;
5682 if (comm->bCartesianPP)
5684 for(i=1; i<DIM; i++)
5686 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5688 if (bDiv[YY] || bDiv[ZZ])
5690 comm->bCartesianPP_PME = TRUE;
5691 /* If we have 2D PME decomposition, which is always in x+y,
5692 * we stack the PME only nodes in z.
5693 * Otherwise we choose the direction that provides the thinnest slab
5694 * of PME only nodes as this will have the least effect
5695 * on the PP communication.
5696 * But for the PME communication the opposite might be better.
5698 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
5700 dd->nc[YY] > dd->nc[ZZ]))
5702 comm->cartpmedim = ZZ;
5706 comm->cartpmedim = YY;
5708 comm->ntot[comm->cartpmedim]
5709 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
5713 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]);
5715 "Will not use a Cartesian communicator for PP <-> PME\n\n");
5720 if (comm->bCartesianPP_PME)
5724 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]);
5727 for(i=0; i<DIM; i++)
5731 MPI_Cart_create(cr->mpi_comm_mysim,DIM,comm->ntot,periods,reorder,
5734 MPI_Comm_rank(comm_cart,&rank);
5735 if (MASTERNODE(cr) && rank != 0)
5737 gmx_fatal(FARGS,"MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
5740 /* With this assigment we loose the link to the original communicator
5741 * which will usually be MPI_COMM_WORLD, unless have multisim.
5743 cr->mpi_comm_mysim = comm_cart;
5744 cr->sim_nodeid = rank;
5746 MPI_Cart_coords(cr->mpi_comm_mysim,cr->sim_nodeid,DIM,dd->ci);
5750 fprintf(fplog,"Cartesian nodeid %d, coordinates %d %d %d\n\n",
5751 cr->sim_nodeid,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5754 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
5758 if (cr->npmenodes == 0 ||
5759 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
5761 cr->duty = DUTY_PME;
5764 /* Split the sim communicator into PP and PME only nodes */
5765 MPI_Comm_split(cr->mpi_comm_mysim,
5767 dd_index(comm->ntot,dd->ci),
5768 &cr->mpi_comm_mygroup);
5772 switch (dd_node_order)
5777 fprintf(fplog,"Order of the nodes: PP first, PME last\n");
5780 case ddnoINTERLEAVE:
5781 /* Interleave the PP-only and PME-only nodes,
5782 * as on clusters with dual-core machines this will double
5783 * the communication bandwidth of the PME processes
5784 * and thus speed up the PP <-> PME and inter PME communication.
5788 fprintf(fplog,"Interleaving PP and PME nodes\n");
5790 comm->pmenodes = dd_pmenodes(cr);
5795 gmx_fatal(FARGS,"Unknown dd_node_order=%d",dd_node_order);
5798 if (dd_simnode2pmenode(cr,cr->sim_nodeid) == -1)
5800 cr->duty = DUTY_PME;
5807 /* Split the sim communicator into PP and PME only nodes */
5808 MPI_Comm_split(cr->mpi_comm_mysim,
5811 &cr->mpi_comm_mygroup);
5812 MPI_Comm_rank(cr->mpi_comm_mygroup,&cr->nodeid);
5818 fprintf(fplog,"This is a %s only node\n\n",
5819 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
5823 void make_dd_communicators(FILE *fplog,t_commrec *cr,int dd_node_order)
5826 gmx_domdec_comm_t *comm;
5832 copy_ivec(dd->nc,comm->ntot);
5834 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
5835 comm->bCartesianPP_PME = FALSE;
5837 /* Reorder the nodes by default. This might change the MPI ranks.
5838 * Real reordering is only supported on very few architectures,
5839 * Blue Gene is one of them.
5841 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
5843 if (cr->npmenodes > 0)
5845 /* Split the communicator into a PP and PME part */
5846 split_communicator(fplog,cr,dd_node_order,CartReorder);
5847 if (comm->bCartesianPP_PME)
5849 /* We (possibly) reordered the nodes in split_communicator,
5850 * so it is no longer required in make_pp_communicator.
5852 CartReorder = FALSE;
5857 /* All nodes do PP and PME */
5859 /* We do not require separate communicators */
5860 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
5864 if (cr->duty & DUTY_PP)
5866 /* Copy or make a new PP communicator */
5867 make_pp_communicator(fplog,cr,CartReorder);
5871 receive_ddindex2simnodeid(cr);
5874 if (!(cr->duty & DUTY_PME))
5876 /* Set up the commnuication to our PME node */
5877 dd->pme_nodeid = dd_simnode2pmenode(cr,cr->sim_nodeid);
5878 dd->pme_receive_vir_ener = receive_vir_ener(cr);
5881 fprintf(debug,"My pme_nodeid %d receive ener %d\n",
5882 dd->pme_nodeid,dd->pme_receive_vir_ener);
5887 dd->pme_nodeid = -1;
5892 dd->ma = init_gmx_domdec_master_t(dd,
5894 comm->cgs_gl.index[comm->cgs_gl.nr]);
5898 static real *get_slb_frac(FILE *fplog,const char *dir,int nc,const char *size_string)
5905 if (nc > 1 && size_string != NULL)
5909 fprintf(fplog,"Using static load balancing for the %s direction\n",
5914 for (i=0; i<nc; i++)
5917 sscanf(size_string,"%lf%n",&dbl,&n);
5920 gmx_fatal(FARGS,"Incorrect or not enough DD cell size entries for direction %s: '%s'",dir,size_string);
5929 fprintf(fplog,"Relative cell sizes:");
5931 for (i=0; i<nc; i++)
5936 fprintf(fplog," %5.3f",slb_frac[i]);
5941 fprintf(fplog,"\n");
5948 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
5951 gmx_mtop_ilistloop_t iloop;
5955 iloop = gmx_mtop_ilistloop_init(mtop);
5956 while (gmx_mtop_ilistloop_next(iloop,&il,&nmol))
5958 for(ftype=0; ftype<F_NRE; ftype++)
5960 if ((interaction_function[ftype].flags & IF_BOND) &&
5963 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
5971 static int dd_nst_env(FILE *fplog,const char *env_var,int def)
5977 val = getenv(env_var);
5980 if (sscanf(val,"%d",&nst) <= 0)
5986 fprintf(fplog,"Found env.var. %s = %s, using value %d\n",
5994 static void dd_warning(t_commrec *cr,FILE *fplog,const char *warn_string)
5998 fprintf(stderr,"\n%s\n",warn_string);
6002 fprintf(fplog,"\n%s\n",warn_string);
6006 static void check_dd_restrictions(t_commrec *cr,gmx_domdec_t *dd,
6007 t_inputrec *ir,FILE *fplog)
6009 if (ir->ePBC == epbcSCREW &&
6010 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6012 gmx_fatal(FARGS,"With pbc=%s can only do domain decomposition in the x-direction",epbc_names[ir->ePBC]);
6015 if (ir->ns_type == ensSIMPLE)
6017 gmx_fatal(FARGS,"Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6020 if (ir->nstlist == 0)
6022 gmx_fatal(FARGS,"Domain decomposition does not work with nstlist=0");
6025 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6027 dd_warning(cr,fplog,"comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6031 static real average_cellsize_min(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
6036 r = ddbox->box_size[XX];
6037 for(di=0; di<dd->ndim; di++)
6040 /* Check using the initial average cell size */
6041 r = min(r,ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6047 static int check_dlb_support(FILE *fplog,t_commrec *cr,
6048 const char *dlb_opt,gmx_bool bRecordLoad,
6049 unsigned long Flags,t_inputrec *ir)
6057 case 'a': eDLB = edlbAUTO; break;
6058 case 'n': eDLB = edlbNO; break;
6059 case 'y': eDLB = edlbYES; break;
6060 default: gmx_incons("Unknown dlb_opt");
6063 if (Flags & MD_RERUN)
6068 if (!EI_DYNAMICS(ir->eI))
6070 if (eDLB == edlbYES)
6072 sprintf(buf,"NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n",EI(ir->eI));
6073 dd_warning(cr,fplog,buf);
6081 dd_warning(cr,fplog,"NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6086 if (Flags & MD_REPRODUCIBLE)
6093 dd_warning(cr,fplog,"NOTE: reproducibility requested, will not use dynamic load balancing\n");
6097 dd_warning(cr,fplog,"WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6100 gmx_fatal(FARGS,"Death horror: undefined case (%d) for load balancing choice",eDLB);
6108 static void set_dd_dim(FILE *fplog,gmx_domdec_t *dd)
6113 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6115 /* Decomposition order z,y,x */
6118 fprintf(fplog,"Using domain decomposition order z, y, x\n");
6120 for(dim=DIM-1; dim>=0; dim--)
6122 if (dd->nc[dim] > 1)
6124 dd->dim[dd->ndim++] = dim;
6130 /* Decomposition order x,y,z */
6131 for(dim=0; dim<DIM; dim++)
6133 if (dd->nc[dim] > 1)
6135 dd->dim[dd->ndim++] = dim;
6141 static gmx_domdec_comm_t *init_dd_comm()
6143 gmx_domdec_comm_t *comm;
6147 snew(comm->cggl_flag,DIM*2);
6148 snew(comm->cgcm_state,DIM*2);
6149 for(i=0; i<DIM*2; i++)
6151 comm->cggl_flag_nalloc[i] = 0;
6152 comm->cgcm_state_nalloc[i] = 0;
6155 comm->nalloc_int = 0;
6156 comm->buf_int = NULL;
6158 vec_rvec_init(&comm->vbuf);
6160 comm->n_load_have = 0;
6161 comm->n_load_collect = 0;
6163 for(i=0; i<ddnatNR-ddnatZONE; i++)
6165 comm->sum_nat[i] = 0;
6169 comm->load_step = 0;
6172 clear_ivec(comm->load_lim);
6179 gmx_domdec_t *init_domain_decomposition(FILE *fplog,t_commrec *cr,
6180 unsigned long Flags,
6182 real comm_distance_min,real rconstr,
6183 const char *dlb_opt,real dlb_scale,
6184 const char *sizex,const char *sizey,const char *sizez,
6185 gmx_mtop_t *mtop,t_inputrec *ir,
6188 int *npme_x,int *npme_y)
6191 gmx_domdec_comm_t *comm;
6194 real r_2b,r_mb,r_bonded=-1,r_bonded_limit=-1,limit,acs;
6201 "\nInitializing Domain Decomposition on %d nodes\n",cr->nnodes);
6206 dd->comm = init_dd_comm();
6208 snew(comm->cggl_flag,DIM*2);
6209 snew(comm->cgcm_state,DIM*2);
6211 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6212 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6214 dd->bSendRecv2 = dd_nst_env(fplog,"GMX_DD_SENDRECV2",0);
6215 comm->dlb_scale_lim = dd_nst_env(fplog,"GMX_DLB_MAX",10);
6216 comm->eFlop = dd_nst_env(fplog,"GMX_DLB_FLOP",0);
6217 recload = dd_nst_env(fplog,"GMX_DD_LOAD",1);
6218 comm->nstSortCG = dd_nst_env(fplog,"GMX_DD_SORT",1);
6219 comm->nstDDDump = dd_nst_env(fplog,"GMX_DD_DUMP",0);
6220 comm->nstDDDumpGrid = dd_nst_env(fplog,"GMX_DD_DUMP_GRID",0);
6221 comm->DD_debug = dd_nst_env(fplog,"GMX_DD_DEBUG",0);
6223 dd->pme_recv_f_alloc = 0;
6224 dd->pme_recv_f_buf = NULL;
6226 if (dd->bSendRecv2 && fplog)
6228 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");
6234 fprintf(fplog,"Will load balance based on FLOP count\n");
6236 if (comm->eFlop > 1)
6238 srand(1+cr->nodeid);
6240 comm->bRecordLoad = TRUE;
6244 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6248 comm->eDLB = check_dlb_support(fplog,cr,dlb_opt,comm->bRecordLoad,Flags,ir);
6250 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6253 fprintf(fplog,"Dynamic load balancing: %s\n",edlb_names[comm->eDLB]);
6255 dd->bGridJump = comm->bDynLoadBal;
6257 if (comm->nstSortCG)
6261 if (comm->nstSortCG == 1)
6263 fprintf(fplog,"Will sort the charge groups at every domain (re)decomposition\n");
6267 fprintf(fplog,"Will sort the charge groups every %d steps\n",
6277 fprintf(fplog,"Will not sort the charge groups\n");
6281 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6282 if (comm->bInterCGBondeds)
6284 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6288 comm->bInterCGMultiBody = FALSE;
6291 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6293 if (ir->rlistlong == 0)
6295 /* Set the cut-off to some very large value,
6296 * so we don't need if statements everywhere in the code.
6297 * We use sqrt, since the cut-off is squared in some places.
6299 comm->cutoff = GMX_CUTOFF_INF;
6303 comm->cutoff = ir->rlistlong;
6305 comm->cutoff_mbody = 0;
6307 comm->cellsize_limit = 0;
6308 comm->bBondComm = FALSE;
6310 if (comm->bInterCGBondeds)
6312 if (comm_distance_min > 0)
6314 comm->cutoff_mbody = comm_distance_min;
6315 if (Flags & MD_DDBONDCOMM)
6317 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6321 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6323 r_bonded_limit = comm->cutoff_mbody;
6325 else if (ir->bPeriodicMols)
6327 /* Can not easily determine the required cut-off */
6328 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");
6329 comm->cutoff_mbody = comm->cutoff/2;
6330 r_bonded_limit = comm->cutoff_mbody;
6336 dd_bonded_cg_distance(fplog,dd,mtop,ir,x,box,
6337 Flags & MD_DDBONDCHECK,&r_2b,&r_mb);
6339 gmx_bcast(sizeof(r_2b),&r_2b,cr);
6340 gmx_bcast(sizeof(r_mb),&r_mb,cr);
6342 /* We use an initial margin of 10% for the minimum cell size,
6343 * except when we are just below the non-bonded cut-off.
6345 if (Flags & MD_DDBONDCOMM)
6347 if (max(r_2b,r_mb) > comm->cutoff)
6349 r_bonded = max(r_2b,r_mb);
6350 r_bonded_limit = 1.1*r_bonded;
6351 comm->bBondComm = TRUE;
6356 r_bonded_limit = min(1.1*r_bonded,comm->cutoff);
6358 /* We determine cutoff_mbody later */
6362 /* No special bonded communication,
6363 * simply increase the DD cut-off.
6365 r_bonded_limit = 1.1*max(r_2b,r_mb);
6366 comm->cutoff_mbody = r_bonded_limit;
6367 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6370 comm->cellsize_limit = max(comm->cellsize_limit,r_bonded_limit);
6374 "Minimum cell size due to bonded interactions: %.3f nm\n",
6375 comm->cellsize_limit);
6379 if (dd->bInterCGcons && rconstr <= 0)
6381 /* There is a cell size limit due to the constraints (P-LINCS) */
6382 rconstr = constr_r_max(fplog,mtop,ir);
6386 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6388 if (rconstr > comm->cellsize_limit)
6390 fprintf(fplog,"This distance will limit the DD cell size, you can override this with -rcon\n");
6394 else if (rconstr > 0 && fplog)
6396 /* Here we do not check for dd->bInterCGcons,
6397 * because one can also set a cell size limit for virtual sites only
6398 * and at this point we don't know yet if there are intercg v-sites.
6401 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6404 comm->cellsize_limit = max(comm->cellsize_limit,rconstr);
6406 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6410 copy_ivec(nc,dd->nc);
6411 set_dd_dim(fplog,dd);
6412 set_ddbox_cr(cr,&dd->nc,ir,box,&comm->cgs_gl,x,ddbox);
6414 if (cr->npmenodes == -1)
6418 acs = average_cellsize_min(dd,ddbox);
6419 if (acs < comm->cellsize_limit)
6423 fprintf(fplog,"ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n",acs,comm->cellsize_limit);
6425 gmx_fatal_collective(FARGS,cr,NULL,
6426 "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",
6427 acs,comm->cellsize_limit);
6432 set_ddbox_cr(cr,NULL,ir,box,&comm->cgs_gl,x,ddbox);
6434 /* We need to choose the optimal DD grid and possibly PME nodes */
6435 limit = dd_choose_grid(fplog,cr,dd,ir,mtop,box,ddbox,
6436 comm->eDLB!=edlbNO,dlb_scale,
6437 comm->cellsize_limit,comm->cutoff,
6438 comm->bInterCGBondeds,comm->bInterCGMultiBody);
6440 if (dd->nc[XX] == 0)
6442 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6443 sprintf(buf,"Change the number of nodes or mdrun option %s%s%s",
6444 !bC ? "-rdd" : "-rcon",
6445 comm->eDLB!=edlbNO ? " or -dds" : "",
6446 bC ? " or your LINCS settings" : "");
6448 gmx_fatal_collective(FARGS,cr,NULL,
6449 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6451 "Look in the log file for details on the domain decomposition",
6452 cr->nnodes-cr->npmenodes,limit,buf);
6454 set_dd_dim(fplog,dd);
6460 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6461 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],cr->npmenodes);
6464 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6465 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6467 gmx_fatal_collective(FARGS,cr,NULL,
6468 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6469 dd->nnodes,cr->nnodes - cr->npmenodes,cr->nnodes);
6471 if (cr->npmenodes > dd->nnodes)
6473 gmx_fatal_collective(FARGS,cr,NULL,
6474 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.",cr->npmenodes,dd->nnodes);
6476 if (cr->npmenodes > 0)
6478 comm->npmenodes = cr->npmenodes;
6482 comm->npmenodes = dd->nnodes;
6485 if (EEL_PME(ir->coulombtype))
6487 /* The following choices should match those
6488 * in comm_cost_est in domdec_setup.c.
6489 * Note that here the checks have to take into account
6490 * that the decomposition might occur in a different order than xyz
6491 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6492 * in which case they will not match those in comm_cost_est,
6493 * but since that is mainly for testing purposes that's fine.
6495 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6496 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6497 getenv("GMX_PMEONEDD") == NULL)
6499 comm->npmedecompdim = 2;
6500 comm->npmenodes_x = dd->nc[XX];
6501 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6505 /* In case nc is 1 in both x and y we could still choose to
6506 * decompose pme in y instead of x, but we use x for simplicity.
6508 comm->npmedecompdim = 1;
6509 if (dd->dim[0] == YY)
6511 comm->npmenodes_x = 1;
6512 comm->npmenodes_y = comm->npmenodes;
6516 comm->npmenodes_x = comm->npmenodes;
6517 comm->npmenodes_y = 1;
6522 fprintf(fplog,"PME domain decomposition: %d x %d x %d\n",
6523 comm->npmenodes_x,comm->npmenodes_y,1);
6528 comm->npmedecompdim = 0;
6529 comm->npmenodes_x = 0;
6530 comm->npmenodes_y = 0;
6533 /* Technically we don't need both of these,
6534 * but it simplifies code not having to recalculate it.
6536 *npme_x = comm->npmenodes_x;
6537 *npme_y = comm->npmenodes_y;
6539 snew(comm->slb_frac,DIM);
6540 if (comm->eDLB == edlbNO)
6542 comm->slb_frac[XX] = get_slb_frac(fplog,"x",dd->nc[XX],sizex);
6543 comm->slb_frac[YY] = get_slb_frac(fplog,"y",dd->nc[YY],sizey);
6544 comm->slb_frac[ZZ] = get_slb_frac(fplog,"z",dd->nc[ZZ],sizez);
6547 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6549 if (comm->bBondComm || comm->eDLB != edlbNO)
6551 /* Set the bonded communication distance to halfway
6552 * the minimum and the maximum,
6553 * since the extra communication cost is nearly zero.
6555 acs = average_cellsize_min(dd,ddbox);
6556 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6557 if (comm->eDLB != edlbNO)
6559 /* Check if this does not limit the scaling */
6560 comm->cutoff_mbody = min(comm->cutoff_mbody,dlb_scale*acs);
6562 if (!comm->bBondComm)
6564 /* Without bBondComm do not go beyond the n.b. cut-off */
6565 comm->cutoff_mbody = min(comm->cutoff_mbody,comm->cutoff);
6566 if (comm->cellsize_limit >= comm->cutoff)
6568 /* We don't loose a lot of efficieny
6569 * when increasing it to the n.b. cut-off.
6570 * It can even be slightly faster, because we need
6571 * less checks for the communication setup.
6573 comm->cutoff_mbody = comm->cutoff;
6576 /* Check if we did not end up below our original limit */
6577 comm->cutoff_mbody = max(comm->cutoff_mbody,r_bonded_limit);
6579 if (comm->cutoff_mbody > comm->cellsize_limit)
6581 comm->cellsize_limit = comm->cutoff_mbody;
6584 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6589 fprintf(debug,"Bonded atom communication beyond the cut-off: %d\n"
6590 "cellsize limit %f\n",
6591 comm->bBondComm,comm->cellsize_limit);
6596 check_dd_restrictions(cr,dd,ir,fplog);
6599 comm->globalcomm_step = INT_MIN;
6602 clear_dd_cycle_counts(dd);
6607 static void set_dlb_limits(gmx_domdec_t *dd)
6612 for(d=0; d<dd->ndim; d++)
6614 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6615 dd->comm->cellsize_min[dd->dim[d]] =
6616 dd->comm->cellsize_min_dlb[dd->dim[d]];
6621 static void turn_on_dlb(FILE *fplog,t_commrec *cr,gmx_large_int_t step)
6624 gmx_domdec_comm_t *comm;
6634 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);
6637 cellsize_min = comm->cellsize_min[dd->dim[0]];
6638 for(d=1; d<dd->ndim; d++)
6640 cellsize_min = min(cellsize_min,comm->cellsize_min[dd->dim[d]]);
6643 if (cellsize_min < comm->cellsize_limit*1.05)
6645 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");
6647 /* Change DLB from "auto" to "no". */
6648 comm->eDLB = edlbNO;
6653 dd_warning(cr,fplog,"NOTE: Turning on dynamic load balancing\n");
6654 comm->bDynLoadBal = TRUE;
6655 dd->bGridJump = TRUE;
6659 /* We can set the required cell size info here,
6660 * so we do not need to communicate this.
6661 * The grid is completely uniform.
6663 for(d=0; d<dd->ndim; d++)
6667 comm->load[d].sum_m = comm->load[d].sum;
6669 nc = dd->nc[dd->dim[d]];
6672 comm->root[d]->cell_f[i] = i/(real)nc;
6675 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6676 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6679 comm->root[d]->cell_f[nc] = 1.0;
6684 static char *init_bLocalCG(gmx_mtop_t *mtop)
6689 ncg = ncg_mtop(mtop);
6691 for(cg=0; cg<ncg; cg++)
6693 bLocalCG[cg] = FALSE;
6699 void dd_init_bondeds(FILE *fplog,
6700 gmx_domdec_t *dd,gmx_mtop_t *mtop,
6701 gmx_vsite_t *vsite,gmx_constr_t constr,
6702 t_inputrec *ir,gmx_bool bBCheck,cginfo_mb_t *cginfo_mb)
6704 gmx_domdec_comm_t *comm;
6708 dd_make_reverse_top(fplog,dd,mtop,vsite,constr,ir,bBCheck);
6712 if (comm->bBondComm)
6714 /* Communicate atoms beyond the cut-off for bonded interactions */
6717 comm->cglink = make_charge_group_links(mtop,dd,cginfo_mb);
6719 comm->bLocalCG = init_bLocalCG(mtop);
6723 /* Only communicate atoms based on cut-off */
6724 comm->cglink = NULL;
6725 comm->bLocalCG = NULL;
6729 static void print_dd_settings(FILE *fplog,gmx_domdec_t *dd,
6731 gmx_bool bDynLoadBal,real dlb_scale,
6734 gmx_domdec_comm_t *comm;
6749 fprintf(fplog,"The maximum number of communication pulses is:");
6750 for(d=0; d<dd->ndim; d++)
6752 fprintf(fplog," %c %d",dim2char(dd->dim[d]),comm->cd[d].np_dlb);
6754 fprintf(fplog,"\n");
6755 fprintf(fplog,"The minimum size for domain decomposition cells is %.3f nm\n",comm->cellsize_limit);
6756 fprintf(fplog,"The requested allowed shrink of DD cells (option -dds) is: %.2f\n",dlb_scale);
6757 fprintf(fplog,"The allowed shrink of domain decomposition cells is:");
6758 for(d=0; d<DIM; d++)
6762 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
6769 comm->cellsize_min_dlb[d]/
6770 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6772 fprintf(fplog," %c %.2f",dim2char(d),shrink);
6775 fprintf(fplog,"\n");
6779 set_dd_cell_sizes_slb(dd,ddbox,FALSE,np);
6780 fprintf(fplog,"The initial number of communication pulses is:");
6781 for(d=0; d<dd->ndim; d++)
6783 fprintf(fplog," %c %d",dim2char(dd->dim[d]),np[dd->dim[d]]);
6785 fprintf(fplog,"\n");
6786 fprintf(fplog,"The initial domain decomposition cell size is:");
6787 for(d=0; d<DIM; d++) {
6790 fprintf(fplog," %c %.2f nm",
6791 dim2char(d),dd->comm->cellsize_min[d]);
6794 fprintf(fplog,"\n\n");
6797 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
6799 fprintf(fplog,"The maximum allowed distance for charge groups involved in interactions is:\n");
6800 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6801 "non-bonded interactions","",comm->cutoff);
6805 limit = dd->comm->cellsize_limit;
6809 if (dynamic_dd_box(ddbox,ir))
6811 fprintf(fplog,"(the following are initial values, they could change due to box deformation)\n");
6813 limit = dd->comm->cellsize_min[XX];
6814 for(d=1; d<DIM; d++)
6816 limit = min(limit,dd->comm->cellsize_min[d]);
6820 if (comm->bInterCGBondeds)
6822 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6823 "two-body bonded interactions","(-rdd)",
6824 max(comm->cutoff,comm->cutoff_mbody));
6825 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6826 "multi-body bonded interactions","(-rdd)",
6827 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff,limit));
6831 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6832 "virtual site constructions","(-rcon)",limit);
6834 if (dd->constraint_comm)
6836 sprintf(buf,"atoms separated by up to %d constraints",
6838 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6839 buf,"(-rcon)",limit);
6841 fprintf(fplog,"\n");
6847 void set_dd_parameters(FILE *fplog,gmx_domdec_t *dd,real dlb_scale,
6848 t_inputrec *ir,t_forcerec *fr,
6851 gmx_domdec_comm_t *comm;
6852 int d,dim,npulse,npulse_d_max,npulse_d;
6859 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
6861 if (EEL_PME(ir->coulombtype))
6863 init_ddpme(dd,&comm->ddpme[0],0);
6864 if (comm->npmedecompdim >= 2)
6866 init_ddpme(dd,&comm->ddpme[1],1);
6871 comm->npmenodes = 0;
6872 if (dd->pme_nodeid >= 0)
6874 gmx_fatal_collective(FARGS,NULL,dd,
6875 "Can not have separate PME nodes without PME electrostatics");
6879 /* If each molecule is a single charge group
6880 * or we use domain decomposition for each periodic dimension,
6881 * we do not need to take pbc into account for the bonded interactions.
6883 if (fr->ePBC == epbcNONE || !comm->bInterCGBondeds ||
6884 (dd->nc[XX]>1 && dd->nc[YY]>1 && (dd->nc[ZZ]>1 || fr->ePBC==epbcXY)))
6886 fr->bMolPBC = FALSE;
6895 fprintf(debug,"The DD cut-off is %f\n",comm->cutoff);
6897 if (comm->eDLB != edlbNO)
6899 /* Determine the maximum number of comm. pulses in one dimension */
6901 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
6903 /* Determine the maximum required number of grid pulses */
6904 if (comm->cellsize_limit >= comm->cutoff)
6906 /* Only a single pulse is required */
6909 else if (!bNoCutOff && comm->cellsize_limit > 0)
6911 /* We round down slightly here to avoid overhead due to the latency
6912 * of extra communication calls when the cut-off
6913 * would be only slightly longer than the cell size.
6914 * Later cellsize_limit is redetermined,
6915 * so we can not miss interactions due to this rounding.
6917 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
6921 /* There is no cell size limit */
6922 npulse = max(dd->nc[XX]-1,max(dd->nc[YY]-1,dd->nc[ZZ]-1));
6925 if (!bNoCutOff && npulse > 1)
6927 /* See if we can do with less pulses, based on dlb_scale */
6929 for(d=0; d<dd->ndim; d++)
6932 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
6933 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
6934 npulse_d_max = max(npulse_d_max,npulse_d);
6936 npulse = min(npulse,npulse_d_max);
6939 /* This env var can override npulse */
6940 d = dd_nst_env(fplog,"GMX_DD_NPULSE",0);
6947 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
6948 for(d=0; d<dd->ndim; d++)
6950 comm->cd[d].np_dlb = min(npulse,dd->nc[dd->dim[d]]-1);
6951 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
6952 snew(comm->cd[d].ind,comm->cd[d].np_nalloc);
6953 comm->maxpulse = max(comm->maxpulse,comm->cd[d].np_dlb);
6954 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
6956 comm->bVacDLBNoLimit = FALSE;
6960 /* cellsize_limit is set for LINCS in init_domain_decomposition */
6961 if (!comm->bVacDLBNoLimit)
6963 comm->cellsize_limit = max(comm->cellsize_limit,
6964 comm->cutoff/comm->maxpulse);
6966 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
6967 /* Set the minimum cell size for each DD dimension */
6968 for(d=0; d<dd->ndim; d++)
6970 if (comm->bVacDLBNoLimit ||
6971 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
6973 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
6977 comm->cellsize_min_dlb[dd->dim[d]] =
6978 comm->cutoff/comm->cd[d].np_dlb;
6981 if (comm->cutoff_mbody <= 0)
6983 comm->cutoff_mbody = min(comm->cutoff,comm->cellsize_limit);
6985 if (comm->bDynLoadBal)
6991 print_dd_settings(fplog,dd,ir,comm->bDynLoadBal,dlb_scale,ddbox);
6992 if (comm->eDLB == edlbAUTO)
6996 fprintf(fplog,"When dynamic load balancing gets turned on, these settings will change to:\n");
6998 print_dd_settings(fplog,dd,ir,TRUE,dlb_scale,ddbox);
7001 if (ir->ePBC == epbcNONE)
7003 vol_frac = 1 - 1/(double)dd->nnodes;
7008 (1 + comm_box_frac(dd->nc,comm->cutoff,ddbox))/(double)dd->nnodes;
7012 fprintf(debug,"Volume fraction for all DD zones: %f\n",vol_frac);
7014 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7016 dd->ga2la = ga2la_init(natoms_tot,vol_frac*natoms_tot);
7019 static void merge_cg_buffers(int ncell,
7020 gmx_domdec_comm_dim_t *cd, int pulse,
7022 int *index_gl, int *recv_i,
7023 rvec *cg_cm, rvec *recv_vr,
7025 cginfo_mb_t *cginfo_mb,int *cginfo)
7027 gmx_domdec_ind_t *ind,*ind_p;
7028 int p,cell,c,cg,cg0,cg1,cg_gl,nat;
7031 ind = &cd->ind[pulse];
7033 /* First correct the already stored data */
7034 shift = ind->nrecv[ncell];
7035 for(cell=ncell-1; cell>=0; cell--)
7037 shift -= ind->nrecv[cell];
7040 /* Move the cg's present from previous grid pulses */
7041 cg0 = ncg_cell[ncell+cell];
7042 cg1 = ncg_cell[ncell+cell+1];
7043 cgindex[cg1+shift] = cgindex[cg1];
7044 for(cg=cg1-1; cg>=cg0; cg--)
7046 index_gl[cg+shift] = index_gl[cg];
7047 copy_rvec(cg_cm[cg],cg_cm[cg+shift]);
7048 cgindex[cg+shift] = cgindex[cg];
7049 cginfo[cg+shift] = cginfo[cg];
7051 /* Correct the already stored send indices for the shift */
7052 for(p=1; p<=pulse; p++)
7054 ind_p = &cd->ind[p];
7056 for(c=0; c<cell; c++)
7058 cg0 += ind_p->nsend[c];
7060 cg1 = cg0 + ind_p->nsend[cell];
7061 for(cg=cg0; cg<cg1; cg++)
7063 ind_p->index[cg] += shift;
7069 /* Merge in the communicated buffers */
7073 for(cell=0; cell<ncell; cell++)
7075 cg1 = ncg_cell[ncell+cell+1] + shift;
7078 /* Correct the old cg indices */
7079 for(cg=ncg_cell[ncell+cell]; cg<cg1; cg++)
7081 cgindex[cg+1] += shift_at;
7084 for(cg=0; cg<ind->nrecv[cell]; cg++)
7086 /* Copy this charge group from the buffer */
7087 index_gl[cg1] = recv_i[cg0];
7088 copy_rvec(recv_vr[cg0],cg_cm[cg1]);
7089 /* Add it to the cgindex */
7090 cg_gl = index_gl[cg1];
7091 cginfo[cg1] = ddcginfo(cginfo_mb,cg_gl);
7092 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7093 cgindex[cg1+1] = cgindex[cg1] + nat;
7098 shift += ind->nrecv[cell];
7099 ncg_cell[ncell+cell+1] = cg1;
7103 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7104 int nzone,int cg0,const int *cgindex)
7108 /* Store the atom block boundaries for easy copying of communication buffers
7111 for(zone=0; zone<nzone; zone++)
7113 for(p=0; p<cd->np; p++) {
7114 cd->ind[p].cell2at0[zone] = cgindex[cg];
7115 cg += cd->ind[p].nrecv[zone];
7116 cd->ind[p].cell2at1[zone] = cgindex[cg];
7121 static gmx_bool missing_link(t_blocka *link,int cg_gl,char *bLocalCG)
7127 for(i=link->index[cg_gl]; i<link->index[cg_gl+1]; i++)
7129 if (!bLocalCG[link->a[i]])
7138 static void setup_dd_communication(gmx_domdec_t *dd,
7139 matrix box,gmx_ddbox_t *ddbox,t_forcerec *fr)
7141 int dim_ind,dim,dim0,dim1=-1,dim2=-1,dimd,p,nat_tot;
7142 int nzone,nzone_send,zone,zonei,cg0,cg1;
7143 int c,i,j,cg,cg_gl,nrcg;
7144 int *zone_cg_range,pos_cg,*index_gl,*cgindex,*recv_i;
7145 gmx_domdec_comm_t *comm;
7146 gmx_domdec_zones_t *zones;
7147 gmx_domdec_comm_dim_t *cd;
7148 gmx_domdec_ind_t *ind;
7149 cginfo_mb_t *cginfo_mb;
7150 gmx_bool bBondComm,bDist2B,bDistMB,bDistMB_pulse,bDistBonded,bScrew;
7151 real r_mb,r_comm2,r_scomm2,r_bcomm2,r,r_0,r_1,r2,rb2,r2inc,inv_ncg,tric_sh;
7153 real corner[DIM][4],corner_round_0=0,corner_round_1[4];
7154 real bcorner[DIM],bcorner_round_1=0;
7156 rvec *cg_cm,*normal,*v_d,*v_0=NULL,*v_1=NULL,*recv_vr;
7157 real skew_fac2_d,skew_fac_01;
7163 fprintf(debug,"Setting up DD communication\n");
7169 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7171 dim = dd->dim[dim_ind];
7173 /* Check if we need to use triclinic distances */
7174 tric_dist[dim_ind] = 0;
7175 for(i=0; i<=dim_ind; i++)
7177 if (ddbox->tric_dir[dd->dim[i]])
7179 tric_dist[dim_ind] = 1;
7184 bBondComm = comm->bBondComm;
7186 /* Do we need to determine extra distances for multi-body bondeds? */
7187 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
7189 /* Do we need to determine extra distances for only two-body bondeds? */
7190 bDist2B = (bBondComm && !bDistMB);
7192 r_comm2 = sqr(comm->cutoff);
7193 r_bcomm2 = sqr(comm->cutoff_mbody);
7197 fprintf(debug,"bBondComm %d, r_bc %f\n",bBondComm,sqrt(r_bcomm2));
7200 zones = &comm->zones;
7203 /* The first dimension is equal for all cells */
7204 corner[0][0] = comm->cell_x0[dim0];
7207 bcorner[0] = corner[0][0];
7212 /* This cell row is only seen from the first row */
7213 corner[1][0] = comm->cell_x0[dim1];
7214 /* All rows can see this row */
7215 corner[1][1] = comm->cell_x0[dim1];
7218 corner[1][1] = max(comm->cell_x0[dim1],comm->zone_d1[1].mch0);
7221 /* For the multi-body distance we need the maximum */
7222 bcorner[1] = max(comm->cell_x0[dim1],comm->zone_d1[1].p1_0);
7225 /* Set the upper-right corner for rounding */
7226 corner_round_0 = comm->cell_x1[dim0];
7233 corner[2][j] = comm->cell_x0[dim2];
7237 /* Use the maximum of the i-cells that see a j-cell */
7238 for(i=0; i<zones->nizone; i++)
7240 for(j=zones->izone[i].j0; j<zones->izone[i].j1; j++)
7246 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7252 /* For the multi-body distance we need the maximum */
7253 bcorner[2] = comm->cell_x0[dim2];
7258 bcorner[2] = max(bcorner[2],
7259 comm->zone_d2[i][j].p1_0);
7265 /* Set the upper-right corner for rounding */
7266 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7267 * Only cell (0,0,0) can see cell 7 (1,1,1)
7269 corner_round_1[0] = comm->cell_x1[dim1];
7270 corner_round_1[3] = comm->cell_x1[dim1];
7273 corner_round_1[0] = max(comm->cell_x1[dim1],
7274 comm->zone_d1[1].mch1);
7277 /* For the multi-body distance we need the maximum */
7278 bcorner_round_1 = max(comm->cell_x1[dim1],
7279 comm->zone_d1[1].p1_1);
7285 /* Triclinic stuff */
7286 normal = ddbox->normal;
7290 v_0 = ddbox->v[dim0];
7291 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
7293 /* Determine the coupling coefficient for the distances
7294 * to the cell planes along dim0 and dim1 through dim2.
7295 * This is required for correct rounding.
7298 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
7301 fprintf(debug,"\nskew_fac_01 %f\n",skew_fac_01);
7307 v_1 = ddbox->v[dim1];
7310 zone_cg_range = zones->cg_range;
7311 index_gl = dd->index_gl;
7312 cgindex = dd->cgindex;
7313 cginfo_mb = fr->cginfo_mb;
7315 zone_cg_range[0] = 0;
7316 zone_cg_range[1] = dd->ncg_home;
7317 comm->zone_ncg1[0] = dd->ncg_home;
7318 pos_cg = dd->ncg_home;
7320 nat_tot = dd->nat_home;
7322 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7324 dim = dd->dim[dim_ind];
7325 cd = &comm->cd[dim_ind];
7327 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
7329 /* No pbc in this dimension, the first node should not comm. */
7337 bScrew = (dd->bScrewPBC && dim == XX);
7339 v_d = ddbox->v[dim];
7340 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
7342 cd->bInPlace = TRUE;
7343 for(p=0; p<cd->np; p++)
7345 /* Only atoms communicated in the first pulse are used
7346 * for multi-body bonded interactions or for bBondComm.
7348 bDistBonded = ((bDistMB || bDist2B) && p == 0);
7349 bDistMB_pulse = (bDistMB && bDistBonded);
7354 for(zone=0; zone<nzone_send; zone++)
7356 if (tric_dist[dim_ind] && dim_ind > 0)
7358 /* Determine slightly more optimized skew_fac's
7360 * This reduces the number of communicated atoms
7361 * by about 10% for 3D DD of rhombic dodecahedra.
7363 for(dimd=0; dimd<dim; dimd++)
7365 sf2_round[dimd] = 1;
7366 if (ddbox->tric_dir[dimd])
7368 for(i=dd->dim[dimd]+1; i<DIM; i++)
7370 /* If we are shifted in dimension i
7371 * and the cell plane is tilted forward
7372 * in dimension i, skip this coupling.
7374 if (!(zones->shift[nzone+zone][i] &&
7375 ddbox->v[dimd][i][dimd] >= 0))
7378 sqr(ddbox->v[dimd][i][dimd]);
7381 sf2_round[dimd] = 1/sf2_round[dimd];
7386 zonei = zone_perm[dim_ind][zone];
7389 /* Here we permutate the zones to obtain a convenient order
7390 * for neighbor searching
7392 cg0 = zone_cg_range[zonei];
7393 cg1 = zone_cg_range[zonei+1];
7397 /* Look only at the cg's received in the previous grid pulse
7399 cg1 = zone_cg_range[nzone+zone+1];
7400 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
7402 ind->nsend[zone] = 0;
7403 for(cg=cg0; cg<cg1; cg++)
7407 if (tric_dist[dim_ind] == 0)
7409 /* Rectangular direction, easy */
7410 r = cg_cm[cg][dim] - corner[dim_ind][zone];
7417 r = cg_cm[cg][dim] - bcorner[dim_ind];
7423 /* Rounding gives at most a 16% reduction
7424 * in communicated atoms
7426 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7428 r = cg_cm[cg][dim0] - corner_round_0;
7429 /* This is the first dimension, so always r >= 0 */
7436 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7438 r = cg_cm[cg][dim1] - corner_round_1[zone];
7445 r = cg_cm[cg][dim1] - bcorner_round_1;
7455 /* Triclinic direction, more complicated */
7458 /* Rounding, conservative as the skew_fac multiplication
7459 * will slightly underestimate the distance.
7461 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7463 rn[dim0] = cg_cm[cg][dim0] - corner_round_0;
7464 for(i=dim0+1; i<DIM; i++)
7466 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7468 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7471 rb[dim0] = rn[dim0];
7474 /* Take care that the cell planes along dim0 might not
7475 * be orthogonal to those along dim1 and dim2.
7477 for(i=1; i<=dim_ind; i++)
7480 if (normal[dim0][dimd] > 0)
7482 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7485 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7490 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7492 rn[dim1] += cg_cm[cg][dim1] - corner_round_1[zone];
7494 for(i=dim1+1; i<DIM; i++)
7496 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7498 rn[dim1] += tric_sh;
7501 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7502 /* Take care of coupling of the distances
7503 * to the planes along dim0 and dim1 through dim2.
7505 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7506 /* Take care that the cell planes along dim1
7507 * might not be orthogonal to that along dim2.
7509 if (normal[dim1][dim2] > 0)
7511 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7517 cg_cm[cg][dim1] - bcorner_round_1 + tric_sh;
7520 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7521 /* Take care of coupling of the distances
7522 * to the planes along dim0 and dim1 through dim2.
7524 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7525 /* Take care that the cell planes along dim1
7526 * might not be orthogonal to that along dim2.
7528 if (normal[dim1][dim2] > 0)
7530 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7535 /* The distance along the communication direction */
7536 rn[dim] += cg_cm[cg][dim] - corner[dim_ind][zone];
7538 for(i=dim+1; i<DIM; i++)
7540 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7545 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7546 /* Take care of coupling of the distances
7547 * to the planes along dim0 and dim1 through dim2.
7549 if (dim_ind == 1 && zonei == 1)
7551 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7557 rb[dim] += cg_cm[cg][dim] - bcorner[dim_ind] + tric_sh;
7560 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7561 /* Take care of coupling of the distances
7562 * to the planes along dim0 and dim1 through dim2.
7564 if (dim_ind == 1 && zonei == 1)
7566 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7574 ((bDistMB && rb2 < r_bcomm2) ||
7575 (bDist2B && r2 < r_bcomm2)) &&
7577 (GET_CGINFO_BOND_INTER(fr->cginfo[cg]) &&
7578 missing_link(comm->cglink,index_gl[cg],
7581 /* Make an index to the local charge groups */
7582 if (nsend+1 > ind->nalloc)
7584 ind->nalloc = over_alloc_large(nsend+1);
7585 srenew(ind->index,ind->nalloc);
7587 if (nsend+1 > comm->nalloc_int)
7589 comm->nalloc_int = over_alloc_large(nsend+1);
7590 srenew(comm->buf_int,comm->nalloc_int);
7592 ind->index[nsend] = cg;
7593 comm->buf_int[nsend] = index_gl[cg];
7595 vec_rvec_check_alloc(&comm->vbuf,nsend+1);
7597 if (dd->ci[dim] == 0)
7599 /* Correct cg_cm for pbc */
7600 rvec_add(cg_cm[cg],box[dim],comm->vbuf.v[nsend]);
7603 comm->vbuf.v[nsend][YY] =
7604 box[YY][YY]-comm->vbuf.v[nsend][YY];
7605 comm->vbuf.v[nsend][ZZ] =
7606 box[ZZ][ZZ]-comm->vbuf.v[nsend][ZZ];
7611 copy_rvec(cg_cm[cg],comm->vbuf.v[nsend]);
7614 nat += cgindex[cg+1] - cgindex[cg];
7618 /* Clear the counts in case we do not have pbc */
7619 for(zone=nzone_send; zone<nzone; zone++)
7621 ind->nsend[zone] = 0;
7623 ind->nsend[nzone] = nsend;
7624 ind->nsend[nzone+1] = nat;
7625 /* Communicate the number of cg's and atoms to receive */
7626 dd_sendrecv_int(dd, dim_ind, dddirBackward,
7627 ind->nsend, nzone+2,
7628 ind->nrecv, nzone+2);
7630 /* The rvec buffer is also required for atom buffers of size nsend
7631 * in dd_move_x and dd_move_f.
7633 vec_rvec_check_alloc(&comm->vbuf,ind->nsend[nzone+1]);
7637 /* We can receive in place if only the last zone is not empty */
7638 for(zone=0; zone<nzone-1; zone++)
7640 if (ind->nrecv[zone] > 0)
7642 cd->bInPlace = FALSE;
7647 /* The int buffer is only required here for the cg indices */
7648 if (ind->nrecv[nzone] > comm->nalloc_int2)
7650 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
7651 srenew(comm->buf_int2,comm->nalloc_int2);
7653 /* The rvec buffer is also required for atom buffers
7654 * of size nrecv in dd_move_x and dd_move_f.
7656 i = max(cd->ind[0].nrecv[nzone+1],ind->nrecv[nzone+1]);
7657 vec_rvec_check_alloc(&comm->vbuf2,i);
7661 /* Make space for the global cg indices */
7662 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
7663 || dd->cg_nalloc == 0)
7665 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
7666 srenew(index_gl,dd->cg_nalloc);
7667 srenew(cgindex,dd->cg_nalloc+1);
7669 /* Communicate the global cg indices */
7672 recv_i = index_gl + pos_cg;
7676 recv_i = comm->buf_int2;
7678 dd_sendrecv_int(dd, dim_ind, dddirBackward,
7679 comm->buf_int, nsend,
7680 recv_i, ind->nrecv[nzone]);
7682 /* Make space for cg_cm */
7683 if (pos_cg + ind->nrecv[nzone] > fr->cg_nalloc)
7685 dd_realloc_fr_cg(fr,pos_cg + ind->nrecv[nzone]);
7688 /* Communicate cg_cm */
7691 recv_vr = cg_cm + pos_cg;
7695 recv_vr = comm->vbuf2.v;
7697 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
7698 comm->vbuf.v, nsend,
7699 recv_vr, ind->nrecv[nzone]);
7701 /* Make the charge group index */
7704 zone = (p == 0 ? 0 : nzone - 1);
7705 while (zone < nzone)
7707 for(cg=0; cg<ind->nrecv[zone]; cg++)
7709 cg_gl = index_gl[pos_cg];
7710 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb,cg_gl);
7711 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
7712 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
7715 /* Update the charge group presence,
7716 * so we can use it in the next pass of the loop.
7718 comm->bLocalCG[cg_gl] = TRUE;
7724 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
7727 zone_cg_range[nzone+zone] = pos_cg;
7732 /* This part of the code is never executed with bBondComm. */
7733 merge_cg_buffers(nzone,cd,p,zone_cg_range,
7734 index_gl,recv_i,cg_cm,recv_vr,
7735 cgindex,fr->cginfo_mb,fr->cginfo);
7736 pos_cg += ind->nrecv[nzone];
7738 nat_tot += ind->nrecv[nzone+1];
7742 /* Store the atom block for easy copying of communication buffers */
7743 make_cell2at_index(cd,nzone,zone_cg_range[nzone],cgindex);
7747 dd->index_gl = index_gl;
7748 dd->cgindex = cgindex;
7750 dd->ncg_tot = zone_cg_range[zones->n];
7751 dd->nat_tot = nat_tot;
7752 comm->nat[ddnatHOME] = dd->nat_home;
7753 for(i=ddnatZONE; i<ddnatNR; i++)
7755 comm->nat[i] = dd->nat_tot;
7760 /* We don't need to update cginfo, since that was alrady done above.
7761 * So we pass NULL for the forcerec.
7763 dd_set_cginfo(dd->index_gl,dd->ncg_home,dd->ncg_tot,
7764 NULL,comm->bLocalCG);
7769 fprintf(debug,"Finished setting up DD communication, zones:");
7770 for(c=0; c<zones->n; c++)
7772 fprintf(debug," %d",zones->cg_range[c+1]-zones->cg_range[c]);
7774 fprintf(debug,"\n");
7778 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
7782 for(c=0; c<zones->nizone; c++)
7784 zones->izone[c].cg1 = zones->cg_range[c+1];
7785 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
7786 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
7790 static int comp_cgsort(const void *a,const void *b)
7794 gmx_cgsort_t *cga,*cgb;
7795 cga = (gmx_cgsort_t *)a;
7796 cgb = (gmx_cgsort_t *)b;
7798 comp = cga->nsc - cgb->nsc;
7801 comp = cga->ind_gl - cgb->ind_gl;
7807 static void order_int_cg(int n,gmx_cgsort_t *sort,
7812 /* Order the data */
7815 buf[i] = a[sort[i].ind];
7818 /* Copy back to the original array */
7825 static void order_vec_cg(int n,gmx_cgsort_t *sort,
7830 /* Order the data */
7833 copy_rvec(v[sort[i].ind],buf[i]);
7836 /* Copy back to the original array */
7839 copy_rvec(buf[i],v[i]);
7843 static void order_vec_atom(int ncg,int *cgindex,gmx_cgsort_t *sort,
7846 int a,atot,cg,cg0,cg1,i;
7848 /* Order the data */
7850 for(cg=0; cg<ncg; cg++)
7852 cg0 = cgindex[sort[cg].ind];
7853 cg1 = cgindex[sort[cg].ind+1];
7854 for(i=cg0; i<cg1; i++)
7856 copy_rvec(v[i],buf[a]);
7862 /* Copy back to the original array */
7863 for(a=0; a<atot; a++)
7865 copy_rvec(buf[a],v[a]);
7869 static void ordered_sort(int nsort2,gmx_cgsort_t *sort2,
7870 int nsort_new,gmx_cgsort_t *sort_new,
7871 gmx_cgsort_t *sort1)
7875 /* The new indices are not very ordered, so we qsort them */
7876 qsort_threadsafe(sort_new,nsort_new,sizeof(sort_new[0]),comp_cgsort);
7878 /* sort2 is already ordered, so now we can merge the two arrays */
7882 while(i2 < nsort2 || i_new < nsort_new)
7886 sort1[i1++] = sort_new[i_new++];
7888 else if (i_new == nsort_new)
7890 sort1[i1++] = sort2[i2++];
7892 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
7893 (sort2[i2].nsc == sort_new[i_new].nsc &&
7894 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
7896 sort1[i1++] = sort2[i2++];
7900 sort1[i1++] = sort_new[i_new++];
7905 static void dd_sort_state(gmx_domdec_t *dd,int ePBC,
7906 rvec *cgcm,t_forcerec *fr,t_state *state,
7909 gmx_domdec_sort_t *sort;
7910 gmx_cgsort_t *cgsort,*sort_i;
7911 int ncg_new,nsort2,nsort_new,i,cell_index,*ibuf,cgsize;
7914 sort = dd->comm->sort;
7916 if (dd->ncg_home > sort->sort_nalloc)
7918 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
7919 srenew(sort->sort1,sort->sort_nalloc);
7920 srenew(sort->sort2,sort->sort_nalloc);
7923 if (ncg_home_old >= 0)
7925 /* The charge groups that remained in the same ns grid cell
7926 * are completely ordered. So we can sort efficiently by sorting
7927 * the charge groups that did move into the stationary list.
7932 for(i=0; i<dd->ncg_home; i++)
7934 /* Check if this cg did not move to another node */
7935 cell_index = fr->ns.grid->cell_index[i];
7936 if (cell_index != 4*fr->ns.grid->ncells)
7938 if (i >= ncg_home_old || cell_index != sort->sort1[i].nsc)
7940 /* This cg is new on this node or moved ns grid cell */
7941 if (nsort_new >= sort->sort_new_nalloc)
7943 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
7944 srenew(sort->sort_new,sort->sort_new_nalloc);
7946 sort_i = &(sort->sort_new[nsort_new++]);
7950 /* This cg did not move */
7951 sort_i = &(sort->sort2[nsort2++]);
7953 /* Sort on the ns grid cell indices
7954 * and the global topology index
7956 sort_i->nsc = cell_index;
7957 sort_i->ind_gl = dd->index_gl[i];
7964 fprintf(debug,"ordered sort cgs: stationary %d moved %d\n",
7967 /* Sort efficiently */
7968 ordered_sort(nsort2,sort->sort2,nsort_new,sort->sort_new,sort->sort1);
7972 cgsort = sort->sort1;
7974 for(i=0; i<dd->ncg_home; i++)
7976 /* Sort on the ns grid cell indices
7977 * and the global topology index
7979 cgsort[i].nsc = fr->ns.grid->cell_index[i];
7980 cgsort[i].ind_gl = dd->index_gl[i];
7982 if (cgsort[i].nsc != 4*fr->ns.grid->ncells)
7989 fprintf(debug,"qsort cgs: %d new home %d\n",dd->ncg_home,ncg_new);
7991 /* Determine the order of the charge groups using qsort */
7992 qsort_threadsafe(cgsort,dd->ncg_home,sizeof(cgsort[0]),comp_cgsort);
7994 cgsort = sort->sort1;
7996 /* We alloc with the old size, since cgindex is still old */
7997 vec_rvec_check_alloc(&dd->comm->vbuf,dd->cgindex[dd->ncg_home]);
7998 vbuf = dd->comm->vbuf.v;
8000 /* Remove the charge groups which are no longer at home here */
8001 dd->ncg_home = ncg_new;
8003 /* Reorder the state */
8004 for(i=0; i<estNR; i++)
8006 if (EST_DISTR(i) && (state->flags & (1<<i)))
8011 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->x,vbuf);
8014 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->v,vbuf);
8017 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->sd_X,vbuf);
8020 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->cg_p,vbuf);
8024 case estDISRE_INITF:
8025 case estDISRE_RM3TAV:
8026 case estORIRE_INITF:
8028 /* No ordering required */
8031 gmx_incons("Unknown state entry encountered in dd_sort_state");
8037 order_vec_cg(dd->ncg_home,cgsort,cgcm,vbuf);
8039 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8041 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8042 srenew(sort->ibuf,sort->ibuf_nalloc);
8045 /* Reorder the global cg index */
8046 order_int_cg(dd->ncg_home,cgsort,dd->index_gl,ibuf);
8047 /* Reorder the cginfo */
8048 order_int_cg(dd->ncg_home,cgsort,fr->cginfo,ibuf);
8049 /* Rebuild the local cg index */
8051 for(i=0; i<dd->ncg_home; i++)
8053 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8054 ibuf[i+1] = ibuf[i] + cgsize;
8056 for(i=0; i<dd->ncg_home+1; i++)
8058 dd->cgindex[i] = ibuf[i];
8060 /* Set the home atom number */
8061 dd->nat_home = dd->cgindex[dd->ncg_home];
8063 /* Copy the sorted ns cell indices back to the ns grid struct */
8064 for(i=0; i<dd->ncg_home; i++)
8066 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
8068 fr->ns.grid->nr = dd->ncg_home;
8071 static void add_dd_statistics(gmx_domdec_t *dd)
8073 gmx_domdec_comm_t *comm;
8078 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8080 comm->sum_nat[ddnat-ddnatZONE] +=
8081 comm->nat[ddnat] - comm->nat[ddnat-1];
8086 void reset_dd_statistics_counters(gmx_domdec_t *dd)
8088 gmx_domdec_comm_t *comm;
8093 /* Reset all the statistics and counters for total run counting */
8094 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8096 comm->sum_nat[ddnat-ddnatZONE] = 0;
8100 comm->load_step = 0;
8103 clear_ivec(comm->load_lim);
8108 void print_dd_statistics(t_commrec *cr,t_inputrec *ir,FILE *fplog)
8110 gmx_domdec_comm_t *comm;
8114 comm = cr->dd->comm;
8116 gmx_sumd(ddnatNR-ddnatZONE,comm->sum_nat,cr);
8123 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");
8125 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8127 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
8132 " av. #atoms communicated per step for force: %d x %.1f\n",
8136 if (cr->dd->vsite_comm)
8139 " av. #atoms communicated per step for vsites: %d x %.1f\n",
8140 (EEL_PME(ir->coulombtype) || ir->coulombtype==eelEWALD) ? 3 : 2,
8145 if (cr->dd->constraint_comm)
8148 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
8149 1 + ir->nLincsIter,av);
8153 gmx_incons(" Unknown type for DD statistics");
8156 fprintf(fplog,"\n");
8158 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
8160 print_dd_load_av(fplog,cr->dd);
8164 void dd_partition_system(FILE *fplog,
8165 gmx_large_int_t step,
8167 gmx_bool bMasterState,
8169 t_state *state_global,
8170 gmx_mtop_t *top_global,
8172 t_state *state_local,
8175 gmx_localtop_t *top_local,
8178 gmx_shellfc_t shellfc,
8179 gmx_constr_t constr,
8181 gmx_wallcycle_t wcycle,
8185 gmx_domdec_comm_t *comm;
8186 gmx_ddbox_t ddbox={0};
8188 gmx_large_int_t step_pcoupl;
8189 rvec cell_ns_x0,cell_ns_x1;
8190 int i,j,n,cg0=0,ncg_home_old=-1,nat_f_novirsum;
8191 gmx_bool bBoxChanged,bNStGlobalComm,bDoDLB,bCheckDLB,bTurnOnDLB,bLogLoad;
8192 gmx_bool bRedist,bSortCG,bResortAll;
8200 bBoxChanged = (bMasterState || DEFORM(*ir));
8201 if (ir->epc != epcNO)
8203 /* With nstpcouple > 1 pressure coupling happens.
8204 * one step after calculating the pressure.
8205 * Box scaling happens at the end of the MD step,
8206 * after the DD partitioning.
8207 * We therefore have to do DLB in the first partitioning
8208 * after an MD step where P-coupling occured.
8209 * We need to determine the last step in which p-coupling occurred.
8210 * MRS -- need to validate this for vv?
8215 step_pcoupl = step - 1;
8219 step_pcoupl = ((step - 1)/n)*n + 1;
8221 if (step_pcoupl >= comm->globalcomm_step)
8227 bNStGlobalComm = (step >= comm->globalcomm_step + nstglobalcomm);
8229 if (!comm->bDynLoadBal)
8235 /* Should we do dynamic load balacing this step?
8236 * Since it requires (possibly expensive) global communication,
8237 * we might want to do DLB less frequently.
8239 if (bBoxChanged || ir->epc != epcNO)
8241 bDoDLB = bBoxChanged;
8245 bDoDLB = bNStGlobalComm;
8249 /* Check if we have recorded loads on the nodes */
8250 if (comm->bRecordLoad && dd_load_count(comm))
8252 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
8254 /* Check if we should use DLB at the second partitioning
8255 * and every 100 partitionings,
8256 * so the extra communication cost is negligible.
8258 n = max(100,nstglobalcomm);
8259 bCheckDLB = (comm->n_load_collect == 0 ||
8260 comm->n_load_have % n == n-1);
8267 /* Print load every nstlog, first and last step to the log file */
8268 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
8269 comm->n_load_collect == 0 ||
8271 (step + ir->nstlist > ir->init_step + ir->nsteps)));
8273 /* Avoid extra communication due to verbose screen output
8274 * when nstglobalcomm is set.
8276 if (bDoDLB || bLogLoad || bCheckDLB ||
8277 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
8279 get_load_distribution(dd,wcycle);
8284 dd_print_load(fplog,dd,step-1);
8288 dd_print_load_verbose(dd);
8291 comm->n_load_collect++;
8294 /* Since the timings are node dependent, the master decides */
8298 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
8301 fprintf(debug,"step %s, imb loss %f\n",
8302 gmx_step_str(step,sbuf),
8303 dd_force_imb_perf_loss(dd));
8306 dd_bcast(dd,sizeof(bTurnOnDLB),&bTurnOnDLB);
8309 turn_on_dlb(fplog,cr,step);
8314 comm->n_load_have++;
8317 cgs_gl = &comm->cgs_gl;
8322 /* Clear the old state */
8323 clear_dd_indices(dd,0,0);
8325 set_ddbox(dd,bMasterState,cr,ir,state_global->box,
8326 TRUE,cgs_gl,state_global->x,&ddbox);
8328 get_cg_distribution(fplog,step,dd,cgs_gl,
8329 state_global->box,&ddbox,state_global->x);
8331 dd_distribute_state(dd,cgs_gl,
8332 state_global,state_local,f);
8334 dd_make_local_cgs(dd,&top_local->cgs);
8336 if (dd->ncg_home > fr->cg_nalloc)
8338 dd_realloc_fr_cg(fr,dd->ncg_home);
8340 calc_cgcm(fplog,0,dd->ncg_home,
8341 &top_local->cgs,state_local->x,fr->cg_cm);
8343 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
8345 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
8349 else if (state_local->ddp_count != dd->ddp_count)
8351 if (state_local->ddp_count > dd->ddp_count)
8353 gmx_fatal(FARGS,"Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)",state_local->ddp_count,dd->ddp_count);
8356 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
8358 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);
8361 /* Clear the old state */
8362 clear_dd_indices(dd,0,0);
8364 /* Build the new indices */
8365 rebuild_cgindex(dd,cgs_gl->index,state_local);
8366 make_dd_indices(dd,cgs_gl->index,0);
8368 /* Redetermine the cg COMs */
8369 calc_cgcm(fplog,0,dd->ncg_home,
8370 &top_local->cgs,state_local->x,fr->cg_cm);
8372 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
8374 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
8376 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
8377 TRUE,&top_local->cgs,state_local->x,&ddbox);
8379 bRedist = comm->bDynLoadBal;
8383 /* We have the full state, only redistribute the cgs */
8385 /* Clear the non-home indices */
8386 clear_dd_indices(dd,dd->ncg_home,dd->nat_home);
8388 /* Avoid global communication for dim's without pbc and -gcom */
8389 if (!bNStGlobalComm)
8391 copy_rvec(comm->box0 ,ddbox.box0 );
8392 copy_rvec(comm->box_size,ddbox.box_size);
8394 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
8395 bNStGlobalComm,&top_local->cgs,state_local->x,&ddbox);
8400 /* For dim's without pbc and -gcom */
8401 copy_rvec(ddbox.box0 ,comm->box0 );
8402 copy_rvec(ddbox.box_size,comm->box_size);
8404 set_dd_cell_sizes(dd,&ddbox,dynamic_dd_box(&ddbox,ir),bMasterState,bDoDLB,
8407 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
8409 write_dd_grid_pdb("dd_grid",step,dd,state_local->box,&ddbox);
8412 /* Check if we should sort the charge groups */
8413 if (comm->nstSortCG > 0)
8415 bSortCG = (bMasterState ||
8416 (bRedist && (step % comm->nstSortCG == 0)));
8423 ncg_home_old = dd->ncg_home;
8427 cg0 = dd_redistribute_cg(fplog,step,dd,ddbox.tric_dir,
8428 state_local,f,fr,mdatoms,
8432 get_nsgrid_boundaries(fr->ns.grid,dd,
8433 state_local->box,&ddbox,&comm->cell_x0,&comm->cell_x1,
8434 dd->ncg_home,fr->cg_cm,
8435 cell_ns_x0,cell_ns_x1,&grid_density);
8439 comm_dd_ns_cell_sizes(dd,&ddbox,cell_ns_x0,cell_ns_x1,step);
8442 copy_ivec(fr->ns.grid->n,ncells_old);
8443 grid_first(fplog,fr->ns.grid,dd,&ddbox,fr->ePBC,
8444 state_local->box,cell_ns_x0,cell_ns_x1,
8445 fr->rlistlong,grid_density);
8446 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
8447 copy_ivec(ddbox.tric_dir,comm->tric_dir);
8451 /* Sort the state on charge group position.
8452 * This enables exact restarts from this step.
8453 * It also improves performance by about 15% with larger numbers
8454 * of atoms per node.
8457 /* Fill the ns grid with the home cell,
8458 * so we can sort with the indices.
8460 set_zones_ncg_home(dd);
8461 fill_grid(fplog,&comm->zones,fr->ns.grid,dd->ncg_home,
8462 0,dd->ncg_home,fr->cg_cm);
8464 /* Check if we can user the old order and ns grid cell indices
8465 * of the charge groups to sort the charge groups efficiently.
8467 bResortAll = (bMasterState ||
8468 fr->ns.grid->n[XX] != ncells_old[XX] ||
8469 fr->ns.grid->n[YY] != ncells_old[YY] ||
8470 fr->ns.grid->n[ZZ] != ncells_old[ZZ]);
8474 fprintf(debug,"Step %s, sorting the %d home charge groups\n",
8475 gmx_step_str(step,sbuf),dd->ncg_home);
8477 dd_sort_state(dd,ir->ePBC,fr->cg_cm,fr,state_local,
8478 bResortAll ? -1 : ncg_home_old);
8479 /* Rebuild all the indices */
8481 ga2la_clear(dd->ga2la);
8484 /* Setup up the communication and communicate the coordinates */
8485 setup_dd_communication(dd,state_local->box,&ddbox,fr);
8487 /* Set the indices */
8488 make_dd_indices(dd,cgs_gl->index,cg0);
8490 /* Set the charge group boundaries for neighbor searching */
8491 set_cg_boundaries(&comm->zones);
8494 write_dd_pdb("dd_home",step,"dump",top_global,cr,
8495 -1,state_local->x,state_local->box);
8498 /* Extract a local topology from the global topology */
8499 for(i=0; i<dd->ndim; i++)
8501 np[dd->dim[i]] = comm->cd[i].np;
8503 dd_make_local_top(fplog,dd,&comm->zones,dd->npbcdim,state_local->box,
8504 comm->cellsize_min,np,
8505 fr,vsite,top_global,top_local);
8507 /* Set up the special atom communication */
8508 n = comm->nat[ddnatZONE];
8509 for(i=ddnatZONE+1; i<ddnatNR; i++)
8514 if (vsite && vsite->n_intercg_vsite)
8516 n = dd_make_local_vsites(dd,n,top_local->idef.il);
8520 if (dd->bInterCGcons)
8522 /* Only for inter-cg constraints we need special code */
8523 n = dd_make_local_constraints(dd,n,top_global,
8524 constr,ir->nProjOrder,
8525 &top_local->idef.il[F_CONSTR]);
8529 gmx_incons("Unknown special atom type setup");
8534 /* Make space for the extra coordinates for virtual site
8535 * or constraint communication.
8537 state_local->natoms = comm->nat[ddnatNR-1];
8538 if (state_local->natoms > state_local->nalloc)
8540 dd_realloc_state(state_local,f,state_local->natoms);
8543 if (fr->bF_NoVirSum)
8545 if (vsite && vsite->n_intercg_vsite)
8547 nat_f_novirsum = comm->nat[ddnatVSITE];
8551 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
8553 nat_f_novirsum = dd->nat_tot;
8557 nat_f_novirsum = dd->nat_home;
8566 /* Set the number of atoms required for the force calculation.
8567 * Forces need to be constrained when using a twin-range setup
8568 * or with energy minimization. For simple simulations we could
8569 * avoid some allocation, zeroing and copying, but this is
8570 * probably not worth the complications ande checking.
8572 forcerec_set_ranges(fr,dd->ncg_home,dd->ncg_tot,
8573 dd->nat_tot,comm->nat[ddnatCON],nat_f_novirsum);
8575 /* We make the all mdatoms up to nat_tot_con.
8576 * We could save some work by only setting invmass
8577 * between nat_tot and nat_tot_con.
8579 /* This call also sets the new number of home particles to dd->nat_home */
8580 atoms2md(top_global,ir,
8581 comm->nat[ddnatCON],dd->gatindex,0,dd->nat_home,mdatoms);
8583 /* Now we have the charges we can sort the FE interactions */
8584 dd_sort_local_top(dd,mdatoms,top_local);
8588 /* Make the local shell stuff, currently no communication is done */
8589 make_local_shells(cr,mdatoms,shellfc);
8592 if (ir->implicit_solvent)
8594 make_local_gb(cr,fr->born,ir->gb_algorithm);
8597 if (!(cr->duty & DUTY_PME))
8599 /* Send the charges to our PME only node */
8600 gmx_pme_send_q(cr,mdatoms->nChargePerturbed,
8601 mdatoms->chargeA,mdatoms->chargeB,
8602 dd_pme_maxshift_x(dd),dd_pme_maxshift_y(dd));
8607 set_constraints(constr,top_local,ir,mdatoms,cr);
8610 if (ir->ePull != epullNO)
8612 /* Update the local pull groups */
8613 dd_make_local_pull_groups(dd,ir->pull,mdatoms);
8618 /* Update the local rotation groups */
8619 dd_make_local_rotation_groups(dd,ir->rot);
8623 add_dd_statistics(dd);
8625 /* Make sure we only count the cycles for this DD partitioning */
8626 clear_dd_cycle_counts(dd);
8628 /* Because the order of the atoms might have changed since
8629 * the last vsite construction, we need to communicate the constructing
8630 * atom coordinates again (for spreading the forces this MD step).
8632 dd_move_x_vsites(dd,state_local->box,state_local->x);
8634 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
8636 dd_move_x(dd,state_local->box,state_local->x);
8637 write_dd_pdb("dd_dump",step,"dump",top_global,cr,
8638 -1,state_local->x,state_local->box);
8643 /* Store the global communication step */
8644 comm->globalcomm_step = step;
8647 /* Increase the DD partitioning counter */
8649 /* The state currently matches this DD partitioning count, store it */
8650 state_local->ddp_count = dd->ddp_count;
8653 /* The DD master node knows the complete cg distribution,
8654 * store the count so we can possibly skip the cg info communication.
8656 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
8659 if (comm->DD_debug > 0)
8661 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
8662 check_index_consistency(dd,top_global->natoms,ncg_mtop(top_global),
8663 "after partitioning");