* \author Prashanth Kanduri <kanduri@cscs.ch>
* \author Sebastian Keller <keller@cscs.ch>
*/
+#include "nblib/exception.h"
#include "nblib/forcecalculator.h"
#include "nblib/gmxcalculator.h"
#include "nblib/gmxsetup.h"
void ForceCalculator::compute(gmx::ArrayRef<const Vec3> coordinates, gmx::ArrayRef<Vec3> forces)
{
- return gmxForceCalculator_->compute(coordinates, forces);
+ if (coordinates.size() != forces.size())
+ {
+ throw InputException("Coordinates array and force buffer size mismatch");
+ }
+
+ gmxForceCalculator_->compute(coordinates, forces);
}
void ForceCalculator::updatePairList(gmx::ArrayRef<const int> particleInfoAllVdW,
* costly to create this object since much of the SimulationState and NBKernelOptions has to be
* passed to the gromacs backend. However, once constructed, compute can be called repeatedly only
* paying the cost of the actual nonbonded force calculation. Repeated calls to compute on the same
- * coordinated will always return the same forces (within precision), so the user must update the
+ * coordinates will always return the same forces (within precision), so the user must update the
* positions using the forces generated here to advance a simulation. If the coordinates move
* sufficiently far from their positions at construction time, the efficiency of the calculation
* will suffer. To alleviate this, the user can call updatePairList.
// update the coordinates in the backend
nbv_->convertCoordinates(gmx::AtomLocality::Local, false, coordinateInput);
- // set forces to zero
- std::fill(forceOutput.begin(), forceOutput.end(), gmx::RVec{ 0, 0, 0 });
-
nbv_->dispatchNonbondedKernel(gmx::InteractionLocality::Local, *interactionConst_, *stepWork_,
enbvClearFYes, *forcerec_, enerd_.get(), nrnb_.get());
class SimulationState;
struct NBKernelOptions;
+/*! \brief GROMACS non-bonded force calculation backend
+ *
+ * This class encapsulates the various GROMACS data structures and their interplay
+ * from the NBLIB user. The class is a private member of the ForceCalculator and
+ * is not intended for the public interface.
+ *
+ * Handles the task of storing the simulation problem description using the internal
+ * representation used within GROMACS. It currently supports short range non-bonded
+ * interactions (PP) on a single node CPU.
+ *
+ */
+
class GmxForceCalculator final
{
public:
const Box& box);
private:
+ //! Friend to allow setting up private members in this class
friend class NbvSetupUtil;
//! Non-Bonded Verlet object for force calculation
// Put everything together
auto nbv = std::make_unique<nonbonded_verlet_t>(std::move(pairlistSets), std::move(pairSearch),
- std::move(atomData), kernelSetup, nullptr, nullptr);
+ std::move(atomData), kernelSetup, nullptr,
+ nullWallcycle);
// Needs to be called with the number of unique ParticleTypes
nbnxn_atomdata_init(gmx::MDLogger(), nbv->nbat.get(), kernelSetup.kernelType, combinationRule,
namespace nblib
{
+/*! \brief Sets up the GROMACS data structures for the non-bonded force calculator
+ *
+ * This data structure initializes the GmxForceCalculator object which internally
+ * contains various objects needed to perform non-bonded force calculations using
+ * the internal representation for the problem as required for GROMACS.
+ *
+ * The public functions of this class basically translate the problem description
+ * specified by the user in NBLIB. This ultimately returns the GmxForceCalculator
+ * object which is used by the ForceCalculator object in the user-facing library.
+ *
+ */
class NbvSetupUtil final
{
public:
//! Sets up t_forcerec object on the GmxForceCalculator
void setupForceRec(const matrix& box);
+ //! Returns a unique pointer a GmxForceCalculator object
std::unique_ptr<GmxForceCalculator> getGmxForceCalculator()
{
return std::move(gmxForceCalculator_);
std::unique_ptr<GmxForceCalculator> gmxForceCalculator_;
};
+/*! \brief Calls the setup utilities needed to initialize a GmxForceCalculator object
+ *
+ * The GmxSetupDirector encapsulates the multi-stage setup of the GmxForceCalculator which
+ * is done using the public functions of the NbvSetupUtil. This separation ensures that the
+ * NbvSetupUtil object is temporary in scope. The function definition makes it easy for the
+ * developers to follow the sequence of calls and the dataflow involved in setting up
+ * the non-bonded force calculation backend. This is the only function needed to be called
+ * from the ForceCalculator during construction.
+ *
+ */
class GmxSetupDirector
{
public:
namespace nblib
{
-// NOLINTNEXTLINE(performance-unnecessary-value-param)
+
LeapFrog::LeapFrog(const Topology& topology, const Box& box) : box_(box)
{
inverseMasses_.resize(topology.numParticles());
}
}
-void Molecule::addExclusion(std::tuple<std::string, std::string> particle,
- std::tuple<std::string, std::string> particleToExclude)
+void Molecule::addExclusion(std::tuple<ParticleName, ResidueName> particle,
+ std::tuple<ParticleName, ResidueName> 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<1>(particleToExclude)));
}
-void Molecule::addExclusion(const std::string& particleName, const std::string& particleNameToExclude)
+void Molecule::addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude)
{
- addExclusion(std::make_tuple(particleName, name_), std::make_tuple(particleNameToExclude, name_));
+ addExclusion(std::make_tuple(particleName, ResidueName(name_)),
+ std::make_tuple(particleNameToExclude, ResidueName(name_)));
}
const ParticleType& Molecule::at(const std::string& particleTypeName) const
namespace nblib
{
+class TopologyBuilder;
+
//! Named type for unique identifier for a particle in a molecule
using ParticleName = StrongType<std::string, struct ParticleNameParameter>;
void addExclusion(int particleIndex, int particleIndexToExclude);
//! Specify an exclusion with particle and residue names that have been added to molecule
- void addExclusion(std::tuple<std::string, std::string> particle,
- std::tuple<std::string, std::string> particleToExclude);
+ void addExclusion(std::tuple<ParticleName, ResidueName> particle,
+ std::tuple<ParticleName, ResidueName> particleToExclude);
//! Specify an exclusion with particle names that have been added to molecule
- void addExclusion(const std::string& particleName, const std::string& particleNameToExclude);
+ void addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude);
//! The number of molecules
int numParticlesInMolecule() const;
box_(box),
topology_(std::move(topology))
{
- if (!checkNumericValues(coordinates))
+ auto numParticles = topology_.numParticles();
+
+ if (int(coordinates.size()) != numParticles)
+ {
+ throw InputException("Coordinates array size mismatch");
+ }
+
+ if (int(velocities.size()) != numParticles)
+ {
+ throw InputException("Velocities array size mismatch");
+ }
+
+ if (int(forces.size()) != numParticles)
+ {
+ throw InputException("Force buffer array size mismatch");
+ }
+
+ if (!isRealValued(coordinates))
{
throw InputException("Input coordinates has at least one NaN");
}
coordinates_ = coordinates;
- if (!checkNumericValues(velocities))
+ if (!isRealValued(velocities))
{
throw InputException("Input velocities has at least one NaN");
}
Molecule water = waterMolecule.waterMoleculeWithoutExclusions();
//! Add the exclusions
- water.addExclusion("H1", "Oxygen");
- water.addExclusion("H2", "Oxygen");
+ water.addExclusion(ParticleName("H1"), ParticleName("Oxygen"));
+ water.addExclusion(ParticleName("H2"), ParticleName("Oxygen"));
water.addExclusion(1, 2);
std::vector<std::tuple<int, int>> exclusions = water.getExclusions();
gmx::ArrayRef<Vec3> forces(simState.forces());
forceCalculator.compute(simState.coordinates(), simState.forces());
+ // copy computed forces to another array
std::vector<Vec3> forces_1(forces.size());
std::copy(forces.begin(), forces.end(), begin(forces_1));
+ // zero original force buffer
+ zeroCartesianArray(forces);
+
// check if forces change without update step
forceCalculator.compute(simState.coordinates(), forces);
// update
integrator.integrate(1.0, simState.coordinates(), simState.velocities(), simState.forces());
+ // zero original force buffer
+ zeroCartesianArray(forces);
+
// step 2
- forceCalculator.compute(simState.coordinates(), simState.forces());
+ forceCalculator.compute(simState.coordinates(), forces);
std::vector<Vec3> forces_2(forces.size());
std::copy(forces.begin(), forces.end(), begin(forces_2));
void WaterMoleculeBuilder::addExclusionsFromNames()
{
- water_.addExclusion("H1", "Oxygen");
- water_.addExclusion("H2", "Oxygen");
- water_.addExclusion("H1", "H2");
+ water_.addExclusion(ParticleName("H1"), ParticleName("Oxygen"));
+ water_.addExclusion(ParticleName("H2"), ParticleName("Oxygen"));
+ water_.addExclusion(ParticleName("H1"), ParticleName("H2"));
}
MethanolMoleculeBuilder::MethanolMoleculeBuilder() : methanol_(MoleculeName("MeOH"))
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");
+ methanol_.addExclusion(ParticleName("Me1"), ParticleName("O2"));
+ methanol_.addExclusion(ParticleName("Me1"), ParticleName("H3"));
+ methanol_.addExclusion(ParticleName("H3"), ParticleName("O2"));
}
Molecule MethanolMoleculeBuilder::methanolMolecule()
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()));
+ // Note this is a programming error as all particles should exclude at least themselves and empty topologies are not allowed.
+ const std::string message =
+ "No exclusions found in the " + molecule.name().value() + " molecule.";
+ assert((!exclusions.empty() && message.c_str()));
std::vector<gmx::ExclusionBlock> exclusionBlockPerMolecule =
detail::toGmxExclusionBlock(exclusions);
Topology TopologyBuilder::buildTopology()
{
+ assert((!(numParticles_ < 0) && "It should not be possible to have negative particles"));
+ if (numParticles_ == 0)
+ {
+ throw InputException("You cannot build a topology with no particles");
+ }
topology_.numParticles_ = numParticles_;
topology_.exclusions_ = createExclusionsListOfLists();
return std::get<0>(tup1) < std::get<0>(tup2);
};
+ // Note this is a programming error as all particles should exclude at least themselves and empty topologies are not allowed.
+ // Note also that this is also checked in the parent function with a more informative error message.
+ assert((!tupleList.empty() && "No exclusions found.\n"));
+
// 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;
#include <type_traits>
#include <vector>
-#include "nblib/basicdefinitions.h"
-#include "nblib/vector.h"
namespace nblib
{
template<auto...>
using void_value_t = void;
+template<class... Tuples>
+using tuple_cat_t = decltype(std::tuple_cat(Tuples{}...));
+
template<class T, class = void>
struct HasValueMember : std::false_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]
+/*! \brief Comparison meta function that compares T to Tuple[N]
+ *
+ * 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>>,
{
};
-//! 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>{}>{}>
+//! \brief Recursion to check the next field N+1
+template<int N, class T, class Tuple, template<int, class, class> class Comparison, class Match = void>
+struct MatchField_ : std::integral_constant<size_t, MatchField_<N + 1, T, Tuple, Comparison>{}>
{
};
-//! recursion stop when Comparison<N, T, Tuple>::value is true
+//! \brief 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>
+struct MatchField_<N, T, Tuple, Comparison, std::enable_if_t<Comparison<N, T, Tuple>{}>> :
+ 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
+
+/*! \brief Meta function to return the first 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, the size of Tuple will be returned.
+ * Note that the default comparison operation supplied here also matches if the type member Tuple[N]::type matches T
*/
-template<class T, class Tuple, template<int, class, class> class Comparison>
-struct MatchField : detail::MatchField_<0, T, Tuple, Comparison, Comparison<0, T, Tuple>{}>
+template<typename T, class TL, template<int, class, class> class Comparison = detail::MatchTypeOrTypeMember>
+struct FindIndex
{
};
-/*! \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
+/*! \brief Specialization to only enable this trait if TL has template parameters
+ *
+ * \tparam T a type to look for in the template parameters of TL
+ * \tparam TL a template template parameter, e.g. std::tuple or nblib::TypeList
+ * \tparam Ts template parameters of TL
+ * \tparam Comparison comparison operation
+ *
+ * Note that \a T is added to \a TL as a sentinel to terminate the recursion
+ * and prevent an out of bounds tuple access compiler error.
*/
-template<typename T, typename Tuple, template<int, class, class> class Comparison = detail::MatchTypeOrTypeMember>
-struct FindIndex : std::integral_constant<size_t, MatchField<T, Tuple, Comparison>{}>
+template<typename T, template<class...> class TL, class... Ts, template<int, class, class> class Comparison>
+struct FindIndex<T, TL<Ts...>, Comparison> : detail::MatchField_<0, T, std::tuple<Ts..., T>, 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
+/*! \brief Meta function to return the element in Tuple whose type matches T
+ *
+ * If there are more than one, the first occurrence will be returned
+ * If there is no such that, a compiler error is generated due to accessing
+ * the tuple out of bounds
+ */
template<typename T, typename Tuple>
decltype(auto) pickType(Tuple& tup)
{
- return std::get<FindIndex<T, Tuple>{}>(tup);
+ return std::get<FindIndex<T, std::decay_t<Tuple>>{}>(tup);
}
+//! \brief template meta function to determine whether T is contained in TL
+template<class T, class TL>
+struct Contains
+{
+};
+
+//! this formatting must be a bug in clang-format... should be:
+// struct Contains<T, TL<Ts...>> : std::bool_constant<FindIndex<T, TL<Ts...>>{} < sizeof...(Ts)>
+template<class T, template<class...> class TL, class... Ts>
+ struct Contains<T, TL<Ts...>> : std::bool_constant < FindIndex<T, TL<Ts...>>{}<sizeof...(Ts)>
+{
+};
+
+
//! Utility to call function with each element in tuple_
template<class F, class... Ts>
void for_each_tuple(F&& func, std::tuple<Ts...>& tuple_)
gmx_add_gtest_executable(
${exename}
CPP_SOURCE_FILES
+ internal.cpp
user.cpp
)
target_link_libraries(${exename} PRIVATE nblib_test_infrastructure 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 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 "nblib/tests/testhelpers.h"
+#include "nblib/util/internal.h"
+#include "nblib/util/user.h"
+
+namespace nblib
+{
+
+TEST(NblibInternalUtils, FindIndexTuple1)
+{
+ using TupleType = std::tuple<float>;
+
+ constexpr int floatIndex = FindIndex<float, TupleType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, TupleType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(1, outOfRange);
+}
+
+TEST(NblibInternalUtils, FindIndexTuple2)
+{
+ using TupleType = std::tuple<float, int>;
+
+ constexpr int floatIndex = FindIndex<float, TupleType>{};
+ constexpr int intIndex = FindIndex<int, TupleType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, TupleType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(1, intIndex);
+ EXPECT_EQ(2, outOfRange);
+}
+
+TEST(NblibInternalUtils, FindIndexTypeList1)
+{
+ using ListType = TypeList<float>;
+
+ constexpr int floatIndex = FindIndex<float, ListType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, ListType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(1, outOfRange);
+}
+
+TEST(NblibInternalUtils, FindIndexTypeList2)
+{
+ using ListType = TypeList<float, int>;
+
+ constexpr int floatIndex = FindIndex<float, ListType>{};
+ constexpr int intIndex = FindIndex<int, ListType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, ListType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(1, intIndex);
+ EXPECT_EQ(2, outOfRange);
+}
+
+
+TEST(NblibInternalUtils, Contains)
+{
+ using ListType = TypeList<float, int>;
+
+ constexpr bool hasFloat = Contains<float, ListType>{};
+ constexpr bool hasInt = Contains<int, ListType>{};
+ constexpr bool hasUint = Contains<unsigned, ListType>{};
+
+ EXPECT_TRUE(hasFloat);
+ EXPECT_TRUE(hasInt);
+ EXPECT_FALSE(hasUint);
+}
+
+TEST(NblibInternalUtils, FindIndexTupleRepeated)
+{
+ using TupleType = std::tuple<float, float, int>;
+
+ constexpr int floatIndex = FindIndex<float, TupleType>{};
+
+ constexpr int intIndex = FindIndex<int, TupleType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, TupleType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(2, intIndex);
+ EXPECT_EQ(3, outOfRange);
+}
+
+TEST(NblibInternalUtils, FindIndexTypeListRepeated)
+{
+ using TupleType = TypeList<float, float, int>;
+
+ constexpr int floatIndex = FindIndex<float, TupleType>{};
+
+ constexpr int intIndex = FindIndex<int, TupleType>{};
+
+ constexpr int outOfRange = FindIndex<unsigned, TupleType>{};
+
+ EXPECT_EQ(0, floatIndex);
+ EXPECT_EQ(2, intIndex);
+ EXPECT_EQ(3, outOfRange);
+}
+
+
+} // namespace nblib
namespace
{
-TEST(NBlibTest, checkNumericValues)
+TEST(NBlibTest, isRealValued)
{
std::vector<Vec3> vec;
vec.emplace_back(1., 1., 1.);
vec.emplace_back(2., 2., 2.);
- bool ret = checkNumericValues(vec);
+ bool ret = isRealValued(vec);
EXPECT_EQ(ret, true);
}
vec.emplace_back(NAN, NAN, NAN);
- bool ret = checkNumericValues(vec);
+ bool ret = isRealValued(vec);
EXPECT_EQ(ret, false);
}
vec.emplace_back(INFINITY, INFINITY, INFINITY);
- bool ret = checkNumericValues(vec);
+ bool ret = isRealValued(vec);
EXPECT_EQ(ret, false);
}
constexpr int N = 10;
std::vector<real> masses(N, 1.0);
auto out = generateVelocity(300.0, 1, masses);
- bool ret = checkNumericValues(out);
+ bool ret = isRealValued(out);
EXPECT_EQ(ret, true);
}
#include "nblib/util/user.h"
#include "gromacs/random/tabulatednormaldistribution.h"
#include "gromacs/random/threefry.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/fatalerror.h"
namespace nblib
}
//! Check within the container of gmx::RVecs for a NaN or inf
-bool checkNumericValues(const std::vector<Vec3>& values)
+bool isRealValued(gmx::ArrayRef<const Vec3> values)
{
for (auto val : values)
{
return true;
}
+void zeroCartesianArray(gmx::ArrayRef<Vec3> cartesianArray)
+{
+ std::fill(cartesianArray.begin(), cartesianArray.end(), Vec3{ 0, 0, 0 });
+}
+
} // namespace nblib
#include "nblib/basicdefinitions.h"
#include "nblib/vector.h"
+namespace gmx
+{
+template<typename T>
+class ArrayRef;
+} // namespace gmx
+
namespace nblib
{
std::vector<Vec3> generateVelocity(real Temperature, unsigned int seed, std::vector<real> const& masses);
//! Check within the container of gmx::RVecs for a NaN or inf
-bool checkNumericValues(const std::vector<Vec3>& values);
+bool isRealValued(gmx::ArrayRef<const Vec3> values);
+
+//! Zero a cartesian buffer
+void zeroCartesianArray(gmx::ArrayRef<Vec3> cartesianArray);
//! Used to ignore unused arguments of a lambda functions
inline void ignore_unused() {}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff