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, int dim_ind,ivec loc)
5216 gmx_domdec_root_t *root;
5217 gmx_bool bPartOfGroup = FALSE;
5219 dim = dd->dim[dim_ind];
5220 copy_ivec(loc,loc_c);
5221 for(i=0; i<dd->nc[dim]; i++)
5224 rank = dd_index(dd->nc,loc_c);
5225 if (rank == dd->rank)
5227 /* This process is part of the group */
5228 bPartOfGroup = TRUE;
5231 MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup?0:MPI_UNDEFINED, dd->rank,
5235 dd->comm->mpi_comm_load[dim_ind] = c_row;
5236 if (dd->comm->eDLB != edlbNO)
5238 if (dd->ci[dim] == dd->master_ci[dim])
5240 /* This is the root process of this row */
5241 snew(dd->comm->root[dim_ind],1);
5242 root = dd->comm->root[dim_ind];
5243 snew(root->cell_f,DD_CELL_F_SIZE(dd,dim_ind));
5244 snew(root->old_cell_f,dd->nc[dim]+1);
5245 snew(root->bCellMin,dd->nc[dim]);
5248 snew(root->cell_f_max0,dd->nc[dim]);
5249 snew(root->cell_f_min1,dd->nc[dim]);
5250 snew(root->bound_min,dd->nc[dim]);
5251 snew(root->bound_max,dd->nc[dim]);
5253 snew(root->buf_ncd,dd->nc[dim]);
5257 /* This is not a root process, we only need to receive cell_f */
5258 snew(dd->comm->cell_f_row,DD_CELL_F_SIZE(dd,dim_ind));
5261 if (dd->ci[dim] == dd->master_ci[dim])
5263 snew(dd->comm->load[dim_ind].load,dd->nc[dim]*DD_NLOAD_MAX);
5269 static void make_load_communicators(gmx_domdec_t *dd)
5276 fprintf(debug,"Making load communicators\n");
5278 snew(dd->comm->load,dd->ndim);
5279 snew(dd->comm->mpi_comm_load,dd->ndim);
5282 make_load_communicator(dd,0,loc);
5285 for(i=0; i<dd->nc[dim0]; i++) {
5287 make_load_communicator(dd,1,loc);
5292 for(i=0; i<dd->nc[dim0]; i++) {
5295 for(j=0; j<dd->nc[dim1]; j++) {
5297 make_load_communicator(dd,2,loc);
5303 fprintf(debug,"Finished making load communicators\n");
5307 void setup_dd_grid(FILE *fplog,gmx_domdec_t *dd)
5313 ivec dd_zp[DD_MAXIZONE];
5314 gmx_domdec_zones_t *zones;
5315 gmx_domdec_ns_ranges_t *izone;
5317 for(d=0; d<dd->ndim; d++)
5320 copy_ivec(dd->ci,tmp);
5321 tmp[dim] = (tmp[dim] + 1) % dd->nc[dim];
5322 dd->neighbor[d][0] = ddcoord2ddnodeid(dd,tmp);
5323 copy_ivec(dd->ci,tmp);
5324 tmp[dim] = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
5325 dd->neighbor[d][1] = ddcoord2ddnodeid(dd,tmp);
5328 fprintf(debug,"DD rank %d neighbor ranks in dir %d are + %d - %d\n",
5331 dd->neighbor[d][1]);
5337 fprintf(stderr,"Making %dD domain decomposition %d x %d x %d\n",
5338 dd->ndim,dd->nc[XX],dd->nc[YY],dd->nc[ZZ]);
5342 fprintf(fplog,"\nMaking %dD domain decomposition grid %d x %d x %d, home cell index %d %d %d\n\n",
5344 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],
5345 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5352 for(i=0; i<nzonep; i++)
5354 copy_ivec(dd_zp3[i],dd_zp[i]);
5360 for(i=0; i<nzonep; i++)
5362 copy_ivec(dd_zp2[i],dd_zp[i]);
5368 for(i=0; i<nzonep; i++)
5370 copy_ivec(dd_zp1[i],dd_zp[i]);
5374 gmx_fatal(FARGS,"Can only do 1, 2 or 3D domain decomposition");
5379 zones = &dd->comm->zones;
5381 for(i=0; i<nzone; i++)
5384 clear_ivec(zones->shift[i]);
5385 for(d=0; d<dd->ndim; d++)
5387 zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
5392 for(i=0; i<nzone; i++)
5394 for(d=0; d<DIM; d++)
5396 s[d] = dd->ci[d] - zones->shift[i][d];
5401 else if (s[d] >= dd->nc[d])
5407 zones->nizone = nzonep;
5408 for(i=0; i<zones->nizone; i++)
5410 if (dd_zp[i][0] != i)
5412 gmx_fatal(FARGS,"Internal inconsistency in the dd grid setup");
5414 izone = &zones->izone[i];
5415 izone->j0 = dd_zp[i][1];
5416 izone->j1 = dd_zp[i][2];
5417 for(dim=0; dim<DIM; dim++)
5419 if (dd->nc[dim] == 1)
5421 /* All shifts should be allowed */
5422 izone->shift0[dim] = -1;
5423 izone->shift1[dim] = 1;
5428 izone->shift0[d] = 0;
5429 izone->shift1[d] = 0;
5430 for(j=izone->j0; j<izone->j1; j++) {
5431 if (dd->shift[j][d] > dd->shift[i][d])
5432 izone->shift0[d] = -1;
5433 if (dd->shift[j][d] < dd->shift[i][d])
5434 izone->shift1[d] = 1;
5440 /* Assume the shift are not more than 1 cell */
5441 izone->shift0[dim] = 1;
5442 izone->shift1[dim] = -1;
5443 for(j=izone->j0; j<izone->j1; j++)
5445 shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
5446 if (shift_diff < izone->shift0[dim])
5448 izone->shift0[dim] = shift_diff;
5450 if (shift_diff > izone->shift1[dim])
5452 izone->shift1[dim] = shift_diff;
5459 if (dd->comm->eDLB != edlbNO)
5461 snew(dd->comm->root,dd->ndim);
5464 if (dd->comm->bRecordLoad)
5466 make_load_communicators(dd);
5470 static void make_pp_communicator(FILE *fplog,t_commrec *cr,int reorder)
5473 gmx_domdec_comm_t *comm;
5484 if (comm->bCartesianPP)
5486 /* Set up cartesian communication for the particle-particle part */
5489 fprintf(fplog,"Will use a Cartesian communicator: %d x %d x %d\n",
5490 dd->nc[XX],dd->nc[YY],dd->nc[ZZ]);
5493 for(i=0; i<DIM; i++)
5497 MPI_Cart_create(cr->mpi_comm_mygroup,DIM,dd->nc,periods,reorder,
5499 /* We overwrite the old communicator with the new cartesian one */
5500 cr->mpi_comm_mygroup = comm_cart;
5503 dd->mpi_comm_all = cr->mpi_comm_mygroup;
5504 MPI_Comm_rank(dd->mpi_comm_all,&dd->rank);
5506 if (comm->bCartesianPP_PME)
5508 /* Since we want to use the original cartesian setup for sim,
5509 * and not the one after split, we need to make an index.
5511 snew(comm->ddindex2ddnodeid,dd->nnodes);
5512 comm->ddindex2ddnodeid[dd_index(dd->nc,dd->ci)] = dd->rank;
5513 gmx_sumi(dd->nnodes,comm->ddindex2ddnodeid,cr);
5514 /* Get the rank of the DD master,
5515 * above we made sure that the master node is a PP node.
5525 MPI_Allreduce(&rank,&dd->masterrank,1,MPI_INT,MPI_SUM,dd->mpi_comm_all);
5527 else if (comm->bCartesianPP)
5529 if (cr->npmenodes == 0)
5531 /* The PP communicator is also
5532 * the communicator for this simulation
5534 cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
5536 cr->nodeid = dd->rank;
5538 MPI_Cart_coords(dd->mpi_comm_all,dd->rank,DIM,dd->ci);
5540 /* We need to make an index to go from the coordinates
5541 * to the nodeid of this simulation.
5543 snew(comm->ddindex2simnodeid,dd->nnodes);
5544 snew(buf,dd->nnodes);
5545 if (cr->duty & DUTY_PP)
5547 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5549 /* Communicate the ddindex to simulation nodeid index */
5550 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5551 cr->mpi_comm_mysim);
5554 /* Determine the master coordinates and rank.
5555 * The DD master should be the same node as the master of this sim.
5557 for(i=0; i<dd->nnodes; i++)
5559 if (comm->ddindex2simnodeid[i] == 0)
5561 ddindex2xyz(dd->nc,i,dd->master_ci);
5562 MPI_Cart_rank(dd->mpi_comm_all,dd->master_ci,&dd->masterrank);
5567 fprintf(debug,"The master rank is %d\n",dd->masterrank);
5572 /* No Cartesian communicators */
5573 /* We use the rank in dd->comm->all as DD index */
5574 ddindex2xyz(dd->nc,dd->rank,dd->ci);
5575 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
5577 clear_ivec(dd->master_ci);
5584 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5585 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5590 "Domain decomposition nodeid %d, coordinates %d %d %d\n\n",
5591 dd->rank,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5595 static void receive_ddindex2simnodeid(t_commrec *cr)
5599 gmx_domdec_comm_t *comm;
5606 if (!comm->bCartesianPP_PME && comm->bCartesianPP)
5608 snew(comm->ddindex2simnodeid,dd->nnodes);
5609 snew(buf,dd->nnodes);
5610 if (cr->duty & DUTY_PP)
5612 buf[dd_index(dd->nc,dd->ci)] = cr->sim_nodeid;
5615 /* Communicate the ddindex to simulation nodeid index */
5616 MPI_Allreduce(buf,comm->ddindex2simnodeid,dd->nnodes,MPI_INT,MPI_SUM,
5617 cr->mpi_comm_mysim);
5624 static gmx_domdec_master_t *init_gmx_domdec_master_t(gmx_domdec_t *dd,
5627 gmx_domdec_master_t *ma;
5632 snew(ma->ncg,dd->nnodes);
5633 snew(ma->index,dd->nnodes+1);
5635 snew(ma->nat,dd->nnodes);
5636 snew(ma->ibuf,dd->nnodes*2);
5637 snew(ma->cell_x,DIM);
5638 for(i=0; i<DIM; i++)
5640 snew(ma->cell_x[i],dd->nc[i]+1);
5643 if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
5649 snew(ma->vbuf,natoms);
5655 static void split_communicator(FILE *fplog,t_commrec *cr,int dd_node_order,
5659 gmx_domdec_comm_t *comm;
5670 if (comm->bCartesianPP)
5672 for(i=1; i<DIM; i++)
5674 bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
5676 if (bDiv[YY] || bDiv[ZZ])
5678 comm->bCartesianPP_PME = TRUE;
5679 /* If we have 2D PME decomposition, which is always in x+y,
5680 * we stack the PME only nodes in z.
5681 * Otherwise we choose the direction that provides the thinnest slab
5682 * of PME only nodes as this will have the least effect
5683 * on the PP communication.
5684 * But for the PME communication the opposite might be better.
5686 if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
5688 dd->nc[YY] > dd->nc[ZZ]))
5690 comm->cartpmedim = ZZ;
5694 comm->cartpmedim = YY;
5696 comm->ntot[comm->cartpmedim]
5697 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
5701 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]);
5703 "Will not use a Cartesian communicator for PP <-> PME\n\n");
5708 if (comm->bCartesianPP_PME)
5712 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]);
5715 for(i=0; i<DIM; i++)
5719 MPI_Cart_create(cr->mpi_comm_mysim,DIM,comm->ntot,periods,reorder,
5722 MPI_Comm_rank(comm_cart,&rank);
5723 if (MASTERNODE(cr) && rank != 0)
5725 gmx_fatal(FARGS,"MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
5728 /* With this assigment we loose the link to the original communicator
5729 * which will usually be MPI_COMM_WORLD, unless have multisim.
5731 cr->mpi_comm_mysim = comm_cart;
5732 cr->sim_nodeid = rank;
5734 MPI_Cart_coords(cr->mpi_comm_mysim,cr->sim_nodeid,DIM,dd->ci);
5738 fprintf(fplog,"Cartesian nodeid %d, coordinates %d %d %d\n\n",
5739 cr->sim_nodeid,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
5742 if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
5746 if (cr->npmenodes == 0 ||
5747 dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
5749 cr->duty = DUTY_PME;
5752 /* Split the sim communicator into PP and PME only nodes */
5753 MPI_Comm_split(cr->mpi_comm_mysim,
5755 dd_index(comm->ntot,dd->ci),
5756 &cr->mpi_comm_mygroup);
5760 switch (dd_node_order)
5765 fprintf(fplog,"Order of the nodes: PP first, PME last\n");
5768 case ddnoINTERLEAVE:
5769 /* Interleave the PP-only and PME-only nodes,
5770 * as on clusters with dual-core machines this will double
5771 * the communication bandwidth of the PME processes
5772 * and thus speed up the PP <-> PME and inter PME communication.
5776 fprintf(fplog,"Interleaving PP and PME nodes\n");
5778 comm->pmenodes = dd_pmenodes(cr);
5783 gmx_fatal(FARGS,"Unknown dd_node_order=%d",dd_node_order);
5786 if (dd_simnode2pmenode(cr,cr->sim_nodeid) == -1)
5788 cr->duty = DUTY_PME;
5795 /* Split the sim communicator into PP and PME only nodes */
5796 MPI_Comm_split(cr->mpi_comm_mysim,
5799 &cr->mpi_comm_mygroup);
5800 MPI_Comm_rank(cr->mpi_comm_mygroup,&cr->nodeid);
5806 fprintf(fplog,"This is a %s only node\n\n",
5807 (cr->duty & DUTY_PP) ? "particle-particle" : "PME-mesh");
5811 void make_dd_communicators(FILE *fplog,t_commrec *cr,int dd_node_order)
5814 gmx_domdec_comm_t *comm;
5820 copy_ivec(dd->nc,comm->ntot);
5822 comm->bCartesianPP = (dd_node_order == ddnoCARTESIAN);
5823 comm->bCartesianPP_PME = FALSE;
5825 /* Reorder the nodes by default. This might change the MPI ranks.
5826 * Real reordering is only supported on very few architectures,
5827 * Blue Gene is one of them.
5829 CartReorder = (getenv("GMX_NO_CART_REORDER") == NULL);
5831 if (cr->npmenodes > 0)
5833 /* Split the communicator into a PP and PME part */
5834 split_communicator(fplog,cr,dd_node_order,CartReorder);
5835 if (comm->bCartesianPP_PME)
5837 /* We (possibly) reordered the nodes in split_communicator,
5838 * so it is no longer required in make_pp_communicator.
5840 CartReorder = FALSE;
5845 /* All nodes do PP and PME */
5847 /* We do not require separate communicators */
5848 cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
5852 if (cr->duty & DUTY_PP)
5854 /* Copy or make a new PP communicator */
5855 make_pp_communicator(fplog,cr,CartReorder);
5859 receive_ddindex2simnodeid(cr);
5862 if (!(cr->duty & DUTY_PME))
5864 /* Set up the commnuication to our PME node */
5865 dd->pme_nodeid = dd_simnode2pmenode(cr,cr->sim_nodeid);
5866 dd->pme_receive_vir_ener = receive_vir_ener(cr);
5869 fprintf(debug,"My pme_nodeid %d receive ener %d\n",
5870 dd->pme_nodeid,dd->pme_receive_vir_ener);
5875 dd->pme_nodeid = -1;
5880 dd->ma = init_gmx_domdec_master_t(dd,
5882 comm->cgs_gl.index[comm->cgs_gl.nr]);
5886 static real *get_slb_frac(FILE *fplog,const char *dir,int nc,const char *size_string)
5893 if (nc > 1 && size_string != NULL)
5897 fprintf(fplog,"Using static load balancing for the %s direction\n",
5902 for (i=0; i<nc; i++)
5905 sscanf(size_string,"%lf%n",&dbl,&n);
5908 gmx_fatal(FARGS,"Incorrect or not enough DD cell size entries for direction %s: '%s'",dir,size_string);
5917 fprintf(fplog,"Relative cell sizes:");
5919 for (i=0; i<nc; i++)
5924 fprintf(fplog," %5.3f",slb_frac[i]);
5929 fprintf(fplog,"\n");
5936 static int multi_body_bondeds_count(gmx_mtop_t *mtop)
5939 gmx_mtop_ilistloop_t iloop;
5943 iloop = gmx_mtop_ilistloop_init(mtop);
5944 while (gmx_mtop_ilistloop_next(iloop,&il,&nmol))
5946 for(ftype=0; ftype<F_NRE; ftype++)
5948 if ((interaction_function[ftype].flags & IF_BOND) &&
5951 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
5959 static int dd_nst_env(FILE *fplog,const char *env_var,int def)
5965 val = getenv(env_var);
5968 if (sscanf(val,"%d",&nst) <= 0)
5974 fprintf(fplog,"Found env.var. %s = %s, using value %d\n",
5982 static void dd_warning(t_commrec *cr,FILE *fplog,const char *warn_string)
5986 fprintf(stderr,"\n%s\n",warn_string);
5990 fprintf(fplog,"\n%s\n",warn_string);
5994 static void check_dd_restrictions(t_commrec *cr,gmx_domdec_t *dd,
5995 t_inputrec *ir,FILE *fplog)
5997 if (ir->ePBC == epbcSCREW &&
5998 (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
6000 gmx_fatal(FARGS,"With pbc=%s can only do domain decomposition in the x-direction",epbc_names[ir->ePBC]);
6003 if (ir->ns_type == ensSIMPLE)
6005 gmx_fatal(FARGS,"Domain decomposition does not support simple neighbor searching, use grid searching or use particle decomposition");
6008 if (ir->nstlist == 0)
6010 gmx_fatal(FARGS,"Domain decomposition does not work with nstlist=0");
6013 if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
6015 dd_warning(cr,fplog,"comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
6019 static real average_cellsize_min(gmx_domdec_t *dd,gmx_ddbox_t *ddbox)
6024 r = ddbox->box_size[XX];
6025 for(di=0; di<dd->ndim; di++)
6028 /* Check using the initial average cell size */
6029 r = min(r,ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6035 static int check_dlb_support(FILE *fplog,t_commrec *cr,
6036 const char *dlb_opt,gmx_bool bRecordLoad,
6037 unsigned long Flags,t_inputrec *ir)
6045 case 'a': eDLB = edlbAUTO; break;
6046 case 'n': eDLB = edlbNO; break;
6047 case 'y': eDLB = edlbYES; break;
6048 default: gmx_incons("Unknown dlb_opt");
6051 if (Flags & MD_RERUN)
6056 if (!EI_DYNAMICS(ir->eI))
6058 if (eDLB == edlbYES)
6060 sprintf(buf,"NOTE: dynamic load balancing is only supported with dynamics, not with integrator '%s'\n",EI(ir->eI));
6061 dd_warning(cr,fplog,buf);
6069 dd_warning(cr,fplog,"NOTE: Cycle counting is not supported on this architecture, will not use dynamic load balancing\n");
6074 if (Flags & MD_REPRODUCIBLE)
6081 dd_warning(cr,fplog,"NOTE: reproducibility requested, will not use dynamic load balancing\n");
6085 dd_warning(cr,fplog,"WARNING: reproducibility requested with dynamic load balancing, the simulation will NOT be binary reproducible\n");
6088 gmx_fatal(FARGS,"Death horror: undefined case (%d) for load balancing choice",eDLB);
6096 static void set_dd_dim(FILE *fplog,gmx_domdec_t *dd)
6101 if (getenv("GMX_DD_ORDER_ZYX") != NULL)
6103 /* Decomposition order z,y,x */
6106 fprintf(fplog,"Using domain decomposition order z, y, x\n");
6108 for(dim=DIM-1; dim>=0; dim--)
6110 if (dd->nc[dim] > 1)
6112 dd->dim[dd->ndim++] = dim;
6118 /* Decomposition order x,y,z */
6119 for(dim=0; dim<DIM; dim++)
6121 if (dd->nc[dim] > 1)
6123 dd->dim[dd->ndim++] = dim;
6129 static gmx_domdec_comm_t *init_dd_comm()
6131 gmx_domdec_comm_t *comm;
6135 snew(comm->cggl_flag,DIM*2);
6136 snew(comm->cgcm_state,DIM*2);
6137 for(i=0; i<DIM*2; i++)
6139 comm->cggl_flag_nalloc[i] = 0;
6140 comm->cgcm_state_nalloc[i] = 0;
6143 comm->nalloc_int = 0;
6144 comm->buf_int = NULL;
6146 vec_rvec_init(&comm->vbuf);
6148 comm->n_load_have = 0;
6149 comm->n_load_collect = 0;
6151 for(i=0; i<ddnatNR-ddnatZONE; i++)
6153 comm->sum_nat[i] = 0;
6157 comm->load_step = 0;
6160 clear_ivec(comm->load_lim);
6167 gmx_domdec_t *init_domain_decomposition(FILE *fplog,t_commrec *cr,
6168 unsigned long Flags,
6170 real comm_distance_min,real rconstr,
6171 const char *dlb_opt,real dlb_scale,
6172 const char *sizex,const char *sizey,const char *sizez,
6173 gmx_mtop_t *mtop,t_inputrec *ir,
6176 int *npme_x,int *npme_y)
6179 gmx_domdec_comm_t *comm;
6182 real r_2b,r_mb,r_bonded=-1,r_bonded_limit=-1,limit,acs;
6189 "\nInitializing Domain Decomposition on %d nodes\n",cr->nnodes);
6194 dd->comm = init_dd_comm();
6196 snew(comm->cggl_flag,DIM*2);
6197 snew(comm->cgcm_state,DIM*2);
6199 dd->npbcdim = ePBC2npbcdim(ir->ePBC);
6200 dd->bScrewPBC = (ir->ePBC == epbcSCREW);
6202 dd->bSendRecv2 = dd_nst_env(fplog,"GMX_DD_SENDRECV2",0);
6203 comm->dlb_scale_lim = dd_nst_env(fplog,"GMX_DLB_MAX",10);
6204 comm->eFlop = dd_nst_env(fplog,"GMX_DLB_FLOP",0);
6205 recload = dd_nst_env(fplog,"GMX_DD_LOAD",1);
6206 comm->nstSortCG = dd_nst_env(fplog,"GMX_DD_SORT",1);
6207 comm->nstDDDump = dd_nst_env(fplog,"GMX_DD_DUMP",0);
6208 comm->nstDDDumpGrid = dd_nst_env(fplog,"GMX_DD_DUMP_GRID",0);
6209 comm->DD_debug = dd_nst_env(fplog,"GMX_DD_DEBUG",0);
6211 dd->pme_recv_f_alloc = 0;
6212 dd->pme_recv_f_buf = NULL;
6214 if (dd->bSendRecv2 && fplog)
6216 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");
6222 fprintf(fplog,"Will load balance based on FLOP count\n");
6224 if (comm->eFlop > 1)
6226 srand(1+cr->nodeid);
6228 comm->bRecordLoad = TRUE;
6232 comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
6236 comm->eDLB = check_dlb_support(fplog,cr,dlb_opt,comm->bRecordLoad,Flags,ir);
6238 comm->bDynLoadBal = (comm->eDLB == edlbYES);
6241 fprintf(fplog,"Dynamic load balancing: %s\n",edlb_names[comm->eDLB]);
6243 dd->bGridJump = comm->bDynLoadBal;
6245 if (comm->nstSortCG)
6249 if (comm->nstSortCG == 1)
6251 fprintf(fplog,"Will sort the charge groups at every domain (re)decomposition\n");
6255 fprintf(fplog,"Will sort the charge groups every %d steps\n",
6265 fprintf(fplog,"Will not sort the charge groups\n");
6269 comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
6270 if (comm->bInterCGBondeds)
6272 comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
6276 comm->bInterCGMultiBody = FALSE;
6279 dd->bInterCGcons = inter_charge_group_constraints(mtop);
6281 if (ir->rlistlong == 0)
6283 /* Set the cut-off to some very large value,
6284 * so we don't need if statements everywhere in the code.
6285 * We use sqrt, since the cut-off is squared in some places.
6287 comm->cutoff = GMX_CUTOFF_INF;
6291 comm->cutoff = ir->rlistlong;
6293 comm->cutoff_mbody = 0;
6295 comm->cellsize_limit = 0;
6296 comm->bBondComm = FALSE;
6298 if (comm->bInterCGBondeds)
6300 if (comm_distance_min > 0)
6302 comm->cutoff_mbody = comm_distance_min;
6303 if (Flags & MD_DDBONDCOMM)
6305 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
6309 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6311 r_bonded_limit = comm->cutoff_mbody;
6313 else if (ir->bPeriodicMols)
6315 /* Can not easily determine the required cut-off */
6316 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");
6317 comm->cutoff_mbody = comm->cutoff/2;
6318 r_bonded_limit = comm->cutoff_mbody;
6324 dd_bonded_cg_distance(fplog,dd,mtop,ir,x,box,
6325 Flags & MD_DDBONDCHECK,&r_2b,&r_mb);
6327 gmx_bcast(sizeof(r_2b),&r_2b,cr);
6328 gmx_bcast(sizeof(r_mb),&r_mb,cr);
6330 /* We use an initial margin of 10% for the minimum cell size,
6331 * except when we are just below the non-bonded cut-off.
6333 if (Flags & MD_DDBONDCOMM)
6335 if (max(r_2b,r_mb) > comm->cutoff)
6337 r_bonded = max(r_2b,r_mb);
6338 r_bonded_limit = 1.1*r_bonded;
6339 comm->bBondComm = TRUE;
6344 r_bonded_limit = min(1.1*r_bonded,comm->cutoff);
6346 /* We determine cutoff_mbody later */
6350 /* No special bonded communication,
6351 * simply increase the DD cut-off.
6353 r_bonded_limit = 1.1*max(r_2b,r_mb);
6354 comm->cutoff_mbody = r_bonded_limit;
6355 comm->cutoff = max(comm->cutoff,comm->cutoff_mbody);
6358 comm->cellsize_limit = max(comm->cellsize_limit,r_bonded_limit);
6362 "Minimum cell size due to bonded interactions: %.3f nm\n",
6363 comm->cellsize_limit);
6367 if (dd->bInterCGcons && rconstr <= 0)
6369 /* There is a cell size limit due to the constraints (P-LINCS) */
6370 rconstr = constr_r_max(fplog,mtop,ir);
6374 "Estimated maximum distance required for P-LINCS: %.3f nm\n",
6376 if (rconstr > comm->cellsize_limit)
6378 fprintf(fplog,"This distance will limit the DD cell size, you can override this with -rcon\n");
6382 else if (rconstr > 0 && fplog)
6384 /* Here we do not check for dd->bInterCGcons,
6385 * because one can also set a cell size limit for virtual sites only
6386 * and at this point we don't know yet if there are intercg v-sites.
6389 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
6392 comm->cellsize_limit = max(comm->cellsize_limit,rconstr);
6394 comm->cgs_gl = gmx_mtop_global_cgs(mtop);
6398 copy_ivec(nc,dd->nc);
6399 set_dd_dim(fplog,dd);
6400 set_ddbox_cr(cr,&dd->nc,ir,box,&comm->cgs_gl,x,ddbox);
6402 if (cr->npmenodes == -1)
6406 acs = average_cellsize_min(dd,ddbox);
6407 if (acs < comm->cellsize_limit)
6411 fprintf(fplog,"ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n",acs,comm->cellsize_limit);
6413 gmx_fatal_collective(FARGS,cr,NULL,
6414 "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",
6415 acs,comm->cellsize_limit);
6420 set_ddbox_cr(cr,NULL,ir,box,&comm->cgs_gl,x,ddbox);
6422 /* We need to choose the optimal DD grid and possibly PME nodes */
6423 limit = dd_choose_grid(fplog,cr,dd,ir,mtop,box,ddbox,
6424 comm->eDLB!=edlbNO,dlb_scale,
6425 comm->cellsize_limit,comm->cutoff,
6426 comm->bInterCGBondeds,comm->bInterCGMultiBody);
6428 if (dd->nc[XX] == 0)
6430 bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
6431 sprintf(buf,"Change the number of nodes or mdrun option %s%s%s",
6432 !bC ? "-rdd" : "-rcon",
6433 comm->eDLB!=edlbNO ? " or -dds" : "",
6434 bC ? " or your LINCS settings" : "");
6436 gmx_fatal_collective(FARGS,cr,NULL,
6437 "There is no domain decomposition for %d nodes that is compatible with the given box and a minimum cell size of %g nm\n"
6439 "Look in the log file for details on the domain decomposition",
6440 cr->nnodes-cr->npmenodes,limit,buf);
6442 set_dd_dim(fplog,dd);
6448 "Domain decomposition grid %d x %d x %d, separate PME nodes %d\n",
6449 dd->nc[XX],dd->nc[YY],dd->nc[ZZ],cr->npmenodes);
6452 dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
6453 if (cr->nnodes - dd->nnodes != cr->npmenodes)
6455 gmx_fatal_collective(FARGS,cr,NULL,
6456 "The size of the domain decomposition grid (%d) does not match the number of nodes (%d). The total number of nodes is %d",
6457 dd->nnodes,cr->nnodes - cr->npmenodes,cr->nnodes);
6459 if (cr->npmenodes > dd->nnodes)
6461 gmx_fatal_collective(FARGS,cr,NULL,
6462 "The number of separate PME nodes (%d) is larger than the number of PP nodes (%d), this is not supported.",cr->npmenodes,dd->nnodes);
6464 if (cr->npmenodes > 0)
6466 comm->npmenodes = cr->npmenodes;
6470 comm->npmenodes = dd->nnodes;
6473 if (EEL_PME(ir->coulombtype))
6475 /* The following choices should match those
6476 * in comm_cost_est in domdec_setup.c.
6477 * Note that here the checks have to take into account
6478 * that the decomposition might occur in a different order than xyz
6479 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
6480 * in which case they will not match those in comm_cost_est,
6481 * but since that is mainly for testing purposes that's fine.
6483 if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
6484 comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
6485 getenv("GMX_PMEONEDD") == NULL)
6487 comm->npmedecompdim = 2;
6488 comm->npmenodes_x = dd->nc[XX];
6489 comm->npmenodes_y = comm->npmenodes/comm->npmenodes_x;
6493 /* In case nc is 1 in both x and y we could still choose to
6494 * decompose pme in y instead of x, but we use x for simplicity.
6496 comm->npmedecompdim = 1;
6497 if (dd->dim[0] == YY)
6499 comm->npmenodes_x = 1;
6500 comm->npmenodes_y = comm->npmenodes;
6504 comm->npmenodes_x = comm->npmenodes;
6505 comm->npmenodes_y = 1;
6510 fprintf(fplog,"PME domain decomposition: %d x %d x %d\n",
6511 comm->npmenodes_x,comm->npmenodes_y,1);
6516 comm->npmedecompdim = 0;
6517 comm->npmenodes_x = 0;
6518 comm->npmenodes_y = 0;
6521 /* Technically we don't need both of these,
6522 * but it simplifies code not having to recalculate it.
6524 *npme_x = comm->npmenodes_x;
6525 *npme_y = comm->npmenodes_y;
6527 snew(comm->slb_frac,DIM);
6528 if (comm->eDLB == edlbNO)
6530 comm->slb_frac[XX] = get_slb_frac(fplog,"x",dd->nc[XX],sizex);
6531 comm->slb_frac[YY] = get_slb_frac(fplog,"y",dd->nc[YY],sizey);
6532 comm->slb_frac[ZZ] = get_slb_frac(fplog,"z",dd->nc[ZZ],sizez);
6535 if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
6537 if (comm->bBondComm || comm->eDLB != edlbNO)
6539 /* Set the bonded communication distance to halfway
6540 * the minimum and the maximum,
6541 * since the extra communication cost is nearly zero.
6543 acs = average_cellsize_min(dd,ddbox);
6544 comm->cutoff_mbody = 0.5*(r_bonded + acs);
6545 if (comm->eDLB != edlbNO)
6547 /* Check if this does not limit the scaling */
6548 comm->cutoff_mbody = min(comm->cutoff_mbody,dlb_scale*acs);
6550 if (!comm->bBondComm)
6552 /* Without bBondComm do not go beyond the n.b. cut-off */
6553 comm->cutoff_mbody = min(comm->cutoff_mbody,comm->cutoff);
6554 if (comm->cellsize_limit >= comm->cutoff)
6556 /* We don't loose a lot of efficieny
6557 * when increasing it to the n.b. cut-off.
6558 * It can even be slightly faster, because we need
6559 * less checks for the communication setup.
6561 comm->cutoff_mbody = comm->cutoff;
6564 /* Check if we did not end up below our original limit */
6565 comm->cutoff_mbody = max(comm->cutoff_mbody,r_bonded_limit);
6567 if (comm->cutoff_mbody > comm->cellsize_limit)
6569 comm->cellsize_limit = comm->cutoff_mbody;
6572 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
6577 fprintf(debug,"Bonded atom communication beyond the cut-off: %d\n"
6578 "cellsize limit %f\n",
6579 comm->bBondComm,comm->cellsize_limit);
6584 check_dd_restrictions(cr,dd,ir,fplog);
6587 comm->globalcomm_step = INT_MIN;
6590 clear_dd_cycle_counts(dd);
6595 static void set_dlb_limits(gmx_domdec_t *dd)
6600 for(d=0; d<dd->ndim; d++)
6602 dd->comm->cd[d].np = dd->comm->cd[d].np_dlb;
6603 dd->comm->cellsize_min[dd->dim[d]] =
6604 dd->comm->cellsize_min_dlb[dd->dim[d]];
6609 static void turn_on_dlb(FILE *fplog,t_commrec *cr,gmx_large_int_t step)
6612 gmx_domdec_comm_t *comm;
6622 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);
6625 cellsize_min = comm->cellsize_min[dd->dim[0]];
6626 for(d=1; d<dd->ndim; d++)
6628 cellsize_min = min(cellsize_min,comm->cellsize_min[dd->dim[d]]);
6631 if (cellsize_min < comm->cellsize_limit*1.05)
6633 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");
6635 /* Change DLB from "auto" to "no". */
6636 comm->eDLB = edlbNO;
6641 dd_warning(cr,fplog,"NOTE: Turning on dynamic load balancing\n");
6642 comm->bDynLoadBal = TRUE;
6643 dd->bGridJump = TRUE;
6647 /* We can set the required cell size info here,
6648 * so we do not need to communicate this.
6649 * The grid is completely uniform.
6651 for(d=0; d<dd->ndim; d++)
6655 comm->load[d].sum_m = comm->load[d].sum;
6657 nc = dd->nc[dd->dim[d]];
6660 comm->root[d]->cell_f[i] = i/(real)nc;
6663 comm->root[d]->cell_f_max0[i] = i /(real)nc;
6664 comm->root[d]->cell_f_min1[i] = (i+1)/(real)nc;
6667 comm->root[d]->cell_f[nc] = 1.0;
6672 static char *init_bLocalCG(gmx_mtop_t *mtop)
6677 ncg = ncg_mtop(mtop);
6679 for(cg=0; cg<ncg; cg++)
6681 bLocalCG[cg] = FALSE;
6687 void dd_init_bondeds(FILE *fplog,
6688 gmx_domdec_t *dd,gmx_mtop_t *mtop,
6689 gmx_vsite_t *vsite,gmx_constr_t constr,
6690 t_inputrec *ir,gmx_bool bBCheck,cginfo_mb_t *cginfo_mb)
6692 gmx_domdec_comm_t *comm;
6696 dd_make_reverse_top(fplog,dd,mtop,vsite,constr,ir,bBCheck);
6700 if (comm->bBondComm)
6702 /* Communicate atoms beyond the cut-off for bonded interactions */
6705 comm->cglink = make_charge_group_links(mtop,dd,cginfo_mb);
6707 comm->bLocalCG = init_bLocalCG(mtop);
6711 /* Only communicate atoms based on cut-off */
6712 comm->cglink = NULL;
6713 comm->bLocalCG = NULL;
6717 static void print_dd_settings(FILE *fplog,gmx_domdec_t *dd,
6719 gmx_bool bDynLoadBal,real dlb_scale,
6722 gmx_domdec_comm_t *comm;
6737 fprintf(fplog,"The maximum number of communication pulses is:");
6738 for(d=0; d<dd->ndim; d++)
6740 fprintf(fplog," %c %d",dim2char(dd->dim[d]),comm->cd[d].np_dlb);
6742 fprintf(fplog,"\n");
6743 fprintf(fplog,"The minimum size for domain decomposition cells is %.3f nm\n",comm->cellsize_limit);
6744 fprintf(fplog,"The requested allowed shrink of DD cells (option -dds) is: %.2f\n",dlb_scale);
6745 fprintf(fplog,"The allowed shrink of domain decomposition cells is:");
6746 for(d=0; d<DIM; d++)
6750 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
6757 comm->cellsize_min_dlb[d]/
6758 (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
6760 fprintf(fplog," %c %.2f",dim2char(d),shrink);
6763 fprintf(fplog,"\n");
6767 set_dd_cell_sizes_slb(dd,ddbox,FALSE,np);
6768 fprintf(fplog,"The initial number of communication pulses is:");
6769 for(d=0; d<dd->ndim; d++)
6771 fprintf(fplog," %c %d",dim2char(dd->dim[d]),np[dd->dim[d]]);
6773 fprintf(fplog,"\n");
6774 fprintf(fplog,"The initial domain decomposition cell size is:");
6775 for(d=0; d<DIM; d++) {
6778 fprintf(fplog," %c %.2f nm",
6779 dim2char(d),dd->comm->cellsize_min[d]);
6782 fprintf(fplog,"\n\n");
6785 if (comm->bInterCGBondeds || dd->vsite_comm || dd->constraint_comm)
6787 fprintf(fplog,"The maximum allowed distance for charge groups involved in interactions is:\n");
6788 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6789 "non-bonded interactions","",comm->cutoff);
6793 limit = dd->comm->cellsize_limit;
6797 if (dynamic_dd_box(ddbox,ir))
6799 fprintf(fplog,"(the following are initial values, they could change due to box deformation)\n");
6801 limit = dd->comm->cellsize_min[XX];
6802 for(d=1; d<DIM; d++)
6804 limit = min(limit,dd->comm->cellsize_min[d]);
6808 if (comm->bInterCGBondeds)
6810 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6811 "two-body bonded interactions","(-rdd)",
6812 max(comm->cutoff,comm->cutoff_mbody));
6813 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6814 "multi-body bonded interactions","(-rdd)",
6815 (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff,limit));
6819 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6820 "virtual site constructions","(-rcon)",limit);
6822 if (dd->constraint_comm)
6824 sprintf(buf,"atoms separated by up to %d constraints",
6826 fprintf(fplog,"%40s %-7s %6.3f nm\n",
6827 buf,"(-rcon)",limit);
6829 fprintf(fplog,"\n");
6835 void set_dd_parameters(FILE *fplog,gmx_domdec_t *dd,real dlb_scale,
6836 t_inputrec *ir,t_forcerec *fr,
6839 gmx_domdec_comm_t *comm;
6840 int d,dim,npulse,npulse_d_max,npulse_d;
6847 bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
6849 if (EEL_PME(ir->coulombtype))
6851 init_ddpme(dd,&comm->ddpme[0],0);
6852 if (comm->npmedecompdim >= 2)
6854 init_ddpme(dd,&comm->ddpme[1],1);
6859 comm->npmenodes = 0;
6860 if (dd->pme_nodeid >= 0)
6862 gmx_fatal_collective(FARGS,NULL,dd,
6863 "Can not have separate PME nodes without PME electrostatics");
6867 /* If each molecule is a single charge group
6868 * or we use domain decomposition for each periodic dimension,
6869 * we do not need to take pbc into account for the bonded interactions.
6871 if (fr->ePBC == epbcNONE || !comm->bInterCGBondeds ||
6872 (dd->nc[XX]>1 && dd->nc[YY]>1 && (dd->nc[ZZ]>1 || fr->ePBC==epbcXY)))
6874 fr->bMolPBC = FALSE;
6883 fprintf(debug,"The DD cut-off is %f\n",comm->cutoff);
6885 if (comm->eDLB != edlbNO)
6887 /* Determine the maximum number of comm. pulses in one dimension */
6889 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
6891 /* Determine the maximum required number of grid pulses */
6892 if (comm->cellsize_limit >= comm->cutoff)
6894 /* Only a single pulse is required */
6897 else if (!bNoCutOff && comm->cellsize_limit > 0)
6899 /* We round down slightly here to avoid overhead due to the latency
6900 * of extra communication calls when the cut-off
6901 * would be only slightly longer than the cell size.
6902 * Later cellsize_limit is redetermined,
6903 * so we can not miss interactions due to this rounding.
6905 npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
6909 /* There is no cell size limit */
6910 npulse = max(dd->nc[XX]-1,max(dd->nc[YY]-1,dd->nc[ZZ]-1));
6913 if (!bNoCutOff && npulse > 1)
6915 /* See if we can do with less pulses, based on dlb_scale */
6917 for(d=0; d<dd->ndim; d++)
6920 npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
6921 /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
6922 npulse_d_max = max(npulse_d_max,npulse_d);
6924 npulse = min(npulse,npulse_d_max);
6927 /* This env var can override npulse */
6928 d = dd_nst_env(fplog,"GMX_DD_NPULSE",0);
6935 comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
6936 for(d=0; d<dd->ndim; d++)
6938 comm->cd[d].np_dlb = min(npulse,dd->nc[dd->dim[d]]-1);
6939 comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
6940 snew(comm->cd[d].ind,comm->cd[d].np_nalloc);
6941 comm->maxpulse = max(comm->maxpulse,comm->cd[d].np_dlb);
6942 if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
6944 comm->bVacDLBNoLimit = FALSE;
6948 /* cellsize_limit is set for LINCS in init_domain_decomposition */
6949 if (!comm->bVacDLBNoLimit)
6951 comm->cellsize_limit = max(comm->cellsize_limit,
6952 comm->cutoff/comm->maxpulse);
6954 comm->cellsize_limit = max(comm->cellsize_limit,comm->cutoff_mbody);
6955 /* Set the minimum cell size for each DD dimension */
6956 for(d=0; d<dd->ndim; d++)
6958 if (comm->bVacDLBNoLimit ||
6959 comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
6961 comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
6965 comm->cellsize_min_dlb[dd->dim[d]] =
6966 comm->cutoff/comm->cd[d].np_dlb;
6969 if (comm->cutoff_mbody <= 0)
6971 comm->cutoff_mbody = min(comm->cutoff,comm->cellsize_limit);
6973 if (comm->bDynLoadBal)
6979 print_dd_settings(fplog,dd,ir,comm->bDynLoadBal,dlb_scale,ddbox);
6980 if (comm->eDLB == edlbAUTO)
6984 fprintf(fplog,"When dynamic load balancing gets turned on, these settings will change to:\n");
6986 print_dd_settings(fplog,dd,ir,TRUE,dlb_scale,ddbox);
6989 if (ir->ePBC == epbcNONE)
6991 vol_frac = 1 - 1/(double)dd->nnodes;
6996 (1 + comm_box_frac(dd->nc,comm->cutoff,ddbox))/(double)dd->nnodes;
7000 fprintf(debug,"Volume fraction for all DD zones: %f\n",vol_frac);
7002 natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
7004 dd->ga2la = ga2la_init(natoms_tot,vol_frac*natoms_tot);
7007 static void merge_cg_buffers(int ncell,
7008 gmx_domdec_comm_dim_t *cd, int pulse,
7010 int *index_gl, int *recv_i,
7011 rvec *cg_cm, rvec *recv_vr,
7013 cginfo_mb_t *cginfo_mb,int *cginfo)
7015 gmx_domdec_ind_t *ind,*ind_p;
7016 int p,cell,c,cg,cg0,cg1,cg_gl,nat;
7019 ind = &cd->ind[pulse];
7021 /* First correct the already stored data */
7022 shift = ind->nrecv[ncell];
7023 for(cell=ncell-1; cell>=0; cell--)
7025 shift -= ind->nrecv[cell];
7028 /* Move the cg's present from previous grid pulses */
7029 cg0 = ncg_cell[ncell+cell];
7030 cg1 = ncg_cell[ncell+cell+1];
7031 cgindex[cg1+shift] = cgindex[cg1];
7032 for(cg=cg1-1; cg>=cg0; cg--)
7034 index_gl[cg+shift] = index_gl[cg];
7035 copy_rvec(cg_cm[cg],cg_cm[cg+shift]);
7036 cgindex[cg+shift] = cgindex[cg];
7037 cginfo[cg+shift] = cginfo[cg];
7039 /* Correct the already stored send indices for the shift */
7040 for(p=1; p<=pulse; p++)
7042 ind_p = &cd->ind[p];
7044 for(c=0; c<cell; c++)
7046 cg0 += ind_p->nsend[c];
7048 cg1 = cg0 + ind_p->nsend[cell];
7049 for(cg=cg0; cg<cg1; cg++)
7051 ind_p->index[cg] += shift;
7057 /* Merge in the communicated buffers */
7061 for(cell=0; cell<ncell; cell++)
7063 cg1 = ncg_cell[ncell+cell+1] + shift;
7066 /* Correct the old cg indices */
7067 for(cg=ncg_cell[ncell+cell]; cg<cg1; cg++)
7069 cgindex[cg+1] += shift_at;
7072 for(cg=0; cg<ind->nrecv[cell]; cg++)
7074 /* Copy this charge group from the buffer */
7075 index_gl[cg1] = recv_i[cg0];
7076 copy_rvec(recv_vr[cg0],cg_cm[cg1]);
7077 /* Add it to the cgindex */
7078 cg_gl = index_gl[cg1];
7079 cginfo[cg1] = ddcginfo(cginfo_mb,cg_gl);
7080 nat = GET_CGINFO_NATOMS(cginfo[cg1]);
7081 cgindex[cg1+1] = cgindex[cg1] + nat;
7086 shift += ind->nrecv[cell];
7087 ncg_cell[ncell+cell+1] = cg1;
7091 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
7092 int nzone,int cg0,const int *cgindex)
7096 /* Store the atom block boundaries for easy copying of communication buffers
7099 for(zone=0; zone<nzone; zone++)
7101 for(p=0; p<cd->np; p++) {
7102 cd->ind[p].cell2at0[zone] = cgindex[cg];
7103 cg += cd->ind[p].nrecv[zone];
7104 cd->ind[p].cell2at1[zone] = cgindex[cg];
7109 static gmx_bool missing_link(t_blocka *link,int cg_gl,char *bLocalCG)
7115 for(i=link->index[cg_gl]; i<link->index[cg_gl+1]; i++)
7117 if (!bLocalCG[link->a[i]])
7126 static void setup_dd_communication(gmx_domdec_t *dd,
7127 matrix box,gmx_ddbox_t *ddbox,t_forcerec *fr)
7129 int dim_ind,dim,dim0,dim1=-1,dim2=-1,dimd,p,nat_tot;
7130 int nzone,nzone_send,zone,zonei,cg0,cg1;
7131 int c,i,j,cg,cg_gl,nrcg;
7132 int *zone_cg_range,pos_cg,*index_gl,*cgindex,*recv_i;
7133 gmx_domdec_comm_t *comm;
7134 gmx_domdec_zones_t *zones;
7135 gmx_domdec_comm_dim_t *cd;
7136 gmx_domdec_ind_t *ind;
7137 cginfo_mb_t *cginfo_mb;
7138 gmx_bool bBondComm,bDist2B,bDistMB,bDistMB_pulse,bDistBonded,bScrew;
7139 real r_mb,r_comm2,r_scomm2,r_bcomm2,r,r_0,r_1,r2,rb2,r2inc,inv_ncg,tric_sh;
7141 real corner[DIM][4],corner_round_0=0,corner_round_1[4];
7142 real bcorner[DIM],bcorner_round_1=0;
7144 rvec *cg_cm,*normal,*v_d,*v_0=NULL,*v_1=NULL,*recv_vr;
7145 real skew_fac2_d,skew_fac_01;
7151 fprintf(debug,"Setting up DD communication\n");
7157 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7159 dim = dd->dim[dim_ind];
7161 /* Check if we need to use triclinic distances */
7162 tric_dist[dim_ind] = 0;
7163 for(i=0; i<=dim_ind; i++)
7165 if (ddbox->tric_dir[dd->dim[i]])
7167 tric_dist[dim_ind] = 1;
7172 bBondComm = comm->bBondComm;
7174 /* Do we need to determine extra distances for multi-body bondeds? */
7175 bDistMB = (comm->bInterCGMultiBody && dd->bGridJump && dd->ndim > 1);
7177 /* Do we need to determine extra distances for only two-body bondeds? */
7178 bDist2B = (bBondComm && !bDistMB);
7180 r_comm2 = sqr(comm->cutoff);
7181 r_bcomm2 = sqr(comm->cutoff_mbody);
7185 fprintf(debug,"bBondComm %d, r_bc %f\n",bBondComm,sqrt(r_bcomm2));
7188 zones = &comm->zones;
7191 /* The first dimension is equal for all cells */
7192 corner[0][0] = comm->cell_x0[dim0];
7195 bcorner[0] = corner[0][0];
7200 /* This cell row is only seen from the first row */
7201 corner[1][0] = comm->cell_x0[dim1];
7202 /* All rows can see this row */
7203 corner[1][1] = comm->cell_x0[dim1];
7206 corner[1][1] = max(comm->cell_x0[dim1],comm->zone_d1[1].mch0);
7209 /* For the multi-body distance we need the maximum */
7210 bcorner[1] = max(comm->cell_x0[dim1],comm->zone_d1[1].p1_0);
7213 /* Set the upper-right corner for rounding */
7214 corner_round_0 = comm->cell_x1[dim0];
7221 corner[2][j] = comm->cell_x0[dim2];
7225 /* Use the maximum of the i-cells that see a j-cell */
7226 for(i=0; i<zones->nizone; i++)
7228 for(j=zones->izone[i].j0; j<zones->izone[i].j1; j++)
7234 comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
7240 /* For the multi-body distance we need the maximum */
7241 bcorner[2] = comm->cell_x0[dim2];
7246 bcorner[2] = max(bcorner[2],
7247 comm->zone_d2[i][j].p1_0);
7253 /* Set the upper-right corner for rounding */
7254 /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
7255 * Only cell (0,0,0) can see cell 7 (1,1,1)
7257 corner_round_1[0] = comm->cell_x1[dim1];
7258 corner_round_1[3] = comm->cell_x1[dim1];
7261 corner_round_1[0] = max(comm->cell_x1[dim1],
7262 comm->zone_d1[1].mch1);
7265 /* For the multi-body distance we need the maximum */
7266 bcorner_round_1 = max(comm->cell_x1[dim1],
7267 comm->zone_d1[1].p1_1);
7273 /* Triclinic stuff */
7274 normal = ddbox->normal;
7278 v_0 = ddbox->v[dim0];
7279 if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
7281 /* Determine the coupling coefficient for the distances
7282 * to the cell planes along dim0 and dim1 through dim2.
7283 * This is required for correct rounding.
7286 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
7289 fprintf(debug,"\nskew_fac_01 %f\n",skew_fac_01);
7295 v_1 = ddbox->v[dim1];
7298 zone_cg_range = zones->cg_range;
7299 index_gl = dd->index_gl;
7300 cgindex = dd->cgindex;
7301 cginfo_mb = fr->cginfo_mb;
7303 zone_cg_range[0] = 0;
7304 zone_cg_range[1] = dd->ncg_home;
7305 comm->zone_ncg1[0] = dd->ncg_home;
7306 pos_cg = dd->ncg_home;
7308 nat_tot = dd->nat_home;
7310 for(dim_ind=0; dim_ind<dd->ndim; dim_ind++)
7312 dim = dd->dim[dim_ind];
7313 cd = &comm->cd[dim_ind];
7315 if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
7317 /* No pbc in this dimension, the first node should not comm. */
7325 bScrew = (dd->bScrewPBC && dim == XX);
7327 v_d = ddbox->v[dim];
7328 skew_fac2_d = sqr(ddbox->skew_fac[dim]);
7330 cd->bInPlace = TRUE;
7331 for(p=0; p<cd->np; p++)
7333 /* Only atoms communicated in the first pulse are used
7334 * for multi-body bonded interactions or for bBondComm.
7336 bDistBonded = ((bDistMB || bDist2B) && p == 0);
7337 bDistMB_pulse = (bDistMB && bDistBonded);
7342 for(zone=0; zone<nzone_send; zone++)
7344 if (tric_dist[dim_ind] && dim_ind > 0)
7346 /* Determine slightly more optimized skew_fac's
7348 * This reduces the number of communicated atoms
7349 * by about 10% for 3D DD of rhombic dodecahedra.
7351 for(dimd=0; dimd<dim; dimd++)
7353 sf2_round[dimd] = 1;
7354 if (ddbox->tric_dir[dimd])
7356 for(i=dd->dim[dimd]+1; i<DIM; i++)
7358 /* If we are shifted in dimension i
7359 * and the cell plane is tilted forward
7360 * in dimension i, skip this coupling.
7362 if (!(zones->shift[nzone+zone][i] &&
7363 ddbox->v[dimd][i][dimd] >= 0))
7366 sqr(ddbox->v[dimd][i][dimd]);
7369 sf2_round[dimd] = 1/sf2_round[dimd];
7374 zonei = zone_perm[dim_ind][zone];
7377 /* Here we permutate the zones to obtain a convenient order
7378 * for neighbor searching
7380 cg0 = zone_cg_range[zonei];
7381 cg1 = zone_cg_range[zonei+1];
7385 /* Look only at the cg's received in the previous grid pulse
7387 cg1 = zone_cg_range[nzone+zone+1];
7388 cg0 = cg1 - cd->ind[p-1].nrecv[zone];
7390 ind->nsend[zone] = 0;
7391 for(cg=cg0; cg<cg1; cg++)
7395 if (tric_dist[dim_ind] == 0)
7397 /* Rectangular direction, easy */
7398 r = cg_cm[cg][dim] - corner[dim_ind][zone];
7405 r = cg_cm[cg][dim] - bcorner[dim_ind];
7411 /* Rounding gives at most a 16% reduction
7412 * in communicated atoms
7414 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7416 r = cg_cm[cg][dim0] - corner_round_0;
7417 /* This is the first dimension, so always r >= 0 */
7424 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7426 r = cg_cm[cg][dim1] - corner_round_1[zone];
7433 r = cg_cm[cg][dim1] - bcorner_round_1;
7443 /* Triclinic direction, more complicated */
7446 /* Rounding, conservative as the skew_fac multiplication
7447 * will slightly underestimate the distance.
7449 if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
7451 rn[dim0] = cg_cm[cg][dim0] - corner_round_0;
7452 for(i=dim0+1; i<DIM; i++)
7454 rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
7456 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
7459 rb[dim0] = rn[dim0];
7462 /* Take care that the cell planes along dim0 might not
7463 * be orthogonal to those along dim1 and dim2.
7465 for(i=1; i<=dim_ind; i++)
7468 if (normal[dim0][dimd] > 0)
7470 rn[dimd] -= rn[dim0]*normal[dim0][dimd];
7473 rb[dimd] -= rb[dim0]*normal[dim0][dimd];
7478 if (dim_ind == 2 && (zonei == 2 || zonei == 3))
7480 rn[dim1] += cg_cm[cg][dim1] - corner_round_1[zone];
7482 for(i=dim1+1; i<DIM; i++)
7484 tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
7486 rn[dim1] += tric_sh;
7489 r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
7490 /* Take care of coupling of the distances
7491 * to the planes along dim0 and dim1 through dim2.
7493 r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
7494 /* Take care that the cell planes along dim1
7495 * might not be orthogonal to that along dim2.
7497 if (normal[dim1][dim2] > 0)
7499 rn[dim2] -= rn[dim1]*normal[dim1][dim2];
7505 cg_cm[cg][dim1] - bcorner_round_1 + tric_sh;
7508 rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
7509 /* Take care of coupling of the distances
7510 * to the planes along dim0 and dim1 through dim2.
7512 rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
7513 /* Take care that the cell planes along dim1
7514 * might not be orthogonal to that along dim2.
7516 if (normal[dim1][dim2] > 0)
7518 rb[dim2] -= rb[dim1]*normal[dim1][dim2];
7523 /* The distance along the communication direction */
7524 rn[dim] += cg_cm[cg][dim] - corner[dim_ind][zone];
7526 for(i=dim+1; i<DIM; i++)
7528 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
7533 r2 += rn[dim]*rn[dim]*skew_fac2_d;
7534 /* Take care of coupling of the distances
7535 * to the planes along dim0 and dim1 through dim2.
7537 if (dim_ind == 1 && zonei == 1)
7539 r2 -= rn[dim0]*rn[dim]*skew_fac_01;
7545 rb[dim] += cg_cm[cg][dim] - bcorner[dim_ind] + tric_sh;
7548 rb2 += rb[dim]*rb[dim]*skew_fac2_d;
7549 /* Take care of coupling of the distances
7550 * to the planes along dim0 and dim1 through dim2.
7552 if (dim_ind == 1 && zonei == 1)
7554 rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
7562 ((bDistMB && rb2 < r_bcomm2) ||
7563 (bDist2B && r2 < r_bcomm2)) &&
7565 (GET_CGINFO_BOND_INTER(fr->cginfo[cg]) &&
7566 missing_link(comm->cglink,index_gl[cg],
7569 /* Make an index to the local charge groups */
7570 if (nsend+1 > ind->nalloc)
7572 ind->nalloc = over_alloc_large(nsend+1);
7573 srenew(ind->index,ind->nalloc);
7575 if (nsend+1 > comm->nalloc_int)
7577 comm->nalloc_int = over_alloc_large(nsend+1);
7578 srenew(comm->buf_int,comm->nalloc_int);
7580 ind->index[nsend] = cg;
7581 comm->buf_int[nsend] = index_gl[cg];
7583 vec_rvec_check_alloc(&comm->vbuf,nsend+1);
7585 if (dd->ci[dim] == 0)
7587 /* Correct cg_cm for pbc */
7588 rvec_add(cg_cm[cg],box[dim],comm->vbuf.v[nsend]);
7591 comm->vbuf.v[nsend][YY] =
7592 box[YY][YY]-comm->vbuf.v[nsend][YY];
7593 comm->vbuf.v[nsend][ZZ] =
7594 box[ZZ][ZZ]-comm->vbuf.v[nsend][ZZ];
7599 copy_rvec(cg_cm[cg],comm->vbuf.v[nsend]);
7602 nat += cgindex[cg+1] - cgindex[cg];
7606 /* Clear the counts in case we do not have pbc */
7607 for(zone=nzone_send; zone<nzone; zone++)
7609 ind->nsend[zone] = 0;
7611 ind->nsend[nzone] = nsend;
7612 ind->nsend[nzone+1] = nat;
7613 /* Communicate the number of cg's and atoms to receive */
7614 dd_sendrecv_int(dd, dim_ind, dddirBackward,
7615 ind->nsend, nzone+2,
7616 ind->nrecv, nzone+2);
7618 /* The rvec buffer is also required for atom buffers of size nsend
7619 * in dd_move_x and dd_move_f.
7621 vec_rvec_check_alloc(&comm->vbuf,ind->nsend[nzone+1]);
7625 /* We can receive in place if only the last zone is not empty */
7626 for(zone=0; zone<nzone-1; zone++)
7628 if (ind->nrecv[zone] > 0)
7630 cd->bInPlace = FALSE;
7635 /* The int buffer is only required here for the cg indices */
7636 if (ind->nrecv[nzone] > comm->nalloc_int2)
7638 comm->nalloc_int2 = over_alloc_dd(ind->nrecv[nzone]);
7639 srenew(comm->buf_int2,comm->nalloc_int2);
7641 /* The rvec buffer is also required for atom buffers
7642 * of size nrecv in dd_move_x and dd_move_f.
7644 i = max(cd->ind[0].nrecv[nzone+1],ind->nrecv[nzone+1]);
7645 vec_rvec_check_alloc(&comm->vbuf2,i);
7649 /* Make space for the global cg indices */
7650 if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
7651 || dd->cg_nalloc == 0)
7653 dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
7654 srenew(index_gl,dd->cg_nalloc);
7655 srenew(cgindex,dd->cg_nalloc+1);
7657 /* Communicate the global cg indices */
7660 recv_i = index_gl + pos_cg;
7664 recv_i = comm->buf_int2;
7666 dd_sendrecv_int(dd, dim_ind, dddirBackward,
7667 comm->buf_int, nsend,
7668 recv_i, ind->nrecv[nzone]);
7670 /* Make space for cg_cm */
7671 if (pos_cg + ind->nrecv[nzone] > fr->cg_nalloc)
7673 dd_realloc_fr_cg(fr,pos_cg + ind->nrecv[nzone]);
7676 /* Communicate cg_cm */
7679 recv_vr = cg_cm + pos_cg;
7683 recv_vr = comm->vbuf2.v;
7685 dd_sendrecv_rvec(dd, dim_ind, dddirBackward,
7686 comm->vbuf.v, nsend,
7687 recv_vr, ind->nrecv[nzone]);
7689 /* Make the charge group index */
7692 zone = (p == 0 ? 0 : nzone - 1);
7693 while (zone < nzone)
7695 for(cg=0; cg<ind->nrecv[zone]; cg++)
7697 cg_gl = index_gl[pos_cg];
7698 fr->cginfo[pos_cg] = ddcginfo(cginfo_mb,cg_gl);
7699 nrcg = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
7700 cgindex[pos_cg+1] = cgindex[pos_cg] + nrcg;
7703 /* Update the charge group presence,
7704 * so we can use it in the next pass of the loop.
7706 comm->bLocalCG[cg_gl] = TRUE;
7712 comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
7715 zone_cg_range[nzone+zone] = pos_cg;
7720 /* This part of the code is never executed with bBondComm. */
7721 merge_cg_buffers(nzone,cd,p,zone_cg_range,
7722 index_gl,recv_i,cg_cm,recv_vr,
7723 cgindex,fr->cginfo_mb,fr->cginfo);
7724 pos_cg += ind->nrecv[nzone];
7726 nat_tot += ind->nrecv[nzone+1];
7730 /* Store the atom block for easy copying of communication buffers */
7731 make_cell2at_index(cd,nzone,zone_cg_range[nzone],cgindex);
7735 dd->index_gl = index_gl;
7736 dd->cgindex = cgindex;
7738 dd->ncg_tot = zone_cg_range[zones->n];
7739 dd->nat_tot = nat_tot;
7740 comm->nat[ddnatHOME] = dd->nat_home;
7741 for(i=ddnatZONE; i<ddnatNR; i++)
7743 comm->nat[i] = dd->nat_tot;
7748 /* We don't need to update cginfo, since that was alrady done above.
7749 * So we pass NULL for the forcerec.
7751 dd_set_cginfo(dd->index_gl,dd->ncg_home,dd->ncg_tot,
7752 NULL,comm->bLocalCG);
7757 fprintf(debug,"Finished setting up DD communication, zones:");
7758 for(c=0; c<zones->n; c++)
7760 fprintf(debug," %d",zones->cg_range[c+1]-zones->cg_range[c]);
7762 fprintf(debug,"\n");
7766 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
7770 for(c=0; c<zones->nizone; c++)
7772 zones->izone[c].cg1 = zones->cg_range[c+1];
7773 zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
7774 zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
7778 static int comp_cgsort(const void *a,const void *b)
7782 gmx_cgsort_t *cga,*cgb;
7783 cga = (gmx_cgsort_t *)a;
7784 cgb = (gmx_cgsort_t *)b;
7786 comp = cga->nsc - cgb->nsc;
7789 comp = cga->ind_gl - cgb->ind_gl;
7795 static void order_int_cg(int n,gmx_cgsort_t *sort,
7800 /* Order the data */
7803 buf[i] = a[sort[i].ind];
7806 /* Copy back to the original array */
7813 static void order_vec_cg(int n,gmx_cgsort_t *sort,
7818 /* Order the data */
7821 copy_rvec(v[sort[i].ind],buf[i]);
7824 /* Copy back to the original array */
7827 copy_rvec(buf[i],v[i]);
7831 static void order_vec_atom(int ncg,int *cgindex,gmx_cgsort_t *sort,
7834 int a,atot,cg,cg0,cg1,i;
7836 /* Order the data */
7838 for(cg=0; cg<ncg; cg++)
7840 cg0 = cgindex[sort[cg].ind];
7841 cg1 = cgindex[sort[cg].ind+1];
7842 for(i=cg0; i<cg1; i++)
7844 copy_rvec(v[i],buf[a]);
7850 /* Copy back to the original array */
7851 for(a=0; a<atot; a++)
7853 copy_rvec(buf[a],v[a]);
7857 static void ordered_sort(int nsort2,gmx_cgsort_t *sort2,
7858 int nsort_new,gmx_cgsort_t *sort_new,
7859 gmx_cgsort_t *sort1)
7863 /* The new indices are not very ordered, so we qsort them */
7864 qsort_threadsafe(sort_new,nsort_new,sizeof(sort_new[0]),comp_cgsort);
7866 /* sort2 is already ordered, so now we can merge the two arrays */
7870 while(i2 < nsort2 || i_new < nsort_new)
7874 sort1[i1++] = sort_new[i_new++];
7876 else if (i_new == nsort_new)
7878 sort1[i1++] = sort2[i2++];
7880 else if (sort2[i2].nsc < sort_new[i_new].nsc ||
7881 (sort2[i2].nsc == sort_new[i_new].nsc &&
7882 sort2[i2].ind_gl < sort_new[i_new].ind_gl))
7884 sort1[i1++] = sort2[i2++];
7888 sort1[i1++] = sort_new[i_new++];
7893 static void dd_sort_state(gmx_domdec_t *dd,int ePBC,
7894 rvec *cgcm,t_forcerec *fr,t_state *state,
7897 gmx_domdec_sort_t *sort;
7898 gmx_cgsort_t *cgsort,*sort_i;
7899 int ncg_new,nsort2,nsort_new,i,cell_index,*ibuf,cgsize;
7902 sort = dd->comm->sort;
7904 if (dd->ncg_home > sort->sort_nalloc)
7906 sort->sort_nalloc = over_alloc_dd(dd->ncg_home);
7907 srenew(sort->sort1,sort->sort_nalloc);
7908 srenew(sort->sort2,sort->sort_nalloc);
7911 if (ncg_home_old >= 0)
7913 /* The charge groups that remained in the same ns grid cell
7914 * are completely ordered. So we can sort efficiently by sorting
7915 * the charge groups that did move into the stationary list.
7920 for(i=0; i<dd->ncg_home; i++)
7922 /* Check if this cg did not move to another node */
7923 cell_index = fr->ns.grid->cell_index[i];
7924 if (cell_index != 4*fr->ns.grid->ncells)
7926 if (i >= ncg_home_old || cell_index != sort->sort1[i].nsc)
7928 /* This cg is new on this node or moved ns grid cell */
7929 if (nsort_new >= sort->sort_new_nalloc)
7931 sort->sort_new_nalloc = over_alloc_dd(nsort_new+1);
7932 srenew(sort->sort_new,sort->sort_new_nalloc);
7934 sort_i = &(sort->sort_new[nsort_new++]);
7938 /* This cg did not move */
7939 sort_i = &(sort->sort2[nsort2++]);
7941 /* Sort on the ns grid cell indices
7942 * and the global topology index
7944 sort_i->nsc = cell_index;
7945 sort_i->ind_gl = dd->index_gl[i];
7952 fprintf(debug,"ordered sort cgs: stationary %d moved %d\n",
7955 /* Sort efficiently */
7956 ordered_sort(nsort2,sort->sort2,nsort_new,sort->sort_new,sort->sort1);
7960 cgsort = sort->sort1;
7962 for(i=0; i<dd->ncg_home; i++)
7964 /* Sort on the ns grid cell indices
7965 * and the global topology index
7967 cgsort[i].nsc = fr->ns.grid->cell_index[i];
7968 cgsort[i].ind_gl = dd->index_gl[i];
7970 if (cgsort[i].nsc != 4*fr->ns.grid->ncells)
7977 fprintf(debug,"qsort cgs: %d new home %d\n",dd->ncg_home,ncg_new);
7979 /* Determine the order of the charge groups using qsort */
7980 qsort_threadsafe(cgsort,dd->ncg_home,sizeof(cgsort[0]),comp_cgsort);
7982 cgsort = sort->sort1;
7984 /* We alloc with the old size, since cgindex is still old */
7985 vec_rvec_check_alloc(&dd->comm->vbuf,dd->cgindex[dd->ncg_home]);
7986 vbuf = dd->comm->vbuf.v;
7988 /* Remove the charge groups which are no longer at home here */
7989 dd->ncg_home = ncg_new;
7991 /* Reorder the state */
7992 for(i=0; i<estNR; i++)
7994 if (EST_DISTR(i) && (state->flags & (1<<i)))
7999 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->x,vbuf);
8002 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->v,vbuf);
8005 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->sd_X,vbuf);
8008 order_vec_atom(dd->ncg_home,dd->cgindex,cgsort,state->cg_p,vbuf);
8012 case estDISRE_INITF:
8013 case estDISRE_RM3TAV:
8014 case estORIRE_INITF:
8016 /* No ordering required */
8019 gmx_incons("Unknown state entry encountered in dd_sort_state");
8025 order_vec_cg(dd->ncg_home,cgsort,cgcm,vbuf);
8027 if (dd->ncg_home+1 > sort->ibuf_nalloc)
8029 sort->ibuf_nalloc = over_alloc_dd(dd->ncg_home+1);
8030 srenew(sort->ibuf,sort->ibuf_nalloc);
8033 /* Reorder the global cg index */
8034 order_int_cg(dd->ncg_home,cgsort,dd->index_gl,ibuf);
8035 /* Reorder the cginfo */
8036 order_int_cg(dd->ncg_home,cgsort,fr->cginfo,ibuf);
8037 /* Rebuild the local cg index */
8039 for(i=0; i<dd->ncg_home; i++)
8041 cgsize = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
8042 ibuf[i+1] = ibuf[i] + cgsize;
8044 for(i=0; i<dd->ncg_home+1; i++)
8046 dd->cgindex[i] = ibuf[i];
8048 /* Set the home atom number */
8049 dd->nat_home = dd->cgindex[dd->ncg_home];
8051 /* Copy the sorted ns cell indices back to the ns grid struct */
8052 for(i=0; i<dd->ncg_home; i++)
8054 fr->ns.grid->cell_index[i] = cgsort[i].nsc;
8056 fr->ns.grid->nr = dd->ncg_home;
8059 static void add_dd_statistics(gmx_domdec_t *dd)
8061 gmx_domdec_comm_t *comm;
8066 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8068 comm->sum_nat[ddnat-ddnatZONE] +=
8069 comm->nat[ddnat] - comm->nat[ddnat-1];
8074 void reset_dd_statistics_counters(gmx_domdec_t *dd)
8076 gmx_domdec_comm_t *comm;
8081 /* Reset all the statistics and counters for total run counting */
8082 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8084 comm->sum_nat[ddnat-ddnatZONE] = 0;
8088 comm->load_step = 0;
8091 clear_ivec(comm->load_lim);
8096 void print_dd_statistics(t_commrec *cr,t_inputrec *ir,FILE *fplog)
8098 gmx_domdec_comm_t *comm;
8102 comm = cr->dd->comm;
8104 gmx_sumd(ddnatNR-ddnatZONE,comm->sum_nat,cr);
8111 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");
8113 for(ddnat=ddnatZONE; ddnat<ddnatNR; ddnat++)
8115 av = comm->sum_nat[ddnat-ddnatZONE]/comm->ndecomp;
8120 " av. #atoms communicated per step for force: %d x %.1f\n",
8124 if (cr->dd->vsite_comm)
8127 " av. #atoms communicated per step for vsites: %d x %.1f\n",
8128 (EEL_PME(ir->coulombtype) || ir->coulombtype==eelEWALD) ? 3 : 2,
8133 if (cr->dd->constraint_comm)
8136 " av. #atoms communicated per step for LINCS: %d x %.1f\n",
8137 1 + ir->nLincsIter,av);
8141 gmx_incons(" Unknown type for DD statistics");
8144 fprintf(fplog,"\n");
8146 if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
8148 print_dd_load_av(fplog,cr->dd);
8152 void dd_partition_system(FILE *fplog,
8153 gmx_large_int_t step,
8155 gmx_bool bMasterState,
8157 t_state *state_global,
8158 gmx_mtop_t *top_global,
8160 t_state *state_local,
8163 gmx_localtop_t *top_local,
8166 gmx_shellfc_t shellfc,
8167 gmx_constr_t constr,
8169 gmx_wallcycle_t wcycle,
8173 gmx_domdec_comm_t *comm;
8174 gmx_ddbox_t ddbox={0};
8176 gmx_large_int_t step_pcoupl;
8177 rvec cell_ns_x0,cell_ns_x1;
8178 int i,j,n,cg0=0,ncg_home_old=-1,nat_f_novirsum;
8179 gmx_bool bBoxChanged,bNStGlobalComm,bDoDLB,bCheckDLB,bTurnOnDLB,bLogLoad;
8180 gmx_bool bRedist,bSortCG,bResortAll;
8188 bBoxChanged = (bMasterState || DEFORM(*ir));
8189 if (ir->epc != epcNO)
8191 /* With nstpcouple > 1 pressure coupling happens.
8192 * one step after calculating the pressure.
8193 * Box scaling happens at the end of the MD step,
8194 * after the DD partitioning.
8195 * We therefore have to do DLB in the first partitioning
8196 * after an MD step where P-coupling occured.
8197 * We need to determine the last step in which p-coupling occurred.
8198 * MRS -- need to validate this for vv?
8203 step_pcoupl = step - 1;
8207 step_pcoupl = ((step - 1)/n)*n + 1;
8209 if (step_pcoupl >= comm->globalcomm_step)
8215 bNStGlobalComm = (step >= comm->globalcomm_step + nstglobalcomm);
8217 if (!comm->bDynLoadBal)
8223 /* Should we do dynamic load balacing this step?
8224 * Since it requires (possibly expensive) global communication,
8225 * we might want to do DLB less frequently.
8227 if (bBoxChanged || ir->epc != epcNO)
8229 bDoDLB = bBoxChanged;
8233 bDoDLB = bNStGlobalComm;
8237 /* Check if we have recorded loads on the nodes */
8238 if (comm->bRecordLoad && dd_load_count(comm))
8240 if (comm->eDLB == edlbAUTO && !comm->bDynLoadBal)
8242 /* Check if we should use DLB at the second partitioning
8243 * and every 100 partitionings,
8244 * so the extra communication cost is negligible.
8246 n = max(100,nstglobalcomm);
8247 bCheckDLB = (comm->n_load_collect == 0 ||
8248 comm->n_load_have % n == n-1);
8255 /* Print load every nstlog, first and last step to the log file */
8256 bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
8257 comm->n_load_collect == 0 ||
8259 (step + ir->nstlist > ir->init_step + ir->nsteps)));
8261 /* Avoid extra communication due to verbose screen output
8262 * when nstglobalcomm is set.
8264 if (bDoDLB || bLogLoad || bCheckDLB ||
8265 (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
8267 get_load_distribution(dd,wcycle);
8272 dd_print_load(fplog,dd,step-1);
8276 dd_print_load_verbose(dd);
8279 comm->n_load_collect++;
8282 /* Since the timings are node dependent, the master decides */
8286 (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
8289 fprintf(debug,"step %s, imb loss %f\n",
8290 gmx_step_str(step,sbuf),
8291 dd_force_imb_perf_loss(dd));
8294 dd_bcast(dd,sizeof(bTurnOnDLB),&bTurnOnDLB);
8297 turn_on_dlb(fplog,cr,step);
8302 comm->n_load_have++;
8305 cgs_gl = &comm->cgs_gl;
8310 /* Clear the old state */
8311 clear_dd_indices(dd,0,0);
8313 set_ddbox(dd,bMasterState,cr,ir,state_global->box,
8314 TRUE,cgs_gl,state_global->x,&ddbox);
8316 get_cg_distribution(fplog,step,dd,cgs_gl,
8317 state_global->box,&ddbox,state_global->x);
8319 dd_distribute_state(dd,cgs_gl,
8320 state_global,state_local,f);
8322 dd_make_local_cgs(dd,&top_local->cgs);
8324 if (dd->ncg_home > fr->cg_nalloc)
8326 dd_realloc_fr_cg(fr,dd->ncg_home);
8328 calc_cgcm(fplog,0,dd->ncg_home,
8329 &top_local->cgs,state_local->x,fr->cg_cm);
8331 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
8333 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
8337 else if (state_local->ddp_count != dd->ddp_count)
8339 if (state_local->ddp_count > dd->ddp_count)
8341 gmx_fatal(FARGS,"Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%d)",state_local->ddp_count,dd->ddp_count);
8344 if (state_local->ddp_count_cg_gl != state_local->ddp_count)
8346 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);
8349 /* Clear the old state */
8350 clear_dd_indices(dd,0,0);
8352 /* Build the new indices */
8353 rebuild_cgindex(dd,cgs_gl->index,state_local);
8354 make_dd_indices(dd,cgs_gl->index,0);
8356 /* Redetermine the cg COMs */
8357 calc_cgcm(fplog,0,dd->ncg_home,
8358 &top_local->cgs,state_local->x,fr->cg_cm);
8360 inc_nrnb(nrnb,eNR_CGCM,dd->nat_home);
8362 dd_set_cginfo(dd->index_gl,0,dd->ncg_home,fr,comm->bLocalCG);
8364 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
8365 TRUE,&top_local->cgs,state_local->x,&ddbox);
8367 bRedist = comm->bDynLoadBal;
8371 /* We have the full state, only redistribute the cgs */
8373 /* Clear the non-home indices */
8374 clear_dd_indices(dd,dd->ncg_home,dd->nat_home);
8376 /* Avoid global communication for dim's without pbc and -gcom */
8377 if (!bNStGlobalComm)
8379 copy_rvec(comm->box0 ,ddbox.box0 );
8380 copy_rvec(comm->box_size,ddbox.box_size);
8382 set_ddbox(dd,bMasterState,cr,ir,state_local->box,
8383 bNStGlobalComm,&top_local->cgs,state_local->x,&ddbox);
8388 /* For dim's without pbc and -gcom */
8389 copy_rvec(ddbox.box0 ,comm->box0 );
8390 copy_rvec(ddbox.box_size,comm->box_size);
8392 set_dd_cell_sizes(dd,&ddbox,dynamic_dd_box(&ddbox,ir),bMasterState,bDoDLB,
8395 if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
8397 write_dd_grid_pdb("dd_grid",step,dd,state_local->box,&ddbox);
8400 /* Check if we should sort the charge groups */
8401 if (comm->nstSortCG > 0)
8403 bSortCG = (bMasterState ||
8404 (bRedist && (step % comm->nstSortCG == 0)));
8411 ncg_home_old = dd->ncg_home;
8415 cg0 = dd_redistribute_cg(fplog,step,dd,ddbox.tric_dir,
8416 state_local,f,fr,mdatoms,
8420 get_nsgrid_boundaries(fr->ns.grid,dd,
8421 state_local->box,&ddbox,&comm->cell_x0,&comm->cell_x1,
8422 dd->ncg_home,fr->cg_cm,
8423 cell_ns_x0,cell_ns_x1,&grid_density);
8427 comm_dd_ns_cell_sizes(dd,&ddbox,cell_ns_x0,cell_ns_x1,step);
8430 copy_ivec(fr->ns.grid->n,ncells_old);
8431 grid_first(fplog,fr->ns.grid,dd,&ddbox,fr->ePBC,
8432 state_local->box,cell_ns_x0,cell_ns_x1,
8433 fr->rlistlong,grid_density);
8434 /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
8435 copy_ivec(ddbox.tric_dir,comm->tric_dir);
8439 /* Sort the state on charge group position.
8440 * This enables exact restarts from this step.
8441 * It also improves performance by about 15% with larger numbers
8442 * of atoms per node.
8445 /* Fill the ns grid with the home cell,
8446 * so we can sort with the indices.
8448 set_zones_ncg_home(dd);
8449 fill_grid(fplog,&comm->zones,fr->ns.grid,dd->ncg_home,
8450 0,dd->ncg_home,fr->cg_cm);
8452 /* Check if we can user the old order and ns grid cell indices
8453 * of the charge groups to sort the charge groups efficiently.
8455 bResortAll = (bMasterState ||
8456 fr->ns.grid->n[XX] != ncells_old[XX] ||
8457 fr->ns.grid->n[YY] != ncells_old[YY] ||
8458 fr->ns.grid->n[ZZ] != ncells_old[ZZ]);
8462 fprintf(debug,"Step %s, sorting the %d home charge groups\n",
8463 gmx_step_str(step,sbuf),dd->ncg_home);
8465 dd_sort_state(dd,ir->ePBC,fr->cg_cm,fr,state_local,
8466 bResortAll ? -1 : ncg_home_old);
8467 /* Rebuild all the indices */
8469 ga2la_clear(dd->ga2la);
8472 /* Setup up the communication and communicate the coordinates */
8473 setup_dd_communication(dd,state_local->box,&ddbox,fr);
8475 /* Set the indices */
8476 make_dd_indices(dd,cgs_gl->index,cg0);
8478 /* Set the charge group boundaries for neighbor searching */
8479 set_cg_boundaries(&comm->zones);
8482 write_dd_pdb("dd_home",step,"dump",top_global,cr,
8483 -1,state_local->x,state_local->box);
8486 /* Extract a local topology from the global topology */
8487 for(i=0; i<dd->ndim; i++)
8489 np[dd->dim[i]] = comm->cd[i].np;
8491 dd_make_local_top(fplog,dd,&comm->zones,dd->npbcdim,state_local->box,
8492 comm->cellsize_min,np,
8493 fr,vsite,top_global,top_local);
8495 /* Set up the special atom communication */
8496 n = comm->nat[ddnatZONE];
8497 for(i=ddnatZONE+1; i<ddnatNR; i++)
8502 if (vsite && vsite->n_intercg_vsite)
8504 n = dd_make_local_vsites(dd,n,top_local->idef.il);
8508 if (dd->bInterCGcons)
8510 /* Only for inter-cg constraints we need special code */
8511 n = dd_make_local_constraints(dd,n,top_global,
8512 constr,ir->nProjOrder,
8513 &top_local->idef.il[F_CONSTR]);
8517 gmx_incons("Unknown special atom type setup");
8522 /* Make space for the extra coordinates for virtual site
8523 * or constraint communication.
8525 state_local->natoms = comm->nat[ddnatNR-1];
8526 if (state_local->natoms > state_local->nalloc)
8528 dd_realloc_state(state_local,f,state_local->natoms);
8531 if (fr->bF_NoVirSum)
8533 if (vsite && vsite->n_intercg_vsite)
8535 nat_f_novirsum = comm->nat[ddnatVSITE];
8539 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
8541 nat_f_novirsum = dd->nat_tot;
8545 nat_f_novirsum = dd->nat_home;
8554 /* Set the number of atoms required for the force calculation.
8555 * Forces need to be constrained when using a twin-range setup
8556 * or with energy minimization. For simple simulations we could
8557 * avoid some allocation, zeroing and copying, but this is
8558 * probably not worth the complications ande checking.
8560 forcerec_set_ranges(fr,dd->ncg_home,dd->ncg_tot,
8561 dd->nat_tot,comm->nat[ddnatCON],nat_f_novirsum);
8563 /* We make the all mdatoms up to nat_tot_con.
8564 * We could save some work by only setting invmass
8565 * between nat_tot and nat_tot_con.
8567 /* This call also sets the new number of home particles to dd->nat_home */
8568 atoms2md(top_global,ir,
8569 comm->nat[ddnatCON],dd->gatindex,0,dd->nat_home,mdatoms);
8571 /* Now we have the charges we can sort the FE interactions */
8572 dd_sort_local_top(dd,mdatoms,top_local);
8576 /* Make the local shell stuff, currently no communication is done */
8577 make_local_shells(cr,mdatoms,shellfc);
8580 if (ir->implicit_solvent)
8582 make_local_gb(cr,fr->born,ir->gb_algorithm);
8585 if (!(cr->duty & DUTY_PME))
8587 /* Send the charges to our PME only node */
8588 gmx_pme_send_q(cr,mdatoms->nChargePerturbed,
8589 mdatoms->chargeA,mdatoms->chargeB,
8590 dd_pme_maxshift_x(dd),dd_pme_maxshift_y(dd));
8595 set_constraints(constr,top_local,ir,mdatoms,cr);
8598 if (ir->ePull != epullNO)
8600 /* Update the local pull groups */
8601 dd_make_local_pull_groups(dd,ir->pull,mdatoms);
8606 /* Update the local rotation groups */
8607 dd_make_local_rotation_groups(dd,ir->rot);
8611 add_dd_statistics(dd);
8613 /* Make sure we only count the cycles for this DD partitioning */
8614 clear_dd_cycle_counts(dd);
8616 /* Because the order of the atoms might have changed since
8617 * the last vsite construction, we need to communicate the constructing
8618 * atom coordinates again (for spreading the forces this MD step).
8620 dd_move_x_vsites(dd,state_local->box,state_local->x);
8622 if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
8624 dd_move_x(dd,state_local->box,state_local->x);
8625 write_dd_pdb("dd_dump",step,"dump",top_global,cr,
8626 -1,state_local->x,state_local->box);
8631 /* Store the global communication step */
8632 comm->globalcomm_step = step;
8635 /* Increase the DD partitioning counter */
8637 /* The state currently matches this DD partitioning count, store it */
8638 state_local->ddp_count = dd->ddp_count;
8641 /* The DD master node knows the complete cg distribution,
8642 * store the count so we can possibly skip the cg info communication.
8644 comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
8647 if (comm->DD_debug > 0)
8649 /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
8650 check_index_consistency(dd,top_global->natoms,ncg_mtop(top_global),
8651 "after partitioning");