From e2eeab4aea5be826cabe858a5584830b4c1fe641 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Mon, 9 Nov 2020 14:55:16 +0000 Subject: [PATCH] Add backend for nblib listed forces API This adds the back-end for the nblib API defined in listed_forces/bondtypes.h and listed_forces/definitions.h. Additionally, header and source are added for the private implementations of these files. These are implementation details, and are not visible to or interactable with by the user, so they are not part of the public API. This backend is a ground up refactoring of gromacs, copying the force and energy kernels but rewriting the disipatch. There are tests that compare forces and energies to gromacs. Each function in this commit has at least one test for correctness, though in most cases there are many more since the functionality is cumulative. All of this functionality has undergone extensive design review and code review already, so it should not be necessary to go through the whole change with a fine-tooth comb. Suggestions on how implementations could be improved are welcome. However, they will likely need to be addressed as follow-up issues unless they are trivial to fix. --- api/nblib/listed_forces/CMakeLists.txt | 6 + api/nblib/listed_forces/calculator.cpp | 171 +++++ api/nblib/listed_forces/conversions.hpp | 278 +++++++++ api/nblib/listed_forces/dataflow.hpp | 417 +++++++++++++ api/nblib/listed_forces/gmxcalculator.cpp | 129 ++++ api/nblib/listed_forces/gmxcalculator.h | 116 ++++ api/nblib/listed_forces/helpers.hpp | 194 ++++++ api/nblib/listed_forces/kernels.hpp | 587 ++++++++++++++++++ api/nblib/listed_forces/tests/CMakeLists.txt | 6 + api/nblib/listed_forces/tests/calculator.cpp | 266 ++++++++ api/nblib/listed_forces/tests/conversions.cpp | 105 ++++ api/nblib/listed_forces/tests/helpers.cpp | 148 +++++ api/nblib/listed_forces/tests/kernels.cpp | 64 ++ .../tests/linear_chain_input.hpp | 134 ++++ ...ForcesProperDihedralTest_CheckListed_0.xml | 28 + ...ForcesProperDihedralTest_CheckListed_1.xml | 28 + ...ForcesProperDihedralTest_CheckListed_2.xml | 28 + .../listed_forces/tests/transformations.cpp | 94 +++ api/nblib/listed_forces/tests/typetests.cpp | 123 ++++ api/nblib/listed_forces/transformations.cpp | 64 ++ api/nblib/listed_forces/transformations.h | 148 +++++ api/nblib/molecules.cpp | 97 +++ api/nblib/pbc.hpp | 95 +++ api/nblib/tests/CMakeLists.txt | 1 + api/nblib/tests/molecules.cpp | 28 + api/nblib/tests/pbcholder.cpp | 73 +++ api/nblib/tests/testsystems.cpp | 13 + api/nblib/tests/topology.cpp | 249 ++++++++ api/nblib/topology.cpp | 60 ++ api/nblib/topologyhelpers.cpp | 151 +++++ api/nblib/util/user.h | 1 + docs/nblib/listed-dev.rst | 51 +- 32 files changed, 3939 insertions(+), 14 deletions(-) create mode 100644 api/nblib/listed_forces/calculator.cpp create mode 100644 api/nblib/listed_forces/conversions.hpp create mode 100644 api/nblib/listed_forces/dataflow.hpp create mode 100644 api/nblib/listed_forces/gmxcalculator.cpp create mode 100644 api/nblib/listed_forces/gmxcalculator.h create mode 100644 api/nblib/listed_forces/helpers.hpp create mode 100644 api/nblib/listed_forces/kernels.hpp create mode 100644 api/nblib/listed_forces/tests/calculator.cpp create mode 100644 api/nblib/listed_forces/tests/conversions.cpp create mode 100644 api/nblib/listed_forces/tests/helpers.cpp create mode 100644 api/nblib/listed_forces/tests/kernels.cpp create mode 100644 api/nblib/listed_forces/tests/linear_chain_input.hpp create mode 100644 api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_0.xml create mode 100644 api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_1.xml create mode 100644 api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_2.xml create mode 100644 api/nblib/listed_forces/tests/transformations.cpp create mode 100644 api/nblib/listed_forces/tests/typetests.cpp create mode 100644 api/nblib/listed_forces/transformations.cpp create mode 100644 api/nblib/listed_forces/transformations.h create mode 100644 api/nblib/pbc.hpp create mode 100644 api/nblib/tests/pbcholder.cpp diff --git a/api/nblib/listed_forces/CMakeLists.txt b/api/nblib/listed_forces/CMakeLists.txt index c4e8f33e63..55507984d5 100644 --- a/api/nblib/listed_forces/CMakeLists.txt +++ b/api/nblib/listed_forces/CMakeLists.txt @@ -38,6 +38,12 @@ # \author Sebastian Keller # +target_sources(nblib + PRIVATE + calculator.cpp + transformations.cpp + ) + if(GMX_INSTALL_NBLIB_API) install(FILES bondtypes.h diff --git a/api/nblib/listed_forces/calculator.cpp b/api/nblib/listed_forces/calculator.cpp new file mode 100644 index 0000000000..09e305d7b2 --- /dev/null +++ b/api/nblib/listed_forces/calculator.cpp @@ -0,0 +1,171 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ +#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(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>( + masterForceBuffer_.data(), rangeStart, rangeEnd)); + } +} + +void ListedForceCalculator::computeForcesAndEnergies(gmx::ArrayRef x, bool usePbc) +{ + energyBuffer_.fill(0); + std::vector::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 + if (thisBuffer.inRange(index)) + { + auto force = outlier.second; + masterForceBuffer_[index] += force; + } + } + } + } +} + +void ListedForceCalculator::compute(gmx::ArrayRef coordinates, gmx::ArrayRef 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 coordinates, + gmx::ArrayRef forces, + EnergyType& energies, + bool usePbc) +{ + compute(coordinates, forces, usePbc); + + energies = energyBuffer_; +} + +} // namespace nblib diff --git a/api/nblib/listed_forces/conversions.hpp b/api/nblib/listed_forces/conversions.hpp new file mode 100644 index 0000000000..1385963e93 --- /dev/null +++ b/api/nblib/listed_forces/conversions.hpp @@ -0,0 +1,278 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ + +#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 +struct ListedIndex +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template<> +struct ListedIndex : std::integral_constant +{ +}; + +template +struct ListedTypeIsImplemented : std::bool_constant>::value> +{ +}; + +namespace detail +{ + +template +void transferParameters([[maybe_unused]] const InteractionData& interactionData, + [[maybe_unused]] gmx_ffparams_t& gmx_params) +{ +} + +template<> +void transferParameters(const TwoCenterData& 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& 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& 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 +void transferIndicesImpl(const TwoCenterData& interactions, InteractionDefinitions& idef, int offset) +{ + for (const auto& index : interactions.indices) + { + int parameterIndex = index[2] + offset; + idef.il[ListedIndex::value].iatoms.push_back(parameterIndex); + idef.il[ListedIndex::value].iatoms.push_back(index[0]); + idef.il[ListedIndex::value].iatoms.push_back(index[1]); + } +} + +template +void transferIndicesImpl(const ThreeCenterData& interactions, InteractionDefinitions& idef, int offset) +{ + for (const auto& index : interactions.indices) + { + int parameterIndex = index[3] + offset; + idef.il[ListedIndex::value].iatoms.push_back(parameterIndex); + idef.il[ListedIndex::value].iatoms.push_back(index[0]); + idef.il[ListedIndex::value].iatoms.push_back(index[1]); + idef.il[ListedIndex::value].iatoms.push_back(index[2]); + } +} + +template +void transferIndicesImpl(const FourCenterData& interactions, InteractionDefinitions& idef, int offset) +{ + for (const auto& index : interactions.indices) + { + int parameterIndex = index[4] + offset; + idef.il[ListedIndex::value].iatoms.push_back(parameterIndex); + idef.il[ListedIndex::value].iatoms.push_back(index[0]); + idef.il[ListedIndex::value].iatoms.push_back(index[1]); + idef.il[ListedIndex::value].iatoms.push_back(index[2]); + idef.il[ListedIndex::value].iatoms.push_back(index[3]); + } +} + +template +void transferIndicesImpl(const FiveCenterData& interactions, InteractionDefinitions& idef, int offset) +{ + for (const auto& index : interactions.indices) + { + int parameterIndex = index[5] + offset; + idef.il[ListedIndex::value].iatoms.push_back(parameterIndex); + idef.il[ListedIndex::value].iatoms.push_back(index[0]); + idef.il[ListedIndex::value].iatoms.push_back(index[1]); + idef.il[ListedIndex::value].iatoms.push_back(index[2]); + idef.il[ListedIndex::value].iatoms.push_back(index[3]); + idef.il[ListedIndex::value].iatoms.push_back(index[4]); + } +} + +template class Container, class InteractionType> +void transferIndices(const Container& interactionData, + InteractionDefinitions& idef, + [[maybe_unused]] int offset) +{ + if constexpr (ListedTypeIsImplemented{}) + { + 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> +createFFparams(const ListedInteractionData& interactions); + +std::tuple, std::unique_ptr> +createFFparams(const ListedInteractionData& interactions) +{ + std::unique_ptr ffparamsHolder = std::make_unique(); + std::unique_ptr idefHolder = std::make_unique(*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> indexOffsets{0}; + auto extractNIndices = [&indexOffsets, &acc](const auto& interactionElement) + { + constexpr int elementIndex = FindIndex, 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, 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 diff --git a/api/nblib/listed_forces/dataflow.hpp b/api/nblib/listed_forces/dataflow.hpp new file mode 100644 index 0000000000..79b3676353 --- /dev/null +++ b/api/nblib/listed_forces/dataflow.hpp @@ -0,0 +1,417 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ + +#ifndef NBLIB_LISTEDFORCES_DATAFLOW_HPP +#define NBLIB_LISTEDFORCES_DATAFLOW_HPP + +#include +#include + +#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 +inline NBLIB_ALWAYS_INLINE +auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj) +{ + using ValueType = BasicVectorValueType_t; + + 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 +inline NBLIB_ALWAYS_INLINE +auto dispatchInteraction(const TwoCenterInteractionIndex& index, + const std::vector& bondInstances, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + KernelEnergy> 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 +inline NBLIB_ALWAYS_INLINE +std::enable_if_t::value, BasicVectorValueType_t> +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 +inline NBLIB_ALWAYS_INLINE +std::enable_if_t::value, BasicVectorValueType_t> +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 +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; + //! 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 +inline NBLIB_ALWAYS_INLINE +auto dispatchInteraction(const ThreeCenterInteractionIndex& index, + const std::vector& parameters, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + KernelEnergy> 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 +inline NBLIB_ALWAYS_INLINE +std::enable_if_t::value, BasicVectorValueType_t> +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; + 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 +inline NBLIB_ALWAYS_INLINE +std::enable_if_t::value, BasicVectorValueType_t> +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 +inline NBLIB_ALWAYS_INLINE +auto dispatchInteraction(const FourCenterInteractionIndex& index, + const std::vector& parameters, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + KernelEnergy> 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 +inline NBLIB_ALWAYS_INLINE +auto dispatchInteraction(const FiveCenterInteractionIndex& index, + const std::vector& parameters, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + KernelEnergy> 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 +auto computeForces(const std::vector& indices, + const std::vector& interactionParameters, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + KernelEnergy> 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 +auto reduceListedForces(const ListedInteractionData& interactions, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + using ValueType = BasicVectorValueType_t; + std::array::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::type; + + KernelEnergy energy = computeForces(interactionElement.indices, interactionElement.parameters, x, forces, pbc); + + energies[CarrierIndex{}] += energy.carrier(); + energies[TwoCenterAggregateIndex{}] += energy.twoCenterAggregate(); + energies[ThreeCenterAggregateIndex{}] += 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 diff --git a/api/nblib/listed_forces/gmxcalculator.cpp b/api/nblib/listed_forces/gmxcalculator.cpp new file mode 100644 index 0000000000..9aa215339d --- /dev/null +++ b/api/nblib/listed_forces/gmxcalculator.cpp @@ -0,0 +1,129 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#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(&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& x, + std::vector& 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::type; + if constexpr (ListedTypeIsImplemented{}) + { + constexpr int index = FindIndex::value; + energies[index] = this->enerd.term[ListedIndex::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 diff --git a/api/nblib/listed_forces/gmxcalculator.h b/api/nblib/listed_forces/gmxcalculator.h new file mode 100644 index 0000000000..cd6a9ed27d --- /dev/null +++ b/api/nblib/listed_forces/gmxcalculator.h @@ -0,0 +1,116 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ + +#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& x, + std::vector& forces, + ListedForceCalculator::EnergyType& energies); + + [[nodiscard]] const InteractionDefinitions& getIdef() const; + +private: + int numParticles; + int numThreads; + + + std::unique_ptr idef; + std::unique_ptr ffparams; + + bonded_threading_t bondedThreading; + + std::vector shiftBuffer; + std::vector 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 lambdaBuffer; +}; + +} // namespace nblib + +#endif // NBLIB_LISTEDFORCES_GMXCALCULATOR_H diff --git a/api/nblib/listed_forces/helpers.hpp b/api/nblib/listed_forces/helpers.hpp new file mode 100644 index 0000000000..4415e9a6ca --- /dev/null +++ b/api/nblib/listed_forces/helpers.hpp @@ -0,0 +1,194 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ + +#ifndef NBLIB_LISTEDFORCSES_HELPERS_HPP +#define NBLIB_LISTEDFORCSES_HELPERS_HPP + +#include + +#include "nblib/pbc.hpp" +#include "definitions.h" +#include "nblib/util/internal.h" + +namespace nblib +{ + +namespace detail +{ +template +inline void gmxRVecZeroWorkaround([[maybe_unused]] T& value) +{ +} + +template<> +inline void gmxRVecZeroWorkaround(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 ForceBuffer +{ + using HashMap = std::unordered_map; + +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 splitListedWork(const ListedInteractionData& interactions, + int totalRange, + int nSplits) +{ + std::vector workDivision(nSplits); + + auto splitOneElement = [totalRange, nSplits, &workDivision](const auto& inputElement) { + // the index of inputElement in the ListedInteractionsTuple + constexpr int elementIndex = + FindIndex, ListedInteractionData>{}; + + // for now, copy all parameters to each split + // Todo: extract only the parameters needed for this split + for (auto& workDivisionSplit : workDivision) + { + std::get(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(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 diff --git a/api/nblib/listed_forces/kernels.hpp b/api/nblib/listed_forces/kernels.hpp new file mode 100644 index 0000000000..425820bb10 --- /dev/null +++ b/api/nblib/listed_forces/kernels.hpp @@ -0,0 +1,587 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ +#ifndef NBLIB_LISTEDFORCES_KERNELS_HPP +#define NBLIB_LISTEDFORCES_KERNELS_HPP + +#include + +#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 + */ +template +inline std::tuple 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 + */ +template +inline std::tuple 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 +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 + */ +template +inline std::tuple 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 + */ +template +inline std::tuple 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 +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 + */ +template +inline std::tuple 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 + */ +template +inline std::tuple 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 +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 + */ +template +inline std::tuple 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 +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 + */ +template +inline std::tuple 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 +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 + */ +template +inline std::tuple 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 + */ +template +inline std::tuple 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 +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 +inline auto threeCenterKernel(T dr, const DefaultAngle& angle) +{ + return harmonicScalarForce(angle.forceConstant(), angle.equilDistance(), dr); +} + + +//! \brief Computes and returns the proper dihedral force +template +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 +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 +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 +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 +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 +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 diff --git a/api/nblib/listed_forces/tests/CMakeLists.txt b/api/nblib/listed_forces/tests/CMakeLists.txt index 2f0bc8efea..af67386835 100644 --- a/api/nblib/listed_forces/tests/CMakeLists.txt +++ b/api/nblib/listed_forces/tests/CMakeLists.txt @@ -49,6 +49,12 @@ gmx_add_gtest_executable( 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) diff --git a/api/nblib/listed_forces/tests/calculator.cpp b/api/nblib/listed_forces/tests/calculator.cpp new file mode 100644 index 0000000000..677ea63d28 --- /dev/null +++ b/api/nblib/listed_forces/tests/calculator.cpp @@ -0,0 +1,266 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include "gmxpre.h" + +// includesorter/check-source want this include here. IMO a bug + + +#include + +#include + +#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 +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 bonds{ bond1, bond2 }; + // one bond between atoms 0-1 with bond1 parameters and another between atoms 1-2 with bond2 parameters + std::vector> bondIndices{ { 0, 1, 0 }, { 1, 2, 1 } }; + + DefaultAngle angle(Degrees(108.53), 397.5); + std::vector angles{ angle }; + std::vector> angleIndices{ { 0, 1, 2, 0 } }; + + pickType(interactions).indices = bondIndices; + pickType(interactions).parameters = bonds; + + pickType(interactions).indices = angleIndices; + pickType(interactions).parameters = angles; + + // initial position for the methanol atoms from the spc-water example + x = std::vector{ { 1.97, 1.46, 1.209 }, { 1.978, 1.415, 1.082 }, { 1.905, 1.46, 1.03 } }; + forces = std::vector(3, gmx::RVec{ 0, 0, 0 }); + + refBondForcesFloat = + std::valarray>{ { -22.8980637, 128.801575, 363.505951 }, + { -43.2698593, -88.0130997, -410.639252 }, + { 66.167923, -40.788475, 47.1333084 } }; + refAngleForcesFloat = + std::valarray>{ { 54.7276611, -40.1688995, 17.6805191 }, + { -81.8118973, 86.1988525, 60.1752243 }, + { 27.0842342, -46.0299492, -77.8557434 } }; + + refBondForcesDouble = std::valarray>{ + { -22.89764839974935, 128.79927224858977, 363.50016834602064 }, + { -43.24622441913251, -88.025652017772231, -410.61635172385434 }, + { 66.14387281888186, -40.773620230817542, 47.116183377833721 } + }; + refAngleForcesDouble = std::valarray>{ + { 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 x; + std::vector forces; + + ListedInteractionData interactions; + + std::shared_ptr box; + std::shared_ptr pbc; + + // reference values + std::valarray> refBondForcesFloat, refAngleForcesFloat; + std::valarray> refBondForcesDouble, refAngleForcesDouble; + + float refBondEnergyFloat, refAngleEnergyFloat; + double refBondEnergyDouble, refAngleEnergyDouble; +}; + +TEST_F(ListedExampleData, ComputeHarmonicBondForces) +{ + auto indices = pickType(interactions).indices; + auto bonds = pickType(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(interactions).indices; + auto angles = pickType(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& forces) const + { + compareVectors(forces, refForces, refForces); + } + + std::vector x; + ListedInteractionData interactions; + std::shared_ptr box; + // std::shared_ptr pbc; + +private: + std::vector refForces; + ListedForceCalculator::EnergyType refEnergies; +}; + +TEST_F(LinearChainDataFixture, Multithreading) +{ + ListedForceCalculator lfCalculator(interactions, x.size(), 4, *box); + + std::vector forces(x.size(), Vec3{ 0, 0, 0 }); + ListedForceCalculator::EnergyType energies; + lfCalculator.compute(x, forces, energies); + + testEnergies(energies); + testForces(forces); +} + + +} // namespace +} // namespace test +} // namespace nblib diff --git a/api/nblib/listed_forces/tests/conversions.cpp b/api/nblib/listed_forces/tests/conversions.cpp new file mode 100644 index 0000000000..71f8e3d077 --- /dev/null +++ b/api/nblib/listed_forces/tests/conversions.cpp @@ -0,0 +1,105 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include + +#include + +#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 bonds{ bond1, bond2 }; + pickType(interactions).parameters = bonds; + + DefaultAngle angle1(Degrees(100), 100); + DefaultAngle angle2(Degrees(101), 200); + std::vector angles{ angle1, angle2 }; + pickType(interactions).parameters = angles; + + std::vector> bondIndices{ { 0, 1, 0 }, { 1, 2, 0 }, { 2, 3, 1 } }; + pickType(interactions).indices = std::move(bondIndices); + + std::vector> angleIndices{ { 0, 1, 2, 0 }, { 1, 2, 3, 1 } }; + pickType(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(interactions).parameters[0].equilDistance()); + EXPECT_REAL_EQ_TOL(gmx_params->iparams[2].harmonic.rA, + pickType(interactions).parameters[0].equilDistance() / DEG2RAD, + gmx::test::defaultRealTolerance()); + + EXPECT_EQ(idef->il[F_BONDS].iatoms.size(), 9); + std::vector bondIatoms{ 0, 0, 1, 0, 1, 2, 1, 2, 3 }; + EXPECT_EQ(idef->il[F_BONDS].iatoms, bondIatoms); + std::vector angleIatoms{ 2, 0, 1, 2, 3, 1, 2, 3 }; + EXPECT_EQ(idef->il[F_ANGLES].iatoms, angleIatoms); + idef->clear(); +} + +} // namespace +} // namespace test +} // namespace nblib diff --git a/api/nblib/listed_forces/tests/helpers.cpp b/api/nblib/listed_forces/tests/helpers.cpp new file mode 100644 index 0000000000..e265bf7f33 --- /dev/null +++ b/api/nblib/listed_forces/tests/helpers.cpp @@ -0,0 +1,148 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include "gmxpre.h" + +#include + +#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> bondIndices{ + { 0, 1, 0 }, { 0, 6, 0 }, { 11, 12, 0 }, { 18, 19, 0 } + }; + std::vector> angleIndices{ + { 0, 1, 2, 0 }, { 0, 6, 7, 0 }, { 11, 12, 13, 0 }, { 17, 19, 18, 0 } + }; + + pickType(interactions).indices = bondIndices; + pickType(interactions).indices = angleIndices; + + std::vector splitInteractions = + splitListedWork(interactions, largestIndex, nSplits); + + std::vector> refBondIndices0{ { 0, 1, 0 }, { 0, 6, 0 } }; + std::vector> refAngleIndices0{ { 0, 1, 2, 0 }, { 0, 6, 7, 0 } }; + std::vector> refBondIndices1{ { 11, 12, 0 } }; + std::vector> refAngleIndices1{ { 11, 12, 13, 0 } }; + std::vector> refBondIndices2{ { 18, 19, 0 } }; + std::vector> refAngleIndices2{ { 17, 19, 18, 0 } }; + + EXPECT_EQ(refBondIndices0, pickType(splitInteractions[0]).indices); + EXPECT_EQ(refBondIndices1, pickType(splitInteractions[1]).indices); + EXPECT_EQ(refBondIndices2, pickType(splitInteractions[2]).indices); + + EXPECT_EQ(refAngleIndices0, pickType(splitInteractions[0]).indices); + EXPECT_EQ(refAngleIndices1, pickType(splitInteractions[1]).indices); + EXPECT_EQ(refAngleIndices2, pickType(splitInteractions[2]).indices); +} + + +TEST(NBlibTest, ListedForceBuffer) +{ + using T = gmx::RVec; + int ncoords = 20; + + T vzero{ 0, 0, 0 }; + std::vector 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 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 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 diff --git a/api/nblib/listed_forces/tests/kernels.cpp b/api/nblib/listed_forces/tests/kernels.cpp new file mode 100644 index 0000000000..317f7ba861 --- /dev/null +++ b/api/nblib/listed_forces/tests/kernels.cpp @@ -0,0 +1,64 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#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 diff --git a/api/nblib/listed_forces/tests/linear_chain_input.hpp b/api/nblib/listed_forces/tests/linear_chain_input.hpp new file mode 100644 index 0000000000..211d4edd23 --- /dev/null +++ b/api/nblib/listed_forces/tests/linear_chain_input.hpp @@ -0,0 +1,134 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ + +#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 bonds{ bond1, bond2 }; + pickType(interactions).parameters = bonds; + + DefaultAngle angle(Degrees(179.9), 397.5); + std::vector angles{ angle }; + pickType(interactions).parameters = angles; + + std::vector> bondIndices; + for (int i = 0; i < nParticles - 1; ++i) + { + // third index: alternate between the two bond parameters + bondIndices.push_back(InteractionIndex{ i, i + 1, i % 2 }); + } + pickType(interactions).indices = bondIndices; + + std::vector> angleIndices; + for (int i = 0; i < nParticles - 2; ++i) + { + angleIndices.push_back(InteractionIndex{ i, i + 1, i + 2, 0 }); + } + pickType(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(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(interactions).parameters.push_back(dummyBond); + int bondIndex = pickType(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(interactions) + .indices.push_back(InteractionIndex{ from, to, bondIndex }); + } + } + } + + int nParticles; + + std::vector x; + std::vector forces; + + ListedInteractionData interactions; + + std::shared_ptr box; +}; + +} // namespace nblib + +#endif // NBLIB_LINEAR_CHAIN_DATA_HPP + diff --git a/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_0.xml b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_0.xml new file mode 100644 index 0000000000..0c97923ca0 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_0.xml @@ -0,0 +1,28 @@ + + + + 15.644075054627638 + + 4 + + 0 + 3000.9773757382213 + 0 + + + -294.15662928338577 + 2953.652678308878 + -14.707831464169288 + + + 589.93598821678029 + -5922.0943245644257 + 29.49679941083901 + + + -295.77935893339446 + -32.535729482673389 + -14.788967946669722 + + + diff --git a/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_1.xml b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_1.xml new file mode 100644 index 0000000000..41155d21f8 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_1.xml @@ -0,0 +1,28 @@ + + + + 11.117714323462188 + + 4 + + 273.20508075688775 + 0 + 0 + + + -858.64453952164718 + 0 + 0 + + + 957.99184161506105 + 0 + 0 + + + -372.55238285030157 + 0 + 0 + + + diff --git a/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_2.xml b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_2.xml new file mode 100644 index 0000000000..3a4f134948 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/FourCenter_ListedForcesProperDihedralTest_CheckListed_2.xml @@ -0,0 +1,28 @@ + + + + 11.117714323462197 + + 4 + + 0 + 0 + 253.81712522292264 + + + 0 + 0 + -884.99823613097033 + + + 0 + 0 + 299.4154280846883 + + + 0 + 0 + 331.76568282335938 + + + diff --git a/api/nblib/listed_forces/tests/transformations.cpp b/api/nblib/listed_forces/tests/transformations.cpp new file mode 100644 index 0000000000..1e02da05f2 --- /dev/null +++ b/api/nblib/listed_forces/tests/transformations.cpp @@ -0,0 +1,94 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include + +#include + +#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> bondIndices{ { 0, 2, 0 }, { 0, 1, 0 } }; + pickType(interactions).indices = std::move(bondIndices); + + std::vector> angleIndices{ { 0, 1, 2, 0 }, { 1, 0, 2, 0 } }; + pickType(interactions).indices = std::move(angleIndices); + + std::vector> dihedralIndices{ { 0, 2, 1, 3, 0 }, { 0, 1, 2, 3, 0 } }; + pickType(interactions).indices = std::move(dihedralIndices); + + return interactions; +} + +TEST(ListedTransformations, SortInteractionIndices) +{ + ListedInteractionData interactions = unsortedInteractions(); + sortInteractions(interactions); + + std::vector> refBondIndices{ { 0, 1, 0 }, { 0, 2, 0 } }; + std::vector> refAngleIndices{ { 1, 0, 2, 0 }, { 0, 1, 2, 0 } }; + std::vector> refDihedralIndices{ { 0, 1, 2, 3, 0 }, + { 0, 2, 1, 3, 0 } }; + + EXPECT_EQ(pickType(interactions).indices, refBondIndices); + EXPECT_EQ(pickType(interactions).indices, refAngleIndices); + EXPECT_EQ(pickType(interactions).indices, refDihedralIndices); +} + +} // namespace +} // namespace test +} // namespace nblib diff --git a/api/nblib/listed_forces/tests/typetests.cpp b/api/nblib/listed_forces/tests/typetests.cpp new file mode 100644 index 0000000000..75a094b948 --- /dev/null +++ b/api/nblib/listed_forces/tests/typetests.cpp @@ -0,0 +1,123 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#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> 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> 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> c_InputDihs = { { RyckaertBellemanDihedral({ -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 }) } }; + +template +class ListedForcesBase +{ +public: + std::vector input_; + std::vector x_; + std::vector> indices_; + PbcHolder pbcHolder_; + gmx::test::TestReferenceData refData_; + gmx::test::TestReferenceChecker checker_; + std::vector forces_; + real energy_; + + ListedForcesBase(std::vector input, + std::vector coordinates, + std::vector> 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, + public testing::TestWithParam, std::vector>> +{ + using Base = ListedForcesBase; + +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 diff --git a/api/nblib/listed_forces/transformations.cpp b/api/nblib/listed_forces/transformations.cpp new file mode 100644 index 0000000000..fc03957593 --- /dev/null +++ b/api/nblib/listed_forces/transformations.cpp @@ -0,0 +1,64 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ +#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 diff --git a/api/nblib/listed_forces/transformations.h b/api/nblib/listed_forces/transformations.h new file mode 100644 index 0000000000..ac1471d59c --- /dev/null +++ b/api/nblib/listed_forces/transformations.h @@ -0,0 +1,148 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ +#ifndef NBLIB_LISTEDFORCES_TRANSFORMATIONS_H +#define NBLIB_LISTEDFORCES_TRANSFORMATIONS_H + +#include + +#include "definitions.h" + +namespace nblib +{ + +namespace detail +{ + +/*! \brief tuple ordering for two center interactions + * + * \param t input array + * \return ordered array + */ +inline std::array nblibOrdering(const std::array& 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{ id1, id2 }; +} + +/*! \brief tuple ordering for three center interactions + * + * \param t input array + * \return ordered array + */ +inline std::array nblibOrdering(const std::array& 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{ id1, std::get<1>(t), id3 }; +} + +/*! \brief tuple ordering for four center interactions + * + * \param t input array + * \return ordered array + */ +inline std::array nblibOrdering(const std::array& t) +{ + return t; +} + +/*! \brief tuple ordering for five center interactions + * + * \param t input array + * \return ordered array + */ +inline std::array nblibOrdering(const std::array& 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 diff --git a/api/nblib/molecules.cpp b/api/nblib/molecules.cpp index c2491e7b55..d9e672c0a9 100644 --- a/api/nblib/molecules.cpp +++ b/api/nblib/molecules.cpp @@ -115,6 +115,98 @@ Molecule& Molecule::addParticle(const ParticleName& particleName, const Particle return *this; } +//! Two-particle interactions such as bonds and LJ1-4 +template +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(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 +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 +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(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 +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(); @@ -146,6 +238,11 @@ void Molecule::addExclusion(const ParticleName& particleName, const ParticleName 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); diff --git a/api/nblib/pbc.hpp b/api/nblib/pbc.hpp new file mode 100644 index 0000000000..95c0c43db3 --- /dev/null +++ b/api/nblib/pbc.hpp @@ -0,0 +1,95 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + * \author Artem Zhmurov + */ +#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 diff --git a/api/nblib/tests/CMakeLists.txt b/api/nblib/tests/CMakeLists.txt index 366027339c..4707595b4d 100644 --- a/api/nblib/tests/CMakeLists.txt +++ b/api/nblib/tests/CMakeLists.txt @@ -57,6 +57,7 @@ gmx_add_gtest_executable( box.cpp interactions.cpp particletype.cpp + pbcholder.cpp molecules.cpp topology.cpp ) diff --git a/api/nblib/tests/molecules.cpp b/api/nblib/tests/molecules.cpp index b906e097b5..009f18f246 100644 --- a/api/nblib/tests/molecules.cpp +++ b/api/nblib/tests/molecules.cpp @@ -206,6 +206,34 @@ TEST(NBlibTest, MoleculeNoThrowsSameParticleTypeName) 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(interactionData).interactions_.size(), 2); + //! cubic bonds + EXPECT_EQ(pickType(interactionData).interactions_.size(), 1); + //! angular interactions + EXPECT_EQ(pickType(interactionData).interactions_.size(), 1); +} + } // namespace } // namespace test } // namespace nblib diff --git a/api/nblib/tests/pbcholder.cpp b/api/nblib/tests/pbcholder.cpp new file mode 100644 index 0000000000..af53143aaa --- /dev/null +++ b/api/nblib/tests/pbcholder.cpp @@ -0,0 +1,73 @@ +/* + * 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 + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include + +#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 diff --git a/api/nblib/tests/testsystems.cpp b/api/nblib/tests/testsystems.cpp index 41045b08c4..cb8f8bc236 100644 --- a/api/nblib/tests/testsystems.cpp +++ b/api/nblib/tests/testsystems.cpp @@ -118,6 +118,10 @@ WaterMoleculeBuilder::WaterMoleculeBuilder() : water_(MoleculeName("SOL")) 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() @@ -151,6 +155,15 @@ MethanolMoleculeBuilder::MethanolMoleculeBuilder() : methanol_(MoleculeName("MeO 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() diff --git a/api/nblib/tests/topology.cpp b/api/nblib/tests/topology.cpp index 63c326afd4..bf9bc58727 100644 --- a/api/nblib/tests/topology.cpp +++ b/api/nblib/tests/topology.cpp @@ -192,6 +192,255 @@ TEST(NBlibTest, TopologyHasSequencing) 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> molecules{ std::make_tuple(water, 2), + std::make_tuple(methanol, 1) }; + std::vector bonds; + std::vector bondsExpansion; + std::tie(bondsExpansion, bonds) = detail::collectInteractions(molecules); + + std::vector 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 waterBonds = + pickType(water.interactionData()).interactionTypes_; + std::vector methanolBonds = + pickType(methanol.interactionData()).interactionTypes_; + + std::vector 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> molecules{ std::make_tuple(water, 2), + std::make_tuple(methanol, 1) }; + detail::ParticleSequencer particleSequencer; + particleSequencer.build(molecules); + auto pairs = detail::sequenceIDs(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> 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> molecules{ std::make_tuple(water, 2), + std::make_tuple(methanol, 1) }; + detail::ParticleSequencer particleSequencer; + particleSequencer.build(molecules); + auto pairs = detail::sequenceIDs(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 bonds{ b2, b2, b3, b1, b2, b1, b3, b3 }; + + // expected output + std::vector uniqueBondsReference{ b1, b2, b3 }; + std::vector indicesReference{ 1, 1, 2, 0, 1, 0, 2, 2 }; + + std::tuple, std::vector> 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(interactionData); + + auto& indices = harmonicBonds.indices; + auto& bonds = harmonicBonds.parameters; + + std::map, 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, 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 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(interactionData); + auto& cubicBonds = pickType(interactionData); + auto& angleInteractions = pickType(interactionData); + + HarmonicBondType ohBond(1., 1.); + HarmonicBondType ohBondMethanol(1.01, 1.02); + HarmonicBondType ometBond(1.1, 1.2); + std::vector 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 cubicBondsReference{ testBond }; + std::vector> 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 angleReference{ testAngle, methanolAngle }; + std::vector> 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 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> testInput{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 }, diff --git a/api/nblib/topology.cpp b/api/nblib/topology.cpp index 4e2e1e4d8c..e6b6acd07b 100644 --- a/api/nblib/topology.cpp +++ b/api/nblib/topology.cpp @@ -101,6 +101,59 @@ gmx::ListOfLists TopologyBuilder::createExclusionsListOfLists() const 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(interactionData)); + + auto create = [this, &particleSequencer](auto& interactionDataElement) { + using InteractionType = typename std::decay_t::type; + + // first compression stage: each bond per molecule listed once, + // eliminates duplicates from multiple identical molecules + auto compressedDataStage1 = detail::collectInteractions(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 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 + auto coordinateIndices = detail::sequenceIDs(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 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 std::vector TopologyBuilder::extractParticleTypeQuantity(Extractor&& extractor) { @@ -164,6 +217,8 @@ Topology TopologyBuilder::buildTopology() 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_) @@ -252,6 +307,11 @@ NonBondedInteractionMap Topology::getNonBondedInteractionMap() const return nonBondedInteractionMap_; } +ListedInteractionData Topology::getInteractionData() const +{ + return interactionData_; +} + CombinationRule Topology::getCombinationRule() const { return combinationRule_; diff --git a/api/nblib/topologyhelpers.cpp b/api/nblib/topologyhelpers.cpp index d2cdd63d66..47dee71d09 100644 --- a/api/nblib/topologyhelpers.cpp +++ b/api/nblib/topologyhelpers.cpp @@ -47,6 +47,7 @@ #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" @@ -156,6 +157,148 @@ void ParticleSequencer::build(const std::vector>& mole } } +template +std::tuple, std::vector> +collectInteractions(const std::vector>& molecules) +{ + std::vector collectedBonds; + std::vector expansionArray; + for (auto& molNumberTuple : molecules) + { + const Molecule& molecule = std::get<0>(molNumberTuple); + size_t numMols = std::get<1>(molNumberTuple); + + auto& interactions = pickType(molecule.interactionData()).interactionTypes_; + + std::vector 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> collectInteractions( \ + const std::vector>&); +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 +auto stringsToIndices_impl(const Tuple& tuple, std::index_sequence 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{ 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 + * + */ +template +auto stringsToIndices(const Tuple& tuple, F&& f, Args... args) +{ + auto is = std::make_index_sequence::value / 2>{}; + return stringsToIndices_impl(tuple, is, std::forward(f), args...); +} + +//! Tuple ordering for two center interactions +[[maybe_unused]] static std::array nblibOrdering(const std::array& 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{ id1, id2 }; +} + +//! Tuple ordering for three center interactions +[[maybe_unused]] static std::array nblibOrdering(const std::array& 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{ id1, std::get<1>(t), id3 }; +} + +//! Tuple ordering for four center interactions +[[maybe_unused]] static std::array nblibOrdering(const std::array& t) +{ + return t; +} + +//! Tuple ordering for five center interactions +[[maybe_unused]] static std::array nblibOrdering(const std::array& t) +{ + return t; +} + +} // namespace sequence_detail + +//! For each interaction, translate particle identifiers (moleculeName, nr, residueName, +//! particleName) to particle coordinate indices +template +std::vector> sequenceIDs(const std::vector>& molecules, + const detail::ParticleSequencer& particleSequencer) +{ + std::vector> 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(molecule.interactionData()).interactions_; + for (const auto& interactionString : interactions) + { + CoordinateIndex 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> sequenceIDs( \ + const std::vector>&, const detail::ParticleSequencer&); +MAP(SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES) +#undef SEQUENCE_PAIR_ID_INSTANTIATE_TEMPLATE +/// \endcond template std::tuple, std::vector> eliminateDuplicateInteractions(const std::vector& aggregatedInteractions) @@ -211,6 +354,14 @@ std::tuple, std::vector> eliminateDuplicateInteractions(c return make_tuple(uniqueIndices, uniquInteractionsInstances); } +/// \cond DO_NOT_DOCUMENT +#define ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE(x) \ + template std::tuple, std::vector> eliminateDuplicateInteractions( \ + const std::vector& aggregatedBonds); +MAP(ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE, SUPPORTED_LISTED_TYPES) +#undef ELIMINATE_DUPLICATE_BONDS_INSTANTIATE_TEMPLATE +/// \endcond + } // namespace detail } // namespace nblib diff --git a/api/nblib/util/user.h b/api/nblib/util/user.h index 6d782fb6bd..148e114051 100644 --- a/api/nblib/util/user.h +++ b/api/nblib/util/user.h @@ -55,6 +55,7 @@ #include #include "nblib/basicdefinitions.h" +#include "nblib/ppmap.h" #include "nblib/vector.h" namespace gmx diff --git a/docs/nblib/listed-dev.rst b/docs/nblib/listed-dev.rst index ab8f18cfe8..2fd7490725 100644 --- a/docs/nblib/listed-dev.rst +++ b/docs/nblib/listed-dev.rst @@ -2,7 +2,7 @@ Adding New Listed-Interaction Types in NB-LIB ============================================= NB-LIB currently has code paths for listed interactions that occur between two, three, four and five different particles. -To extend NB-LIB to support more types of particle interactions, modify the following three files. +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)``. @@ -20,11 +20,40 @@ Accepting these constraints, it is possible to add a new kernel by modifying the 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 + { + NewBondType() = default; + NewBondType(ForceConstant f, EquilDistance d, ScaleFactor s) : + std::tuple{ f, d, s } + { + } + + NAMED_MEMBERS(forceConstant, equilDistance, scaleFactor) + }; .. _definitions.h: @@ -46,16 +75,10 @@ Assuming that the only other two center interaction is called ``DefaultBond``, t 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 @@ -91,7 +114,7 @@ The kernel return type is always an ``std::tuple`` of the force and the potentia } template - inline std::tuple bondKernel(T dr, const NewBondType& bond) + inline auto bondKernel(T dr, const NewBondType& bond) { return newBondForce(bond.forceConstant(), bond.equilDistance(), bond.scaleFactor(), dr); } -- 2.22.0