Bonds get lost when deleting an uninvolved particle

[jonaslandsgesell]

Starting from commit a2bb7e1 from 31st October 2019, there is unexpected behaviour: deleting a particle without bond results in removing some bonds from particles which have bonds.
The reaction ensemble inserts and deletes particles, with the deletion taking place in function ReactionAlgorithm::delete_particle(int p_id) in reaction_ensemble.cpp.
Belwo you find a minimal working example of the problem - the script throws if some bonds got unintendetly deleted.

Using the n squared cell system removes the wrong behaviour.
First reported by Oleg Rud @helvrud :

########################################################
# The srip illustrationg the bonds disappearence issue #
########################################################
from __future__ import print_function
import espressomd
from espressomd import interactions, electrostatics, diamond, reaction_ensemble
import numpy as np
import pprint
pprint = pprint.PrettyPrinter(indent = 4).pprint

type_na = 0
type_cl = 1
type_chains = 3

MPC = 10 # monomers per chain

# length for Kremer-Grest chain
bond_length = 0.966

# The physical distance between nodes such that a line of monomers "fit" needs to be worked out.
# This is done via the unit diamond lattice size parameter "a".
a = (MPC + 1) * bond_length / (0.25 * np.sqrt(3))

# Lastly, the created periodic connections requires a specific simulation box.

system = espressomd.System(box_l = [a]*3)
system.time_step = 0.005
#system.cell_system.set_n_square()

# set the bonded interaction between the gel beads
bond = interactions.FeneBond(k=10.0, d_r_max=2., r_0 = 0.)
system.bonded_inter.add(bond)

from espressomd import diamond
# The physical distance between nodes such that a line of monomers "fit" needs to be worked out.
# This is done via the unit diamond lattice size parameter "a".
diamond.Diamond(a=a, bond_length=bond_length, MPC=MPC)

system.part[:].type=type_chains

# Set up salt ion insertions
RE = reaction_ensemble.ReactionEnsemble(temperature=1., seed = 1, exclusion_radius=1.)
gamma = np.exp( -3 )    
RE.add_reaction(gamma=gamma,
    reactant_types=[],
    reactant_coefficients=[],
    product_types=[0,1],
    product_coefficients=[1, 1],
    default_charges={0:1,1:-1})



num_gel_particles=16*MPC+8
bonds_from_a=[]
bonds_to_a=[]
for i in range(num_gel_particles):
    for j in range(len(system.part[i].bonds)):
        bonds_from_a.append(i)
        bonds_to_a.append(system.part[i].bonds[j][1])

import copy
a=copy.copy(system.part[0:num_gel_particles-1].bonds)

sim_steps = 100
print("simulating...")
#system.integrator.run(sim_steps)
for i in range(10000):
    RE.reaction(1)

    b=copy.copy(system.part[0:num_gel_particles-1].bonds)

    bonds_from_b=[]
    bonds_to_b=[]
    for i in range(num_gel_particles):
        for j in range(len(system.part[i].bonds)):
            bonds_from_b.append(i)
            bonds_to_b.append(system.part[i].bonds[j][1])
               
    print("hashes", hash(str(a)), hash(str(b)) )
    print("nr bonds", len(bonds_from_a), len(bonds_from_b))
    if (len(system.part.select(type=type_chains))==num_gel_particles):
        print("all gel particles still exist")            
    if (len(bonds_from_a) != len(bonds_from_b)):
        raise RuntimeError("Bond lost")