# \author Sebastian Keller <keller@cscs.ch>
#
+target_sources(nblib
+ PRIVATE
+ calculator.cpp
+ transformations.cpp
+ )
+
if(GMX_INSTALL_NBLIB_API)
install(FILES
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.
+ */
+/*! \internal \file
+ * \brief
+ * Implements a bonded force calculator
+ *
+ * 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>
+ */
+#include "nblib/box.h"
+#include "nblib/exception.h"
+#include "nblib/pbc.hpp"
+#include "nblib/listed_forces/calculator.h"
+#include "nblib/listed_forces/dataflow.hpp"
+#include "nblib/listed_forces/helpers.hpp"
+
+namespace nblib
+{
+
+ListedForceCalculator::~ListedForceCalculator() = default;
+
+ListedForceCalculator::ListedForceCalculator(const ListedInteractionData& interactions,
+ size_t bufferSize,
+ int nthr,
+ const Box& box) :
+ numThreads(nthr),
+ masterForceBuffer_(bufferSize, Vec3{ 0, 0, 0 }),
+ pbcHolder_(std::make_unique<PbcHolder>(box))
+{
+ // split up the work
+ threadedInteractions_ = splitListedWork(interactions, bufferSize, numThreads);
+
+ // set up the reduction buffers
+ int threadRange = bufferSize / numThreads;
+ for (int i = 0; i < numThreads; ++i)
+ {
+ int rangeStart = i * threadRange;
+ int rangeEnd = rangeStart + threadRange;
+ // the last range is a bit bigger due to integer division
+ if (i == numThreads - 1)
+ {
+ rangeEnd = bufferSize;
+ }
+
+ threadedForceBuffers_.push_back(std::make_unique<ForceBuffer<Vec3>>(
+ masterForceBuffer_.data(), rangeStart, rangeEnd));
+ }
+}
+
+void ListedForceCalculator::computeForcesAndEnergies(gmx::ArrayRef<const Vec3> x, bool usePbc)
+{
+ energyBuffer_.fill(0);
+ std::vector<std::array<real, std::tuple_size<ListedInteractionData>::value>> energiesPerThread(numThreads);
+
+#pragma omp parallel for num_threads(numThreads) schedule(static)
+ for (int thread = 0; thread < numThreads; ++thread)
+ {
+ if (usePbc)
+ {
+ energiesPerThread[thread] = reduceListedForces(
+ threadedInteractions_[thread], x, threadedForceBuffers_[thread].get(), *pbcHolder_);
+ }
+ else
+ {
+ energiesPerThread[thread] = reduceListedForces(
+ threadedInteractions_[thread], x, threadedForceBuffers_[thread].get(), NoPbc{});
+ }
+ }
+
+ // reduce energies
+ for (int thread = 0; thread < numThreads; ++thread)
+ {
+ for (int type = 0; type < int(energyBuffer_.size()); ++type)
+ {
+ energyBuffer_[type] += energiesPerThread[thread][type];
+ }
+ }
+
+ // reduce forces
+#pragma omp parallel for num_threads(numThreads) schedule(static)
+ for (int thread = 0; thread < numThreads; ++thread)
+ {
+ auto& thisBuffer = *threadedForceBuffers_[thread];
+ // access outliers from other threads
+ for (int otherThread = 0; otherThread < numThreads; ++otherThread)
+ {
+ auto& otherBuffer = *threadedForceBuffers_[otherThread];
+ for (const auto& outlier : otherBuffer)
+ {
+ int index = outlier.first;
+ // check whether this outlier falls within the range of <thread>
+ if (thisBuffer.inRange(index))
+ {
+ auto force = outlier.second;
+ masterForceBuffer_[index] += force;
+ }
+ }
+ }
+ }
+}
+
+void ListedForceCalculator::compute(gmx::ArrayRef<const Vec3> coordinates, gmx::ArrayRef<Vec3> forces, bool usePbc)
+{
+ if (coordinates.size() != forces.size())
+ {
+ throw InputException("Coordinates array and force buffer size mismatch");
+ }
+
+ // check if the force buffers have the same size
+ if (masterForceBuffer_.size() != forces.size())
+ {
+ throw InputException("Input force buffer size mismatch with listed forces buffer");
+ }
+
+ // compute forces and fill in local buffers
+ computeForcesAndEnergies(coordinates, usePbc);
+
+ // add forces to output force buffers
+ for (int pIndex = 0; pIndex < int(forces.size()); pIndex++)
+ {
+ forces[pIndex] += masterForceBuffer_[pIndex];
+ }
+}
+void ListedForceCalculator::compute(gmx::ArrayRef<const Vec3> coordinates,
+ gmx::ArrayRef<Vec3> forces,
+ EnergyType& energies,
+ bool usePbc)
+{
+ compute(coordinates, forces, usePbc);
+
+ energies = energyBuffer_;
+}
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib 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>
+ */
+
+#ifndef NBLIB_LISTEDFORCES_CONVERSION_HPP
+#define NBLIB_LISTEDFORCES_CONVERSION_HPP
+
+#include "gromacs/topology/forcefieldparameters.h"
+#include "gromacs/topology/idef.h"
+#include "nblib/listed_forces/traits.h"
+
+namespace nblib
+{
+
+/*! this trait maps nblib InteractionTypes to the corresponding gmx enum
+ *
+ * @tparam InteractionType
+ */
+template<class InteractionType>
+struct ListedIndex
+{
+};
+
+template<>
+struct ListedIndex<HarmonicBondType> : std::integral_constant<int, F_BONDS>
+{
+};
+
+template<>
+struct ListedIndex<MorseBondType> : std::integral_constant<int, F_MORSE>
+{
+};
+
+template<>
+struct ListedIndex<FENEBondType> : std::integral_constant<int, F_FENEBONDS>
+{
+};
+
+template<>
+struct ListedIndex<CubicBondType> : std::integral_constant<int, F_CUBICBONDS>
+{
+};
+
+template<>
+struct ListedIndex<G96BondType> : std::integral_constant<int, F_G96BONDS>
+{
+};
+
+template<>
+struct ListedIndex<DefaultAngle> : std::integral_constant<int, F_ANGLES>
+{
+};
+
+template<>
+struct ListedIndex<ProperDihedral> : std::integral_constant<int, F_PDIHS>
+{
+};
+
+template<>
+struct ListedIndex<ImproperDihedral> : std::integral_constant<int, F_IDIHS>
+{
+};
+
+template<>
+struct ListedIndex<RyckaertBellemanDihedral> : std::integral_constant<int, F_RBDIHS>
+{
+};
+
+template<class InteractionType>
+struct ListedTypeIsImplemented : std::bool_constant<detail::HasValueMember<ListedIndex<InteractionType>>::value>
+{
+};
+
+namespace detail
+{
+
+template<class InteractionData>
+void transferParameters([[maybe_unused]] const InteractionData& interactionData,
+ [[maybe_unused]] gmx_ffparams_t& gmx_params)
+{
+}
+
+template<>
+void transferParameters(const TwoCenterData<HarmonicBondType>& interactions, gmx_ffparams_t& gmx_params)
+{
+ for (const auto& hbond : interactions.parameters)
+ {
+ t_iparams param;
+ param.harmonic.krA = hbond.forceConstant();
+ param.harmonic.rA = hbond.equilDistance();
+ gmx_params.iparams.push_back(param);
+ }
+}
+
+template<>
+void transferParameters(const ThreeCenterData<DefaultAngle>& interactions, gmx_ffparams_t& gmx_params)
+{
+ for (const auto& angle : interactions.parameters)
+ {
+ t_iparams param;
+ param.harmonic.krA = angle.forceConstant();
+ param.harmonic.rA = angle.equilDistance() / DEG2RAD;
+ gmx_params.iparams.push_back(param);
+ }
+}
+
+template<>
+void transferParameters(const FourCenterData<ProperDihedral>& interactions, gmx_ffparams_t& gmx_params)
+{
+ for (const auto& dihedral : interactions.parameters)
+ {
+ t_iparams param;
+ param.pdihs.phiA = dihedral.equilDistance() / DEG2RAD;
+ param.pdihs.cpA = dihedral.forceConstant();
+ param.pdihs.mult = dihedral.multiplicity();
+ gmx_params.iparams.push_back(param);
+ }
+}
+
+template<class TwoCenterType>
+void transferIndicesImpl(const TwoCenterData<TwoCenterType>& interactions, InteractionDefinitions& idef, int offset)
+{
+ for (const auto& index : interactions.indices)
+ {
+ int parameterIndex = index[2] + offset;
+ idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(parameterIndex);
+ idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(index[0]);
+ idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(index[1]);
+ }
+}
+
+template<class ThreeCenterType>
+void transferIndicesImpl(const ThreeCenterData<ThreeCenterType>& interactions, InteractionDefinitions& idef, int offset)
+{
+ for (const auto& index : interactions.indices)
+ {
+ int parameterIndex = index[3] + offset;
+ idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(parameterIndex);
+ idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[0]);
+ idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[1]);
+ idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[2]);
+ }
+}
+
+template<class FourCenterType>
+void transferIndicesImpl(const FourCenterData<FourCenterType>& interactions, InteractionDefinitions& idef, int offset)
+{
+ for (const auto& index : interactions.indices)
+ {
+ int parameterIndex = index[4] + offset;
+ idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(parameterIndex);
+ idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[0]);
+ idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[1]);
+ idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[2]);
+ idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[3]);
+ }
+}
+
+template<class FiveCenterType>
+void transferIndicesImpl(const FiveCenterData<FiveCenterType>& interactions, InteractionDefinitions& idef, int offset)
+{
+ for (const auto& index : interactions.indices)
+ {
+ int parameterIndex = index[5] + offset;
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(parameterIndex);
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[0]);
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[1]);
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[2]);
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[3]);
+ idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[4]);
+ }
+}
+
+template<template<class> class Container, class InteractionType>
+void transferIndices(const Container<InteractionType>& interactionData,
+ InteractionDefinitions& idef,
+ [[maybe_unused]] int offset)
+{
+ if constexpr (ListedTypeIsImplemented<InteractionType>{})
+ {
+ transferIndicesImpl(interactionData, idef, offset);
+ }
+}
+
+} // namespace detail
+
+/*! \brief function to convert from nblib-ListedInteractionData to gmx-InteractionDefinitions
+ *
+ * Currently only supports harmonic bonds and angles, other types are ignored
+ *
+ * \param interactions
+ * \return
+ */
+std::tuple<std::unique_ptr<InteractionDefinitions>, std::unique_ptr<gmx_ffparams_t>>
+createFFparams(const ListedInteractionData& interactions);
+
+std::tuple<std::unique_ptr<InteractionDefinitions>, std::unique_ptr<gmx_ffparams_t>>
+createFFparams(const ListedInteractionData& interactions)
+{
+ std::unique_ptr<gmx_ffparams_t> ffparamsHolder = std::make_unique<gmx_ffparams_t>();
+ std::unique_ptr<InteractionDefinitions> idefHolder = std::make_unique<InteractionDefinitions>(*ffparamsHolder);
+
+ gmx_ffparams_t& ffparams = *ffparamsHolder;
+ InteractionDefinitions& idef = *idefHolder;
+
+ auto copyParamsOneType = [&ffparams](const auto& interactionElement)
+ {
+ detail::transferParameters(interactionElement, ffparams);
+ };
+ for_each_tuple(copyParamsOneType, interactions);
+
+ // since gmx_ffparams_t.iparams is a flattened vector over all interaction types,
+ // we need to compute offsets for each type to know where the parameters for each type start
+ // in the flattened iparams vectors
+ int acc = 0;
+ std::array<int, std::tuple_size_v<ListedInteractionData>> indexOffsets{0};
+ auto extractNIndices = [&indexOffsets, &acc](const auto& interactionElement)
+ {
+ constexpr int elementIndex = FindIndex<std::decay_t<decltype(interactionElement)>, ListedInteractionData>::value;
+ indexOffsets[elementIndex] = acc;
+ acc += interactionElement.parameters.size();
+ };
+ for_each_tuple(extractNIndices, interactions);
+
+ auto copyIndicesOneType = [&idef, &indexOffsets](const auto& interactionElement)
+ {
+ constexpr int elementIndex = FindIndex<std::decay_t<decltype(interactionElement)>, ListedInteractionData>::value;
+ detail::transferIndices(interactionElement, idef, indexOffsets[elementIndex]);
+ };
+ for_each_tuple(copyIndicesOneType, interactions);
+
+ return std::make_tuple(std::move(idefHolder), std::move(ffparamsHolder));
+}
+
+
+} // namespace nblib
+
+#endif // NBLIB_LISTEDFORCES_CONVERSION_HPP
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements the data flow from ListedInteractionData and coordinates
+ * down to the individual interaction type kernels
+ *
+ * Intended for internal use inside ListedCalculator only.
+ *
+ * \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_DATAFLOW_HPP
+#define NBLIB_LISTEDFORCES_DATAFLOW_HPP
+
+#include <tuple>
+#include <vector>
+
+#include "nblib/listed_forces/traits.h"
+#include "nblib/listed_forces/kernels.hpp"
+#include "nblib/pbc.hpp"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/arrayref.h"
+
+namespace nblib
+{
+
+template<class TwoCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+
+ ValueType dr2 = dot(dx, dx);
+ ValueType dr = std::sqrt(dr2);
+ auto [force, energy] = bondKernel(dr, parameters);
+
+ // avoid division by 0
+ if (dr2 != 0.0)
+ {
+ force /= dr;
+ spreadTwoCenterForces(force, dx, fi, fj);
+ }
+
+ return energy;
+}
+
+/*! \brief calculate two-center interactions
+ *
+ * \tparam BondType
+ * \param index
+ * \param bondInstances
+ * \param x
+ * \param forces
+ * \return
+ */
+template <class Buffer, class TwoCenterType, class BasicVector, class Pbc>
+inline NBLIB_ALWAYS_INLINE
+auto dispatchInteraction(const TwoCenterInteractionIndex& index,
+ const std::vector<TwoCenterType>& bondInstances,
+ gmx::ArrayRef<const BasicVector> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
+
+ 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);
+
+ energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j]);
+ return energy;
+}
+
+
+template<class ThreeCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+std::enable_if_t<HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
+addTwoCenterAggregate(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
+ BasicVector* fi, BasicVector* fj, BasicVector* fk)
+{
+ if (parameters.manifest == ThreeCenterType::Cargo::ij)
+ {
+ // i-j bond
+ return computeTwoCenter(parameters.twoCenter(), rij, fi, fj);
+ }
+ if (parameters.manifest == ThreeCenterType::Cargo::jk)
+ {
+ // j-k bond
+ return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj);
+ }
+
+ // aggregate is empty
+ return 0.0;
+};
+
+template<class ThreeCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+std::enable_if_t<!HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
+addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters,
+ [[maybe_unused]] const BasicVector& rij,
+ [[maybe_unused]] const BasicVector& rkj,
+ [[maybe_unused]] BasicVector* fi,
+ [[maybe_unused]] BasicVector* fj,
+ [[maybe_unused]] BasicVector* fk)
+{
+ return 0.0;
+};
+
+template<class ThreeCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
+ BasicVector* fi, BasicVector* fj, BasicVector* fk)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+ //! calculate 3-center common quantities: angle between x1-x2 and x2-x3
+ //! Todo: after sufficient evaluation, switch over to atan2 based algorithm
+ ValueType costh = cos_angle(rij, rkj); /* 25 */
+ ValueType theta = std::acos(costh); /* 10 */
+
+ //! call type-specific angle kernel, e.g. harmonic, linear, quartic, etc.
+ auto [force, energy] = threeCenterKernel(theta, parameters);
+
+ spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk);
+
+ return energy;
+}
+
+/*! \brief calculate three-center interactions
+ *
+ * \tparam BondType
+ * \param index
+ * \param bondInstances
+ * \param x
+ * \param forces
+ * \return
+ */
+template <class Buffer, class ThreeCenterType, class BasicVector, class Pbc>
+inline NBLIB_ALWAYS_INLINE
+auto dispatchInteraction(const ThreeCenterInteractionIndex& index,
+ const std::vector<ThreeCenterType>& parameters,
+ gmx::ArrayRef<const BasicVector> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
+
+ //! fetch input data: position vectors x1-x3 and interaction parameters
+ int i = std::get<0>(index);
+ int j = std::get<1>(index);
+ int k = std::get<2>(index);
+ const gmx::RVec& xi = x[i];
+ const gmx::RVec& xj = x[j];
+ const gmx::RVec& xk = x[k];
+ const ThreeCenterType& threeCenterParameters = parameters[std::get<3>(index)];
+
+ gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
+
+ gmx::RVec rij, rkj;
+ pbc.dxAiuc(xi, xj, rij); /* 3 */
+ pbc.dxAiuc(xk, xj, rkj); /* 3 */
+
+ energy.carrier() = computeThreeCenter(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
+ energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
+
+ (*forces)[i] += fi;
+ (*forces)[j] += fj;
+ (*forces)[k] += fk;
+
+ return energy;
+}
+
+template<class FourCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+std::enable_if_t<HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
+addThreeCenterAggregate(const FourCenterType& parameters,
+ const BasicVector& rij, const BasicVector& rkj, const BasicVector& rkl,
+ BasicVector* fi, BasicVector* fj, BasicVector* fk, BasicVector* fl)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+ if (parameters.manifest == FourCenterType::Cargo::j)
+ {
+ return computeThreeCenter(parameters.threeCenter(), rij, rkj, fi, fj, fk);
+ }
+ if (parameters.manifest == FourCenterType::Cargo::k)
+ {
+ return computeThreeCenter(parameters.threeCenter(), ValueType(-1.0)*rkj, ValueType(-1.0)*rkl, fj, fk, fl);
+ //return computeThreeCenter(parameters.threeCenter(), rkj, rkl, fj, fk, fl);
+ }
+
+ // aggregate is empty
+ return 0.0;
+};
+
+template<class FourCenterType, class BasicVector>
+inline NBLIB_ALWAYS_INLINE
+std::enable_if_t<!HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
+addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
+ [[maybe_unused]] const BasicVector& rij,
+ [[maybe_unused]] const BasicVector& rkj,
+ [[maybe_unused]] const BasicVector& rkl,
+ [[maybe_unused]] BasicVector* fi,
+ [[maybe_unused]] BasicVector* fj,
+ [[maybe_unused]] BasicVector* fk,
+ [[maybe_unused]] BasicVector* fl)
+{
+return 0.0;
+};
+
+/*! \brief calculate four-center interactions
+ *
+ * \tparam FourCenterType The bond type to compute; used for type deduction
+ * \param[in] index The atom and parameter indices used for computing the interaction
+ * \param[in] parameters The full type-specific interaction list
+ * \param[in] x The coordinates
+ * \param[in/out] forces The forces
+ * \param[in] pbc Object used for computing distances accounting for PBC's
+ * \return
+ */
+template <class Buffer, class FourCenterType, class BasicVector, class Pbc>
+inline NBLIB_ALWAYS_INLINE
+auto dispatchInteraction(const FourCenterInteractionIndex& index,
+ const std::vector<FourCenterType>& parameters,
+ gmx::ArrayRef<const BasicVector> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
+
+ int i = std::get<0>(index);
+ int j = std::get<1>(index);
+ int k = std::get<2>(index);
+ int l = std::get<3>(index);
+
+ const gmx::RVec& xi = x[i];
+ const gmx::RVec& xj = x[j];
+ const gmx::RVec& xk = x[k];
+ const gmx::RVec& xl = x[l];
+
+ gmx::RVec fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
+
+ gmx::RVec dxIJ, dxKJ, dxKL;
+ pbc.dxAiuc(xi, xj, dxIJ);
+ pbc.dxAiuc(xk, xj, dxKJ);
+ pbc.dxAiuc(xk, xl, dxKL);
+
+ const FourCenterType& fourCenterTypeParams = parameters[std::get<4>(index)];
+
+ rvec m, n;
+ real phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
+
+ auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
+
+ energy.carrier() = kernelEnergy;
+ energy.threeCenterAggregate() = addThreeCenterAggregate(fourCenterTypeParams, dxIJ, dxKJ, dxKL, &fi, &fj, &fk, &fl);
+
+ spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl);
+
+ (*forces)[i] += fi;
+ (*forces)[j] += fj;
+ (*forces)[k] += fk;
+ (*forces)[l] += fl;
+
+ return energy;
+}
+
+/*! \brief calculate five-center interactions
+ *
+ * \tparam BondType
+ * \param index
+ * \param bondInstances
+ * \param x
+ * \param forces
+ * \return
+ */
+template <class Buffer, class FiveCenterType, class BasicVector, class Pbc>
+inline NBLIB_ALWAYS_INLINE
+auto dispatchInteraction(const FiveCenterInteractionIndex& index,
+ const std::vector<FiveCenterType>& parameters,
+ gmx::ArrayRef<const BasicVector> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
+
+ int i = std::get<0>(index);
+ int j = std::get<1>(index);
+ int k = std::get<2>(index);
+ int l = std::get<3>(index);
+ int m = std::get<4>(index);
+
+ const gmx::RVec& xi = x[i];
+ const gmx::RVec& xj = x[j];
+ const gmx::RVec& xk = x[k];
+ const gmx::RVec& xl = x[l];
+ const gmx::RVec& xm = x[m];
+
+ gmx::RVec dxIJ, dxJK, dxKL, dxLM;
+ pbc.dxAiuc(xi, xj, dxIJ);
+ pbc.dxAiuc(xj, xk, dxJK);
+ pbc.dxAiuc(xk, xl, dxKL);
+ pbc.dxAiuc(xl, xm, dxLM);
+
+ const FiveCenterType& fiveCenterTypeParams = parameters[std::get<5>(index)];
+
+ ignore_unused(x, forces, fiveCenterTypeParams);
+ return energy;
+}
+
+
+/*! \brief implement a loop over bonds for a given BondType and Kernel
+ * corresponds to e.g. the "bonds" function at Gromacs:bonded.cpp@450
+ *
+ * \tparam BondType
+ * \tparam Kernel unused for now
+ * \param indices interaction atom pair indices + bond parameter index
+ * \param bondInstances bond parameters
+ * \param x coordinate input
+ * \param kernel unused for now
+ * \return
+ */
+template <class Index, class InteractionType, class Buffer, class Pbc>
+auto computeForces(const std::vector<Index>& indices,
+ const std::vector<InteractionType>& interactionParameters,
+ gmx::ArrayRef<const Vec3> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ KernelEnergy<BasicVectorValueType_t<Vec3>> energy;
+
+ for (const auto& index : indices)
+ {
+ energy += dispatchInteraction(index, interactionParameters, x, forces, pbc);
+ }
+
+ return energy;
+}
+
+/*! \brief implement a loop over bond types and accumulate their force contributions
+ *
+ * \param interactions interaction pairs and bond parameters
+ * \param x coordinate input
+ * \param forces output force buffer
+ */
+template<class Buffer, class Pbc>
+auto reduceListedForces(const ListedInteractionData& interactions,
+ gmx::ArrayRef<const Vec3> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ using ValueType = BasicVectorValueType_t<Vec3>;
+ std::array<ValueType, std::tuple_size<ListedInteractionData>::value> energies{0};
+ energies.fill(0);
+
+ // calculate one bond type
+ auto computeForceType = [forces, &x, &energies, &pbc](const auto& interactionElement) {
+ using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
+
+ KernelEnergy<ValueType> energy = computeForces(interactionElement.indices, interactionElement.parameters, x, forces, pbc);
+
+ energies[CarrierIndex<InteractionType>{}] += energy.carrier();
+ energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();
+ energies[ThreeCenterAggregateIndex<InteractionType>{}] += energy.threeCenterAggregate();
+ // energies[FepIndex{}] += energy.freeEnergyDerivative();
+ };
+
+ // calculate all bond types, returns a tuple with the energies for each type
+ for_each_tuple(computeForceType, interactions);
+
+ return energies;
+}
+
+} // namespace nblib
+
+#endif // NBLIB_LISTEDFORCES_DATAFLOW_HPP
--- /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/exception.h"
+#include "nblib/listed_forces/gmxcalculator.h"
+#include "gromacs/listed_forces/listed_forces.cpp"
+#include "gromacs/listed_forces/manage_threading.h"
+#include "gromacs/timing/wallcycle.h"
+#include "conversions.hpp"
+
+namespace nblib
+{
+
+ListedGmxCalculator::ListedGmxCalculator(const ListedInteractionData& interactions,
+ int nP,
+ int nThr,
+ const Box& box) :
+ numParticles(nP),
+ numThreads(nThr),
+ bondedThreading(numThreads, 1, nullptr),
+ shiftBuffer(SHIFTS), // this is gromacs setup code so here use SHIFTS instead of numShiftVectors
+ forceBuffer(2 * numParticles, gmx::RVec{ 0, 0, 0 }),
+ shiftProxy(gmx::ArrayRefWithPadding<gmx::RVec>(&forceBuffer[0],
+ &forceBuffer[numParticles],
+ &forceBuffer[2 * numParticles]),
+ false,
+ shiftBuffer),
+ virialProxy(forceBuffer, false),
+ forceOutputs(shiftProxy, false, virialProxy),
+ enerd(1, 0),
+ lambdaBuffer(numParticles) // just something big enough
+{
+ std::tie(idef, ffparams) = createFFparams(interactions);
+ idef->ilsort = ilsortNO_FE;
+
+ setup_bonded_threading(&bondedThreading, numParticles, false, *idef);
+
+ wcycle = wallcycle_init(nullptr, 0, &cr);
+ set_pbc(&pbc, PbcType::Xyz, box.legacyMatrix());
+
+ stepWork.computeDhdl = false;
+ stepWork.computeVirial = false;
+ stepWork.computeEnergy = true;
+
+ fr.natoms_force = numParticles;
+}
+
+void ListedGmxCalculator::compute(const std::vector<gmx::RVec>& x,
+ std::vector<gmx::RVec>& forces,
+ ListedForceCalculator::EnergyType& energies)
+{
+ if (forces.size() != x.size() || forces.size() >= forceBuffer.size())
+ {
+ throw InputException("Provided force and/or coordinate buffers inconsistent");
+ }
+
+ const rvec* xdata = &(x[0].as_vec());
+
+ energies.fill(0);
+ std::fill(enerd.term, enerd.term + F_NRE, 0.0);
+
+ calc_listed(wcycle, *idef, &bondedThreading, xdata, &forceOutputs, &fr, &pbc, &enerd, &nrnb,
+ lambdaBuffer.data(), nullptr, nullptr, nullptr, stepWork);
+
+ auto transferEnergy = [&energies, this](auto& interactionElement) {
+ using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
+ if constexpr (ListedTypeIsImplemented<InteractionType>{})
+ {
+ constexpr int index = FindIndex<InteractionType, ListedInteractionData>::value;
+ energies[index] = this->enerd.term[ListedIndex<InteractionType>::value];
+ }
+ };
+ for_each_tuple(transferEnergy, ListedInteractionData{});
+
+ // add forces to output force buffers
+ for (int pIndex = 0; pIndex < forces.size(); pIndex++)
+ {
+ forces[pIndex] += forceBuffer[pIndex];
+ }
+}
+
+const InteractionDefinitions& ListedGmxCalculator::getIdef() const
+{
+ return *idef;
+}
+
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements a fixture for calling calc_listed in gromacs
+ * with nblib interaction data
+ *
+ * \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>
+ */
+
+#ifndef NBLIB_LISTEDFORCES_GMXCALCULATOR_H
+#define NBLIB_LISTEDFORCES_GMXCALCULATOR_H
+
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/listed_forces/listed_forces.h"
+#include "gromacs/listed_forces/listed_internal.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forceoutput.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/simulation_workload.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/topology/forcefieldparameters.h"
+#include "gromacs/topology/idef.h"
+#include "nblib/box.h"
+#include "calculator.h"
+
+namespace nblib
+{
+
+/* an encapsulation class for gmx calc_listed
+ *
+ * Holds all the necessary data to call calc_listed
+ * same ctor signature and behavior as the corresponding nblib
+ * ListedForceCalculator
+ */
+class ListedGmxCalculator
+{
+public:
+ ListedGmxCalculator(const ListedInteractionData& interactions, int nP, int nThr, const Box& box);
+
+ void compute(const std::vector<gmx::RVec>& x,
+ std::vector<gmx::RVec>& forces,
+ ListedForceCalculator::EnergyType& energies);
+
+ [[nodiscard]] const InteractionDefinitions& getIdef() const;
+
+private:
+ int numParticles;
+ int numThreads;
+
+
+ std::unique_ptr<InteractionDefinitions> idef;
+ std::unique_ptr<gmx_ffparams_t> ffparams;
+
+ bonded_threading_t bondedThreading;
+
+ std::vector<gmx::RVec> shiftBuffer;
+ std::vector<gmx::RVec> forceBuffer;
+
+ gmx::ForceWithShiftForces shiftProxy;
+ gmx::ForceWithVirial virialProxy;
+ gmx::ForceOutputs forceOutputs; // yet another proxy
+
+ t_forcerec fr;
+ t_fcdata fcdata; // unused
+
+ t_pbc pbc;
+ gmx_wallcycle_t wcycle;
+ gmx_enerdata_t enerd;
+ gmx::StepWorkload stepWork;
+
+ t_nrnb nrnb;
+ t_commrec cr;
+ std::vector<real> lambdaBuffer;
+};
+
+} // namespace nblib
+
+#endif // NBLIB_LISTEDFORCES_GMXCALCULATOR_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Helper data structures and utility functions for the nblib force calculator.
+ * Intended for internal use.
+ *
+ * \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_LISTEDFORCSES_HELPERS_HPP
+#define NBLIB_LISTEDFORCSES_HELPERS_HPP
+
+#include <unordered_map>
+
+#include "nblib/pbc.hpp"
+#include "definitions.h"
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+template<class T>
+inline void gmxRVecZeroWorkaround([[maybe_unused]] T& value)
+{
+}
+
+template<>
+inline void gmxRVecZeroWorkaround<gmx::RVec>(gmx::RVec& value)
+{
+ for (int i = 0; i < dimSize; ++i)
+ {
+ value[i] = 0;
+ }
+}
+} // namespace detail
+
+/*! \internal \brief object to store forces for multithreaded listed forces computation
+ *
+ */
+template<class T>
+class ForceBuffer
+{
+ using HashMap = std::unordered_map<int, T>;
+
+public:
+ ForceBuffer() : rangeStart(0), rangeEnd(0) { }
+
+ ForceBuffer(T* mbuf, int rs, int re) :
+ masterForceBuffer(mbuf),
+ rangeStart(rs),
+ rangeEnd(re)
+ {
+ }
+
+ void clear() { outliers.clear(); }
+
+ inline NBLIB_ALWAYS_INLINE T& operator[](int i)
+ {
+ if (i >= rangeStart && i < rangeEnd)
+ {
+ return masterForceBuffer[i];
+ }
+ else
+ {
+ if (outliers.count(i) == 0)
+ {
+ T zero = T();
+ // if T = gmx::RVec, need to explicitly initialize it to zeros
+ detail::gmxRVecZeroWorkaround(zero);
+ outliers[i] = zero;
+ }
+ return outliers[i];
+ }
+ }
+
+ typename HashMap::const_iterator begin() { return outliers.begin(); }
+ typename HashMap::const_iterator end() { return outliers.end(); }
+
+ [[nodiscard]] bool inRange(int index) const { return (index >= rangeStart && index < rangeEnd); }
+
+private:
+ T* masterForceBuffer;
+ int rangeStart;
+ int rangeEnd;
+
+ HashMap outliers;
+};
+
+namespace detail
+{
+
+static int computeChunkIndex(int index, int totalRange, int nSplits)
+{
+ if (totalRange < nSplits)
+ {
+ // if there's more threads than particles
+ return index;
+ }
+
+ int splitLength = totalRange / nSplits;
+ return std::min(index / splitLength, nSplits - 1);
+}
+
+} // namespace detail
+
+
+/*! \internal \brief splits an interaction tuple into nSplits interaction tuples
+ *
+ * \param interactions
+ * \param totalRange the number of particle sequence coordinates
+ * \param nSplits number to divide the total work by
+ * \return
+ */
+inline
+std::vector<ListedInteractionData> splitListedWork(const ListedInteractionData& interactions,
+ int totalRange,
+ int nSplits)
+{
+ std::vector<ListedInteractionData> workDivision(nSplits);
+
+ auto splitOneElement = [totalRange, nSplits, &workDivision](const auto& inputElement) {
+ // the index of inputElement in the ListedInteractionsTuple
+ constexpr int elementIndex =
+ FindIndex<std::decay_t<decltype(inputElement)>, ListedInteractionData>{};
+
+ // for now, copy all parameters to each split
+ // Todo: extract only the parameters needed for this split
+ for (auto& workDivisionSplit : workDivision)
+ {
+ std::get<elementIndex>(workDivisionSplit).parameters = inputElement.parameters;
+ }
+
+ // loop over all interactions in inputElement
+ for (const auto& interactionIndex : inputElement.indices)
+ {
+ // each interaction has multiple coordinate indices
+ // we must pick one of them to assign this interaction to one of the output index ranges
+ // Todo: count indices outside the current split range in order to minimize the buffer size
+ int representativeIndex =
+ *std::min_element(begin(interactionIndex), end(interactionIndex) - 1);
+ int splitIndex = detail::computeChunkIndex(representativeIndex, totalRange, nSplits);
+
+ std::get<elementIndex>(workDivision[splitIndex]).indices.push_back(interactionIndex);
+ }
+ };
+
+ // split each interaction type in the input interaction tuple
+ for_each_tuple(splitOneElement, interactions);
+
+ return workDivision;
+}
+
+} // namespace nblib
+
+#endif // NBLIB_LISTEDFORCSES_HELPERS_HPP
--- /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 kernels for nblib supported bondtypes
+ *
+ * \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_KERNELS_HPP
+#define NBLIB_LISTEDFORCES_KERNELS_HPP
+
+#include <tuple>
+
+#include "gromacs/math/vec.h"
+#include "nblib/basicdefinitions.h"
+#include "bondtypes.h"
+
+namespace nblib
+{
+
+/*! \brief kernel to calculate the scalar part for simple harmonic bond forces
+ * for lambda = 0
+ *
+ * \param k spring constant
+ * \param x0 equilibrium distance
+ * \param x input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> harmonicScalarForce(T k, T x0, T x)
+{
+ T dx = x - x0;
+ T dx2 = dx * dx;
+
+ T force = -k * dx;
+ T epot = 0.5 * k * dx2;
+
+ return std::make_tuple(force, epot);
+
+ /* That was 6 flops */
+}
+
+/*! \brief kernel to calculate the scalar part for simple harmonic bond forces
+ * for non-zero lambda to interpolate between A and B states
+ *
+ * \param kA spring constant state A
+ * \param kB spring constant state B
+ * \param xA equilibrium distance state A
+ * \param xB equilibrium distance state B
+ * \param x input bond length
+ * \param lambda interpolation factor between A and B state
+ *
+ * \return tuple<force, potential energy, lambda-interpolated energy>
+ */
+template <class T>
+inline std::tuple<T, T, T> harmonicScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
+{
+ // code unchanged relative to Gromacs
+
+ T L1 = 1.0 - lambda;
+ T kk = L1 * kA + lambda * kB;
+ T x0 = L1 * xA + lambda * xB;
+
+ T dx = x - x0;
+ T dx2 = dx * dx;
+
+ T force = -kk * dx;
+ T epot = 0.5 * kk * dx2;
+ T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
+
+ return std::make_tuple(force, epot, dvdlambda);
+
+ /* That was 19 flops */
+}
+
+//! abstraction layer for different 2-center bonds
+template <class T>
+inline auto bondKernel(T dr, const HarmonicBondType& bond)
+{
+ return harmonicScalarForce(bond.forceConstant(), bond.equilDistance(), dr);
+}
+
+
+/*! \brief kernel to calculate the scalar part for the forth power pontential bond forces
+ * for lambda = 0
+ *
+ * \param k spring constant
+ * \param x0 squared equilibrium distance
+ * \param x squared input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> g96ScalarForce(T k, T x0, T x)
+{
+ T dx = x - x0;
+ T dx2 = dx * dx;
+
+ T force = -k * dx * x;
+ T epot = 0.25 * k * dx2;
+
+ return std::make_tuple(force, epot);
+
+ /* That was 7 flops */
+}
+
+/*! \brief kernel to calculate the scalar part for forth power pontential bond forces
+ * for non-zero lambda to interpolate between A and B states
+ *
+ * \param kA spring constant state A
+ * \param kB spring constant state B
+ * \param xA squared equilibrium distance state A
+ * \param xB squared equilibrium distance state B
+ * \param x squared input bond length
+ * \param lambda interpolation factor between A and B state
+ *
+ * \return tuple<force, potential energy, lambda-interpolated energy>
+ */
+template <class T>
+inline std::tuple<T, T, T> g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
+{
+ T L1 = 1.0 - lambda;
+ T kk = L1 * kA + lambda * kB;
+ T x0 = L1 * xA + lambda * xB;
+
+ T dx = x - x0;
+ T dx2 = dx * dx;
+
+ T force = -kk * dx * x;
+ T epot = 0.25 * kk * dx2;
+ // TODO: Check if this is 1/2 or 1/4
+ T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
+
+ return std::make_tuple(force, epot, dvdlambda);
+
+ /* That was 21 flops */
+}
+
+//! Abstraction layer for different 2-center bonds. Forth power case
+template <class T>
+inline auto bondKernel(T dr, const G96BondType& bond)
+{
+ // NOTE: Not assuming GROMACS' convention of storing squared bond.equilDistance() for this type
+ return g96ScalarForce(bond.forceConstant(), bond.equilDistance() * bond.equilDistance(), dr * dr);
+}
+
+
+/*! \brief kernel to calculate the scalar part for the morse pontential bond forces
+ * for lambda = 0
+ *
+ * \param k force constant
+ * \param beta beta exponent
+ * \param x0 equilibrium distance
+ * \param x input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> morseScalarForce(T k, T beta, T x0, T x)
+{
+ T exponent = std::exp(-beta * (x - x0)); /* 12 */
+ T omexp = 1.0 - exponent; /* 1 */
+ T kexp = k * omexp; /* 1 */
+
+ T epot = kexp * omexp; /* 1 */
+ T force = -2.0 * beta * exponent * omexp; /* 4 */
+
+ return std::make_tuple(force, epot);
+
+ /* That was 23 flops */
+}
+
+/*! \brief kernel to calculate the scalar part for morse potential bond forces
+ * for non-zero lambda to interpolate between A and B states
+ *
+ * \param kA force constant state A
+ * \param kB force constant state B
+ * \param betaA beta exponent state A
+ * \param betaB beta exponent state B
+ * \param xA equilibrium distance state A
+ * \param xB equilibrium distance state B
+ * \param x input bond length
+ * \param lambda interpolation factor between A and B state
+ *
+ * \return tuple<force, potential energy, lambda-interpolated energy>
+ */
+template <class T>
+inline std::tuple<T, T, T> morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, T xB, T x, T lambda)
+{
+ T L1 = 1.0 - lambda; /* 1 */
+ T x0 = L1 * xA + lambda * xB; /* 3 */
+ T beta = L1 * betaA + lambda * betaB; /* 3 */
+ T k = L1 * kA + lambda * kB; /* 3 */
+
+ T exponent = std::exp(-beta * (x - x0)); /* 12 */
+ T omexp = 1.0 - exponent; /* 1 */
+ T kexp = k * omexp; /* 1 */
+
+ T epot = kexp * omexp; /* 1 */
+ T force = -2.0 * beta * exponent * omexp; /* 4 */
+
+ T dvdlambda = (kB - kA) * omexp * omexp
+ - (2.0 - 2.0 * omexp) * omexp * k
+ * ((xB - xA) * beta - (betaB - betaA) * (x - x0)); /* 15 */
+
+ return std::make_tuple(force, epot, dvdlambda);
+
+ /* That was 44 flops */
+}
+
+//! Abstraction layer for different 2-center bonds. Morse case
+template <class T>
+inline auto bondKernel(T dr, const MorseBondType& bond)
+{
+ return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
+}
+
+
+/*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
+ * for lambda = 0
+ *
+ * \param k spring constant
+ * \param x0 equilibrium distance
+ * \param x input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> FENEScalarForce(T k, T x0, T x)
+{
+ T x02 = x0 * x0;
+ T x2 = x * x;
+
+ T omx2_ox02 = 1.0 - (x2 / x02);
+
+ T epot = -0.5 * k * x02 * std::log(omx2_ox02);
+ T force = -k * x / omx2_ox02;
+
+ return std::make_tuple(force, epot);
+
+ /* That was 24 flops */
+}
+
+// TODO: Implement the free energy version of FENE (finitely extensible nonlinear elastic) bond types
+
+//! Abstraction layer for different 2-center bonds. FENE case
+template <class T>
+inline auto bondKernel(T dr, const FENEBondType& bond)
+{
+ return FENEScalarForce(bond.forceConstant(), bond.equilDistance(), dr);
+}
+
+
+/*! \brief kernel to calculate the scalar part for cubic potential bond forces
+ * for lambda = 0
+ *
+ * \param kc cubic spring constant
+ * \param kq quadratic spring constant
+ * \param x0 equilibrium distance
+ * \param x input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> cubicScalarForce(T kc, T kq, T x0, T x)
+{
+ T dx = x - x0;
+
+ T kdist = kq * dx;
+ T kdist2 = kdist * dx;
+
+ T epot = kdist2 + (kc * kdist2 * dx);
+ T force = -((2.0 * kdist) + (3.0 * kdist2 * kc));
+
+ return std::make_tuple(force, epot);
+
+ /* That was 16 flops */
+}
+
+// TODO: Implement the free energy version of Cubic bond types
+
+template <class T>
+inline auto bondKernel(T dr, const CubicBondType& bond)
+{
+ return cubicScalarForce(bond.cubicForceConstant(), bond.quadraticForceConstant(), bond.equilDistance(), dr);
+}
+
+
+/*! \brief kernel to calculate the scalar part for half attractive potential bond forces
+ * for lambda = 0
+ *
+ * \param k spring constant
+ * \param x0 equilibrium distance
+ * \param x input bond length
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> halfAttractiveScalarForce(T k, T x0, T x)
+{
+ T dx = x - x0;
+ T dx2 = dx * dx;
+ T dx3 = dx2 * dx;
+ T dx4 = dx2 * dx2;
+
+ T epot = -0.5 * k * dx4;
+ T force = -2.0 * k * dx3;
+
+ return std::make_tuple(force, epot);
+
+ /* That was 10 flops */
+}
+
+/*! \brief kernel to calculate the scalar part for half attractive potential bond forces
+ * for non-zero lambda to interpolate between A and B states
+ *
+ * \param kA spring constant state A
+ * \param kB spring constant state B
+ * \param xA equilibrium distance state A
+ * \param xB equilibrium distance state B
+ * \param x input bond length
+ * \param lambda interpolation factor between A and B state
+ *
+ * \return tuple<force, potential energy, lambda-interpolated energy>
+ */
+template <class T>
+inline std::tuple<T, T, T> halfAttractiveScalarForce(T kA, T kB, T xA, T xB, T x, T lambda)
+{
+ T L1 = 1.0 - lambda;
+ T kk = L1 * kA + lambda * kB;
+ T x0 = L1 * xA + lambda * xB;
+
+ T dx = x - x0;
+ T dx2 = dx * dx;
+ T dx3 = dx2 * dx;
+ T dx4 = dx2 * dx2;
+
+ T epot = -0.5 * kk * dx4;
+ T force = -2.0 * kk * dx3;
+ T dvdlambda = 0.5 * (kB - kA) * dx4 + (2.0 * (xA - xB) * kk * dx3);
+
+ return std::make_tuple(force, epot, dvdlambda);
+
+ /* That was 29 flops */
+}
+
+template <class T>
+inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond)
+{
+ return halfAttractiveScalarForce(bond.forceConstant(), bond.equilDistance(), dr);
+}
+
+
+//! Three-center interaction type kernels
+
+// linear angles go here
+// quartic angles go here
+
+//! Three-center interaction type dispatch
+
+template <class T>
+inline auto threeCenterKernel(T dr, const DefaultAngle& angle)
+{
+ return harmonicScalarForce(angle.forceConstant(), angle.equilDistance(), dr);
+}
+
+
+//! \brief Computes and returns the proper dihedral force
+template <class T>
+inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
+{
+ const T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
+ const T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
+ const T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
+ return std::make_tuple(force, ePot);
+}
+
+
+//! \brief Ensure that a geometric quantity lies in (-pi, pi)
+static inline void makeAnglePeriodic(real& angle)
+{
+ if (angle >= M_PI)
+ {
+ angle -= 2 * M_PI;
+ }
+ else if (angle < -M_PI)
+ {
+ angle += 2 * M_PI;
+ }
+}
+
+//! \brief Computes and returns a dihedral phi angle
+static inline real dihedralPhi(rvec dxIJ, rvec dxKJ, rvec dxKL, rvec m, rvec n)
+{
+ cprod(dxIJ, dxKJ, m);
+ cprod(dxKJ, dxKL, n);
+ real phi = gmx_angle(m, n);
+ real ipr = iprod(dxIJ, n);
+ real sign = (ipr < 0.0) ? -1.0 : 1.0;
+ phi = sign * phi;
+ return phi;
+}
+
+//! \brief Computes and returns the improper dihedral force
+template <class T>
+inline auto fourCenterKernel(T phi, const ImproperDihedral& improperDihedral)
+{
+ T deltaPhi = phi - improperDihedral.equilDistance();
+ /* deltaPhi cannot be outside (-pi,pi) */
+ makeAnglePeriodic(deltaPhi);
+ const T force = -improperDihedral.forceConstant() * deltaPhi;
+ const T ePot = 0.5 * improperDihedral.forceConstant() * deltaPhi * deltaPhi;
+ return std::make_tuple(force, ePot);
+}
+
+//! \brief Computes and returns the Ryckaert-Belleman dihedral force
+template <class T>
+inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBellemanDihedral)
+{
+ /* Change to polymer convention */
+ const T localPhi = (phi < 0) ? (phi += M_PI) : (phi -= M_PI);
+ T cos_phi = std::cos(localPhi);
+ T ePot = ryckaertBellemanDihedral[0];
+ T force = 0;
+ T cosineFactor = 1;
+
+ for (int i = 1; i < int(ryckaertBellemanDihedral.size()); i++)
+ {
+ force += ryckaertBellemanDihedral[i] * cosineFactor * i;
+ cosineFactor *= cos_phi;
+ ePot += cosineFactor * ryckaertBellemanDihedral[i];
+ }
+ /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
+ force = -force * std::sin(localPhi);
+ return std::make_tuple(force, ePot);
+}
+
+//! Two-center category common
+
+/*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
+ *
+ * \p shiftIndex is used as the periodic shift.
+ */
+template <class T>
+inline void spreadTwoCenterForces(const T bondForce,
+ const gmx::RVec& dx,
+ gmx::RVec* force_i,
+ gmx::RVec* force_j)
+{
+ for (int m = 0; m < dimSize; m++) /* 15 */
+ {
+ const T fij = bondForce * dx[m];
+ (*force_i)[m] += fij;
+ (*force_j)[m] -= fij;
+ }
+}
+
+//! Three-center category common
+
+/*! \brief spread force to 3 centers based on scalar force and angle
+ *
+ * @tparam T
+ * @param cos_theta
+ * @param force
+ * @param r_ij
+ * @param r_kj
+ * @param force_i
+ * @param force_j
+ * @param force_k
+ */
+template <class T>
+inline void spreadThreeCenterForces(T cos_theta,
+ T force,
+ const gmx::RVec& r_ij,
+ const gmx::RVec& r_kj,
+ gmx::RVec* force_i,
+ gmx::RVec* force_j,
+ gmx::RVec* force_k)
+{
+ T cos_theta2 = cos_theta * cos_theta;
+ if (cos_theta2 < 1)
+ {
+ T st = force / std::sqrt(1 - cos_theta2); /* 12 */
+ T sth = st * cos_theta; /* 1 */
+ T nrij2 = dot(r_ij, r_ij); /* 5 */
+ T nrkj2 = dot(r_kj, r_kj); /* 5 */
+
+ T nrij_1 = 1.0 / std::sqrt(nrij2); /* 10 */
+ T nrkj_1 = 1.0 / std::sqrt(nrkj2); /* 10 */
+
+ T cik = st * nrij_1 * nrkj_1; /* 2 */
+ T cii = sth * nrij_1 * nrij_1; /* 2 */
+ T ckk = sth * nrkj_1 * nrkj_1; /* 2 */
+
+ gmx::RVec f_i{0, 0, 0};
+ gmx::RVec f_j{0, 0, 0};
+ gmx::RVec f_k{0, 0, 0};
+ for (int m = 0; m < dimSize; m++)
+ { /* 39 */
+ f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
+ f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
+ f_j[m] = -f_i[m] - f_k[m];
+ (*force_i)[m] += f_i[m];
+ (*force_j)[m] += f_j[m];
+ (*force_k)[m] += f_k[m];
+ }
+ } /* 161 TOTAL */
+}
+
+//! Four-center category common
+template <class T>
+inline void spreadFourCenterForces(T force, rvec dxIJ, rvec dxJK, rvec dxKL, rvec m, rvec n,
+ gmx::RVec* force_i,
+ gmx::RVec* force_j,
+ gmx::RVec* force_k,
+ gmx::RVec* force_l)
+{
+ rvec f_i, f_j, f_k, f_l;
+ rvec uvec, vvec, svec;
+ T iprm = iprod(m, m); /* 5 */
+ T iprn = iprod(n, n); /* 5 */
+ T nrkj2 = iprod(dxJK, dxJK); /* 5 */
+ T toler = nrkj2 * GMX_REAL_EPS;
+ if ((iprm > toler) && (iprn > toler))
+ {
+ T nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
+ T nrkj_2 = nrkj_1 * nrkj_1; /* 1 */
+ T nrkj = nrkj2 * nrkj_1; /* 1 */
+ T a = -force * nrkj / iprm; /* 11 */
+ svmul(a, m, f_i); /* 3 */
+ T b = force * nrkj / iprn; /* 11 */
+ svmul(b, n, f_l); /* 3 */
+ T p = iprod(dxIJ, dxJK); /* 5 */
+ p *= nrkj_2; /* 1 */
+ T q = iprod(dxKL, dxJK); /* 5 */
+ q *= nrkj_2; /* 1 */
+ svmul(p, f_i, uvec); /* 3 */
+ svmul(q, f_l, vvec); /* 3 */
+ rvec_sub(uvec, vvec, svec); /* 3 */
+ rvec_sub(f_i, svec, f_j); /* 3 */
+ rvec_add(f_l, svec, f_k); /* 3 */
+ rvec_inc(force_i->as_vec(), f_i); /* 3 */
+ rvec_dec(force_j->as_vec(), f_j); /* 3 */
+ rvec_dec(force_k->as_vec(), f_k); /* 3 */
+ rvec_inc(force_l->as_vec(), f_l); /* 3 */
+ }
+}
+
+} // namespace nblib
+#endif // NBLIB_LISTEDFORCES_KERNELS_HPP
CPP_SOURCE_FILES
# files with code for tests
bondtypes.cpp
+ helpers.cpp
+ kernels.cpp
+ typetests.cpp
+ calculator.cpp
+ conversions.cpp
+ transformations.cpp
)
target_link_libraries(${exename} PRIVATE nblib_test_infrastructure nblib)
gmx_register_gtest_test(${testname} ${exename} INTEGRATION_TEST)
--- /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 "gmxpre.h"
+
+// includesorter/check-source want this include here. IMO a bug
+
+
+#include <valarray>
+
+#include <gtest/gtest.h>
+
+#include "nblib/listed_forces/calculator.h"
+#include "nblib/listed_forces/dataflow.hpp"
+#include "nblib/listed_forces/tests/linear_chain_input.hpp"
+#include "nblib/tests/testhelpers.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+TEST(NBlibTest, ListedForceCalculatorCanConstruct)
+{
+ ListedInteractionData interactions;
+ Box box(1, 1, 1);
+ EXPECT_NO_THROW(ListedForceCalculator listedForceCalculator(interactions, 2, 1, box));
+}
+
+template<class TestSeq, class SeqFloat, class SeqDouble>
+void compareVectors(const TestSeq& forces,
+ [[maybe_unused]] const SeqFloat& refForcesFloat,
+ [[maybe_unused]] const SeqDouble& refForcesDouble)
+{
+ for (size_t i = 0; i < forces.size(); ++i)
+ {
+ for (int m = 0; m < dimSize; ++m)
+ {
+ EXPECT_FLOAT_DOUBLE_EQ_TOL(
+ forces[i][m], refForcesFloat[i][m], refForcesDouble[i][m],
+ // Todo: why does the tolerance need to be so low?
+ gmx::test::relativeToleranceAsFloatingPoint(refForcesDouble[i][m], 5e-5));
+ }
+ }
+}
+
+class ListedExampleData : public ::testing::Test
+{
+protected:
+ void SetUp() override
+ {
+ // methanol-spc data
+ HarmonicBondType bond1{ 376560, 0.136 };
+ HarmonicBondType bond2{ 313800, 0.1 };
+ std::vector<HarmonicBondType> bonds{ bond1, bond2 };
+ // one bond between atoms 0-1 with bond1 parameters and another between atoms 1-2 with bond2 parameters
+ std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 1, 0 }, { 1, 2, 1 } };
+
+ DefaultAngle angle(Degrees(108.53), 397.5);
+ std::vector<DefaultAngle> angles{ angle };
+ std::vector<InteractionIndex<DefaultAngle>> angleIndices{ { 0, 1, 2, 0 } };
+
+ pickType<HarmonicBondType>(interactions).indices = bondIndices;
+ pickType<HarmonicBondType>(interactions).parameters = bonds;
+
+ pickType<DefaultAngle>(interactions).indices = angleIndices;
+ pickType<DefaultAngle>(interactions).parameters = angles;
+
+ // initial position for the methanol atoms from the spc-water example
+ x = std::vector<gmx::RVec>{ { 1.97, 1.46, 1.209 }, { 1.978, 1.415, 1.082 }, { 1.905, 1.46, 1.03 } };
+ forces = std::vector<gmx::RVec>(3, gmx::RVec{ 0, 0, 0 });
+
+ refBondForcesFloat =
+ std::valarray<gmx::BasicVector<float>>{ { -22.8980637, 128.801575, 363.505951 },
+ { -43.2698593, -88.0130997, -410.639252 },
+ { 66.167923, -40.788475, 47.1333084 } };
+ refAngleForcesFloat =
+ std::valarray<gmx::BasicVector<float>>{ { 54.7276611, -40.1688995, 17.6805191 },
+ { -81.8118973, 86.1988525, 60.1752243 },
+ { 27.0842342, -46.0299492, -77.8557434 } };
+
+ refBondForcesDouble = std::valarray<gmx::BasicVector<double>>{
+ { -22.89764839974935, 128.79927224858977, 363.50016834602064 },
+ { -43.24622441913251, -88.025652017772231, -410.61635172385434 },
+ { 66.14387281888186, -40.773620230817542, 47.116183377833721 }
+ };
+ refAngleForcesDouble = std::valarray<gmx::BasicVector<double>>{
+ { 54.726206806506234, -40.167809526198099, 17.680008528590257 },
+ { -81.809781666748606, 86.196545126117257, 60.173723525141448 },
+ { 27.083574860242372, -46.028735599919159, -77.853732053731704 }
+ };
+
+ refBondEnergyFloat = 0.2113433;
+ refAngleEnergyFloat = 0.112774156;
+
+ refBondEnergyDouble = 0.2113273434867636;
+ refAngleEnergyDouble = 0.11276812148357591;
+
+ box.reset(new Box(3, 3, 3));
+ pbc.reset(new PbcHolder(*box));
+ }
+
+ std::vector<gmx::RVec> x;
+ std::vector<gmx::RVec> forces;
+
+ ListedInteractionData interactions;
+
+ std::shared_ptr<Box> box;
+ std::shared_ptr<PbcHolder> pbc;
+
+ // reference values
+ std::valarray<gmx::BasicVector<float>> refBondForcesFloat, refAngleForcesFloat;
+ std::valarray<gmx::BasicVector<double>> refBondForcesDouble, refAngleForcesDouble;
+
+ float refBondEnergyFloat, refAngleEnergyFloat;
+ double refBondEnergyDouble, refAngleEnergyDouble;
+};
+
+TEST_F(ListedExampleData, ComputeHarmonicBondForces)
+{
+ auto indices = pickType<HarmonicBondType>(interactions).indices;
+ auto bonds = pickType<HarmonicBondType>(interactions).parameters;
+ real energy = computeForces(indices, bonds, x, &forces, *pbc);
+
+ EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refBondEnergyFloat, refBondEnergyDouble,
+ gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
+
+ compareVectors(forces, refBondForcesFloat, refBondForcesDouble);
+}
+
+TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
+{
+ auto indices = pickType<DefaultAngle>(interactions).indices;
+ auto angles = pickType<DefaultAngle>(interactions).parameters;
+ real energy = computeForces(indices, angles, x, &forces, *pbc);
+
+ EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refAngleEnergyFloat, refAngleEnergyDouble,
+ gmx::test::relativeToleranceAsFloatingPoint(refAngleEnergyDouble, 1e-5));
+
+ compareVectors(forces, refAngleForcesFloat, refAngleForcesDouble);
+}
+
+TEST_F(ListedExampleData, CanReduceForces)
+{
+ auto energies = reduceListedForces(interactions, x, &forces, *pbc);
+ real totalEnergy = std::accumulate(begin(energies), end(energies), 0.0);
+
+ EXPECT_FLOAT_DOUBLE_EQ_TOL(totalEnergy, refBondEnergyFloat + refAngleEnergyFloat,
+ refBondEnergyDouble + refAngleEnergyDouble,
+ gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
+
+ compareVectors(forces, refBondForcesFloat + refAngleForcesFloat,
+ refBondForcesDouble + refAngleForcesDouble);
+}
+
+
+void compareArray(const ListedForceCalculator::EnergyType& energies,
+ const ListedForceCalculator::EnergyType& refEnergies)
+{
+ for (size_t i = 0; i < energies.size(); ++i)
+ {
+ EXPECT_REAL_EQ_TOL(energies[i], refEnergies[i],
+ gmx::test::relativeToleranceAsFloatingPoint(refEnergies[i], 1e-5));
+ }
+}
+
+
+//! \brief sets up an interaction tuple for a linear chain with nParticles
+class LinearChainDataFixture : public ::testing::Test
+{
+protected:
+ void SetUp() override
+ {
+ LinearChainData data(430);
+
+ x = data.x;
+ interactions = data.interactions;
+ box = data.box;
+ refForces = data.forces;
+ // pbc.reset(new PbcHolder(*box));
+
+ refEnergies = reduceListedForces(interactions, x, &refForces, NoPbc{});
+ }
+
+ void testEnergies(const ListedForceCalculator::EnergyType& energies) const
+ {
+ compareArray(energies, refEnergies);
+ }
+
+ void testForces(const std::vector<gmx::RVec>& forces) const
+ {
+ compareVectors(forces, refForces, refForces);
+ }
+
+ std::vector<gmx::RVec> x;
+ ListedInteractionData interactions;
+ std::shared_ptr<Box> box;
+ // std::shared_ptr<PbcHolder> pbc;
+
+private:
+ std::vector<gmx::RVec> refForces;
+ ListedForceCalculator::EnergyType refEnergies;
+};
+
+TEST_F(LinearChainDataFixture, Multithreading)
+{
+ ListedForceCalculator lfCalculator(interactions, x.size(), 4, *box);
+
+ std::vector<Vec3> forces(x.size(), Vec3{ 0, 0, 0 });
+ ListedForceCalculator::EnergyType energies;
+ lfCalculator.compute(x, forces, energies);
+
+ testEnergies(energies);
+ testForces(forces);
+}
+
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib 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 <numeric>
+
+#include <gtest/gtest.h>
+
+#include "nblib/listed_forces/conversions.hpp"
+
+#include "testutils/testasserts.h"
+
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+ListedInteractionData someBondsAndAngles()
+{
+ ListedInteractionData interactions;
+ HarmonicBondType bond1{ 10, 0.1 };
+ HarmonicBondType bond2{ 20, 0.2 };
+ std::vector<HarmonicBondType> bonds{ bond1, bond2 };
+ pickType<HarmonicBondType>(interactions).parameters = bonds;
+
+ DefaultAngle angle1(Degrees(100), 100);
+ DefaultAngle angle2(Degrees(101), 200);
+ std::vector<DefaultAngle> angles{ angle1, angle2 };
+ pickType<DefaultAngle>(interactions).parameters = angles;
+
+ std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 1, 0 }, { 1, 2, 0 }, { 2, 3, 1 } };
+ pickType<HarmonicBondType>(interactions).indices = std::move(bondIndices);
+
+ std::vector<InteractionIndex<DefaultAngle>> angleIndices{ { 0, 1, 2, 0 }, { 1, 2, 3, 1 } };
+ pickType<DefaultAngle>(interactions).indices = std::move(angleIndices);
+
+ return interactions;
+}
+
+TEST(ListedShims, ParameterConversion)
+{
+ ListedInteractionData interactions = someBondsAndAngles();
+
+ auto [idef, gmx_params] = createFFparams(interactions);
+
+ EXPECT_EQ(gmx_params->iparams.size(), 4);
+ EXPECT_EQ(gmx_params->iparams[0].harmonic.rA,
+ pickType<HarmonicBondType>(interactions).parameters[0].equilDistance());
+ EXPECT_REAL_EQ_TOL(gmx_params->iparams[2].harmonic.rA,
+ pickType<DefaultAngle>(interactions).parameters[0].equilDistance() / DEG2RAD,
+ gmx::test::defaultRealTolerance());
+
+ EXPECT_EQ(idef->il[F_BONDS].iatoms.size(), 9);
+ std::vector<int> bondIatoms{ 0, 0, 1, 0, 1, 2, 1, 2, 3 };
+ EXPECT_EQ(idef->il[F_BONDS].iatoms, bondIatoms);
+ std::vector<int> angleIatoms{ 2, 0, 1, 2, 3, 1, 2, 3 };
+ EXPECT_EQ(idef->il[F_ANGLES].iatoms, angleIatoms);
+ idef->clear();
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib 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 "gmxpre.h"
+
+#include <gtest/gtest.h>
+
+#include "nblib/listed_forces/helpers.hpp"
+#include "nblib/listed_forces/traits.h"
+
+#include "testutils/testasserts.h"
+
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+TEST(NBlibTest, CanSplitListedWork)
+{
+ ListedInteractionData interactions;
+
+ DefaultAngle angle(Degrees(1), 1);
+ HarmonicBondType bond(1, 1);
+
+ int largestIndex = 20;
+ int nSplits = 3; // split ranges: [0,5], [6,11], [12, 19]
+
+ std::vector<InteractionIndex<HarmonicBondType>> bondIndices{
+ { 0, 1, 0 }, { 0, 6, 0 }, { 11, 12, 0 }, { 18, 19, 0 }
+ };
+ std::vector<InteractionIndex<DefaultAngle>> angleIndices{
+ { 0, 1, 2, 0 }, { 0, 6, 7, 0 }, { 11, 12, 13, 0 }, { 17, 19, 18, 0 }
+ };
+
+ pickType<HarmonicBondType>(interactions).indices = bondIndices;
+ pickType<DefaultAngle>(interactions).indices = angleIndices;
+
+ std::vector<ListedInteractionData> splitInteractions =
+ splitListedWork(interactions, largestIndex, nSplits);
+
+ std::vector<InteractionIndex<HarmonicBondType>> refBondIndices0{ { 0, 1, 0 }, { 0, 6, 0 } };
+ std::vector<InteractionIndex<DefaultAngle>> refAngleIndices0{ { 0, 1, 2, 0 }, { 0, 6, 7, 0 } };
+ std::vector<InteractionIndex<HarmonicBondType>> refBondIndices1{ { 11, 12, 0 } };
+ std::vector<InteractionIndex<DefaultAngle>> refAngleIndices1{ { 11, 12, 13, 0 } };
+ std::vector<InteractionIndex<HarmonicBondType>> refBondIndices2{ { 18, 19, 0 } };
+ std::vector<InteractionIndex<DefaultAngle>> refAngleIndices2{ { 17, 19, 18, 0 } };
+
+ EXPECT_EQ(refBondIndices0, pickType<HarmonicBondType>(splitInteractions[0]).indices);
+ EXPECT_EQ(refBondIndices1, pickType<HarmonicBondType>(splitInteractions[1]).indices);
+ EXPECT_EQ(refBondIndices2, pickType<HarmonicBondType>(splitInteractions[2]).indices);
+
+ EXPECT_EQ(refAngleIndices0, pickType<DefaultAngle>(splitInteractions[0]).indices);
+ EXPECT_EQ(refAngleIndices1, pickType<DefaultAngle>(splitInteractions[1]).indices);
+ EXPECT_EQ(refAngleIndices2, pickType<DefaultAngle>(splitInteractions[2]).indices);
+}
+
+
+TEST(NBlibTest, ListedForceBuffer)
+{
+ using T = gmx::RVec;
+ int ncoords = 20;
+
+ T vzero{ 0, 0, 0 };
+ std::vector<T> masterBuffer(ncoords, vzero);
+
+ // the ForceBuffer is going to access indices [10-15) through the masterBuffer
+ // and the outliers internally
+ int rangeStart = 10;
+ int rangeEnd = 15;
+
+ ForceBuffer<T> forceBuffer(masterBuffer.data(), rangeStart, rangeEnd);
+
+ // in range
+ T internal1{ 1, 2, 3 };
+ T internal2{ 4, 5, 6 };
+ forceBuffer[10] = internal1;
+ forceBuffer[14] = internal2;
+
+ // outliers
+ T outlier{ 0, 1, 2 };
+ forceBuffer[5] = outlier;
+
+ std::vector<T> refMasterBuffer(ncoords, vzero);
+ refMasterBuffer[10] = internal1;
+ refMasterBuffer[14] = internal2;
+
+ for (size_t i = 0; i < masterBuffer.size(); ++i)
+ {
+ for (size_t m = 0; m < dimSize; ++m)
+ {
+ EXPECT_REAL_EQ_TOL(refMasterBuffer[i][m], masterBuffer[i][m],
+ gmx::test::defaultRealTolerance());
+ }
+ }
+
+ for (size_t m = 0; m < dimSize; ++m)
+ {
+ EXPECT_REAL_EQ_TOL(outlier[m], forceBuffer[5][m], gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(vzero[m], forceBuffer[4][m], gmx::test::defaultRealTolerance());
+ }
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements kernel 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/kernels.hpp"
+
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+
+TEST(Kernels, HarmonicScalarKernelCanCompute)
+{
+ real k = 1.1;
+ real x0 = 1.0;
+ real x = 1.2;
+
+ real force, epot;
+ std::tie(force, epot) = harmonicScalarForce(k, x0, x);
+
+ EXPECT_REAL_EQ_TOL(-k * (x - x0), force, gmx::test::absoluteTolerance(1e-4));
+ EXPECT_REAL_EQ_TOL(0.5 * k * (x - x0) * (x - x0), epot, gmx::test::absoluteTolerance(1e-4));
+}
+
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib 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>
+ */
+
+#ifndef NBLIB_LINEAR_CHAIN_DATA_HPP
+#define NBLIB_LINEAR_CHAIN_DATA_HPP
+
+#include "nblib/box.h"
+#include "nblib/topologyhelpers.h"
+#include "nblib/listed_forces/traits.h"
+
+namespace nblib
+{
+
+//! \brief sets up an interaction tuples for a linear chain with nParticles
+class LinearChainData
+{
+public:
+ LinearChainData(int nP, float outlierRatio = 0.02) : nParticles(nP)
+ {
+ HarmonicBondType bond1{ 376560, real(0.1001 * std::sqrt(3)) };
+ HarmonicBondType bond2{ 313800, real(0.1001 * std::sqrt(3)) };
+ std::vector<HarmonicBondType> bonds{ bond1, bond2 };
+ pickType<HarmonicBondType>(interactions).parameters = bonds;
+
+ DefaultAngle angle(Degrees(179.9), 397.5);
+ std::vector<DefaultAngle> angles{ angle };
+ pickType<DefaultAngle>(interactions).parameters = angles;
+
+ std::vector<InteractionIndex<HarmonicBondType>> bondIndices;
+ for (int i = 0; i < nParticles - 1; ++i)
+ {
+ // third index: alternate between the two bond parameters
+ bondIndices.push_back(InteractionIndex<HarmonicBondType>{ i, i + 1, i % 2 });
+ }
+ pickType<HarmonicBondType>(interactions).indices = bondIndices;
+
+ std::vector<InteractionIndex<DefaultAngle>> angleIndices;
+ for (int i = 0; i < nParticles - 2; ++i)
+ {
+ angleIndices.push_back(InteractionIndex<DefaultAngle>{ i, i + 1, i + 2, 0 });
+ }
+ pickType<DefaultAngle>(interactions).indices = angleIndices;
+
+ // initialize coordinates
+ x.resize(nParticles);
+ for (int i = 0; i < nParticles; ++i)
+ {
+ x[i] = real(i) * gmx::RVec{ 0.1, 0.1, 0.1 };
+ }
+
+ forces = std::vector<gmx::RVec>(nParticles, gmx::RVec{ 0, 0, 0 });
+
+ box.reset(new Box(0.1 * (nParticles + 1), 0.1 * (nParticles + 1), 0.1 * (nParticles + 1)));
+
+ addOutliers(outlierRatio);
+ }
+
+ // void createAggregates() { aggregateTransformations(interactions); }
+
+ void addOutliers(float outlierRatio)
+ {
+ HarmonicBondType dummyBond{ 1e-6, real(0.1001 * std::sqrt(3)) };
+ pickType<HarmonicBondType>(interactions).parameters.push_back(dummyBond);
+ int bondIndex = pickType<HarmonicBondType>(interactions).parameters.size() - 1;
+
+ int nOutliers = nParticles * outlierRatio;
+ srand(42);
+ for (int s = 0; s < nOutliers; ++s)
+ {
+ int from = rand() % nParticles;
+ int to = rand() % nParticles;
+ if (from != to)
+ {
+ pickType<HarmonicBondType>(interactions)
+ .indices.push_back(InteractionIndex<HarmonicBondType>{ from, to, bondIndex });
+ }
+ }
+ }
+
+ int nParticles;
+
+ std::vector<gmx::RVec> x;
+ std::vector<gmx::RVec> forces;
+
+ ListedInteractionData interactions;
+
+ std::shared_ptr<Box> box;
+};
+
+} // namespace nblib
+
+#endif // NBLIB_LINEAR_CHAIN_DATA_HPP
+
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">15.644075054627638</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">4</Int>
+ <Vector>
+ <Real Name="X">0</Real>
+ <Real Name="Y">3000.9773757382213</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-294.15662928338577</Real>
+ <Real Name="Y">2953.652678308878</Real>
+ <Real Name="Z">-14.707831464169288</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">589.93598821678029</Real>
+ <Real Name="Y">-5922.0943245644257</Real>
+ <Real Name="Z">29.49679941083901</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-295.77935893339446</Real>
+ <Real Name="Y">-32.535729482673389</Real>
+ <Real Name="Z">-14.788967946669722</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">11.117714323462188</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">4</Int>
+ <Vector>
+ <Real Name="X">273.20508075688775</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-858.64453952164718</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">957.99184161506105</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-372.55238285030157</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">11.117714323462197</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">4</Int>
+ <Vector>
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">253.81712522292264</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">-884.99823613097033</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">299.4154280846883</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">331.76568282335938</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib utility tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+#include <numeric>
+
+#include <gtest/gtest.h>
+
+#include "nblib/listed_forces/traits.h"
+#include "nblib/listed_forces/transformations.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+namespace test
+{
+namespace
+{
+
+ListedInteractionData unsortedInteractions()
+{
+ ListedInteractionData interactions;
+
+ std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 2, 0 }, { 0, 1, 0 } };
+ pickType<HarmonicBondType>(interactions).indices = std::move(bondIndices);
+
+ std::vector<InteractionIndex<DefaultAngle>> angleIndices{ { 0, 1, 2, 0 }, { 1, 0, 2, 0 } };
+ pickType<DefaultAngle>(interactions).indices = std::move(angleIndices);
+
+ std::vector<InteractionIndex<ProperDihedral>> dihedralIndices{ { 0, 2, 1, 3, 0 }, { 0, 1, 2, 3, 0 } };
+ pickType<ProperDihedral>(interactions).indices = std::move(dihedralIndices);
+
+ return interactions;
+}
+
+TEST(ListedTransformations, SortInteractionIndices)
+{
+ ListedInteractionData interactions = unsortedInteractions();
+ sortInteractions(interactions);
+
+ std::vector<InteractionIndex<HarmonicBondType>> refBondIndices{ { 0, 1, 0 }, { 0, 2, 0 } };
+ std::vector<InteractionIndex<DefaultAngle>> refAngleIndices{ { 1, 0, 2, 0 }, { 0, 1, 2, 0 } };
+ std::vector<InteractionIndex<ProperDihedral>> refDihedralIndices{ { 0, 1, 2, 3, 0 },
+ { 0, 2, 1, 3, 0 } };
+
+ EXPECT_EQ(pickType<HarmonicBondType>(interactions).indices, refBondIndices);
+ EXPECT_EQ(pickType<DefaultAngle>(interactions).indices, refAngleIndices);
+ EXPECT_EQ(pickType<ProperDihedral>(interactions).indices, refDihedralIndices);
+}
+
+} // namespace
+} // namespace test
+} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This implements basic nblib 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/dataflow.hpp"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace nblib
+{
+
+//! Number of atoms used in these tests.
+constexpr int c_numAtoms = 4;
+
+//! Coordinates for testing
+std::vector<std::vector<gmx::RVec>> c_coordinatesForTests = {
+ { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 } },
+ { { 0.5, 0.0, 0.0 }, { 0.5, 0.0, 0.15 }, { 0.5, 0.07, 0.22 }, { 0.5, 0.18, 0.22 } },
+ { { -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 } }
+};
+
+//! Function types for testing dihedrals. Add new terms at the end.
+std::vector<std::vector<ProperDihedral>> c_InputDihs = { { { ProperDihedral(Degrees(-105.0), 15.0, 2) } } /*, { ImproperDihedral(100.0, 50.0) }*/ };
+// Todo: update test setup to allow more than one interaction type and add the following to the inputs
+// std::vector<std::vector<RyckaertBellemanDihedral>> c_InputDihs = { { RyckaertBellemanDihedral({ -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 }) } };
+
+template<class Interaction>
+class ListedForcesBase
+{
+public:
+ std::vector<Interaction> input_;
+ std::vector<gmx::RVec> x_;
+ std::vector<InteractionIndex<Interaction>> indices_;
+ PbcHolder pbcHolder_;
+ gmx::test::TestReferenceData refData_;
+ gmx::test::TestReferenceChecker checker_;
+ std::vector<gmx::RVec> forces_;
+ real energy_;
+
+ ListedForcesBase(std::vector<Interaction> input,
+ std::vector<gmx::RVec> coordinates,
+ std::vector<InteractionIndex<Interaction>> indices) :
+ input_(std::move(input)),
+ x_(std::move(coordinates)),
+ indices_(std::move(indices)),
+ pbcHolder_(Box(1.5)),
+ checker_(refData_.rootChecker()),
+ forces_(c_numAtoms, gmx::RVec{ 0, 0, 0 })
+ {
+ energy_ = computeForces(indices_, input_, x_, &forces_, pbcHolder_);
+ }
+
+ void checkForcesAndEnergies()
+ {
+ checker_.checkReal(energy_, "Epot");
+ checker_.checkSequence(std::begin(forces_), std::end(forces_), "forces");
+ }
+};
+
+class ListedForcesProperDihedralTest :
+ public ListedForcesBase<ProperDihedral>,
+ public testing::TestWithParam<std::tuple<std::vector<ProperDihedral>, std::vector<gmx::RVec>>>
+{
+ using Base = ListedForcesBase<ProperDihedral>;
+
+public:
+ ListedForcesProperDihedralTest() :
+ Base(std::get<0>(GetParam()), std::get<1>(GetParam()), { { 0, 1, 2, 3, 0 } })
+ {
+ }
+};
+
+TEST_P(ListedForcesProperDihedralTest, CheckListed)
+{
+ checkForcesAndEnergies();
+}
+
+INSTANTIATE_TEST_CASE_P(FourCenter,
+ ListedForcesProperDihedralTest,
+ ::testing::Combine(::testing::ValuesIn(c_InputDihs),
+ ::testing::ValuesIn(c_coordinatesForTests)));
+
+} // 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
+ *
+ * Implement transformations and manipulations of ListedInteractionData,
+ * such as sorting and splitting
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+#include "nblib/listed_forces/transformations.h"
+
+#include "nblib/util/internal.h"
+
+namespace nblib
+{
+
+void sortInteractions(ListedInteractionData& interactions)
+{
+ auto sortKeyObj = [](const auto& lhs, const auto& rhs) { return interactionSortKey(lhs, rhs); };
+ auto sortOneElement = [comparison = sortKeyObj](auto& interactionElement) {
+ std::sort(begin(interactionElement.indices), end(interactionElement.indices), comparison);
+ };
+
+ for_each_tuple(sortOneElement, interactions);
+}
+
+} // 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
+ *
+ * Implement transformations and manipulations of ListedInteractionData,
+ * such as sorting and splitting
+ *
+ * \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_TRANSFORMATIONS_H
+#define NBLIB_LISTEDFORCES_TRANSFORMATIONS_H
+
+#include <algorithm>
+
+#include "definitions.h"
+
+namespace nblib
+{
+
+namespace detail
+{
+
+/*! \brief tuple ordering for two center interactions
+ *
+ * \param t input array
+ * \return ordered array
+ */
+inline std::array<int, 2> nblibOrdering(const std::array<int, 2>& t)
+{
+ // for bonds (two coordinate indices),
+ // we choose to store the lower sequence ID first. this allows for better unit tests
+ // that are agnostic to how the input was set up
+ int id1 = std::min(std::get<0>(t), std::get<1>(t));
+ int id2 = std::max(std::get<0>(t), std::get<1>(t));
+
+ return std::array<int, 2>{ id1, id2 };
+}
+
+/*! \brief tuple ordering for three center interactions
+ *
+ * \param t input array
+ * \return ordered array
+ */
+inline std::array<int, 3> nblibOrdering(const std::array<int, 3>& t)
+{
+ // for angles (three coordinate indices),
+ // we choose to store the two non-center coordinate indices sorted.
+ // such that ret[0] < ret[2] always (ret = returned tuple)
+ int id1 = std::min(std::get<0>(t), std::get<2>(t));
+ int id3 = std::max(std::get<0>(t), std::get<2>(t));
+
+ return std::array<int, 3>{ id1, std::get<1>(t), id3 };
+}
+
+/*! \brief tuple ordering for four center interactions
+ *
+ * \param t input array
+ * \return ordered array
+ */
+inline std::array<int, 4> nblibOrdering(const std::array<int, 4>& t)
+{
+ return t;
+}
+
+/*! \brief tuple ordering for five center interactions
+ *
+ * \param t input array
+ * \return ordered array
+ */
+inline std::array<int, 5> nblibOrdering(const std::array<int, 5>& t)
+{
+ return t;
+}
+
+} // namespace detail
+
+//! \brief sort key function object to sort 2-center interactions
+inline bool interactionSortKey(const TwoCenterInteractionIndex& lhs, const TwoCenterInteractionIndex& rhs)
+{
+ return lhs < rhs;
+}
+
+//! \brief sort key function object to sort 3-center interactions
+inline bool interactionSortKey(const ThreeCenterInteractionIndex& lhs, const ThreeCenterInteractionIndex& rhs)
+{
+ // position [1] is the center atom of the angle and is the only sort key
+ // to allow use of std::equal_range to obtain a range of all angles with a given central atom
+ return lhs[1] < rhs[1];
+}
+
+//! \brief sort key function object to sort 4-center interactions
+inline bool interactionSortKey(const FourCenterInteractionIndex& lhs, const FourCenterInteractionIndex& rhs)
+{
+ // we only take the first center-axis-particle into account
+ // this allows use of std::equal_range to find all four-center interactions with a given j-index
+ // Note that there exists a second version that compares the k-index which is used in
+ // aggregate transformation
+ return lhs[1] < rhs[1];
+}
+
+//! \brief sort key function object to sort 5-center interactions
+inline bool interactionSortKey(const FiveCenterInteractionIndex& lhs, const FiveCenterInteractionIndex& rhs)
+{
+ return lhs < rhs;
+}
+
+//! \brief sorts all interaction indices according to the keys defined in the implementation
+void sortInteractions(ListedInteractionData& interactions);
+
+} // namespace nblib
+#endif // NBLIB_LISTEDFORCES_TRANSFORMATIONS_H
return *this;
}
+//! Two-particle interactions such as bonds and LJ1-4
+template<class Interaction>
+void Molecule::addInteraction(const ParticleName& particleNameI,
+ const ResidueName& residueNameI,
+ const ParticleName& particleNameJ,
+ const ResidueName& residueNameJ,
+ const Interaction& interaction)
+{
+ if (particleNameI == particleNameJ and residueNameI == residueNameJ)
+ {
+ throw InputException(std::string("Cannot add interaction of particle ")
+ + particleNameI.value() + " with itself in molecule " + name_.value());
+ }
+
+ auto& interactionContainer = pickType<Interaction>(interactionData_);
+ interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ, residueNameJ);
+ interactionContainer.interactionTypes_.push_back(interaction);
+}
+
+//! \cond DO_NOT_DOCUMENT
+#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
+ template void Molecule::addInteraction( \
+ const ParticleName& particleNameI, const ResidueName& residueNameI, \
+ const ParticleName& particleNameJ, const ResidueName& residueNameJ, const x& interaction);
+MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
+#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
+//! \endcond
+
+// add interactions with default residue name
+template<class Interaction>
+void Molecule::addInteraction(const ParticleName& particleNameI,
+ const ParticleName& particleNameJ,
+ const Interaction& interaction)
+{
+ addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_), interaction);
+}
+
+//! \cond DO_NOT_DOCUMENT
+#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
+ template void Molecule::addInteraction(const ParticleName& particleNameI, \
+ const ParticleName& particleNameJ, const x& interaction);
+MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
+#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
+
+//! 3-particle interactions such as angles
+template<class Interaction>
+void Molecule::addInteraction(const ParticleName& particleNameI,
+ const ResidueName& residueNameI,
+ const ParticleName& particleNameJ,
+ const ResidueName& residueNameJ,
+ const ParticleName& particleNameK,
+ const ResidueName& residueNameK,
+ const Interaction& interaction)
+{
+ if (particleNameI == particleNameJ and particleNameJ == particleNameK)
+ {
+ throw InputException(std::string("Cannot add interaction of particle ")
+ + particleNameI.value() + " with itself in molecule " + name_.value());
+ }
+
+ auto& interactionContainer = pickType<Interaction>(interactionData_);
+ interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ,
+ residueNameJ, particleNameK, residueNameK);
+ interactionContainer.interactionTypes_.push_back(interaction);
+}
+
+#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
+ 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_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
+#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
+
+template<class Interaction>
+void Molecule::addInteraction(const ParticleName& particleNameI,
+ const ParticleName& particleNameJ,
+ const ParticleName& particleNameK,
+ const Interaction& interaction)
+{
+ addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_),
+ particleNameK, ResidueName(name_), interaction);
+}
+
+#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
+ template void Molecule::addInteraction(const ParticleName& particleNameI, \
+ const ParticleName& particleNameJ, \
+ const ParticleName& particleNameK, const x& interaction);
+MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
+#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
+//! \endcond
+
int Molecule::numParticlesInMolecule() const
{
return particles_.size();
std::make_tuple(particleNameToExclude, ResidueName(name_)));
}
+const Molecule::InteractionTuple& Molecule::interactionData() const
+{
+ return interactionData_;
+}
+
const ParticleType& Molecule::at(const std::string& particleTypeName) const
{
return particleTypes_.at(particleTypeName);
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements a RAII encapsulation for PbcAiuc
+ * note this is a header implementation to make inlining easier if needed
+ *
+ * \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_PBC_HPP
+#define NBLIB_PBC_HPP
+
+#include "nblib/box.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc_aiuc.h"
+
+struct PbcAiuc;
+
+namespace nblib
+{
+
+/*! \brief life-time manager for PbcAiuc
+ *
+ */
+class PbcHolder
+{
+public:
+ //! \brief ctor
+ explicit PbcHolder(const Box& box)
+ {
+ setPbcAiuc(3, box.legacyMatrix(), &m_pbc);
+ }
+
+ //! \brief calculate pbc-aware r1-r2
+ inline void dxAiuc(const rvec& r1, const rvec& r2, rvec dr) const
+ {
+ ::pbcDxAiuc(m_pbc, r1, r2, dr);
+ }
+
+private:
+ PbcAiuc m_pbc;
+};
+
+/*! \brief dummy class used to turn off Pbc in listed forces
+ *
+ */
+class NoPbc
+{
+public:
+
+ //! \brief calculate r1-r2, ignore pbc
+ inline void dxAiuc(const rvec& r1, const rvec& r2, rvec dr) const
+ {
+ rvec_sub(r1, r2, dr);
+ }
+};
+
+} // namespace nblib
+#endif // NBLIB_PBC_HPP
box.cpp
interactions.cpp
particletype.cpp
+ pbcholder.cpp
molecules.cpp
topology.cpp
)
EXPECT_NO_THROW(molecule.addParticle(ParticleName("U2"), atom2));
}
+TEST(NBlibTest, CanAddInteractions)
+{
+ Molecule molecule(MoleculeName("BondTest"));
+ ParticleType O(ParticleTypeName("Ow"), Mass(1));
+ ParticleType H(ParticleTypeName("Hw"), Mass(1));
+ molecule.addParticle(ParticleName("O"), O);
+ molecule.addParticle(ParticleName("H1"), H);
+ molecule.addParticle(ParticleName("H2"), H);
+
+ HarmonicBondType hb(1, 2);
+ CubicBondType cub(1, 2, 3);
+ DefaultAngle ang(Degrees(1), 1);
+
+ molecule.addInteraction(ParticleName("O"), ParticleName("H1"), hb);
+ molecule.addInteraction(ParticleName("O"), ParticleName("H2"), hb);
+ molecule.addInteraction(ParticleName("H1"), ParticleName("H2"), cub);
+ molecule.addInteraction(ParticleName("H1"), ParticleName("O"), ParticleName("H2"), ang);
+
+ const auto& interactionData = molecule.interactionData();
+
+ //! harmonic bonds
+ EXPECT_EQ(pickType<HarmonicBondType>(interactionData).interactions_.size(), 2);
+ //! cubic bonds
+ EXPECT_EQ(pickType<CubicBondType>(interactionData).interactions_.size(), 1);
+ //! angular interactions
+ EXPECT_EQ(pickType<DefaultAngle>(interactionData).interactions_.size(), 1);
+}
+
} // namespace
} // namespace test
} // namespace nblib
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * This PbcHolder 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 <cmath>
+
+#include "nblib/pbc.hpp"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+using gmx::test::defaultRealTolerance;
+
+namespace nblib
+{
+
+TEST(NBlibTest, PbcHolderWorks)
+{
+ Box box(10, 10, 10);
+
+ PbcHolder pbcHolder(box);
+
+ rvec x1{ 1, 1, 1 }, x2{ 9, 9, 9 };
+ rvec dx;
+
+ pbcHolder.dxAiuc(x1, x2, dx);
+ rvec ref{ 2, 2, 2 };
+
+ EXPECT_REAL_EQ_TOL(ref[0], dx[0], gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(ref[1], dx[1], gmx::test::defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(ref[2], dx[2], gmx::test::defaultRealTolerance());
+}
+
+} // namespace nblib
water_.addParticle(ParticleName("Oxygen"), Charges.at("Ow"), plib.type("Ow"));
water_.addParticle(ParticleName("H1"), Charges.at("Hw"), plib.type("H"));
water_.addParticle(ParticleName("H2"), Charges.at("Hw"), plib.type("H"));
+
+ HarmonicBondType ohBond(1., 1.);
+ water_.addInteraction(ParticleName("Oxygen"), ParticleName("H1"), ohBond);
+ water_.addInteraction(ParticleName("Oxygen"), ParticleName("H2"), ohBond);
}
Molecule WaterMoleculeBuilder::waterMolecule()
methanol_.addExclusion(ParticleName("Me1"), ParticleName("O2"));
methanol_.addExclusion(ParticleName("Me1"), ParticleName("H3"));
methanol_.addExclusion(ParticleName("H3"), ParticleName("O2"));
+
+ HarmonicBondType ohBond(1.01, 1.02);
+ methanol_.addInteraction(ParticleName("O2"), ParticleName("H3"), ohBond);
+
+ HarmonicBondType ometBond(1.1, 1.2);
+ methanol_.addInteraction(ParticleName("O2"), ParticleName("Me1"), ometBond);
+
+ DefaultAngle ochAngle(Degrees(108.52), 397.5);
+ methanol_.addInteraction(ParticleName("O2"), ParticleName("Me1"), ParticleName("H3"), ochAngle);
}
Molecule MethanolMoleculeBuilder::methanolMolecule()
EXPECT_EQ(5, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H2")));
}
+TEST(NBlibTest, TopologyCanAggregateBonds)
+{
+ Molecule water = WaterMoleculeBuilder{}.waterMolecule();
+ Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
+
+ std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
+ std::make_tuple(methanol, 1) };
+ std::vector<HarmonicBondType> bonds;
+ std::vector<size_t> bondsExpansion;
+ std::tie(bondsExpansion, bonds) = detail::collectInteractions<HarmonicBondType>(molecules);
+
+ std::vector<HarmonicBondType> bondsTest;
+ // use the expansionArray (bondsExpansion) to expand to the full list if bonds
+ std::transform(begin(bondsExpansion), end(bondsExpansion), std::back_inserter(bondsTest),
+ [&bonds](size_t i) { return bonds[i]; });
+
+ std::vector<HarmonicBondType> waterBonds =
+ pickType<HarmonicBondType>(water.interactionData()).interactionTypes_;
+ std::vector<HarmonicBondType> methanolBonds =
+ pickType<HarmonicBondType>(methanol.interactionData()).interactionTypes_;
+
+ std::vector<HarmonicBondType> bondsReference;
+ std::copy(begin(waterBonds), end(waterBonds), std::back_inserter(bondsReference));
+ std::copy(begin(waterBonds), end(waterBonds), std::back_inserter(bondsReference));
+ std::copy(begin(methanolBonds), end(methanolBonds), std::back_inserter(bondsReference));
+
+ EXPECT_EQ(bondsTest, bondsReference);
+}
+
+TEST(NBlibTest, TopologyCanSequencePairIDs)
+{
+ Molecule water = WaterMoleculeBuilder{}.waterMolecule();
+ Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
+
+ std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
+ std::make_tuple(methanol, 1) };
+ detail::ParticleSequencer particleSequencer;
+ particleSequencer.build(molecules);
+ auto pairs = detail::sequenceIDs<HarmonicBondType>(molecules, particleSequencer);
+
+ int Ow1 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
+ int H11 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
+ int H12 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
+ int Ow2 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("Oxygen"));
+ int H21 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H1"));
+ int H22 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H2"));
+
+ int Me = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
+ int MeO = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
+ int MeH = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
+
+ /// \cond DO_NOT_DOCUMENT
+#define SORT(i, j) (i < j) ? i : j, (i < j) ? j : i
+
+ std::vector<CoordinateIndex<HarmonicBondType>> pairs_reference{
+ { SORT(Ow1, H11) }, { SORT(Ow1, H12) }, { SORT(Ow2, H21) },
+ { SORT(Ow2, H22) }, { SORT(MeO, MeH) }, { SORT(MeO, Me) }
+ };
+#undef SORT
+ /// \endcond
+
+ EXPECT_EQ(pairs, pairs_reference);
+}
+
+TEST(NBlibTest, TopologySequenceIdThrows)
+{
+ Molecule water = WaterMoleculeBuilder{}.waterMolecule();
+ Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
+
+ std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
+ std::make_tuple(methanol, 1) };
+ detail::ParticleSequencer particleSequencer;
+ particleSequencer.build(molecules);
+ auto pairs = detail::sequenceIDs<HarmonicBondType>(molecules, particleSequencer);
+
+ // Input error: no particle called O-Atom in molecule "water"
+ EXPECT_THROW(particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("O-Atom")),
+ InputException);
+}
+
+TEST(NBlibTest, TopologyCanEliminateDuplicateBonds)
+{
+ HarmonicBondType b1(1.0, 2.0);
+ HarmonicBondType b2(1.1, 2.1);
+ HarmonicBondType b3(1.2, 2.2);
+
+ // can be compressed to {b1,b2,b3} + {1,1,2,0,1,0,2,2}
+ std::vector<HarmonicBondType> bonds{ b2, b2, b3, b1, b2, b1, b3, b3 };
+
+ // expected output
+ std::vector<HarmonicBondType> uniqueBondsReference{ b1, b2, b3 };
+ std::vector<size_t> indicesReference{ 1, 1, 2, 0, 1, 0, 2, 2 };
+
+ std::tuple<std::vector<size_t>, std::vector<HarmonicBondType>> bondData =
+ detail::eliminateDuplicateInteractions(bonds);
+
+ auto indices = std::get<0>(bondData);
+ auto uniqueBonds = std::get<1>(bondData);
+
+ EXPECT_EQ(uniqueBondsReference, uniqueBonds);
+ EXPECT_EQ(indicesReference, indices);
+}
+
+TEST(NBlibTest, TopologyListedInteractions)
+{
+ // Todo: add angles to SpcMethanolTopologyBuilder and extend interactions_reference below to
+ // Todo: include angles
+
+ Topology spcTopology = SpcMethanolTopologyBuilder{}.buildTopology(1, 2);
+
+ auto interactionData = spcTopology.getInteractionData();
+ auto& harmonicBonds = pickType<HarmonicBondType>(interactionData);
+
+ auto& indices = harmonicBonds.indices;
+ auto& bonds = harmonicBonds.parameters;
+
+ std::map<std::tuple<int, int>, HarmonicBondType> interactions_test;
+ for (auto& ituple : indices)
+ {
+ interactions_test[std::make_tuple(std::get<0>(ituple), std::get<1>(ituple))] =
+ bonds[std::get<2>(ituple)];
+ }
+
+ // there should be 3 unique HarmonicBondType instances
+ EXPECT_EQ(bonds.size(), 3);
+ // and 6 interaction pairs (bonds)
+ EXPECT_EQ(indices.size(), 6);
+
+ HarmonicBondType ohBond(1., 1.);
+ HarmonicBondType ohBondMethanol(1.01, 1.02);
+ HarmonicBondType ometBond(1.1, 1.2);
+
+ std::map<std::tuple<int, int>, HarmonicBondType> interactions_reference;
+
+ int Ow = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
+ int H1 = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
+ int H2 = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
+
+ int Me1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
+ int MeO1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
+ int MeH1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
+
+ int Me2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("Me1"));
+ int MeO2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("O2"));
+ int MeH2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("H3"));
+
+ /// \cond DO_NOT_DOCUMENT
+#define SORT(i, j) (i < j) ? i : j, (i < j) ? j : i
+ interactions_reference[std::make_tuple(SORT(Ow, H1))] = ohBond;
+ interactions_reference[std::make_tuple(SORT(Ow, H2))] = ohBond;
+ interactions_reference[std::make_tuple(SORT(MeO1, MeH1))] = ohBondMethanol;
+ interactions_reference[std::make_tuple(SORT(MeO1, Me1))] = ometBond;
+ interactions_reference[std::make_tuple(SORT(MeO2, MeH2))] = ohBondMethanol;
+ interactions_reference[std::make_tuple(SORT(MeO2, Me2))] = ometBond;
+#undef SORT
+ /// \endcond
+
+ EXPECT_TRUE(std::equal(begin(interactions_reference), end(interactions_reference),
+ begin(interactions_test)));
+}
+
+TEST(NBlibTest, TopologyListedInteractionsMultipleTypes)
+{
+ // Todo: add an angle type here
+
+ Molecule water = WaterMoleculeBuilder{}.waterMolecule();
+ Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
+
+ CubicBondType testBond(1., 1., 1.);
+ DefaultAngle testAngle(Degrees(1), 1);
+
+ water.addInteraction(ParticleName("H1"), ParticleName("H2"), testBond);
+ water.addInteraction(ParticleName("H1"), ParticleName("Oxygen"), ParticleName("H2"), testAngle);
+
+ ParticleTypesInteractions nbInteractions;
+ std::vector<std::string> particleTypeNames = { "Ow", "H", "OMet", "CMet" };
+ for (const auto& name : particleTypeNames)
+ {
+ nbInteractions.add(ParticleTypeName(name), C6(0), C12(0));
+ }
+
+ TopologyBuilder topologyBuilder;
+ topologyBuilder.addMolecule(water, 1);
+ topologyBuilder.addMolecule(methanol, 1);
+ topologyBuilder.addParticleTypesInteractions(nbInteractions);
+
+ Topology topology = topologyBuilder.buildTopology();
+
+ auto interactionData = topology.getInteractionData();
+ auto& harmonicBonds = pickType<HarmonicBondType>(interactionData);
+ auto& cubicBonds = pickType<CubicBondType>(interactionData);
+ auto& angleInteractions = pickType<DefaultAngle>(interactionData);
+
+ HarmonicBondType ohBond(1., 1.);
+ HarmonicBondType ohBondMethanol(1.01, 1.02);
+ HarmonicBondType ometBond(1.1, 1.2);
+ std::vector<HarmonicBondType> harmonicBondsReference{ ohBond, ohBondMethanol, ometBond };
+
+ EXPECT_EQ(harmonicBonds.parameters, harmonicBondsReference);
+
+ int H1 = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
+ int H2 = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
+ int Ow = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
+
+ int Me1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
+ int MeO1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
+ int MeH1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
+
+ std::vector<CubicBondType> cubicBondsReference{ testBond };
+ std::vector<InteractionIndex<CubicBondType>> cubicIndicesReference{ { std::min(H1, H2),
+ std::max(H1, H2), 0 } };
+ EXPECT_EQ(cubicBondsReference, cubicBonds.parameters);
+ EXPECT_EQ(cubicIndicesReference, cubicBonds.indices);
+
+ DefaultAngle methanolAngle(Degrees(108.52), 397.5);
+ std::vector<DefaultAngle> angleReference{ testAngle, methanolAngle };
+ std::vector<InteractionIndex<DefaultAngle>> angleIndicesReference{
+ { std::min(H1, H2), Ow, std::max(H1, H2), 0 }, { std::min(MeH1, MeO1), Me1, std::max(MeO1, MeH1), 1 }
+ };
+ EXPECT_EQ(angleReference, angleInteractions.parameters);
+ EXPECT_EQ(angleIndicesReference, angleInteractions.indices);
+}
+
+TEST(NBlibTest, TopologyInvalidParticleInInteractionThrows)
+{
+ Molecule water = WaterMoleculeBuilder{}.waterMolecule();
+ Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
+
+ HarmonicBondType testBond(1., 1.);
+
+ // Invalid input: no particle named "Iron" in molecule water
+ water.addInteraction(ParticleName("H1"), ParticleName("Iron"), testBond);
+
+ ParticleTypesInteractions nbInteractions;
+ std::vector<std::string> particleTypeNames = { "Ow", "H", "OMet", "CMet" };
+ for (const auto& name : particleTypeNames)
+ {
+ nbInteractions.add(ParticleTypeName(name), C6(0), C12(0));
+ }
+
+ TopologyBuilder topologyBuilder;
+ topologyBuilder.addMolecule(water, 1);
+ topologyBuilder.addMolecule(methanol, 1);
+ topologyBuilder.addParticleTypesInteractions(nbInteractions);
+
+ EXPECT_THROW(topologyBuilder.buildTopology(), InputException);
+}
+
+
TEST(NBlibTest, toGmxExclusionBlockWorks)
{
std::vector<std::tuple<int, int>> testInput{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
return exclusionsListOfListsGlobal;
}
+ListedInteractionData TopologyBuilder::createInteractionData(const detail::ParticleSequencer& particleSequencer)
+{
+ ListedInteractionData interactionData;
+
+ // this code is doing the compile time equivalent of
+ // for (int i = 0; i < interactionData.size(); ++i)
+ // create(get<i>(interactionData));
+
+ auto create = [this, &particleSequencer](auto& interactionDataElement) {
+ using InteractionType = typename std::decay_t<decltype(interactionDataElement)>::type;
+
+ // first compression stage: each bond per molecule listed once,
+ // eliminates duplicates from multiple identical molecules
+ auto compressedDataStage1 = detail::collectInteractions<InteractionType>(this->molecules_);
+ auto& expansionArrayStage1 = std::get<0>(compressedDataStage1);
+ auto& moleculeInteractions = std::get<1>(compressedDataStage1);
+
+ // second compression stage: recognize bond duplicates among bonds from all molecules put together
+ auto compressedDataStage2 = detail::eliminateDuplicateInteractions(moleculeInteractions);
+ auto& expansionArrayStage2 = std::get<0>(compressedDataStage2);
+ auto& uniqueInteractionInstances = std::get<1>(compressedDataStage2);
+
+ // combine stage 1 + 2 expansion arrays
+ std::vector<size_t> expansionArray(expansionArrayStage1.size());
+ std::transform(begin(expansionArrayStage1), end(expansionArrayStage1), begin(expansionArray),
+ [& S2 = expansionArrayStage2](size_t S1Element) { return S2[S1Element]; });
+
+ // add data about InteractionType instances
+ interactionDataElement.parameters = std::move(uniqueInteractionInstances);
+
+ interactionDataElement.indices.resize(expansionArray.size());
+ // coordinateIndices contains the particle sequence IDs of all interaction coordinates of type <BondType>
+ auto coordinateIndices = detail::sequenceIDs<InteractionType>(this->molecules_, particleSequencer);
+ // zip coordinateIndices(i,j,...) + expansionArray(k) -> interactionDataElement.indices(i,j,...,k)
+ std::transform(begin(coordinateIndices), end(coordinateIndices), begin(expansionArray),
+ begin(interactionDataElement.indices),
+ [](auto coordinateIndex, auto interactionIndex) {
+ std::array<int, coordinateIndex.size() + 1> ret;
+ for (int i = 0; i < int(coordinateIndex.size()); ++i)
+ {
+ ret[i] = coordinateIndex[i];
+ }
+ ret[coordinateIndex.size()] = interactionIndex;
+ return ret;
+ // return std::tuple_cat(coordinateIndex, std::make_tuple(interactionIndex));
+ });
+ };
+
+ for_each_tuple(create, interactionData);
+
+ return interactionData;
+}
+
template<typename T, class Extractor>
std::vector<T> TopologyBuilder::extractParticleTypeQuantity(Extractor&& extractor)
{
topology_.combinationRule_ = particleTypesInteractions_.getCombinationRule();
topology_.nonBondedInteractionMap_ = particleTypesInteractions_.generateTable();
+ topology_.interactionData_ = createInteractionData(topology_.particleSequencer_);
+
// Check whether there is any missing term in the particleTypesInteractions compared to the
// list of particletypes
for (const auto& particleType1 : particleTypes_)
return nonBondedInteractionMap_;
}
+ListedInteractionData Topology::getInteractionData() const
+{
+ return interactionData_;
+}
+
CombinationRule Topology::getCombinationRule() const
{
return combinationRule_;
#include "gromacs/topology/exclusionblocks.h"
#include "gromacs/utility/smalloc.h"
#include "nblib/exception.h"
+#include "nblib/listed_forces/transformations.h"
#include "nblib/topologyhelpers.h"
#include "nblib/util/internal.h"
}
}
+template<class I>
+std::tuple<std::vector<size_t>, std::vector<I>>
+collectInteractions(const std::vector<std::tuple<Molecule, int>>& molecules)
+{
+ std::vector<I> collectedBonds;
+ std::vector<size_t> expansionArray;
+ for (auto& molNumberTuple : molecules)
+ {
+ const Molecule& molecule = std::get<0>(molNumberTuple);
+ size_t numMols = std::get<1>(molNumberTuple);
+
+ auto& interactions = pickType<I>(molecule.interactionData()).interactionTypes_;
+
+ std::vector<size_t> moleculeExpansion(interactions.size());
+ // assign indices to the bonds in the current molecule, continue counting from
+ // the number of bonds seen so far (=collectedBonds.size())
+ std::iota(begin(moleculeExpansion), end(moleculeExpansion), collectedBonds.size());
+
+ std::copy(begin(interactions), end(interactions), std::back_inserter(collectedBonds));
+
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ std::copy(begin(moleculeExpansion), end(moleculeExpansion), std::back_inserter(expansionArray));
+ }
+ }
+ return std::make_tuple(expansionArray, collectedBonds);
+}
+
+/// \cond DO_NOT_DOCUMENT
+#define COLLECT_BONDS_INSTANTIATE_TEMPLATE(x) \
+ template std::tuple<std::vector<size_t>, std::vector<x>> collectInteractions( \
+ const std::vector<std::tuple<Molecule, int>>&);
+MAP(COLLECT_BONDS_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
+/// \endcond
+
+namespace sequence_detail
+{
+
+//! Helper function to convert a tuple of strings into a particle index sequence
+template<class Tuple, class F, class... Args, size_t... Is>
+auto stringsToIndices_impl(const Tuple& tuple, std::index_sequence<Is...> is, F&& f, Args... args)
+{
+ // return std::make_tuple(f(args..., std::get<2 * Is + 1>(tuple), std::get<2 * Is>(tuple))...);
+ ignore_unused(is);
+ return std::array<int, sizeof...(Is)>{ f(args..., std::get<2 * Is + 1>(tuple),
+ std::get<2 * Is>(tuple))... };
+}
+
+/*! \brief
+ * This takes a tuple<(string, string) * nCenter> from molecule
+ * where nCenter = 2 for bonds, 3 for angles and 4 for dihedrals
+ * each (ResidueName, ParticleName)-pair is converted to a particle sequence index
+ * by calling the supplied function object f, containing the particleSequencer at the call site
+ * Therefore, the return type is tuple<int * nCenter>
+ *
+ */
+template<class Tuple, class F, class... Args>
+auto stringsToIndices(const Tuple& tuple, F&& f, Args... args)
+{
+ auto is = std::make_index_sequence<std::tuple_size<Tuple>::value / 2>{};
+ return stringsToIndices_impl(tuple, is, std::forward<F>(f), args...);
+}
+
+//! Tuple ordering for two center interactions
+[[maybe_unused]] static std::array<int, 2> nblibOrdering(const std::array<int, 2>& t)
+{
+ // for bonds (two coordinate indices),
+ // we choose to store the lower sequence ID first. this allows for better unit tests
+ // that are agnostic to how the input was set up
+ int id1 = std::min(std::get<0>(t), std::get<1>(t));
+ int id2 = std::max(std::get<0>(t), std::get<1>(t));
+
+ return std::array<int, 2>{ id1, id2 };
+}
+
+//! Tuple ordering for three center interactions
+[[maybe_unused]] static std::array<int, 3> nblibOrdering(const std::array<int, 3>& t)
+{
+ // for angles (three coordinate indices),
+ // we choose to store the two non-center coordinate indices sorted.
+ // such that ret[0] < ret[2] always (ret = returned tuple)
+ int id1 = std::min(std::get<0>(t), std::get<2>(t));
+ int id3 = std::max(std::get<0>(t), std::get<2>(t));
+
+ return std::array<int, 3>{ id1, std::get<1>(t), id3 };
+}
+
+//! Tuple ordering for four center interactions
+[[maybe_unused]] static std::array<int, 4> nblibOrdering(const std::array<int, 4>& t)
+{
+ return t;
+}
+
+//! Tuple ordering for five center interactions
+[[maybe_unused]] static std::array<int, 5> nblibOrdering(const std::array<int, 5>& t)
+{
+ return t;
+}
+
+} // namespace sequence_detail
+
+//! For each interaction, translate particle identifiers (moleculeName, nr, residueName,
+//! particleName) to particle coordinate indices
+template<class B>
+std::vector<CoordinateIndex<B>> sequenceIDs(const std::vector<std::tuple<Molecule, int>>& molecules,
+ const detail::ParticleSequencer& particleSequencer)
+{
+ std::vector<CoordinateIndex<B>> coordinateIndices;
+
+ auto callSequencer = [&particleSequencer](const MoleculeName& moleculeName, int i,
+ const ResidueName& residueName,
+ const ParticleName& particleName) {
+ return particleSequencer(moleculeName, i, residueName, particleName);
+ };
+
+ // loop over all molecules
+ for (const auto& molNumberTuple : molecules)
+ {
+ const Molecule& molecule = std::get<0>(molNumberTuple);
+ size_t numMols = std::get<1>(molNumberTuple);
+
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ auto& interactions = pickType<B>(molecule.interactionData()).interactions_;
+ for (const auto& interactionString : interactions)
+ {
+ CoordinateIndex<B> index = sequence_detail::stringsToIndices(
+ interactionString, callSequencer, molecule.name(), i);
+ coordinateIndices.push_back(nblibOrdering(index));
+ }
+ }
+ }
+ return coordinateIndices;
+}
+
+/// \cond DO_NOT_DOCUMENT
+#define SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE(x) \
+ template std::vector<CoordinateIndex<x>> sequenceIDs<x>( \
+ const std::vector<std::tuple<Molecule, int>>&, const detail::ParticleSequencer&);
+MAP(SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
+#undef SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE
+/// \endcond
template<class I>
std::tuple<std::vector<size_t>, std::vector<I>> eliminateDuplicateInteractions(const std::vector<I>& aggregatedInteractions)
return make_tuple(uniqueIndices, uniquInteractionsInstances);
}
+/// \cond DO_NOT_DOCUMENT
+#define ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE(x) \
+ template std::tuple<std::vector<size_t>, std::vector<x>> eliminateDuplicateInteractions( \
+ const std::vector<x>& aggregatedBonds);
+MAP(ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES)
+#undef ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE
+/// \endcond
+
} // namespace detail
} // namespace nblib
#include <vector>
#include "nblib/basicdefinitions.h"
+#include "nblib/ppmap.h"
#include "nblib/vector.h"
namespace gmx
=============================================
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.
+It is easy to extend NB-LIB to support novel formulations of particle interactions by modifying 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)``.
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.
+This file contains one ``struct`` for each interaction type parameter set.
+New interaction types are added here as separate structs. There
+are no content requirements, but for convenience, the existing ``NAMED_MEBERS``
+macro in combination with inheriting from a ``std::tuple`` or ``std::array``
+may be used. The macro can be used to define the
+parameter names for the corresponding setters and getters.
+For example, ``NAMED_MEMBERS(forceConstant, equilDistance)`` will expand to
+
+.. code:: cpp
+
+ inline auto& forceConstant() { return std::get<0>(*this); }
+ inline auto& equilDistance() { return std::get<1>(*this); }
+ inline const auto& forceConstant() const { return std::get<0>(*this); }
+ inline const auto& equilDistance() const { return std::get<1>(*this); }
+
+Putting everything together, one could define the complete parameter set for a new interaction type as follows.
+
+.. code:: cpp
+
+ /*! \brief new bond type
+ *
+ * V(r; forceConstant, equilDistance, scaleFactor)
+ * = forceConstant * exp( (r - equilDistance) / scaleFactor)
+ */
+ struct NewBondType : public std::tuple<real, real, int>
+ {
+ NewBondType() = default;
+ NewBondType(ForceConstant f, EquilDistance d, ScaleFactor s) :
+ std::tuple<real, real, int>{ f, d, s }
+ {
+ }
+
+ NAMED_MEMBERS(forceConstant, equilDistance, scaleFactor)
+ };
.. _definitions.h:
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.
+provides.
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
}
template <class T>
- inline std::tuple<T, T> bondKernel(T dr, const NewBondType& bond)
+ inline auto bondKernel(T dr, const NewBondType& bond)
{
return newBondForce(bond.forceConstant(), bond.equilDistance(), bond.scaleFactor(), dr);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff