Particle image box is wrong with Verlet lists
jngrad opened this issue · 2 comments
When Verlet lists are used, and the Verlet skin is too large, particles are no longer resorted and their image box is miscounted.
MWE:
import espressomd
import numpy as np
import matplotlib.pyplot as plt
system = espressomd.System(box_l=[10.0, 10.0, 10.0])
system.time_step = 0.01
system.cell_system.skin = 0.4 # bug with 0.4, ok with 0.04 and 0.004
system.periodicity = [True, True, True]
x0 = np.array([9.95, 0.0, 0.0])
v0 = np.array([0.1, 0.0, 0.0])
p = system.part.add(pos=x0, v=v0)
pos_unfolded_ref = []
pos_unfolded = []
pos_folded = []
image_box = []
for i in range(10):
pos_unfolded_ref.append(x0 + system.time * v0)
pos_unfolded.append(p.pos)
pos_folded.append(p.pos_folded)
image_box.append(p.image_box)
system.integrator.run(10)
pos_unfolded_ref = np.array(pos_unfolded_ref)
pos_unfolded = np.array(pos_unfolded)
pos_folded = np.array(pos_folded)
image_box = np.array(image_box)
plt.plot(pos_unfolded_ref[:,0], image_box[:,0], "o-")
plt.show()
Three skin values can be used in this MWE:
skin=0.4
: no particle resort happens besides the first time step (local resort)skin=0.04
: particle resort (local resort) happens in one cellskin=0.004
: particle resort (local resort) happens in multiple cells
Depending on the choice of initial conditions, e.g. with x0=[9.50, 0.0, 0.0]
and skin=0.4
, the image box is updated, but with a random time lag, i.e. it gets incremented when the position is anywhere between box_l
and box_l + skin
. The only particle folding test we have doesn't use Verlet lists.
This bug is reproducible on the Python branch, and the 4.2.2 and 4.1.4 releases. Other releases were not checked.
Many thanks to Xin Yuan for reporting this issue with a MWE on the mailing list thread Position folding issue.
I think this needs to be fixed in the getter in the script interface.
There, currently we just return the iage box. Instead, fold_positoin should be used to calculate the updated image box from position and image box on the particle. Cf. the pos_folded getter.
Our hdf5 interface is also affected by the bug.
Applying @RudolfWeeber's suggestion in both ScriptInterface::Particles::ParticleHandle::ParticleHandle()
and
Writer::H5md::File::write()
seems to fix the bug, with and without Lees-Edwards boundary conditions.