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 Topology helpers
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>
47 #include "gromacs/topology/exclusionblocks.h"
48 #include "gromacs/utility/smalloc.h"
49 #include "nblib/exception.h"
50 #include "nblib/listed_forces/transformations.h"
51 #include "nblib/topologyhelpers.h"
52 #include "nblib/util/internal.h"
60 std::vector<gmx::ExclusionBlock> toGmxExclusionBlock(const std::vector<std::tuple<int, int>>& tupleList)
62 std::vector<gmx::ExclusionBlock> ret;
64 auto firstLowerThan = [](auto const& tup1, auto const& tup2) {
65 return std::get<0>(tup1) < std::get<0>(tup2);
68 // Note this is a programming error as all particles should exclude at least themselves and empty topologies are not allowed.
69 // Note also that this is also checked in the parent function with a more informative error message.
70 assert((!tupleList.empty() && "No exclusions found.\n"));
72 // initialize pair of iterators delimiting the range of exclusions for
73 // the first particle in the list
74 auto range = std::equal_range(std::begin(tupleList), std::end(tupleList), tupleList[0], firstLowerThan);
75 auto it1 = range.first;
76 auto it2 = range.second;
78 // loop over all exclusions in molecule, linear in tupleList.size()
79 while (it1 != std::end(tupleList))
81 gmx::ExclusionBlock localBlock;
82 // loop over all exclusions for current particle
83 for (; it1 != it2; ++it1)
85 localBlock.atomNumber.push_back(std::get<1>(*it1));
88 ret.push_back(localBlock);
90 // update the upper bound of the range for the next particle
91 if (it1 != end(tupleList))
93 it2 = std::upper_bound(it1, std::end(tupleList), *it1, firstLowerThan);
100 std::vector<gmx::ExclusionBlock> offsetGmxBlock(std::vector<gmx::ExclusionBlock> inBlock, int offset)
102 // shift particle numbers by offset
103 for (auto& localBlock : inBlock)
105 std::transform(std::begin(localBlock.atomNumber), std::end(localBlock.atomNumber),
106 std::begin(localBlock.atomNumber), [offset](auto i) { return i + offset; });
112 int ParticleSequencer::operator()(const MoleculeName& moleculeName,
114 const ResidueName& residueName,
115 const ParticleName& particleName) const
119 return data_.at(moleculeName).at(moleculeNr).at(residueName).at(particleName);
121 catch (const std::out_of_range& outOfRange)
123 // TODO: use string format function once we have it
124 if (moleculeName.value() == residueName.value())
126 printf("No particle %s in residue %s in molecule %s found\n", particleName.value().c_str(),
127 residueName.value().c_str(), moleculeName.value().c_str());
131 printf("No particle %s in molecule %s found\n", particleName.value().c_str(),
132 moleculeName.value().c_str());
135 throw InputException(outOfRange.what());
139 void ParticleSequencer::build(const std::vector<std::tuple<Molecule, int>>& moleculesList)
142 for (auto& molNumberTuple : moleculesList)
144 const Molecule& molecule = std::get<0>(molNumberTuple);
145 const size_t numMols = std::get<1>(molNumberTuple);
147 auto& moleculeMap = data_[molecule.name()];
149 for (size_t i = 0; i < numMols; ++i)
151 auto& moleculeNrMap = moleculeMap[i];
152 for (int j = 0; j < molecule.numParticlesInMolecule(); ++j)
154 moleculeNrMap[molecule.residueName(j)][molecule.particleName(j)] = currentID++;
161 std::tuple<std::vector<size_t>, std::vector<I>>
162 collectInteractions(const std::vector<std::tuple<Molecule, int>>& molecules)
164 std::vector<I> collectedBonds;
165 std::vector<size_t> expansionArray;
166 for (auto& molNumberTuple : molecules)
168 const Molecule& molecule = std::get<0>(molNumberTuple);
169 size_t numMols = std::get<1>(molNumberTuple);
171 auto& interactions = pickType<I>(molecule.interactionData()).interactionTypes_;
173 std::vector<size_t> moleculeExpansion(interactions.size());
174 // assign indices to the bonds in the current molecule, continue counting from
175 // the number of bonds seen so far (=collectedBonds.size())
176 std::iota(begin(moleculeExpansion), end(moleculeExpansion), collectedBonds.size());
178 std::copy(begin(interactions), end(interactions), std::back_inserter(collectedBonds));
180 for (size_t i = 0; i < numMols; ++i)
182 std::copy(begin(moleculeExpansion), end(moleculeExpansion), std::back_inserter(expansionArray));
185 return std::make_tuple(expansionArray, collectedBonds);
188 /// \cond DO_NOT_DOCUMENT
189 #define COLLECT_BONDS_INSTANTIATE_TEMPLATE(x) \
190 template std::tuple<std::vector<size_t>, std::vector<x>> collectInteractions( \
191 const std::vector<std::tuple<Molecule, int>>&);
192 MAP(COLLECT_BONDS_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
195 namespace sequence_detail
198 //! Helper function to convert a tuple of strings into a particle index sequence
199 template<class Tuple, class F, class... Args, size_t... Is>
200 auto stringsToIndices_impl(const Tuple& tuple, std::index_sequence<Is...> is, F&& f, Args... args)
202 // return std::make_tuple(f(args..., std::get<2 * Is + 1>(tuple), std::get<2 * Is>(tuple))...);
204 return std::array<int, sizeof...(Is)>{ f(args..., std::get<2 * Is + 1>(tuple),
205 std::get<2 * Is>(tuple))... };
209 * This takes a tuple<(string, string) * nCenter> from molecule
210 * where nCenter = 2 for bonds, 3 for angles and 4 for dihedrals
211 * each (ResidueName, ParticleName)-pair is converted to a particle sequence index
212 * by calling the supplied function object f, containing the particleSequencer at the call site
213 * Therefore, the return type is tuple<int * nCenter>
216 template<class Tuple, class F, class... Args>
217 auto stringsToIndices(const Tuple& tuple, F&& f, Args... args)
219 auto is = std::make_index_sequence<std::tuple_size<Tuple>::value / 2>{};
220 return stringsToIndices_impl(tuple, is, std::forward<F>(f), args...);
223 //! Tuple ordering for two center interactions
224 [[maybe_unused]] static std::array<int, 2> nblibOrdering(const std::array<int, 2>& t)
226 // for bonds (two coordinate indices),
227 // we choose to store the lower sequence ID first. this allows for better unit tests
228 // that are agnostic to how the input was set up
229 int id1 = std::min(std::get<0>(t), std::get<1>(t));
230 int id2 = std::max(std::get<0>(t), std::get<1>(t));
232 return std::array<int, 2>{ id1, id2 };
235 //! Tuple ordering for three center interactions
236 [[maybe_unused]] static std::array<int, 3> nblibOrdering(const std::array<int, 3>& t)
238 // for angles (three coordinate indices),
239 // we choose to store the two non-center coordinate indices sorted.
240 // such that ret[0] < ret[2] always (ret = returned tuple)
241 int id1 = std::min(std::get<0>(t), std::get<2>(t));
242 int id3 = std::max(std::get<0>(t), std::get<2>(t));
244 return std::array<int, 3>{ id1, std::get<1>(t), id3 };
247 //! Tuple ordering for four center interactions
248 [[maybe_unused]] static std::array<int, 4> nblibOrdering(const std::array<int, 4>& t)
253 //! Tuple ordering for five center interactions
254 [[maybe_unused]] static std::array<int, 5> nblibOrdering(const std::array<int, 5>& t)
259 } // namespace sequence_detail
261 //! For each interaction, translate particle identifiers (moleculeName, nr, residueName,
262 //! particleName) to particle coordinate indices
264 std::vector<CoordinateIndex<B>> sequenceIDs(const std::vector<std::tuple<Molecule, int>>& molecules,
265 const detail::ParticleSequencer& particleSequencer)
267 std::vector<CoordinateIndex<B>> coordinateIndices;
269 auto callSequencer = [&particleSequencer](const MoleculeName& moleculeName, int i,
270 const ResidueName& residueName,
271 const ParticleName& particleName) {
272 return particleSequencer(moleculeName, i, residueName, particleName);
275 // loop over all molecules
276 for (const auto& molNumberTuple : molecules)
278 const Molecule& molecule = std::get<0>(molNumberTuple);
279 size_t numMols = std::get<1>(molNumberTuple);
281 for (size_t i = 0; i < numMols; ++i)
283 auto& interactions = pickType<B>(molecule.interactionData()).interactions_;
284 for (const auto& interactionString : interactions)
286 CoordinateIndex<B> index = sequence_detail::stringsToIndices(
287 interactionString, callSequencer, molecule.name(), i);
288 coordinateIndices.push_back(nblibOrdering(index));
292 return coordinateIndices;
295 /// \cond DO_NOT_DOCUMENT
296 #define SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE(x) \
297 template std::vector<CoordinateIndex<x>> sequenceIDs<x>( \
298 const std::vector<std::tuple<Molecule, int>>&, const detail::ParticleSequencer&);
299 MAP(SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
300 #undef SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE
304 std::tuple<std::vector<size_t>, std::vector<I>> eliminateDuplicateInteractions(const std::vector<I>& aggregatedInteractions)
306 std::vector<size_t> uniqueIndices(aggregatedInteractions.size());
307 std::vector<I> uniquInteractionsInstances;
308 // if there are no interactions of type B we're done now
309 if (aggregatedInteractions.empty())
311 return std::make_tuple(uniqueIndices, uniquInteractionsInstances);
314 // create 0,1,2,... sequence
315 std::iota(begin(uniqueIndices), end(uniqueIndices), 0);
317 std::vector<std::tuple<I, size_t>> enumeratedBonds(aggregatedInteractions.size());
318 // append each interaction with its index
319 std::transform(begin(aggregatedInteractions), end(aggregatedInteractions), begin(uniqueIndices),
320 begin(enumeratedBonds), [](I b, size_t i) { return std::make_tuple(b, i); });
322 auto sortKey = [](const auto& t1, const auto& t2) { return std::get<0>(t1) < std::get<0>(t2); };
323 // sort w.r.t bonds. the result will contain contiguous segments of identical bond instances
324 // the associated int indicates the original index of each BondType instance in the input vector
325 std::sort(begin(enumeratedBonds), end(enumeratedBonds), sortKey);
327 // initialize it1 and it2 to delimit first range of equal BondType instances
328 auto range = std::equal_range(begin(enumeratedBonds), end(enumeratedBonds), enumeratedBonds[0], sortKey);
329 auto it1 = range.first;
330 auto it2 = range.second;
332 // number of unique instances of BondType B = number of contiguous segments in enumeratedBonds =
333 // number of iterations in the outer while loop below
334 while (it1 != end(enumeratedBonds))
336 uniquInteractionsInstances.push_back(std::get<0>(*it1));
338 // loop over all identical BondType instances;
339 for (; it1 != it2; ++it1)
341 // we note down that the BondType instance at index <interactionIndex>
342 // can be found in the uniqueBondInstances container at index <uniqueBondInstances.size()>
343 int interactionIndex = std::get<1>(*it1);
344 uniqueIndices[interactionIndex] = uniquInteractionsInstances.size() - 1;
347 // Note it1 has been incremented and is now equal to it2
348 if (it1 != end(enumeratedBonds))
350 it2 = std::upper_bound(it1, end(enumeratedBonds), *it1, sortKey);
354 return make_tuple(uniqueIndices, uniquInteractionsInstances);
357 /// \cond DO_NOT_DOCUMENT
358 #define ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE(x) \
359 template std::tuple<std::vector<size_t>, std::vector<x>> eliminateDuplicateInteractions( \
360 const std::vector<x>& aggregatedBonds);
361 MAP(ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
362 #undef ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE
365 } // namespace detail