target_sources(nblib
PRIVATE
box.cpp
+ interactions.cpp
+ molecules.cpp
+ particletype.cpp
+ topologyhelpers.cpp
+ topology.cpp
)
+
set_target_properties(nblib
PROPERTIES
LINKER_LANGUAGE CXX
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib particle-types interactions
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include <cmath>
+#include <set>
+
+#include "nblib/exception.h"
+#include "nblib/interactions.h"
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+C6 NonBondedInteractionMap::getC6(const ParticleTypeName& first, const ParticleTypeName& second) const
+{
+ return std::get<0>(interactionMap_.at(std::make_tuple(first, second)));
+}
+
+C12 NonBondedInteractionMap::getC12(const ParticleTypeName& first, const ParticleTypeName& second) const
+{
+ return std::get<1>(interactionMap_.at(std::make_tuple(first, second)));
+}
+
+void NonBondedInteractionMap::setInteractions(const ParticleTypeName& first,
+ const ParticleTypeName& second,
+ const C6 c6_combo,
+ const C12 c12_combo)
+{
+ auto interactionKey = std::make_tuple(first, second);
+ interactionMap_[interactionKey] = std::make_tuple(c6_combo, c12_combo);
+}
+
+size_t NonBondedInteractionMap::count(const NonBondedInteractionMap::NamePairTuple& namePairTuple)
+{
+ return interactionMap_.count(namePairTuple);
+}
+
+namespace detail
+{
+
+//! Combines the non-bonded parameters from two particles for pairwise interactions
+real combineNonbondedParameters(real v, real w, CombinationRule combinationRule)
+{
+ if (combinationRule == CombinationRule::Geometric)
+ {
+ return std::sqrt(v * w);
+ }
+ else
+ {
+ throw InputException("unknown LJ Combination rule specified\n");
+ }
+}
+
+} // namespace detail
+
+ParticleTypesInteractions::ParticleTypesInteractions(CombinationRule cr) : combinationRule_(cr) {}
+
+ParticleTypesInteractions& ParticleTypesInteractions::add(const ParticleTypeName& particleTypeName,
+ C6 c6,
+ C12 c12)
+{
+ auto insertLocation = singleParticleInteractionsMap_.insert(
+ std::make_pair(particleTypeName, std::make_tuple(c6, c12)));
+
+ if (!insertLocation.second) // if particleTypeName already existed
+ {
+ if (std::get<0>(insertLocation.first->second) != c6
+ || std::get<1>(insertLocation.first->second) != c12)
+ {
+ std::string message = formatString(
+ "Attempting to add nonbonded interaction parameters for particle "
+ "type {} twice",
+ particleTypeName.value());
+ throw InputException(message);
+ }
+ }
+ return *this;
+}
+
+ParticleTypesInteractions& ParticleTypesInteractions::add(const ParticleTypeName& particleTypeName1,
+ const ParticleTypeName& particleTypeName2,
+ C6 c6,
+ C12 c12)
+{
+ auto interactionKey = std::make_tuple(particleTypeName1, particleTypeName2);
+ auto possibleInteractionKey = std::make_tuple(particleTypeName2, particleTypeName1);
+
+ auto insertLocation = twoParticlesInteractionsMap_.insert(
+ std::make_pair(interactionKey, std::make_tuple(c6, c12)));
+ twoParticlesInteractionsMap_.insert(std::make_pair(possibleInteractionKey, std::make_tuple(c6, c12)));
+
+ if (!insertLocation.second) // if particleTypeName already existed
+ {
+ if (std::get<0>(insertLocation.first->second) != c6
+ || std::get<1>(insertLocation.first->second) != c12)
+ {
+ std::string message = formatString(
+ "Attempting to add nonbonded interaction parameters between the particle types "
+ "{} {} twice",
+ particleTypeName1.value(), particleTypeName2.value());
+ throw InputException(message);
+ }
+ }
+ return *this;
+}
+
+NonBondedInteractionMap ParticleTypesInteractions::generateTable() const
+{
+ NonBondedInteractionMap nonbondedParameters_;
+
+ // creating the combination rule based interaction matrix
+ for (const auto& particleType1 : singleParticleInteractionsMap_)
+ {
+ C6 c6_1 = std::get<0>(particleType1.second);
+ C12 c12_1 = std::get<1>(particleType1.second);
+
+ for (const auto& particleType2 : singleParticleInteractionsMap_)
+ {
+ C6 c6_2 = std::get<0>(particleType2.second);
+ C12 c12_2 = std::get<1>(particleType2.second);
+
+ C6 c6_combo{ detail::combineNonbondedParameters(c6_1, c6_2, combinationRule_) };
+ C12 c12_combo{ detail::combineNonbondedParameters(c12_1, c12_2, combinationRule_) };
+
+ nonbondedParameters_.setInteractions(particleType1.first, particleType2.first, c6_combo,
+ c12_combo);
+ }
+ }
+
+ // updating the interaction matrix based on the user fine tuned parameters
+ for (const auto& particleTypeTuple : twoParticlesInteractionsMap_)
+ {
+ const auto& first = std::get<0>(particleTypeTuple.first);
+ const auto& second = std::get<1>(particleTypeTuple.first);
+
+ C6 c6_combo = std::get<0>(particleTypeTuple.second);
+ C12 c12_combo = std::get<1>(particleTypeTuple.second);
+
+ nonbondedParameters_.setInteractions(first, second, c6_combo, c12_combo);
+ }
+
+ std::set<ParticleTypeName> particleTypes;
+ for (auto const& typeKey : nonbondedParameters_)
+ { // we don't need to get<1> because the list is guaranteed to be symmetric
+ particleTypes.insert(std::get<0>(typeKey.first));
+ }
+
+ // check whether there is any missing interaction
+ for (const ParticleTypeName& particleTypeName1 : particleTypes)
+ {
+ for (const ParticleTypeName& particleTypeName2 : particleTypes)
+ {
+ auto interactionKey = std::make_tuple(particleTypeName1, particleTypeName2);
+ if (nonbondedParameters_.count(interactionKey) == 0)
+ {
+ std::string message = formatString("Missing interaction between {} {}",
+ particleTypeName1.value(), particleTypeName2.value());
+ throw InputException(message);
+ }
+ }
+ }
+ return nonbondedParameters_;
+}
+
+CombinationRule ParticleTypesInteractions::getCombinationRule() const
+{
+ return combinationRule_;
+}
+
+void ParticleTypesInteractions::merge(const ParticleTypesInteractions& other)
+{
+ for (const auto& keyval : other.singleParticleInteractionsMap_)
+ {
+ add(keyval.first, std::get<0>(keyval.second), std::get<1>(keyval.second));
+ }
+
+ for (const auto& keyval : other.twoParticlesInteractionsMap_)
+ {
+ add(std::get<0>(keyval.first), std::get<1>(keyval.first), std::get<0>(keyval.second),
+ std::get<1>(keyval.second));
+ }
+}
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib Molecule
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include <algorithm>
+#include <tuple>
+
+#include "nblib/exception.h"
+#include "nblib/molecules.h"
+#include "nblib/particletype.h"
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+
+Molecule::Molecule(MoleculeName moleculeName) : name_(std::move(moleculeName)) {}
+
+MoleculeName Molecule::name() const
+{
+ return name_;
+}
+
+Molecule& Molecule::addParticle(const ParticleName& particleName,
+ const ResidueName& residueName,
+ const Charge& charge,
+ ParticleType const& particleType)
+{
+ auto found = particleTypes_.find(particleType.name());
+ if (found == particleTypes_.end())
+ {
+ particleTypes_.insert(std::make_pair(particleType.name(), particleType));
+ }
+ else
+ {
+ if (!(found->second == particleType))
+ {
+ throw InputException(
+ "Differing ParticleTypes with identical names encountered in the same "
+ "molecule.");
+ }
+ }
+
+ particles_.emplace_back(ParticleData{ particleName, residueName, particleType.name(), charge });
+
+ // Add self exclusion. We just added the particle, so we know its index and that the exclusion doesn't exist yet
+ std::size_t id = particles_.size() - 1;
+ exclusions_.emplace_back(id, id);
+
+ return *this;
+}
+
+Molecule& Molecule::addParticle(const ParticleName& particleName,
+ const ResidueName& residueName,
+ ParticleType const& particleType)
+{
+ addParticle(particleName, residueName, Charge(0), particleType);
+
+ return *this;
+}
+
+Molecule& Molecule::addParticle(const ParticleName& particleName,
+ const Charge& charge,
+ ParticleType const& particleType)
+{
+ addParticle(particleName, ResidueName(name_), charge, particleType);
+
+ return *this;
+}
+
+Molecule& Molecule::addParticle(const ParticleName& particleName, const ParticleType& particleType)
+{
+ addParticle(particleName, ResidueName(name_), Charge(0), particleType);
+
+ return *this;
+}
+
+int Molecule::numParticlesInMolecule() const
+{
+ return particles_.size();
+}
+
+void Molecule::addExclusion(const int particleIndex, const int particleIndexToExclude)
+{
+ // We do not need to add exclusion in case the particle indexes are the same
+ // because self exclusion are added by addParticle
+ if (particleIndex != particleIndexToExclude)
+ {
+ exclusions_.emplace_back(particleIndex, particleIndexToExclude);
+ exclusions_.emplace_back(particleIndexToExclude, particleIndex);
+ }
+}
+
+void Molecule::addExclusion(std::tuple<std::string, std::string> particle,
+ std::tuple<std::string, std::string> particleToExclude)
+{
+ // duplication for the swapped pair happens in getExclusions()
+ exclusionsByName_.emplace_back(std::make_tuple(std::get<0>(particle), std::get<1>(particle),
+ std::get<0>(particleToExclude),
+ std::get<1>(particleToExclude)));
+}
+
+void Molecule::addExclusion(const std::string& particleName, const std::string& particleNameToExclude)
+{
+ addExclusion(std::make_tuple(particleName, name_), std::make_tuple(particleNameToExclude, name_));
+}
+
+const ParticleType& Molecule::at(const std::string& particleTypeName) const
+{
+ return particleTypes_.at(particleTypeName);
+}
+
+ParticleName Molecule::particleName(int i) const
+{
+ return ParticleName(particles_[i].particleName_);
+}
+
+ResidueName Molecule::residueName(int i) const
+{
+ return ResidueName(particles_[i].residueName_);
+}
+
+std::vector<std::tuple<int, int>> Molecule::getExclusions() const
+{
+ // tuples of (particleName, residueName, index)
+ std::vector<std::tuple<std::string, std::string, int>> indexKey;
+ indexKey.reserve(numParticlesInMolecule());
+
+ for (int i = 0; i < numParticlesInMolecule(); ++i)
+ {
+ indexKey.emplace_back(particles_[i].particleName_, particles_[i].residueName_, i);
+ }
+
+ std::sort(std::begin(indexKey), std::end(indexKey));
+
+ std::vector<std::tuple<int, int>> ret = exclusions_;
+ ret.reserve(exclusions_.size() + exclusionsByName_.size());
+
+ // normal operator<, except ignore third element
+ auto sortKey = [](const auto& tup1, const auto& tup2) {
+ if (std::get<0>(tup1) < std::get<0>(tup2))
+ {
+ return true;
+ }
+ else
+ {
+ return std::get<1>(tup1) < std::get<1>(tup2);
+ }
+ };
+
+ // convert exclusions given by names to indices and append
+ for (auto& tup : exclusionsByName_)
+ {
+ const std::string& particleName1 = std::get<0>(tup);
+ const std::string& residueName1 = std::get<1>(tup);
+ const std::string& particleName2 = std::get<2>(tup);
+ const std::string& residueName2 = std::get<3>(tup);
+
+ // look up first index (binary search)
+ auto it1 = std::lower_bound(std::begin(indexKey), std::end(indexKey),
+ std::make_tuple(particleName1, residueName2, 0), sortKey);
+
+ // make sure we have the (particleName,residueName) combo
+ if (it1 == std::end(indexKey) or std::get<0>(*it1) != particleName1 or std::get<1>(*it1) != residueName1)
+ {
+ throw std::runtime_error(
+ (std::string("Particle ") += particleName1 + std::string(" in residue ") +=
+ residueName1 + std::string(" not found in list of particles\n")));
+ }
+
+ int firstIndex = std::get<2>(*it1);
+
+ // look up second index (binary search)
+ auto it2 = std::lower_bound(std::begin(indexKey), std::end(indexKey),
+ std::make_tuple(particleName2, residueName2, 0), sortKey);
+
+ // make sure we have the (particleName,residueName) combo
+ if (it2 == std::end(indexKey) or std::get<0>(*it2) != particleName2 or std::get<1>(*it2) != residueName2)
+ {
+ throw std::runtime_error(
+ (std::string("Particle ") += particleName2 + std::string(" in residue ") +=
+ residueName2 + std::string(" not found in list of particles\n")));
+ }
+
+ int secondIndex = std::get<2>(*it2);
+
+ ret.emplace_back(firstIndex, secondIndex);
+ ret.emplace_back(secondIndex, firstIndex);
+ }
+
+ std::sort(std::begin(ret), std::end(ret));
+
+ auto uniqueEnd = std::unique(std::begin(ret), std::end(ret));
+ if (uniqueEnd != std::end(ret))
+ {
+ printf("[nblib] Warning: exclusionList for molecule %s contained duplicates",
+ name_.value().c_str());
+ }
+
+ ret.erase(uniqueEnd, std::end(ret));
+ return ret;
+}
+
+std::unordered_map<std::string, ParticleType> Molecule::particleTypesMap() const
+{
+ return particleTypes_;
+}
+
+std::vector<ParticleData> Molecule::particleData() const
+{
+ return particles_;
+}
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib ParticleType
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include "nblib/particletype.h"
+
+namespace nblib
+{
+
+ParticleType::ParticleType(ParticleTypeName name, Mass mass) : name_(std::move(name)), mass_(mass)
+{
+}
+
+ParticleTypeName ParticleType::name() const
+{
+ return name_;
+}
+
+Mass ParticleType::mass() const
+{
+ return mass_;
+}
+
+bool operator==(const ParticleType& a, const ParticleType& b)
+{
+ return a.name() == b.name() && a.mass() == b.mass();
+}
+
+} // namespace nblib
# in multiple test executables across the repository.
gmx_add_unit_test_library(nblib_test_infrastructure
testhelpers.cpp
- )
+ testsystems.cpp
+ )
target_include_directories(nblib_test_infrastructure PRIVATE ${PROJECT_SOURCE_DIR}/api)
target_include_directories(nblib_test_infrastructure SYSTEM PRIVATE ${PROJECT_SOURCE_DIR}/src/external)
-
set(testname "NbLibSetupTests")
set(exename "nblib-setup-test")
${exename}
CPP_SOURCE_FILES
# files with code for tests
- box.cpp
+ box.cpp
+ interactions.cpp
+ particletype.cpp
+ molecules.cpp
+ topology.cpp
)
target_link_libraries(${exename} PRIVATE nblib_test_infrastructure nblib)
target_include_directories(${exename} PRIVATE ${PROJECT_SOURCE_DIR}/api)
Box test = Box(length);
for (int i = 0; i < dimSize; ++i)
+ {
for (int j = 0; j < dimSize; ++j)
+ {
EXPECT_REAL_EQ_TOL(ref[i][j], test.legacyMatrix()[i][j], defaultRealTolerance());
+ }
+ }
}
} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements molecule setup tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include "nblib/interactions.h"
+#include "nblib/exception.h"
+#include "nblib/particletype.h"
+#include "nblib/tests/testsystems.h"
+
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+TEST(NBlibTest, NonBondedForceParamsCorrect)
+{
+ ParticleType atom1(ParticleTypeName("a1"), Mass(1));
+ ParticleType atom2(ParticleTypeName("a2"), Mass(1));
+ ParticleType atom3(ParticleTypeName("a3"), Mass(1));
+
+ ParticleTypesInteractions interactions;
+
+ auto c6_1 = C6(1.6);
+ auto c12_1 = C12(1.12);
+ C6 c6_2{ 2.6 };
+ C12 c12_2{ 2.12 };
+ C6 c6_3 = C6(3.6);
+ C12 c12_3 = C12(3.12);
+
+ auto c6comb = C6{ 40 };
+ auto c12comb = C12{ 41 };
+
+ interactions.add(atom1.name(), c6_1, c12_1);
+ interactions.add(atom2.name(), c6_2, c12_2);
+ interactions.add(atom3.name(), c6_3, c12_3);
+ interactions.add(atom2.name(), atom3.name(), c6comb, c12comb);
+
+ auto nbfp = interactions.generateTable();
+
+ //! self interaction for c6
+ EXPECT_REAL_EQ_TOL(c6_1, nbfp.getC6(atom1.name(), atom1.name()), gmx::test::defaultRealTolerance());
+ //! + symmetric pair
+ EXPECT_REAL_EQ_TOL(c12_1, nbfp.getC12(atom1.name(), atom1.name()), gmx::test::defaultRealTolerance());
+
+ //! geometric comb rule for c6
+ EXPECT_REAL_EQ_TOL(std::sqrt(c6_1 * c6_2), nbfp.getC6(atom1.name(), atom2.name()),
+ gmx::test::defaultRealTolerance());
+ //! + symmetric pair
+ EXPECT_REAL_EQ_TOL(std::sqrt(c6_1 * c6_2), nbfp.getC6(atom2.name(), atom1.name()),
+ gmx::test::defaultRealTolerance());
+
+ //! geometric comb rule for c12
+ EXPECT_REAL_EQ_TOL(std::sqrt(c12_1 * c12_2), nbfp.getC12(atom1.name(), atom2.name()),
+ gmx::test::defaultRealTolerance());
+
+ //! + symmetric par
+ EXPECT_REAL_EQ_TOL(std::sqrt(c12_1 * c12_2), nbfp.getC12(atom2.name(), atom1.name()),
+ gmx::test::defaultRealTolerance());
+
+ //! explicit pairwise interaction c6
+ EXPECT_REAL_EQ_TOL(c6comb, nbfp.getC6(atom2.name(), atom3.name()), gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(c6comb, nbfp.getC6(atom3.name(), atom2.name()), gmx::test::defaultRealTolerance());
+
+ //! explicit pairwise interaction c12
+ EXPECT_REAL_EQ_TOL(c12comb, nbfp.getC12(atom2.name(), atom3.name()),
+ gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(c12comb, nbfp.getC12(atom3.name(), atom2.name()),
+ gmx::test::defaultRealTolerance());
+
+ ParticleType atom4(ParticleTypeName("a4"), Mass(1));
+ interactions.add(atom3.name(), atom4.name(), C6(1), C12(2));
+ //! need to have self-interaction for all particles
+ EXPECT_THROW(auto interactionsTable = interactions.generateTable(), InputException);
+}
+
+TEST(NBlibTest, CanMergeInteractions)
+{
+ ParticleType atom1(ParticleTypeName("a1"), Mass(1));
+ ParticleType atom2(ParticleTypeName("a2"), Mass(1));
+ ParticleType atom3(ParticleTypeName("a3"), Mass(1));
+
+ ParticleTypesInteractions interactions;
+
+ auto c6_1 = C6(1.6);
+ auto c12_1 = C12(1.12);
+ C6 c6_2{ 2.6 };
+ C12 c12_2{ 2.12 };
+ C6 c6_3 = C6(3.6);
+ C12 c12_3 = C12(3.12);
+
+ auto c6comb = C6{ 40 };
+ auto c12comb = C12{ 41 };
+
+ interactions.add(atom1.name(), c6_1, c12_1);
+ interactions.add(atom2.name(), c6_2, c12_2);
+ interactions.add(atom3.name(), c6_3, c12_3);
+ interactions.add(atom2.name(), atom3.name(), c6comb, c12comb);
+
+ ParticleType atom4(ParticleTypeName("a4"), Mass(1));
+ ParticleType atom5(ParticleTypeName("a5"), Mass(1));
+ auto c6_4 = C6{ 4.6 };
+ auto c12_4 = C12{ 4.12 };
+
+ C6 c6_5{ 5.6 };
+ C12 c12_5{ 5.12 };
+
+ C6 c6_override = C6(45.6);
+ C12 c12_override = C12(45.12);
+
+ ParticleTypesInteractions otherInteractions;
+ otherInteractions.add(atom4.name(), c6_4, c12_4);
+ otherInteractions.add(atom5.name(), c6_5, c12_5);
+ otherInteractions.add(atom4.name(), atom5.name(), c6_override, c12_override);
+
+ interactions.merge(otherInteractions);
+
+ auto nbfp = interactions.generateTable();
+
+ EXPECT_REAL_EQ_TOL(std::sqrt(c6_3 * c6_4), nbfp.getC6(atom3.name(), atom4.name()),
+ gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(std::sqrt(c12_3 * c12_4), nbfp.getC12(atom3.name(), atom4.name()),
+ gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(c6_override, nbfp.getC6(atom4.name(), atom5.name()),
+ gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(c12_override, nbfp.getC12(atom4.name(), atom5.name()),
+ gmx::test::defaultRealTolerance());
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements molecule setup tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include "nblib/molecules.h"
+#include "nblib/exception.h"
+#include "nblib/particletype.h"
+#include "nblib/tests/testsystems.h"
+
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+TEST(NBlibTest, CanConstructMoleculeWithoutChargeOrResidueName)
+{
+ ArAtom arAtom;
+ ParticleType Ar(arAtom.particleTypeName, arAtom.mass);
+ Molecule argon(arAtom.moleculeName);
+ EXPECT_NO_THROW(argon.addParticle(arAtom.particleName, Ar));
+}
+
+TEST(NBlibTest, CanConstructMoleculeWithChargeWithoutResidueName)
+{
+ ArAtom arAtom;
+ ParticleType Ar(arAtom.particleTypeName, arAtom.mass);
+ Molecule argon(arAtom.moleculeName);
+ EXPECT_NO_THROW(argon.addParticle(arAtom.particleName, Charge(0), Ar));
+}
+
+TEST(NBlibTest, CanConstructMoleculeWithoutChargeWithResidueName)
+{
+ ArAtom arAtom;
+ ParticleType Ar(arAtom.particleTypeName, arAtom.mass);
+ Molecule argon(arAtom.moleculeName);
+ EXPECT_NO_THROW(argon.addParticle(arAtom.particleName, ResidueName("ar2"), Ar));
+}
+
+TEST(NBlibTest, CanConstructMoleculeWithChargeWithResidueName)
+{
+ ArAtom arAtom;
+ ParticleType Ar(arAtom.particleTypeName, arAtom.mass);
+ Molecule argon(arAtom.moleculeName);
+ EXPECT_NO_THROW(argon.addParticle(arAtom.particleName, ResidueName("ar2"), Charge(0), Ar));
+}
+
+TEST(NBlibTest, CanGetNumParticlesInMolecule)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMolecule();
+ auto numParticles = water.numParticlesInMolecule();
+
+ EXPECT_EQ(3, numParticles);
+}
+
+TEST(NBlibTest, CanConstructExclusionListFromNames)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMolecule();
+
+ std::vector<std::tuple<int, int>> exclusions = water.getExclusions();
+
+ std::vector<std::tuple<int, int>> reference{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
+ { 1, 2 }, { 2, 0 }, { 2, 1 }, { 2, 2 } };
+
+ ASSERT_EQ(exclusions.size(), 9);
+ for (std::size_t i = 0; i < exclusions.size(); ++i)
+ {
+ EXPECT_EQ(exclusions[i], reference[i]);
+ }
+}
+
+TEST(NBlibTest, CanConstructExclusionListFromIndices)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMoleculeWithoutExclusions();
+
+ //! Add the exclusions
+ water.addExclusion(1, 0);
+ water.addExclusion(2, 0);
+ water.addExclusion(1, 2);
+
+ std::vector<std::tuple<int, int>> exclusions = water.getExclusions();
+
+ std::vector<std::tuple<int, int>> reference{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
+ { 1, 2 }, { 2, 0 }, { 2, 1 }, { 2, 2 } };
+
+ ASSERT_EQ(exclusions.size(), 9);
+ for (std::size_t i = 0; i < exclusions.size(); ++i)
+ {
+ EXPECT_EQ(exclusions[i], reference[i]);
+ }
+}
+
+TEST(NBlibTest, CanConstructExclusionListFromNamesAndIndicesMixed)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMoleculeWithoutExclusions();
+
+ //! Add the exclusions
+ water.addExclusion("H1", "Oxygen");
+ water.addExclusion("H2", "Oxygen");
+ water.addExclusion(1, 2);
+
+ std::vector<std::tuple<int, int>> exclusions = water.getExclusions();
+
+ std::vector<std::tuple<int, int>> reference{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
+ { 1, 2 }, { 2, 0 }, { 2, 1 }, { 2, 2 } };
+
+ ASSERT_EQ(exclusions.size(), 9);
+ for (std::size_t i = 0; i < exclusions.size(); ++i)
+ {
+ EXPECT_EQ(exclusions[i], reference[i]);
+ }
+}
+
+TEST(NBlibTest, AtWorks)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMolecule();
+ EXPECT_NO_THROW(water.at("Ow"));
+ EXPECT_NO_THROW(water.at("H"));
+}
+
+TEST(NBlibTest, AtThrows)
+{
+ WaterMoleculeBuilder waterMolecule;
+ Molecule water = waterMolecule.waterMolecule();
+ EXPECT_THROW_GMX(water.at("Hw"), std::out_of_range);
+}
+
+TEST(NBlibTest, MoleculeThrowsSameParticleTypeNameDifferentMass)
+{
+ //! User error: Two different ParticleTypes with the same name
+ ParticleType atom1(ParticleTypeName("Atom"), Mass(1));
+ ParticleType atom2(ParticleTypeName("Atom"), Mass(2));
+
+ Molecule molecule(MoleculeName("UraniumDimer"));
+ EXPECT_NO_THROW(molecule.addParticle(ParticleName("U1"), atom1));
+ EXPECT_THROW(molecule.addParticle(ParticleName("U2"), atom2), InputException);
+}
+
+TEST(NBlibTest, MoleculeDontThrowsSameParticleTypeNameDifferentMass)
+{
+ //! User error: Two different ParticleTypes with the same name
+ ParticleType atom1(ParticleTypeName("Atom"), Mass(1));
+ ParticleType atom2(ParticleTypeName("Atom"), Mass(1));
+
+ Molecule molecule(MoleculeName("UraniumDimer"));
+ EXPECT_NO_THROW(molecule.addParticle(ParticleName("U1"), atom1));
+ EXPECT_NO_THROW(molecule.addParticle(ParticleName("U2"), atom2));
+}
+
+TEST(NBlibTest, MoleculeNoThrowsSameParticleTypeName)
+{
+ //! User error: Two different ParticleTypes with the same name
+ ParticleType atom1(ParticleTypeName("Atom"), Mass(1));
+ ParticleType atom2(ParticleTypeName("Atom"), Mass(1));
+
+ Molecule molecule(MoleculeName("UraniumDimer"));
+ EXPECT_NO_THROW(molecule.addParticle(ParticleName("U1"), atom1));
+ EXPECT_NO_THROW(molecule.addParticle(ParticleName("U2"), atom2));
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib AtomType tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ *
+ */
+#include <cmath>
+
+#include "nblib/particletype.h"
+#include "nblib/tests/testsystems.h"
+
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+
+TEST(NBlibTest, ParticleTypeNameCanBeConstructed)
+{
+ ArAtom arAtom;
+ ParticleType argonAtom(arAtom.particleTypeName, arAtom.mass);
+ EXPECT_EQ(ParticleTypeName(argonAtom.name()), arAtom.particleTypeName);
+}
+
+TEST(NBlibTest, ParticleTypeMassCanBeConstructed)
+{
+ ArAtom arAtom;
+ ParticleType argonAtom(arAtom.particleTypeName, arAtom.mass);
+ EXPECT_EQ(argonAtom.mass(), arAtom.mass);
+}
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements nblib test systems
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include "nblib/tests/testsystems.h"
+#include "nblib/exception.h"
+
+namespace nblib
+{
+
+//! User class to hold ParticleTypes
+//! Note: this is not part of NBLIB, users should write their own
+class ParticleLibrary
+{
+public:
+ ParticleLibrary()
+ {
+ ParticleType Ow(ParticleTypeName("Ow"), Mass(15.99940));
+ ParticleType H(ParticleTypeName("H"), Mass(1.008));
+ ParticleType OMet(ParticleTypeName("OMet"), Mass(15.999));
+ ParticleType CMet(ParticleTypeName("CMet"), Mass(15.035));
+ ParticleType Ar(ParticleTypeName("Ar"), Mass(39.94800));
+
+ particles_.insert(std::make_pair(Ow.name(), Ow));
+ particles_.insert(std::make_pair(H.name(), H));
+ particles_.insert(std::make_pair(OMet.name(), OMet));
+ particles_.insert(std::make_pair(CMet.name(), CMet));
+ particles_.insert(std::make_pair(Ar.name(), Ar));
+
+ c6_[Ow.name()] = 0.0026173456;
+ c6_[H.name()] = 0;
+ c6_[OMet.name()] = 0.0022619536;
+ c6_[CMet.name()] = 0.0088755241;
+ c6_[Ar.name()] = 0.0062647225;
+
+ c12_[Ow.name()] = 2.634129e-06;
+ c12_[H.name()] = 0;
+ c12_[OMet.name()] = 1.505529e-06;
+ c12_[CMet.name()] = 2.0852922e-05;
+ c12_[Ar.name()] = 9.847044e-06;
+ }
+
+ //! Get particle type using the string identifier
+ [[nodiscard]] ParticleType type(const std::string& particleTypeName) const
+ {
+ return particles_.at(ParticleTypeName(particleTypeName));
+ }
+
+ //! Get C6 parameter of a given particle
+ [[nodiscard]] C6 c6(const ParticleName& particleName) const
+ {
+ return c6_.at(ParticleTypeName(particleName.value()));
+ }
+
+ //! Get C12 parameter of a given particle
+ [[nodiscard]] C12 c12(const ParticleName& particleName) const
+ {
+ return c12_.at(ParticleTypeName(particleName.value()));
+ }
+
+private:
+ std::map<ParticleTypeName, ParticleType> particles_;
+ std::map<ParticleTypeName, C6> c6_;
+ std::map<ParticleTypeName, C12> c12_;
+};
+
+std::unordered_map<std::string, Charge> Charges{ { "Ow", Charge(-0.82) },
+ { "Hw", Charge(+0.41) },
+ { "OMet", Charge(-0.574) },
+ { "CMet", Charge(+0.176) },
+ { "HMet", Charge(+0.398) } };
+
+WaterMoleculeBuilder::WaterMoleculeBuilder() : water_(MoleculeName("SOL"))
+{
+ ParticleLibrary plib;
+
+ //! Add the particles
+ water_.addParticle(ParticleName("Oxygen"), Charges.at("Ow"), plib.type("Ow"));
+ water_.addParticle(ParticleName("H1"), Charges.at("Hw"), plib.type("H"));
+ water_.addParticle(ParticleName("H2"), Charges.at("Hw"), plib.type("H"));
+}
+
+Molecule WaterMoleculeBuilder::waterMolecule()
+{
+ addExclusionsFromNames();
+ return water_;
+}
+
+Molecule WaterMoleculeBuilder::waterMoleculeWithoutExclusions()
+{
+ return water_;
+}
+
+void WaterMoleculeBuilder::addExclusionsFromNames()
+{
+ water_.addExclusion("H1", "Oxygen");
+ water_.addExclusion("H2", "Oxygen");
+ water_.addExclusion("H1", "H2");
+}
+
+MethanolMoleculeBuilder::MethanolMoleculeBuilder() : methanol_(MoleculeName("MeOH"))
+{
+ ParticleLibrary library;
+
+ //! Add the particles
+ methanol_.addParticle(ParticleName("Me1"), Charges.at("CMet"), library.type("CMet"));
+ methanol_.addParticle(ParticleName("O2"), Charges.at("OMet"), library.type("OMet"));
+ methanol_.addParticle(ParticleName("H3"), Charges.at("HMet"), library.type("H"));
+
+ // Add the exclusions
+ methanol_.addExclusion("Me1", "O2");
+ methanol_.addExclusion("Me1", "H3");
+ methanol_.addExclusion("H3", "O2");
+}
+
+Molecule MethanolMoleculeBuilder::methanolMolecule()
+{
+ return methanol_;
+}
+
+
+Topology WaterTopologyBuilder::buildTopology(int numMolecules)
+{
+ ParticleLibrary library;
+
+ ParticleTypesInteractions interactions;
+ std::vector<std::string> typeNames = { "Ow", "H" };
+ for (const auto& name : typeNames)
+ {
+ interactions.add(ParticleTypeName(name), library.c6(ParticleName(name)),
+ library.c12(ParticleName(name)));
+ }
+
+ // Add some molecules to the topology
+ TopologyBuilder topologyBuilder;
+ topologyBuilder.addMolecule(water(), numMolecules);
+
+ // Add non-bonded interaction information
+ topologyBuilder.addParticleTypesInteractions(interactions);
+
+ Topology topology = topologyBuilder.buildTopology();
+ return topology;
+}
+
+Molecule WaterTopologyBuilder::water()
+{
+ return waterMolecule_.waterMolecule();
+}
+
+Topology SpcMethanolTopologyBuilder::buildTopology(int numWater, int numMethanol)
+{
+ ParticleLibrary library;
+
+ ParticleTypesInteractions interactions;
+ std::vector<std::string> typeNames = { "Ow", "H", "OMet", "CMet" };
+ for (const auto& name : typeNames)
+ {
+ interactions.add(ParticleTypeName(name), library.c6(ParticleName(name)),
+ library.c12(ParticleName(name)));
+ }
+
+ // Add some molecules to the topology
+ TopologyBuilder topologyBuilder;
+ topologyBuilder.addMolecule(methanol(), numMethanol);
+ topologyBuilder.addMolecule(water(), numWater);
+
+ // Add non-bonded interaction information
+ topologyBuilder.addParticleTypesInteractions(interactions);
+
+ Topology topology = topologyBuilder.buildTopology();
+ return topology;
+}
+
+Molecule SpcMethanolTopologyBuilder::methanol()
+{
+ return methanolMolecule_.methanolMolecule();
+}
+
+Molecule SpcMethanolTopologyBuilder::water()
+{
+ return waterMolecule_.waterMolecule();
+}
+
+ArgonTopologyBuilder::ArgonTopologyBuilder(const int& numParticles)
+{
+ ParticleLibrary library;
+
+ ParticleTypesInteractions nbinteractions;
+ nbinteractions.add(ParticleTypeName("Ar"), library.c6(ParticleName("Ar")),
+ library.c12(ParticleName("Ar")));
+
+ Molecule argonMolecule(MoleculeName("AR"));
+ argonMolecule.addParticle(ParticleName("AR"), library.type("Ar"));
+
+ topologyBuilder_.addMolecule(argonMolecule, numParticles);
+ topologyBuilder_.addParticleTypesInteractions((nbinteractions));
+}
+
+Topology ArgonTopologyBuilder::argonTopology()
+{
+ return topologyBuilder_.buildTopology();
+}
+
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib test systems
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#ifndef NBLIB_TESTSYSTEMS_H
+#define NBLIB_TESTSYSTEMS_H
+
+#include <cmath>
+
+#include "nblib/box.h"
+#include "nblib/molecules.h"
+#include "nblib/particletype.h"
+#include "nblib/simulationstate.h"
+#include "nblib/topology.h"
+
+namespace nblib
+{
+
+//! \internal \brief Parameters from gromos43A1
+struct ArAtom
+{
+ //! Argon Particle Name
+ ParticleName particleName = ParticleName("Ar");
+ //! Argon particle type name
+ ParticleTypeName particleTypeName = ParticleTypeName("Ar");
+ //! Argon molecule name
+ MoleculeName moleculeName = MoleculeName("Ar");
+ //! Argon Particle Mass
+ Mass mass = Mass(39.94800);
+ //! Argon C6 parameter
+ C6 c6{ 0.0062647225 };
+ //! Argon C12 parameter
+ C12 c12{ 9.847044e-06 };
+};
+
+//! Lookup table for charges needed for building topologies
+extern std::unordered_map<std::string, Charge> Charges;
+
+//! \internal \brief Make an SPC water molecule with parameters from gromos43A1
+class WaterMoleculeBuilder
+{
+public:
+ // There is no default ctor for a Molecule so it must be initialized
+ WaterMoleculeBuilder();
+
+ //! Return the initialized water Molecule, with exclusions
+ Molecule waterMolecule();
+
+ //! Return the initialized water Molecule, without exclusions
+ Molecule waterMoleculeWithoutExclusions();
+
+private:
+ //! The molecule
+ Molecule water_;
+
+ //! Add the exclusions from particle names. Private to prevent multiple calls
+ void addExclusionsFromNames();
+};
+
+//! \internal \brief Make a methanol molecule with parameters from gromos43A1
+class MethanolMoleculeBuilder
+{
+public:
+ // There is no default ctor for a Molecule so it must be initialized
+ MethanolMoleculeBuilder();
+
+ //! Return the initialized water Molecule, with exclusions
+ Molecule methanolMolecule();
+
+private:
+ //! The molecule
+ Molecule methanol_;
+};
+
+//! \internal \brief Build topology of water molecules of a specified number
+class WaterTopologyBuilder
+{
+public:
+ //! Return a topology with specified SPC water molecules
+ Topology buildTopology(int numMolecules);
+
+ //! Return the actual water Molecule used in the topology
+ Molecule water();
+
+private:
+ WaterMoleculeBuilder waterMolecule_;
+};
+
+//! \internal \brief Build topology of methanol+water molecules from specified numbers
+class SpcMethanolTopologyBuilder
+{
+public:
+ //! Return a topology with specified methanol molecules
+ Topology buildTopology(int numWater, int numMethanol);
+
+ //! Return the actual methanol Molecule used in the topology
+ Molecule methanol();
+
+ //! Return the actual water Molecule used in the topology
+ Molecule water();
+
+private:
+ MethanolMoleculeBuilder methanolMolecule_;
+ WaterMoleculeBuilder waterMolecule_;
+};
+
+//! \internal \brief Build topology of argon molecules of a specified number
+class ArgonTopologyBuilder
+{
+public:
+ //! Build a topology with specified argon molecules
+ ArgonTopologyBuilder(const int& numParticles);
+
+ //! Get the topology with specified argon molecules
+ Topology argonTopology();
+
+private:
+ TopologyBuilder topologyBuilder_;
+};
+
+} // namespace nblib
+#endif // NBLIB_TESTSYSTEMS_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements topology setup tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "gromacs/topology/exclusionblocks.h"
+#include "nblib/exception.h"
+#include "nblib/particletype.h"
+#include "nblib/tests/testsystems.h"
+#include "nblib/topology.h"
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+using ::testing::Eq;
+using ::testing::Pointwise;
+
+//! Compares all element between two lists of lists
+//! Todo: unify this with the identical function in nbkernelsystem test make this a method
+//! of ListOfLists<>
+template<typename T>
+void compareLists(const gmx::ListOfLists<T>& list, const std::vector<std::vector<T>>& v)
+{
+ ASSERT_EQ(list.size(), v.size());
+ for (std::size_t i = 0; i < list.size(); i++)
+ {
+ ASSERT_EQ(list[i].size(), v[i].size());
+ EXPECT_THAT(list[i], Pointwise(Eq(), v[i]));
+ }
+}
+
+// This is defined in src/gromacs/mdtypes/forcerec.h but there is also a
+// legacy C6 macro defined there that conflicts with the nblib C6 type.
+// Todo: Once that C6 has been refactored into a regular function, this
+// file can just include forcerec.h
+//! Macro to marks particles to have Van der Waals interactions
+#define SET_CGINFO_HAS_VDW(cgi) (cgi) = ((cgi) | (1 << 23))
+
+TEST(NBlibTest, TopologyHasNumParticles)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+ const int test = watersTopology.numParticles();
+ const int ref = 6;
+ EXPECT_EQ(ref, test);
+}
+
+TEST(NBlibTest, TopologyHasCharges)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+ const std::vector<real>& test = watersTopology.getCharges();
+ const std::vector<real>& ref = { Charges.at("Ow"), Charges.at("Hw"), Charges.at("Hw"),
+ Charges.at("Ow"), Charges.at("Hw"), Charges.at("Hw") };
+ EXPECT_EQ(ref, test);
+}
+
+TEST(NBlibTest, TopologyHasMasses)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+
+ const Mass refOwMass = waters.water().at("Ow").mass();
+ const Mass refHwMass = waters.water().at("H").mass();
+ const std::vector<Mass> ref = { refOwMass, refHwMass, refHwMass, refOwMass, refHwMass, refHwMass };
+ const std::vector<Mass> test = expandQuantity(watersTopology, &ParticleType::mass);
+ EXPECT_EQ(ref, test);
+}
+
+TEST(NBlibTest, TopologyHasParticleTypes)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+ const std::vector<ParticleType>& test = watersTopology.getParticleTypes();
+ const ParticleType refOw = waters.water().at("Ow");
+ const ParticleType refHw = waters.water().at("H");
+ const std::vector<ParticleType>& ref = { refOw, refHw };
+ const std::vector<ParticleType>& ref2 = { refHw, refOw };
+ EXPECT_TRUE(ref == test || ref2 == test);
+}
+
+TEST(NBlibTest, TopologyHasParticleTypeIds)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+
+ const std::vector<int>& testIds = watersTopology.getParticleTypeIdOfAllParticles();
+ const std::vector<ParticleType>& testTypes = watersTopology.getParticleTypes();
+
+ std::vector<ParticleType> testTypesExpanded;
+ testTypesExpanded.reserve(testTypes.size());
+ for (int i : testIds)
+ {
+ testTypesExpanded.push_back(testTypes[i]);
+ }
+
+ const ParticleType refOw = waters.water().at("Ow");
+ const ParticleType refHw = waters.water().at("H");
+ const std::vector<ParticleType> ref = { refOw, refHw, refHw, refOw, refHw, refHw };
+
+ EXPECT_TRUE(ref == testTypesExpanded);
+}
+
+TEST(NBlibTest, TopologyThrowsIdenticalParticleType)
+{
+ //! User error: Two different ParticleTypes with the same name
+ ParticleType U235(ParticleTypeName("Uranium"), Mass(235));
+ ParticleType U238(ParticleTypeName("Uranium"), Mass(238));
+
+ Molecule ud235(MoleculeName("UraniumDimer235"));
+ ud235.addParticle(ParticleName("U1"), U235);
+ ud235.addParticle(ParticleName("U2"), U235);
+
+ Molecule ud238(MoleculeName("UraniumDimer238"));
+ ud238.addParticle(ParticleName("U1"), U238);
+ ud238.addParticle(ParticleName("U2"), U238);
+
+ TopologyBuilder topologyBuilder;
+ topologyBuilder.addMolecule(ud235, 1);
+ EXPECT_THROW(topologyBuilder.addMolecule(ud238, 1), InputException);
+}
+
+TEST(NBlibTest, TopologyHasExclusions)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+ const gmx::ListOfLists<int>& testExclusions = watersTopology.getGmxExclusions();
+
+ const std::vector<std::vector<int>>& refExclusions = { { 0, 1, 2 }, { 0, 1, 2 }, { 0, 1, 2 },
+ { 3, 4, 5 }, { 3, 4, 5 }, { 3, 4, 5 } };
+
+ compareLists(testExclusions, refExclusions);
+}
+
+TEST(NBlibTest, TopologyHasSequencing)
+{
+ WaterTopologyBuilder waters;
+ Topology watersTopology = waters.buildTopology(2);
+
+ EXPECT_EQ(0, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"),
+ ParticleName("Oxygen")));
+ EXPECT_EQ(1, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1")));
+ EXPECT_EQ(2, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2")));
+ EXPECT_EQ(3, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"),
+ ParticleName("Oxygen")));
+ EXPECT_EQ(4, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H1")));
+ EXPECT_EQ(5, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H2")));
+}
+
+TEST(NBlibTest, toGmxExclusionBlockWorks)
+{
+ std::vector<std::tuple<int, int>> testInput{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
+ { 1, 2 }, { 2, 0 }, { 2, 1 }, { 2, 2 } };
+
+ std::vector<gmx::ExclusionBlock> reference;
+
+ gmx::ExclusionBlock localBlock;
+ localBlock.atomNumber.push_back(0);
+ localBlock.atomNumber.push_back(1);
+ localBlock.atomNumber.push_back(2);
+
+ reference.push_back(localBlock);
+ reference.push_back(localBlock);
+ reference.push_back(localBlock);
+
+ std::vector<gmx::ExclusionBlock> probe = detail::toGmxExclusionBlock(testInput);
+
+ ASSERT_EQ(reference.size(), probe.size());
+ for (size_t i = 0; i < reference.size(); ++i)
+ {
+ ASSERT_EQ(reference[i].atomNumber.size(), probe[i].atomNumber.size());
+ for (size_t j = 0; j < reference[i].atomNumber.size(); ++j)
+ {
+ EXPECT_EQ(reference[i].atomNumber[j], probe[i].atomNumber[j]);
+ }
+ }
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib Topology and TopologyBuilder
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include <numeric>
+
+#include "gromacs/topology/exclusionblocks.h"
+#include "gromacs/utility/listoflists.h"
+#include "gromacs/utility/smalloc.h"
+#include "nblib/exception.h"
+#include "nblib/particletype.h"
+#include "nblib/topology.h"
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+TopologyBuilder::TopologyBuilder() : numParticles_(0) {}
+
+gmx::ListOfLists<int> TopologyBuilder::createExclusionsListOfLists() const
+{
+ const auto& moleculesList = molecules_;
+
+ std::vector<gmx::ExclusionBlock> exclusionBlockGlobal;
+ exclusionBlockGlobal.reserve(numParticles_);
+
+ size_t particleNumberOffset = 0;
+ for (const auto& molNumberTuple : moleculesList)
+ {
+ const Molecule& molecule = std::get<0>(molNumberTuple);
+ size_t numMols = std::get<1>(molNumberTuple);
+ const auto& exclusions = molecule.getExclusions();
+
+ assert((!exclusions.empty()
+ && std::string("No exclusions found in the " + molecule.name().value() + " molecule.")
+ .c_str()));
+
+ std::vector<gmx::ExclusionBlock> exclusionBlockPerMolecule =
+ detail::toGmxExclusionBlock(exclusions);
+
+ // duplicate the exclusionBlockPerMolecule for the number of Molecules of (numMols)
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ auto offsetExclusions =
+ detail::offsetGmxBlock(exclusionBlockPerMolecule, particleNumberOffset);
+
+ std::copy(std::begin(offsetExclusions), std::end(offsetExclusions),
+ std::back_inserter(exclusionBlockGlobal));
+
+ particleNumberOffset += molecule.numParticlesInMolecule();
+ }
+ }
+
+ gmx::ListOfLists<int> exclusionsListOfListsGlobal;
+ for (const auto& block : exclusionBlockGlobal)
+ {
+ exclusionsListOfListsGlobal.pushBack(block.atomNumber);
+ }
+
+ return exclusionsListOfListsGlobal;
+}
+
+template<typename T, class Extractor>
+std::vector<T> TopologyBuilder::extractParticleTypeQuantity(Extractor&& extractor)
+{
+ auto& moleculesList = molecules_;
+
+ // returned object
+ std::vector<T> ret;
+ ret.reserve(numParticles_);
+
+ for (auto& molNumberTuple : moleculesList)
+ {
+ Molecule& molecule = std::get<0>(molNumberTuple);
+ size_t numMols = std::get<1>(molNumberTuple);
+
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ for (auto& particleData : molecule.particleData())
+ {
+ auto particleTypesMap = molecule.particleTypesMap();
+ ret.push_back(extractor(particleData, particleTypesMap));
+ }
+ }
+ }
+
+ return ret;
+}
+
+Topology TopologyBuilder::buildTopology()
+{
+ topology_.numParticles_ = numParticles_;
+
+ topology_.exclusions_ = createExclusionsListOfLists();
+ topology_.charges_ = extractParticleTypeQuantity<real>([](const auto& data, auto& map) {
+ ignore_unused(map);
+ return data.charge_;
+ });
+
+ // map unique ParticleTypes to IDs
+ std::unordered_map<std::string, int> nameToId;
+ for (auto& name_particleType_tuple : particleTypes_)
+ {
+ topology_.particleTypes_.push_back(name_particleType_tuple.second);
+ nameToId[name_particleType_tuple.first] = nameToId.size();
+ }
+
+ topology_.particleTypeIdOfAllParticles_ =
+ extractParticleTypeQuantity<int>([&nameToId](const auto& data, auto& map) {
+ ignore_unused(map);
+ return nameToId[data.particleTypeName_];
+ });
+
+ detail::ParticleSequencer particleSequencer;
+ particleSequencer.build(molecules_);
+ topology_.particleSequencer_ = std::move(particleSequencer);
+
+ topology_.combinationRule_ = particleTypesInteractions_.getCombinationRule();
+ topology_.nonBondedInteractionMap_ = particleTypesInteractions_.generateTable();
+
+ // Check whether there is any missing term in the particleTypesInteractions compared to the
+ // list of particletypes
+ for (const auto& particleType1 : particleTypes_)
+ {
+ for (const auto& particleType2 : particleTypes_)
+ {
+ auto interactionKey = std::make_tuple(ParticleTypeName(particleType1.first),
+ ParticleTypeName(particleType2.first));
+ if (topology_.nonBondedInteractionMap_.count(interactionKey) == 0)
+ {
+ std::string message =
+ formatString("Missing nonbonded interaction parameters for pair {} {}",
+ particleType1.first, particleType2.first);
+ throw InputException(message);
+ }
+ }
+ }
+
+ return topology_;
+}
+
+TopologyBuilder& TopologyBuilder::addMolecule(const Molecule& molecule, const int nMolecules)
+{
+ /*
+ * 1. Push-back a tuple of molecule type and nMolecules
+ * 2. Append exclusion list into the data structure
+ */
+
+ molecules_.emplace_back(molecule, nMolecules);
+ numParticles_ += nMolecules * molecule.numParticlesInMolecule();
+
+ auto particleTypesInMolecule = molecule.particleTypesMap();
+
+ for (const auto& name_type_tuple : particleTypesInMolecule)
+ {
+ // If we already have the particleType, we need to make
+ // sure that the type's parameters are actually the same
+ // otherwise we would overwrite them
+ if (particleTypes_.count(name_type_tuple.first) > 0)
+ {
+ if (!(particleTypes_.at(name_type_tuple.first) == name_type_tuple.second))
+ {
+ throw InputException("Differing ParticleTypes with identical names encountered");
+ }
+ }
+ }
+
+ // Note: insert does nothing if the key already exists
+ particleTypes_.insert(particleTypesInMolecule.begin(), particleTypesInMolecule.end());
+
+ return *this;
+}
+
+void TopologyBuilder::addParticleTypesInteractions(const ParticleTypesInteractions& particleTypesInteractions)
+{
+ particleTypesInteractions_.merge(particleTypesInteractions);
+}
+
+int Topology::numParticles() const
+{
+ return numParticles_;
+}
+
+std::vector<real> Topology::getCharges() const
+{
+ return charges_;
+}
+
+std::vector<ParticleType> Topology::getParticleTypes() const
+{
+ return particleTypes_;
+}
+
+std::vector<int> Topology::getParticleTypeIdOfAllParticles() const
+{
+ return particleTypeIdOfAllParticles_;
+}
+
+int Topology::sequenceID(MoleculeName moleculeName, int moleculeNr, ResidueName residueName, ParticleName particleName) const
+{
+ return particleSequencer_(moleculeName, moleculeNr, residueName, particleName);
+}
+
+NonBondedInteractionMap Topology::getNonBondedInteractionMap() const
+{
+ return nonBondedInteractionMap_;
+}
+
+CombinationRule Topology::getCombinationRule() const
+{
+ return combinationRule_;
+}
+
+gmx::ListOfLists<int> Topology::getGmxExclusions() const
+{
+ return exclusions_;
+}
+
+} // namespace nblib
ParticleTypesInteractions particleTypesInteractions_;
};
+//! utility function to extract Particle quantities and expand them to the full
+//! array of length numParticles()
+template<class F>
+inline auto expandQuantity(const Topology& topology, F&& particleTypeExtractor)
+{
+ using ValueType = decltype((std::declval<ParticleType>().*std::declval<F>())());
+
+ std::vector<ValueType> ret;
+ ret.reserve(topology.numParticles());
+
+ const std::vector<ParticleType>& particleTypes = topology.getParticleTypes();
+
+ for (size_t id : topology.getParticleTypeIdOfAllParticles())
+ {
+ ret.push_back((particleTypes[id].*particleTypeExtractor)());
+ }
+
+ return ret;
+}
+
} // namespace nblib
#endif // NBLIB_TOPOLOGY_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib Topology helpers
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include <numeric>
+
+#include "gromacs/topology/exclusionblocks.h"
+#include "gromacs/utility/smalloc.h"
+#include "nblib/exception.h"
+#include "nblib/topologyhelpers.h"
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+
+std::vector<gmx::ExclusionBlock> toGmxExclusionBlock(const std::vector<std::tuple<int, int>>& tupleList)
+{
+ std::vector<gmx::ExclusionBlock> ret;
+
+ auto firstLowerThan = [](auto const& tup1, auto const& tup2) {
+ return std::get<0>(tup1) < std::get<0>(tup2);
+ };
+
+ // initialize pair of iterators delimiting the range of exclusions for
+ // the first particle in the list
+ assert((!tupleList.empty() && "tupleList must not be empty\n"));
+ auto range = std::equal_range(std::begin(tupleList), std::end(tupleList), tupleList[0], firstLowerThan);
+ auto it1 = range.first;
+ auto it2 = range.second;
+
+ // loop over all exclusions in molecule, linear in tupleList.size()
+ while (it1 != std::end(tupleList))
+ {
+ gmx::ExclusionBlock localBlock;
+ // loop over all exclusions for current particle
+ for (; it1 != it2; ++it1)
+ {
+ localBlock.atomNumber.push_back(std::get<1>(*it1));
+ }
+
+ ret.push_back(localBlock);
+
+ // update the upper bound of the range for the next particle
+ if (it1 != end(tupleList))
+ {
+ it2 = std::upper_bound(it1, std::end(tupleList), *it1, firstLowerThan);
+ }
+ }
+
+ return ret;
+}
+
+std::vector<gmx::ExclusionBlock> offsetGmxBlock(std::vector<gmx::ExclusionBlock> inBlock, int offset)
+{
+ // shift particle numbers by offset
+ for (auto& localBlock : inBlock)
+ {
+ std::transform(std::begin(localBlock.atomNumber), std::end(localBlock.atomNumber),
+ std::begin(localBlock.atomNumber), [offset](auto i) { return i + offset; });
+ }
+
+ return inBlock;
+}
+
+int ParticleSequencer::operator()(const MoleculeName& moleculeName,
+ int moleculeNr,
+ const ResidueName& residueName,
+ const ParticleName& particleName) const
+{
+ try
+ {
+ return data_.at(moleculeName).at(moleculeNr).at(residueName).at(particleName);
+ }
+ catch (const std::out_of_range& outOfRange)
+ {
+ // TODO: use string format function once we have it
+ if (moleculeName.value() == residueName.value())
+ {
+ printf("No particle %s in residue %s in molecule %s found\n", particleName.value().c_str(),
+ residueName.value().c_str(), moleculeName.value().c_str());
+ }
+ else
+ {
+ printf("No particle %s in molecule %s found\n", particleName.value().c_str(),
+ moleculeName.value().c_str());
+ }
+
+ throw InputException(outOfRange.what());
+ };
+}
+
+void ParticleSequencer::build(const std::vector<std::tuple<Molecule, int>>& moleculesList)
+{
+ int currentID = 0;
+ for (auto& molNumberTuple : moleculesList)
+ {
+ const Molecule& molecule = std::get<0>(molNumberTuple);
+ const size_t numMols = std::get<1>(molNumberTuple);
+
+ auto& moleculeMap = data_[molecule.name()];
+
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ auto& moleculeNrMap = moleculeMap[i];
+ for (int j = 0; j < molecule.numParticlesInMolecule(); ++j)
+ {
+ moleculeNrMap[molecule.residueName(j)][molecule.particleName(j)] = currentID++;
+ }
+ }
+ }
+}
+
+
+template<class I>
+std::tuple<std::vector<size_t>, std::vector<I>> eliminateDuplicateInteractions(const std::vector<I>& aggregatedInteractions)
+{
+ std::vector<size_t> uniqueIndices(aggregatedInteractions.size());
+ std::vector<I> uniquInteractionsInstances;
+ // if there are no interactions of type B we're done now
+ if (aggregatedInteractions.empty())
+ {
+ return std::make_tuple(uniqueIndices, uniquInteractionsInstances);
+ }
+
+ // create 0,1,2,... sequence
+ std::iota(begin(uniqueIndices), end(uniqueIndices), 0);
+
+ std::vector<std::tuple<I, size_t>> enumeratedBonds(aggregatedInteractions.size());
+ // append each interaction with its index
+ std::transform(begin(aggregatedInteractions), end(aggregatedInteractions), begin(uniqueIndices),
+ begin(enumeratedBonds), [](I b, size_t i) { return std::make_tuple(b, i); });
+
+ auto sortKey = [](const auto& t1, const auto& t2) { return std::get<0>(t1) < std::get<0>(t2); };
+ // sort w.r.t bonds. the result will contain contiguous segments of identical bond instances
+ // the associated int indicates the original index of each BondType instance in the input vector
+ std::sort(begin(enumeratedBonds), end(enumeratedBonds), sortKey);
+
+ // initialize it1 and it2 to delimit first range of equal BondType instances
+ auto range = std::equal_range(begin(enumeratedBonds), end(enumeratedBonds), enumeratedBonds[0], sortKey);
+ auto it1 = range.first;
+ auto it2 = range.second;
+
+ // number of unique instances of BondType B = number of contiguous segments in enumeratedBonds =
+ // number of iterations in the outer while loop below
+ while (it1 != end(enumeratedBonds))
+ {
+ uniquInteractionsInstances.push_back(std::get<0>(*it1));
+
+ // loop over all identical BondType instances;
+ for (; it1 != it2; ++it1)
+ {
+ // we note down that the BondType instance at index <interactionIndex>
+ // can be found in the uniqueBondInstances container at index <uniqueBondInstances.size()>
+ int interactionIndex = std::get<1>(*it1);
+ uniqueIndices[interactionIndex] = uniquInteractionsInstances.size() - 1;
+ }
+
+ // Note it1 has been incremented and is now equal to it2
+ if (it1 != end(enumeratedBonds))
+ {
+ it2 = std::upper_bound(it1, end(enumeratedBonds), *it1, sortKey);
+ }
+ }
+
+ return make_tuple(uniqueIndices, uniquInteractionsInstances);
+}
+
+} // namespace detail
+
+} // namespace nblib
# \author Sebastian Keller <keller@cscs.ch>
#
+
+target_sources(nblib
+ PRIVATE
+ user.cpp
+ )
+
if(GMX_INSTALL_NBLIB_API)
install(FILES
user.h
DESTINATION include/nblib/util)
endif()
+
+if(BUILD_TESTING)
+ add_subdirectory(tests)
+endif()
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \inpublicapi \file
+ * \brief
+ * Implements nblib utilities
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+
+#ifndef NBLIB_UTIL_INTERNAL_H
+#define NBLIB_UTIL_INTERNAL_H
+
+#include <cassert>
+
+#include <sstream>
+#include <string>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+#include "nblib/basicdefinitions.h"
+#include "nblib/vector.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+//! Format strings for use in error messages
+std::string next_token(std::string& s, const std::string& delimiter);
+
+// Like std::void_t but for values
+template<auto...>
+using void_value_t = void;
+
+template<class T, class = void>
+struct HasValueMember : std::false_type
+{
+};
+
+template<class T>
+struct HasValueMember<T, void_value_t<T::value>> : std::true_type
+{
+};
+
+template<class T, class = void>
+struct AccessTypeMemberIfPresent
+{
+ typedef T type;
+};
+
+template<class T>
+struct AccessTypeMemberIfPresent<T, typename std::void_t<typename T::type>>
+{
+ typedef typename T::type type;
+};
+
+template<class T>
+using AccessTypeMemberIfPresent_t = typename AccessTypeMemberIfPresent<T>::type;
+
+//! this trait evaluates to std::true_type if T is the same as Tuple[N]
+//! OR if T is the same as the type member of Tuple[N]
+template<int N, typename T, typename Tuple>
+struct MatchTypeOrTypeMember :
+ std::disjunction<std::is_same<T, std::tuple_element_t<N, Tuple>>,
+ std::is_same<T, AccessTypeMemberIfPresent_t<std::tuple_element_t<N, Tuple>>>>
+{
+};
+
+//! recursion to check the next field N+1
+template<int N, class T, class Tuple, template<int, class, class> class Comparison, bool Match = false>
+struct MatchField_ :
+ std::integral_constant<size_t, MatchField_<N + 1, T, Tuple, Comparison, Comparison<N + 1, T, Tuple>{}>{}>
+{
+};
+
+//! recursion stop when Comparison<N, T, Tuple>::value is true
+template<int N, class T, class Tuple, template<int, class, class> class Comparison>
+struct MatchField_<N, T, Tuple, Comparison, true> : std::integral_constant<size_t, N>
+{
+};
+
+} // namespace detail
+
+/*! \brief The value member of this struct evaluates to the integral constant N for which
+ * the value member of Comparison<N, T, Tuple> is true
+ * and generates a compiler error if there is no such N
+ */
+template<class T, class Tuple, template<int, class, class> class Comparison>
+struct MatchField : detail::MatchField_<0, T, Tuple, Comparison, Comparison<0, T, Tuple>{}>
+{
+};
+
+/*! \brief Function to return the index in Tuple whose type matches T
+ * - If there are more than one, the first occurrence will be returned
+ * - If there is no such type, a compiler error from accessing a tuple out of range is generated
+ * Note that the default comparison operation supplied here also matches if the type member of Tuple[N] matches T
+ */
+template<typename T, typename Tuple, template<int, class, class> class Comparison = detail::MatchTypeOrTypeMember>
+struct FindIndex : std::integral_constant<size_t, MatchField<T, Tuple, Comparison>{}>
+{
+};
+
+//! Function to return the element in Tuple whose type matches T
+//! Note: if there are more than one, the first occurrence will be returned
+template<typename T, typename Tuple>
+decltype(auto) pickType(Tuple& tup)
+{
+ return std::get<FindIndex<T, Tuple>{}>(tup);
+}
+
+//! Utility to call function with each element in tuple_
+template<class F, class... Ts>
+void for_each_tuple(F&& func, std::tuple<Ts...>& tuple_)
+{
+ std::apply(
+ [f = func](auto&... args) {
+ [[maybe_unused]] auto list = std::initializer_list<int>{ (f(args), 0)... };
+ },
+ tuple_);
+}
+
+//! Utility to call function with each element in tuple_ with const guarantee
+template<class F, class... Ts>
+void for_each_tuple(F&& func, const std::tuple<Ts...>& tuple_)
+{
+ std::apply(
+ [f = func](auto&... args) {
+ [[maybe_unused]] auto list = std::initializer_list<int>{ (f(args), 0)... };
+ },
+ tuple_);
+}
+
+//! Format strings for use in error messages
+template<class... Args>
+std::string formatString(std::string fmt, Args... args)
+{
+ std::ostringstream os;
+ std::string delimiter = "{}";
+
+ std::initializer_list<int> unused{ 0, (os << detail::next_token(fmt, delimiter) << args, 0)... };
+ static_cast<void>(unused); // unused is not actually used
+
+ os << detail::next_token(fmt, delimiter);
+
+ return os.str();
+}
+
+} // namespace nblib
+
+#endif // NBLIB_UTIL_INTERNAL_H
--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2020, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+#
+# \author Victor Holanda <victor.holanda@cscs.ch>
+# \author Joe Jordan <ejjordan@kth.se>
+# \author Prashanth Kanduri <kanduri@cscs.ch>
+# \author Sebastian Keller <keller@cscs.ch>
+#
+
+set(testname "NbLibUtilTests")
+set(exename "nblib-util-test")
+
+gmx_add_gtest_executable(
+ ${exename}
+ CPP_SOURCE_FILES
+ user.cpp
+)
+target_link_libraries(${exename} PRIVATE nblib_test_infrastructure nblib)
+gmx_register_gtest_test(${testname} ${exename})
+add_dependencies(check-nblib ${exename})
+
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Sequence Name="generated-velocities">
+ <Int Name="Length">10</Int>
+ <Vector>
+ <Real Name="X">2.6981033066448674</Real>
+ <Real Name="Y">0.5539709822574701</Real>
+ <Real Name="Z">-0.66999586711644377</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">1.2082614830533767</Real>
+ <Real Name="Y">-0.94750311926855424</Real>
+ <Real Name="Z">-0.3939451998729408</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-0.25639656387655385</Real>
+ <Real Name="Y">0.15094422478273684</Real>
+ <Real Name="Z">-1.9023008166840518</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-2.6653384151271498</Real>
+ <Real Name="Y">1.0284873953993694</Real>
+ <Real Name="Z">1.8633553870427515</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-0.51905872128588526</Real>
+ <Real Name="Y">-1.5808470201670317</Real>
+ <Real Name="Z">0.59660484580158024</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-1.5358923264099191</Real>
+ <Real Name="Y">-4.0045494036796079</Real>
+ <Real Name="Z">2.3295413128004117</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">2.0461371078863495</Real>
+ <Real Name="Y">-0.65718768377791814</Real>
+ <Real Name="Z">-0.84789603041054018</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0.52471564906472223</Real>
+ <Real Name="Y">2.0471783903038139</Real>
+ <Real Name="Z">1.0757782548241039</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-0.53067549253173119</Real>
+ <Real Name="Y">1.0085627678362028</Real>
+ <Real Name="Z">1.5091817455269219</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0.71045843460233404</Real>
+ <Real Name="Y">-1.4262273469861517</Real>
+ <Real Name="Z">2.2175722402808948</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib utility tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+#include <vector>
+
+#include "nblib/tests/testhelpers.h"
+#include "nblib/util/user.h"
+
+#include "testutils/testasserts.h"
+
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+TEST(NBlibTest, checkNumericValues)
+{
+ std::vector<Vec3> vec;
+ vec.emplace_back(1., 1., 1.);
+ vec.emplace_back(2., 2., 2.);
+
+ bool ret = checkNumericValues(vec);
+ EXPECT_EQ(ret, true);
+}
+
+TEST(NBlibTest, checkNumericValuesHasNan)
+{
+ std::vector<Vec3> vec;
+ vec.emplace_back(1., 1., 1.);
+ vec.emplace_back(2., 2., 2.);
+
+ vec.emplace_back(NAN, NAN, NAN);
+
+ bool ret = checkNumericValues(vec);
+ EXPECT_EQ(ret, false);
+}
+
+TEST(NBlibTest, checkNumericValuesHasInf)
+{
+ std::vector<Vec3> vec;
+ vec.emplace_back(1., 1., 1.);
+ vec.emplace_back(2., 2., 2.);
+
+ vec.emplace_back(INFINITY, INFINITY, INFINITY);
+
+ bool ret = checkNumericValues(vec);
+ EXPECT_EQ(ret, false);
+}
+
+
+TEST(NBlibTest, GeneratedVelocitiesAreCorrect)
+{
+ constexpr size_t N = 10;
+ std::vector<real> masses(N, 1.0);
+ std::vector<Vec3> velocities;
+ velocities = generateVelocity(300.0, 1, masses);
+
+ Vector3DTest velocitiesTest;
+ velocitiesTest.testVectors(velocities, "generated-velocities");
+}
+TEST(NBlibTest, generateVelocitySize)
+{
+ constexpr int N = 10;
+ std::vector<real> masses(N, 1.0);
+ auto out = generateVelocity(300.0, 1, masses);
+ EXPECT_EQ(out.size(), N);
+}
+
+TEST(NBlibTest, generateVelocityCheckNumbers)
+{
+ constexpr int N = 10;
+ std::vector<real> masses(N, 1.0);
+ auto out = generateVelocity(300.0, 1, masses);
+ bool ret = checkNumericValues(out);
+ EXPECT_EQ(ret, true);
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements nblib utility functions
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+
+#include "nblib/util/user.h"
+#include "gromacs/random/tabulatednormaldistribution.h"
+#include "gromacs/random/threefry.h"
+#include "gromacs/utility/fatalerror.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+
+std::string next_token(std::string& s, const std::string& delimiter)
+{
+ std::string token = s.substr(0, s.find(delimiter));
+
+ std::size_t next = s.find(delimiter);
+ if (next == std::string::npos)
+ s.clear();
+ else
+ s.erase(0, next + delimiter.length());
+
+ return token;
+}
+
+} // namespace detail
+
+//! Generates an array of particle velocities based on the Maxwell-Boltzmann distribution
+//! using temperature, masses and a random number generator
+static std::vector<Vec3> low_mspeed(real tempi, std::vector<real> const& masses, gmx::ThreeFry2x64<>* rng)
+{
+ int nrdf;
+ real boltz;
+ real ekin, temp;
+ gmx::TabulatedNormalDistribution<real> normalDist;
+
+ std::vector<Vec3> velocities(masses.size());
+
+ boltz = BOLTZ * tempi;
+ ekin = 0.0;
+ nrdf = 0;
+ for (size_t i = 0; i < masses.size(); i++)
+ {
+ real mass = masses[i];
+ if (mass > 0)
+ {
+ rng->restart(i, 0);
+ real sd = std::sqrt(boltz / mass);
+ for (int m = 0; (m < dimSize); m++)
+ {
+ velocities[i][m] = sd * normalDist(*rng);
+ ekin += 0.5 * mass * velocities[i][m] * velocities[i][m];
+ }
+ nrdf += dimSize;
+ }
+ }
+ temp = (2.0 * ekin) / (nrdf * BOLTZ);
+ if (temp > 0)
+ {
+ real scal = std::sqrt(tempi / temp);
+ for (auto& vel : velocities)
+ {
+ for (int m = 0; (m < dimSize); m++)
+ {
+ vel[m] *= scal;
+ }
+ }
+ }
+ fprintf(stderr, "Velocities were taken from a Maxwell distribution at %g K\n", tempi);
+ if (debug)
+ {
+ fprintf(debug,
+ "Velocities were taken from a Maxwell distribution\n"
+ "Initial generated temperature: %12.5e (scaled to: %12.5e)\n",
+ temp, tempi);
+ }
+
+ return velocities;
+}
+
+//! Generate velocities from a Maxwell Boltzmann distribution, masses should be the
+//! same as the ones specified for the Topology object
+std::vector<Vec3> generateVelocity(real tempi, unsigned int seed, std::vector<real> const& masses)
+{
+
+ if (seed == 0)
+ {
+ seed = static_cast<int>(gmx::makeRandomSeed());
+ fprintf(stderr, "Using random seed %u for generating velocities\n", seed);
+ }
+ gmx::ThreeFry2x64<> rng(seed, gmx::RandomDomain::MaxwellVelocities);
+
+ return low_mspeed(tempi, masses, &rng);
+}
+
+//! Check within the container of gmx::RVecs for a NaN or inf
+bool checkNumericValues(const std::vector<Vec3>& values)
+{
+ for (auto val : values)
+ {
+ for (int m = 0; (m < dimSize); m++)
+ {
+ if (std::isnan(val[m]) or std::isinf(val[m]))
+ {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+} // namespace nblib
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff