1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2006,2007,2008,2009,2010 by the GROMACS development team.
5 * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
6 * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37
38 /*! \internal \file
39 *
40 * \brief This file implements functions for domdec to use
41 * while managing inter-atomic constraints.
42 *
43 * \author Berk Hess <hess@kth.se>
44 * \ingroup module_domdec
45 */
46
47 #include "gmxpre.h"
48
49 #include "domdec_constraints.h"
50
51 #include <cassert>
52
53 #include <algorithm>
54 #include <memory>
55
56 #include "gromacs/domdec/dlbtiming.h"
57 #include "gromacs/domdec/domdec.h"
58 #include "gromacs/domdec/domdec_struct.h"
59 #include "gromacs/domdec/ga2la.h"
60 #include "gromacs/domdec/hashedmap.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdlib/constr.h"
63 #include "gromacs/mdlib/gmx_omp_nthreads.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/forcerec.h" // only for GET_CGINFO_*
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/topology/ifunc.h"
68 #include "gromacs/topology/mtop_lookup.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/gmxassert.h"
72 #include "gromacs/utility/listoflists.h"
73
74 #include "domdec_internal.h"
75 #include "domdec_specatomcomm.h"
76
77 using gmx::ListOfLists;
78
79 /*! \brief Struct used during constraint setup with domain decomposition */
80 struct gmx_domdec_constraints_t
81 {
82 //! @cond Doxygen_Suppress
83 std::vector<int> molb_con_offset; /**< Offset in the constraint array for each molblock */
84 std::vector<int> molb_ncon_mol; /**< The number of constraints per molecule for each molblock */
85
86 int ncon; /**< The fully local and conneced constraints */
87 /* The global constraint number, only required for clearing gc_req */
88 std::vector<int> con_gl; /**< Global constraint indices for local constraints */
89 std::vector<int> con_nlocat; /**< Number of local atoms (2/1/0) for each constraint */
90
91 std::vector<bool> gc_req; /**< Boolean that tells if a global constraint index has been requested; note: size global #constraints */
92
93 /* Hash table for keeping track of requests */
94 std::unique_ptr<gmx::HashedMap<int>> ga2la; /**< Global to local communicated constraint atom only index */
95
96 /* Multi-threading stuff */
97 int nthread; /**< Number of threads used for DD constraint setup */
98 std::vector<InteractionList> ils; /**< Constraint ilist working arrays, size \p nthread */
99
100 /* Buffers for requesting atoms */
101 std::vector<std::vector<int>> requestedGlobalAtomIndices; /**< Buffers for requesting global atom indices, one per thread */
102
103 //! @endcond
104 };
105
dd_move_x_constraints(gmx_domdec_t * dd,const matrix box,gmx::ArrayRef<gmx::RVec> x0,gmx::ArrayRef<gmx::RVec> x1,gmx_bool bX1IsCoord)106 void dd_move_x_constraints(gmx_domdec_t* dd,
107 const matrix box,
108 gmx::ArrayRef<gmx::RVec> x0,
109 gmx::ArrayRef<gmx::RVec> x1,
110 gmx_bool bX1IsCoord)
111 {
112 if (dd->constraint_comm)
113 {
114 dd_move_x_specat(dd, dd->constraint_comm, box, as_rvec_array(x0.data()),
115 as_rvec_array(x1.data()), bX1IsCoord);
116
117 ddReopenBalanceRegionCpu(dd);
118 }
119 }
120
dd_constraints_nlocalatoms(const gmx_domdec_t * dd)121 gmx::ArrayRef<const int> dd_constraints_nlocalatoms(const gmx_domdec_t* dd)
122 {
123 if (dd && dd->constraints)
124 {
125 return dd->constraints->con_nlocat;
126 }
127 else
128 {
129 return {};
130 }
131 }
132
dd_clear_local_constraint_indices(gmx_domdec_t * dd)133 void dd_clear_local_constraint_indices(gmx_domdec_t* dd)
134 {
135 gmx_domdec_constraints_t* dc = dd->constraints;
136
137 std::fill(dc->gc_req.begin(), dc->gc_req.end(), false);
138
139 if (dd->constraint_comm)
140 {
141 dc->ga2la->clearAndResizeHashTable();
142 }
143 }
144
145 /*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */
walk_out(int con,int con_offset,int a,int offset,int nrec,gmx::ArrayRef<const int> ia1,gmx::ArrayRef<const int> ia2,const ListOfLists<int> & at2con,const gmx_ga2la_t & ga2la,gmx_bool bHomeConnect,gmx_domdec_constraints_t * dc,gmx_domdec_specat_comm_t * dcc,InteractionList * il_local,std::vector<int> * ireq)146 static void walk_out(int con,
147 int con_offset,
148 int a,
149 int offset,
150 int nrec,
151 gmx::ArrayRef<const int> ia1,
152 gmx::ArrayRef<const int> ia2,
153 const ListOfLists<int>& at2con,
154 const gmx_ga2la_t& ga2la,
155 gmx_bool bHomeConnect,
156 gmx_domdec_constraints_t* dc,
157 gmx_domdec_specat_comm_t* dcc,
158 InteractionList* il_local,
159 std::vector<int>* ireq)
160 {
161 if (!dc->gc_req[con_offset + con])
162 {
163 /* Add this non-home constraint to the list */
164 dc->con_gl.push_back(con_offset + con);
165 dc->con_nlocat.push_back(bHomeConnect ? 1 : 0);
166 dc->gc_req[con_offset + con] = true;
167 const int* iap = constr_iatomptr(ia1, ia2, con);
168 const int parameterType = iap[0];
169 const int a1_gl = offset + iap[1];
170 const int a2_gl = offset + iap[2];
171 std::array<int, 2> atoms;
172 /* The following indexing code can probably be optizimed */
173 if (const int* a_loc = ga2la.findHome(a1_gl))
174 {
175 atoms[0] = *a_loc;
176 }
177 else
178 {
179 /* We set this index later */
180 atoms[0] = -a1_gl - 1;
181 }
182 if (const int* a_loc = ga2la.findHome(a2_gl))
183 {
184 atoms[1] = *a_loc;
185 }
186 else
187 {
188 /* We set this index later */
189 atoms[1] = -a2_gl - 1;
190 }
191 il_local->push_back(parameterType, atoms);
192 dc->ncon++;
193 }
194 /* Check to not ask for the same atom more than once */
195 if (!dc->ga2la->find(offset + a))
196 {
197 assert(dcc);
198 /* Add this non-home atom to the list */
199 ireq->push_back(offset + a);
200 /* Temporarily mark with -2, we get the index later */
201 dc->ga2la->insert(offset + a, -2);
202 }
203
204 if (nrec > 0)
205 {
206 /* Loop over the constraint connected to atom a */
207 for (const int coni : at2con[a])
208 {
209 if (coni != con)
210 {
211 /* Walk further */
212 const int* iap = constr_iatomptr(ia1, ia2, coni);
213 int b;
214 if (a == iap[1])
215 {
216 b = iap[2];
217 }
218 else
219 {
220 b = iap[1];
221 }
222 if (!ga2la.findHome(offset + b))
223 {
224 walk_out(coni, con_offset, b, offset, nrec - 1, ia1, ia2, at2con, ga2la, FALSE,
225 dc, dcc, il_local, ireq);
226 }
227 }
228 }
229 }
230 }
231
232 /*! \brief Looks up SETTLE constraints for a range of charge-groups */
atoms_to_settles(gmx_domdec_t * dd,const gmx_mtop_t * mtop,const int * cginfo,gmx::ArrayRef<const std::vector<int>> at2settle_mt,int cg_start,int cg_end,InteractionList * ils_local,std::vector<int> * ireq)233 static void atoms_to_settles(gmx_domdec_t* dd,
234 const gmx_mtop_t* mtop,
235 const int* cginfo,
236 gmx::ArrayRef<const std::vector<int>> at2settle_mt,
237 int cg_start,
238 int cg_end,
239 InteractionList* ils_local,
240 std::vector<int>* ireq)
241 {
242 const gmx_ga2la_t& ga2la = *dd->ga2la;
243 int nral = NRAL(F_SETTLE);
244
245 int mb = 0;
246 for (int a = cg_start; a < cg_end; a++)
247 {
248 if (GET_CGINFO_SETTLE(cginfo[a]))
249 {
250 int a_gl = dd->globalAtomIndices[a];
251 int a_mol;
252 mtopGetMolblockIndex(mtop, a_gl, &mb, nullptr, &a_mol);
253
254 const gmx_molblock_t* molb = &mtop->molblock[mb];
255 int settle = at2settle_mt[molb->type][a_mol];
256
257 if (settle >= 0)
258 {
259 int offset = a_gl - a_mol;
260
261 const int* ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms.data();
262
263 int a_gls[3];
264 gmx_bool bAssign = FALSE;
265 int nlocal = 0;
266 for (int sa = 0; sa < nral; sa++)
267 {
268 int a_glsa = offset + ia1[settle * (1 + nral) + 1 + sa];
269 a_gls[sa] = a_glsa;
270 if (ga2la.findHome(a_glsa))
271 {
272 if (nlocal == 0 && a_gl == a_glsa)
273 {
274 bAssign = TRUE;
275 }
276 nlocal++;
277 }
278 }
279
280 if (bAssign)
281 {
282 const int parameterType = ia1[settle * 4];
283 std::array<int, 3> atoms;
284 for (int sa = 0; sa < nral; sa++)
285 {
286 if (const int* a_loc = ga2la.findHome(a_gls[sa]))
287 {
288 atoms[sa] = *a_loc;
289 }
290 else
291 {
292 atoms[sa] = -a_gls[sa] - 1;
293 /* Add this non-home atom to the list */
294 ireq->push_back(a_gls[sa]);
295 /* A check on double atom requests is
296 * not required for settle.
297 */
298 }
299 }
300 ils_local->push_back(parameterType, atoms);
301 }
302 }
303 }
304 }
305 }
306
307 /*! \brief Looks up constraint for the local atoms */
atoms_to_constraints(gmx_domdec_t * dd,const gmx_mtop_t * mtop,const int * cginfo,gmx::ArrayRef<const ListOfLists<int>> at2con_mt,int nrec,InteractionList * ilc_local,std::vector<int> * ireq)308 static void atoms_to_constraints(gmx_domdec_t* dd,
309 const gmx_mtop_t* mtop,
310 const int* cginfo,
311 gmx::ArrayRef<const ListOfLists<int>> at2con_mt,
312 int nrec,
313 InteractionList* ilc_local,
314 std::vector<int>* ireq)
315 {
316 gmx_domdec_constraints_t* dc = dd->constraints;
317 gmx_domdec_specat_comm_t* dcc = dd->constraint_comm;
318
319 const gmx_ga2la_t& ga2la = *dd->ga2la;
320
321 dc->con_gl.clear();
322 dc->con_nlocat.clear();
323
324 int mb = 0;
325 int nhome = 0;
326 for (int a = 0; a < dd->ncg_home; a++)
327 {
328 if (GET_CGINFO_CONSTR(cginfo[a]))
329 {
330 int a_gl = dd->globalAtomIndices[a];
331 int molnr, a_mol;
332 mtopGetMolblockIndex(mtop, a_gl, &mb, &molnr, &a_mol);
333
334 const gmx_molblock_t& molb = mtop->molblock[mb];
335
336 gmx::ArrayRef<const int> ia1 = mtop->moltype[molb.type].ilist[F_CONSTR].iatoms;
337 gmx::ArrayRef<const int> ia2 = mtop->moltype[molb.type].ilist[F_CONSTRNC].iatoms;
338
339 /* Calculate the global constraint number offset for the molecule.
340 * This is only required for the global index to make sure
341 * that we use each constraint only once.
342 */
343 const int con_offset = dc->molb_con_offset[mb] + molnr * dc->molb_ncon_mol[mb];
344
345 /* The global atom number offset for this molecule */
346 const int offset = a_gl - a_mol;
347 /* Loop over the constraints connected to atom a_mol in the molecule */
348 const auto& at2con = at2con_mt[molb.type];
349 for (const int con : at2con[a_mol])
350 {
351 const int* iap = constr_iatomptr(ia1, ia2, con);
352 int b_mol;
353 if (a_mol == iap[1])
354 {
355 b_mol = iap[2];
356 }
357 else
358 {
359 b_mol = iap[1];
360 }
361 if (const int* a_loc = ga2la.findHome(offset + b_mol))
362 {
363 /* Add this fully home constraint at the first atom */
364 if (a_mol < b_mol)
365 {
366 dc->con_gl.push_back(con_offset + con);
367 dc->con_nlocat.push_back(2);
368 const int b_lo = *a_loc;
369 const int parameterType = iap[0];
370 std::array<int, 2> atoms;
371 atoms[0] = (a_gl == iap[1] ? a : b_lo);
372 atoms[1] = (a_gl == iap[1] ? b_lo : a);
373 ilc_local->push_back(parameterType, atoms);
374 dc->ncon++;
375 nhome++;
376 }
377 }
378 else
379 {
380 /* We need the nrec constraints coupled to this constraint,
381 * so we need to walk out of the home cell by nrec+1 atoms,
382 * since already atom bg is not locally present.
383 * Therefore we call walk_out with nrec recursions to go
384 * after this first call.
385 */
386 walk_out(con, con_offset, b_mol, offset, nrec, ia1, ia2, at2con, ga2la, TRUE,
387 dc, dcc, ilc_local, ireq);
388 }
389 }
390 }
391 }
392
393 GMX_ASSERT(dc->con_gl.size() == static_cast<size_t>(dc->ncon),
394 "con_gl size should match the number of constraints");
395 GMX_ASSERT(dc->con_nlocat.size() == static_cast<size_t>(dc->ncon),
396 "con_nlocat size should match the number of constraints");
397
398 if (debug)
399 {
400 fprintf(debug, "Constraints: home %3d border %3d atoms: %3zu\n", nhome, dc->ncon - nhome,
401 dd->constraint_comm ? ireq->size() : 0);
402 }
403 }
404
dd_make_local_constraints(gmx_domdec_t * dd,int at_start,const struct gmx_mtop_t * mtop,const int * cginfo,gmx::Constraints * constr,int nrec,gmx::ArrayRef<InteractionList> il_local)405 int dd_make_local_constraints(gmx_domdec_t* dd,
406 int at_start,
407 const struct gmx_mtop_t* mtop,
408 const int* cginfo,
409 gmx::Constraints* constr,
410 int nrec,
411 gmx::ArrayRef<InteractionList> il_local)
412 {
413 gmx_domdec_constraints_t* dc;
414 InteractionList * ilc_local, *ils_local;
415 gmx::HashedMap<int>* ga2la_specat;
416 int at_end, i, j;
417
418 // This code should not be called unless this condition is true,
419 // because that's the only time init_domdec_constraints is
420 // called...
421 GMX_RELEASE_ASSERT(dd->comm->systemInfo.haveSplitConstraints || dd->comm->systemInfo.haveSplitSettles,
422 "dd_make_local_constraints called when there are no local constraints");
423 // ... and init_domdec_constraints always sets
424 // dd->constraint_comm...
425 GMX_RELEASE_ASSERT(
426 dd->constraint_comm,
427 "Invalid use of dd_make_local_constraints before construction of constraint_comm");
428 // ... which static analysis needs to be reassured about, because
429 // otherwise, when dd->splitSettles is
430 // true. dd->constraint_comm is unilaterally dereferenced before
431 // the call to atoms_to_settles.
432
433 dc = dd->constraints;
434
435 ilc_local = &il_local[F_CONSTR];
436 ils_local = &il_local[F_SETTLE];
437
438 dc->ncon = 0;
439 gmx::ArrayRef<const ListOfLists<int>> at2con_mt;
440 std::vector<int>* ireq = nullptr;
441 ilc_local->clear();
442 if (dd->constraint_comm)
443 {
444 // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
445 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
446 at2con_mt = constr->atom2constraints_moltype();
447 ireq = &dc->requestedGlobalAtomIndices[0];
448 ireq->clear();
449 }
450
451 gmx::ArrayRef<const std::vector<int>> at2settle_mt;
452 /* When settle works inside charge groups, we assigned them already */
453 if (dd->comm->systemInfo.haveSplitSettles)
454 {
455 // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
456 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
457 at2settle_mt = constr->atom2settle_moltype();
458 ils_local->clear();
459 }
460
461 if (at2settle_mt.empty())
462 {
463 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec, ilc_local, ireq);
464 }
465 else
466 {
467 int t0_set;
468
469 /* Do the constraints, if present, on the first thread.
470 * Do the settles on all other threads.
471 */
472 t0_set = ((!at2con_mt.empty() && dc->nthread > 1) ? 1 : 0);
473
474 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
475 for (int thread = 0; thread < dc->nthread; thread++)
476 {
477 try
478 {
479 if (!at2con_mt.empty() && thread == 0)
480 {
481 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec, ilc_local, ireq);
482 }
483
484 if (thread >= t0_set)
485 {
486 int cg0, cg1;
487 InteractionList* ilst;
488
489 /* Distribute the settle check+assignments over
490 * dc->nthread or dc->nthread-1 threads.
491 */
492 cg0 = (dd->ncg_home * (thread - t0_set)) / (dc->nthread - t0_set);
493 cg1 = (dd->ncg_home * (thread - t0_set + 1)) / (dc->nthread - t0_set);
494
495 if (thread == t0_set)
496 {
497 ilst = ils_local;
498 }
499 else
500 {
501 ilst = &dc->ils[thread];
502 }
503 ilst->clear();
504
505 std::vector<int>& ireqt = dc->requestedGlobalAtomIndices[thread];
506 if (thread > 0)
507 {
508 ireqt.clear();
509 }
510
511 atoms_to_settles(dd, mtop, cginfo, at2settle_mt, cg0, cg1, ilst, &ireqt);
512 }
513 }
514 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
515 }
516
517 /* Combine the generate settles and requested indices */
518 for (int thread = 1; thread < dc->nthread; thread++)
519 {
520 if (thread > t0_set)
521 {
522 ils_local->append(dc->ils[thread]);
523 }
524
525 const std::vector<int>& ireqt = dc->requestedGlobalAtomIndices[thread];
526 ireq->insert(ireq->end(), ireqt.begin(), ireqt.end());
527 }
528
529 if (debug)
530 {
531 fprintf(debug, "Settles: total %3d\n", ils_local->size() / 4);
532 }
533 }
534
535 if (dd->constraint_comm)
536 {
537 int nral1;
538
539 at_end = setup_specat_communication(dd, ireq, dd->constraint_comm, dd->constraints->ga2la.get(),
540 at_start, 2, "constraint", " or lincs-order");
541
542 /* Fill in the missing indices */
543 ga2la_specat = dd->constraints->ga2la.get();
544
545 nral1 = 1 + NRAL(F_CONSTR);
546 for (i = 0; i < ilc_local->size(); i += nral1)
547 {
548 int* iap = ilc_local->iatoms.data() + i;
549 for (j = 1; j < nral1; j++)
550 {
551 if (iap[j] < 0)
552 {
553 const int* a = ga2la_specat->find(-iap[j] - 1);
554 GMX_ASSERT(a, "We have checked before that this atom index has been set");
555 iap[j] = *a;
556 }
557 }
558 }
559
560 nral1 = 1 + NRAL(F_SETTLE);
561 for (i = 0; i < ils_local->size(); i += nral1)
562 {
563 int* iap = ils_local->iatoms.data() + i;
564 for (j = 1; j < nral1; j++)
565 {
566 if (iap[j] < 0)
567 {
568 const int* a = ga2la_specat->find(-iap[j] - 1);
569 GMX_ASSERT(a, "We have checked before that this atom index has been set");
570 iap[j] = *a;
571 }
572 }
573 }
574 }
575 else
576 {
577 // Currently unreachable
578 at_end = at_start;
579 }
580
581 return at_end;
582 }
583
init_domdec_constraints(gmx_domdec_t * dd,const gmx_mtop_t * mtop)584 void init_domdec_constraints(gmx_domdec_t* dd, const gmx_mtop_t* mtop)
585 {
586 gmx_domdec_constraints_t* dc;
587 const gmx_molblock_t* molb;
588
589 if (debug)
590 {
591 fprintf(debug, "Begin init_domdec_constraints\n");
592 }
593
594 dd->constraints = new gmx_domdec_constraints_t;
595 dc = dd->constraints;
596
597 dc->molb_con_offset.resize(mtop->molblock.size());
598 dc->molb_ncon_mol.resize(mtop->molblock.size());
599
600 int ncon = 0;
601 for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
602 {
603 molb = &mtop->molblock[mb];
604 dc->molb_con_offset[mb] = ncon;
605 dc->molb_ncon_mol[mb] = mtop->moltype[molb->type].ilist[F_CONSTR].size() / 3
606 + mtop->moltype[molb->type].ilist[F_CONSTRNC].size() / 3;
607 ncon += molb->nmol * dc->molb_ncon_mol[mb];
608 }
609
610 if (ncon > 0)
611 {
612 dc->gc_req.resize(ncon);
613 }
614
615 /* Use a hash table for the global to local index.
616 * The number of keys is a rough estimate, it will be optimized later.
617 */
618 int numKeysEstimate = std::min(mtop->natoms / 20, mtop->natoms / (2 * dd->nnodes));
619 dc->ga2la = std::make_unique<gmx::HashedMap<int>>(numKeysEstimate);
620
621 dc->nthread = gmx_omp_nthreads_get(emntDomdec);
622 dc->ils.resize(dc->nthread);
623
624 dd->constraint_comm = new gmx_domdec_specat_comm_t;
625
626 dc->requestedGlobalAtomIndices.resize(dc->nthread);
627 }
628