api/nblib/molecules.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 Molecule
  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 <algorithm>
  46 #include <tuple>
  47
  48 #include "nblib/exception.h"
  49 #include "nblib/molecules.h"
  50 #include "nblib/particletype.h"
  51 #include "nblib/util/internal.h"
  52
  53 namespace nblib
  54 {
  55
  56
  57 Molecule::Molecule(MoleculeName moleculeName) : name_(std::move(moleculeName)) {}
  58
  59 MoleculeName Molecule::name() const
  60 {
  61     return name_;
  62 }
  63
  64 Molecule& Molecule::addParticle(const ParticleName& particleName,
  65                                 const ResidueName&  residueName,
  66                                 const Charge&       charge,
  67                                 ParticleType const& particleType)
  68 {
  69     auto found = particleTypes_.find(particleType.name());
  70     if (found == particleTypes_.end())
  71     {
  72         particleTypes_.insert(std::make_pair(particleType.name(), particleType));
  73     }
  74     else
  75     {
  76         if (!(found->second == particleType))
  77         {
  78             throw InputException(
  79                     "Differing ParticleTypes with identical names encountered in the same "
  80                     "molecule.");
  81         }
  82     }
  83
  84     particles_.emplace_back(ParticleData{ particleName, residueName, particleType.name(), charge });
  85
  86     // Add self exclusion. We just added the particle, so we know its index and that the exclusion doesn't exist yet
  87     std::size_t id = particles_.size() - 1;
  88     exclusions_.emplace_back(id, id);
  89
  90     return *this;
  91 }
  92
  93 Molecule& Molecule::addParticle(const ParticleName& particleName,
  94                                 const ResidueName&  residueName,
  95                                 ParticleType const& particleType)
  96 {
  97     addParticle(particleName, residueName, Charge(0), particleType);
  98
  99     return *this;
 100 }
 101
 102 Molecule& Molecule::addParticle(const ParticleName& particleName,
 103                                 const Charge&       charge,
 104                                 ParticleType const& particleType)
 105 {
 106     addParticle(particleName, ResidueName(name_), charge, particleType);
 107
 108     return *this;
 109 }
 110
 111 Molecule& Molecule::addParticle(const ParticleName& particleName, const ParticleType& particleType)
 112 {
 113     addParticle(particleName, ResidueName(name_), Charge(0), particleType);
 114
 115     return *this;
 116 }
 117
 118 //! Two-particle interactions such as bonds and LJ1-4
 119 template<class Interaction>
 120 void Molecule::addInteraction(const ParticleName& particleNameI,
 121                               const ResidueName&  residueNameI,
 122                               const ParticleName& particleNameJ,
 123                               const ResidueName&  residueNameJ,
 124                               const Interaction&  interaction)
 125 {
 126     if (particleNameI == particleNameJ and residueNameI == residueNameJ)
 127     {
 128         throw InputException(std::string("Cannot add interaction of particle ")
 129                              + particleNameI.value() + " with itself in molecule " + name_.value());
 130     }
 131
 132     auto& interactionContainer = pickType<Interaction>(interactionData_);
 133     interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ, residueNameJ);
 134     interactionContainer.interactionTypes_.push_back(interaction);
 135 }
 136
 137 //! \cond DO_NOT_DOCUMENT
 138 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
 139     template void Molecule::addInteraction(const ParticleName& particleNameI, \
 140                                            const ResidueName&  residueNameI,  \
 141                                            const ParticleName& particleNameJ, \
 142                                            const ResidueName&  residueNameJ,  \
 143                                            const x&            interaction);
 144 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
 145 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
 146 //! \endcond
 147
 148 // add interactions with default residue name
 149 template<class Interaction>
 150 void Molecule::addInteraction(const ParticleName& particleNameI,
 151                               const ParticleName& particleNameJ,
 152                               const Interaction&  interaction)
 153 {
 154     addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_), interaction);
 155 }
 156
 157 //! \cond DO_NOT_DOCUMENT
 158 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
 159     template void Molecule::addInteraction(     \
 160             const ParticleName& particleNameI, const ParticleName& particleNameJ, const x& interaction);
 161 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
 162 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
 163
 164 //! 3-particle interactions such as angles
 165 template<class Interaction>
 166 void Molecule::addInteraction(const ParticleName& particleNameI,
 167                               const ResidueName&  residueNameI,
 168                               const ParticleName& particleNameJ,
 169                               const ResidueName&  residueNameJ,
 170                               const ParticleName& particleNameK,
 171                               const ResidueName&  residueNameK,
 172                               const Interaction&  interaction)
 173 {
 174     if (particleNameI == particleNameJ and particleNameJ == particleNameK)
 175     {
 176         throw InputException(std::string("Cannot add interaction of particle ")
 177                              + particleNameI.value() + " with itself in molecule " + name_.value());
 178     }
 179
 180     auto& interactionContainer = pickType<Interaction>(interactionData_);
 181     interactionContainer.interactions_.emplace_back(
 182             particleNameI, residueNameI, particleNameJ, residueNameJ, particleNameK, residueNameK);
 183     interactionContainer.interactionTypes_.push_back(interaction);
 184 }
 185
 186 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
 187     template void Molecule::addInteraction(const ParticleName& particleNameI, \
 188                                            const ResidueName&  residueNameI,  \
 189                                            const ParticleName& particleNameJ, \
 190                                            const ResidueName&  residueNameJ,  \
 191                                            const ParticleName& particleNameK, \
 192                                            const ResidueName&  residueNameK,  \
 193                                            const x&            interaction);
 194 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
 195 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
 196
 197 template<class Interaction>
 198 void Molecule::addInteraction(const ParticleName& particleNameI,
 199                               const ParticleName& particleNameJ,
 200                               const ParticleName& particleNameK,
 201                               const Interaction&  interaction)
 202 {
 203     addInteraction(particleNameI,
 204                    ResidueName(name_),
 205                    particleNameJ,
 206                    ResidueName(name_),
 207                    particleNameK,
 208                    ResidueName(name_),
 209                    interaction);
 210 }
 211
 212 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
 213     template void Molecule::addInteraction(const ParticleName& particleNameI, \
 214                                            const ParticleName& particleNameJ, \
 215                                            const ParticleName& particleNameK, \
 216                                            const x&            interaction);
 217 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
 218 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
 219 //! \endcond
 220
 221 int Molecule::numParticlesInMolecule() const
 222 {
 223     return particles_.size();
 224 }
 225
 226 void Molecule::addExclusion(const int particleIndex, const int particleIndexToExclude)
 227 {
 228     // We do not need to add exclusion in case the particle indexes are the same
 229     // because self exclusion are added by addParticle
 230     if (particleIndex != particleIndexToExclude)
 231     {
 232         exclusions_.emplace_back(particleIndex, particleIndexToExclude);
 233         exclusions_.emplace_back(particleIndexToExclude, particleIndex);
 234     }
 235 }
 236
 237 void Molecule::addExclusion(std::tuple<ParticleName, ResidueName> particle,
 238                             std::tuple<ParticleName, ResidueName> particleToExclude)
 239 {
 240     // duplication for the swapped pair happens in getExclusions()
 241     exclusionsByName_.emplace_back(std::make_tuple(std::get<0>(particle),
 242                                                    std::get<1>(particle),
 243                                                    std::get<0>(particleToExclude),
 244                                                    std::get<1>(particleToExclude)));
 245 }
 246
 247 void Molecule::addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude)
 248 {
 249     addExclusion(std::make_tuple(particleName, ResidueName(name_)),
 250                  std::make_tuple(particleNameToExclude, ResidueName(name_)));
 251 }
 252
 253 const Molecule::InteractionTuple& Molecule::interactionData() const
 254 {
 255     return interactionData_;
 256 }
 257
 258 const ParticleType& Molecule::at(const std::string& particleTypeName) const
 259 {
 260     return particleTypes_.at(particleTypeName);
 261 }
 262
 263 ParticleName Molecule::particleName(int i) const
 264 {
 265     return ParticleName(particles_[i].particleName_);
 266 }
 267
 268 ResidueName Molecule::residueName(int i) const
 269 {
 270     return ResidueName(particles_[i].residueName_);
 271 }
 272
 273 std::vector<std::tuple<int, int>> Molecule::getExclusions() const
 274 {
 275     // tuples of (particleName, residueName, index)
 276     std::vector<std::tuple<std::string, std::string, int>> indexKey;
 277     indexKey.reserve(numParticlesInMolecule());
 278
 279     for (int i = 0; i < numParticlesInMolecule(); ++i)
 280     {
 281         indexKey.emplace_back(particles_[i].particleName_, particles_[i].residueName_, i);
 282     }
 283
 284     std::sort(std::begin(indexKey), std::end(indexKey));
 285
 286     std::vector<std::tuple<int, int>> ret = exclusions_;
 287     ret.reserve(exclusions_.size() + exclusionsByName_.size());
 288
 289     // normal operator<, except ignore third element
 290     auto sortKey = [](const auto& tup1, const auto& tup2) {
 291         if (std::get<0>(tup1) < std::get<0>(tup2))
 292         {
 293             return true;
 294         }
 295         else
 296         {
 297             return std::get<1>(tup1) < std::get<1>(tup2);
 298         }
 299     };
 300
 301     // convert exclusions given by names to indices and append
 302     for (auto& tup : exclusionsByName_)
 303     {
 304         const std::string& particleName1 = std::get<0>(tup);
 305         const std::string& residueName1  = std::get<1>(tup);
 306         const std::string& particleName2 = std::get<2>(tup);
 307         const std::string& residueName2  = std::get<3>(tup);
 308
 309         // look up first index (binary search)
 310         auto it1 = std::lower_bound(std::begin(indexKey),
 311                                     std::end(indexKey),
 312                                     std::make_tuple(particleName1, residueName2, 0),
 313                                     sortKey);
 314
 315         // make sure we have the (particleName,residueName) combo
 316         if (it1 == std::end(indexKey) or std::get<0>(*it1) != particleName1 or std::get<1>(*it1) != residueName1)
 317         {
 318             throw std::runtime_error(
 319                     (std::string("Particle ") += particleName1 + std::string(" in residue ") +=
 320                      residueName1 + std::string(" not found in list of particles\n")));
 321         }
 322
 323         int firstIndex = std::get<2>(*it1);
 324
 325         // look up second index (binary search)
 326         auto it2 = std::lower_bound(std::begin(indexKey),
 327                                     std::end(indexKey),
 328                                     std::make_tuple(particleName2, residueName2, 0),
 329                                     sortKey);
 330
 331         // make sure we have the (particleName,residueName) combo
 332         if (it2 == std::end(indexKey) or std::get<0>(*it2) != particleName2 or std::get<1>(*it2) != residueName2)
 333         {
 334             throw std::runtime_error(
 335                     (std::string("Particle ") += particleName2 + std::string(" in residue ") +=
 336                      residueName2 + std::string(" not found in list of particles\n")));
 337         }
 338
 339         int secondIndex = std::get<2>(*it2);
 340
 341         ret.emplace_back(firstIndex, secondIndex);
 342         ret.emplace_back(secondIndex, firstIndex);
 343     }
 344
 345     std::sort(std::begin(ret), std::end(ret));
 346
 347     auto uniqueEnd = std::unique(std::begin(ret), std::end(ret));
 348     if (uniqueEnd != std::end(ret))
 349     {
 350         printf("[nblib] Warning: exclusionList for molecule %s contained duplicates",
 351                name_.value().c_str());
 352     }
 353
 354     ret.erase(uniqueEnd, std::end(ret));
 355     return ret;
 356 }
 357
 358 std::unordered_map<std::string, ParticleType> Molecule::particleTypesMap() const
 359 {
 360     return particleTypes_;
 361 }
 362
 363 std::vector<ParticleData> Molecule::particleData() const
 364 {
 365     return particles_;
 366 }
 367
 368 } // namespace nblib