2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2006 - 2014, The GROMACS development team.
5 * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
39 * \brief This file defines functions used by the domdec module
40 * while managing the construction, use and error checking for
41 * topologies local to a DD rank.
43 * \author Berk Hess <hess@kth.se>
44 * \ingroup module_domdec
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/domdec/domdec_network.h"
60 #include "gromacs/domdec/ga2la.h"
61 #include "gromacs/domdec/localtopologychecker.h"
62 #include "gromacs/domdec/options.h"
63 #include "gromacs/domdec/reversetopology.h"
64 #include "gromacs/gmxlib/network.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/mdlib/forcerec.h"
67 #include "gromacs/mdlib/gmx_omp_nthreads.h"
68 #include "gromacs/mdlib/vsite.h"
69 #include "gromacs/mdtypes/commrec.h"
70 #include "gromacs/mdtypes/forcerec.h"
71 #include "gromacs/mdtypes/inputrec.h"
72 #include "gromacs/mdtypes/md_enums.h"
73 #include "gromacs/mdtypes/mdatom.h"
74 #include "gromacs/mdtypes/state.h"
75 #include "gromacs/pbcutil/mshift.h"
76 #include "gromacs/pbcutil/pbc.h"
77 #include "gromacs/topology/mtop_util.h"
78 #include "gromacs/topology/topsort.h"
79 #include "gromacs/utility/cstringutil.h"
80 #include "gromacs/utility/exceptions.h"
81 #include "gromacs/utility/fatalerror.h"
82 #include "gromacs/utility/gmxassert.h"
83 #include "gromacs/utility/logger.h"
84 #include "gromacs/utility/smalloc.h"
85 #include "gromacs/utility/strconvert.h"
86 #include "gromacs/utility/stringstream.h"
87 #include "gromacs/utility/stringutil.h"
88 #include "gromacs/utility/textwriter.h"
90 #include "domdec_constraints.h"
91 #include "domdec_internal.h"
92 #include "domdec_vsite.h"
96 using gmx::DDBondedChecking;
99 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
100 #define NITEM_DD_INIT_LOCAL_STATE 5
102 void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local)
104 int buf[NITEM_DD_INIT_LOCAL_STATE];
108 buf[0] = state_global->flags;
109 buf[1] = state_global->ngtc;
110 buf[2] = state_global->nnhpres;
111 buf[3] = state_global->nhchainlength;
112 buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
114 dd_bcast(&dd, NITEM_DD_INIT_LOCAL_STATE * sizeof(int), buf);
116 init_gtc_state(state_local, buf[1], buf[2], buf[3]);
117 init_dfhist_state(state_local, buf[4]);
118 state_local->flags = buf[0];
121 /*! \brief Check if a link is stored in \p link between charge groups \p cg_gl and \p cg_gl_j and if not so, store a link */
122 static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
125 for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
127 GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
128 if (link->a[k] == cg_gl_j)
135 GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
136 "Inconsistent allocation of link");
137 /* Add this charge group link */
138 if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
140 link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
141 srenew(link->a, link->nalloc_a);
143 link->a[link->index[cg_gl + 1]] = cg_gl_j;
144 link->index[cg_gl + 1]++;
148 void makeBondedLinks(gmx_domdec_t* dd,
149 const gmx_mtop_t& mtop,
150 gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock)
153 if (!dd->comm->systemInfo.filterBondedCommunication)
155 /* Only communicate atoms based on cut-off */
156 dd->comm->bondedLinks = nullptr;
160 t_blocka* link = nullptr;
162 /* For each atom make a list of other atoms in the system
163 * that a linked to it via bonded interactions
164 * which are also stored in reverse_top.
167 reverse_ilist_t ril_intermol;
168 if (mtop.bIntermolecularInteractions)
172 atoms.nr = mtop.natoms;
173 atoms.atom = nullptr;
175 GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
176 "We should have an ilist when intermolecular interactions are on");
178 ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
180 *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
184 snew(link->index, mtop.natoms + 1);
189 int indexOfFirstAtomInMolecule = 0;
190 int numLinkedAtoms = 0;
191 for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
193 const gmx_molblock_t& molb = mtop.molblock[mb];
198 const gmx_moltype_t& molt = mtop.moltype[molb.type];
199 /* Make a reverse ilist in which the interactions are linked
200 * to all atoms, not only the first atom as in gmx_reverse_top.
201 * The constraints are discarded here.
203 ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
205 make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
207 AtomInfoWithinMoleculeBlock* atomInfoOfMoleculeBlock = &atomInfoForEachMoleculeBlock[mb];
210 for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
212 for (int a = 0; a < molt.atoms.nr; a++)
214 int atomIndex = indexOfFirstAtomInMolecule + a;
215 link->index[atomIndex + 1] = link->index[atomIndex];
216 int i = ril.index[a];
217 while (i < ril.index[a + 1])
219 int ftype = ril.il[i++];
220 int nral = NRAL(ftype);
221 /* Skip the ifunc index */
223 for (int j = 0; j < nral; j++)
225 int aj = ril.il[i + j];
228 check_link(link, atomIndex, indexOfFirstAtomInMolecule + aj);
234 if (mtop.bIntermolecularInteractions)
236 int i = ril_intermol.index[atomIndex];
237 while (i < ril_intermol.index[atomIndex + 1])
239 int ftype = ril_intermol.il[i++];
240 int nral = NRAL(ftype);
241 /* Skip the ifunc index */
243 for (int j = 0; j < nral; j++)
245 /* Here we assume we have no charge groups;
246 * this has been checked above.
248 int aj = ril_intermol.il[i + j];
249 check_link(link, atomIndex, aj);
254 if (link->index[atomIndex + 1] - link->index[atomIndex] > 0)
256 SET_CGINFO_BOND_INTER(atomInfoOfMoleculeBlock->atomInfo[a]);
261 indexOfFirstAtomInMolecule += molt.atoms.nr;
263 int nlink_mol = link->index[indexOfFirstAtomInMolecule]
264 - link->index[indexOfFirstAtomInMolecule - molt.atoms.nr];
269 "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
277 /* Copy the data for the rest of the molecules in this block */
278 link->nalloc_a += (molb.nmol - mol) * nlink_mol;
279 srenew(link->a, link->nalloc_a);
280 for (; mol < molb.nmol; mol++)
282 for (int a = 0; a < molt.atoms.nr; a++)
284 int atomIndex = indexOfFirstAtomInMolecule + a;
285 link->index[atomIndex + 1] = link->index[atomIndex + 1 - molt.atoms.nr] + nlink_mol;
286 for (int j = link->index[atomIndex]; j < link->index[atomIndex + 1]; j++)
288 link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
290 if (link->index[atomIndex + 1] - link->index[atomIndex] > 0
291 && atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock
292 < gmx::ssize(atomInfoOfMoleculeBlock->atomInfo))
294 SET_CGINFO_BOND_INTER(
295 atomInfoOfMoleculeBlock
296 ->atomInfo[atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock]);
300 indexOfFirstAtomInMolecule += molt.atoms.nr;
307 fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, numLinkedAtoms);
310 dd->comm->bondedLinks = link;
321 /*! \brief Compare distance^2 \p r2 against the distance in \p bd and if larger store it along with \p ftype and atom indices \p a1 and \p a2 */
322 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, bonded_distance_t* bd)
333 /*! \brief Set the distance, function type and atom indices for the longest distance between atoms of molecule type \p molt for two-body and multi-body bonded interactions */
334 static void bonded_cg_distance_mol(const gmx_moltype_t* molt,
335 const DDBondedChecking ddBondedChecking,
337 ArrayRef<const RVec> x,
338 bonded_distance_t* bd_2b,
339 bonded_distance_t* bd_mb)
341 const ReverseTopOptions rtOptions(ddBondedChecking);
343 for (int ftype = 0; ftype < F_NRE; ftype++)
345 if (dd_check_ftype(ftype, rtOptions))
347 const auto& il = molt->ilist[ftype];
348 int nral = NRAL(ftype);
351 for (int i = 0; i < il.size(); i += 1 + nral)
353 for (int ai = 0; ai < nral; ai++)
355 int atomI = il.iatoms[i + 1 + ai];
356 for (int aj = ai + 1; aj < nral; aj++)
358 int atomJ = il.iatoms[i + 1 + aj];
361 real rij2 = distance2(x[atomI], x[atomJ]);
363 update_max_bonded_distance(
364 rij2, ftype, atomI, atomJ, (nral == 2) ? bd_2b : bd_mb);
374 const auto& excls = molt->excls;
375 for (gmx::index ai = 0; ai < excls.ssize(); ai++)
377 for (const int aj : excls[ai])
381 real rij2 = distance2(x[ai], x[aj]);
383 /* There is no function type for exclusions, use -1 */
384 update_max_bonded_distance(rij2, -1, ai, aj, bd_2b);
391 /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
392 static void bonded_distance_intermol(const InteractionLists& ilists_intermol,
393 const DDBondedChecking ddBondedChecking,
394 ArrayRef<const RVec> x,
397 bonded_distance_t* bd_2b,
398 bonded_distance_t* bd_mb)
402 set_pbc(&pbc, pbcType, box);
404 const ReverseTopOptions rtOptions(ddBondedChecking);
406 for (int ftype = 0; ftype < F_NRE; ftype++)
408 if (dd_check_ftype(ftype, rtOptions))
410 const auto& il = ilists_intermol[ftype];
411 int nral = NRAL(ftype);
413 /* No nral>1 check here, since intermol interactions always
414 * have nral>=2 (and the code is also correct for nral=1).
416 for (int i = 0; i < il.size(); i += 1 + nral)
418 for (int ai = 0; ai < nral; ai++)
420 int atom_i = il.iatoms[i + 1 + ai];
422 for (int aj = ai + 1; aj < nral; aj++)
426 int atom_j = il.iatoms[i + 1 + aj];
428 pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
430 const real rij2 = norm2(dx);
432 update_max_bonded_distance(rij2, ftype, atom_i, atom_j, (nral == 2) ? bd_2b : bd_mb);
440 //! Returns whether \p molt has at least one virtual site
441 static bool moltypeHasVsite(const gmx_moltype_t& molt)
443 bool hasVsite = false;
444 for (int i = 0; i < F_NRE; i++)
446 if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
455 //! Returns coordinates not broken over PBC for a molecule
456 static void getWholeMoleculeCoordinates(const gmx_moltype_t* molt,
457 const gmx_ffparams_t* ffparams,
461 ArrayRef<const RVec> x,
464 if (pbcType != PbcType::No)
466 mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(x.data()));
468 shift_x(graph, box, as_rvec_array(x.data()), as_rvec_array(xs.data()));
469 /* By doing an extra mk_mshift the molecules that are broken
470 * because they were e.g. imported from another software
471 * will be made whole again. Such are the healing powers
474 mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(xs.data()));
478 /* We copy the coordinates so the original coordinates remain
479 * unchanged, just to be 100% sure that we do not affect
480 * binary reproducibility of simulations.
482 for (int i = 0; i < molt->atoms.nr; i++)
484 copy_rvec(x[i], xs[i]);
488 if (moltypeHasVsite(*molt))
490 gmx::constructVirtualSites(xs, ffparams->iparams, molt->ilist);
494 void dd_bonded_cg_distance(const gmx::MDLogger& mdlog,
495 const gmx_mtop_t& mtop,
496 const t_inputrec& inputrec,
497 ArrayRef<const RVec> x,
499 const DDBondedChecking ddBondedChecking,
503 bonded_distance_t bd_2b = { 0, -1, -1, -1 };
504 bonded_distance_t bd_mb = { 0, -1, -1, -1 };
506 bool bExclRequired = inputrecExclForces(&inputrec);
511 for (const gmx_molblock_t& molb : mtop.molblock)
513 const gmx_moltype_t& molt = mtop.moltype[molb.type];
514 if (molt.atoms.nr == 1 || molb.nmol == 0)
516 at_offset += molb.nmol * molt.atoms.nr;
521 if (inputrec.pbcType != PbcType::No)
523 graph = mk_graph_moltype(molt);
526 std::vector<RVec> xs(molt.atoms.nr);
527 for (int mol = 0; mol < molb.nmol; mol++)
529 getWholeMoleculeCoordinates(&molt,
534 x.subArray(at_offset, molt.atoms.nr),
537 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
538 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
540 bonded_cg_distance_mol(&molt, ddBondedChecking, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb);
542 /* Process the mol data adding the atom index offset */
543 update_max_bonded_distance(bd_mol_2b.r2,
545 at_offset + bd_mol_2b.a1,
546 at_offset + bd_mol_2b.a2,
548 update_max_bonded_distance(bd_mol_mb.r2,
550 at_offset + bd_mol_mb.a1,
551 at_offset + bd_mol_mb.a2,
554 at_offset += molt.atoms.nr;
559 if (mtop.bIntermolecularInteractions)
561 GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
562 "We should have an ilist when intermolecular interactions are on");
564 bonded_distance_intermol(
565 *mtop.intermolecular_ilist, ddBondedChecking, x, inputrec.pbcType, box, &bd_2b, &bd_mb);
568 *r_2b = sqrt(bd_2b.r2);
569 *r_mb = sqrt(bd_mb.r2);
571 if (*r_2b > 0 || *r_mb > 0)
573 GMX_LOG(mdlog.info).appendText("Initial maximum distances in bonded interactions:");
577 .appendTextFormatted(
578 " two-body bonded interactions: %5.3f nm, %s, atoms %d %d",
580 (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
587 .appendTextFormatted(
588 " multi-body bonded interactions: %5.3f nm, %s, atoms %d %d",
590 interaction_function[bd_mb.ftype].longname,