docs/nblib/listed-dev.rst

   1 Adding New Listed-Interaction Types in NB-LIB
   2 =============================================
   3
   4 NB-LIB currently has code paths for listed interactions that occur between two, three, four and five different particles.
   5 It is easy to extend NB-LIB to support novel formulations of particle interactions by modifying the following three files.
   6
   7 Two center interactions must use the distance between the centers as an input to the force kernel.
   8 Three center interactions take the form ``(particleI, particleJ, ParticleK)``.
   9 In this case, the middle particle, ``particleJ`` is taken as the center around which the angle is computed.
  10 This angle must be an input to a three center force kernel.
  11 Likewise for four center interactions, the dihedral angle phi must be an input to the force kernel.
  12 Accepting these constraints, it is possible to add a new kernel by modifying the following three files.
  13
  14 1) bondtypes.h_
  15 2) definitions.h_
  16 3) kernels.hpp_
  17
  18 .. _bondtypes.h:
  19
  20 1) bondtypes.h
  21 ---------------
  22
  23 This file contains one ``struct`` for each interaction type parameter set.
  24 New interaction types are added here as separate structs. There
  25 are no content requirements, but for convenience, the existing ``NAMED_MEBERS``
  26 macro in combination with inheriting from a ``std::tuple`` or ``std::array``
  27 may be used. The macro can be used to define the
  28 parameter names for the corresponding setters and getters.
  29 For example, ``NAMED_MEMBERS(forceConstant, equilDistance)`` will expand to
  30
  31 .. code:: cpp
  32
  33    inline auto& forceConstant() { return std::get<0>(*this); }
  34    inline auto& equilDistance() { return std::get<1>(*this); }
  35    inline const auto& forceConstant() const { return std::get<0>(*this); }
  36    inline const auto& equilDistance() const { return std::get<1>(*this); }
  37
  38 Putting everything together, one could define the complete parameter set for a new interaction type as follows.
  39
  40 .. code:: cpp
  41
  42    /*! \brief new bond type
  43     *
  44     * V(r; forceConstant, equilDistance, scaleFactor)
  45     *      = forceConstant * exp( (r - equilDistance) / scaleFactor)
  46     */
  47    struct NewBondType : public std::tuple<real, real, int>
  48    {
  49        NewBondType() = default;
  50        NewBondType(ForceConstant f, EquilDistance d, ScaleFactor s) :
  51            std::tuple<real, real, int>{ f, d, s }
  52        {
  53        }
  54
  55        NAMED_MEMBERS(forceConstant, equilDistance, scaleFactor)
  56    };
  57
  58 .. _definitions.h:
  59
  60 2) definitions.h
  61 ------------------------
  62
  63 This file begins with pre-processor macro lists that classify concrete interaction types into two, three, four and five center types.
  64 To add a new type, the user must add the interaction type parameter struct name to the macro of the correct center number.
  65 In this case, ``NewBondType`` is an example of a two center interaction.
  66 As such it would get added to the ``SUPPORTED_TWO_CENTER_TYPES`` macro.
  67 Assuming that the only other two center interaction is called ``DefaultBond``, the result would look like the following snippet.
  68
  69 .. code:: cpp
  70
  71     #define SUPPORTED_TWO_CENTER_TYPES DefaultBond, NewBondType
  72
  73 .. _kernels.hpp:
  74
  75 Adding ``NewBondType`` to this macro ensures that the NBLIB ``molecule``
  76 class ``addInteraction`` function supports adding the new bond type
  77 and includes it in the listed interaction data that the ``topology`` class
  78 provides.
  79
  80 Note that, as of C++17, there's no alternative to preprocessor macros for adding
  81 the required template instantiations controlled through the macros described here.
  82 In NBLIB, the design decision we took, was that we did not want to expose a templated
  83 interface in a user header and it is for this reason that we explicitly need
  84 to instantiate the interface with all the supported listed interaction types defined
  85 in this macro.
  86
  87 3) kernels.hpp
  88 ---------------------
  89
  90 In this file the actual force kernels for each interaction type are implemented.
  91 Each kernel call is templated to allow various precisions and is
  92 accessed through an overload ``bondKernel`` that extracts the relevant
  93 parameters from a ``const NewBondType&`` argument.
  94 The kernel return type is always an ``std::tuple`` of the force and the potential.
  95
  96 .. code:: cpp
  97
  98    /*! \brief kernel to calculate the new bond type force
  99     *
 100     * \param k     Force constant
 101     * \param x0    Equilibrium distance
 102     * \param scale The scaling factor
 103     * \param x     Input bond length
 104     *
 105     * \return tuple<force, potential energy>
 106     */
 107    template <class T>
 108    std::tuple<T, T> newBondForce(T k, T x0, T scale, T x)
 109    {
 110        real exponent = std::exp( (x - x0) / scale);
 111        real epot = k * exponent;
 112        real force =  epot / scale;
 113        return std::make_tuple(force, epot);
 114    }
 115
 116   template <class T>
 117   inline auto bondKernel(T dr, const NewBondType& bond)
 118   {
 119       return newBondForce(bond.forceConstant(), bond.equilDistance(), bond.scaleFactor(), dr);
 120   }
 121