src/gromacs/mdlib/updategroupscog.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019,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 /* \internal \file
  36  *
  37  * \brief Implements the UpdateGroupsCog class.
  38  *
  39  * \author Berk Hess <hess@kth.se>
  40  * \ingroup module_mdlib
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "updategroupscog.h"
  46
  47 #include "gromacs/mdlib/updategroups.h"
  48 #include "gromacs/topology/block.h"
  49 #include "gromacs/topology/mtop_lookup.h"
  50 #include "gromacs/topology/topology.h"
  51
  52 namespace gmx
  53 {
  54
  55 UpdateGroupsCog::UpdateGroupsCog(const gmx_mtop_t&                           mtop,
  56                                  gmx::ArrayRef<const gmx::RangePartitioning> updateGroupingsPerMoleculeType,
  57                                  real                                        temperature,
  58                                  int                                         numHomeAtoms) :
  59     globalToLocalMap_(numHomeAtoms), mtop_(mtop)
  60 {
  61     int firstUpdateGroupInMolecule = 0;
  62     for (const auto& molblock : mtop.molblock)
  63     {
  64         const auto& updateGrouping = updateGroupingsPerMoleculeType[molblock.type];
  65         indicesPerMoleculeblock_.push_back({ firstUpdateGroupInMolecule, updateGrouping.numBlocks(), {} });
  66         auto& groupIndex = indicesPerMoleculeblock_.back().groupIndex_;
  67
  68         for (int block = 0; block < updateGrouping.numBlocks(); block++)
  69         {
  70             groupIndex.insert(groupIndex.end(), updateGrouping.block(block).size(), block);
  71         }
  72
  73         firstUpdateGroupInMolecule += molblock.nmol * updateGrouping.numBlocks();
  74     }
  75
  76     maxUpdateGroupRadius_ = computeMaxUpdateGroupRadius(mtop, updateGroupingsPerMoleculeType, temperature);
  77 }
  78
  79 void UpdateGroupsCog::addCogs(gmx::ArrayRef<const int>       globalAtomIndices,
  80                               gmx::ArrayRef<const gmx::RVec> coordinates)
  81 {
  82     const int    localAtomBegin = cogIndices_.size();
  83     const size_t cogBegin       = cogs_.size();
  84
  85     GMX_RELEASE_ASSERT(globalAtomIndices.ssize() >= localAtomBegin,
  86                        "addCogs should only be called to add COGs to the list that is already "
  87                        "present (which could be empty)");
  88
  89     cogIndices_.reserve(globalAtomIndices.size());
  90
  91     int moleculeBlock = 0;
  92     for (gmx::index localAtom = localAtomBegin; localAtom < globalAtomIndices.ssize(); localAtom++)
  93     {
  94         const int globalAtom = globalAtomIndices[localAtom];
  95         int       moleculeIndex;
  96         int       atomIndexInMolecule;
  97         mtopGetMolblockIndex(mtop_, globalAtom, &moleculeBlock, &moleculeIndex, &atomIndexInMolecule);
  98         const auto& indicesForBlock        = indicesPerMoleculeblock_[moleculeBlock];
  99         int         globalUpdateGroupIndex = indicesForBlock.groupStart_
 100                                      + moleculeIndex * indicesForBlock.numGroupsPerMolecule_
 101                                      + indicesForBlock.groupIndex_[atomIndexInMolecule];
 102
 103         if (const int* cogIndexPtr = globalToLocalMap_.find(globalUpdateGroupIndex))
 104         {
 105             GMX_ASSERT(static_cast<size_t>(*cogIndexPtr) >= cogBegin,
 106                        "Added atoms should not be part of previously present groups");
 107
 108             cogIndices_.push_back(*cogIndexPtr);
 109
 110             cogs_[*cogIndexPtr] += coordinates[localAtom];
 111             numAtomsPerCog_[*cogIndexPtr]++;
 112         }
 113         else
 114         {
 115             const int cogIndex = cogs_.size();
 116
 117             globalToLocalMap_.insert(globalUpdateGroupIndex, cogIndex);
 118             cogIndices_.push_back(cogIndex);
 119             cogs_.push_back(coordinates[localAtom]);
 120             numAtomsPerCog_.push_back(1);
 121         }
 122     }
 123
 124     /* Divide sum of coordinates for each COG by the number of atoms */
 125     for (size_t i = cogBegin; i < cogs_.size(); i++)
 126     {
 127         const int numAtoms = numAtomsPerCog_[i];
 128         if (numAtoms > 1)
 129         {
 130             cogs_[i] /= numAtoms;
 131         }
 132     }
 133 }
 134
 135 void UpdateGroupsCog::clear()
 136 {
 137     cogIndices_.clear();
 138     cogs_.clear();
 139     numAtomsPerCog_.clear();
 140     globalToLocalMap_.clearAndResizeHashTable();
 141 }
 142
 143 } // namespace gmx