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> updateGroupsPerMoleculetype,
  57                                  real                                        temperature,
  58                                  int                                         numHomeAtoms) :
  59     globalToLocalMap_(numHomeAtoms),
  60     mtop_(mtop)
  61 {
  62     int firstUpdateGroupInMolecule = 0;
  63     for (const auto& molblock : mtop.molblock)
  64     {
  65         const auto& updateGroups = updateGroupsPerMoleculetype[molblock.type];
  66         indicesPerMoleculeblock_.push_back({ firstUpdateGroupInMolecule, updateGroups.numBlocks(), {} });
  67         auto& groupIndex = indicesPerMoleculeblock_.back().groupIndex_;
  68
  69         for (int block = 0; block < updateGroups.numBlocks(); block++)
  70         {
  71             groupIndex.insert(groupIndex.end(), updateGroups.block(block).size(), block);
  72         }
  73
  74         firstUpdateGroupInMolecule += molblock.nmol * updateGroups.numBlocks();
  75     }
  76
  77     maxUpdateGroupRadius_ = computeMaxUpdateGroupRadius(mtop, updateGroupsPerMoleculetype, temperature);
  78 }
  79
  80 void UpdateGroupsCog::addCogs(gmx::ArrayRef<const int>       globalAtomIndices,
  81                               gmx::ArrayRef<const gmx::RVec> coordinates)
  82 {
  83     const int    localAtomBegin = cogIndices_.size();
  84     const size_t cogBegin       = cogs_.size();
  85
  86     GMX_RELEASE_ASSERT(globalAtomIndices.ssize() >= localAtomBegin,
  87                        "addCogs should only be called to add COGs to the list that is already "
  88                        "present (which could be empty)");
  89
  90     cogIndices_.reserve(globalAtomIndices.size());
  91
  92     int moleculeBlock = 0;
  93     for (gmx::index localAtom = localAtomBegin; localAtom < globalAtomIndices.ssize(); localAtom++)
  94     {
  95         const int globalAtom = globalAtomIndices[localAtom];
  96         int       moleculeIndex;
  97         int       atomIndexInMolecule;
  98         mtopGetMolblockIndex(&mtop_, globalAtom, &moleculeBlock, &moleculeIndex, &atomIndexInMolecule);
  99         const auto& indicesForBlock        = indicesPerMoleculeblock_[moleculeBlock];
 100         int         globalUpdateGroupIndex = indicesForBlock.groupStart_
 101                                      + moleculeIndex * indicesForBlock.numGroupsPerMolecule_
 102                                      + indicesForBlock.groupIndex_[atomIndexInMolecule];
 103
 104         if (const int* cogIndexPtr = globalToLocalMap_.find(globalUpdateGroupIndex))
 105         {
 106             GMX_ASSERT(static_cast<size_t>(*cogIndexPtr) >= cogBegin,
 107                        "Added atoms should not be part of previously present groups");
 108
 109             cogIndices_.push_back(*cogIndexPtr);
 110
 111             cogs_[*cogIndexPtr] += coordinates[localAtom];
 112             numAtomsPerCog_[*cogIndexPtr]++;
 113         }
 114         else
 115         {
 116             const int cogIndex = cogs_.size();
 117
 118             globalToLocalMap_.insert(globalUpdateGroupIndex, cogIndex);
 119             cogIndices_.push_back(cogIndex);
 120             cogs_.push_back(coordinates[localAtom]);
 121             numAtomsPerCog_.push_back(1);
 122         }
 123     }
 124
 125     /* Divide sum of coordinates for each COG by the number of atoms */
 126     for (size_t i = cogBegin; i < cogs_.size(); i++)
 127     {
 128         const int numAtoms = numAtomsPerCog_[i];
 129         if (numAtoms > 1)
 130         {
 131             cogs_[i] /= numAtoms;
 132         }
 133     }
 134 }
 135
 136 void UpdateGroupsCog::clear()
 137 {
 138     cogIndices_.clear();
 139     cogs_.clear();
 140     numAtomsPerCog_.clear();
 141     globalToLocalMap_.clearAndResizeHashTable();
 142 }
 143
 144 } // namespace gmx