2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2006,2007,2008,2009,2010,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.
38 * \brief This file implements functions for domdec to use
39 * while managing inter-atomic constraints.
41 * \author Berk Hess <hess@kth.se>
42 * \ingroup module_domdec
47 #include "domdec_constraints.h"
54 #include "gromacs/domdec/dlbtiming.h"
55 #include "gromacs/domdec/domdec.h"
56 #include "gromacs/domdec/domdec_struct.h"
57 #include "gromacs/domdec/ga2la.h"
58 #include "gromacs/domdec/hashedmap.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/constr.h"
61 #include "gromacs/mdlib/gmx_omp_nthreads.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/forcerec.h" // only for GET_CGINFO_*
64 #include "gromacs/pbcutil/ishift.h"
65 #include "gromacs/topology/ifunc.h"
66 #include "gromacs/topology/mtop_lookup.h"
67 #include "gromacs/utility/exceptions.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/gmxassert.h"
70 #include "gromacs/utility/smalloc.h"
72 #include "domdec_internal.h"
73 #include "domdec_specatomcomm.h"
75 /*! \brief Struct used during constraint setup with domain decomposition */
76 struct gmx_domdec_constraints_t
78 //! @cond Doxygen_Suppress
79 std::vector
<int> molb_con_offset
; /**< Offset in the constraint array for each molblock */
80 std::vector
<int> molb_ncon_mol
; /**< The number of constraints per molecule for each molblock */
82 int ncon
; /**< The fully local and conneced constraints */
83 /* The global constraint number, only required for clearing gc_req */
84 std::vector
<int> con_gl
; /**< Global constraint indices for local constraints */
85 std::vector
<int> con_nlocat
; /**< Number of local atoms (2/1/0) for each constraint */
87 std::vector
<bool> gc_req
; /**< Boolean that tells if a global constraint index has been requested; note: size global #constraints */
89 /* Hash table for keeping track of requests */
90 std::unique_ptr
< gmx::HashedMap
< int>> ga2la
; /**< Global to local communicated constraint atom only index */
92 /* Multi-threading stuff */
93 int nthread
; /**< Number of threads used for DD constraint setup */
94 std::vector
<t_ilist
> ils
; /**< Constraint ilist working arrays, size \p nthread */
96 /* Buffers for requesting atoms */
97 std::vector
< std::vector
< int>> requestedGlobalAtomIndices
; /**< Buffers for requesting global atom indices, one per thread */
102 void dd_move_x_constraints(gmx_domdec_t
*dd
, matrix box
,
103 rvec
*x0
, rvec
*x1
, gmx_bool bX1IsCoord
)
105 if (dd
->constraint_comm
)
107 dd_move_x_specat(dd
, dd
->constraint_comm
, box
, x0
, x1
, bX1IsCoord
);
109 ddReopenBalanceRegionCpu(dd
);
113 gmx::ArrayRef
<const int> dd_constraints_nlocalatoms(const gmx_domdec_t
*dd
)
115 if (dd
&& dd
->constraints
)
117 return dd
->constraints
->con_nlocat
;
125 void dd_clear_local_constraint_indices(gmx_domdec_t
*dd
)
127 gmx_domdec_constraints_t
*dc
= dd
->constraints
;
129 std::fill(dc
->gc_req
.begin(), dc
->gc_req
.end(), false);
131 if (dd
->constraint_comm
)
137 /*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */
138 static void walk_out(int con
, int con_offset
, int a
, int offset
, int nrec
,
139 gmx::ArrayRef
<const int> ia1
,
140 gmx::ArrayRef
<const int> ia2
,
141 const t_blocka
*at2con
,
142 const gmx_ga2la_t
&ga2la
, gmx_bool bHomeConnect
,
143 gmx_domdec_constraints_t
*dc
,
144 gmx_domdec_specat_comm_t
*dcc
,
146 std::vector
<int> *ireq
)
148 int a1_gl
, a2_gl
, i
, coni
, b
;
151 if (!dc
->gc_req
[con_offset
+ con
])
153 /* Add this non-home constraint to the list */
154 dc
->con_gl
.push_back(con_offset
+ con
);
155 dc
->con_nlocat
.push_back(bHomeConnect
? 1 : 0);
156 dc
->gc_req
[con_offset
+ con
] = true;
157 if (il_local
->nr
+ 3 > il_local
->nalloc
)
159 il_local
->nalloc
= over_alloc_dd(il_local
->nr
+3);
160 srenew(il_local
->iatoms
, il_local
->nalloc
);
162 iap
= constr_iatomptr(ia1
, ia2
, con
);
163 il_local
->iatoms
[il_local
->nr
++] = iap
[0];
164 a1_gl
= offset
+ iap
[1];
165 a2_gl
= offset
+ iap
[2];
166 /* The following indexing code can probably be optizimed */
167 if (const int *a_loc
= ga2la
.findHome(a1_gl
))
169 il_local
->iatoms
[il_local
->nr
++] = *a_loc
;
173 /* We set this index later */
174 il_local
->iatoms
[il_local
->nr
++] = -a1_gl
- 1;
176 if (const int *a_loc
= ga2la
.findHome(a2_gl
))
178 il_local
->iatoms
[il_local
->nr
++] = *a_loc
;
182 /* We set this index later */
183 il_local
->iatoms
[il_local
->nr
++] = -a2_gl
- 1;
187 /* Check to not ask for the same atom more than once */
188 if (!dc
->ga2la
->find(offset
+ a
))
191 /* Add this non-home atom to the list */
192 ireq
->push_back(offset
+ a
);
193 /* Temporarily mark with -2, we get the index later */
194 dc
->ga2la
->insert(offset
+ a
, -2);
199 for (i
= at2con
->index
[a
]; i
< at2con
->index
[a
+1]; i
++)
205 iap
= constr_iatomptr(ia1
, ia2
, coni
);
214 if (!ga2la
.findHome(offset
+ b
))
216 walk_out(coni
, con_offset
, b
, offset
, nrec
-1,
218 ga2la
, FALSE
, dc
, dcc
, il_local
, ireq
);
225 /*! \brief Looks up SETTLE constraints for a range of charge-groups */
226 static void atoms_to_settles(gmx_domdec_t
*dd
,
227 const gmx_mtop_t
*mtop
,
229 gmx::ArrayRef
< const std::vector
< int>> at2settle_mt
,
230 int cg_start
, int cg_end
,
232 std::vector
<int> *ireq
)
234 const gmx_ga2la_t
&ga2la
= *dd
->ga2la
;
235 int nral
= NRAL(F_SETTLE
);
238 for (int a
= cg_start
; a
< cg_end
; a
++)
240 if (GET_CGINFO_SETTLE(cginfo
[a
]))
242 int a_gl
= dd
->globalAtomIndices
[a
];
244 mtopGetMolblockIndex(mtop
, a_gl
, &mb
, nullptr, &a_mol
);
246 const gmx_molblock_t
*molb
= &mtop
->molblock
[mb
];
247 int settle
= at2settle_mt
[molb
->type
][a_mol
];
251 int offset
= a_gl
- a_mol
;
253 const int *ia1
= mtop
->moltype
[molb
->type
].ilist
[F_SETTLE
].iatoms
.data();
256 gmx_bool bAssign
= FALSE
;
258 for (int sa
= 0; sa
< nral
; sa
++)
260 int a_glsa
= offset
+ ia1
[settle
*(1+nral
)+1+sa
];
262 if (ga2la
.findHome(a_glsa
))
264 if (nlocal
== 0 && a_gl
== a_glsa
)
274 if (ils_local
->nr
+1+nral
> ils_local
->nalloc
)
276 ils_local
->nalloc
= over_alloc_dd(ils_local
->nr
+1+nral
);
277 srenew(ils_local
->iatoms
, ils_local
->nalloc
);
280 ils_local
->iatoms
[ils_local
->nr
++] = ia1
[settle
*4];
282 for (int sa
= 0; sa
< nral
; sa
++)
284 if (const int *a_loc
= ga2la
.findHome(a_gls
[sa
]))
286 ils_local
->iatoms
[ils_local
->nr
++] = *a_loc
;
290 ils_local
->iatoms
[ils_local
->nr
++] = -a_gls
[sa
] - 1;
291 /* Add this non-home atom to the list */
292 ireq
->push_back(a_gls
[sa
]);
293 /* A check on double atom requests is
294 * not required for settle.
304 /*! \brief Looks up constraint for the local atoms */
305 static void atoms_to_constraints(gmx_domdec_t
*dd
,
306 const gmx_mtop_t
*mtop
,
308 gmx::ArrayRef
<const t_blocka
> at2con_mt
, int nrec
,
310 std::vector
<int> *ireq
)
312 const t_blocka
*at2con
;
313 int b_lo
, offset
, b_mol
, i
, con
, con_offset
;
315 gmx_domdec_constraints_t
*dc
= dd
->constraints
;
316 gmx_domdec_specat_comm_t
*dcc
= dd
->constraint_comm
;
318 const gmx_ga2la_t
&ga2la
= *dd
->ga2la
;
321 dc
->con_nlocat
.clear();
325 for (int a
= 0; a
< dd
->ncg_home
; a
++)
327 if (GET_CGINFO_CONSTR(cginfo
[a
]))
329 int a_gl
= dd
->globalAtomIndices
[a
];
331 mtopGetMolblockIndex(mtop
, a_gl
, &mb
, &molnr
, &a_mol
);
333 const gmx_molblock_t
&molb
= mtop
->molblock
[mb
];
335 gmx::ArrayRef
<const int> ia1
= mtop
->moltype
[molb
.type
].ilist
[F_CONSTR
].iatoms
;
336 gmx::ArrayRef
<const int> ia2
= mtop
->moltype
[molb
.type
].ilist
[F_CONSTRNC
].iatoms
;
338 /* Calculate the global constraint number offset for the molecule.
339 * This is only required for the global index to make sure
340 * that we use each constraint only once.
343 dc
->molb_con_offset
[mb
] + molnr
*dc
->molb_ncon_mol
[mb
];
345 /* The global atom number offset for this molecule */
346 offset
= a_gl
- a_mol
;
347 at2con
= &at2con_mt
[molb
.type
];
348 for (i
= at2con
->index
[a_mol
]; i
< at2con
->index
[a_mol
+1]; i
++)
351 const int *iap
= constr_iatomptr(ia1
, ia2
, con
);
360 if (const int *a_loc
= ga2la
.findHome(offset
+ b_mol
))
362 /* Add this fully home constraint at the first atom */
365 dc
->con_gl
.push_back(con_offset
+ con
);
366 dc
->con_nlocat
.push_back(2);
367 if (ilc_local
->nr
+ 3 > ilc_local
->nalloc
)
369 ilc_local
->nalloc
= over_alloc_dd(ilc_local
->nr
+ 3);
370 srenew(ilc_local
->iatoms
, ilc_local
->nalloc
);
373 ilc_local
->iatoms
[ilc_local
->nr
++] = iap
[0];
374 ilc_local
->iatoms
[ilc_local
->nr
++] = (a_gl
== iap
[1] ? a
: b_lo
);
375 ilc_local
->iatoms
[ilc_local
->nr
++] = (a_gl
== iap
[1] ? b_lo
: a
);
382 /* We need the nrec constraints coupled to this constraint,
383 * so we need to walk out of the home cell by nrec+1 atoms,
384 * since already atom bg is not locally present.
385 * Therefore we call walk_out with nrec recursions to go
386 * after this first call.
388 walk_out(con
, con_offset
, b_mol
, offset
, nrec
,
390 ga2la
, TRUE
, dc
, dcc
, ilc_local
, ireq
);
396 GMX_ASSERT(dc
->con_gl
.size() == static_cast<size_t>(dc
->ncon
), "con_gl size should match the number of constraints");
397 GMX_ASSERT(dc
->con_nlocat
.size() == static_cast<size_t>(dc
->ncon
), "con_nlocat size should match the number of constraints");
402 "Constraints: home %3d border %3d atoms: %3zu\n",
403 nhome
, dc
->ncon
-nhome
,
404 dd
->constraint_comm
? ireq
->size() : 0);
408 int dd_make_local_constraints(gmx_domdec_t
*dd
, int at_start
,
409 const struct gmx_mtop_t
*mtop
,
411 gmx::Constraints
*constr
, int nrec
,
414 gmx_domdec_constraints_t
*dc
;
415 t_ilist
*ilc_local
, *ils_local
;
416 std::vector
<int> *ireq
;
417 gmx::ArrayRef
<const t_blocka
> at2con_mt
;
418 gmx::HashedMap
<int> *ga2la_specat
;
422 // This code should not be called unless this condition is true,
423 // because that's the only time init_domdec_constraints is
425 GMX_RELEASE_ASSERT(dd
->comm
->systemInfo
.haveSplitConstraints
||
426 dd
->comm
->systemInfo
.haveSplitSettles
,
427 "dd_make_local_constraints called when there are no local constraints");
428 // ... and init_domdec_constraints always sets
429 // dd->constraint_comm...
430 GMX_RELEASE_ASSERT(dd
->constraint_comm
, "Invalid use of dd_make_local_constraints before construction of constraint_comm");
431 // ... which static analysis needs to be reassured about, because
432 // otherwise, when dd->splitSettles is
433 // true. dd->constraint_comm is unilaterally dereferenced before
434 // the call to atoms_to_settles.
436 dc
= dd
->constraints
;
438 ilc_local
= &il_local
[F_CONSTR
];
439 ils_local
= &il_local
[F_SETTLE
];
443 if (dd
->constraint_comm
)
445 // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
446 GMX_RELEASE_ASSERT(constr
!= nullptr, "Must have valid constraints object");
447 at2con_mt
= constr
->atom2constraints_moltype();
448 ireq
= &dc
->requestedGlobalAtomIndices
[0];
453 // Currently unreachable
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
,
477 /* Do the constraints, if present, on the first thread.
478 * Do the settles on all other threads.
480 t0_set
= ((!at2con_mt
.empty() && dc
->nthread
> 1) ? 1 : 0);
482 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
483 for (int thread
= 0; thread
< dc
->nthread
; thread
++)
487 if (!at2con_mt
.empty() && thread
== 0)
489 atoms_to_constraints(dd
, mtop
, cginfo
, at2con_mt
, nrec
,
493 if (thread
>= t0_set
)
498 /* Distribute the settle check+assignments over
499 * dc->nthread or dc->nthread-1 threads.
501 cg0
= (dd
->ncg_home
*(thread
-t0_set
))/(dc
->nthread
-t0_set
);
502 cg1
= (dd
->ncg_home
*(thread
-t0_set
+1))/(dc
->nthread
-t0_set
);
504 if (thread
== t0_set
)
510 ilst
= &dc
->ils
[thread
];
514 std::vector
<int> &ireqt
= dc
->requestedGlobalAtomIndices
[thread
];
520 atoms_to_settles(dd
, mtop
, cginfo
, at2settle_mt
,
525 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
528 /* Combine the generate settles and requested indices */
529 for (int thread
= 1; thread
< dc
->nthread
; thread
++)
536 ilst
= &dc
->ils
[thread
];
537 if (ils_local
->nr
+ ilst
->nr
> ils_local
->nalloc
)
539 ils_local
->nalloc
= over_alloc_large(ils_local
->nr
+ ilst
->nr
);
540 srenew(ils_local
->iatoms
, ils_local
->nalloc
);
542 for (ia
= 0; ia
< ilst
->nr
; ia
++)
544 ils_local
->iatoms
[ils_local
->nr
+ia
] = ilst
->iatoms
[ia
];
546 ils_local
->nr
+= ilst
->nr
;
549 const std::vector
<int> &ireqt
= dc
->requestedGlobalAtomIndices
[thread
];
550 ireq
->insert(ireq
->end(), ireqt
.begin(), ireqt
.end());
555 fprintf(debug
, "Settles: total %3d\n", ils_local
->nr
/4);
559 if (dd
->constraint_comm
)
564 setup_specat_communication(dd
, ireq
, dd
->constraint_comm
,
565 dd
->constraints
->ga2la
.get(),
567 "constraint", " or lincs-order");
569 /* Fill in the missing indices */
570 ga2la_specat
= dd
->constraints
->ga2la
.get();
572 nral1
= 1 + NRAL(F_CONSTR
);
573 for (i
= 0; i
< ilc_local
->nr
; i
+= nral1
)
575 iap
= ilc_local
->iatoms
+ i
;
576 for (j
= 1; j
< nral1
; j
++)
580 const int *a
= ga2la_specat
->find(-iap
[j
] - 1);
581 GMX_ASSERT(a
, "We have checked before that this atom index has been set");
587 nral1
= 1 + NRAL(F_SETTLE
);
588 for (i
= 0; i
< ils_local
->nr
; i
+= nral1
)
590 iap
= ils_local
->iatoms
+ i
;
591 for (j
= 1; j
< nral1
; j
++)
595 const int *a
= ga2la_specat
->find(-iap
[j
] - 1);
596 GMX_ASSERT(a
, "We have checked before that this atom index has been set");
604 // Currently unreachable
611 void init_domdec_constraints(gmx_domdec_t
*dd
,
612 const gmx_mtop_t
*mtop
)
614 gmx_domdec_constraints_t
*dc
;
615 const gmx_molblock_t
*molb
;
619 fprintf(debug
, "Begin init_domdec_constraints\n");
622 dd
->constraints
= new gmx_domdec_constraints_t
;
623 dc
= dd
->constraints
;
625 dc
->molb_con_offset
.resize(mtop
->molblock
.size());
626 dc
->molb_ncon_mol
.resize(mtop
->molblock
.size());
629 for (size_t mb
= 0; mb
< mtop
->molblock
.size(); mb
++)
631 molb
= &mtop
->molblock
[mb
];
632 dc
->molb_con_offset
[mb
] = ncon
;
633 dc
->molb_ncon_mol
[mb
] =
634 mtop
->moltype
[molb
->type
].ilist
[F_CONSTR
].size()/3 +
635 mtop
->moltype
[molb
->type
].ilist
[F_CONSTRNC
].size()/3;
636 ncon
+= molb
->nmol
*dc
->molb_ncon_mol
[mb
];
641 dc
->gc_req
.resize(ncon
);
644 /* Use a hash table for the global to local index.
645 * The number of keys is a rough estimate, it will be optimized later.
647 int numKeysEstimate
= std::min(mtop
->natoms
/20,
648 mtop
->natoms
/(2*dd
->nnodes
));
649 dc
->ga2la
= std::make_unique
< gmx::HashedMap
< int>>(numKeysEstimate
);
651 dc
->nthread
= gmx_omp_nthreads_get(emntDomdec
);
652 dc
->ils
.resize(dc
->nthread
);
654 dd
->constraint_comm
= new gmx_domdec_specat_comm_t
;
656 dc
->requestedGlobalAtomIndices
.resize(dc
->nthread
);