Finish breaking up domdec_topology.cpp
authorMark Abraham <mark.j.abraham@gmail.com>
Mon, 31 May 2021 15:32:14 +0000 (15:32 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Mon, 31 May 2021 15:32:14 +0000 (15:32 +0000)
src/gromacs/domdec/computemultibodycutoffs.cpp [moved from src/gromacs/domdec/domdec_topology.cpp with 57% similarity]
src/gromacs/domdec/computemultibodycutoffs.h [new file with mode: 0644]
src/gromacs/domdec/domdec.cpp
src/gromacs/domdec/domdec.h
src/gromacs/domdec/makebondedlinks.cpp [new file with mode: 0644]
src/gromacs/domdec/makebondedlinks.h [new file with mode: 0644]
src/gromacs/mdrun/runner.cpp

similarity index 57%
rename from src/gromacs/domdec/domdec_topology.cpp
rename to src/gromacs/domdec/computemultibodycutoffs.cpp
index 289a0598712a806d1b45c75f51502a022c944e0c..24c63e8f7167ece696d57431725ad1d37ff01d76 100644 (file)
@@ -1,8 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2006 - 2014, The GROMACS development team.
- * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
+ * Copyright (c) 2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 
 #include "gmxpre.h"
 
-#include <cassert>
-#include <cstdlib>
-#include <cstring>
+#include "gromacs/domdec/computemultibodycutoffs.h"
 
-#include <algorithm>
-#include <memory>
-#include <optional>
-#include <string>
+#include <vector>
 
-#include "gromacs/domdec/domdec.h"
-#include "gromacs/domdec/domdec_network.h"
-#include "gromacs/domdec/ga2la.h"
-#include "gromacs/domdec/localtopologychecker.h"
 #include "gromacs/domdec/options.h"
 #include "gromacs/domdec/reversetopology.h"
-#include "gromacs/gmxlib/network.h"
 #include "gromacs/math/vec.h"
-#include "gromacs/mdlib/forcerec.h"
-#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/vsite.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/mshift.h"
 #include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/mtop_util.h"
-#include "gromacs/topology/topsort.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/exceptions.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/logger.h"
-#include "gromacs/utility/smalloc.h"
-#include "gromacs/utility/strconvert.h"
-#include "gromacs/utility/stringstream.h"
-#include "gromacs/utility/stringutil.h"
-#include "gromacs/utility/textwriter.h"
-
-#include "domdec_constraints.h"
-#include "domdec_internal.h"
-#include "domdec_vsite.h"
-#include "dump.h"
+#include "gromacs/topology/mtop_util.h"
 
 using gmx::ArrayRef;
 using gmx::DDBondedChecking;
 using gmx::RVec;
 
-/*! \brief The number of integer item in the local state, used for broadcasting of the state */
-#define NITEM_DD_INIT_LOCAL_STATE 5
-
-void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local)
-{
-    int buf[NITEM_DD_INIT_LOCAL_STATE];
-
-    if (DDMASTER(dd))
-    {
-        buf[0] = state_global->flags;
-        buf[1] = state_global->ngtc;
-        buf[2] = state_global->nnhpres;
-        buf[3] = state_global->nhchainlength;
-        buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
-    }
-    dd_bcast(&dd, NITEM_DD_INIT_LOCAL_STATE * sizeof(int), buf);
-
-    init_gtc_state(state_local, buf[1], buf[2], buf[3]);
-    init_dfhist_state(state_local, buf[4]);
-    state_local->flags = buf[0];
-}
-
-/*! \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 */
-static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
-{
-    bool bFound = false;
-    for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
-    {
-        GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
-        if (link->a[k] == cg_gl_j)
-        {
-            bFound = TRUE;
-        }
-    }
-    if (!bFound)
-    {
-        GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
-                           "Inconsistent allocation of link");
-        /* Add this charge group link */
-        if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
-        {
-            link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
-            srenew(link->a, link->nalloc_a);
-        }
-        link->a[link->index[cg_gl + 1]] = cg_gl_j;
-        link->index[cg_gl + 1]++;
-    }
-}
-
-void makeBondedLinks(gmx_domdec_t*                              dd,
-                     const gmx_mtop_t&                          mtop,
-                     gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock)
-{
-
-    if (!dd->comm->systemInfo.filterBondedCommunication)
-    {
-        /* Only communicate atoms based on cut-off */
-        dd->comm->bondedLinks = nullptr;
-        return;
-    }
-
-    t_blocka* link = nullptr;
-
-    /* For each atom make a list of other atoms in the system
-     * that a linked to it via bonded interactions
-     * which are also stored in reverse_top.
-     */
-
-    reverse_ilist_t ril_intermol;
-    if (mtop.bIntermolecularInteractions)
-    {
-        t_atoms atoms;
-
-        atoms.nr   = mtop.natoms;
-        atoms.atom = nullptr;
-
-        GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
-                           "We should have an ilist when intermolecular interactions are on");
-
-        ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
-        make_reverse_ilist(
-                *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
-    }
-
-    snew(link, 1);
-    snew(link->index, mtop.natoms + 1);
-    link->nalloc_a = 0;
-    link->a        = nullptr;
-
-    link->index[0]                 = 0;
-    int indexOfFirstAtomInMolecule = 0;
-    int numLinkedAtoms             = 0;
-    for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
-    {
-        const gmx_molblock_t& molb = mtop.molblock[mb];
-        if (molb.nmol == 0)
-        {
-            continue;
-        }
-        const gmx_moltype_t& molt = mtop.moltype[molb.type];
-        /* Make a reverse ilist in which the interactions are linked
-         * to all atoms, not only the first atom as in gmx_reverse_top.
-         * The constraints are discarded here.
-         */
-        ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
-        reverse_ilist_t   ril;
-        make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
-
-        AtomInfoWithinMoleculeBlock* atomInfoOfMoleculeBlock = &atomInfoForEachMoleculeBlock[mb];
-
-        int mol = 0;
-        for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
-        {
-            for (int a = 0; a < molt.atoms.nr; a++)
-            {
-                int atomIndex              = indexOfFirstAtomInMolecule + a;
-                link->index[atomIndex + 1] = link->index[atomIndex];
-                int i                      = ril.index[a];
-                while (i < ril.index[a + 1])
-                {
-                    int ftype = ril.il[i++];
-                    int nral  = NRAL(ftype);
-                    /* Skip the ifunc index */
-                    i++;
-                    for (int j = 0; j < nral; j++)
-                    {
-                        int aj = ril.il[i + j];
-                        if (aj != a)
-                        {
-                            check_link(link, atomIndex, indexOfFirstAtomInMolecule + aj);
-                        }
-                    }
-                    i += nral_rt(ftype);
-                }
-
-                if (mtop.bIntermolecularInteractions)
-                {
-                    int i = ril_intermol.index[atomIndex];
-                    while (i < ril_intermol.index[atomIndex + 1])
-                    {
-                        int ftype = ril_intermol.il[i++];
-                        int nral  = NRAL(ftype);
-                        /* Skip the ifunc index */
-                        i++;
-                        for (int j = 0; j < nral; j++)
-                        {
-                            /* Here we assume we have no charge groups;
-                             * this has been checked above.
-                             */
-                            int aj = ril_intermol.il[i + j];
-                            check_link(link, atomIndex, aj);
-                        }
-                        i += nral_rt(ftype);
-                    }
-                }
-                if (link->index[atomIndex + 1] - link->index[atomIndex] > 0)
-                {
-                    SET_CGINFO_BOND_INTER(atomInfoOfMoleculeBlock->atomInfo[a]);
-                    numLinkedAtoms++;
-                }
-            }
-
-            indexOfFirstAtomInMolecule += molt.atoms.nr;
-        }
-        int nlink_mol = link->index[indexOfFirstAtomInMolecule]
-                        - link->index[indexOfFirstAtomInMolecule - molt.atoms.nr];
-
-        if (debug)
-        {
-            fprintf(debug,
-                    "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
-                    *molt.name,
-                    molt.atoms.nr,
-                    nlink_mol);
-        }
-
-        if (molb.nmol > mol)
-        {
-            /* Copy the data for the rest of the molecules in this block */
-            link->nalloc_a += (molb.nmol - mol) * nlink_mol;
-            srenew(link->a, link->nalloc_a);
-            for (; mol < molb.nmol; mol++)
-            {
-                for (int a = 0; a < molt.atoms.nr; a++)
-                {
-                    int atomIndex              = indexOfFirstAtomInMolecule + a;
-                    link->index[atomIndex + 1] = link->index[atomIndex + 1 - molt.atoms.nr] + nlink_mol;
-                    for (int j = link->index[atomIndex]; j < link->index[atomIndex + 1]; j++)
-                    {
-                        link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
-                    }
-                    if (link->index[atomIndex + 1] - link->index[atomIndex] > 0
-                        && atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock
-                                   < gmx::ssize(atomInfoOfMoleculeBlock->atomInfo))
-                    {
-                        SET_CGINFO_BOND_INTER(
-                                atomInfoOfMoleculeBlock
-                                        ->atomInfo[atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock]);
-                        numLinkedAtoms++;
-                    }
-                }
-                indexOfFirstAtomInMolecule += molt.atoms.nr;
-            }
-        }
-    }
-
-    if (debug)
-    {
-        fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, numLinkedAtoms);
-    }
-
-    dd->comm->bondedLinks = link;
-}
-
 typedef struct
 {
     real r2;
diff --git a/src/gromacs/domdec/computemultibodycutoffs.h b/src/gromacs/domdec/computemultibodycutoffs.h
new file mode 100644 (file)
index 0000000..ddc1447
--- /dev/null
@@ -0,0 +1,72 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \libinternal \file
+ *
+ * \brief This file declares the function for computing the required
+ * cutoff distance for inter-domain multi-body interactions, when
+ * those exist.
+ *
+ * \inlibraryapi
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_COMPUTEMULTIBODYCUTOFFS_H
+#define GMX_DOMDEC_COMPUTEMULTIBODYCUTOFFS_H
+
+#include "gromacs/math/vectypes.h"
+
+struct gmx_mtop_t;
+struct t_inputrec;
+
+namespace gmx
+{
+template<typename>
+class ArrayRef;
+class MDLogger;
+enum class DDBondedChecking : bool;
+} // namespace gmx
+
+/*! \brief Calculate the maximum distance involved in 2-body and multi-body bonded interactions */
+void dd_bonded_cg_distance(const gmx::MDLogger&           mdlog,
+                           const gmx_mtop_t&              mtop,
+                           const t_inputrec&              ir,
+                           gmx::ArrayRef<const gmx::RVec> x,
+                           const matrix                   box,
+                           gmx::DDBondedChecking          ddBondedChecking,
+                           real*                          r_2b,
+                           real*                          r_mb);
+
+#endif
index b3532bd189f612e16ce969139cdd161438b381d0..42770b651fc5ad5c2e7fe096fdddf61f48eb9bc3 100644 (file)
@@ -55,6 +55,7 @@
 
 #include "gromacs/domdec/builder.h"
 #include "gromacs/domdec/collect.h"
+#include "gromacs/domdec/computemultibodycutoffs.h"
 #include "gromacs/domdec/dlb.h"
 #include "gromacs/domdec/dlbtiming.h"
 #include "gromacs/domdec/domdec_network.h"
@@ -3217,3 +3218,22 @@ void communicateGpuHaloForces(const t_commrec& cr, bool accumulateForces)
         }
     }
 }
+
+void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local)
+{
+    std::array<int, 5> buf;
+
+    if (DDMASTER(dd))
+    {
+        buf[0] = state_global->flags;
+        buf[1] = state_global->ngtc;
+        buf[2] = state_global->nnhpres;
+        buf[3] = state_global->nhchainlength;
+        buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
+    }
+    dd_bcast(&dd, buf.size() * sizeof(int), buf.data());
+
+    init_gtc_state(state_local, buf[1], buf[2], buf[3]);
+    init_dfhist_state(state_local, buf[4]);
+    state_local->flags = buf[0];
+}
index e6564bf752b844e74f8d40ce8f52a649a83dfd99..301c3c82b64d2d20a2a48eea67f02921e99bb2e5 100644 (file)
@@ -244,24 +244,6 @@ gmx::ArrayRef<const int> dd_constraints_nlocalatoms(const gmx_domdec_t* dd);
 /*! \brief Construct local state */
 void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* local_state);
 
-/*! \brief Generate a list of links between atoms that are linked by bonded interactions
- *
- * Also stores whether atoms are linked in \p atomInfoForEachMoleculeBlock.
- */
-void makeBondedLinks(gmx_domdec_t*                              dd,
-                     const gmx_mtop_t&                          mtop,
-                     gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock);
-
-/*! \brief Calculate the maximum distance involved in 2-body and multi-body bonded interactions */
-void dd_bonded_cg_distance(const gmx::MDLogger&           mdlog,
-                           const gmx_mtop_t&              mtop,
-                           const t_inputrec&              ir,
-                           gmx::ArrayRef<const gmx::RVec> x,
-                           const matrix                   box,
-                           gmx::DDBondedChecking          ddBondedChecking,
-                           real*                          r_2b,
-                           real*                          r_mb);
-
 /*! \brief Construct the GPU halo exchange object(s).
  *
  * \param[in] mdlog               The logger object.
diff --git a/src/gromacs/domdec/makebondedlinks.cpp b/src/gromacs/domdec/makebondedlinks.cpp
new file mode 100644 (file)
index 0000000..265832f
--- /dev/null
@@ -0,0 +1,251 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief Defines a function that makes the list of links between
+ * atoms connected by bonded interactions.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_domdec
+ */
+
+#include "gmxpre.h"
+
+#include "gromacs/domdec/makebondedlinks.h"
+
+#include "gromacs/domdec/domdec_internal.h"
+#include "gromacs/domdec/options.h"
+#include "gromacs/domdec/reversetopology.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/mtop_util.h"
+
+using gmx::ArrayRef;
+using gmx::DDBondedChecking;
+
+/*! \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 */
+static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
+{
+    bool bFound = false;
+    for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
+    {
+        GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
+        if (link->a[k] == cg_gl_j)
+        {
+            bFound = TRUE;
+        }
+    }
+    if (!bFound)
+    {
+        GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
+                           "Inconsistent allocation of link");
+        /* Add this charge group link */
+        if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
+        {
+            link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
+            srenew(link->a, link->nalloc_a);
+        }
+        link->a[link->index[cg_gl + 1]] = cg_gl_j;
+        link->index[cg_gl + 1]++;
+    }
+}
+
+void makeBondedLinks(gmx_domdec_t*                              dd,
+                     const gmx_mtop_t&                          mtop,
+                     gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock)
+{
+
+    if (!dd->comm->systemInfo.filterBondedCommunication)
+    {
+        /* Only communicate atoms based on cut-off */
+        dd->comm->bondedLinks = nullptr;
+        return;
+    }
+
+    t_blocka* link = nullptr;
+
+    /* For each atom make a list of other atoms in the system
+     * that a linked to it via bonded interactions
+     * which are also stored in reverse_top.
+     */
+
+    reverse_ilist_t ril_intermol;
+    if (mtop.bIntermolecularInteractions)
+    {
+        t_atoms atoms;
+
+        atoms.nr   = mtop.natoms;
+        atoms.atom = nullptr;
+
+        GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
+                           "We should have an ilist when intermolecular interactions are on");
+
+        ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
+        make_reverse_ilist(
+                *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
+    }
+
+    snew(link, 1);
+    snew(link->index, mtop.natoms + 1);
+    link->nalloc_a = 0;
+    link->a        = nullptr;
+
+    link->index[0]                 = 0;
+    int indexOfFirstAtomInMolecule = 0;
+    int numLinkedAtoms             = 0;
+    for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
+    {
+        const gmx_molblock_t& molb = mtop.molblock[mb];
+        if (molb.nmol == 0)
+        {
+            continue;
+        }
+        const gmx_moltype_t& molt = mtop.moltype[molb.type];
+        /* Make a reverse ilist in which the interactions are linked
+         * to all atoms, not only the first atom as in gmx_reverse_top.
+         * The constraints are discarded here.
+         */
+        ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
+        reverse_ilist_t   ril;
+        make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
+
+        AtomInfoWithinMoleculeBlock* atomInfoOfMoleculeBlock = &atomInfoForEachMoleculeBlock[mb];
+
+        int mol = 0;
+        for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
+        {
+            for (int a = 0; a < molt.atoms.nr; a++)
+            {
+                int atomIndex              = indexOfFirstAtomInMolecule + a;
+                link->index[atomIndex + 1] = link->index[atomIndex];
+                int i                      = ril.index[a];
+                while (i < ril.index[a + 1])
+                {
+                    int ftype = ril.il[i++];
+                    int nral  = NRAL(ftype);
+                    /* Skip the ifunc index */
+                    i++;
+                    for (int j = 0; j < nral; j++)
+                    {
+                        int aj = ril.il[i + j];
+                        if (aj != a)
+                        {
+                            check_link(link, atomIndex, indexOfFirstAtomInMolecule + aj);
+                        }
+                    }
+                    i += nral_rt(ftype);
+                }
+
+                if (mtop.bIntermolecularInteractions)
+                {
+                    int i = ril_intermol.index[atomIndex];
+                    while (i < ril_intermol.index[atomIndex + 1])
+                    {
+                        int ftype = ril_intermol.il[i++];
+                        int nral  = NRAL(ftype);
+                        /* Skip the ifunc index */
+                        i++;
+                        for (int j = 0; j < nral; j++)
+                        {
+                            /* Here we assume we have no charge groups;
+                             * this has been checked above.
+                             */
+                            int aj = ril_intermol.il[i + j];
+                            check_link(link, atomIndex, aj);
+                        }
+                        i += nral_rt(ftype);
+                    }
+                }
+                if (link->index[atomIndex + 1] - link->index[atomIndex] > 0)
+                {
+                    SET_CGINFO_BOND_INTER(atomInfoOfMoleculeBlock->atomInfo[a]);
+                    numLinkedAtoms++;
+                }
+            }
+
+            indexOfFirstAtomInMolecule += molt.atoms.nr;
+        }
+        int nlink_mol = link->index[indexOfFirstAtomInMolecule]
+                        - link->index[indexOfFirstAtomInMolecule - molt.atoms.nr];
+
+        if (debug)
+        {
+            fprintf(debug,
+                    "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
+                    *molt.name,
+                    molt.atoms.nr,
+                    nlink_mol);
+        }
+
+        if (molb.nmol > mol)
+        {
+            /* Copy the data for the rest of the molecules in this block */
+            link->nalloc_a += (molb.nmol - mol) * nlink_mol;
+            srenew(link->a, link->nalloc_a);
+            for (; mol < molb.nmol; mol++)
+            {
+                for (int a = 0; a < molt.atoms.nr; a++)
+                {
+                    int atomIndex              = indexOfFirstAtomInMolecule + a;
+                    link->index[atomIndex + 1] = link->index[atomIndex + 1 - molt.atoms.nr] + nlink_mol;
+                    for (int j = link->index[atomIndex]; j < link->index[atomIndex + 1]; j++)
+                    {
+                        link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
+                    }
+                    if (link->index[atomIndex + 1] - link->index[atomIndex] > 0
+                        && atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock
+                                   < gmx::ssize(atomInfoOfMoleculeBlock->atomInfo))
+                    {
+                        SET_CGINFO_BOND_INTER(
+                                atomInfoOfMoleculeBlock
+                                        ->atomInfo[atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock]);
+                        numLinkedAtoms++;
+                    }
+                }
+                indexOfFirstAtomInMolecule += molt.atoms.nr;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, numLinkedAtoms);
+    }
+
+    dd->comm->bondedLinks = link;
+}
diff --git a/src/gromacs/domdec/makebondedlinks.h b/src/gromacs/domdec/makebondedlinks.h
new file mode 100644 (file)
index 0000000..01bc7a8
--- /dev/null
@@ -0,0 +1,65 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \libinternal \file
+ *
+ * \brief Declares a function that makes the list of links between
+ * atoms connected by bonded interactions.
+ *
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_MAKEBONDEDLINKS_H
+#define GMX_DOMDEC_MAKEBONDEDLINKS_H
+
+struct AtomInfoWithinMoleculeBlock;
+struct gmx_domdec_t;
+struct gmx_mtop_t;
+
+namespace gmx
+{
+template<typename>
+class ArrayRef;
+}
+
+/*! \brief Generate a list of links between atoms that are linked by bonded interactions
+ *
+ * Also stores whether atoms are linked in \p atomInfoForEachMoleculeBlock.
+ */
+void makeBondedLinks(gmx_domdec_t*                              dd,
+                     const gmx_mtop_t&                          mtop,
+                     gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock);
+
+#endif
index 022642cbfc8b2316905f1cc4bc5b0145abfa2b10..20c0c3e136ce036e67b3a38f57125d3bb83a1791 100644 (file)
@@ -62,6 +62,7 @@
 #include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/domdec/gpuhaloexchange.h"
 #include "gromacs/domdec/localatomsetmanager.h"
+#include "gromacs/domdec/makebondedlinks.h"
 #include "gromacs/domdec/partition.h"
 #include "gromacs/domdec/reversetopology.h"
 #include "gromacs/ewald/ewald_utils.h"