2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements nblib Molecule
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>
48 #include "nblib/exception.h"
49 #include "nblib/molecules.h"
50 #include "nblib/particletype.h"
51 #include "nblib/util/internal.h"
57 Molecule::Molecule(MoleculeName moleculeName) : name_(std::move(moleculeName)) {}
59 MoleculeName Molecule::name() const
64 Molecule& Molecule::addParticle(const ParticleName& particleName,
65 const ResidueName& residueName,
67 ParticleType const& particleType)
69 auto found = particleTypes_.find(particleType.name());
70 if (found == particleTypes_.end())
72 particleTypes_.insert(std::make_pair(particleType.name(), particleType));
76 if (!(found->second == particleType))
79 "Differing ParticleTypes with identical names encountered in the same "
84 particles_.emplace_back(ParticleData{ particleName, residueName, particleType.name(), charge });
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);
93 Molecule& Molecule::addParticle(const ParticleName& particleName,
94 const ResidueName& residueName,
95 ParticleType const& particleType)
97 addParticle(particleName, residueName, Charge(0), particleType);
102 Molecule& Molecule::addParticle(const ParticleName& particleName,
103 const Charge& charge,
104 ParticleType const& particleType)
106 addParticle(particleName, ResidueName(name_), charge, particleType);
111 Molecule& Molecule::addParticle(const ParticleName& particleName, const ParticleType& particleType)
113 addParticle(particleName, ResidueName(name_), Charge(0), particleType);
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)
126 if (particleNameI == particleNameJ and residueNameI == residueNameJ)
128 throw InputException(std::string("Cannot add interaction of particle ")
129 + particleNameI.value() + " with itself in molecule " + name_.value());
132 auto& interactionContainer = pickType<Interaction>(interactionData_);
133 interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ, residueNameJ);
134 interactionContainer.interactionTypes_.push_back(interaction);
137 //! \cond DO_NOT_DOCUMENT
138 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
139 template void Molecule::addInteraction( \
140 const ParticleName& particleNameI, const ResidueName& residueNameI, \
141 const ParticleName& particleNameJ, const ResidueName& residueNameJ, const x& interaction);
142 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
143 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
146 // add interactions with default residue name
147 template<class Interaction>
148 void Molecule::addInteraction(const ParticleName& particleNameI,
149 const ParticleName& particleNameJ,
150 const Interaction& interaction)
152 addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_), interaction);
155 //! \cond DO_NOT_DOCUMENT
156 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
157 template void Molecule::addInteraction(const ParticleName& particleNameI, \
158 const ParticleName& particleNameJ, const x& interaction);
159 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
160 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
162 //! 3-particle interactions such as angles
163 template<class Interaction>
164 void Molecule::addInteraction(const ParticleName& particleNameI,
165 const ResidueName& residueNameI,
166 const ParticleName& particleNameJ,
167 const ResidueName& residueNameJ,
168 const ParticleName& particleNameK,
169 const ResidueName& residueNameK,
170 const Interaction& interaction)
172 if (particleNameI == particleNameJ and particleNameJ == particleNameK)
174 throw InputException(std::string("Cannot add interaction of particle ")
175 + particleNameI.value() + " with itself in molecule " + name_.value());
178 auto& interactionContainer = pickType<Interaction>(interactionData_);
179 interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ,
180 residueNameJ, particleNameK, residueNameK);
181 interactionContainer.interactionTypes_.push_back(interaction);
184 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
185 template void Molecule::addInteraction( \
186 const ParticleName& particleNameI, const ResidueName& residueNameI, \
187 const ParticleName& particleNameJ, const ResidueName& residueNameJ, \
188 const ParticleName& particleNameK, const ResidueName& residueNameK, const x& interaction);
189 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
190 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
192 template<class Interaction>
193 void Molecule::addInteraction(const ParticleName& particleNameI,
194 const ParticleName& particleNameJ,
195 const ParticleName& particleNameK,
196 const Interaction& interaction)
198 addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_),
199 particleNameK, ResidueName(name_), interaction);
202 #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
203 template void Molecule::addInteraction(const ParticleName& particleNameI, \
204 const ParticleName& particleNameJ, \
205 const ParticleName& particleNameK, const x& interaction);
206 MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
207 #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
210 int Molecule::numParticlesInMolecule() const
212 return particles_.size();
215 void Molecule::addExclusion(const int particleIndex, const int particleIndexToExclude)
217 // We do not need to add exclusion in case the particle indexes are the same
218 // because self exclusion are added by addParticle
219 if (particleIndex != particleIndexToExclude)
221 exclusions_.emplace_back(particleIndex, particleIndexToExclude);
222 exclusions_.emplace_back(particleIndexToExclude, particleIndex);
226 void Molecule::addExclusion(std::tuple<ParticleName, ResidueName> particle,
227 std::tuple<ParticleName, ResidueName> particleToExclude)
229 // duplication for the swapped pair happens in getExclusions()
230 exclusionsByName_.emplace_back(std::make_tuple(std::get<0>(particle), std::get<1>(particle),
231 std::get<0>(particleToExclude),
232 std::get<1>(particleToExclude)));
235 void Molecule::addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude)
237 addExclusion(std::make_tuple(particleName, ResidueName(name_)),
238 std::make_tuple(particleNameToExclude, ResidueName(name_)));
241 const Molecule::InteractionTuple& Molecule::interactionData() const
243 return interactionData_;
246 const ParticleType& Molecule::at(const std::string& particleTypeName) const
248 return particleTypes_.at(particleTypeName);
251 ParticleName Molecule::particleName(int i) const
253 return ParticleName(particles_[i].particleName_);
256 ResidueName Molecule::residueName(int i) const
258 return ResidueName(particles_[i].residueName_);
261 std::vector<std::tuple<int, int>> Molecule::getExclusions() const
263 // tuples of (particleName, residueName, index)
264 std::vector<std::tuple<std::string, std::string, int>> indexKey;
265 indexKey.reserve(numParticlesInMolecule());
267 for (int i = 0; i < numParticlesInMolecule(); ++i)
269 indexKey.emplace_back(particles_[i].particleName_, particles_[i].residueName_, i);
272 std::sort(std::begin(indexKey), std::end(indexKey));
274 std::vector<std::tuple<int, int>> ret = exclusions_;
275 ret.reserve(exclusions_.size() + exclusionsByName_.size());
277 // normal operator<, except ignore third element
278 auto sortKey = [](const auto& tup1, const auto& tup2) {
279 if (std::get<0>(tup1) < std::get<0>(tup2))
285 return std::get<1>(tup1) < std::get<1>(tup2);
289 // convert exclusions given by names to indices and append
290 for (auto& tup : exclusionsByName_)
292 const std::string& particleName1 = std::get<0>(tup);
293 const std::string& residueName1 = std::get<1>(tup);
294 const std::string& particleName2 = std::get<2>(tup);
295 const std::string& residueName2 = std::get<3>(tup);
297 // look up first index (binary search)
298 auto it1 = std::lower_bound(std::begin(indexKey), std::end(indexKey),
299 std::make_tuple(particleName1, residueName2, 0), sortKey);
301 // make sure we have the (particleName,residueName) combo
302 if (it1 == std::end(indexKey) or std::get<0>(*it1) != particleName1 or std::get<1>(*it1) != residueName1)
304 throw std::runtime_error(
305 (std::string("Particle ") += particleName1 + std::string(" in residue ") +=
306 residueName1 + std::string(" not found in list of particles\n")));
309 int firstIndex = std::get<2>(*it1);
311 // look up second index (binary search)
312 auto it2 = std::lower_bound(std::begin(indexKey), std::end(indexKey),
313 std::make_tuple(particleName2, residueName2, 0), sortKey);
315 // make sure we have the (particleName,residueName) combo
316 if (it2 == std::end(indexKey) or std::get<0>(*it2) != particleName2 or std::get<1>(*it2) != residueName2)
318 throw std::runtime_error(
319 (std::string("Particle ") += particleName2 + std::string(" in residue ") +=
320 residueName2 + std::string(" not found in list of particles\n")));
323 int secondIndex = std::get<2>(*it2);
325 ret.emplace_back(firstIndex, secondIndex);
326 ret.emplace_back(secondIndex, firstIndex);
329 std::sort(std::begin(ret), std::end(ret));
331 auto uniqueEnd = std::unique(std::begin(ret), std::end(ret));
332 if (uniqueEnd != std::end(ret))
334 printf("[nblib] Warning: exclusionList for molecule %s contained duplicates",
335 name_.value().c_str());
338 ret.erase(uniqueEnd, std::end(ret));
342 std::unordered_map<std::string, ParticleType> Molecule::particleTypesMap() const
344 return particleTypes_;
347 std::vector<ParticleData> Molecule::particleData() const