bonded_interactions.dox

/*
 * Copyright (C) 2019 The ESPResSo project
 *
 * This file is part of ESPResSo.
 *
 * ESPResSo is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * ESPResSo 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/** @page bondedIA Bonded interactions
 *
 *  @tableofcontents
 *
 *  @section bondedIA_new Adding new interactions
 *
 *  To add a new bonded interaction, the following steps are required:
 *  * core:
 *   - define a data structure for the interaction parameters (prefactors, etc.)
 *   - write a setter function, which takes the parameters of the interactions
 *     and stores them in the interaction data structure
 *   - write functions for calculating force and energy
 *   - add calls to the force and energy calculation functions to the force
 *     calculation in the integration loop as well as to energy and
 *     pressure tensor analysis
 *  * Python interface:
 *   - import the definition of the interaction struct from the core
 *   - implement a class for the bonded interaction derived from the Python
 *     BondedInteraction base class
 *
 *  @subsection bondedIA_new_struct Defining the interaction data structure
 *
 *  The data structures for bonded interactions reside in
 *  bonded_interaction_data.hpp.
 *
 *  * Add your interaction to the @ref BondedInteraction enum, it is used to
 *    identify different bonded interactions.
 *  * Create a new struct containing the parameters of the interaction, using
 *    for example @ref Fene_bond_parameters as a guide.
 *  * Add the struct to the @ref Bond_parameters union
 *
 *  @subsection bondedIA_new_fun Functions for calculating force and energy, \
 *                               and for setting parameters
 *
 *  Every interaction resides in its own source .cpp and .hpp files. A simple
 *  example for a bonded interaction is the FENE bond in @ref fene.hpp and
 *  @ref fene.cpp. Use these two files as templates for your interaction.
 *  In particular, use the same function names: \c NAME_pair_force,
 *  \c NAME_pair_energy, \c NAME_set_params, \c angle_NAME_force,
 *  \c angle_NAME_energy, \c angle_NAME_set_params, etc. with \c NAME your
 *  interaction name. In the checklist below, we'll use \c NAME = \c fene.
 *
 *  * The recommended signatures of the force and energy functions are:
 *    @code{.cpp}
 *    inline boost::optional<Utils::Vector3d>
 *    fene_pair_force(Bonded_ia_parameters const &iaparams,
 *                    Utils::Vector3d const &dx);
 *    inline boost::optional<double>
 *    fene_pair_energy(Bonded_ia_parameters const &iaparams,
 *                     Utils::Vector3d const &dx);
 *    @endcode
 *  * The setter function gets a \p bond_type which is a numerical id
 *    identifying the number of the bond type in the simulation.
 *    It DOES NOT determine the type of the bond potential (harmonic vs FENE).
 *    The signature of the setter function has to contain the \c bond_type, the
 *    remaining parameters are specific to the interaction:
 *    @code{.cpp}
 *    int fene_set_params(int bond_type, double k, double drmax, double r0);
 *    @endcode
 *    @ref ES_OK is returned on success, @ref ES_ERROR on error, e.g. when
 *    parameters are invalid.
 *  * The setter function must call @ref make_bond_type_exist() with
 *    \p bond_type to allocate the memory for storing the parameters.
 *  * Afterwards, the bond parameters can be stored in the global variable
 *    @ref bonded_ia_params "bonded_ia_params[bond_type]"
 *  * \c bonded_ia_params[bond_type].num is the number of particles involved
 *    in the bond, minus 1, i.e. 1 for a pairwise bonded potential such as the
 *    FENE bond.
 *  * The parameters for the individual bonded interaction go to the member of
 *    @ref Bond_parameters for your interaction defined in the previous step.
 *    For the FENE bond, this would be:
 *    @code{.cpp}
 *    bonded_ia_params[bond_type].p.fene
 *    @endcode
 *  * At the end of the parameter setter function, do not forget the call to
 *    @ref mpi_bcast_ia_params(), which will sync the parameters just set to
 *    other compute nodes in a parallel simulation.
 *  * The functions calculating force(s) and energy return their values in a
 *    \c boost::optional container if the bond is breakable. If the bond is
 *    broken, the returned object is empty; this will stop the integrator
 *    with a runtime error.
 *  * The functions calculating force and energy can make use of a
 *    pre-calculated distance vector (\p dx) pointing from particle 2 to
 *    particle 1, that takes periodic boundary conditions into account.
 *  * The force and energy functions cannot alter the particle states.
 *
 *  @subsection bondedIA_new_integration Including the bonded interaction in \
 *                                       the force calculation and the energy \
 *                                       and pressure analysis
 *
 *  * In bonded_interaction_data.cpp:
 *    - in function @ref maximal_cutoff_bonded(): add a case for the new
 *      interaction which makes sure that the cutoff is as large as
 *      the interaction range of the new interaction. This is only relevant to
 *      pairwise bonds and dihedral bonds. The range can be calculated by a
 *      formula, so it is not strictly necessary that the maximal interaction
 *      range be stored explicitly in the struct.
 *  * In forces_inline.hpp:
 *    - include the header file of the new interaction
 *    - in function @ref calc_bond_pair_force(): add the new interaction to the
 *      switch statement if it is a pairwise bond. For the FENE bond, the
 *      code looks like this:
 *      @code{.cpp}
 *      boost::optional<Utils::Vector3d> result;
 *      // ...
 *      case BONDED_IA_FENE:
 *        result = fene_pair_force(iaparams, dx);
 *        break;
 *      @endcode
 *    - in functions called by @ref add_bonded_force(): add the new interaction
 *      to the switch statements if it is not a pairwise bond. For the harmonic
 *      angle, the code looks like this:
 *      @code{.cpp}
 *      case BONDED_IA_ANGLE_HARMONIC:
 *        return angle_harmonic_force(p1.r.p, p2.r.p, p3.r.p, iaparams);
 *      @endcode
 *  * In energy_inline.hpp:
 *    - include the header file of the new interaction
 *    - in function @ref calc_bonded_energy(): add the new interaction to the
 *      switch statement. For the FENE bond, e.g., the code looks like this:
 *      @code{.cpp}
 *      case BONDED_IA_FENE:
 *        return fene_pair_energy(iaparams, dx);
 *      @endcode
 *  * Pressure tensor and virial calculation: if your bonded
 *    interaction is not a pair bond or modifies the particles involved,
 *    you have to implement a custom solution for virial calculation.
 *    The pressure occurs twice, once for the parallelized isotropic pressure
 *    and once for the tensorial pressure calculation. For pair forces, the
 *    pressure is calculated using the virials, for many body interactions
 *    currently no pressure is calculated.
 *
 *  @subsection bondedIA_new_interface Adding the interaction in the Python interface
 *
 *  Please note that the following is Cython code (www.cython.org), rather than
 *  pure Python.
 *
 *  * In file <tt>src/python/espressomd/interactions.pxd</tt>:
 *    - Import the parameter data structure from the C++ header file for the new
 *      interaction. For the FENE bond, this looks like:
 *      @code{.py}
 *      cdef extern from "bonded_interactions/bonded_interaction_data.hpp":
 *          cdef struct Fene_bond_parameters:
 *              double k
 *              double drmax
 *              double r0
 *              double drmax2
 *              double drmax2i
 *      @endcode
 *    - Add the bonded interaction to the Cython copy of the
 *      @ref BondedInteraction enum analogous to the one in the core
 *      The spelling has to match the one in the C++ enum exactly:
 *      @code{.py}
 *      cdef extern from "bonded_interactions/bonded_interaction_data.hpp":
 *          cdef enum enum_bonded_interaction "BondedInteraction":
 *              BONDED_IA_NONE = -1,
 *              BONDED_IA_FENE,
 *              [...]
 *      @endcode
 *    - Add the bonded interaction to the Cython copy of the
 *      @ref Bond_parameters union analogous to the C++ core.
 *      The member name has to match the one in C++ exactly:
 *      @code{.py}
 *      cdef union Bond_parameters:
 *          Fene_bond_parameters fene
 *          [...]
 *      @endcode
 *    - Import the declaration of the setter function implemented in the core.
 *      For the FENE bond, this looks like:
 *      @code{.py}
 *      cdef extern from "bonded_interactions/fene.hpp":
 *          int fene_set_params(int bond_type, double k, double drmax, double r0)
 *      @endcode
 *  * In file <tt>src/python/espressomd/interactions.pyx</tt>:
 *    - Implement the Cython class for the bonded interaction, using the one
 *      for the FENE bond as template. Please use pep8 naming convention.
 *  * In file <tt>testsuite/python/interactions_bonded_interface.py</tt>:
 *    - Add a test case to verify that parameters set and gotten from the
 *      interaction are consistent.
 *  * In file <tt>testsuite/python/interactions_bonded.py</tt> or
 *    <tt>testsuite/python/interactions_bond_angle.py</tt> or
 *    <tt>testsuite/python/interactions_dihedral.py</tt>:
 *    - Add a test case to verify the forces and energies are correct.
 *
 *  @section bondedIA_bond_angles Bond angle potentials
 *
 *  @subsection bondedIA_angle_force General expressions for the forces
 *
 *  This section uses the particle force expressions derived in @cite swope92a.
 *
 *  The gradient of the potential at particle @f$ i @f$ is given by the chain
 *  rule in equation 6:
 *
 *  @f{equation}{
 *      \label{eq:Swope-eq-6}
 *      \nabla_i U(\theta_{ijk})
 *          = \left(
 *                  \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *            \right)
 *            \left(
 *                  \frac{\mathrm{d}\theta_{ijk}}{\mathrm{d}\cos(\theta_{ijk})}
 *            \right)
 *            \left(
 *                  \nabla_i \cos(\theta_{ijk})
 *            \right)
 *  @f}
 *
 *  with
 *
 *  @f[
 *      \left(
 *            \frac{\mathrm{d}\theta_{ijk}}{\mathrm{d}\cos(\theta_{ijk})}
 *      \right)
 *          = \left(
 *                  \frac{-1}{\sin(\theta_{ijk})}
 *            \right)
 *  @f]
 *
 *  and @f$ \theta_{ijk} @f$ the angle formed by the three particles,
 *  @f$ U(\theta_{ijk}) @f$ the bond angle potential, @f$ \vec{r_{ij}} @f$
 *  the vector from particle @f$ j @f$ to particle @f$ i @f$ and
 *  @f$ \nabla_i @f$ the gradient in the direction @f$ \vec{r_{ij}} @f$.
 *
 *  The expression for @f$ \cos(\theta_{ijk}) @f$ is given by equation 4:
 *
 *  @f{equation}{
 *      \label{eq:Swope-eq-4}
 *      \cos(\theta_{ijk})
 *          = \frac{\vec{r_{ij}}\cdot\vec{r_{kj}}}
 *                 {\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|}
 *  @f}
 *
 *  The expression for its gradient is given by equation 9:
 *
 *  @f{equation}{
 *      \label{eq:Swope-eq-9}
 *      \nabla_i \cos(\theta_{ijk})
 *          = \vec{e_x}\frac{\partial\cos(\theta_{ijk})}{\partial x_{ij}}
 *          + \vec{e_y}\frac{\partial\cos(\theta_{ijk})}{\partial y_{ij}}
 *          + \vec{e_z}\frac{\partial\cos(\theta_{ijk})}{\partial z_{ij}}
 *  @f}
 *
 *  with @f$ \left(\vec{e_x}, \vec{e_y}, \vec{e_z}\right) @f$ the unit vectors
 *  of the reference coordinate system and
 *  @f$ \vec{r_{ij}} = \left(x_{ij}, y_{ij}, z_{ij}\right) @f$.
 *
 *  Applying the quotient rule:
 *
 *  @f[
 *      \frac{\partial\cos(\theta_{ijk})}{\partial x_{ij}}
 *          = \frac{\partial}{\partial x_{ij}}
 *            \left(
 *                  \frac{\vec{r_{ij}}\cdot\vec{r_{kj}}}
 *                       {\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|}
 *            \right)
 *          = \frac{\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|
 *                  \partial \left(\vec{r_{ij}}\cdot\vec{r_{kj}}\right)
 *                /\, \partial x_{ij}
 *                - \vec{r_{ij}}\cdot\vec{r_{kj}}\cdot
 *                  \partial
 *                  \left(
 *                       \left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|
 *                  \right)
 *                /\, \partial x_{ij}}
 *                 {\left\|\vec{r_{ij}}\right\|^2\left\|\vec{r_{kj}}\right\|^2}
 *  @f]
 *
 *  with
 *
 *  @f[
 *      \frac{\partial \left(\vec{r_{ij}}\cdot\vec{r_{kj}}\right)}
 *           {\partial x_{ij}}
 *          = \frac{\partial \left(x_{ij} \cdot x_{kj} + y_{ij} \cdot y_{kj} + z_{ij} \cdot z_{kj}\right)}
 *           {\partial x_{ij}}
 *          = x_{kj}
 *  @f]
 *
 *  and
 *
 *  @f[
 *      \frac{\partial \left(\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|\right)}
 *           {\partial x_{ij}}
 *          = \left\|\vec{r_{kj}}\right\|
 *            \frac{\partial}{\partial x_{ij}}
 *            \sqrt{x_{ij}^2 + y_{ij}^2 + z_{ij}^2}
 *          = \left\|\vec{r_{kj}}\right\|
 *            \frac{0.5 \cdot 2 \cdot x_{ij}}
 *                 {\sqrt{x_{ij}^2 + y_{ij}^2 + z_{ij}^2}}
 *          = x_{ij}
 *            \frac{\left\|\vec{r_{kj}}\right\|}
 *                 {\left\|\vec{r_{ij}}\right\|}
 *  @f]
 *
 *  leading to equation 12:
 *
 *  @f{align*}{
 *      \label{eq:Swope-eq-12}
 *      \frac{\partial\cos(\theta_{ijk})}{\partial x_{ij}}
 *         &= \frac{\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|x_{kj}
 *                  - \vec{r_{ij}}\cdot\vec{r_{kj}}\cdot x_{ij}
 *                    \left\|\vec{r_{kj}}\right\|\left\|\vec{r_{ij}}\right\|^{-1}}
 *                 {\left\|\vec{r_{ij}}\right\|^2\left\|\vec{r_{kj}}\right\|^2}
 *      \\
 *         &= \frac{x_{kj}}
 *                 {\left\|\vec{r_{ij}}\right\|\left\|\vec{r_{kj}}\right\|}
 *          - \frac{\vec{r_{ij}}\cdot\vec{r_{kj}}\cdot x_{ij}}
 *                 {\left\|\vec{r_{ij}}\right\|^3\left\|\vec{r_{kj}}\right\|}
 *      \\
 *         &= \frac{1}{\left\|\vec{r_{ij}}\right\|}
 *                 \left(
 *                    \frac{x_{kj}}{\left\|\vec{r_{kj}}\right\|}
 *                  - \frac{x_{ij}}{\left\|\vec{r_{ij}}\right\|}
 *                    \cos(\theta_{ijk})
 *                 \right)
 *  @f}
 *
 *  Applying these steps to equations 9-11 leads to the force equations
 *  for all three particles:
 *
 *  @f{align*}{
 *      \vec{F_i}
 *          &= - K(\theta_{ijk})
 *                 \frac{1}{\left\|\vec{r_{ij}}\right\|}
 *                 \left(
 *                    \frac{\vec{r_{kj}}}{\left\|\vec{r_{kj}}\right\|}
 *                  - \frac{\vec{r_{ij}}}{\left\|\vec{r_{ij}}\right\|}
 *                    \cos(\theta_{ijk})
 *                 \right)
 *      \\
 *      \vec{F_k}
 *          &= - K(\theta_{ijk})
 *                 \frac{1}{\left\lVert\vec{r_{kj}}\right\rVert}
 *                 \left(
 *                    \frac{\vec{r_{ij}}}{\left\|\vec{r_{ij}}\right\|}
 *                  - \frac{\vec{r_{kj}}}{\left\|\vec{r_{kj}}\right\|}
 *                    \cos(\theta_{ijk})
 *                 \right)
 *      \\
 *      \vec{F_j} &= -\left(\vec{F_i} + \vec{F_k}\right)
 *  @f}
 *
 *  with @f$ K(\theta_{ijk}) @f$ the angle force term, which depends on the
 *  expression used for the angle potential. Forces @f$ \vec{F_i} @f$ and
 *  @f$ \vec{F_k} @f$ are perpendicular to the displacement vectors
 *  @f$ \vec{r_{ij}} @f$ resp. @f$ \vec{r_{kj}} @f$ and their magnitude
 *  are proportional to the potential gradient normalized by the displacement
 *  vectors:
 *
 *  @f{align*}{
 *      \left\|\vec{F_i}\right\|
 *          &= \left(
 *                  \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *             \right)
 *             \frac{1}{\left\|\vec{r_{ij}}\right\|}
 *      \\
 *      \left\|\vec{F_k}\right\|
 *          &= \left(
 *                  \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *             \right)
 *             \frac{1}{\left\|\vec{r_{kj}}\right\|}
 *  @f}
 *
 *
 *  @subsection bondedIA_angle_potentials Available potentials
 *
 *
 *  @subsubsection bondedIA_angle_harmonic Harmonic angle potential
 *
 *  The harmonic angle potential takes the form:
 *
 *  @f{equation}{
 *      \label{eq:harmonic-angle-pot}
 *      U(\theta_{ijk})
 *          = \frac{1}{2}k_{ijk}\left[\theta_{ijk} - \theta_{ijk}^0\right]^2
 *  @f}
 *
 *  with @f$ \theta_{ijk} @f$ the angle formed by the three particles,
 *  @f$ \theta_{ijk}^0 @f$ the equilibrium angle and @f$ k_{ijk} @f$
 *  the bond angle force constant.
 *
 *  The derivative with respect to the angle is:
 *
 *  @f[
 *      \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *          = k_{ijk}\left[\theta_{ijk} - \theta_{ijk}^0\right]
 *  @f]
 *
 *  resulting in the following angle force term:
 *
 *  @f{equation}{
 *      \label{eq:harmonic-angle-pot-angle-term}
 *      K(\theta_{ijk})
 *          = -k_{ijk}\frac{\theta_{ijk} - \theta_{ijk}^0}
 *                         {\sin(\theta_{ijk})}
 *  @f}
 *
 *  which can lead to numerical instability at @f$ \theta_{ijk} = 0 @f$ and
 *  @f$ \theta_{ijk} = \pi @f$.
 *
 *
 *  @subsubsection bondedIA_angle_cossquare Harmonic cosine potential
 *
 *  The harmonic cosine potential takes the form:
 *
 *  @f{equation}{
 *      \label{eq:harmonic-cosine-pot}
 *      U(\theta_{ijk})
 *          = \frac{1}{2}
 *            k_{ijk}\left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right]^2
 *  @f}
 *
 *  with @f$ \theta_{ijk} @f$ the angle formed by the three particles,
 *  @f$ \theta_{ijk}^0 @f$ the equilibrium angle and @f$ k_{ijk} @f$
 *  the bond angle force constant.
 *
 *  The derivative with respect to the angle is:
 *
 *  @f[
 *      \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *          = -k_{ijk}\sin(\theta_{ijk})
 *            \left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right]
 *  @f]
 *
 *  resulting in the following angle force term:
 *
 *  @f{equation}{
 *      \label{eq:harmonic-cosine-pot-angle-term}
 *      K(\theta_{ijk})
 *          = k_{ijk}\left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right]
 *  @f}
 *
 *  which does not suffer from numerical instability.
 *
 *
 *  @subsubsection bondedIA_angle_cosine Cosine potential
 *
 *  The cosine potential takes the form:
 *
 *  @f{equation}{
 *      \label{eq:cosine-pot}
 *      U(\theta_{ijk})
 *          = k_{ijk}\left[1 - \cos(\theta_{ijk} - \theta_{ijk}^0)\right]
 *  @f}
 *
 *  with @f$ \theta_{ijk} @f$ the angle formed by the three particles,
 *  @f$ \theta_{ijk}^0 @f$ the equilibrium angle and @f$ k_{ijk} @f$
 *  the bond angle force constant.
 *
 *  The derivative with respect to the angle is:
 *
 *  @f[
 *      \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *          = k_{ijk}\sin(\theta_{ijk} - \theta_{ijk}^0)
 *  @f]
 *
 *  resulting in the following angle force term:
 *
 *  @f{equation}{
 *      \label{eq:cosine-pot-angle-term}
 *      K(\theta_{ijk})
 *          = -k_{ijk}\frac{\sin(\theta_{ijk} - \theta_{ijk}^0)}
 *                         {\sin(\theta_{ijk})}
 *  @f}
 *
 *  which can lead to numerical instability at @f$ \theta_{ijk} = 0 @f$ and
 *  @f$ \theta_{ijk} = \pi @f$.
 *
 *
 *  @subsubsection bondedIA_angle_tab Tabulated potential
 *
 *  The tabulated potential and its derivative with respect to the angle are
 *  provided by the user. The angle force term takes the form:
 *
 *  @f{equation}{
 *      \label{eq:tabulated-pot-angle-term}
 *      K(\theta_{ijk})
 *          = \frac{-1}{\sin(\theta_{ijk})}
 *            \left(
 *                  \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}}
 *            \right)
 *  @f}
 *
 *  which can lead to numerical instability at @f$ \theta_{ijk} = 0 @f$ and
 *  @f$ \theta_{ijk} = \pi @f$.
 *
 */