src/gromacs/domdec/makebondedlinks.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief Defines a function that makes the list of links between
  39  * atoms connected by bonded interactions.
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_domdec
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "gromacs/domdec/makebondedlinks.h"
  48
  49 #include "gromacs/domdec/domdec_internal.h"
  50 #include "gromacs/domdec/options.h"
  51 #include "gromacs/domdec/reversetopology.h"
  52 #include "gromacs/mdtypes/atominfo.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/smalloc.h"
  55 #include "gromacs/topology/block.h"
  56 #include "gromacs/topology/mtop_util.h"
  57
  58 using gmx::ArrayRef;
  59 using gmx::DDBondedChecking;
  60
  61 /*! \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 */
  62 static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
  63 {
  64     bool bFound = false;
  65     for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
  66     {
  67         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
  68         if (link->a[k] == cg_gl_j)
  69         {
  70             bFound = TRUE;
  71         }
  72     }
  73     if (!bFound)
  74     {
  75         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
  76                            "Inconsistent allocation of link");
  77         /* Add this charge group link */
  78         if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
  79         {
  80             link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
  81             srenew(link->a, link->nalloc_a);
  82         }
  83         link->a[link->index[cg_gl + 1]] = cg_gl_j;
  84         link->index[cg_gl + 1]++;
  85     }
  86 }
  87
  88 void makeBondedLinks(gmx_domdec_t*                                   dd,
  89                      const gmx_mtop_t&                               mtop,
  90                      gmx::ArrayRef<gmx::AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock)
  91 {
  92
  93     if (!dd->comm->systemInfo.filterBondedCommunication)
  94     {
  95         /* Only communicate atoms based on cut-off */
  96         dd->comm->bondedLinks = nullptr;
  97         return;
  98     }
  99
 100     t_blocka* link = nullptr;
 101
 102     /* For each atom make a list of other atoms in the system
 103      * that a linked to it via bonded interactions
 104      * which are also stored in reverse_top.
 105      */
 106
 107     reverse_ilist_t ril_intermol;
 108     if (mtop.bIntermolecularInteractions)
 109     {
 110         t_atoms atoms;
 111
 112         atoms.nr   = mtop.natoms;
 113         atoms.atom = nullptr;
 114
 115         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
 116                            "We should have an ilist when intermolecular interactions are on");
 117
 118         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
 119         make_reverse_ilist(
 120                 *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
 121     }
 122
 123     snew(link, 1);
 124     snew(link->index, mtop.natoms + 1);
 125     link->nalloc_a = 0;
 126     link->a        = nullptr;
 127
 128     link->index[0]                 = 0;
 129     int indexOfFirstAtomInMolecule = 0;
 130     int numLinkedAtoms             = 0;
 131     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
 132     {
 133         const gmx_molblock_t& molb = mtop.molblock[mb];
 134         if (molb.nmol == 0)
 135         {
 136             continue;
 137         }
 138         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 139         /* Make a reverse ilist in which the interactions are linked
 140          * to all atoms, not only the first atom as in gmx_reverse_top.
 141          * The constraints are discarded here.
 142          */
 143         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
 144         reverse_ilist_t   ril;
 145         make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
 146
 147         gmx::AtomInfoWithinMoleculeBlock* atomInfoOfMoleculeBlock = &atomInfoForEachMoleculeBlock[mb];
 148
 149         int mol = 0;
 150         for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
 151         {
 152             for (int a = 0; a < molt.atoms.nr; a++)
 153             {
 154                 int atomIndex              = indexOfFirstAtomInMolecule + a;
 155                 link->index[atomIndex + 1] = link->index[atomIndex];
 156                 int i                      = ril.index[a];
 157                 while (i < ril.index[a + 1])
 158                 {
 159                     int ftype = ril.il[i++];
 160                     int nral  = NRAL(ftype);
 161                     /* Skip the ifunc index */
 162                     i++;
 163                     for (int j = 0; j < nral; j++)
 164                     {
 165                         int aj = ril.il[i + j];
 166                         if (aj != a)
 167                         {
 168                             check_link(link, atomIndex, indexOfFirstAtomInMolecule + aj);
 169                         }
 170                     }
 171                     i += nral_rt(ftype);
 172                 }
 173
 174                 if (mtop.bIntermolecularInteractions)
 175                 {
 176                     int i = ril_intermol.index[atomIndex];
 177                     while (i < ril_intermol.index[atomIndex + 1])
 178                     {
 179                         int ftype = ril_intermol.il[i++];
 180                         int nral  = NRAL(ftype);
 181                         /* Skip the ifunc index */
 182                         i++;
 183                         for (int j = 0; j < nral; j++)
 184                         {
 185                             /* Here we assume we have no charge groups;
 186                              * this has been checked above.
 187                              */
 188                             int aj = ril_intermol.il[i + j];
 189                             check_link(link, atomIndex, aj);
 190                         }
 191                         i += nral_rt(ftype);
 192                     }
 193                 }
 194                 if (link->index[atomIndex + 1] - link->index[atomIndex] > 0)
 195                 {
 196                     atomInfoOfMoleculeBlock->atomInfo[a] |= gmx::sc_atomInfo_BondCommunication;
 197                     numLinkedAtoms++;
 198                 }
 199             }
 200
 201             indexOfFirstAtomInMolecule += molt.atoms.nr;
 202         }
 203         int nlink_mol = link->index[indexOfFirstAtomInMolecule]
 204                         - link->index[indexOfFirstAtomInMolecule - molt.atoms.nr];
 205
 206         if (debug)
 207         {
 208             fprintf(debug,
 209                     "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
 210                     *molt.name,
 211                     molt.atoms.nr,
 212                     nlink_mol);
 213         }
 214
 215         if (molb.nmol > mol)
 216         {
 217             /* Copy the data for the rest of the molecules in this block */
 218             link->nalloc_a += (molb.nmol - mol) * nlink_mol;
 219             srenew(link->a, link->nalloc_a);
 220             for (; mol < molb.nmol; mol++)
 221             {
 222                 for (int a = 0; a < molt.atoms.nr; a++)
 223                 {
 224                     int atomIndex              = indexOfFirstAtomInMolecule + a;
 225                     link->index[atomIndex + 1] = link->index[atomIndex + 1 - molt.atoms.nr] + nlink_mol;
 226                     for (int j = link->index[atomIndex]; j < link->index[atomIndex + 1]; j++)
 227                     {
 228                         link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
 229                     }
 230                     if (link->index[atomIndex + 1] - link->index[atomIndex] > 0
 231                         && atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock
 232                                    < gmx::ssize(atomInfoOfMoleculeBlock->atomInfo))
 233                     {
 234                         atomInfoOfMoleculeBlock->atomInfo[atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock] |=
 235                                 gmx::sc_atomInfo_BondCommunication;
 236                         numLinkedAtoms++;
 237                     }
 238                 }
 239                 indexOfFirstAtomInMolecule += molt.atoms.nr;
 240             }
 241         }
 242     }
 243
 244     if (debug)
 245     {
 246         fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, numLinkedAtoms);
 247     }
 248
 249     dd->comm->bondedLinks = link;
 250 }