Bonded interactions

Table of Contents

• Adding new interactions
  • Defining the new interaction
  • Integrating the new bond type into the C++ core
  • Registering the new interaction in the ScriptInterface
  • Adding the interaction in the Python interface
• Bond angle potentials
  • General expressions for the forces
  • Available potentials
    • Harmonic angle potential
    • Harmonic cosine potential
    • Cosine potential
    • Tabulated potential

Adding new interactions

To add a new bonded interaction, the following steps are required:

• C++ core:
  • Create a new header file for your new bond type, preferably inside the folder src/core/bonded_interactions.
  • Define a data structure for the interaction parameters (prefactors, etc.).
  • Write a constructor and force/energy kernels.
  • Add calls to these kernels in the force, energy and pressure calculation functions.
  • Register the new bond type.
• ScriptInterface:
  • Define the ScriptInterface class of the new bond type, which serves as the connection between the C++ core and the Python representation of the bond.
• 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

Defining the new interaction

Every interaction resides in its own source .hpp file. A simple example for a bonded interaction is the FENE bond in fene.hpp. Use this file as template for your new interaction.

The first step is to create a new struct which represents your new bond type inside the .hpp file. It needs to have the following members:

• static constexpr int num = 1;

  This is the number of particles involved in the bond, minus 1, i.e. 1 for a pairwise bonded potential such as the FENE bond.

• double cutoff() const { return r0 + drmax; }

  The return value of cutoff() should be as large as the interaction range of the new interaction. This is only relevant to pairwise bonds. In all other cases, the return value should be 0, namely angle bonds, dihedral bonds as well as other bonds that don't have an interaction range. The VirtualBond is the exception to this rule; its range is BONDED_INACTIVE_CUTOFF to ensure that it is always skipped by the short-range loop.

• std::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;

  This function returns the bond force. If it is a bond involving three or four particles, a std::tuple with three or four force vectors has to be returned, respectively.
  • The returned value is in a std::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 function can make use of a pre-calculated distance vector (dx) pointing from particle 2 to particle 1, that takes periodic boundary conditions into account.
  • The function name and signature may be different, e.g. for angle bonds. To determine the right signature for the new interaction, you can have a look at forces_inline.hpp to see where this function will be called and which other variables may be available for your calculation.

• std::optional<double> energy(Utils::Vector3d const &dx) const;

  This function returns the bond energy. The same information as given for the force calculation above applies here. This function will be called from energy_inline.hpp .
• A constructor which is used to set all parameters and to do preparatory calculations if necessary, for example

  FeneBond(double k, double drmax, double r0);

  All values the bond needs to function properly should be passed as arguments to this constructor.

Integrating the new bond type into the C++ core

In most cases, there are three files that need to be updated to integrate the new bond type into the core, namely

• bonded_interaction_data.hpp
• forces_inline.hpp
• energy_inline.hpp

In some cases, you may also need to modify pressure_inline.hpp.

• In bonded_interaction_data.cpp:
  • Include the header file containing the new bond type.
  • Add the new bond type to Bonded_IA_Parameters at the end of the types list.
  • If by doing this, the length of the list in Bonded_IA_Parameters passes over a multiple of 10, you may have to update ESPResSo's top level CMakeLists.txt:

    # enable boost::variant with more than 20 types
    target_compile_options(
      espresso_cpp_flags INTERFACE -DBOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
                                   -DBOOST_MPL_LIMIT_LIST_SIZE=40)

• In forces_inline.hpp:
  • A call to the new bond's force calculation needs to be placed in either of the functions calc_bond_pair_force(), calc_bonded_three_body_force() or calc_bonded_four_body_force(), depending on how many bond partners there are.
  • Add the new entry to the if - else chain, like in the following example

    // ...
    else if (auto const *iap = boost::get<QuarticBond>(&iaparams)) {
      return iap->force(dx);
    }
    // ...

• In energy_inline.hpp:
  • A call to the new bond's force calculation needs to be placed in calc_bonded_energy. Find the if - else chain that corresponds to the correct number of bond partners.
  • Add the new entry to the if - else chain, like in the following example

    // ...
    else if (auto const *iap = boost::get<QuarticBond>(&iaparams)) {
      return iap->energy(dx);
    }
    // ...

• Pressure tensor and virial calculation (pressure_inline.hpp): 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.

Note that the force and energy functions cannot alter the particle states.

Registering the new interaction in the ScriptInterface

• In src/script_interface/interactions/BondedInteraction.hpp: Add a new class representing your new bond type in the ScriptInterface.
  • We recommend that the new class has the same name as the interaction in the core.
  • You can use ScriptInterface::Interactions::FeneBond as a template.
  • The class must be derived from ScriptInterface::Interactions::BondedInteraction.
  • It is recommended to include the statement

    using CoreBondedInteraction = ::YourNewBond;

    where YourNewBond is the core type you defined.
  • Implement a member function with the signature

    void construct_bond(VariantMap const &params) override { /* ... */ }

    In this function, the member m_bonded_ia shall be initialized using the parameters that are given in params. Use the construct_bond() method of FeneBond as a template. An instance of the core type YourNewBond should be initialized, which is then used to initialize a std::shared_ptr to a Bonded_IA_Parameters, which is then assigned to m_bonded_ia. The values of the parameters are extracted from params using

    get_value<parameter_type>(params, "parameter_name")

    where parameter_type is the type of the parameter, e.g. double or int or even std::string, and "parameter_name" must be replaced by the name of the respective parameter. This name must be the same as in the Python interface, but may differ from the name in the core interaction type. It is, however, recommended to use the same names for both the Python interface and the ESPResSo core for consistency whenever possible.
  • Implement the constructor. We recommend to adapt it from FeneBond. All it needs to do is to register its parameters so they can be set from Python. For this purpose, call

    add_parameters(/* ... */);

    inside the constructor. It expects a vector of ScriptInterface::AutoParameter. Usually, this vector is initialized using an initializer list, each element of which is in itself a list which initializes one instance of AutoParameter (see AutoParameter.hpp). One of many ways to initialize these parameters is to pass the parameter name as a string, a custom setter and a custom getter function. The parameters are typically made read-only by passing AutoParameter::read_only instead of a setter function. The getter function can be a lambda function, which is directly passed to the constructor.
• In src/script_interface/interactions/initialize.cpp: Your new interaction type needs to be registered here so that the ScriptInterface can find it by its name. In the function initialize add a new line of the form

  om->register_new<YourNewBond>("Interactions::YourNewBond");

  where YourNewBond must be replaced by the name of your new bond type. The string is used to match the ScriptInterface class with the Python class (see below).

Adding the interaction in the Python interface

• In file src/python/espressomd/interactions.py:
  • Add the bonded interaction to BONDED_IA. The order of the enum values must match the order of types in Bonded_IA_Parameters exactly:

    class BONDED_IA(enum.IntEnum):
        NONE = 0
        FENE = enum.auto()
        HARMONIC = enum.auto()
        # ...

  • Implement the Python class for the bonded interaction, using the one for the FENE bond as template. Please use pep8 naming convention.
  • Set the class' members

    _so_name = "Interactions::YourNewBond"
    _type_number = BONDED_IA.YOURNEWBOND

    where you put the name of your bond instead of YourNewBond. This connects the ScriptInterface class with the Python class. The type number is the new enum value from BONDED_IA.
• In file testsuite/python/interactions_bonded_interface.py:
  • Add a test case to verify that parameters set and gotten from the interaction are consistent.
• In file testsuite/python/interactions_bonded.py or testsuite/python/interactions_bond_angle.py or testsuite/python/interactions_dihedral.py:
  • Add a test case to verify the forces and energies are correct.

Bond angle potentials

General expressions for the forces

This section uses the particle force expressions derived in [40].

The gradient of the potential at particle \( i \) is given by the chain rule in equation 6:

\begin{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) \end{equation}

with

\[ \left( \frac{\mathrm{d}\theta_{ijk}}{\mathrm{d}\cos(\theta_{ijk})} \right) = \left( \frac{-1}{\sin(\theta_{ijk})} \right) \]

and \( \theta_{ijk} \) the angle formed by the three particles, \( U(\theta_{ijk}) \) the bond angle potential, \( \vec{r_{ij}} \) the vector from particle \( j \) to particle \( i \) and \( \nabla_i \) the gradient in the direction \( \vec{r_{ij}} \).

The expression for \( \cos(\theta_{ijk}) \) is given by equation 4:

\begin{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\|} \end{equation}

The expression for its gradient is given by equation 9:

\begin{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}} \end{equation}

with \( \left(\vec{e_x}, \vec{e_y}, \vec{e_z}\right) \) the unit vectors of the reference coordinate system and \( \vec{r_{ij}} = \left(x_{ij}, y_{ij}, z_{ij}\right) \).

Applying the quotient rule:

\[ \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} \]

with

\[ \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} \]

and

\[ \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\|} \]

leading to equation 12:

\begin{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) \end{align*}

Applying these steps to equations 9-11 leads to the force equations for all three particles:

\begin{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) \end{align*}

with \( K(\theta_{ijk}) \) the angle force term, which depends on the expression used for the angle potential. Forces \( \vec{F_i} \) and \( \vec{F_k} \) are perpendicular to the displacement vectors \( \vec{r_{ij}} \) resp. \( \vec{r_{kj}} \) and their magnitude are proportional to the potential gradient normalized by the displacement vectors:

\begin{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\|} \end{align*}

Available potentials

Harmonic angle potential

The harmonic angle potential takes the form:

\begin{equation} \label{eq:harmonic-angle-pot} U(\theta_{ijk}) = \frac{1}{2}k_{ijk}\left[\theta_{ijk} - \theta_{ijk}^0\right]^2 \end{equation}

with \( \theta_{ijk} \) the angle formed by the three particles, \( \theta_{ijk}^0 \) the equilibrium angle and \( k_{ijk} \) the bond angle force constant.

The derivative with respect to the angle is:

\[ \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}} = k_{ijk}\left[\theta_{ijk} - \theta_{ijk}^0\right] \]

resulting in the following angle force term:

\begin{equation} \label{eq:harmonic-angle-pot-angle-term} K(\theta_{ijk}) = -k_{ijk}\frac{\theta_{ijk} - \theta_{ijk}^0} {\sin(\theta_{ijk})} \end{equation}

which can lead to numerical instability at \( \theta_{ijk} = 0 \) and \( \theta_{ijk} = \pi \).

Harmonic cosine potential

The harmonic cosine potential takes the form:

\begin{equation} \label{eq:harmonic-cosine-pot} U(\theta_{ijk}) = \frac{1}{2} k_{ijk}\left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right]^2 \end{equation}

with \( \theta_{ijk} \) the angle formed by the three particles, \( \theta_{ijk}^0 \) the equilibrium angle and \( k_{ijk} \) the bond angle force constant.

The derivative with respect to the angle is:

\[ \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}} = -k_{ijk}\sin(\theta_{ijk}) \left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right] \]

resulting in the following angle force term:

\begin{equation} \label{eq:harmonic-cosine-pot-angle-term} K(\theta_{ijk}) = k_{ijk}\left[\cos(\theta_{ijk}) - \cos(\theta_{ijk}^0)\right] \end{equation}

which does not suffer from numerical instability.

Cosine potential

The cosine potential takes the form:

\begin{equation} \label{eq:cosine-pot} U(\theta_{ijk}) = k_{ijk}\left[1 - \cos(\theta_{ijk} - \theta_{ijk}^0)\right] \end{equation}

with \( \theta_{ijk} \) the angle formed by the three particles, \( \theta_{ijk}^0 \) the equilibrium angle and \( k_{ijk} \) the bond angle force constant.

The derivative with respect to the angle is:

\[ \frac{\mathrm{d}U(\theta_{ijk})}{\mathrm{d}\theta_{ijk}} = k_{ijk}\sin(\theta_{ijk} - \theta_{ijk}^0) \]

resulting in the following angle force term:

\begin{equation} \label{eq:cosine-pot-angle-term} K(\theta_{ijk}) = -k_{ijk}\frac{\sin(\theta_{ijk} - \theta_{ijk}^0)} {\sin(\theta_{ijk})} \end{equation}

which can lead to numerical instability at \( \theta_{ijk} = 0 \) and \( \theta_{ijk} = \pi \).

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:

\begin{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) \end{equation}

which can lead to numerical instability at \( \theta_{ijk} = 0 \) and \( \theta_{ijk} = \pi \).