src/gromacs/topology/exclusionblocks.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "exclusionblocks.h"
  40
  41 #include <algorithm>
  42
  43 #include "gromacs/topology/block.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "gromacs/utility/stringutil.h"
  47
  48 namespace gmx
  49 {
  50
  51 void blockaToExclusionBlocks(const t_blocka* b, gmx::ArrayRef<ExclusionBlock> b2)
  52 {
  53     for (int i = 0; (i < b->nr); i++)
  54     {
  55         for (int j = b->index[i]; (j < b->index[i + 1]); j++)
  56         {
  57             b2[i].atomNumber.push_back(b->a[j]);
  58         }
  59     }
  60 }
  61
  62 void exclusionBlocksToBlocka(gmx::ArrayRef<const ExclusionBlock> b2, t_blocka* b)
  63 {
  64     int nra = 0;
  65     int i   = 0;
  66     for (const auto& block : b2)
  67     {
  68         b->index[i] = nra;
  69         int j       = 0;
  70         for (const auto& entry : block.atomNumber)
  71         {
  72             b->a[nra + j] = entry;
  73             j++;
  74         }
  75         nra += block.nra();
  76         i++;
  77     }
  78     /* terminate list */
  79     b->index[i] = nra;
  80 }
  81
  82 void mergeExclusions(t_blocka* excl, gmx::ArrayRef<ExclusionBlock> b2)
  83 {
  84     if (b2.empty())
  85     {
  86         return;
  87     }
  88     GMX_RELEASE_ASSERT(b2.ssize() == excl->nr,
  89                        "Cannot merge exclusions for "
  90                        "blocks that do not describe the same number "
  91                        "of particles");
  92
  93     /* Convert the t_blocka entries to ExclusionBlock form */
  94     blockaToExclusionBlocks(excl, b2);
  95
  96     /* Count and sort the exclusions */
  97     int nra = 0;
  98     for (auto& block : b2)
  99     {
 100         if (block.nra() > 0)
 101         {
 102             /* remove double entries */
 103             std::sort(block.atomNumber.begin(), block.atomNumber.end());
 104             for (auto atom = block.atomNumber.begin() + 1; atom != block.atomNumber.end();)
 105             {
 106                 GMX_RELEASE_ASSERT(atom < block.atomNumber.end(),
 107                                    "Need to stay in range of the size of the blocks");
 108                 auto prev = atom - 1;
 109                 if (*prev == *atom)
 110                 {
 111                     atom = block.atomNumber.erase(atom);
 112                 }
 113                 else
 114                 {
 115                     ++atom;
 116                 }
 117             }
 118             nra += block.nra();
 119         }
 120     }
 121     excl->nra = nra;
 122     srenew(excl->a, excl->nra);
 123
 124     exclusionBlocksToBlocka(b2, excl);
 125 }
 126
 127 } // namespace gmx