2 * This file is part of the GROMACS molecular simulation package.
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, 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.
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.
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.
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.
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.
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.
40 * \brief This file implements functions for domdec to use
41 * while managing inter-atomic constraints.
43 * \author Berk Hess <hess@kth.se>
44 * \ingroup module_domdec
49 #include "domdec_constraints.h"
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"
74 #include "domdec_internal.h"
75 #include "domdec_specatomcomm.h"
77 using gmx::ListOfLists
;
79 /*! \brief Struct used during constraint setup with domain decomposition */
80 struct gmx_domdec_constraints_t
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 */
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 */
91 std::vector
<bool> gc_req
; /**< Boolean that tells if a global constraint index has been requested; note: size global #constraints */
93 /* Hash table for keeping track of requests */
94 std::unique_ptr
<gmx::HashedMap
<int>> ga2la
; /**< Global to local communicated constraint atom only index */
96 /* Multi-threading stuff */
97 int nthread
; /**< Number of threads used for DD constraint setup */
98 std::vector
<t_ilist
> ils
; /**< Constraint ilist working arrays, size \p nthread */
100 /* Buffers for requesting atoms */
101 std::vector
<std::vector
<int>> requestedGlobalAtomIndices
; /**< Buffers for requesting global atom indices, one per thread */
106 void dd_move_x_constraints(gmx_domdec_t
* dd
, const matrix box
, rvec
* x0
, rvec
* x1
, gmx_bool bX1IsCoord
)
108 if (dd
->constraint_comm
)
110 dd_move_x_specat(dd
, dd
->constraint_comm
, box
, x0
, x1
, bX1IsCoord
);
112 ddReopenBalanceRegionCpu(dd
);
116 gmx::ArrayRef
<const int> dd_constraints_nlocalatoms(const gmx_domdec_t
* dd
)
118 if (dd
&& dd
->constraints
)
120 return dd
->constraints
->con_nlocat
;
128 void dd_clear_local_constraint_indices(gmx_domdec_t
* dd
)
130 gmx_domdec_constraints_t
* dc
= dd
->constraints
;
132 std::fill(dc
->gc_req
.begin(), dc
->gc_req
.end(), false);
134 if (dd
->constraint_comm
)
140 /*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */
141 static void walk_out(int con
,
146 gmx::ArrayRef
<const int> ia1
,
147 gmx::ArrayRef
<const int> ia2
,
148 const ListOfLists
<int>& at2con
,
149 const gmx_ga2la_t
& ga2la
,
150 gmx_bool bHomeConnect
,
151 gmx_domdec_constraints_t
* dc
,
152 gmx_domdec_specat_comm_t
* dcc
,
154 std::vector
<int>* ireq
)
156 if (!dc
->gc_req
[con_offset
+ con
])
158 /* Add this non-home constraint to the list */
159 dc
->con_gl
.push_back(con_offset
+ con
);
160 dc
->con_nlocat
.push_back(bHomeConnect
? 1 : 0);
161 dc
->gc_req
[con_offset
+ con
] = true;
162 if (il_local
->nr
+ 3 > il_local
->nalloc
)
164 il_local
->nalloc
= over_alloc_dd(il_local
->nr
+ 3);
165 srenew(il_local
->iatoms
, il_local
->nalloc
);
167 const int* iap
= constr_iatomptr(ia1
, ia2
, con
);
168 il_local
->iatoms
[il_local
->nr
++] = iap
[0];
169 const int a1_gl
= offset
+ iap
[1];
170 const int a2_gl
= offset
+ iap
[2];
171 /* The following indexing code can probably be optizimed */
172 if (const int* a_loc
= ga2la
.findHome(a1_gl
))
174 il_local
->iatoms
[il_local
->nr
++] = *a_loc
;
178 /* We set this index later */
179 il_local
->iatoms
[il_local
->nr
++] = -a1_gl
- 1;
181 if (const int* a_loc
= ga2la
.findHome(a2_gl
))
183 il_local
->iatoms
[il_local
->nr
++] = *a_loc
;
187 /* We set this index later */
188 il_local
->iatoms
[il_local
->nr
++] = -a2_gl
- 1;
192 /* Check to not ask for the same atom more than once */
193 if (!dc
->ga2la
->find(offset
+ a
))
196 /* Add this non-home atom to the list */
197 ireq
->push_back(offset
+ a
);
198 /* Temporarily mark with -2, we get the index later */
199 dc
->ga2la
->insert(offset
+ a
, -2);
204 /* Loop over the constraint connected to atom a */
205 for (const int coni
: at2con
[a
])
210 const int* iap
= constr_iatomptr(ia1
, ia2
, coni
);
220 if (!ga2la
.findHome(offset
+ b
))
222 walk_out(coni
, con_offset
, b
, offset
, nrec
- 1, ia1
, ia2
, at2con
, ga2la
, FALSE
,
223 dc
, dcc
, il_local
, ireq
);
230 /*! \brief Looks up SETTLE constraints for a range of charge-groups */
231 static void atoms_to_settles(gmx_domdec_t
* dd
,
232 const gmx_mtop_t
* mtop
,
234 gmx::ArrayRef
<const std::vector
<int>> at2settle_mt
,
238 std::vector
<int>* ireq
)
240 const gmx_ga2la_t
& ga2la
= *dd
->ga2la
;
241 int nral
= NRAL(F_SETTLE
);
244 for (int a
= cg_start
; a
< cg_end
; a
++)
246 if (GET_CGINFO_SETTLE(cginfo
[a
]))
248 int a_gl
= dd
->globalAtomIndices
[a
];
250 mtopGetMolblockIndex(mtop
, a_gl
, &mb
, nullptr, &a_mol
);
252 const gmx_molblock_t
* molb
= &mtop
->molblock
[mb
];
253 int settle
= at2settle_mt
[molb
->type
][a_mol
];
257 int offset
= a_gl
- a_mol
;
259 const int* ia1
= mtop
->moltype
[molb
->type
].ilist
[F_SETTLE
].iatoms
.data();
262 gmx_bool bAssign
= FALSE
;
264 for (int sa
= 0; sa
< nral
; sa
++)
266 int a_glsa
= offset
+ ia1
[settle
* (1 + nral
) + 1 + sa
];
268 if (ga2la
.findHome(a_glsa
))
270 if (nlocal
== 0 && a_gl
== a_glsa
)
280 if (ils_local
->nr
+ 1 + nral
> ils_local
->nalloc
)
282 ils_local
->nalloc
= over_alloc_dd(ils_local
->nr
+ 1 + nral
);
283 srenew(ils_local
->iatoms
, ils_local
->nalloc
);
286 ils_local
->iatoms
[ils_local
->nr
++] = ia1
[settle
* 4];
288 for (int sa
= 0; sa
< nral
; sa
++)
290 if (const int* a_loc
= ga2la
.findHome(a_gls
[sa
]))
292 ils_local
->iatoms
[ils_local
->nr
++] = *a_loc
;
296 ils_local
->iatoms
[ils_local
->nr
++] = -a_gls
[sa
] - 1;
297 /* Add this non-home atom to the list */
298 ireq
->push_back(a_gls
[sa
]);
299 /* A check on double atom requests is
300 * not required for settle.
310 /*! \brief Looks up constraint for the local atoms */
311 static void atoms_to_constraints(gmx_domdec_t
* dd
,
312 const gmx_mtop_t
* mtop
,
314 gmx::ArrayRef
<const ListOfLists
<int>> at2con_mt
,
317 std::vector
<int>* ireq
)
319 gmx_domdec_constraints_t
* dc
= dd
->constraints
;
320 gmx_domdec_specat_comm_t
* dcc
= dd
->constraint_comm
;
322 const gmx_ga2la_t
& ga2la
= *dd
->ga2la
;
325 dc
->con_nlocat
.clear();
329 for (int a
= 0; a
< dd
->ncg_home
; a
++)
331 if (GET_CGINFO_CONSTR(cginfo
[a
]))
333 int a_gl
= dd
->globalAtomIndices
[a
];
335 mtopGetMolblockIndex(mtop
, a_gl
, &mb
, &molnr
, &a_mol
);
337 const gmx_molblock_t
& molb
= mtop
->molblock
[mb
];
339 gmx::ArrayRef
<const int> ia1
= mtop
->moltype
[molb
.type
].ilist
[F_CONSTR
].iatoms
;
340 gmx::ArrayRef
<const int> ia2
= mtop
->moltype
[molb
.type
].ilist
[F_CONSTRNC
].iatoms
;
342 /* Calculate the global constraint number offset for the molecule.
343 * This is only required for the global index to make sure
344 * that we use each constraint only once.
346 const int con_offset
= dc
->molb_con_offset
[mb
] + molnr
* dc
->molb_ncon_mol
[mb
];
348 /* The global atom number offset for this molecule */
349 const int offset
= a_gl
- a_mol
;
350 /* Loop over the constraints connected to atom a_mol in the molecule */
351 const auto& at2con
= at2con_mt
[molb
.type
];
352 for (const int con
: at2con
[a_mol
])
354 const int* iap
= constr_iatomptr(ia1
, ia2
, con
);
364 if (const int* a_loc
= ga2la
.findHome(offset
+ b_mol
))
366 /* Add this fully home constraint at the first atom */
369 dc
->con_gl
.push_back(con_offset
+ con
);
370 dc
->con_nlocat
.push_back(2);
371 if (ilc_local
->nr
+ 3 > ilc_local
->nalloc
)
373 ilc_local
->nalloc
= over_alloc_dd(ilc_local
->nr
+ 3);
374 srenew(ilc_local
->iatoms
, ilc_local
->nalloc
);
376 const int b_lo
= *a_loc
;
377 ilc_local
->iatoms
[ilc_local
->nr
++] = iap
[0];
378 ilc_local
->iatoms
[ilc_local
->nr
++] = (a_gl
== iap
[1] ? a
: b_lo
);
379 ilc_local
->iatoms
[ilc_local
->nr
++] = (a_gl
== iap
[1] ? b_lo
: a
);
386 /* We need the nrec constraints coupled to this constraint,
387 * so we need to walk out of the home cell by nrec+1 atoms,
388 * since already atom bg is not locally present.
389 * Therefore we call walk_out with nrec recursions to go
390 * after this first call.
392 walk_out(con
, con_offset
, b_mol
, offset
, nrec
, ia1
, ia2
, at2con
, ga2la
, TRUE
,
393 dc
, dcc
, ilc_local
, ireq
);
399 GMX_ASSERT(dc
->con_gl
.size() == static_cast<size_t>(dc
->ncon
),
400 "con_gl size should match the number of constraints");
401 GMX_ASSERT(dc
->con_nlocat
.size() == static_cast<size_t>(dc
->ncon
),
402 "con_nlocat size should match the number of constraints");
406 fprintf(debug
, "Constraints: home %3d border %3d atoms: %3zu\n", nhome
, dc
->ncon
- nhome
,
407 dd
->constraint_comm
? ireq
->size() : 0);
411 int dd_make_local_constraints(gmx_domdec_t
* dd
,
413 const struct gmx_mtop_t
* mtop
,
415 gmx::Constraints
* constr
,
419 gmx_domdec_constraints_t
* dc
;
420 t_ilist
* ilc_local
, *ils_local
;
421 gmx::HashedMap
<int>* ga2la_specat
;
425 // This code should not be called unless this condition is true,
426 // because that's the only time init_domdec_constraints is
428 GMX_RELEASE_ASSERT(dd
->comm
->systemInfo
.haveSplitConstraints
|| dd
->comm
->systemInfo
.haveSplitSettles
,
429 "dd_make_local_constraints called when there are no local constraints");
430 // ... and init_domdec_constraints always sets
431 // dd->constraint_comm...
434 "Invalid use of dd_make_local_constraints before construction of constraint_comm");
435 // ... which static analysis needs to be reassured about, because
436 // otherwise, when dd->splitSettles is
437 // true. dd->constraint_comm is unilaterally dereferenced before
438 // the call to atoms_to_settles.
440 dc
= dd
->constraints
;
442 ilc_local
= &il_local
[F_CONSTR
];
443 ils_local
= &il_local
[F_SETTLE
];
447 gmx::ArrayRef
<const ListOfLists
<int>> at2con_mt
;
448 std::vector
<int>* ireq
= nullptr;
449 if (dd
->constraint_comm
)
451 // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
452 GMX_RELEASE_ASSERT(constr
!= nullptr, "Must have valid constraints object");
453 at2con_mt
= constr
->atom2constraints_moltype();
454 ireq
= &dc
->requestedGlobalAtomIndices
[0];
458 gmx::ArrayRef
<const std::vector
<int>> at2settle_mt
;
459 /* When settle works inside charge groups, we assigned them already */
460 if (dd
->comm
->systemInfo
.haveSplitSettles
)
462 // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
463 GMX_RELEASE_ASSERT(constr
!= nullptr, "Must have valid constraints object");
464 at2settle_mt
= constr
->atom2settle_moltype();
468 if (at2settle_mt
.empty())
470 atoms_to_constraints(dd
, mtop
, cginfo
, at2con_mt
, nrec
, ilc_local
, ireq
);
476 /* Do the constraints, if present, on the first thread.
477 * Do the settles on all other threads.
479 t0_set
= ((!at2con_mt
.empty() && dc
->nthread
> 1) ? 1 : 0);
481 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
482 for (int thread
= 0; thread
< dc
->nthread
; thread
++)
486 if (!at2con_mt
.empty() && thread
== 0)
488 atoms_to_constraints(dd
, mtop
, cginfo
, at2con_mt
, nrec
, ilc_local
, ireq
);
491 if (thread
>= t0_set
)
496 /* Distribute the settle check+assignments over
497 * dc->nthread or dc->nthread-1 threads.
499 cg0
= (dd
->ncg_home
* (thread
- t0_set
)) / (dc
->nthread
- t0_set
);
500 cg1
= (dd
->ncg_home
* (thread
- t0_set
+ 1)) / (dc
->nthread
- t0_set
);
502 if (thread
== t0_set
)
508 ilst
= &dc
->ils
[thread
];
512 std::vector
<int>& ireqt
= dc
->requestedGlobalAtomIndices
[thread
];
518 atoms_to_settles(dd
, mtop
, cginfo
, at2settle_mt
, cg0
, cg1
, ilst
, &ireqt
);
521 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
524 /* Combine the generate settles and requested indices */
525 for (int thread
= 1; thread
< dc
->nthread
; thread
++)
532 ilst
= &dc
->ils
[thread
];
533 if (ils_local
->nr
+ ilst
->nr
> ils_local
->nalloc
)
535 ils_local
->nalloc
= over_alloc_large(ils_local
->nr
+ ilst
->nr
);
536 srenew(ils_local
->iatoms
, ils_local
->nalloc
);
538 for (ia
= 0; ia
< ilst
->nr
; ia
++)
540 ils_local
->iatoms
[ils_local
->nr
+ ia
] = ilst
->iatoms
[ia
];
542 ils_local
->nr
+= ilst
->nr
;
545 const std::vector
<int>& ireqt
= dc
->requestedGlobalAtomIndices
[thread
];
546 ireq
->insert(ireq
->end(), ireqt
.begin(), ireqt
.end());
551 fprintf(debug
, "Settles: total %3d\n", ils_local
->nr
/ 4);
555 if (dd
->constraint_comm
)
559 at_end
= setup_specat_communication(dd
, ireq
, dd
->constraint_comm
, dd
->constraints
->ga2la
.get(),
560 at_start
, 2, "constraint", " or lincs-order");
562 /* Fill in the missing indices */
563 ga2la_specat
= dd
->constraints
->ga2la
.get();
565 nral1
= 1 + NRAL(F_CONSTR
);
566 for (i
= 0; i
< ilc_local
->nr
; i
+= nral1
)
568 iap
= ilc_local
->iatoms
+ i
;
569 for (j
= 1; j
< nral1
; j
++)
573 const int* a
= ga2la_specat
->find(-iap
[j
] - 1);
574 GMX_ASSERT(a
, "We have checked before that this atom index has been set");
580 nral1
= 1 + NRAL(F_SETTLE
);
581 for (i
= 0; i
< ils_local
->nr
; i
+= nral1
)
583 iap
= ils_local
->iatoms
+ i
;
584 for (j
= 1; j
< nral1
; j
++)
588 const int* a
= ga2la_specat
->find(-iap
[j
] - 1);
589 GMX_ASSERT(a
, "We have checked before that this atom index has been set");
597 // Currently unreachable
604 void init_domdec_constraints(gmx_domdec_t
* dd
, const gmx_mtop_t
* mtop
)
606 gmx_domdec_constraints_t
* dc
;
607 const gmx_molblock_t
* molb
;
611 fprintf(debug
, "Begin init_domdec_constraints\n");
614 dd
->constraints
= new gmx_domdec_constraints_t
;
615 dc
= dd
->constraints
;
617 dc
->molb_con_offset
.resize(mtop
->molblock
.size());
618 dc
->molb_ncon_mol
.resize(mtop
->molblock
.size());
621 for (size_t mb
= 0; mb
< mtop
->molblock
.size(); mb
++)
623 molb
= &mtop
->molblock
[mb
];
624 dc
->molb_con_offset
[mb
] = ncon
;
625 dc
->molb_ncon_mol
[mb
] = mtop
->moltype
[molb
->type
].ilist
[F_CONSTR
].size() / 3
626 + mtop
->moltype
[molb
->type
].ilist
[F_CONSTRNC
].size() / 3;
627 ncon
+= molb
->nmol
* dc
->molb_ncon_mol
[mb
];
632 dc
->gc_req
.resize(ncon
);
635 /* Use a hash table for the global to local index.
636 * The number of keys is a rough estimate, it will be optimized later.
638 int numKeysEstimate
= std::min(mtop
->natoms
/ 20, mtop
->natoms
/ (2 * dd
->nnodes
));
639 dc
->ga2la
= std::make_unique
<gmx::HashedMap
<int>>(numKeysEstimate
);
641 dc
->nthread
= gmx_omp_nthreads_get(emntDomdec
);
642 dc
->ils
.resize(dc
->nthread
);
644 dd
->constraint_comm
= new gmx_domdec_specat_comm_t
;
646 dc
->requestedGlobalAtomIndices
.resize(dc
->nthread
);