Bond_parameters should be a variant

[fweik]

This is currently a tagged union, which requires manual construction and destruction when changing the active type (the union contains non-trivial types). This requires boiler plate code for switching the active type and is error prone. This can be made automatic by turning the union into a boost::variant.

[fweik]

Is anybody working on this?

[jngrad]

I'm not aware of anyone working on this.

[jngrad]

I gave it a try yesterday, but the refactoring effort is quite significant because the kernels are written as free functions. Today I tried to break it down into smaller tasks, and came up with the following plan:

• make bond parameter structs serializable (easy)
  • bond communication can then be simplified substantially
• move the force/energy kernels into the corresponding structs (difficult)
  • dependency inversion required in bonded_interaction_data.hpp: the kernels are currently separated from the corresponding structs so that the short range loop has access to the cutoff method but not to the kernels (which bring in too many dependencies)
• make kernel signatures more homogeneous
  • pass positions instead of particles in IBM_*_CalcForce()
  • not sure what to do about bonded_coulomb_pair_{force,energy}() and calc_oif_local()
• create visitor patterns for energy and force kernels

To illustrate this strategy, a small C++ prototype is attached (variant.cpp). The key part is this line:

// calculate the bonded energy for a distance vector `vec` and bond `bond_id`
auto energy = boost::apply_visitor([&vec](auto const &ia){ return ia.energy(vec); },
                                   bonded_ia_params[bond_id]);

The difficulty here is that the number of arguments in the kernels is different for pair, angle and dihedral bonds, and further vary for OIF and coulomb kernels, yet the visitor pattern requires a homogeneous function signature. To work around that, we could define in every bond pair an energy(v1, v2) {return 0;} overload and in every bond angle an energy(v) {return 0;} overload (see variant_angle.cpp), but this is not a clean solution, and we definitively don't want the OIF and Coulomb kernel signatures to bleed into all the other bonds.

Notes for the compiler explorer: the macro INLINE_KERNEL is commented out to make the function calls stand out in the generated assembly. If you use Clang 11 instead of GCC 10.2, make sure to enable Boost 1.71+ to get less boilerplate instructions.

To write templated visitors that return functions, we need to know the function signature during template instantiation. We cannot write visitors that store the arguments in private members to later pass them to the kernels (reification), because the number of arguments is different for OIF and Coulomb. I couldn't find an example of visitor pattern that fulfills these constraints.

To solve this, we could "unwrap" the visitor to pass the arguments directly. This can be done by adapting the switch statements in {energy,forces}_inline.hpp, as shown below. If we go with this option, we can keep the kernels as free functions.

double calc_pair_bond_energy(Bonded_ia_parameters const *iaparams, Particle const &p1, Particle const &p2) {
  auto const dx = get_mi_vector(p1.r.p, p2.r.p, box_geo);
  if (auto* ia = boost::get<Harmonic_bond_parameters>(iaparams)) {
    return ia->energy(dx);
  } else if (auto* ia = boost::get<Coulomb_bond_parameters>(iaparams)) {
    return ia->energy(dx, p1.p.q * p2.p.q);
  }
  return 0;
}

[RudolfWeeber]

I think the plan makes sense. I'd say we go ahead with the pair bonds and the visitor pattern.

The handful of exceptions we can handle in a switch statement as in your 2nd example, initially.

Once the bonds are changed in the core, we should also replace the interface to use ScriptInterface.

There may be more reasons, to split pair bonds from the rest. They can, e.g., be stored symetrically on both particles. That would allow for shared-memory parallel evaluation without a particle index and would make the force ghost communication un-necessary for systems with only non-bonded interactions and pair bonds.

[jngrad]

offline discussion:

• implementing a visitor pattern for the subset of pair interactions that only take the positions as arguments isn't worth it
• using an if ... else if chain to replace the current switch statement should be sufficient for our needs, although it's important to remember to add a new else if for every new bonded interaction (src/core/bonded_interactions/bonded_interactions.dox needs to reflect that)
• the force/energy kernels can take the struct as parameter directly, instead of the current union
• see if it's possible to move the class declarations in the respective .hpp files, so that the kernels are defined inside the class
  • watch out for dependency inversion issues that might arise during this change
  • watch out for header files visibility, e.g. electrostatics and thermostats should not become indirectly visible to the entire codebase
  • currently the kernels are separated from the class because we need them to be inline for kernel consumers ({energy,forces}_inline.hpp), while not being visible to other components of the core (e.g. the short_range_loop() that consumes the bond cutoff()), because the kernels expose a lot of dependencies (electrostatics, ScaFaCoS, etc.)

[jngrad]

Depends on #4079.

[jngrad]

Prototype: jngrad/espresso@fix-3019-prototype. The last commit can be removed, i.e. start from jngrad/espresso@7e7b0ea4e898.