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 To extend NB-LIB to support more types of particle interactions, modify 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 C++ type to store the parameters for each interaction type.
  24 New interaction types are added here as separate C++ types.
  25 The interface of these types is completely unrestricted.
  26 The only requirements are equality and less than comparison, and that the interface be
  27 compatible with the corresponding (user-added) kernel.
  28
  29 .. _definitions.h:
  30
  31 2) definitions.h
  32 ------------------------
  33
  34 This file begins with pre-processor macro lists that classify concrete interaction types into two, three, four and five center types.
  35 To add a new type, the user must add the interaction type parameter struct name to the macro of the correct center number.
  36 In this case, ``NewBondType`` is an example of a two center interaction.
  37 As such it would get added to the ``SUPPORTED_TWO_CENTER_TYPES`` macro.
  38 Assuming that the only other two center interaction is called ``DefaultBond``, the result would look like the following snippet.
  39
  40 .. code:: cpp
  41
  42     #define SUPPORTED_TWO_CENTER_TYPES DefaultBond, NewBondType
  43
  44 .. _kernels.hpp:
  45
  46 Adding ``NewBondType`` to this macro ensures that the NBLIB ``molecule``
  47 class ``addInteraction`` function supports adding the new bond type
  48 and includes it in the listed interaction data that the ``topology`` class
  49 provides. The ``SUPPORTED_TWO_CENTER_TYPES`` macro is immediately converted into a
  50 C++ type list that is implemented as a variadic template. The type list
  51 is then used to define all the dependent data structures. Apart from creating
  52 the type list, the only place where the macro is needed is explicit template instantiation.
  53
  54 Note that, as of C++17, there's no alternative to preprocessor macros for adding
  55 the required template instantiations controlled through the macros described here.
  56 (Other than manually adding the template instantiations, which would require the instantiation list
  57 of several templates to be updated each time a new interaction type is added. Compared to the preprocessor
  58 based solution where just a single macro has to be extended, this would clearly be an inferior solution.)
  59 In NBLIB, the design decision we took, was that we did not want to expose a templated
  60 interface in a user header and it is for this reason that we explicitly need
  61 to instantiate the interface with all the supported listed interaction types defined
  62 in this macro.
  63
  64 3) kernels.hpp
  65 ---------------------
  66
  67 In this file the actual force kernels for each interaction type are implemented.
  68 Each kernel call is templated to allow various precisions and is
  69 accessed through an overload ``bondKernel`` that extracts the relevant
  70 parameters from a ``const NewBondType&`` argument.
  71 The kernel return type is always an ``std::tuple`` of the force and the potential.
  72
  73 .. code:: cpp
  74
  75    /*! \brief kernel to calculate the new bond type force
  76     *
  77     * \param k     Force constant
  78     * \param x0    Equilibrium distance
  79     * \param scale The scaling factor
  80     * \param x     Input bond length
  81     *
  82     * \return tuple<force, potential energy>
  83     */
  84    template <class T>
  85    std::tuple<T, T> newBondForce(T k, T x0, T scale, T x)
  86    {
  87        real exponent = std::exp( (x - x0) / scale);
  88        real epot = k * exponent;
  89        real force =  epot / scale;
  90        return std::make_tuple(force, epot);
  91    }
  92
  93   template <class T>
  94   inline std::tuple<T, T> bondKernel(T dr, const NewBondType& bond)
  95   {
  96       return newBondForce(bond.forceConstant(), bond.equilDistance(), bond.scaleFactor(), dr);
  97   }
  98