api/nblib/topologyhelpers.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, 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  * \brief
  37  * Implements nblib Topology helpers
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  * \author Artem Zhmurov <zhmurov@gmail.com>
  44  */
  45 #include <numeric>
  46
  47 #include "gromacs/topology/exclusionblocks.h"
  48 #include "gromacs/utility/smalloc.h"
  49 #include "nblib/exception.h"
  50 #include "nblib/topologyhelpers.h"
  51 #include "nblib/util/internal.h"
  52
  53 namespace nblib
  54 {
  55
  56 namespace detail
  57 {
  58
  59 std::vector<gmx::ExclusionBlock> toGmxExclusionBlock(const std::vector<std::tuple<int, int>>& tupleList)
  60 {
  61     std::vector<gmx::ExclusionBlock> ret;
  62
  63     auto firstLowerThan = [](auto const& tup1, auto const& tup2) {
  64         return std::get<0>(tup1) < std::get<0>(tup2);
  65     };
  66
  67     // Note this is a programming error as all particles should exclude at least themselves and empty topologies are not allowed.
  68     // Note also that this is also checked in the parent function with a more informative error message.
  69     assert((!tupleList.empty() && "No exclusions found.\n"));
  70
  71     // initialize pair of iterators delimiting the range of exclusions for
  72     // the first particle in the list
  73     auto range = std::equal_range(std::begin(tupleList), std::end(tupleList), tupleList[0], firstLowerThan);
  74     auto it1 = range.first;
  75     auto it2 = range.second;
  76
  77     // loop over all exclusions in molecule, linear in tupleList.size()
  78     while (it1 != std::end(tupleList))
  79     {
  80         gmx::ExclusionBlock localBlock;
  81         // loop over all exclusions for current particle
  82         for (; it1 != it2; ++it1)
  83         {
  84             localBlock.atomNumber.push_back(std::get<1>(*it1));
  85         }
  86
  87         ret.push_back(localBlock);
  88
  89         // update the upper bound of the range for the next particle
  90         if (it1 != end(tupleList))
  91         {
  92             it2 = std::upper_bound(it1, std::end(tupleList), *it1, firstLowerThan);
  93         }
  94     }
  95
  96     return ret;
  97 }
  98
  99 std::vector<gmx::ExclusionBlock> offsetGmxBlock(std::vector<gmx::ExclusionBlock> inBlock, int offset)
 100 {
 101     // shift particle numbers by offset
 102     for (auto& localBlock : inBlock)
 103     {
 104         std::transform(std::begin(localBlock.atomNumber), std::end(localBlock.atomNumber),
 105                        std::begin(localBlock.atomNumber), [offset](auto i) { return i + offset; });
 106     }
 107
 108     return inBlock;
 109 }
 110
 111 int ParticleSequencer::operator()(const MoleculeName& moleculeName,
 112                                   int                 moleculeNr,
 113                                   const ResidueName&  residueName,
 114                                   const ParticleName& particleName) const
 115 {
 116     try
 117     {
 118         return data_.at(moleculeName).at(moleculeNr).at(residueName).at(particleName);
 119     }
 120     catch (const std::out_of_range& outOfRange)
 121     {
 122         // TODO: use string format function once we have it
 123         if (moleculeName.value() == residueName.value())
 124         {
 125             printf("No particle %s in residue %s in molecule %s found\n", particleName.value().c_str(),
 126                    residueName.value().c_str(), moleculeName.value().c_str());
 127         }
 128         else
 129         {
 130             printf("No particle %s in molecule %s found\n", particleName.value().c_str(),
 131                    moleculeName.value().c_str());
 132         }
 133
 134         throw InputException(outOfRange.what());
 135     };
 136 }
 137
 138 void ParticleSequencer::build(const std::vector<std::tuple<Molecule, int>>& moleculesList)
 139 {
 140     int currentID = 0;
 141     for (auto& molNumberTuple : moleculesList)
 142     {
 143         const Molecule& molecule = std::get<0>(molNumberTuple);
 144         const size_t    numMols  = std::get<1>(molNumberTuple);
 145
 146         auto& moleculeMap = data_[molecule.name()];
 147
 148         for (size_t i = 0; i < numMols; ++i)
 149         {
 150             auto& moleculeNrMap = moleculeMap[i];
 151             for (int j = 0; j < molecule.numParticlesInMolecule(); ++j)
 152             {
 153                 moleculeNrMap[molecule.residueName(j)][molecule.particleName(j)] = currentID++;
 154             }
 155         }
 156     }
 157 }
 158
 159
 160 template<class I>
 161 std::tuple<std::vector<size_t>, std::vector<I>> eliminateDuplicateInteractions(const std::vector<I>& aggregatedInteractions)
 162 {
 163     std::vector<size_t> uniqueIndices(aggregatedInteractions.size());
 164     std::vector<I>      uniquInteractionsInstances;
 165     // if there are no interactions of type B we're done now
 166     if (aggregatedInteractions.empty())
 167     {
 168         return std::make_tuple(uniqueIndices, uniquInteractionsInstances);
 169     }
 170
 171     // create 0,1,2,... sequence
 172     std::iota(begin(uniqueIndices), end(uniqueIndices), 0);
 173
 174     std::vector<std::tuple<I, size_t>> enumeratedBonds(aggregatedInteractions.size());
 175     // append each interaction with its index
 176     std::transform(begin(aggregatedInteractions), end(aggregatedInteractions), begin(uniqueIndices),
 177                    begin(enumeratedBonds), [](I b, size_t i) { return std::make_tuple(b, i); });
 178
 179     auto sortKey = [](const auto& t1, const auto& t2) { return std::get<0>(t1) < std::get<0>(t2); };
 180     // sort w.r.t bonds. the result will contain contiguous segments of identical bond instances
 181     // the associated int indicates the original index of each BondType instance in the input vector
 182     std::sort(begin(enumeratedBonds), end(enumeratedBonds), sortKey);
 183
 184     // initialize it1 and it2 to delimit first range of equal BondType instances
 185     auto range = std::equal_range(begin(enumeratedBonds), end(enumeratedBonds), enumeratedBonds[0], sortKey);
 186     auto it1 = range.first;
 187     auto it2 = range.second;
 188
 189     // number of unique instances of BondType B = number of contiguous segments in enumeratedBonds =
 190     //         number of iterations in the outer while loop below
 191     while (it1 != end(enumeratedBonds))
 192     {
 193         uniquInteractionsInstances.push_back(std::get<0>(*it1));
 194
 195         // loop over all identical BondType instances;
 196         for (; it1 != it2; ++it1)
 197         {
 198             // we note down that the BondType instance at index <interactionIndex>
 199             // can be found in the uniqueBondInstances container at index <uniqueBondInstances.size()>
 200             int interactionIndex            = std::get<1>(*it1);
 201             uniqueIndices[interactionIndex] = uniquInteractionsInstances.size() - 1;
 202         }
 203
 204         // Note it1 has been incremented and is now equal to it2
 205         if (it1 != end(enumeratedBonds))
 206         {
 207             it2 = std::upper_bound(it1, end(enumeratedBonds), *it1, sortKey);
 208         }
 209     }
 210
 211     return make_tuple(uniqueIndices, uniquInteractionsInstances);
 212 }
 213
 214 } // namespace detail
 215
 216 } // namespace nblib