Box length change invalidates particle positions

[jngrad]

When manually modifying the box length, the new particle positions are recalculated with the wrong box size.

MWE:

import espressomd
import numpy as np
np.set_printoptions(precision=2)
system = espressomd.System(box_l=[15., 15., 15.])
system.time_step = 0.01
system.cell_system.skin = 0.4
p = system.part.add(pos=[-0.1, 15.1, 30.])
expected = p.pos_folded * 10./15.
print(f"before:   {p.pos_folded}")
system.box_l = [10., 10., 10.]
print(f"after:    {p.pos_folded}")
print(f"expected: {expected}")

Output:

before:   [14.90 0.10 0. ]
after:    [ 4.90 0.10 0. ]
expected: [ 9.93 0.07 0. ]

This bug is reproducible on the python branch and with ESPResSo releases 4.2.1 and 4.1.4.

It is unclear to me, how to gracefully handle a box resize. One could fold the particle position onto the unit cell, apply the box resize via rescaling, and re-use the older image box counter. But then, a jump may become visible in trajectories in unfolded coordinates. One would also need to maintain the old and new value in the BoxGeometry object to apply the folding and then unfolding with the two box sizes. It is also unclear what should happen with Lees-Edwards boundary conditions, where particle positions gain an offset when going through one of the three periodic boundaries.

The easiest course of action is to prevent manually changing the box length if the box contains particles, and rely instead on the system.change_volume_and_rescale_particles(), which implements the correct behavior. Alternatively, one could make the box length setter call the rescale method under-the-hood, but this introduces a new problem: the on_boxl_change event would be called three times, causing long-range solvers to re-tune all solver parameters three times with a different box size. Since the box might no longer be cubic for the first two re-tunes, the tuning code may become extremely slow, have to sample more mesh sizes, and throw an exception if the accuracy cannot be reached.

Note: this bug doesn't affect the NpT integrator, which uses the rescale algorithm.