kerneloptions.h
nblib.h
particletype.h
+ ppmap.h
simulationstate.h
topology.h
topologyhelpers.h
DESTINATION include/nblib)
endif()
+add_subdirectory(listed_forces)
add_subdirectory(samples)
add_subdirectory(util)
--- /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>
+#
+
+if(GMX_INSTALL_NBLIB_API)
+ install(FILES
+ bondtypes.h
+ calculator.h
+ definitions.h
+ DESTINATION include/nblib)
+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 supported bondtypes
+ *
+ * We choose to forward comparison operations to the
+ * corresponding std::tuple comparison operations.
+ * In order to do that without temporary copies,
+ * we employ std::tie, which requires lvalues as input.
+ * For this reason, bond type parameter getters are implemented
+ * with a const lvalue reference return.
+ *
+ * \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_LISTEDFORCES_BONDTYPES_H
+#define NBLIB_LISTEDFORCES_BONDTYPES_H
+
+#include <array>
+
+#include "nblib/particletype.h"
+#include "nblib/ppmap.h"
+#include "nblib/util/user.h"
+
+namespace nblib
+{
+using Name = std::string;
+using ForceConstant = real;
+using EquilDistance = real;
+using Exponent = real;
+
+using Degrees = StrongType<real, struct DegreeParameter>;
+using Radians = StrongType<real, struct RadianParameter>;
+
+/*! \brief Basic template for interactions with 2 parameters named forceConstant and equilDistance
+ *
+ * \tparam Phantom unused template parameter for type distinction
+ *
+ * Distinct bond types can be generated from this template with using declarations
+ * and declared, but undefined structs. For example:
+ * using HarmonicBondType = TwoParameterInteraction<struct HarmonicBondTypeParameter>;
+ * Note that HarmonicBondTypeParameter does not have to be defined.
+ */
+template<class Phantom>
+class TwoParameterInteraction
+{
+public:
+ TwoParameterInteraction() = default;
+ TwoParameterInteraction(ForceConstant f, EquilDistance d) : forceConstant_(f), equilDistance_(d)
+ {
+ }
+
+ [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
+ [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
+
+private:
+ ForceConstant forceConstant_;
+ EquilDistance equilDistance_;
+};
+
+template<class Phantom>
+inline bool operator<(const TwoParameterInteraction<Phantom>& a, const TwoParameterInteraction<Phantom>& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistance())
+ < std::tie(b.forceConstant(), b.equilDistance());
+}
+
+template<class Phantom>
+inline bool operator==(const TwoParameterInteraction<Phantom>& a, const TwoParameterInteraction<Phantom>& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistance())
+ == std::tie(b.forceConstant(), b.equilDistance());
+}
+
+/*! \brief harmonic bond type
+ *
+ * It represents the interaction of the form
+ * V(r; forceConstant, equilDistance) = 0.5 * forceConstant * (r - equilDistance)^2
+ */
+using HarmonicBondType = TwoParameterInteraction<struct HarmonicBondTypeParameter>;
+
+
+/*! \brief GROMOS bond type
+ *
+ * It represents the interaction of the form
+ * V(r; forceConstant, equilDistance) = 0.25 * forceConstant * (r^2 - equilDistance^2)^2
+ */
+using G96BondType = TwoParameterInteraction<struct G96BondTypeParameter>;
+
+
+/*! \brief FENE bond type
+ *
+ * It represents the interaction of the form
+ * V(r; forceConstant, equilDistance) = - 0.5 * forceConstant * equilDistance^2 * log( 1 - (r / equilDistance)^2)
+ */
+using FENEBondType = TwoParameterInteraction<struct FENEBondTypeParameter>;
+
+
+/*! \brief Half-attractive quartic bond type
+ *
+ * It represents the interaction of the form
+ * V(r; forceConstant, equilDistance) = 0.5 * forceConstant * (r - equilDistance)^4
+ */
+using HalfAttractiveQuarticBondType =
+ TwoParameterInteraction<struct HalfAttractiveQuarticBondTypeParameter>;
+
+
+/*! \brief Cubic bond type
+ *
+ * It represents the interaction of the form
+ * V(r; quadraticForceConstant, cubicForceConstant, equilDistance) = quadraticForceConstant * (r -
+ * equilDistance)^2 + quadraticForceConstant * cubicForceConstant * (r - equilDistance)
+ */
+struct CubicBondType
+{
+ CubicBondType() = default;
+ CubicBondType(ForceConstant fq, ForceConstant fc, EquilDistance d) :
+ quadraticForceConstant_(fq),
+ cubicForceConstant_(fc),
+ equilDistance_(d)
+ {
+ }
+
+ [[nodiscard]] const ForceConstant& quadraticForceConstant() const
+ {
+ return quadraticForceConstant_;
+ }
+ [[nodiscard]] const ForceConstant& cubicForceConstant() const { return cubicForceConstant_; }
+ [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
+
+private:
+ ForceConstant quadraticForceConstant_;
+ ForceConstant cubicForceConstant_;
+ EquilDistance equilDistance_;
+};
+
+inline bool operator<(const CubicBondType& a, const CubicBondType& b)
+{
+ return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
+ < std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
+}
+
+inline bool operator==(const CubicBondType& a, const CubicBondType& b)
+{
+ return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
+ == std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
+}
+
+/*! \brief Morse bond type
+ *
+ * It represents the interaction of the form
+ * V(r; forceConstant, exponent, equilDistance) = forceConstant * ( 1 - exp( -exponent * (r - equilDistance))
+ */
+class MorseBondType
+{
+public:
+ MorseBondType() = default;
+ MorseBondType(ForceConstant f, Exponent e, EquilDistance d) :
+ forceConstant_(f),
+ exponent_(e),
+ equilDistance_(d)
+ {
+ }
+
+ [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
+ [[nodiscard]] const Exponent& exponent() const { return exponent_; }
+ [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
+
+private:
+ ForceConstant forceConstant_;
+ Exponent exponent_;
+ EquilDistance equilDistance_;
+};
+
+inline bool operator<(const MorseBondType& a, const MorseBondType& b)
+{
+ return std::tie(a.forceConstant(), a.exponent(), a.equilDistance())
+ < std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
+}
+
+inline bool operator==(const MorseBondType& a, const MorseBondType& b)
+{
+ return std::tie(a.forceConstant(), a.exponent(), a.equilDistance())
+ == std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
+}
+
+
+/*! \brief default angle type
+ *
+ * Note: the angle is always stored as radians internally
+ */
+struct DefaultAngle : public TwoParameterInteraction<struct DefaultAngleParameter>
+{
+ DefaultAngle() = default;
+ //! \brief construct from angle given in radians
+ DefaultAngle(Radians angle, ForceConstant f) :
+ TwoParameterInteraction<struct DefaultAngleParameter>{ f, angle }
+ {
+ }
+
+ //! \brief construct from angle given in degrees
+ DefaultAngle(Degrees angle, ForceConstant f) :
+ TwoParameterInteraction<struct DefaultAngleParameter>{ f, angle * DEG2RAD }
+ {
+ }
+};
+
+/*! \brief Proper Dihedral Implementation
+ */
+class ProperDihedral
+{
+public:
+ using Multiplicity = int;
+
+ ProperDihedral() = default;
+ ProperDihedral(Radians phi, ForceConstant f, Multiplicity m) :
+ phi_(phi),
+ forceConstant_(f),
+ multiplicity_(m)
+ {
+ }
+ ProperDihedral(Degrees phi, ForceConstant f, Multiplicity m) :
+ phi_(phi * DEG2RAD),
+ forceConstant_(f),
+ multiplicity_(m)
+ {
+ }
+
+ [[nodiscard]] const EquilDistance& equilDistance() const { return phi_; }
+ [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
+ [[nodiscard]] const Multiplicity& multiplicity() const { return multiplicity_; }
+
+private:
+ EquilDistance phi_;
+ ForceConstant forceConstant_;
+ Multiplicity multiplicity_;
+};
+
+inline bool operator<(const ProperDihedral& a, const ProperDihedral& b)
+{
+ return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
+ < std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
+}
+
+inline bool operator==(const ProperDihedral& a, const ProperDihedral& b)
+{
+ return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
+ == std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
+}
+
+
+/*! \brief Improper Dihedral Implementation
+ */
+struct ImproperDihedral : public TwoParameterInteraction<struct ImproperDihdedralParameter>
+{
+ ImproperDihedral() = default;
+ ImproperDihedral(Radians phi, ForceConstant f) :
+ TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi }
+ {
+ }
+ ImproperDihedral(Degrees phi, ForceConstant f) :
+ TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi * DEG2RAD }
+ {
+ }
+};
+
+/*! \brief Ryckaert-Belleman Dihedral Implementation
+ */
+class RyckaertBellemanDihedral
+{
+public:
+ RyckaertBellemanDihedral() = default;
+ RyckaertBellemanDihedral(real p1, real p2, real p3, real p4, real p5, real p6) :
+ parameters_{ p1, p2, p3, p4, p5, p6 }
+ {
+ }
+
+ const real& operator[](std::size_t i) const { return parameters_[i]; }
+
+ [[nodiscard]] const std::array<real, 6>& parameters() const { return parameters_; }
+
+ [[nodiscard]] std::size_t size() const { return parameters_.size(); }
+
+private:
+ std::array<real, 6> parameters_;
+};
+
+inline bool operator<(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
+{
+ return a.parameters() < b.parameters();
+}
+
+inline bool operator==(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
+{
+ return a.parameters() == b.parameters();
+}
+
+
+/*! \brief Type for 5-center interaction (C-MAP)
+ *
+ * Note: no kernels currently implemented
+ */
+class Default5Center
+{
+public:
+ Default5Center() = default;
+ Default5Center(Radians phi, Radians psi, ForceConstant fphi, ForceConstant fpsi) :
+ phi_(phi),
+ psi_(psi),
+ fphi_(fphi),
+ fpsi_(fpsi)
+ {
+ }
+
+ [[nodiscard]] const Radians& phi() const { return phi_; }
+ [[nodiscard]] const Radians& psi() const { return psi_; }
+ [[nodiscard]] const ForceConstant& fphi() const { return fphi_; }
+ [[nodiscard]] const ForceConstant& fpsi() const { return fpsi_; }
+
+private:
+ Radians phi_, psi_;
+ ForceConstant fphi_, fpsi_;
+};
+
+inline bool operator<(const Default5Center& a, const Default5Center& b)
+{
+ return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
+ < std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
+}
+
+inline bool operator==(const Default5Center& a, const Default5Center& b)
+{
+ return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
+ == std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
+}
+
+
+} // namespace nblib
+#endif // NBLIB_LISTEDFORCES_BONDTYPES_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.
+ */
+/*! \inpublicapi \file
+ * \brief
+ * Implements a force calculator based on GROMACS data structures.
+ *
+ * Intended for internal use inside the ForceCalculator.
+ *
+ * \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_LISTEDFORCES_CALCULATOR_H
+#define NBLIB_LISTEDFORCES_CALCULATOR_H
+
+#include <memory>
+#include <unordered_map>
+
+#include "nblib/listed_forces/definitions.h"
+
+namespace gmx
+{
+template<typename T>
+class ArrayRef;
+} // namespace gmx
+
+namespace nblib
+{
+class Box;
+class PbcHolder;
+template<class T>
+class ForceBuffer;
+
+/*! \internal \brief object to calculate listed forces
+ *
+ */
+class ListedForceCalculator
+{
+public:
+ using EnergyType = std::array<real, std::tuple_size<ListedInteractionData>::value>;
+
+ ListedForceCalculator(const ListedInteractionData& interactions,
+ size_t bufferSize,
+ int numThreads,
+ const Box& box);
+
+ /*! \brief Dispatch the listed force kernels and reduce the forces
+ *
+ * This function adds the computed listed forces to all values in the passed in forces buffer,
+ * so it can be regarded as an output only param. In case this is being used in a simulation
+ * that uses the same force buffer for both non-bonded and listed forces, this call should be
+ * made only after the compute() call from the non-bonded ForceCalculator
+ *
+ * This function also stores the forces and energies from listed interactions in the internal
+ * buffer of the ListedForceCalculator object
+ *
+ * \param[in] coordinates to be used for the force calculation
+ * \param[out] forces buffer to store the output forces
+ */
+ void compute(gmx::ArrayRef<const Vec3> coordinates, gmx::ArrayRef<Vec3> forces, bool usePbc = false);
+
+ //! Alternative overload with the energies in an output buffer
+ void compute(gmx::ArrayRef<const Vec3> coordinates,
+ gmx::ArrayRef<Vec3> forces,
+ EnergyType& energies,
+ bool usePbc = false);
+
+ /*! \brief We need to declare the destructor here to move the (still default) implementation
+ * to the .cpp file. Omitting this declaration would mean an inline destructor
+ * which can't compile because the unique_ptr dtor needs ~ForceBuffer, which is not available
+ * here because it's incomplete.
+ */
+ ~ListedForceCalculator();
+
+private:
+ int numThreads;
+
+ //! the main buffer to hold the final listed forces
+ std::vector<gmx::RVec> masterForceBuffer_;
+
+ //! holds the array of energies computed
+ EnergyType energyBuffer_;
+
+ //! holds the listed interactions split into groups for multithreading
+ std::vector<ListedInteractionData> threadedInteractions_;
+
+ //! reduction force buffers
+ std::vector<std::unique_ptr<ForceBuffer<gmx::RVec>>> threadedForceBuffers_;
+
+ //! PBC objects
+ std::unique_ptr<PbcHolder> pbcHolder_;
+
+ //! compute listed forces and energies, overwrites the internal buffers
+ void computeForcesAndEnergies(gmx::ArrayRef<const Vec3> x, bool usePbc = false);
+};
+
+} // namespace nblib
+
+#endif // NBLIB_LISTEDFORCES_CALCULATOR_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.
+ */
+/*! \inpublicapi \file
+ * \brief
+ * Definitions for supported nblib listed interaction data, such as bonds, angles, dihedrals, etc
+ *
+ * \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>
+ *
+ * A note on the preprocessor (PP) usage in this file:
+ *
+ * The PP macros defined here are used exclusively to generate
+ * template instantiations declarations of the form "extern template function(X)"
+ * in header files and "template function(X)" in .cpp files.
+ * These declarations do not affect the program logic in any way and neither are they
+ * required to read and understand the behavior of the code as they do not
+ * result in any executable instructions.
+ * In fact, it would even be technically possible to omit these PP generated
+ * declarations in the header files and replace them with an unused static function
+ * in the .cpp file that calls the template function in question
+ * (e.g. Molecule::addInteraction) once with each type from the variadic template
+ * TypeLists declared in this file. This would be enough to create the required instantiations.
+ * It would, however, create more work for the compiler which then has to instantiate the
+ * templates in the header in each translation unit where the header is included.
+ * Doing this results in a compiler warning.
+ *
+ */
+#ifndef NBLIB_LISTEDFORCES_DEFINITIONS_H
+#define NBLIB_LISTEDFORCES_DEFINITIONS_H
+
+#include "nblib/util/user.h"
+#include "bondtypes.h"
+
+namespace nblib
+{
+
+//***********************************************************************************
+
+/*! \brief These macros define what interaction types are supported in
+ * -Molecule
+ * -Topology
+ * -ListedForceCalculator
+ *
+ * To enable force calculation for your new interaction type that you've added to bondtypes.h,
+ * list your new type here under the appropriate category and make sure that you've added
+ * a kernel in kernels.hpp
+ */
+
+#define SUPPORTED_TWO_CENTER_TYPES \
+ HarmonicBondType, G96BondType, CubicBondType, FENEBondType, HalfAttractiveQuarticBondType
+
+#define SUPPORTED_THREE_CENTER_TYPES DefaultAngle
+
+#define SUPPORTED_FOUR_CENTER_TYPES ProperDihedral, ImproperDihedral, RyckaertBellemanDihedral
+
+#define SUPPORTED_FIVE_CENTER_TYPES Default5Center
+
+//***********************************************************************************
+
+#define SUPPORTED_LISTED_TYPES \
+ SUPPORTED_TWO_CENTER_TYPES, SUPPORTED_THREE_CENTER_TYPES, SUPPORTED_FOUR_CENTER_TYPES, \
+ SUPPORTED_FIVE_CENTER_TYPES
+
+#define NBLIB_ALWAYS_INLINE __attribute((always_inline))
+
+//! \brief encodes the number of integers needed to represent 2-center interactions (bonds, pairs)
+using TwoCenterInteractionIndex = std::array<int, 3>;
+//! \brief encodes the number of integers needed to represent 3-center interactions (angles)
+using ThreeCenterInteractionIndex = std::array<int, 4>;
+//! \brief encodes the number of integers needed to represent 4-center interactions (dihedrals)
+using FourCenterInteractionIndex = std::array<int, 5>;
+//! \brief encodes the number of integers needed to represent 5-center interactions (CMAP)
+using FiveCenterInteractionIndex = std::array<int, 6>;
+
+//! \brief data type for pairwise interactions, e.g. bonds
+template<class TwoCenterType>
+struct TwoCenterData
+{
+ using type = TwoCenterType;
+
+ // tuple format: <particleID i, particleID j, TwoCenterInstanceIndex>
+ std::vector<TwoCenterInteractionIndex> indices;
+ // vector of unique TwoCenterType instances
+ std::vector<TwoCenterType> parameters;
+};
+
+//! \brief data type for three-center interactions, e.g. angles
+template<class ThreeCenterType>
+struct ThreeCenterData
+{
+ using type = ThreeCenterType;
+
+ // tuple format: <particleID i, particleID j, particleID k, ThreeCenterInstanceIndex>
+ std::vector<ThreeCenterInteractionIndex> indices;
+ // vector of unique ThreeCenterType instances
+ std::vector<ThreeCenterType> parameters;
+};
+
+//! \brief data type for four-center interactions, e.g. dihedrals
+template<class FourCenterType>
+struct FourCenterData
+{
+ using type = FourCenterType;
+
+ // tuple format: <particleID i, particleID j, particleID k, particleID l, FourCenterInstanceIndex>
+ std::vector<FourCenterInteractionIndex> indices;
+ // vector of unique FiveCenterType instances
+ std::vector<FourCenterType> parameters;
+};
+
+//! \brief data type for five-center interactions, e.g. CMAP
+template<class FiveCenterType>
+struct FiveCenterData
+{
+ using type = FiveCenterType;
+
+ // tuple format: <particleID i, particleID j, particleID k, particleID l, particleID m, FiveCenterInstanceIndex>
+ std::vector<FiveCenterInteractionIndex> indices;
+ // vector of unique FiveCenterType instances
+ std::vector<FiveCenterType> parameters;
+};
+
+
+using SupportedTwoCenterTypes = TypeList<SUPPORTED_TWO_CENTER_TYPES>;
+// std::tuple<TwoCenterData<TwoCenterType1>, ...>
+using TwoCenterInteractionData = Reduce<std::tuple, Map<TwoCenterData, SupportedTwoCenterTypes>>;
+
+using SupportedThreeCenterTypes = TypeList<SUPPORTED_THREE_CENTER_TYPES>;
+// std::tuple<AngleData<ThreeCenterType1>, ...>
+using ThreeCenterInteractionData = Reduce<std::tuple, Map<ThreeCenterData, SupportedThreeCenterTypes>>;
+
+using SupportedFourCenterTypes = TypeList<SUPPORTED_FOUR_CENTER_TYPES>;
+// std::tuple<FourCenterData<FourCenterType1>, ...>
+using FourCenterInteractionData = Reduce<std::tuple, Map<FourCenterData, SupportedFourCenterTypes>>;
+
+using SupportedFiveCenterTypes = TypeList<SUPPORTED_FIVE_CENTER_TYPES>;
+// std::tuple<FiveCenterData<FiveCenterType1>, ...>
+using FiveCenterInteractionData = Reduce<std::tuple, Map<FiveCenterData, SupportedFiveCenterTypes>>;
+
+//! This is the complete type that holds all listed interaction data
+using ListedInteractionData = decltype(std::tuple_cat(TwoCenterInteractionData{},
+ ThreeCenterInteractionData{},
+ FourCenterInteractionData{},
+ FiveCenterInteractionData{}));
+} // namespace nblib
+#endif // NBLIB_LISTEDFORCES_DEFINITIONS_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>
+#
+
+# Make a static library for test infrastructure code that we re-use
+# in multiple test executables across the repository.
+
+set(testname "NbLibListedForcesTests")
+set(exename "nblib-listed-forces-test")
+
+gmx_add_gtest_executable(
+ ${exename}
+ CPP_SOURCE_FILES
+ # files with code for tests
+ bondtypes.cpp
+)
+target_link_libraries(${exename} PRIVATE nblib_test_infrastructure nblib)
+gmx_register_gtest_test(${testname} ${exename} INTEGRATION_TEST)
+add_dependencies(check-nblib ${exename})
--- /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 box 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/listed_forces/bondtypes.h"
+#include "nblib/util/internal.h"
+
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+
+namespace test_detail
+{
+
+template<class B>
+void testTwoParameterBondEquality(const B& deduceType)
+{
+ ignore_unused(deduceType);
+ B a(1, 2);
+ B b(1, 2);
+ EXPECT_TRUE(a == b);
+
+ B c(1, 3);
+ EXPECT_FALSE(a == c);
+}
+
+template<class B>
+void testThreeParameterBondEquality(const B& deduceType)
+{
+ ignore_unused(deduceType);
+ B a(1, 2, 3);
+ B b(1, 2, 3);
+ EXPECT_TRUE(a == b);
+
+ B c(2, 3, 4);
+ EXPECT_FALSE(a == c);
+}
+
+template<class B>
+void testTwoParameterBondLessThan(const B& deduceType)
+{
+ ignore_unused(deduceType);
+ B a(1, 2);
+ B b(1, 3);
+ EXPECT_TRUE(a < b);
+ EXPECT_FALSE(b < a);
+
+ B c(1, 2);
+ B d(1, 2);
+ EXPECT_FALSE(c < d);
+
+ B e(2, 1);
+ B f(3, 1);
+ EXPECT_TRUE(e < f);
+ EXPECT_FALSE(f < e);
+}
+
+template<class B>
+void testThreeParameterBondLessThan(const B& deduceType)
+{
+ ignore_unused(deduceType);
+ B a(1, 2, 1);
+ B b(1, 3, 1);
+ EXPECT_TRUE(a < b);
+ EXPECT_FALSE(b < a);
+
+ B c(1, 2, 3);
+ B d(1, 2, 3);
+ EXPECT_FALSE(c < d);
+
+ B e(4, 1, 3);
+ B f(5, 1, 2);
+ EXPECT_TRUE(e < f);
+ EXPECT_FALSE(f < e);
+}
+
+} // namespace test_detail
+
+TEST(NBlibTest, BondTypesOperatorEqualWorks)
+{
+ auto bondList3 = std::make_tuple(HarmonicBondType(), G96BondType(), FENEBondType(),
+ HalfAttractiveQuarticBondType());
+ for_each_tuple([](const auto& b) { test_detail::testTwoParameterBondEquality(b); }, bondList3);
+
+ auto bondList4 = std::make_tuple(CubicBondType(), MorseBondType());
+ for_each_tuple([](const auto& b) { test_detail::testThreeParameterBondEquality(b); }, bondList4);
+}
+
+TEST(NBlibTest, BondTypesLessThanWorks)
+{
+ auto bondList3 = std::make_tuple(HarmonicBondType(), G96BondType(), FENEBondType(),
+ HalfAttractiveQuarticBondType());
+ for_each_tuple([](const auto& b) { test_detail::testTwoParameterBondLessThan(b); }, bondList3);
+
+ auto bondList4 = std::make_tuple(CubicBondType(), MorseBondType());
+ for_each_tuple([](const auto& b) { test_detail::testThreeParameterBondLessThan(b); }, bondList4);
+}
+
+
+} // 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
+ * These traits defined here for supported nblib listed interaction data types
+ * are used to control the dataflow in dataflow.hpp
+ *
+ * \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_LISTEDFORCES_TRAITS_H
+#define NBLIB_LISTEDFORCES_TRAITS_H
+
+#include <numeric>
+
+#include "nblib/util/internal.h"
+#include "bondtypes.h"
+#include "definitions.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+
+template<class InteractionType, class = void>
+struct CoordinateIndex_
+{
+};
+
+template<class InteractionType>
+struct CoordinateIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedTwoCenterTypes>{}>>
+{
+ typedef std::array<int, 2> type;
+};
+
+template<class InteractionType>
+struct CoordinateIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedThreeCenterTypes>{}>>
+{
+ typedef std::array<int, 3> type;
+};
+
+template<class InteractionType>
+struct CoordinateIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedFourCenterTypes>{}>>
+{
+ typedef std::array<int, 4> type;
+};
+
+template<class InteractionType>
+struct CoordinateIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedFiveCenterTypes>{}>>
+{
+ typedef std::array<int, 5> type;
+};
+
+} // namespace detail
+
+/*! \brief traits class to determine the coordinate index type for InteractionType
+ * \internal
+ *
+ * \tparam InteractionCategory
+ */
+template<class InteractionType>
+using CoordinateIndex = typename detail::CoordinateIndex_<InteractionType>::type;
+
+
+namespace detail
+{
+
+template<class InteractionType, class = void>
+struct InteractionIndex_
+{
+};
+
+template<class InteractionType>
+struct InteractionIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedTwoCenterTypes>{}>>
+{
+ typedef TwoCenterInteractionIndex type;
+};
+
+template<class InteractionType>
+struct InteractionIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedThreeCenterTypes>{}>>
+{
+ typedef ThreeCenterInteractionIndex type;
+};
+
+template<class InteractionType>
+struct InteractionIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedFourCenterTypes>{}>>
+{
+ typedef FourCenterInteractionIndex type;
+};
+
+template<class InteractionType>
+struct InteractionIndex_<InteractionType, std::enable_if_t<Contains<InteractionType, SupportedFiveCenterTypes>{}>>
+{
+ typedef FiveCenterInteractionIndex type;
+};
+
+} // namespace detail
+
+/*! \brief traits class to determine the InteractionIndex type for InteractionType
+ * \internal
+ *
+ * \tparam InteractionType
+ */
+template<class InteractionType>
+using InteractionIndex = typename detail::InteractionIndex_<InteractionType>::type;
+
+
+template<class I, class = void>
+struct HasTwoCenterAggregate : std::false_type
+{
+};
+
+template<class I>
+struct HasTwoCenterAggregate<I, std::void_t<typename I::TwoCenterAggregateType>> : std::true_type
+{
+};
+
+template<class I, class = void>
+struct HasThreeCenterAggregate : std::false_type
+{
+};
+
+template<class I>
+struct HasThreeCenterAggregate<I, std::void_t<typename I::ThreeCenterAggregateType>> : std::true_type
+{
+};
+
+//! \internal \brief determines the energy storage location of the carrier part for InteractionTypes without aggregates
+template<class InteractionType, class = void>
+struct CarrierIndex :
+ std::integral_constant<size_t, FindIndex<InteractionType, ListedInteractionData>{}>
+{
+};
+
+//! \internal \brief determines the energy storage location of the carrier part for InteractionTypes with aggregates
+template<class InteractionType>
+struct CarrierIndex<InteractionType, std::void_t<typename InteractionType::CarrierType>> :
+ std::integral_constant<size_t, FindIndex<typename InteractionType::CarrierType, ListedInteractionData>{}>
+{
+};
+
+//! \internal \brief determines the energy storage location of the 2-C aggregate part for InteractionTypes without aggregates
+template<class InteractionType, class = void>
+struct TwoCenterAggregateIndex : std::integral_constant<size_t, 0>
+{
+};
+
+//! \internal \brief determines the energy storage location of the 2-C aggregate part for InteractionTypes with 2-C aggregates
+template<class InteractionType>
+struct TwoCenterAggregateIndex<InteractionType, std::void_t<typename InteractionType::TwoCenterAggregateType>> :
+ std::integral_constant<size_t, FindIndex<typename InteractionType::TwoCenterAggregateType, ListedInteractionData>{}>
+{
+};
+
+//! \internal \brief determines the energy storage location of the 3-C aggregate part for InteractionTypes without aggregates
+template<class InteractionType, class = void>
+struct ThreeCenterAggregateIndex : std::integral_constant<size_t, 0>
+{
+};
+
+//! \internal \brief determines the energy storage location of the 3-C aggregate part for InteractionTypes with 3-C aggregates
+template<class InteractionType>
+struct ThreeCenterAggregateIndex<InteractionType, std::void_t<typename InteractionType::ThreeCenterAggregateType>> :
+ std::integral_constant<size_t, FindIndex<typename InteractionType::ThreeCenterAggregateType, ListedInteractionData>{}>
+{
+};
+
+/*! \brief return type to hold the energies of the different overloads of "dispatchInteraction"
+ * \internal
+ *
+ * \tparam T
+ */
+template<class T>
+class KernelEnergy
+{
+public:
+ KernelEnergy() : energies_{ 0, 0, 0, 0 } {}
+
+ T& carrier() { return energies_[0]; }
+ const T& carrier() const { return energies_[0]; }
+
+ T& twoCenterAggregate() { return energies_[1]; }
+ const T& twoCenterAggregate() const { return energies_[1]; }
+
+ T& threeCenterAggregate() { return energies_[2]; }
+ const T& threeCenterAggregate() const { return energies_[2]; }
+
+ T& freeEnergyDerivative() { return energies_[3]; }
+ const T& freeEnergyDerivative() const { return energies_[3]; }
+
+ KernelEnergy& operator+=(const KernelEnergy& other)
+ {
+ for (size_t i = 0; i < energies_.size(); ++i)
+ {
+ energies_[i] += other.energies_[i];
+ }
+ return *this;
+ }
+
+ operator T() const { return std::accumulate(begin(energies_), end(energies_), T{}); }
+
+private:
+ std::array<T, 4> energies_;
+};
+
+template<class BasicVector>
+using BasicVectorValueType_t = std::remove_all_extents_t<typename BasicVector::RawArray>;
+
+} // namespace nblib
+#endif // NBLIB_LISTEDFORCES_TRAITS_H
#include <unordered_map>
#include <vector>
+#include "nblib/listed_forces/definitions.h"
#include "nblib/particletype.h"
namespace nblib
class Molecule final
{
+ //! \brief string based listed interaction data type used during construction
+ template<class TwoCenterType>
+ struct TwoCenterData
+ {
+ using type = TwoCenterType;
+
+ std::vector<TwoCenterType> interactionTypes_;
+ std::vector<std::tuple<ParticleName, ResidueName, ParticleName, ResidueName>> interactions_;
+ };
+
+ template<class ThreeCenterType>
+ struct ThreeCenterData
+ {
+ using type = ThreeCenterType;
+
+ std::vector<ThreeCenterType> interactionTypes_;
+ std::vector<std::tuple<ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName>> interactions_;
+ };
+
+ template<class FourCenter>
+ struct FourCenterDataHolder
+ {
+ using type = FourCenter;
+
+ std::vector<FourCenter> interactionTypes_;
+ std::vector<std::tuple<ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName>> interactions_;
+ };
+
+ template<class FiveCenter>
+ struct FiveCenterDataHolder
+ {
+ using type = FiveCenter;
+
+ std::vector<FiveCenter> interactionTypes_;
+ std::vector<std::tuple<ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName, ParticleName, ResidueName>>
+ interactions_;
+ };
+
+ // BondContainerTypes is TypeList<TwoCenterData<HarmonicBondType>, ...>
+ using TwoCenterContainerTypes = Map<TwoCenterData, SupportedTwoCenterTypes>;
+
+ using ThreeCenterContainerTypes = Map<ThreeCenterData, SupportedThreeCenterTypes>;
+
+ using FourCenterContainerTypes = Map<FourCenterDataHolder, SupportedFourCenterTypes>;
+
+ using FiveCenterContainerTypes = Map<FiveCenterDataHolder, SupportedFiveCenterTypes>;
+
+ // InteractionTuple is std::tuple<TwoCenterData<HarmonicBondType>, ...>
+ using InteractionTuple = decltype(std::tuple_cat(Reduce<std::tuple, TwoCenterContainerTypes>{},
+ Reduce<std::tuple, ThreeCenterContainerTypes>{},
+ Reduce<std::tuple, FourCenterContainerTypes>{},
+ Reduce<std::tuple, FiveCenterContainerTypes>{}));
+
public:
explicit Molecule(MoleculeName moleculeName);
//! Specify an exclusion with particle names that have been added to molecule
void addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude);
+ // Add various types of interactions to the molecule
+ // Note: adding an interaction type not listed in SUPPORTED_TWO_CENTER_TYPES results in a compilation error
+
+ //! For 2-particle interactions such as harmonic bonds
+ template<class Interaction>
+ void addInteraction(const ParticleName& particleNameI,
+ const ResidueName& residueNameI,
+ const ParticleName& particleNameJ,
+ const ResidueName& residueNameJ,
+ const Interaction& interaction);
+
+ //! Add 2-particle interactions with the default residue name
+ template<class Interaction>
+ void addInteraction(const ParticleName& particleNameI,
+ const ParticleName& particleNameJ,
+ const Interaction& interaction);
+
+ //! For 3-particle interactions such as angles
+ template<class Interaction>
+ void addInteraction(const ParticleName& particleNameI,
+ const ResidueName& residueNameI,
+ const ParticleName& particleNameJ,
+ const ResidueName& residueNameJ,
+ const ParticleName& particleNameK,
+ const ResidueName& residueNameK,
+ const Interaction& interaction);
+
+ //! Add 3-particle interactions with the default residue name
+ template<class Interaction>
+ void addInteraction(const ParticleName& particleNameI,
+ const ParticleName& particleNameJ,
+ const ParticleName& particleNameK,
+ const Interaction& interaction);
+
//! The number of molecules
int numParticlesInMolecule() const;
//! returns a sorted vector containing no duplicates of particles to exclude by indices
std::vector<std::tuple<int, int>> getExclusions() const;
+ //! Return all interactions stored in Molecule
+ const InteractionTuple& interactionData() const;
+
//! Return name of ith particle
ParticleName particleName(int i) const;
//! we cannot efficiently compute indices during the build-phase
//! so we delay the conversion until TopologyBuilder requests it
std::vector<std::tuple<std::string, std::string, std::string, std::string>> exclusionsByName_;
+
+ //! collection of data for all types of interactions
+ InteractionTuple interactionData_;
};
+//! \cond DO_NOT_DOCUMENT
+#define ADD_INTERACTION_EXTERN_TEMPLATE(x) \
+ extern template void Molecule::addInteraction( \
+ const ParticleName& particleNameI, const ResidueName& residueNameI, \
+ const ParticleName& particleNameJ, const ResidueName& residueNameJ, const x& interaction);
+MAP(ADD_INTERACTION_EXTERN_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
+#undef ADD_INTERACTION_EXTERN_TEMPLATE
+
+#define ADD_INTERACTION_EXTERN_TEMPLATE(x) \
+ extern template void Molecule::addInteraction( \
+ const ParticleName& particleNameI, const ParticleName& particleNameJ, const x& interaction);
+MAP(ADD_INTERACTION_EXTERN_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
+#undef ADD_INTERACTION_EXTERN_TEMPLATE
+
+#define ADD_INTERACTION_EXTERN_TEMPLATE(x) \
+ extern template void Molecule::addInteraction( \
+ const ParticleName& particleNameI, const ResidueName& residueNameI, \
+ const ParticleName& particleNameJ, const ResidueName& residueNameJ, \
+ const ParticleName& particleNameK, const ResidueName& residueNameK, const x& interaction);
+MAP(ADD_INTERACTION_EXTERN_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
+#undef ADD_INTERACTION_EXTERN_TEMPLATE
+
+#define ADD_INTERACTION_EXTERN_TEMPLATE(x) \
+ extern template void Molecule::addInteraction( \
+ const ParticleName& particleNameI, const ParticleName& particleNameJ, \
+ const ParticleName& particleNameK, const x& interaction);
+MAP(ADD_INTERACTION_EXTERN_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
+#undef ADD_INTERACTION_EXTERN_TEMPLATE
+//! \endcond
+
} // namespace nblib
#endif // NBLIB_MOLECULES_H
#include "nblib/integrator.h"
#include "nblib/interactions.h"
#include "nblib/kerneloptions.h"
+#include "nblib/listed_forces/bondtypes.h"
+#include "nblib/listed_forces/calculator.h"
+#include "nblib/listed_forces/definitions.h"
#include "nblib/molecules.h"
#include "nblib/particletype.h"
+#include "nblib/ppmap.h"
#include "nblib/simulationstate.h"
#include "nblib/topology.h"
#include "nblib/topologyhelpers.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.
+ */
+/*
+ * Copyright (C) 2012 William Swanson
+ *
+ * Permission is hereby granted, free of charge, to any person
+ * obtaining a copy of this software and associated documentation
+ * files (the "Software"), to deal in the Software without
+ * restriction, including without limitation the rights to use, copy,
+ * modify, merge, publish, distribute, sublicense, and/or sell copies
+ * of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ * Except as contained in this notice, the names of the authors or
+ * their institutions shall not be used in advertising or otherwise to
+ * promote the sale, use or other dealings in this Software without
+ * prior written authorization from the authors.
+ */
+
+/*! \inpublicapi \file
+ * \brief
+ * Provides MAP and MAP_LIST to apply a macro to a variadic argument list
+ *
+ * The MAP and MAP_LIST macros implement calling a supplied macro with
+ * all of the subsequent arguments. For example:
+ * MAP(macro, x, y, z) expands to macro(x) macro(y) macro(z) while
+ * MAP_LIST(macro, x, y, z) expands to macro(x), macro(y), macro(z)
+ *
+ * Due to the limitations of the preprocessor, it is unfortunately not
+ * possible to implement this functionality in a more straight-forward way.
+ * Since this use-case is not too uncommon, Boost for example implements
+ * BOOST_PP_SEQ_FOR_EACH which provides equivalent functionality implemented
+ * with the same technique, but is more comprehensive in scope,
+ * beyond what's required here.
+ *
+ * A discussion of how and why this macro works can be found here:
+ * https://stackoverflow.com/questions/27765387/distributing-an-argument-in-a-variadic-macro
+ * and the original repository of this implementation is this one:
+ * https://github.com/swansontec/map-macro
+ * It also contains some useful explanations of how it works.
+ */
+
+#ifndef NBLIB_PPMAP_H
+#define NBLIB_PPMAP_H
+
+#define EVAL0(...) __VA_ARGS__
+#define EVAL1(...) EVAL0(EVAL0(EVAL0(__VA_ARGS__)))
+#define EVAL2(...) EVAL1(EVAL1(EVAL1(__VA_ARGS__)))
+#define EVAL3(...) EVAL2(EVAL2(EVAL2(__VA_ARGS__)))
+#define EVAL4(...) EVAL3(EVAL3(EVAL3(__VA_ARGS__)))
+#define EVAL(...) EVAL4(EVAL4(EVAL4(__VA_ARGS__)))
+
+#define MAP_END(...)
+#define MAP_OUT
+#define MAP_COMMA ,
+
+#define MAP_GET_END2() 0, MAP_END
+#define MAP_GET_END1(...) MAP_GET_END2
+#define MAP_GET_END(...) MAP_GET_END1
+#define MAP_NEXT0(test, next, ...) next MAP_OUT
+#define MAP_NEXT1(test, next) MAP_NEXT0(test, next, 0)
+#define MAP_NEXT(test, next) MAP_NEXT1(MAP_GET_END test, next)
+
+#define MAP0(f, x, peek, ...) f(x) MAP_NEXT(peek, MAP1)(f, peek, __VA_ARGS__)
+#define MAP1(f, x, peek, ...) f(x) MAP_NEXT(peek, MAP0)(f, peek, __VA_ARGS__)
+
+#define MAP_LIST_NEXT1(test, next) MAP_NEXT0(test, MAP_COMMA next, 0)
+#define MAP_LIST_NEXT(test, next) MAP_LIST_NEXT1(MAP_GET_END test, next)
+
+#define MAP_LIST0(f, x, peek, ...) f(x) MAP_LIST_NEXT(peek, MAP_LIST1)(f, peek, __VA_ARGS__)
+#define MAP_LIST1(f, x, peek, ...) f(x) MAP_LIST_NEXT(peek, MAP_LIST0)(f, peek, __VA_ARGS__)
+
+/**
+ * Applies the function macro `f` to each of the remaining parameters.
+ */
+#define MAP(f, ...) EVAL(MAP1(f, __VA_ARGS__, ()()(), ()()(), ()()(), 0))
+
+/**
+ * Applies the function macro `f` to each of the remaining parameters and
+ * inserts commas between the results.
+ */
+#define MAP_LIST(f, ...) EVAL(MAP_LIST1(f, __VA_ARGS__, ()()(), ()()(), ()()(), 0))
+
+
+/** The PP_NARG macro returns the number of arguments that have been
+ * passed to it.
+ */
+#define PP_NARG(...) PP_NARG_(__VA_ARGS__, PP_RSEQ_N())
+#define PP_NARG_(...) PP_ARG_N(__VA_ARGS__)
+#define PP_ARG_N(_1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, _16, _17, _18, \
+ _19, _20, _21, _22, _23, _24, _25, _26, _27, _28, _29, _30, _31, _32, _33, _34, \
+ _35, _36, _37, _38, _39, _40, _41, _42, _43, _44, _45, _46, _47, _48, _49, _50, \
+ _51, _52, _53, _54, _55, _56, _57, _58, _59, _60, _61, _62, _63, N, ...) \
+ N
+#define PP_RSEQ_N() \
+ 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, \
+ 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, \
+ 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0
+
+/** MAP_ENUMERATE macro:
+ * MAP_ENUMERATE(action, args...)
+ * like MAP, calls action with each argument, but also forwards the index of the argument to action
+ */
+#define FE_0(WHAT)
+#define FE_1(WHAT, N, X) WHAT(X, N - 1) // NOLINT bugprone-macro-parentheses
+#define FE_2(WHAT, N, X, ...) WHAT(X, N - 2) FE_1(WHAT, N, __VA_ARGS__)
+#define FE_3(WHAT, N, X, ...) WHAT(X, N - 3) FE_2(WHAT, N, __VA_ARGS__)
+#define FE_4(WHAT, N, X, ...) WHAT(X, N - 4) FE_3(WHAT, N, __VA_ARGS__)
+#define FE_5(WHAT, N, X, ...) WHAT(X, N - 5) FE_4(WHAT, N, __VA_ARGS__)
+
+#define GET_MACRO(_0, _1, _2, _3, _4, _5, NAME, ...) NAME
+#define MAP_ENUMERATE(action, ...) \
+ GET_MACRO(_0, __VA_ARGS__, FE_5, FE_4, FE_3, FE_2, FE_1, FE_0, ) \
+ (action, PP_NARG(__VA_ARGS__), __VA_ARGS__)
+
+#endif // NBLIB_PPMAP_H
#include <vector>
#include "nblib/interactions.h"
+#include "nblib/listed_forces/definitions.h"
#include "nblib/molecules.h"
#include "nblib/topologyhelpers.h"
//! Returns a map of non-bonded force parameters indexed by ParticleType names
NonBondedInteractionMap getNonBondedInteractionMap() const;
+ //! Returns the interaction data
+ ListedInteractionData getInteractionData() const;
+
//! Returns the combination rule used to generate the NonBondedInteractionMap
CombinationRule getCombinationRule() const;
detail::ParticleSequencer particleSequencer_;
//! Map that should hold all nonbonded interactions for all particle types
NonBondedInteractionMap nonBondedInteractionMap_;
+ //! data about bonds for all supported types
+ ListedInteractionData interactionData_;
//! Combination Rule used to generate the nonbonded interactions
CombinationRule combinationRule_;
};
//! Builds a GROMACS-compliant performant exclusions list aggregating exclusions from all molecules
gmx::ListOfLists<int> createExclusionsListOfLists() const;
+ //! Gather interaction data from molecules
+ ListedInteractionData createInteractionData(const detail::ParticleSequencer&);
+
//! Helper function to extract quantities like mass, charge, etc from the system
template<typename T, class Extractor>
std::vector<T> extractParticleTypeQuantity(Extractor&& extractor);
#include <vector>
#include "gromacs/utility/listoflists.h"
+#include "nblib/listed_forces/traits.h"
#include "nblib/molecules.h"
namespace gmx
//! Add offset to all indices in inBlock
std::vector<gmx::ExclusionBlock> offsetGmxBlock(std::vector<gmx::ExclusionBlock> inBlock, int offset);
+/*!
+ * \brief
+ * Extract all interactions of type I from a vector of molecules. The second argument tuple element
+ * specifies multiples of the molecule given as first tuple element. Let (S, I) denote the return
+ * value tuple. Then J[i] = I[S[i]] for all i in 0...S.size() is the full sequence of BondType
+ * instances as they occur in the input tuple
+ *
+ */
+template<class I>
+std::tuple<std::vector<size_t>, std::vector<I>>
+collectInteractions(const std::vector<std::tuple<Molecule, int>>&);
+
+#define COLLECT_BONDS_EXTERN_TEMPLATE(x) \
+ extern template std::tuple<std::vector<size_t>, std::vector<x>> collectInteractions( \
+ const std::vector<std::tuple<Molecule, int>>&);
+MAP(COLLECT_BONDS_EXTERN_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
+#undef COLLECT_BONDS_EXTERN_TEMPLATE
+
+/*!
+ * \brief
+ * Return a list of unique BondType instances U and an index list S of size aggregatedBonds.size()
+ * such that the BondType instance at aggregatedBonds[i] is equal to U[S[i]]
+ * returns std::tuple(S, U)
+ *
+ */
+template<class I>
+std::tuple<std::vector<size_t>, std::vector<I>> eliminateDuplicateInteractions(const std::vector<I>& collectedBonds);
+
+/// \cond DO_NOT_DOCUMENT
+#define ELIMINATE_DUPLICATE_EXTERN_TEMPLATE(x) \
+ extern template std::tuple<std::vector<size_t>, std::vector<x>> eliminateDuplicateInteractions( \
+ const std::vector<x>& collectedBonds);
+MAP(ELIMINATE_DUPLICATE_EXTERN_TEMPLATE, SUPPORTED_LISTED_TYPES)
+#undef ELIMINATE_DUPLICATE_EXTERN_TEMPLATE
+/// \endcond
+
//! Helper class for Topology to keep track of particle IDs
class ParticleSequencer
{
DataType data_;
};
+//!
+template<class B>
+std::vector<CoordinateIndex<B>> sequenceIDs(const std::vector<std::tuple<Molecule, int>>&,
+ const detail::ParticleSequencer&);
+
+/// \cond DO_NOT_DOCUMENT
+#define SEQUENCE_PAIR_ID_EXTERN_TEMPLATE(x) \
+ extern template std::vector<CoordinateIndex<x>> sequenceIDs<x>( \
+ const std::vector<std::tuple<Molecule, int>>&, const detail::ParticleSequencer&);
+MAP(SEQUENCE_PAIR_ID_EXTERN_TEMPLATE, SUPPORTED_LISTED_TYPES)
+#undef SEQUENCE_PAIR_ID_EXTERN_TEMPLATE
+/// \endcond
+
} // namespace detail
} // namespace nblib
--- /dev/null
+Design goals and motivation for the data format of bonded forces in NB-LIB
+--------------------------------------------------------------------------
+
+
+The current format for listed forces in GROMACS looks like this:
+
+.. code:: cpp
+
+ struct InteractionDefinitions
+ {
+ std::vector<t_iparams> iparams;
+ std::array<std::vector<int>, F_NRE> il;
+ };
+
+The format covers all interaction types, i.e. \ ``t_iparams`` is a union
+type which can hold the parameters of any type.
+The other member called ``il`` contains the
+indices for each interaction type, where ``F_NRE`` is the number of
+interaction types that GROMACS supports. More precisely, each
+member of ``il``, a ``std::vector<int>``, is a flattened list of all
+interactions for a given interaction type. The vector contains ``N+1`` integer indices
+for each interaction, where ``N`` is the number of particles that are
+involved in the interaction. An additional index is needed to retrieve
+the correct parameters in ``iparams``, hence the total number of indices sums up
+to ``N+1`` per interaction.
+
+The big advantage of storing all types in a union data type is (was),
+that it allows looping over all types with a simple for-loop.
+In pre C++11 and perhaps even pre C++14 times, looping over different
+types was a big hassle and the union data type approach likely was the
+only practicable solution. One downside of this approach, however, is
+that with just a single (union) type, one can't leverage the compiler's
+type system, most importantly static branching, for example with overload resolution.
+As a consequence, only dynamic branching with ``if`` statements remains.
+
+Consider, for instance, the implementation of the top-level
+``calc_listed(const InteractionDefinitions& idef, ...)`` in GROMACS, which in its essence,
+looks like this:
+
+.. code:: cpp
+
+ void calc_listed(const InteractionDefinitions& idef, ...)
+ {
+ // manage timing and multi-threading
+
+ for (int ftype = 0; ftype < F_NRE; ++type)
+ {
+ // branch out and descend stack for 2 intermediate functions based on
+ // the type of interaction that ftype corresponds to
+ // then call a function from a pointer table
+
+ bondFunction* bonded = bondedInteractionFunctions[ftype];
+
+ // compute all forces for ftype
+ bonded(idef.iparams, idef.il[ftype], ...);
+ }
+
+ // reduce thread output
+ }
+
+GROMACS supports a lot of different listed interaction types, such as different
+types of bonds, angles and proper and improper dihedrals. These different types
+require different handling and finally the right force kernel chosen from a table
+of function pointers.
+The handling code required to correctly branch out to all the different cases
+results in quite a deep call stack, a lot of branching logic and ends up accounting
+for a fair part of the overall complexity, which should ideally just consist of
+the type-specific force calculation implementations.
+
+
+A type-aware approach to listed forces
+--------------------------------------
+
+NB-LIB aims to reduce the overall code complexity with a type-aware data format
+where each interaction type is implemented as a separate (C++)-type.
+The format for a given interaction type looks like this:
+
+.. code:: cpp
+
+ template <class Interaction>
+ struct InteractionData
+ {
+ std::vector<Index<Interaction>> indices;
+ std::vector<Interaction> parameters;
+ };
+
+For each type of interaction, we store the interaction indices plus the
+interaction parameters. While the (C++)-types are different, the actual data stored is
+exactly the same: ``N+1`` integer indices per ``N``-center interaction plus the unique parameters.
+An example for ``Interaction`` would be ``HarmonicBond``, the public part of which looks like this:
+
+.. code:: cpp
+
+ class HarmonicBond
+ {
+ public:
+ // return lvalue ref for use with std::tie
+ // in order to leverage std::tuple comparison ops
+ const real& forceConstant();
+ const real& equilDistance();
+ };
+
+The ``Index`` traits class deduces to ``std::array<int, 3>``, because
+for each harmonic bond, we need two ``int``\ s for the coordinate
+indices and a third ``int`` to look up the bond parameters in the
+``parameters`` vector. For angles and dihedrals, the ``Index`` trait
+would add an additional one or two ``int``\ s to hold the additional
+coordinate indices.
+
+Finally, we gather all types of interactions in a
+``std::tuple``, such that the complete definition for listed forces
+in NB-LIB looks like this:
+
+.. code:: cpp
+
+ using ListedInteractions = std::tuple<InteractionData<HarmonicBond>, ..., InteractionData<HarmonicAngle>, ...>;
+
+One important property of ``ListedInteractions`` is that it stores exactly the same information as ``InteractionDefinitions``
+and therefore conversion in either direction is easy to implement.
+
+
+The NB-LIB listed forces pipeline
+---------------------------------
+
+Given the listed interaction data provided in the format described above,
+the steps required to calculate the corresponding forces
+are, in brief:
+
+ * Loop over all interaction types
+ * Loop over all interactions for given type
+ * Call interaction type kernel, store forces and return energy
+
+
+This procedure is identical to the current implementation in GROMACS.
+In actual code, the first step looks like this:
+
+.. code:: cpp
+
+ template<class Buffer, class Pbc>
+ auto reduceListedForces(const ListedInteractions& interactions,
+ const std::vector<gmx::RVec>& x,
+ Buffer* forces,
+ const Pbc& pbc)
+ {
+ std::array<real, std::tuple_size<ListedInteractions>::value> energies;
+
+ // lambda function, will be applied to each type
+ auto computeForceType = [forces, &x, &energies, &pbc](const auto& ielem) {
+ real energy = computeForces(ielem.indices, ielem.parameters, x, forces, pbc);
+ energies[FindIndex<std::decay_t<decltype(ilem)>, ListedInteractions>{}] = energy;
+ };
+
+ // apply the lambda to all bond types
+ for_each_tuple(computeForceType, interactions);
+
+ return energies;
+ }
+
+With the help of a generic lambda and C++17’s ``std::apply`` in the
+one-liner ``for_each_tuple``, we can generate the loop over the
+different types in the tuple quite effortlessly. While
+``reduceListedForces`` implements a loop over the interaction types, the
+next layer, ``computeForces`` implements a loop over all interactions of
+a given type:
+
+.. code:: cpp
+
+ template <class Index, class InteractionType, class Buffer, class Pbc>
+ real computeForces(const std::vector<Index>& indices,
+ const std::vector<InteractionType>& iParams,
+ const std::vector<gmx::RVec>& x,
+ Buffer* forces,
+ const Pbc& pbc)
+ {
+ real Epot = 0.0;
+
+ for (const auto& index : indices)
+ {
+ Epot += dispatchInteraction(index, iParams, x, forces);
+ }
+
+ return Epot;
+ }
+
+Compared to the union data type approach where this loop has been manually
+implemented for all interaction types, in NB-LIB, only a single implementation
+is required.
+
+We’re now down to the level of individual bonds, angles and dihedrals.
+At this point, the next steps depend on the actual type of the
+interaction. But instead of dispatching each harmonic bond, cubic bond,
+harmonic angle and so on to their seperate paths just yet, we just
+differentiate based on the number of interaction centers for now.
+Through overload resolution, the appropriate version
+``dispatchInteraction`` gets called now, such as this one for the case
+of 2-center interactions:
+
+.. code:: cpp
+
+ template <class Buffer, class TwoCenterType, class Pbc>
+ std::enable_if_t<IsTwoCenter<TwoCenterType>::value, real>
+ dispatchInteraction(const InteractionIndex<TwoCenterType>& index,
+ const std::vector<TwoCenterType>& bondInstances,
+ const std::vector<gmx::RVec>& x,
+ Buffer* forces,
+ const Pbc& pbc)
+ {
+ int i = std::get<0>(index);
+ int j = std::get<1>(index);
+ const gmx::RVec& x1 = x[i];
+ const gmx::RVec& x2 = x[j];
+ const TwoCenterType& bond = bondInstances[std::get<2>(index)];
+
+ gmx::RVec dx;
+ // calculate x1 - x2 modulo pbc
+ pbc.dxAiuc(x1, x2, dx);
+ real dr2 = dot(dx, dx);
+ real dr = std::sqrt(dr2);
+
+ auto [force, energy] = bondKernel(dr, bond);
+
+ // avoid division by 0
+ if (dr2 != 0.0)
+ {
+ force /= dr;
+ detail::spreadTwoCenterForces(force, dx, &(*forces)[i], &(*forces)[j]);
+ }
+
+ return energy;
+ }
+
+We can again observe that common parts among different 2-center interaction types
+are reused. The common parts are
+
+ * coordinate retrieval
+ * computation of the scalar distance
+ * spreading of the scalar part of the force to the two centers
+
+The only remaining thing to do now is to call the actual
+kernel to compute the force. Since ``bond`` has a distinct type, we can
+again use overload resolution:
+
+.. code:: cpp
+
+ template <class T>
+ auto bondKernel(T dr, const HarmonicBond& bond)
+ {
+ return harmonicScalarForce(bond.forceConstant(), bond.equilDistance(), dr);
+ }
+
+and call the actual kernel, which in its simplest form for a harmonic
+bond looks like this:
+
+.. code:: cpp
+
+ template <class T>
+ std::tuple<T, T> harmonicScalarForce(T k, T x0, T x)
+ {
+ real dx = x - x0;
+ real dx2 = dx * dx;
+
+ real force = -k * dx;
+ real epot = 0.5 * k * dx2;
+
+ return std::make_tuple(force, epot);
+
+ /* That was 6 flops */
+ }
+
+That’s it! The approach outlined here manages to reuse (between different types)
+a significant part of the code that feeds input data to force kernels.
+Notably, not a single ``if(ftype)`` is required to implement the control flow.
+The remaining parts for a feature complete implementation are
+overloads of ``dispatchInteraction`` for the 3- to 5-center interactions and
+the type-aware wrappers for all the different kernels implemented in
+GROMACS. They have been omitted for brevity.
+
+A note on **multithreading**: multithreading is handled above the top-level
+``reduceListedForces`` described here. For parallel execution, the
+input ``ListedInteractions`` tuple is split into ``nThreads`` parts and a
+``Buffer`` object is set up for each thread. ``reduceListedForces`` is then
+called once by each thread with the assigned fraction of ``ListedInteractions``
+and the ``Buffer`` as argument.
+The lifetime of the ``ListedInteractions`` splits is coupled to the domain decomposition.
+
+Summary
+-------
+
+NB-LIB listed forces employs a (C++)-type aware data format that
+is otherwise equivalent to its counter-part in GROMACS.
+The type-aware data format is then used to simplify the "routing" layer that
+connects data input to the appropriate kernels. Thanks to static branching and polymorphism,
+increased code reuse and simplified branching logic could be achieved.
+**The force kernels themselves do not need to be changed and NB-LIB refers to
+GROMACS for their implementation.**
+
+
+Outlook
+-------
+
+The data flow management for listed forces described here allows further
+improvements to be implemented:
+
+* Aggregate interaction types: fuse interactions of different types into
+ aggregated types. For example, a dihedral interaction and the bonds and angles
+ that are present among the same four particle indices can be combined into a single
+ aggregated interaction. This allows to reuse the particle coordinates loaded from memory
+ for multiple types and also combines the store operations for the forces.
+ Type aggregates also likely simplify an efficient GPU implementation of listed forces.
+
+* Separation of a topology containing both parameter sets for a system state A and B into two
+ separate topologies for the A and B system states.
--- /dev/null
+Adding New Listed-Interaction Types in NB-LIB
+=============================================
+
+NB-LIB currently has code paths for listed interactions that occur between two, three, four and five different particles.
+To extend NB-LIB to support more types of particle interactions, modify the following three files.
+
+Two center interactions must use the distance between the centers as an input to the force kernel.
+Three center interactions take the form ``(particleI, particleJ, ParticleK)``.
+In this case, the middle particle, ``particleJ`` is taken as the center around which the angle is computed.
+This angle must be an input to a three center force kernel.
+Likewise for four center interactions, the dihedral angle phi must be an input to the force kernel.
+Accepting these constraints, it is possible to add a new kernel by modifying the following three files.
+
+1) bondtypes.h_
+2) definitions.h_
+3) kernels.hpp_
+
+.. _bondtypes.h:
+
+1) bondtypes.h
+---------------
+
+This file contains one C++ type to store the parameters for each interaction type.
+New interaction types are added here as separate C++ types.
+The interface of these types is completely unrestricted.
+The only requirements are equality and less than comparison, and that the interface be
+compatible with the corresponding (user-added) kernel.
+
+.. _definitions.h:
+
+2) definitions.h
+------------------------
+
+This file begins with pre-processor macro lists that classify concrete interaction types into two, three, four and five center types.
+To add a new type, the user must add the interaction type parameter struct name to the macro of the correct center number.
+In this case, ``NewBondType`` is an example of a two center interaction.
+As such it would get added to the ``SUPPORTED_TWO_CENTER_TYPES`` macro.
+Assuming that the only other two center interaction is called ``DefaultBond``, the result would look like the following snippet.
+
+.. code:: cpp
+
+ #define SUPPORTED_TWO_CENTER_TYPES DefaultBond, NewBondType
+
+.. _kernels.hpp:
+
+Adding ``NewBondType`` to this macro ensures that the NBLIB ``molecule``
+class ``addInteraction`` function supports adding the new bond type
+and includes it in the listed interaction data that the ``topology`` class
+provides. The ``SUPPORTED_TWO_CENTER_TYPES`` macro is immediately converted into a
+C++ type list that is implemented as a variadic template. The type list
+is then used to define all the dependent data structures. Apart from creating
+the type list, the only place where the macro is needed is explicit template instantiation.
+
+Note that, as of C++17, there's no alternative to preprocessor macros for adding
+the required template instantiations controlled through the macros described here.
+(Other than manually adding the template instantiations, which would require the instantiation list
+of several templates to be updated each time a new interaction type is added. Compared to the preprocessor
+based solution where just a single macro has to be extended, this would clearly be an inferior solution.)
+In NBLIB, the design decision we took, was that we did not want to expose a templated
+interface in a user header and it is for this reason that we explicitly need
+to instantiate the interface with all the supported listed interaction types defined
+in this macro.
+
+3) kernels.hpp
+---------------------
+
+In this file the actual force kernels for each interaction type are implemented.
+Each kernel call is templated to allow various precisions and is
+accessed through an overload ``bondKernel`` that extracts the relevant
+parameters from a ``const NewBondType&`` argument.
+The kernel return type is always an ``std::tuple`` of the force and the potential.
+
+.. code:: cpp
+
+ /*! \brief kernel to calculate the new bond type force
+ *
+ * \param k Force constant
+ * \param x0 Equilibrium distance
+ * \param scale The scaling factor
+ * \param x Input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+ template <class T>
+ std::tuple<T, T> newBondForce(T k, T x0, T scale, T x)
+ {
+ real exponent = std::exp( (x - x0) / scale);
+ real epot = k * exponent;
+ real force = epot / scale;
+ return std::make_tuple(force, epot);
+ }
+
+ template <class T>
+ inline std::tuple<T, T> bondKernel(T dr, const NewBondType& bond)
+ {
+ return newBondForce(bond.forceConstant(), bond.equilDistance(), bond.scaleFactor(), dr);
+ }
+
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff