2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gromacs/domdec/collect.h"
54 #include "gromacs/domdec/dlb.h"
55 #include "gromacs/domdec/dlbtiming.h"
56 #include "gromacs/domdec/domdec_network.h"
57 #include "gromacs/domdec/ga2la.h"
58 #include "gromacs/domdec/options.h"
59 #include "gromacs/domdec/partition.h"
60 #include "gromacs/gmxlib/network.h"
61 #include "gromacs/gmxlib/nrnb.h"
62 #include "gromacs/gpu_utils/gpu_utils.h"
63 #include "gromacs/hardware/hw_info.h"
64 #include "gromacs/listed_forces/manage_threading.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/math/vectypes.h"
67 #include "gromacs/mdlib/calc_verletbuf.h"
68 #include "gromacs/mdlib/constr.h"
69 #include "gromacs/mdlib/constraintrange.h"
70 #include "gromacs/mdlib/updategroups.h"
71 #include "gromacs/mdlib/vsite.h"
72 #include "gromacs/mdtypes/commrec.h"
73 #include "gromacs/mdtypes/forceoutput.h"
74 #include "gromacs/mdtypes/inputrec.h"
75 #include "gromacs/mdtypes/mdrunoptions.h"
76 #include "gromacs/mdtypes/state.h"
77 #include "gromacs/pbcutil/ishift.h"
78 #include "gromacs/pbcutil/pbc.h"
79 #include "gromacs/pulling/pull.h"
80 #include "gromacs/timing/wallcycle.h"
81 #include "gromacs/topology/block.h"
82 #include "gromacs/topology/idef.h"
83 #include "gromacs/topology/ifunc.h"
84 #include "gromacs/topology/mtop_lookup.h"
85 #include "gromacs/topology/mtop_util.h"
86 #include "gromacs/topology/topology.h"
87 #include "gromacs/utility/basedefinitions.h"
88 #include "gromacs/utility/basenetwork.h"
89 #include "gromacs/utility/cstringutil.h"
90 #include "gromacs/utility/exceptions.h"
91 #include "gromacs/utility/fatalerror.h"
92 #include "gromacs/utility/gmxmpi.h"
93 #include "gromacs/utility/logger.h"
94 #include "gromacs/utility/real.h"
95 #include "gromacs/utility/smalloc.h"
96 #include "gromacs/utility/strconvert.h"
97 #include "gromacs/utility/stringstream.h"
98 #include "gromacs/utility/stringutil.h"
99 #include "gromacs/utility/textwriter.h"
101 #include "atomdistribution.h"
103 #include "cellsizes.h"
104 #include "distribute.h"
105 #include "domdec_constraints.h"
106 #include "domdec_internal.h"
107 #include "domdec_vsite.h"
108 #include "redistribute.h"
111 // TODO remove this when moving domdec into gmx namespace
112 using gmx::DdRankOrder
;
113 using gmx::DlbOption
;
114 using gmx::DomdecOptions
;
116 static const char *edlbs_names
[int(DlbState::nr
)] = { "off", "auto", "locked", "on", "on" };
118 /* The size per atom group of the cggl_flag buffer in gmx_domdec_comm_t */
121 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
122 #define DD_FLAG_NRCG 65535
123 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
124 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
126 /* The DD zone order */
127 static const ivec dd_zo
[DD_MAXZONE
] =
128 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
130 /* The non-bonded zone-pair setup for domain decomposition
131 * The first number is the i-zone, the second number the first j-zone seen by
132 * this i-zone, the third number the last+1 j-zone seen by this i-zone.
133 * As is, this is for 3D decomposition, where there are 4 i-zones.
134 * With 2D decomposition use only the first 2 i-zones and a last+1 j-zone of 4.
135 * With 1D decomposition use only the first i-zone and a last+1 j-zone of 2.
138 ddNonbondedZonePairRanges
[DD_MAXIZONE
][3] = {{0, 0, 8},
145 #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
147 static void index2xyz(ivec nc,int ind,ivec xyz)
149 xyz[XX] = ind % nc[XX];
150 xyz[YY] = (ind / nc[XX]) % nc[YY];
151 xyz[ZZ] = ind / (nc[YY]*nc[XX]);
155 static void ddindex2xyz(const ivec nc
, int ind
, ivec xyz
)
157 xyz
[XX
] = ind
/ (nc
[YY
]*nc
[ZZ
]);
158 xyz
[YY
] = (ind
/ nc
[ZZ
]) % nc
[YY
];
159 xyz
[ZZ
] = ind
% nc
[ZZ
];
162 static int ddcoord2ddnodeid(gmx_domdec_t
*dd
, ivec c
)
167 ddindex
= dd_index(dd
->nc
, c
);
168 if (dd
->comm
->bCartesianPP_PME
)
170 ddnodeid
= dd
->comm
->ddindex2ddnodeid
[ddindex
];
172 else if (dd
->comm
->bCartesianPP
)
175 MPI_Cart_rank(dd
->mpi_comm_all
, c
, &ddnodeid
);
186 int ddglatnr(const gmx_domdec_t
*dd
, int i
)
196 if (i
>= dd
->comm
->atomRanges
.numAtomsTotal())
198 gmx_fatal(FARGS
, "glatnr called with %d, which is larger than the local number of atoms (%d)", i
, dd
->comm
->atomRanges
.numAtomsTotal());
200 atnr
= dd
->globalAtomIndices
[i
] + 1;
206 t_block
*dd_charge_groups_global(gmx_domdec_t
*dd
)
208 return &dd
->comm
->cgs_gl
;
211 gmx::ArrayRef
<const gmx::RangePartitioning
> getUpdateGroupingPerMoleculetype(const gmx_domdec_t
&dd
)
213 GMX_RELEASE_ASSERT(dd
.comm
, "Need a valid dd.comm");
214 return dd
.comm
->updateGroupingPerMoleculetype
;
217 void dd_store_state(gmx_domdec_t
*dd
, t_state
*state
)
221 if (state
->ddp_count
!= dd
->ddp_count
)
223 gmx_incons("The MD state does not match the domain decomposition state");
226 state
->cg_gl
.resize(dd
->ncg_home
);
227 for (i
= 0; i
< dd
->ncg_home
; i
++)
229 state
->cg_gl
[i
] = dd
->globalAtomGroupIndices
[i
];
232 state
->ddp_count_cg_gl
= dd
->ddp_count
;
235 gmx_domdec_zones_t
*domdec_zones(gmx_domdec_t
*dd
)
237 return &dd
->comm
->zones
;
240 void dd_get_ns_ranges(const gmx_domdec_t
*dd
, int icg
,
241 int *jcg0
, int *jcg1
, ivec shift0
, ivec shift1
)
243 gmx_domdec_zones_t
*zones
;
246 zones
= &dd
->comm
->zones
;
249 while (icg
>= zones
->izone
[izone
].cg1
)
258 else if (izone
< zones
->nizone
)
260 *jcg0
= zones
->izone
[izone
].jcg0
;
264 gmx_fatal(FARGS
, "DD icg %d out of range: izone (%d) >= nizone (%d)",
265 icg
, izone
, zones
->nizone
);
268 *jcg1
= zones
->izone
[izone
].jcg1
;
270 for (d
= 0; d
< dd
->ndim
; d
++)
273 shift0
[dim
] = zones
->izone
[izone
].shift0
[dim
];
274 shift1
[dim
] = zones
->izone
[izone
].shift1
[dim
];
275 if (dd
->comm
->tric_dir
[dim
] || (isDlbOn(dd
->comm
) && d
> 0))
277 /* A conservative approach, this can be optimized */
284 int dd_numHomeAtoms(const gmx_domdec_t
&dd
)
286 return dd
.comm
->atomRanges
.numHomeAtoms();
289 int dd_natoms_mdatoms(const gmx_domdec_t
*dd
)
291 /* We currently set mdatoms entries for all atoms:
292 * local + non-local + communicated for vsite + constraints
295 return dd
->comm
->atomRanges
.numAtomsTotal();
298 int dd_natoms_vsite(const gmx_domdec_t
*dd
)
300 return dd
->comm
->atomRanges
.end(DDAtomRanges::Type::Vsites
);
303 void dd_get_constraint_range(const gmx_domdec_t
*dd
, int *at_start
, int *at_end
)
305 *at_start
= dd
->comm
->atomRanges
.start(DDAtomRanges::Type::Constraints
);
306 *at_end
= dd
->comm
->atomRanges
.end(DDAtomRanges::Type::Constraints
);
309 void dd_move_x(gmx_domdec_t
*dd
,
311 gmx::ArrayRef
<gmx::RVec
> x
,
312 gmx_wallcycle
*wcycle
)
314 wallcycle_start(wcycle
, ewcMOVEX
);
317 gmx_domdec_comm_t
*comm
;
318 gmx_domdec_comm_dim_t
*cd
;
319 rvec shift
= {0, 0, 0};
320 gmx_bool bPBC
, bScrew
;
325 nat_tot
= comm
->atomRanges
.numHomeAtoms();
326 for (int d
= 0; d
< dd
->ndim
; d
++)
328 bPBC
= (dd
->ci
[dd
->dim
[d
]] == 0);
329 bScrew
= (bPBC
&& dd
->unitCellInfo
.haveScrewPBC
&& dd
->dim
[d
] == XX
);
332 copy_rvec(box
[dd
->dim
[d
]], shift
);
335 for (const gmx_domdec_ind_t
&ind
: cd
->ind
)
337 DDBufferAccess
<gmx::RVec
> sendBufferAccess(comm
->rvecBuffer
, ind
.nsend
[nzone
+ 1]);
338 gmx::ArrayRef
<gmx::RVec
> &sendBuffer
= sendBufferAccess
.buffer
;
342 for (int j
: ind
.index
)
344 sendBuffer
[n
] = x
[j
];
350 for (int j
: ind
.index
)
352 /* We need to shift the coordinates */
353 for (int d
= 0; d
< DIM
; d
++)
355 sendBuffer
[n
][d
] = x
[j
][d
] + shift
[d
];
362 for (int j
: ind
.index
)
365 sendBuffer
[n
][XX
] = x
[j
][XX
] + shift
[XX
];
367 * This operation requires a special shift force
368 * treatment, which is performed in calc_vir.
370 sendBuffer
[n
][YY
] = box
[YY
][YY
] - x
[j
][YY
];
371 sendBuffer
[n
][ZZ
] = box
[ZZ
][ZZ
] - x
[j
][ZZ
];
376 DDBufferAccess
<gmx::RVec
> receiveBufferAccess(comm
->rvecBuffer2
, cd
->receiveInPlace
? 0 : ind
.nrecv
[nzone
+ 1]);
378 gmx::ArrayRef
<gmx::RVec
> receiveBuffer
;
379 if (cd
->receiveInPlace
)
381 receiveBuffer
= gmx::arrayRefFromArray(x
.data() + nat_tot
, ind
.nrecv
[nzone
+ 1]);
385 receiveBuffer
= receiveBufferAccess
.buffer
;
387 /* Send and receive the coordinates */
388 ddSendrecv(dd
, d
, dddirBackward
,
389 sendBuffer
, receiveBuffer
);
391 if (!cd
->receiveInPlace
)
394 for (int zone
= 0; zone
< nzone
; zone
++)
396 for (int i
= ind
.cell2at0
[zone
]; i
< ind
.cell2at1
[zone
]; i
++)
398 x
[i
] = receiveBuffer
[j
++];
402 nat_tot
+= ind
.nrecv
[nzone
+1];
407 wallcycle_stop(wcycle
, ewcMOVEX
);
410 void dd_move_f(gmx_domdec_t
*dd
,
411 gmx::ForceWithShiftForces
*forceWithShiftForces
,
412 gmx_wallcycle
*wcycle
)
414 wallcycle_start(wcycle
, ewcMOVEF
);
416 gmx::ArrayRef
<gmx::RVec
> f
= forceWithShiftForces
->force();
417 gmx::ArrayRef
<gmx::RVec
> fshift
= forceWithShiftForces
->shiftForces();
419 gmx_domdec_comm_t
&comm
= *dd
->comm
;
420 int nzone
= comm
.zones
.n
/2;
421 int nat_tot
= comm
.atomRanges
.end(DDAtomRanges::Type::Zones
);
422 for (int d
= dd
->ndim
-1; d
>= 0; d
--)
424 /* Only forces in domains near the PBC boundaries need to
425 consider PBC in the treatment of fshift */
426 const bool shiftForcesNeedPbc
= (forceWithShiftForces
->computeVirial() && dd
->ci
[dd
->dim
[d
]] == 0);
427 const bool applyScrewPbc
= (shiftForcesNeedPbc
&& dd
->unitCellInfo
.haveScrewPBC
&& dd
->dim
[d
] == XX
);
428 /* Determine which shift vector we need */
429 ivec vis
= { 0, 0, 0 };
431 const int is
= IVEC2IS(vis
);
433 /* Loop over the pulses */
434 const gmx_domdec_comm_dim_t
&cd
= comm
.cd
[d
];
435 for (int p
= cd
.numPulses() - 1; p
>= 0; p
--)
437 const gmx_domdec_ind_t
&ind
= cd
.ind
[p
];
438 DDBufferAccess
<gmx::RVec
> receiveBufferAccess(comm
.rvecBuffer
, ind
.nsend
[nzone
+ 1]);
439 gmx::ArrayRef
<gmx::RVec
> &receiveBuffer
= receiveBufferAccess
.buffer
;
441 nat_tot
-= ind
.nrecv
[nzone
+1];
443 DDBufferAccess
<gmx::RVec
> sendBufferAccess(comm
.rvecBuffer2
, cd
.receiveInPlace
? 0 : ind
.nrecv
[nzone
+ 1]);
445 gmx::ArrayRef
<gmx::RVec
> sendBuffer
;
446 if (cd
.receiveInPlace
)
448 sendBuffer
= gmx::arrayRefFromArray(f
.data() + nat_tot
, ind
.nrecv
[nzone
+ 1]);
452 sendBuffer
= sendBufferAccess
.buffer
;
454 for (int zone
= 0; zone
< nzone
; zone
++)
456 for (int i
= ind
.cell2at0
[zone
]; i
< ind
.cell2at1
[zone
]; i
++)
458 sendBuffer
[j
++] = f
[i
];
462 /* Communicate the forces */
463 ddSendrecv(dd
, d
, dddirForward
,
464 sendBuffer
, receiveBuffer
);
465 /* Add the received forces */
467 if (!shiftForcesNeedPbc
)
469 for (int j
: ind
.index
)
471 for (int d
= 0; d
< DIM
; d
++)
473 f
[j
][d
] += receiveBuffer
[n
][d
];
478 else if (!applyScrewPbc
)
480 for (int j
: ind
.index
)
482 for (int d
= 0; d
< DIM
; d
++)
484 f
[j
][d
] += receiveBuffer
[n
][d
];
486 /* Add this force to the shift force */
487 for (int d
= 0; d
< DIM
; d
++)
489 fshift
[is
][d
] += receiveBuffer
[n
][d
];
496 for (int j
: ind
.index
)
498 /* Rotate the force */
499 f
[j
][XX
] += receiveBuffer
[n
][XX
];
500 f
[j
][YY
] -= receiveBuffer
[n
][YY
];
501 f
[j
][ZZ
] -= receiveBuffer
[n
][ZZ
];
502 if (shiftForcesNeedPbc
)
504 /* Add this force to the shift force */
505 for (int d
= 0; d
< DIM
; d
++)
507 fshift
[is
][d
] += receiveBuffer
[n
][d
];
516 wallcycle_stop(wcycle
, ewcMOVEF
);
519 /* Convenience function for extracting a real buffer from an rvec buffer
521 * To reduce the number of temporary communication buffers and avoid
522 * cache polution, we reuse gmx::RVec buffers for storing reals.
523 * This functions return a real buffer reference with the same number
524 * of elements as the gmx::RVec buffer (so 1/3 of the size in bytes).
526 static gmx::ArrayRef
<real
>
527 realArrayRefFromRvecArrayRef(gmx::ArrayRef
<gmx::RVec
> arrayRef
)
529 return gmx::arrayRefFromArray(as_rvec_array(arrayRef
.data())[0],
533 void dd_atom_spread_real(gmx_domdec_t
*dd
, real v
[])
536 gmx_domdec_comm_t
*comm
;
537 gmx_domdec_comm_dim_t
*cd
;
542 nat_tot
= comm
->atomRanges
.numHomeAtoms();
543 for (int d
= 0; d
< dd
->ndim
; d
++)
546 for (const gmx_domdec_ind_t
&ind
: cd
->ind
)
548 /* Note: We provision for RVec instead of real, so a factor of 3
549 * more than needed. The buffer actually already has this size
550 * and we pass a plain pointer below, so this does not matter.
552 DDBufferAccess
<gmx::RVec
> sendBufferAccess(comm
->rvecBuffer
, ind
.nsend
[nzone
+ 1]);
553 gmx::ArrayRef
<real
> sendBuffer
= realArrayRefFromRvecArrayRef(sendBufferAccess
.buffer
);
555 for (int j
: ind
.index
)
557 sendBuffer
[n
++] = v
[j
];
560 DDBufferAccess
<gmx::RVec
> receiveBufferAccess(comm
->rvecBuffer2
, cd
->receiveInPlace
? 0 : ind
.nrecv
[nzone
+ 1]);
562 gmx::ArrayRef
<real
> receiveBuffer
;
563 if (cd
->receiveInPlace
)
565 receiveBuffer
= gmx::arrayRefFromArray(v
+ nat_tot
, ind
.nrecv
[nzone
+ 1]);
569 receiveBuffer
= realArrayRefFromRvecArrayRef(receiveBufferAccess
.buffer
);
571 /* Send and receive the data */
572 ddSendrecv(dd
, d
, dddirBackward
,
573 sendBuffer
, receiveBuffer
);
574 if (!cd
->receiveInPlace
)
577 for (int zone
= 0; zone
< nzone
; zone
++)
579 for (int i
= ind
.cell2at0
[zone
]; i
< ind
.cell2at1
[zone
]; i
++)
581 v
[i
] = receiveBuffer
[j
++];
585 nat_tot
+= ind
.nrecv
[nzone
+1];
591 void dd_atom_sum_real(gmx_domdec_t
*dd
, real v
[])
594 gmx_domdec_comm_t
*comm
;
595 gmx_domdec_comm_dim_t
*cd
;
599 nzone
= comm
->zones
.n
/2;
600 nat_tot
= comm
->atomRanges
.end(DDAtomRanges::Type::Zones
);
601 for (int d
= dd
->ndim
-1; d
>= 0; d
--)
604 for (int p
= cd
->numPulses() - 1; p
>= 0; p
--)
606 const gmx_domdec_ind_t
&ind
= cd
->ind
[p
];
608 /* Note: We provision for RVec instead of real, so a factor of 3
609 * more than needed. The buffer actually already has this size
610 * and we typecast, so this works as intended.
612 DDBufferAccess
<gmx::RVec
> receiveBufferAccess(comm
->rvecBuffer
, ind
.nsend
[nzone
+ 1]);
613 gmx::ArrayRef
<real
> receiveBuffer
= realArrayRefFromRvecArrayRef(receiveBufferAccess
.buffer
);
614 nat_tot
-= ind
.nrecv
[nzone
+ 1];
616 DDBufferAccess
<gmx::RVec
> sendBufferAccess(comm
->rvecBuffer2
, cd
->receiveInPlace
? 0 : ind
.nrecv
[nzone
+ 1]);
618 gmx::ArrayRef
<real
> sendBuffer
;
619 if (cd
->receiveInPlace
)
621 sendBuffer
= gmx::arrayRefFromArray(v
+ nat_tot
, ind
.nrecv
[nzone
+ 1]);
625 sendBuffer
= realArrayRefFromRvecArrayRef(sendBufferAccess
.buffer
);
627 for (int zone
= 0; zone
< nzone
; zone
++)
629 for (int i
= ind
.cell2at0
[zone
]; i
< ind
.cell2at1
[zone
]; i
++)
631 sendBuffer
[j
++] = v
[i
];
635 /* Communicate the forces */
636 ddSendrecv(dd
, d
, dddirForward
,
637 sendBuffer
, receiveBuffer
);
638 /* Add the received forces */
640 for (int j
: ind
.index
)
642 v
[j
] += receiveBuffer
[n
];
650 real
dd_cutoff_multibody(const gmx_domdec_t
*dd
)
652 gmx_domdec_comm_t
*comm
;
659 if (comm
->haveInterDomainMultiBodyBondeds
)
661 if (comm
->cutoff_mbody
> 0)
663 r
= comm
->cutoff_mbody
;
667 /* cutoff_mbody=0 means we do not have DLB */
668 r
= comm
->cellsize_min
[dd
->dim
[0]];
669 for (di
= 1; di
< dd
->ndim
; di
++)
671 r
= std::min(r
, comm
->cellsize_min
[dd
->dim
[di
]]);
675 r
= std::max(r
, comm
->cutoff_mbody
);
679 r
= std::min(r
, comm
->cutoff
);
687 real
dd_cutoff_twobody(const gmx_domdec_t
*dd
)
691 r_mb
= dd_cutoff_multibody(dd
);
693 return std::max(dd
->comm
->cutoff
, r_mb
);
697 static void dd_cart_coord2pmecoord(const gmx_domdec_t
*dd
, const ivec coord
,
702 nc
= dd
->nc
[dd
->comm
->cartpmedim
];
703 ntot
= dd
->comm
->ntot
[dd
->comm
->cartpmedim
];
704 copy_ivec(coord
, coord_pme
);
705 coord_pme
[dd
->comm
->cartpmedim
] =
706 nc
+ (coord
[dd
->comm
->cartpmedim
]*(ntot
- nc
) + (ntot
- nc
)/2)/nc
;
709 static int ddindex2pmeindex(const gmx_domdec_t
*dd
, int ddindex
)
714 npme
= dd
->comm
->npmenodes
;
716 /* Here we assign a PME node to communicate with this DD node
717 * by assuming that the major index of both is x.
718 * We add cr->npmenodes/2 to obtain an even distribution.
720 return (ddindex
*npme
+ npme
/2)/npp
;
723 static int *dd_interleaved_pme_ranks(const gmx_domdec_t
*dd
)
728 snew(pme_rank
, dd
->comm
->npmenodes
);
730 for (i
= 0; i
< dd
->nnodes
; i
++)
732 p0
= ddindex2pmeindex(dd
, i
);
733 p1
= ddindex2pmeindex(dd
, i
+1);
734 if (i
+1 == dd
->nnodes
|| p1
> p0
)
738 fprintf(debug
, "pme_rank[%d] = %d\n", n
, i
+1+n
);
740 pme_rank
[n
] = i
+ 1 + n
;
748 static int gmx_ddcoord2pmeindex(const t_commrec
*cr
, int x
, int y
, int z
)
756 if (dd->comm->bCartesian) {
757 gmx_ddindex2xyz(dd->nc,ddindex,coords);
758 dd_coords2pmecoords(dd,coords,coords_pme);
759 copy_ivec(dd->ntot,nc);
760 nc[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
761 coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
763 slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
765 slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
771 slab
= ddindex2pmeindex(dd
, dd_index(dd
->nc
, coords
));
776 static int ddcoord2simnodeid(const t_commrec
*cr
, int x
, int y
, int z
)
778 gmx_domdec_comm_t
*comm
;
780 int ddindex
, nodeid
= -1;
787 if (comm
->bCartesianPP_PME
)
790 MPI_Cart_rank(cr
->mpi_comm_mysim
, coords
, &nodeid
);
795 ddindex
= dd_index(cr
->dd
->nc
, coords
);
796 if (comm
->bCartesianPP
)
798 nodeid
= comm
->ddindex2simnodeid
[ddindex
];
804 nodeid
= ddindex
+ gmx_ddcoord2pmeindex(cr
, x
, y
, z
);
816 static int dd_simnode2pmenode(const gmx_domdec_t
*dd
,
817 const t_commrec gmx_unused
*cr
,
822 const gmx_domdec_comm_t
*comm
= dd
->comm
;
824 /* This assumes a uniform x domain decomposition grid cell size */
825 if (comm
->bCartesianPP_PME
)
828 ivec coord
, coord_pme
;
829 MPI_Cart_coords(cr
->mpi_comm_mysim
, sim_nodeid
, DIM
, coord
);
830 if (coord
[comm
->cartpmedim
] < dd
->nc
[comm
->cartpmedim
])
832 /* This is a PP node */
833 dd_cart_coord2pmecoord(dd
, coord
, coord_pme
);
834 MPI_Cart_rank(cr
->mpi_comm_mysim
, coord_pme
, &pmenode
);
838 else if (comm
->bCartesianPP
)
840 if (sim_nodeid
< dd
->nnodes
)
842 pmenode
= dd
->nnodes
+ ddindex2pmeindex(dd
, sim_nodeid
);
847 /* This assumes DD cells with identical x coordinates
848 * are numbered sequentially.
850 if (dd
->comm
->pmenodes
== nullptr)
852 if (sim_nodeid
< dd
->nnodes
)
854 /* The DD index equals the nodeid */
855 pmenode
= dd
->nnodes
+ ddindex2pmeindex(dd
, sim_nodeid
);
861 while (sim_nodeid
> dd
->comm
->pmenodes
[i
])
865 if (sim_nodeid
< dd
->comm
->pmenodes
[i
])
867 pmenode
= dd
->comm
->pmenodes
[i
];
875 NumPmeDomains
getNumPmeDomains(const gmx_domdec_t
*dd
)
880 dd
->comm
->npmenodes_x
, dd
->comm
->npmenodes_y
891 std::vector
<int> get_pme_ddranks(const t_commrec
*cr
, int pmenodeid
)
895 ivec coord
, coord_pme
;
899 std::vector
<int> ddranks
;
900 ddranks
.reserve((dd
->nnodes
+cr
->npmenodes
-1)/cr
->npmenodes
);
902 for (x
= 0; x
< dd
->nc
[XX
]; x
++)
904 for (y
= 0; y
< dd
->nc
[YY
]; y
++)
906 for (z
= 0; z
< dd
->nc
[ZZ
]; z
++)
908 if (dd
->comm
->bCartesianPP_PME
)
913 dd_cart_coord2pmecoord(dd
, coord
, coord_pme
);
914 if (dd
->ci
[XX
] == coord_pme
[XX
] &&
915 dd
->ci
[YY
] == coord_pme
[YY
] &&
916 dd
->ci
[ZZ
] == coord_pme
[ZZ
])
918 ddranks
.push_back(ddcoord2simnodeid(cr
, x
, y
, z
));
923 /* The slab corresponds to the nodeid in the PME group */
924 if (gmx_ddcoord2pmeindex(cr
, x
, y
, z
) == pmenodeid
)
926 ddranks
.push_back(ddcoord2simnodeid(cr
, x
, y
, z
));
935 static gmx_bool
receive_vir_ener(const gmx_domdec_t
*dd
, const t_commrec
*cr
)
937 gmx_bool bReceive
= TRUE
;
939 if (cr
->npmenodes
< dd
->nnodes
)
941 gmx_domdec_comm_t
*comm
= dd
->comm
;
942 if (comm
->bCartesianPP_PME
)
945 int pmenode
= dd_simnode2pmenode(dd
, cr
, cr
->sim_nodeid
);
947 MPI_Cart_coords(cr
->mpi_comm_mysim
, cr
->sim_nodeid
, DIM
, coords
);
948 coords
[comm
->cartpmedim
]++;
949 if (coords
[comm
->cartpmedim
] < dd
->nc
[comm
->cartpmedim
])
952 MPI_Cart_rank(cr
->mpi_comm_mysim
, coords
, &rank
);
953 if (dd_simnode2pmenode(dd
, cr
, rank
) == pmenode
)
955 /* This is not the last PP node for pmenode */
960 GMX_RELEASE_ASSERT(false, "Without MPI we should not have Cartesian PP-PME with #PMEnodes < #DDnodes");
965 int pmenode
= dd_simnode2pmenode(dd
, cr
, cr
->sim_nodeid
);
966 if (cr
->sim_nodeid
+1 < cr
->nnodes
&&
967 dd_simnode2pmenode(dd
, cr
, cr
->sim_nodeid
+1) == pmenode
)
969 /* This is not the last PP node for pmenode */
978 static void set_slb_pme_dim_f(gmx_domdec_t
*dd
, int dim
, real
**dim_f
)
980 gmx_domdec_comm_t
*comm
;
985 snew(*dim_f
, dd
->nc
[dim
]+1);
987 for (i
= 1; i
< dd
->nc
[dim
]; i
++)
989 if (comm
->slb_frac
[dim
])
991 (*dim_f
)[i
] = (*dim_f
)[i
-1] + comm
->slb_frac
[dim
][i
-1];
995 (*dim_f
)[i
] = static_cast<real
>(i
)/static_cast<real
>(dd
->nc
[dim
]);
998 (*dim_f
)[dd
->nc
[dim
]] = 1;
1001 static void init_ddpme(gmx_domdec_t
*dd
, gmx_ddpme_t
*ddpme
, int dimind
)
1003 int pmeindex
, slab
, nso
, i
;
1006 if (dimind
== 0 && dd
->dim
[0] == YY
&& dd
->comm
->npmenodes_x
== 1)
1012 ddpme
->dim
= dimind
;
1014 ddpme
->dim_match
= (ddpme
->dim
== dd
->dim
[dimind
]);
1016 ddpme
->nslab
= (ddpme
->dim
== 0 ?
1017 dd
->comm
->npmenodes_x
:
1018 dd
->comm
->npmenodes_y
);
1020 if (ddpme
->nslab
<= 1)
1025 nso
= dd
->comm
->npmenodes
/ddpme
->nslab
;
1026 /* Determine for each PME slab the PP location range for dimension dim */
1027 snew(ddpme
->pp_min
, ddpme
->nslab
);
1028 snew(ddpme
->pp_max
, ddpme
->nslab
);
1029 for (slab
= 0; slab
< ddpme
->nslab
; slab
++)
1031 ddpme
->pp_min
[slab
] = dd
->nc
[dd
->dim
[dimind
]] - 1;
1032 ddpme
->pp_max
[slab
] = 0;
1034 for (i
= 0; i
< dd
->nnodes
; i
++)
1036 ddindex2xyz(dd
->nc
, i
, xyz
);
1037 /* For y only use our y/z slab.
1038 * This assumes that the PME x grid size matches the DD grid size.
1040 if (dimind
== 0 || xyz
[XX
] == dd
->ci
[XX
])
1042 pmeindex
= ddindex2pmeindex(dd
, i
);
1045 slab
= pmeindex
/nso
;
1049 slab
= pmeindex
% ddpme
->nslab
;
1051 ddpme
->pp_min
[slab
] = std::min(ddpme
->pp_min
[slab
], xyz
[dimind
]);
1052 ddpme
->pp_max
[slab
] = std::max(ddpme
->pp_max
[slab
], xyz
[dimind
]);
1056 set_slb_pme_dim_f(dd
, ddpme
->dim
, &ddpme
->slb_dim_f
);
1059 int dd_pme_maxshift_x(const gmx_domdec_t
*dd
)
1061 if (dd
->comm
->ddpme
[0].dim
== XX
)
1063 return dd
->comm
->ddpme
[0].maxshift
;
1071 int dd_pme_maxshift_y(const gmx_domdec_t
*dd
)
1073 if (dd
->comm
->ddpme
[0].dim
== YY
)
1075 return dd
->comm
->ddpme
[0].maxshift
;
1077 else if (dd
->comm
->npmedecompdim
>= 2 && dd
->comm
->ddpme
[1].dim
== YY
)
1079 return dd
->comm
->ddpme
[1].maxshift
;
1087 void dd_cycles_add(const gmx_domdec_t
*dd
, float cycles
, int ddCycl
)
1089 /* Note that the cycles value can be incorrect, either 0 or some
1090 * extremely large value, when our thread migrated to another core
1091 * with an unsynchronized cycle counter. If this happens less often
1092 * that once per nstlist steps, this will not cause issues, since
1093 * we later subtract the maximum value from the sum over nstlist steps.
1094 * A zero count will slightly lower the total, but that's a small effect.
1095 * Note that the main purpose of the subtraction of the maximum value
1096 * is to avoid throwing off the load balancing when stalls occur due
1097 * e.g. system activity or network congestion.
1099 dd
->comm
->cycl
[ddCycl
] += cycles
;
1100 dd
->comm
->cycl_n
[ddCycl
]++;
1101 if (cycles
> dd
->comm
->cycl_max
[ddCycl
])
1103 dd
->comm
->cycl_max
[ddCycl
] = cycles
;
1108 static void make_load_communicator(gmx_domdec_t
*dd
, int dim_ind
, ivec loc
)
1113 gmx_bool bPartOfGroup
= FALSE
;
1115 dim
= dd
->dim
[dim_ind
];
1116 copy_ivec(loc
, loc_c
);
1117 for (i
= 0; i
< dd
->nc
[dim
]; i
++)
1120 rank
= dd_index(dd
->nc
, loc_c
);
1121 if (rank
== dd
->rank
)
1123 /* This process is part of the group */
1124 bPartOfGroup
= TRUE
;
1127 MPI_Comm_split(dd
->mpi_comm_all
, bPartOfGroup
? 0 : MPI_UNDEFINED
, dd
->rank
,
1131 dd
->comm
->mpi_comm_load
[dim_ind
] = c_row
;
1132 if (!isDlbDisabled(dd
->comm
))
1134 DDCellsizesWithDlb
&cellsizes
= dd
->comm
->cellsizesWithDlb
[dim_ind
];
1136 if (dd
->ci
[dim
] == dd
->master_ci
[dim
])
1138 /* This is the root process of this row */
1139 cellsizes
.rowMaster
= std::make_unique
<RowMaster
>();
1141 RowMaster
&rowMaster
= *cellsizes
.rowMaster
;
1142 rowMaster
.cellFrac
.resize(ddCellFractionBufferSize(dd
, dim_ind
));
1143 rowMaster
.oldCellFrac
.resize(dd
->nc
[dim
] + 1);
1144 rowMaster
.isCellMin
.resize(dd
->nc
[dim
]);
1147 rowMaster
.bounds
.resize(dd
->nc
[dim
]);
1149 rowMaster
.buf_ncd
.resize(dd
->nc
[dim
]);
1153 /* This is not a root process, we only need to receive cell_f */
1154 cellsizes
.fracRow
.resize(ddCellFractionBufferSize(dd
, dim_ind
));
1157 if (dd
->ci
[dim
] == dd
->master_ci
[dim
])
1159 snew(dd
->comm
->load
[dim_ind
].load
, dd
->nc
[dim
]*DD_NLOAD_MAX
);
1165 void dd_setup_dlb_resource_sharing(t_commrec
*cr
,
1169 int physicalnode_id_hash
;
1171 MPI_Comm mpi_comm_pp_physicalnode
;
1173 if (!thisRankHasDuty(cr
, DUTY_PP
) || gpu_id
< 0)
1175 /* Only ranks with short-ranged tasks (currently) use GPUs.
1176 * If we don't have GPUs assigned, there are no resources to share.
1181 physicalnode_id_hash
= gmx_physicalnode_id_hash();
1187 fprintf(debug
, "dd_setup_dd_dlb_gpu_sharing:\n");
1188 fprintf(debug
, "DD PP rank %d physical node hash %d gpu_id %d\n",
1189 dd
->rank
, physicalnode_id_hash
, gpu_id
);
1191 /* Split the PP communicator over the physical nodes */
1192 /* TODO: See if we should store this (before), as it's also used for
1193 * for the nodecomm summation.
1195 // TODO PhysicalNodeCommunicator could be extended/used to handle
1196 // the need for per-node per-group communicators.
1197 MPI_Comm_split(dd
->mpi_comm_all
, physicalnode_id_hash
, dd
->rank
,
1198 &mpi_comm_pp_physicalnode
);
1199 MPI_Comm_split(mpi_comm_pp_physicalnode
, gpu_id
, dd
->rank
,
1200 &dd
->comm
->mpi_comm_gpu_shared
);
1201 MPI_Comm_free(&mpi_comm_pp_physicalnode
);
1202 MPI_Comm_size(dd
->comm
->mpi_comm_gpu_shared
, &dd
->comm
->nrank_gpu_shared
);
1206 fprintf(debug
, "nrank_gpu_shared %d\n", dd
->comm
->nrank_gpu_shared
);
1209 /* Note that some ranks could share a GPU, while others don't */
1211 if (dd
->comm
->nrank_gpu_shared
== 1)
1213 MPI_Comm_free(&dd
->comm
->mpi_comm_gpu_shared
);
1216 GMX_UNUSED_VALUE(cr
);
1217 GMX_UNUSED_VALUE(gpu_id
);
1221 static void make_load_communicators(gmx_domdec_t gmx_unused
*dd
)
1224 int dim0
, dim1
, i
, j
;
1229 fprintf(debug
, "Making load communicators\n");
1232 dd
->comm
->load
= new domdec_load_t
[std::max(dd
->ndim
, 1)];
1233 snew(dd
->comm
->mpi_comm_load
, std::max(dd
->ndim
, 1));
1241 make_load_communicator(dd
, 0, loc
);
1245 for (i
= 0; i
< dd
->nc
[dim0
]; i
++)
1248 make_load_communicator(dd
, 1, loc
);
1254 for (i
= 0; i
< dd
->nc
[dim0
]; i
++)
1258 for (j
= 0; j
< dd
->nc
[dim1
]; j
++)
1261 make_load_communicator(dd
, 2, loc
);
1268 fprintf(debug
, "Finished making load communicators\n");
1273 /*! \brief Sets up the relation between neighboring domains and zones */
1274 static void setup_neighbor_relations(gmx_domdec_t
*dd
)
1276 int d
, dim
, i
, j
, m
;
1278 gmx_domdec_zones_t
*zones
;
1279 gmx_domdec_ns_ranges_t
*izone
;
1280 GMX_ASSERT(dd
->ndim
>= 0, "Must have non-negative number of dimensions for DD");
1282 for (d
= 0; d
< dd
->ndim
; d
++)
1285 copy_ivec(dd
->ci
, tmp
);
1286 tmp
[dim
] = (tmp
[dim
] + 1) % dd
->nc
[dim
];
1287 dd
->neighbor
[d
][0] = ddcoord2ddnodeid(dd
, tmp
);
1288 copy_ivec(dd
->ci
, tmp
);
1289 tmp
[dim
] = (tmp
[dim
] - 1 + dd
->nc
[dim
]) % dd
->nc
[dim
];
1290 dd
->neighbor
[d
][1] = ddcoord2ddnodeid(dd
, tmp
);
1293 fprintf(debug
, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
1296 dd
->neighbor
[d
][1]);
1300 int nzone
= (1 << dd
->ndim
);
1301 GMX_ASSERT(dd
->ndim
< DIM
, "Invalid number of dimensions");
1302 int nizone
= (1 << std::max(dd
->ndim
- 1, 0));
1303 assert(nizone
>= 1 && nizone
<= DD_MAXIZONE
);
1305 zones
= &dd
->comm
->zones
;
1307 for (i
= 0; i
< nzone
; i
++)
1310 clear_ivec(zones
->shift
[i
]);
1311 for (d
= 0; d
< dd
->ndim
; d
++)
1313 zones
->shift
[i
][dd
->dim
[d
]] = dd_zo
[i
][m
++];
1318 for (i
= 0; i
< nzone
; i
++)
1320 for (d
= 0; d
< DIM
; d
++)
1322 s
[d
] = dd
->ci
[d
] - zones
->shift
[i
][d
];
1327 else if (s
[d
] >= dd
->nc
[d
])
1333 zones
->nizone
= nizone
;
1334 for (i
= 0; i
< zones
->nizone
; i
++)
1336 assert(ddNonbondedZonePairRanges
[i
][0] == i
);
1338 izone
= &zones
->izone
[i
];
1339 /* dd_zp3 is for 3D decomposition, for fewer dimensions use only
1340 * j-zones up to nzone.
1342 izone
->j0
= std::min(ddNonbondedZonePairRanges
[i
][1], nzone
);
1343 izone
->j1
= std::min(ddNonbondedZonePairRanges
[i
][2], nzone
);
1344 for (dim
= 0; dim
< DIM
; dim
++)
1346 if (dd
->nc
[dim
] == 1)
1348 /* All shifts should be allowed */
1349 izone
->shift0
[dim
] = -1;
1350 izone
->shift1
[dim
] = 1;
1354 /* Determine the min/max j-zone shift wrt the i-zone */
1355 izone
->shift0
[dim
] = 1;
1356 izone
->shift1
[dim
] = -1;
1357 for (j
= izone
->j0
; j
< izone
->j1
; j
++)
1359 int shift_diff
= zones
->shift
[j
][dim
] - zones
->shift
[i
][dim
];
1360 if (shift_diff
< izone
->shift0
[dim
])
1362 izone
->shift0
[dim
] = shift_diff
;
1364 if (shift_diff
> izone
->shift1
[dim
])
1366 izone
->shift1
[dim
] = shift_diff
;
1373 if (!isDlbDisabled(dd
->comm
))
1375 dd
->comm
->cellsizesWithDlb
.resize(dd
->ndim
);
1378 if (dd
->comm
->bRecordLoad
)
1380 make_load_communicators(dd
);
1384 static void make_pp_communicator(const gmx::MDLogger
&mdlog
,
1386 t_commrec gmx_unused
*cr
,
1387 bool gmx_unused reorder
)
1390 gmx_domdec_comm_t
*comm
;
1397 if (comm
->bCartesianPP
)
1399 /* Set up cartesian communication for the particle-particle part */
1400 GMX_LOG(mdlog
.info
).appendTextFormatted(
1401 "Will use a Cartesian communicator: %d x %d x %d",
1402 dd
->nc
[XX
], dd
->nc
[YY
], dd
->nc
[ZZ
]);
1404 for (int i
= 0; i
< DIM
; i
++)
1408 MPI_Cart_create(cr
->mpi_comm_mygroup
, DIM
, dd
->nc
, periods
, static_cast<int>(reorder
),
1410 /* We overwrite the old communicator with the new cartesian one */
1411 cr
->mpi_comm_mygroup
= comm_cart
;
1414 dd
->mpi_comm_all
= cr
->mpi_comm_mygroup
;
1415 MPI_Comm_rank(dd
->mpi_comm_all
, &dd
->rank
);
1417 if (comm
->bCartesianPP_PME
)
1419 /* Since we want to use the original cartesian setup for sim,
1420 * and not the one after split, we need to make an index.
1422 snew(comm
->ddindex2ddnodeid
, dd
->nnodes
);
1423 comm
->ddindex2ddnodeid
[dd_index(dd
->nc
, dd
->ci
)] = dd
->rank
;
1424 gmx_sumi(dd
->nnodes
, comm
->ddindex2ddnodeid
, cr
);
1425 /* Get the rank of the DD master,
1426 * above we made sure that the master node is a PP node.
1436 MPI_Allreduce(&rank
, &dd
->masterrank
, 1, MPI_INT
, MPI_SUM
, dd
->mpi_comm_all
);
1438 else if (comm
->bCartesianPP
)
1440 if (cr
->npmenodes
== 0)
1442 /* The PP communicator is also
1443 * the communicator for this simulation
1445 cr
->mpi_comm_mysim
= cr
->mpi_comm_mygroup
;
1447 cr
->nodeid
= dd
->rank
;
1449 MPI_Cart_coords(dd
->mpi_comm_all
, dd
->rank
, DIM
, dd
->ci
);
1451 /* We need to make an index to go from the coordinates
1452 * to the nodeid of this simulation.
1454 snew(comm
->ddindex2simnodeid
, dd
->nnodes
);
1455 snew(buf
, dd
->nnodes
);
1456 if (thisRankHasDuty(cr
, DUTY_PP
))
1458 buf
[dd_index(dd
->nc
, dd
->ci
)] = cr
->sim_nodeid
;
1460 /* Communicate the ddindex to simulation nodeid index */
1461 MPI_Allreduce(buf
, comm
->ddindex2simnodeid
, dd
->nnodes
, MPI_INT
, MPI_SUM
,
1462 cr
->mpi_comm_mysim
);
1465 /* Determine the master coordinates and rank.
1466 * The DD master should be the same node as the master of this sim.
1468 for (int i
= 0; i
< dd
->nnodes
; i
++)
1470 if (comm
->ddindex2simnodeid
[i
] == 0)
1472 ddindex2xyz(dd
->nc
, i
, dd
->master_ci
);
1473 MPI_Cart_rank(dd
->mpi_comm_all
, dd
->master_ci
, &dd
->masterrank
);
1478 fprintf(debug
, "The master rank is %d\n", dd
->masterrank
);
1483 /* No Cartesian communicators */
1484 /* We use the rank in dd->comm->all as DD index */
1485 ddindex2xyz(dd
->nc
, dd
->rank
, dd
->ci
);
1486 /* The simulation master nodeid is 0, so the DD master rank is also 0 */
1488 clear_ivec(dd
->master_ci
);
1492 GMX_LOG(mdlog
.info
).appendTextFormatted(
1493 "Domain decomposition rank %d, coordinates %d %d %d\n",
1494 dd
->rank
, dd
->ci
[XX
], dd
->ci
[YY
], dd
->ci
[ZZ
]);
1498 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
1499 dd
->rank
, dd
->ci
[XX
], dd
->ci
[YY
], dd
->ci
[ZZ
]);
1503 static void receive_ddindex2simnodeid(gmx_domdec_t
*dd
,
1507 gmx_domdec_comm_t
*comm
= dd
->comm
;
1509 if (!comm
->bCartesianPP_PME
&& comm
->bCartesianPP
)
1512 snew(comm
->ddindex2simnodeid
, dd
->nnodes
);
1513 snew(buf
, dd
->nnodes
);
1514 if (thisRankHasDuty(cr
, DUTY_PP
))
1516 buf
[dd_index(dd
->nc
, dd
->ci
)] = cr
->sim_nodeid
;
1518 /* Communicate the ddindex to simulation nodeid index */
1519 MPI_Allreduce(buf
, comm
->ddindex2simnodeid
, dd
->nnodes
, MPI_INT
, MPI_SUM
,
1520 cr
->mpi_comm_mysim
);
1524 GMX_UNUSED_VALUE(dd
);
1525 GMX_UNUSED_VALUE(cr
);
1529 static void split_communicator(const gmx::MDLogger
&mdlog
,
1530 t_commrec
*cr
, gmx_domdec_t
*dd
,
1531 DdRankOrder gmx_unused rankOrder
,
1532 bool gmx_unused reorder
)
1534 gmx_domdec_comm_t
*comm
;
1543 if (comm
->bCartesianPP
)
1545 for (i
= 1; i
< DIM
; i
++)
1547 bDiv
[i
] = ((cr
->npmenodes
*dd
->nc
[i
]) % (dd
->nnodes
) == 0);
1549 if (bDiv
[YY
] || bDiv
[ZZ
])
1551 comm
->bCartesianPP_PME
= TRUE
;
1552 /* If we have 2D PME decomposition, which is always in x+y,
1553 * we stack the PME only nodes in z.
1554 * Otherwise we choose the direction that provides the thinnest slab
1555 * of PME only nodes as this will have the least effect
1556 * on the PP communication.
1557 * But for the PME communication the opposite might be better.
1559 if (bDiv
[ZZ
] && (comm
->npmenodes_y
> 1 ||
1561 dd
->nc
[YY
] > dd
->nc
[ZZ
]))
1563 comm
->cartpmedim
= ZZ
;
1567 comm
->cartpmedim
= YY
;
1569 comm
->ntot
[comm
->cartpmedim
]
1570 += (cr
->npmenodes
*dd
->nc
[comm
->cartpmedim
])/dd
->nnodes
;
1574 GMX_LOG(mdlog
.info
).appendTextFormatted(
1575 "Number of PME-only ranks (%d) is not a multiple of nx*ny (%d*%d) or nx*nz (%d*%d)",
1576 cr
->npmenodes
, dd
->nc
[XX
], dd
->nc
[YY
], dd
->nc
[XX
], dd
->nc
[ZZ
]);
1577 GMX_LOG(mdlog
.info
).appendText("Will not use a Cartesian communicator for PP <-> PME\n");
1581 if (comm
->bCartesianPP_PME
)
1587 GMX_LOG(mdlog
.info
).appendTextFormatted(
1588 "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d",
1589 comm
->ntot
[XX
], comm
->ntot
[YY
], comm
->ntot
[ZZ
]);
1591 for (i
= 0; i
< DIM
; i
++)
1595 MPI_Cart_create(cr
->mpi_comm_mysim
, DIM
, comm
->ntot
, periods
, static_cast<int>(reorder
),
1597 MPI_Comm_rank(comm_cart
, &rank
);
1598 if (MASTER(cr
) && rank
!= 0)
1600 gmx_fatal(FARGS
, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
1603 /* With this assigment we loose the link to the original communicator
1604 * which will usually be MPI_COMM_WORLD, unless have multisim.
1606 cr
->mpi_comm_mysim
= comm_cart
;
1607 cr
->sim_nodeid
= rank
;
1609 MPI_Cart_coords(cr
->mpi_comm_mysim
, cr
->sim_nodeid
, DIM
, dd
->ci
);
1611 GMX_LOG(mdlog
.info
).appendTextFormatted(
1612 "Cartesian rank %d, coordinates %d %d %d\n",
1613 cr
->sim_nodeid
, dd
->ci
[XX
], dd
->ci
[YY
], dd
->ci
[ZZ
]);
1615 if (dd
->ci
[comm
->cartpmedim
] < dd
->nc
[comm
->cartpmedim
])
1619 if (cr
->npmenodes
== 0 ||
1620 dd
->ci
[comm
->cartpmedim
] >= dd
->nc
[comm
->cartpmedim
])
1622 cr
->duty
= DUTY_PME
;
1625 /* Split the sim communicator into PP and PME only nodes */
1626 MPI_Comm_split(cr
->mpi_comm_mysim
,
1627 getThisRankDuties(cr
),
1628 dd_index(comm
->ntot
, dd
->ci
),
1629 &cr
->mpi_comm_mygroup
);
1636 case DdRankOrder::pp_pme
:
1637 GMX_LOG(mdlog
.info
).appendText("Order of the ranks: PP first, PME last");
1639 case DdRankOrder::interleave
:
1640 /* Interleave the PP-only and PME-only ranks */
1641 GMX_LOG(mdlog
.info
).appendText("Interleaving PP and PME ranks");
1642 comm
->pmenodes
= dd_interleaved_pme_ranks(dd
);
1644 case DdRankOrder::cartesian
:
1647 gmx_fatal(FARGS
, "Invalid ddRankOrder=%d", static_cast<int>(rankOrder
));
1650 if (dd_simnode2pmenode(dd
, cr
, cr
->sim_nodeid
) == -1)
1652 cr
->duty
= DUTY_PME
;
1659 /* Split the sim communicator into PP and PME only nodes */
1660 MPI_Comm_split(cr
->mpi_comm_mysim
,
1661 getThisRankDuties(cr
),
1663 &cr
->mpi_comm_mygroup
);
1664 MPI_Comm_rank(cr
->mpi_comm_mygroup
, &cr
->nodeid
);
1668 GMX_LOG(mdlog
.info
).appendTextFormatted(
1669 "This rank does only %s work.\n",
1670 thisRankHasDuty(cr
, DUTY_PP
) ? "particle-particle" : "PME-mesh");
1673 /*! \brief Generates the MPI communicators for domain decomposition */
1674 static void make_dd_communicators(const gmx::MDLogger
&mdlog
,
1676 gmx_domdec_t
*dd
, DdRankOrder ddRankOrder
)
1678 gmx_domdec_comm_t
*comm
;
1683 copy_ivec(dd
->nc
, comm
->ntot
);
1685 comm
->bCartesianPP
= (ddRankOrder
== DdRankOrder::cartesian
);
1686 comm
->bCartesianPP_PME
= FALSE
;
1688 /* Reorder the nodes by default. This might change the MPI ranks.
1689 * Real reordering is only supported on very few architectures,
1690 * Blue Gene is one of them.
1692 CartReorder
= getenv("GMX_NO_CART_REORDER") == nullptr;
1694 if (cr
->npmenodes
> 0)
1696 /* Split the communicator into a PP and PME part */
1697 split_communicator(mdlog
, cr
, dd
, ddRankOrder
, CartReorder
);
1698 if (comm
->bCartesianPP_PME
)
1700 /* We (possibly) reordered the nodes in split_communicator,
1701 * so it is no longer required in make_pp_communicator.
1703 CartReorder
= FALSE
;
1708 /* All nodes do PP and PME */
1710 /* We do not require separate communicators */
1711 cr
->mpi_comm_mygroup
= cr
->mpi_comm_mysim
;
1715 if (thisRankHasDuty(cr
, DUTY_PP
))
1717 /* Copy or make a new PP communicator */
1718 make_pp_communicator(mdlog
, dd
, cr
, CartReorder
);
1722 receive_ddindex2simnodeid(dd
, cr
);
1725 if (!thisRankHasDuty(cr
, DUTY_PME
))
1727 /* Set up the commnuication to our PME node */
1728 dd
->pme_nodeid
= dd_simnode2pmenode(dd
, cr
, cr
->sim_nodeid
);
1729 dd
->pme_receive_vir_ener
= receive_vir_ener(dd
, cr
);
1732 fprintf(debug
, "My pme_nodeid %d receive ener %s\n",
1733 dd
->pme_nodeid
, gmx::boolToString(dd
->pme_receive_vir_ener
));
1738 dd
->pme_nodeid
= -1;
1741 /* We can not use DDMASTER(dd), because dd->masterrank is set later */
1744 dd
->ma
= std::make_unique
<AtomDistribution
>(dd
->nc
,
1746 comm
->cgs_gl
.index
[comm
->cgs_gl
.nr
]);
1750 static real
*get_slb_frac(const gmx::MDLogger
&mdlog
,
1751 const char *dir
, int nc
, const char *size_string
)
1753 real
*slb_frac
, tot
;
1758 if (nc
> 1 && size_string
!= nullptr)
1760 GMX_LOG(mdlog
.info
).appendTextFormatted(
1761 "Using static load balancing for the %s direction", dir
);
1764 for (i
= 0; i
< nc
; i
++)
1767 sscanf(size_string
, "%20lf%n", &dbl
, &n
);
1770 gmx_fatal(FARGS
, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir
, size_string
);
1777 std::string relativeCellSizes
= "Relative cell sizes:";
1778 for (i
= 0; i
< nc
; i
++)
1781 relativeCellSizes
+= gmx::formatString(" %5.3f", slb_frac
[i
]);
1783 GMX_LOG(mdlog
.info
).appendText(relativeCellSizes
);
1789 static int multi_body_bondeds_count(const gmx_mtop_t
*mtop
)
1792 gmx_mtop_ilistloop_t iloop
= gmx_mtop_ilistloop_init(mtop
);
1794 while (const InteractionLists
*ilists
= gmx_mtop_ilistloop_next(iloop
, &nmol
))
1796 for (auto &ilist
: extractILists(*ilists
, IF_BOND
))
1798 if (NRAL(ilist
.functionType
) > 2)
1800 n
+= nmol
*(ilist
.iatoms
.size()/ilistStride(ilist
));
1808 static int dd_getenv(const gmx::MDLogger
&mdlog
,
1809 const char *env_var
, int def
)
1815 val
= getenv(env_var
);
1818 if (sscanf(val
, "%20d", &nst
) <= 0)
1822 GMX_LOG(mdlog
.info
).appendTextFormatted(
1823 "Found env.var. %s = %s, using value %d",
1830 static void check_dd_restrictions(const gmx_domdec_t
*dd
,
1831 const t_inputrec
*ir
,
1832 const gmx::MDLogger
&mdlog
)
1834 if (ir
->ePBC
== epbcSCREW
&&
1835 (dd
->nc
[XX
] == 1 || dd
->nc
[YY
] > 1 || dd
->nc
[ZZ
] > 1))
1837 gmx_fatal(FARGS
, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names
[ir
->ePBC
]);
1840 if (ir
->ns_type
== ensSIMPLE
)
1842 gmx_fatal(FARGS
, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
1845 if (ir
->nstlist
== 0)
1847 gmx_fatal(FARGS
, "Domain decomposition does not work with nstlist=0");
1850 if (ir
->comm_mode
== ecmANGULAR
&& ir
->ePBC
!= epbcNONE
)
1852 GMX_LOG(mdlog
.warning
).appendText("comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
1856 static real
average_cellsize_min(gmx_domdec_t
*dd
, gmx_ddbox_t
*ddbox
)
1861 r
= ddbox
->box_size
[XX
];
1862 for (di
= 0; di
< dd
->ndim
; di
++)
1865 /* Check using the initial average cell size */
1866 r
= std::min(r
, ddbox
->box_size
[d
]*ddbox
->skew_fac
[d
]/dd
->nc
[d
]);
1872 /*! \brief Depending on the DLB initial value return the DLB switched off state or issue an error.
1874 static DlbState
forceDlbOffOrBail(DlbState cmdlineDlbState
,
1875 const std::string
&reasonStr
,
1876 const gmx::MDLogger
&mdlog
)
1878 std::string dlbNotSupportedErr
= "Dynamic load balancing requested, but ";
1879 std::string dlbDisableNote
= "NOTE: disabling dynamic load balancing as ";
1881 if (cmdlineDlbState
== DlbState::onUser
)
1883 gmx_fatal(FARGS
, "%s", (dlbNotSupportedErr
+ reasonStr
).c_str());
1885 else if (cmdlineDlbState
== DlbState::offCanTurnOn
)
1887 GMX_LOG(mdlog
.info
).appendText(dlbDisableNote
+ reasonStr
);
1889 return DlbState::offForever
;
1892 /*! \brief Return the dynamic load balancer's initial state based on initial conditions and user inputs.
1894 * This function parses the parameters of "-dlb" command line option setting
1895 * corresponding state values. Then it checks the consistency of the determined
1896 * state with other run parameters and settings. As a result, the initial state
1897 * may be altered or an error may be thrown if incompatibility of options is detected.
1899 * \param [in] mdlog Logger.
1900 * \param [in] dlbOption Enum value for the DLB option.
1901 * \param [in] bRecordLoad True if the load balancer is recording load information.
1902 * \param [in] mdrunOptions Options for mdrun.
1903 * \param [in] ir Pointer mdrun to input parameters.
1904 * \returns DLB initial/startup state.
1906 static DlbState
determineInitialDlbState(const gmx::MDLogger
&mdlog
,
1907 DlbOption dlbOption
, gmx_bool bRecordLoad
,
1908 const gmx::MdrunOptions
&mdrunOptions
,
1909 const t_inputrec
*ir
)
1911 DlbState dlbState
= DlbState::offCanTurnOn
;
1915 case DlbOption::turnOnWhenUseful
: dlbState
= DlbState::offCanTurnOn
; break;
1916 case DlbOption::no
: dlbState
= DlbState::offUser
; break;
1917 case DlbOption::yes
: dlbState
= DlbState::onUser
; break;
1918 default: gmx_incons("Invalid dlbOption enum value");
1921 /* Reruns don't support DLB: bail or override auto mode */
1922 if (mdrunOptions
.rerun
)
1924 std::string reasonStr
= "it is not supported in reruns.";
1925 return forceDlbOffOrBail(dlbState
, reasonStr
, mdlog
);
1928 /* Unsupported integrators */
1929 if (!EI_DYNAMICS(ir
->eI
))
1931 auto reasonStr
= gmx::formatString("it is only supported with dynamics, not with integrator '%s'.", EI(ir
->eI
));
1932 return forceDlbOffOrBail(dlbState
, reasonStr
, mdlog
);
1935 /* Without cycle counters we can't time work to balance on */
1938 std::string reasonStr
= "cycle counters unsupported or not enabled in the operating system kernel.";
1939 return forceDlbOffOrBail(dlbState
, reasonStr
, mdlog
);
1942 if (mdrunOptions
.reproducible
)
1944 std::string reasonStr
= "you started a reproducible run.";
1947 case DlbState::offUser
:
1949 case DlbState::offForever
:
1950 GMX_RELEASE_ASSERT(false, "DlbState::offForever is not a valid initial state");
1952 case DlbState::offCanTurnOn
:
1953 return forceDlbOffOrBail(dlbState
, reasonStr
, mdlog
);
1954 case DlbState::onCanTurnOff
:
1955 GMX_RELEASE_ASSERT(false, "DlbState::offCanTurnOff is not a valid initial state");
1957 case DlbState::onUser
:
1958 return forceDlbOffOrBail(dlbState
, reasonStr
+ " In load balanced runs binary reproducibility cannot be ensured.", mdlog
);
1960 gmx_fatal(FARGS
, "Death horror: undefined case (%d) for load balancing choice", static_cast<int>(dlbState
));
1967 static void set_dd_dim(const gmx::MDLogger
&mdlog
, gmx_domdec_t
*dd
)
1970 if (getenv("GMX_DD_ORDER_ZYX") != nullptr)
1972 /* Decomposition order z,y,x */
1973 GMX_LOG(mdlog
.info
).appendText("Using domain decomposition order z, y, x");
1974 for (int dim
= DIM
-1; dim
>= 0; dim
--)
1976 if (dd
->nc
[dim
] > 1)
1978 dd
->dim
[dd
->ndim
++] = dim
;
1984 /* Decomposition order x,y,z */
1985 for (int dim
= 0; dim
< DIM
; dim
++)
1987 if (dd
->nc
[dim
] > 1)
1989 dd
->dim
[dd
->ndim
++] = dim
;
1996 /* Set dim[0] to avoid extra checks on ndim in several places */
2001 static gmx_domdec_comm_t
*init_dd_comm()
2003 gmx_domdec_comm_t
*comm
= new gmx_domdec_comm_t
;
2005 comm
->n_load_have
= 0;
2006 comm
->n_load_collect
= 0;
2008 comm
->haveTurnedOffDlb
= false;
2010 for (int i
= 0; i
< static_cast<int>(DDAtomRanges::Type::Number
); i
++)
2012 comm
->sum_nat
[i
] = 0;
2016 comm
->load_step
= 0;
2019 clear_ivec(comm
->load_lim
);
2023 /* This should be replaced by a unique pointer */
2024 comm
->balanceRegion
= ddBalanceRegionAllocate();
2029 /* Returns whether mtop contains constraints and/or vsites */
2030 static bool systemHasConstraintsOrVsites(const gmx_mtop_t
&mtop
)
2032 auto ilistLoop
= gmx_mtop_ilistloop_init(mtop
);
2034 while (const InteractionLists
*ilists
= gmx_mtop_ilistloop_next(ilistLoop
, &nmol
))
2036 if (!extractILists(*ilists
, IF_CONSTRAINT
| IF_VSITE
).empty())
2045 static void setupUpdateGroups(const gmx::MDLogger
&mdlog
,
2046 const gmx_mtop_t
&mtop
,
2047 const t_inputrec
&inputrec
,
2049 int numMpiRanksTotal
,
2050 gmx_domdec_comm_t
*comm
)
2052 /* When we have constraints and/or vsites, it is beneficial to use
2053 * update groups (when possible) to allow independent update of groups.
2055 if (!systemHasConstraintsOrVsites(mtop
))
2057 /* No constraints or vsites, atoms can be updated independently */
2061 comm
->updateGroupingPerMoleculetype
= gmx::makeUpdateGroups(mtop
);
2062 comm
->useUpdateGroups
=
2063 (!comm
->updateGroupingPerMoleculetype
.empty() &&
2064 getenv("GMX_NO_UPDATEGROUPS") == nullptr);
2066 if (comm
->useUpdateGroups
)
2068 int numUpdateGroups
= 0;
2069 for (const auto &molblock
: mtop
.molblock
)
2071 numUpdateGroups
+= molblock
.nmol
*comm
->updateGroupingPerMoleculetype
[molblock
.type
].numBlocks();
2074 /* Note: We would like to use dd->nnodes for the atom count estimate,
2075 * but that is not yet available here. But this anyhow only
2076 * affect performance up to the second dd_partition_system call.
2078 int homeAtomCountEstimate
= mtop
.natoms
/numMpiRanksTotal
;
2079 comm
->updateGroupsCog
=
2080 std::make_unique
<gmx::UpdateGroupsCog
>(mtop
,
2081 comm
->updateGroupingPerMoleculetype
,
2082 maxReferenceTemperature(inputrec
),
2083 homeAtomCountEstimate
);
2085 /* To use update groups, the large domain-to-domain cutoff distance
2086 * should be compatible with the box size.
2088 comm
->useUpdateGroups
= (atomToAtomIntoDomainToDomainCutoff(*comm
, 0) < cutoffMargin
);
2090 if (comm
->useUpdateGroups
)
2092 GMX_LOG(mdlog
.info
).appendTextFormatted(
2093 "Using update groups, nr %d, average size %.1f atoms, max. radius %.3f nm\n",
2095 mtop
.natoms
/static_cast<double>(numUpdateGroups
),
2096 comm
->updateGroupsCog
->maxUpdateGroupRadius());
2100 GMX_LOG(mdlog
.info
).appendTextFormatted("The combination of rlist and box size prohibits the use of update groups\n");
2101 comm
->updateGroupingPerMoleculetype
.clear();
2102 comm
->updateGroupsCog
.reset(nullptr);
2107 UnitCellInfo::UnitCellInfo(const t_inputrec
&ir
) :
2108 npbcdim(ePBC2npbcdim(ir
.ePBC
)),
2109 numBoundedDimensions(inputrec2nboundeddim(&ir
)),
2110 ddBoxIsDynamic(numBoundedDimensions
< DIM
|| inputrecDynamicBox(&ir
)),
2111 haveScrewPBC(ir
.ePBC
== epbcSCREW
)
2115 /*! \brief Set the cell size and interaction limits, as well as the DD grid */
2116 static void set_dd_limits_and_grid(const gmx::MDLogger
&mdlog
,
2117 t_commrec
*cr
, gmx_domdec_t
*dd
,
2118 const DomdecOptions
&options
,
2119 const gmx::MdrunOptions
&mdrunOptions
,
2120 const gmx_mtop_t
*mtop
,
2121 const t_inputrec
*ir
,
2123 gmx::ArrayRef
<const gmx::RVec
> xGlobal
,
2127 real r_bonded_limit
= -1;
2128 const real tenPercentMargin
= 1.1;
2129 gmx_domdec_comm_t
*comm
= dd
->comm
;
2131 dd
->pme_recv_f_alloc
= 0;
2132 dd
->pme_recv_f_buf
= nullptr;
2134 /* Initialize to GPU share count to 0, might change later */
2135 comm
->nrank_gpu_shared
= 0;
2137 comm
->dlbState
= determineInitialDlbState(mdlog
, options
.dlbOption
, comm
->bRecordLoad
, mdrunOptions
, ir
);
2138 dd_dlb_set_should_check_whether_to_turn_dlb_on(dd
, TRUE
);
2139 /* To consider turning DLB on after 2*nstlist steps we need to check
2140 * at partitioning count 3. Thus we need to increase the first count by 2.
2142 comm
->ddPartioningCountFirstDlbOff
+= 2;
2144 GMX_LOG(mdlog
.info
).appendTextFormatted(
2145 "Dynamic load balancing: %s", edlbs_names
[int(comm
->dlbState
)]);
2147 comm
->bPMELoadBalDLBLimits
= FALSE
;
2149 /* Allocate the charge group/atom sorting struct */
2150 comm
->sort
= std::make_unique
<gmx_domdec_sort_t
>();
2152 /* We need to decide on update groups early, as this affects communication distances */
2153 comm
->useUpdateGroups
= false;
2154 if (ir
->cutoff_scheme
== ecutsVERLET
)
2156 real cutoffMargin
= std::sqrt(max_cutoff2(ir
->ePBC
, box
)) - ir
->rlist
;
2157 setupUpdateGroups(mdlog
, *mtop
, *ir
, cutoffMargin
, cr
->nnodes
, comm
);
2160 // TODO: Check whether all bondeds are within update groups
2161 comm
->haveInterDomainBondeds
= (mtop
->natoms
> gmx_mtop_num_molecules(*mtop
) ||
2162 mtop
->bIntermolecularInteractions
);
2163 comm
->haveInterDomainMultiBodyBondeds
= (multi_body_bondeds_count(mtop
) > 0);
2165 if (comm
->useUpdateGroups
)
2167 dd
->splitConstraints
= false;
2168 dd
->splitSettles
= false;
2172 dd
->splitConstraints
= gmx::inter_charge_group_constraints(*mtop
);
2173 dd
->splitSettles
= gmx::inter_charge_group_settles(*mtop
);
2178 /* Set the cut-off to some very large value,
2179 * so we don't need if statements everywhere in the code.
2180 * We use sqrt, since the cut-off is squared in some places.
2182 comm
->cutoff
= GMX_CUTOFF_INF
;
2186 comm
->cutoff
= atomToAtomIntoDomainToDomainCutoff(*comm
, ir
->rlist
);
2188 comm
->cutoff_mbody
= 0;
2190 /* Determine the minimum cell size limit, affected by many factors */
2191 comm
->cellsize_limit
= 0;
2192 comm
->bBondComm
= FALSE
;
2194 /* We do not allow home atoms to move beyond the neighboring domain
2195 * between domain decomposition steps, which limits the cell size.
2196 * Get an estimate of cell size limit due to atom displacement.
2197 * In most cases this is a large overestimate, because it assumes
2198 * non-interaction atoms.
2199 * We set the chance to 1 in a trillion steps.
2201 constexpr real c_chanceThatAtomMovesBeyondDomain
= 1e-12;
2202 const real limitForAtomDisplacement
=
2203 minCellSizeForAtomDisplacement(*mtop
, *ir
,
2204 comm
->updateGroupingPerMoleculetype
,
2205 c_chanceThatAtomMovesBeyondDomain
);
2206 GMX_LOG(mdlog
.info
).appendTextFormatted(
2207 "Minimum cell size due to atom displacement: %.3f nm",
2208 limitForAtomDisplacement
);
2210 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
,
2211 limitForAtomDisplacement
);
2213 /* TODO: PME decomposition currently requires atoms not to be more than
2214 * 2/3 of comm->cutoff, which is >=rlist, outside of their domain.
2215 * In nearly all cases, limitForAtomDisplacement will be smaller
2216 * than 2/3*rlist, so the PME requirement is satisfied.
2217 * But it would be better for both correctness and performance
2218 * to use limitForAtomDisplacement instead of 2/3*comm->cutoff.
2219 * Note that we would need to improve the pairlist buffer case.
2222 if (comm
->haveInterDomainBondeds
)
2224 if (options
.minimumCommunicationRange
> 0)
2226 comm
->cutoff_mbody
= atomToAtomIntoDomainToDomainCutoff(*comm
, options
.minimumCommunicationRange
);
2227 if (options
.useBondedCommunication
)
2229 comm
->bBondComm
= (comm
->cutoff_mbody
> comm
->cutoff
);
2233 comm
->cutoff
= std::max(comm
->cutoff
, comm
->cutoff_mbody
);
2235 r_bonded_limit
= comm
->cutoff_mbody
;
2237 else if (ir
->bPeriodicMols
)
2239 /* Can not easily determine the required cut-off */
2240 GMX_LOG(mdlog
.warning
).appendText("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.");
2241 comm
->cutoff_mbody
= comm
->cutoff
/2;
2242 r_bonded_limit
= comm
->cutoff_mbody
;
2250 dd_bonded_cg_distance(mdlog
, mtop
, ir
, as_rvec_array(xGlobal
.data()), box
,
2251 options
.checkBondedInteractions
,
2254 gmx_bcast(sizeof(r_2b
), &r_2b
, cr
);
2255 gmx_bcast(sizeof(r_mb
), &r_mb
, cr
);
2257 /* We use an initial margin of 10% for the minimum cell size,
2258 * except when we are just below the non-bonded cut-off.
2260 if (options
.useBondedCommunication
)
2262 if (std::max(r_2b
, r_mb
) > comm
->cutoff
)
2264 r_bonded
= std::max(r_2b
, r_mb
);
2265 r_bonded_limit
= tenPercentMargin
*r_bonded
;
2266 comm
->bBondComm
= TRUE
;
2271 r_bonded_limit
= std::min(tenPercentMargin
*r_bonded
, comm
->cutoff
);
2273 /* We determine cutoff_mbody later */
2277 /* No special bonded communication,
2278 * simply increase the DD cut-off.
2280 r_bonded_limit
= tenPercentMargin
*std::max(r_2b
, r_mb
);
2281 comm
->cutoff_mbody
= r_bonded_limit
;
2282 comm
->cutoff
= std::max(comm
->cutoff
, comm
->cutoff_mbody
);
2285 GMX_LOG(mdlog
.info
).appendTextFormatted(
2286 "Minimum cell size due to bonded interactions: %.3f nm",
2289 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
, r_bonded_limit
);
2293 if (dd
->splitConstraints
&& options
.constraintCommunicationRange
<= 0)
2295 /* There is a cell size limit due to the constraints (P-LINCS) */
2296 rconstr
= gmx::constr_r_max(mdlog
, mtop
, ir
);
2297 GMX_LOG(mdlog
.info
).appendTextFormatted(
2298 "Estimated maximum distance required for P-LINCS: %.3f nm",
2300 if (rconstr
> comm
->cellsize_limit
)
2302 GMX_LOG(mdlog
.info
).appendText("This distance will limit the DD cell size, you can override this with -rcon");
2305 else if (options
.constraintCommunicationRange
> 0)
2307 /* Here we do not check for dd->splitConstraints.
2308 * because one can also set a cell size limit for virtual sites only
2309 * and at this point we don't know yet if there are intercg v-sites.
2311 GMX_LOG(mdlog
.info
).appendTextFormatted(
2312 "User supplied maximum distance required for P-LINCS: %.3f nm",
2313 options
.constraintCommunicationRange
);
2314 rconstr
= options
.constraintCommunicationRange
;
2316 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
, rconstr
);
2318 comm
->cgs_gl
= gmx_mtop_global_cgs(mtop
);
2320 if (options
.numCells
[XX
] > 0)
2322 copy_ivec(options
.numCells
, dd
->nc
);
2323 set_dd_dim(mdlog
, dd
);
2324 set_ddbox_cr(*cr
, &dd
->nc
, *ir
, box
, xGlobal
, ddbox
);
2326 if (options
.numPmeRanks
>= 0)
2328 cr
->npmenodes
= options
.numPmeRanks
;
2332 /* When the DD grid is set explicitly and -npme is set to auto,
2333 * don't use PME ranks. We check later if the DD grid is
2334 * compatible with the total number of ranks.
2339 real acs
= average_cellsize_min(dd
, ddbox
);
2340 if (acs
< comm
->cellsize_limit
)
2342 gmx_fatal_collective(FARGS
, cr
->mpi_comm_mysim
, MASTER(cr
),
2343 "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",
2344 acs
, comm
->cellsize_limit
);
2349 set_ddbox_cr(*cr
, nullptr, *ir
, box
, xGlobal
, ddbox
);
2351 /* We need to choose the optimal DD grid and possibly PME nodes */
2353 dd_choose_grid(mdlog
, cr
, dd
, ir
, mtop
, box
, ddbox
,
2354 options
.numPmeRanks
,
2355 !isDlbDisabled(comm
),
2357 comm
->cellsize_limit
, comm
->cutoff
,
2358 comm
->haveInterDomainBondeds
);
2360 if (dd
->nc
[XX
] == 0)
2363 gmx_bool bC
= (dd
->splitConstraints
&& rconstr
> r_bonded_limit
);
2364 sprintf(buf
, "Change the number of ranks or mdrun option %s%s%s",
2365 !bC
? "-rdd" : "-rcon",
2366 comm
->dlbState
!= DlbState::offUser
? " or -dds" : "",
2367 bC
? " or your LINCS settings" : "");
2369 gmx_fatal_collective(FARGS
, cr
->mpi_comm_mysim
, MASTER(cr
),
2370 "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
2372 "Look in the log file for details on the domain decomposition",
2373 cr
->nnodes
-cr
->npmenodes
, limit
, buf
);
2375 set_dd_dim(mdlog
, dd
);
2378 GMX_LOG(mdlog
.info
).appendTextFormatted(
2379 "Domain decomposition grid %d x %d x %d, separate PME ranks %d",
2380 dd
->nc
[XX
], dd
->nc
[YY
], dd
->nc
[ZZ
], cr
->npmenodes
);
2382 dd
->nnodes
= dd
->nc
[XX
]*dd
->nc
[YY
]*dd
->nc
[ZZ
];
2383 if (cr
->nnodes
- dd
->nnodes
!= cr
->npmenodes
)
2385 gmx_fatal_collective(FARGS
, cr
->mpi_comm_mysim
, MASTER(cr
),
2386 "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
2387 dd
->nnodes
, cr
->nnodes
- cr
->npmenodes
, cr
->nnodes
);
2389 if (cr
->npmenodes
> dd
->nnodes
)
2391 gmx_fatal_collective(FARGS
, cr
->mpi_comm_mysim
, MASTER(cr
),
2392 "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr
->npmenodes
, dd
->nnodes
);
2394 if (cr
->npmenodes
> 0)
2396 comm
->npmenodes
= cr
->npmenodes
;
2400 comm
->npmenodes
= dd
->nnodes
;
2403 if (EEL_PME(ir
->coulombtype
) || EVDW_PME(ir
->vdwtype
))
2405 /* The following choices should match those
2406 * in comm_cost_est in domdec_setup.c.
2407 * Note that here the checks have to take into account
2408 * that the decomposition might occur in a different order than xyz
2409 * (for instance through the env.var. GMX_DD_ORDER_ZYX),
2410 * in which case they will not match those in comm_cost_est,
2411 * but since that is mainly for testing purposes that's fine.
2413 if (dd
->ndim
>= 2 && dd
->dim
[0] == XX
&& dd
->dim
[1] == YY
&&
2414 comm
->npmenodes
> dd
->nc
[XX
] && comm
->npmenodes
% dd
->nc
[XX
] == 0 &&
2415 getenv("GMX_PMEONEDD") == nullptr)
2417 comm
->npmedecompdim
= 2;
2418 comm
->npmenodes_x
= dd
->nc
[XX
];
2419 comm
->npmenodes_y
= comm
->npmenodes
/comm
->npmenodes_x
;
2423 /* In case nc is 1 in both x and y we could still choose to
2424 * decompose pme in y instead of x, but we use x for simplicity.
2426 comm
->npmedecompdim
= 1;
2427 if (dd
->dim
[0] == YY
)
2429 comm
->npmenodes_x
= 1;
2430 comm
->npmenodes_y
= comm
->npmenodes
;
2434 comm
->npmenodes_x
= comm
->npmenodes
;
2435 comm
->npmenodes_y
= 1;
2438 GMX_LOG(mdlog
.info
).appendTextFormatted(
2439 "PME domain decomposition: %d x %d x %d",
2440 comm
->npmenodes_x
, comm
->npmenodes_y
, 1);
2444 comm
->npmedecompdim
= 0;
2445 comm
->npmenodes_x
= 0;
2446 comm
->npmenodes_y
= 0;
2449 snew(comm
->slb_frac
, DIM
);
2450 if (isDlbDisabled(comm
))
2452 comm
->slb_frac
[XX
] = get_slb_frac(mdlog
, "x", dd
->nc
[XX
], options
.cellSizeX
);
2453 comm
->slb_frac
[YY
] = get_slb_frac(mdlog
, "y", dd
->nc
[YY
], options
.cellSizeY
);
2454 comm
->slb_frac
[ZZ
] = get_slb_frac(mdlog
, "z", dd
->nc
[ZZ
], options
.cellSizeZ
);
2457 if (comm
->haveInterDomainBondeds
&& comm
->cutoff_mbody
== 0)
2459 if (comm
->bBondComm
|| !isDlbDisabled(comm
))
2461 /* Set the bonded communication distance to halfway
2462 * the minimum and the maximum,
2463 * since the extra communication cost is nearly zero.
2465 real acs
= average_cellsize_min(dd
, ddbox
);
2466 comm
->cutoff_mbody
= 0.5*(r_bonded
+ acs
);
2467 if (!isDlbDisabled(comm
))
2469 /* Check if this does not limit the scaling */
2470 comm
->cutoff_mbody
= std::min(comm
->cutoff_mbody
,
2471 options
.dlbScaling
*acs
);
2473 if (!comm
->bBondComm
)
2475 /* Without bBondComm do not go beyond the n.b. cut-off */
2476 comm
->cutoff_mbody
= std::min(comm
->cutoff_mbody
, comm
->cutoff
);
2477 if (comm
->cellsize_limit
>= comm
->cutoff
)
2479 /* We don't loose a lot of efficieny
2480 * when increasing it to the n.b. cut-off.
2481 * It can even be slightly faster, because we need
2482 * less checks for the communication setup.
2484 comm
->cutoff_mbody
= comm
->cutoff
;
2487 /* Check if we did not end up below our original limit */
2488 comm
->cutoff_mbody
= std::max(comm
->cutoff_mbody
, r_bonded_limit
);
2490 if (comm
->cutoff_mbody
> comm
->cellsize_limit
)
2492 comm
->cellsize_limit
= comm
->cutoff_mbody
;
2495 /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
2500 fprintf(debug
, "Bonded atom communication beyond the cut-off: %s\n"
2501 "cellsize limit %f\n",
2502 gmx::boolToString(comm
->bBondComm
), comm
->cellsize_limit
);
2507 check_dd_restrictions(dd
, ir
, mdlog
);
2511 static char *init_bLocalCG(const gmx_mtop_t
*mtop
)
2516 ncg
= ncg_mtop(mtop
);
2517 snew(bLocalCG
, ncg
);
2518 for (cg
= 0; cg
< ncg
; cg
++)
2520 bLocalCG
[cg
] = FALSE
;
2526 void dd_init_bondeds(FILE *fplog
,
2528 const gmx_mtop_t
*mtop
,
2529 const gmx_vsite_t
*vsite
,
2530 const t_inputrec
*ir
,
2531 gmx_bool bBCheck
, cginfo_mb_t
*cginfo_mb
)
2533 gmx_domdec_comm_t
*comm
;
2535 dd_make_reverse_top(fplog
, dd
, mtop
, vsite
, ir
, bBCheck
);
2539 if (comm
->bBondComm
)
2541 /* Communicate atoms beyond the cut-off for bonded interactions */
2544 comm
->cglink
= make_charge_group_links(mtop
, dd
, cginfo_mb
);
2546 comm
->bLocalCG
= init_bLocalCG(mtop
);
2550 /* Only communicate atoms based on cut-off */
2551 comm
->cglink
= nullptr;
2552 comm
->bLocalCG
= nullptr;
2556 static void writeSettings(gmx::TextWriter
*log
,
2558 const gmx_mtop_t
*mtop
,
2559 const t_inputrec
*ir
,
2560 gmx_bool bDynLoadBal
,
2562 const gmx_ddbox_t
*ddbox
)
2564 gmx_domdec_comm_t
*comm
;
2573 log
->writeString("The maximum number of communication pulses is:");
2574 for (d
= 0; d
< dd
->ndim
; d
++)
2576 log
->writeStringFormatted(" %c %d", dim2char(dd
->dim
[d
]), comm
->cd
[d
].np_dlb
);
2578 log
->ensureLineBreak();
2579 log
->writeLineFormatted("The minimum size for domain decomposition cells is %.3f nm", comm
->cellsize_limit
);
2580 log
->writeLineFormatted("The requested allowed shrink of DD cells (option -dds) is: %.2f", dlb_scale
);
2581 log
->writeString("The allowed shrink of domain decomposition cells is:");
2582 for (d
= 0; d
< DIM
; d
++)
2586 if (d
>= ddbox
->npbcdim
&& dd
->nc
[d
] == 2)
2593 comm
->cellsize_min_dlb
[d
]/
2594 (ddbox
->box_size
[d
]*ddbox
->skew_fac
[d
]/dd
->nc
[d
]);
2596 log
->writeStringFormatted(" %c %.2f", dim2char(d
), shrink
);
2599 log
->ensureLineBreak();
2603 set_dd_cell_sizes_slb(dd
, ddbox
, setcellsizeslbPULSE_ONLY
, np
);
2604 log
->writeString("The initial number of communication pulses is:");
2605 for (d
= 0; d
< dd
->ndim
; d
++)
2607 log
->writeStringFormatted(" %c %d", dim2char(dd
->dim
[d
]), np
[dd
->dim
[d
]]);
2609 log
->ensureLineBreak();
2610 log
->writeString("The initial domain decomposition cell size is:");
2611 for (d
= 0; d
< DIM
; d
++)
2615 log
->writeStringFormatted(" %c %.2f nm",
2616 dim2char(d
), dd
->comm
->cellsize_min
[d
]);
2619 log
->ensureLineBreak();
2623 const bool haveInterDomainVsites
=
2624 (countInterUpdategroupVsites(*mtop
, comm
->updateGroupingPerMoleculetype
) != 0);
2626 if (comm
->haveInterDomainBondeds
||
2627 haveInterDomainVsites
||
2628 dd
->splitConstraints
|| dd
->splitSettles
)
2630 std::string decompUnits
;
2631 if (comm
->useUpdateGroups
)
2633 decompUnits
= "atom groups";
2637 decompUnits
= "atoms";
2640 log
->writeLineFormatted("The maximum allowed distance for %s involved in interactions is:", decompUnits
.c_str());
2641 log
->writeLineFormatted("%40s %-7s %6.3f nm", "non-bonded interactions", "", comm
->cutoff
);
2645 limit
= dd
->comm
->cellsize_limit
;
2649 if (dd
->unitCellInfo
.ddBoxIsDynamic
)
2651 log
->writeLine("(the following are initial values, they could change due to box deformation)");
2653 limit
= dd
->comm
->cellsize_min
[XX
];
2654 for (d
= 1; d
< DIM
; d
++)
2656 limit
= std::min(limit
, dd
->comm
->cellsize_min
[d
]);
2660 if (comm
->haveInterDomainBondeds
)
2662 log
->writeLineFormatted("%40s %-7s %6.3f nm",
2663 "two-body bonded interactions", "(-rdd)",
2664 std::max(comm
->cutoff
, comm
->cutoff_mbody
));
2665 log
->writeLineFormatted("%40s %-7s %6.3f nm",
2666 "multi-body bonded interactions", "(-rdd)",
2667 (comm
->bBondComm
|| isDlbOn(dd
->comm
)) ? comm
->cutoff_mbody
: std::min(comm
->cutoff
, limit
));
2669 if (haveInterDomainVsites
)
2671 log
->writeLineFormatted("%40s %-7s %6.3f nm",
2672 "virtual site constructions", "(-rcon)", limit
);
2674 if (dd
->splitConstraints
|| dd
->splitSettles
)
2676 std::string separation
= gmx::formatString("atoms separated by up to %d constraints",
2678 log
->writeLineFormatted("%40s %-7s %6.3f nm\n",
2679 separation
.c_str(), "(-rcon)", limit
);
2681 log
->ensureLineBreak();
2685 static void logSettings(const gmx::MDLogger
&mdlog
,
2687 const gmx_mtop_t
*mtop
,
2688 const t_inputrec
*ir
,
2690 const gmx_ddbox_t
*ddbox
)
2692 gmx::StringOutputStream stream
;
2693 gmx::TextWriter
log(&stream
);
2694 writeSettings(&log
, dd
, mtop
, ir
, isDlbOn(dd
->comm
), dlb_scale
, ddbox
);
2695 if (dd
->comm
->dlbState
== DlbState::offCanTurnOn
)
2698 log
.ensureEmptyLine();
2699 log
.writeLine("When dynamic load balancing gets turned on, these settings will change to:");
2701 writeSettings(&log
, dd
, mtop
, ir
, true, dlb_scale
, ddbox
);
2703 GMX_LOG(mdlog
.info
).asParagraph().appendText(stream
.toString());
2706 static void set_cell_limits_dlb(const gmx::MDLogger
&mdlog
,
2709 const t_inputrec
*ir
,
2710 const gmx_ddbox_t
*ddbox
)
2712 gmx_domdec_comm_t
*comm
;
2713 int d
, dim
, npulse
, npulse_d_max
, npulse_d
;
2718 bNoCutOff
= (ir
->rvdw
== 0 || ir
->rcoulomb
== 0);
2720 /* Determine the maximum number of comm. pulses in one dimension */
2722 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
, comm
->cutoff_mbody
);
2724 /* Determine the maximum required number of grid pulses */
2725 if (comm
->cellsize_limit
>= comm
->cutoff
)
2727 /* Only a single pulse is required */
2730 else if (!bNoCutOff
&& comm
->cellsize_limit
> 0)
2732 /* We round down slightly here to avoid overhead due to the latency
2733 * of extra communication calls when the cut-off
2734 * would be only slightly longer than the cell size.
2735 * Later cellsize_limit is redetermined,
2736 * so we can not miss interactions due to this rounding.
2738 npulse
= static_cast<int>(0.96 + comm
->cutoff
/comm
->cellsize_limit
);
2742 /* There is no cell size limit */
2743 npulse
= std::max(dd
->nc
[XX
]-1, std::max(dd
->nc
[YY
]-1, dd
->nc
[ZZ
]-1));
2746 if (!bNoCutOff
&& npulse
> 1)
2748 /* See if we can do with less pulses, based on dlb_scale */
2750 for (d
= 0; d
< dd
->ndim
; d
++)
2753 npulse_d
= static_cast<int>(1 + dd
->nc
[dim
]*comm
->cutoff
2754 /(ddbox
->box_size
[dim
]*ddbox
->skew_fac
[dim
]*dlb_scale
));
2755 npulse_d_max
= std::max(npulse_d_max
, npulse_d
);
2757 npulse
= std::min(npulse
, npulse_d_max
);
2760 /* This env var can override npulse */
2761 d
= dd_getenv(mdlog
, "GMX_DD_NPULSE", 0);
2768 comm
->bVacDLBNoLimit
= (ir
->ePBC
== epbcNONE
);
2769 for (d
= 0; d
< dd
->ndim
; d
++)
2771 comm
->cd
[d
].np_dlb
= std::min(npulse
, dd
->nc
[dd
->dim
[d
]]-1);
2772 comm
->maxpulse
= std::max(comm
->maxpulse
, comm
->cd
[d
].np_dlb
);
2773 if (comm
->cd
[d
].np_dlb
< dd
->nc
[dd
->dim
[d
]]-1)
2775 comm
->bVacDLBNoLimit
= FALSE
;
2779 /* cellsize_limit is set for LINCS in init_domain_decomposition */
2780 if (!comm
->bVacDLBNoLimit
)
2782 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
,
2783 comm
->cutoff
/comm
->maxpulse
);
2785 comm
->cellsize_limit
= std::max(comm
->cellsize_limit
, comm
->cutoff_mbody
);
2786 /* Set the minimum cell size for each DD dimension */
2787 for (d
= 0; d
< dd
->ndim
; d
++)
2789 if (comm
->bVacDLBNoLimit
||
2790 comm
->cd
[d
].np_dlb
*comm
->cellsize_limit
>= comm
->cutoff
)
2792 comm
->cellsize_min_dlb
[dd
->dim
[d
]] = comm
->cellsize_limit
;
2796 comm
->cellsize_min_dlb
[dd
->dim
[d
]] =
2797 comm
->cutoff
/comm
->cd
[d
].np_dlb
;
2800 if (comm
->cutoff_mbody
<= 0)
2802 comm
->cutoff_mbody
= std::min(comm
->cutoff
, comm
->cellsize_limit
);
2810 gmx_bool
dd_bonded_molpbc(const gmx_domdec_t
*dd
, int ePBC
)
2812 /* If each molecule is a single charge group
2813 * or we use domain decomposition for each periodic dimension,
2814 * we do not need to take pbc into account for the bonded interactions.
2816 return (ePBC
!= epbcNONE
&& dd
->comm
->haveInterDomainBondeds
&&
2819 (dd
->nc
[ZZ
] > 1 || ePBC
== epbcXY
)));
2822 /*! \brief Sets grid size limits and PP-PME setup, prints settings to log */
2823 static void set_ddgrid_parameters(const gmx::MDLogger
&mdlog
,
2824 gmx_domdec_t
*dd
, real dlb_scale
,
2825 const gmx_mtop_t
*mtop
, const t_inputrec
*ir
,
2826 const gmx_ddbox_t
*ddbox
)
2828 gmx_domdec_comm_t
*comm
;
2834 if (EEL_PME(ir
->coulombtype
) || EVDW_PME(ir
->vdwtype
))
2836 init_ddpme(dd
, &comm
->ddpme
[0], 0);
2837 if (comm
->npmedecompdim
>= 2)
2839 init_ddpme(dd
, &comm
->ddpme
[1], 1);
2844 comm
->npmenodes
= 0;
2845 if (dd
->pme_nodeid
>= 0)
2847 gmx_fatal_collective(FARGS
, dd
->mpi_comm_all
, DDMASTER(dd
),
2848 "Can not have separate PME ranks without PME electrostatics");
2854 fprintf(debug
, "The DD cut-off is %f\n", comm
->cutoff
);
2856 if (!isDlbDisabled(comm
))
2858 set_cell_limits_dlb(mdlog
, dd
, dlb_scale
, ir
, ddbox
);
2861 logSettings(mdlog
, dd
, mtop
, ir
, dlb_scale
, ddbox
);
2863 if (ir
->ePBC
== epbcNONE
)
2865 vol_frac
= 1 - 1/static_cast<double>(dd
->nnodes
);
2870 (1 + comm_box_frac(dd
->nc
, comm
->cutoff
, ddbox
))/static_cast<double>(dd
->nnodes
);
2874 fprintf(debug
, "Volume fraction for all DD zones: %f\n", vol_frac
);
2876 natoms_tot
= comm
->cgs_gl
.index
[comm
->cgs_gl
.nr
];
2878 dd
->ga2la
= new gmx_ga2la_t(natoms_tot
,
2879 static_cast<int>(vol_frac
*natoms_tot
));
2882 /*! \brief Set some important DD parameters that can be modified by env.vars */
2883 static void set_dd_envvar_options(const gmx::MDLogger
&mdlog
,
2884 gmx_domdec_t
*dd
, int rank_mysim
)
2886 gmx_domdec_comm_t
*comm
= dd
->comm
;
2888 dd
->bSendRecv2
= (dd_getenv(mdlog
, "GMX_DD_USE_SENDRECV2", 0) != 0);
2889 comm
->dlb_scale_lim
= dd_getenv(mdlog
, "GMX_DLB_MAX_BOX_SCALING", 10);
2890 comm
->eFlop
= dd_getenv(mdlog
, "GMX_DLB_BASED_ON_FLOPS", 0);
2891 int recload
= dd_getenv(mdlog
, "GMX_DD_RECORD_LOAD", 1);
2892 comm
->nstDDDump
= dd_getenv(mdlog
, "GMX_DD_NST_DUMP", 0);
2893 comm
->nstDDDumpGrid
= dd_getenv(mdlog
, "GMX_DD_NST_DUMP_GRID", 0);
2894 comm
->DD_debug
= dd_getenv(mdlog
, "GMX_DD_DEBUG", 0);
2898 GMX_LOG(mdlog
.info
).appendText("Will use two sequential MPI_Sendrecv calls instead of two simultaneous non-blocking MPI_Irecv and MPI_Isend pairs for constraint and vsite communication");
2903 GMX_LOG(mdlog
.info
).appendText("Will load balance based on FLOP count");
2904 if (comm
->eFlop
> 1)
2906 srand(1 + rank_mysim
);
2908 comm
->bRecordLoad
= TRUE
;
2912 comm
->bRecordLoad
= (wallcycle_have_counter() && recload
> 0);
2916 gmx_domdec_t::gmx_domdec_t(const t_inputrec
&ir
) :
2921 gmx_domdec_t
*init_domain_decomposition(const gmx::MDLogger
&mdlog
,
2923 const DomdecOptions
&options
,
2924 const gmx::MdrunOptions
&mdrunOptions
,
2925 const gmx_mtop_t
*mtop
,
2926 const t_inputrec
*ir
,
2928 gmx::ArrayRef
<const gmx::RVec
> xGlobal
,
2929 gmx::LocalAtomSetManager
*atomSets
)
2933 GMX_LOG(mdlog
.info
).appendTextFormatted(
2934 "\nInitializing Domain Decomposition on %d ranks", cr
->nnodes
);
2936 dd
= new gmx_domdec_t(*ir
);
2938 dd
->comm
= init_dd_comm();
2940 /* Initialize DD paritioning counters */
2941 dd
->comm
->partition_step
= INT_MIN
;
2944 set_dd_envvar_options(mdlog
, dd
, cr
->nodeid
);
2946 gmx_ddbox_t ddbox
= {0};
2947 set_dd_limits_and_grid(mdlog
, cr
, dd
, options
, mdrunOptions
,
2952 make_dd_communicators(mdlog
, cr
, dd
, options
.rankOrder
);
2954 if (thisRankHasDuty(cr
, DUTY_PP
))
2956 set_ddgrid_parameters(mdlog
, dd
, options
.dlbScaling
, mtop
, ir
, &ddbox
);
2958 setup_neighbor_relations(dd
);
2961 /* Set overallocation to avoid frequent reallocation of arrays */
2962 set_over_alloc_dd(TRUE
);
2964 clear_dd_cycle_counts(dd
);
2966 dd
->atomSets
= atomSets
;
2971 static gmx_bool
test_dd_cutoff(t_commrec
*cr
,
2972 const t_state
&state
,
2973 real cutoffRequested
)
2983 set_ddbox(*dd
, false, state
.box
, true, state
.x
, &ddbox
);
2987 for (d
= 0; d
< dd
->ndim
; d
++)
2991 inv_cell_size
= DD_CELL_MARGIN
*dd
->nc
[dim
]/ddbox
.box_size
[dim
];
2992 if (dd
->unitCellInfo
.ddBoxIsDynamic
)
2994 inv_cell_size
*= DD_PRES_SCALE_MARGIN
;
2997 np
= 1 + static_cast<int>(cutoffRequested
*inv_cell_size
*ddbox
.skew_fac
[dim
]);
2999 if (!isDlbDisabled(dd
->comm
) && (dim
< ddbox
.npbcdim
) && (dd
->comm
->cd
[d
].np_dlb
> 0))
3001 if (np
> dd
->comm
->cd
[d
].np_dlb
)
3006 /* If a current local cell size is smaller than the requested
3007 * cut-off, we could still fix it, but this gets very complicated.
3008 * Without fixing here, we might actually need more checks.
3010 real cellSizeAlongDim
= (dd
->comm
->cell_x1
[dim
] - dd
->comm
->cell_x0
[dim
])*ddbox
.skew_fac
[dim
];
3011 if (cellSizeAlongDim
*dd
->comm
->cd
[d
].np_dlb
< cutoffRequested
)
3018 if (!isDlbDisabled(dd
->comm
))
3020 /* If DLB is not active yet, we don't need to check the grid jumps.
3021 * Actually we shouldn't, because then the grid jump data is not set.
3023 if (isDlbOn(dd
->comm
) &&
3024 check_grid_jump(0, dd
, cutoffRequested
, &ddbox
, FALSE
))
3029 gmx_sumi(1, &LocallyLimited
, cr
);
3031 if (LocallyLimited
> 0)
3040 gmx_bool
change_dd_cutoff(t_commrec
*cr
,
3041 const t_state
&state
,
3042 real cutoffRequested
)
3044 gmx_bool bCutoffAllowed
;
3046 bCutoffAllowed
= test_dd_cutoff(cr
, state
, cutoffRequested
);
3050 cr
->dd
->comm
->cutoff
= cutoffRequested
;
3053 return bCutoffAllowed
;