api/nblib/molecules.cpp

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