Store bonds on all partners

[RudolfWeeber]

Current situation

Currently, a pair bond between particles p1 and p2 or a 3- or 4-particle bond between p1, p2, p3(, p4) is only stored on p1.

Advantages:

• Save some storage
• Calculate the bond energy/force/stress only once

Disadvantages:

• The bond calculation cannot be shared-memory parallelized without atomic writes on the particle forces (because working on the bonds of p1 involves adding forces to p2)
• In MPI-parallel simulations, a reduction on forces in ghost layers is needed (because a particle might receive a bond force from a bond stored on a particle managed by the neighboring processor)
• For pair bonds, the one-sided storage does not represent the physics (the Hamiltonian being symmetric with respect to the bond partners).

Proposed change

• The bonds will be stored on all involved particles
• The Python interface will reflect this change, particularly with respect to pair bonds which are symmetric with regards to exchanging the two particles

Steps:

• Fully hide the bond storage from the rest of ESPResSo. For read access, this is already done via the BondView class. But for bond creation and deletion, the explicit storage format is still visible in particle_data.cpp. Once this is done, the storage format can be changed
• Alternative A: Store bonds as [bond_type, p1.id, p2.id] on all involved particles. This storage format can handle pair- and 3/4-particle bonds.
• Alternative B: Store pair bonds separately as [bond_id, partner.id], and use the full format [bond_id, p1.id, p2.id, p3.id] on angle and dihedral bonds (because for those, the order matters. It has to be known on each of the involved particles, which one is the tip of the triangle)

Proposed Bonds Python interface

• Bonds will no longer be accessible via the ParticleHandle.
• Instead, they are accessed by system.bonds[bond], where bond is an instance of one of the descendants of BondedInteraction, such as HarmonicBond
• system.bonds[bond] contains tuples of particles which are connected by that bond. It obeys set-semantics
• The elements in system.bonds[bond] have
  • set-like behaviour for pair bonds, i.e., PairBondPartners(p1,p2) == PairBondPartners(P2,P1), reflecting the fact that pair bonds are symmetrical with respect to exchanging the two particles
  • tuple-like behavior for bonds involving 3 or more particles. I.e, MultiBondPartners(p1,p2,3) != MultiBondPartners(p3,p1,p2)

Example

harmonic = HarmonicBonx(...)
angle_harmonic = AngleHarmonic(...)

p1 = system.part.add(...)
p2 = system.part.add(...)
p3 = system.part.add(...)

system.bonds[harmonic].add((p1,p2))
system.bonds[harmonic].add((p2,p3))

assert (p1,p2) in system.bonds[harmonic] 
assert (p2,p1) in system.bonds[harmonic]

system.bonds[angle_harmonic].add((p2,p1,p3)) # p2 at the tip of the angle

[andreashaertel]

My only concern would be the size of additionally needed memory, if it were too large - what, I think, is not the case. Apart from this concern, the proposal sounds reasonable.