api/nblib/listed_forces/conversions.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * This implements conversion utilities between the internal
  38  * representations of the listed forces parameters for NBLIB
  39  * and that of the GROMACS backend
  40  *
  41  * \author Victor Holanda <victor.holanda@cscs.ch>
  42  * \author Joe Jordan <ejjordan@kth.se>
  43  * \author Prashanth Kanduri <kanduri@cscs.ch>
  44  * \author Sebastian Keller <keller@cscs.ch>
  45  */
  46
  47 #ifndef NBLIB_LISTEDFORCES_CONVERSION_HPP
  48 #define NBLIB_LISTEDFORCES_CONVERSION_HPP
  49
  50 #include <memory>
  51
  52 #include "gromacs/topology/forcefieldparameters.h"
  53 #include "gromacs/topology/idef.h"
  54 #include "nblib/listed_forces/traits.h"
  55
  56 namespace nblib
  57 {
  58
  59 /*! this trait maps nblib InteractionTypes to the corresponding gmx enum
  60  *
  61  * \tparam InteractionType
  62  */
  63 template<class InteractionType>
  64 struct ListedIndex
  65 {
  66 };
  67
  68 template<>
  69 struct ListedIndex<HarmonicBondType> : std::integral_constant<int, F_BONDS>
  70 {
  71 };
  72
  73 template<>
  74 struct ListedIndex<MorseBondType> : std::integral_constant<int, F_MORSE>
  75 {
  76 };
  77
  78 template<>
  79 struct ListedIndex<FENEBondType> : std::integral_constant<int, F_FENEBONDS>
  80 {
  81 };
  82
  83 template<>
  84 struct ListedIndex<CubicBondType> : std::integral_constant<int, F_CUBICBONDS>
  85 {
  86 };
  87
  88 template<>
  89 struct ListedIndex<G96BondType> : std::integral_constant<int, F_G96BONDS>
  90 {
  91 };
  92
  93 template<>
  94 struct ListedIndex<HarmonicAngle> : std::integral_constant<int, F_ANGLES>
  95 {
  96 };
  97
  98 template<>
  99 struct ListedIndex<ProperDihedral> : std::integral_constant<int, F_PDIHS>
 100 {
 101 };
 102
 103 template<>
 104 struct ListedIndex<ImproperDihedral> : std::integral_constant<int, F_IDIHS>
 105 {
 106 };
 107
 108 template<>
 109 struct ListedIndex<RyckaertBellemanDihedral> : std::integral_constant<int, F_RBDIHS>
 110 {
 111 };
 112
 113 template<class InteractionType>
 114 struct ListedTypeIsImplemented : std::bool_constant<detail::HasValueMember<ListedIndex<InteractionType>>::value>
 115 {
 116 };
 117
 118 namespace detail
 119 {
 120
 121 template<class InteractionData>
 122 void transferParameters([[maybe_unused]] const InteractionData& interactionData,
 123                         [[maybe_unused]] gmx_ffparams_t&        gmx_params)
 124 {
 125 }
 126
 127 //! \brief Harmonic bond parameter conversion function
 128 template<>
 129 void transferParameters(const ListedTypeData<HarmonicBondType>& interactions, gmx_ffparams_t& gmx_params)
 130 {
 131     for (const auto& hbond : interactions.parameters)
 132     {
 133         t_iparams param;
 134         param.harmonic.krA = hbond.forceConstant();
 135         param.harmonic.rA  = hbond.equilConstant();
 136         gmx_params.iparams.push_back(param);
 137     }
 138 }
 139
 140 //! \brief Harmonic angle parameter conversion function
 141 template<>
 142 void transferParameters(const ListedTypeData<HarmonicAngle>& interactions, gmx_ffparams_t& gmx_params)
 143 {
 144     for (const auto& angle : interactions.parameters)
 145     {
 146         t_iparams param;
 147         param.harmonic.krA = angle.forceConstant();
 148         param.harmonic.rA  = angle.equilConstant() / DEG2RAD;
 149         gmx_params.iparams.push_back(param);
 150     }
 151 }
 152
 153 //! \brief Proper dihedral parameter conversion function
 154 template<>
 155 void transferParameters(const ListedTypeData<ProperDihedral>& interactions, gmx_ffparams_t& gmx_params)
 156 {
 157     for (const auto& dihedral : interactions.parameters)
 158     {
 159         t_iparams param;
 160         param.pdihs.phiA = dihedral.equilDistance() / DEG2RAD;
 161         param.pdihs.cpA  = dihedral.forceConstant();
 162         param.pdihs.mult = dihedral.multiplicity();
 163         gmx_params.iparams.push_back(param);
 164     }
 165 }
 166
 167 template<class TwoCenterType>
 168 std::enable_if_t<Contains<TwoCenterType, SupportedTwoCenterTypes>{}>
 169 transferIndicesImpl(const ListedTypeData<TwoCenterType>& interactions, InteractionDefinitions& idef, int offset)
 170 {
 171     for (const auto& index : interactions.indices)
 172     {
 173         int parameterIndex = index[2] + offset;
 174         idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(parameterIndex);
 175         idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(index[0]);
 176         idef.il[ListedIndex<TwoCenterType>::value].iatoms.push_back(index[1]);
 177     }
 178 }
 179
 180 template<class ThreeCenterType>
 181 std::enable_if_t<Contains<ThreeCenterType, SupportedThreeCenterTypes>{}>
 182 transferIndicesImpl(const ListedTypeData<ThreeCenterType>& interactions, InteractionDefinitions& idef, int offset)
 183 {
 184     for (const auto& index : interactions.indices)
 185     {
 186         int parameterIndex = index[3] + offset;
 187         idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(parameterIndex);
 188         idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[0]);
 189         idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[1]);
 190         idef.il[ListedIndex<ThreeCenterType>::value].iatoms.push_back(index[2]);
 191     }
 192 }
 193
 194 template<class FourCenterType>
 195 std::enable_if_t<Contains<FourCenterType, SupportedFourCenterTypes>{}>
 196 transferIndicesImpl(const ListedTypeData<FourCenterType>& interactions, InteractionDefinitions& idef, int offset)
 197 {
 198     for (const auto& index : interactions.indices)
 199     {
 200         int parameterIndex = index[4] + offset;
 201         idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(parameterIndex);
 202         idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[0]);
 203         idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[1]);
 204         idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[2]);
 205         idef.il[ListedIndex<FourCenterType>::value].iatoms.push_back(index[3]);
 206     }
 207 }
 208
 209 template<class FiveCenterType>
 210 std::enable_if_t<Contains<FiveCenterType, SupportedFiveCenterTypes>{}>
 211 transferIndicesImpl(const ListedTypeData<FiveCenterType>& interactions, InteractionDefinitions& idef, int offset)
 212 {
 213     for (const auto& index : interactions.indices)
 214     {
 215         int parameterIndex = index[5] + offset;
 216         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(parameterIndex);
 217         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[0]);
 218         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[1]);
 219         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[2]);
 220         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[3]);
 221         idef.il[ListedIndex<FiveCenterType>::value].iatoms.push_back(index[4]);
 222     }
 223 }
 224
 225 template<template<class> class Container, class InteractionType>
 226 void transferIndices(const Container<InteractionType>&  interactionData,
 227                      InteractionDefinitions& idef,
 228                      [[maybe_unused]] int offset)
 229 {
 230     if constexpr (ListedTypeIsImplemented<InteractionType>{})
 231     {
 232         transferIndicesImpl(interactionData, idef, offset);
 233     }
 234 }
 235
 236 } // namespace detail
 237
 238 /*! \brief function to convert from nblib-ListedInteractionData to gmx-InteractionDefinitions
 239  *
 240  * Currently only supports harmonic bonds and angles, other types are ignored
 241  *
 242  * \param interactions
 243  * \return
 244  */
 245 std::tuple<std::unique_ptr<InteractionDefinitions>, std::unique_ptr<gmx_ffparams_t>>
 246 createFFparams(const ListedInteractionData& interactions);
 247
 248 std::tuple<std::unique_ptr<InteractionDefinitions>, std::unique_ptr<gmx_ffparams_t>>
 249 createFFparams(const ListedInteractionData& interactions)
 250 {
 251     std::unique_ptr<gmx_ffparams_t> ffparamsHolder = std::make_unique<gmx_ffparams_t>();
 252     std::unique_ptr<InteractionDefinitions> idefHolder = std::make_unique<InteractionDefinitions>(*ffparamsHolder);
 253
 254     gmx_ffparams_t& ffparams = *ffparamsHolder;
 255     InteractionDefinitions& idef = *idefHolder;
 256
 257     auto copyParamsOneType = [&ffparams](const auto& interactionElement)
 258     {
 259         detail::transferParameters(interactionElement, ffparams);
 260     };
 261     for_each_tuple(copyParamsOneType, interactions);
 262
 263     // since gmx_ffparams_t.iparams is a flattened vector over all interaction types,
 264     // we need to compute offsets for each type to know where the parameters for each type start
 265     // in the flattened iparams vectors
 266     int acc = 0;
 267     std::array<int, std::tuple_size_v<ListedInteractionData>> indexOffsets{0};
 268     auto extractNIndices = [&indexOffsets, &acc](const auto& interactionElement)
 269     {
 270         constexpr int elementIndex = FindIndex<std::decay_t<decltype(interactionElement)>, ListedInteractionData>::value;
 271         indexOffsets[elementIndex] = acc;
 272         acc += interactionElement.parameters.size();
 273     };
 274     for_each_tuple(extractNIndices, interactions);
 275
 276     auto copyIndicesOneType = [&idef, &indexOffsets](const auto& interactionElement)
 277     {
 278         constexpr int elementIndex = FindIndex<std::decay_t<decltype(interactionElement)>, ListedInteractionData>::value;
 279         detail::transferIndices(interactionElement, idef, indexOffsets[elementIndex]);
 280     };
 281     for_each_tuple(copyIndicesOneType, interactions);
 282
 283     return std::make_tuple(std::move(idefHolder), std::move(ffparamsHolder));
 284 }
 285
 286
 287 } // namespace nblib
 288
 289 #endif // NBLIB_LISTEDFORCES_CONVERSION_HPP