espressomd / espresso

The ESPResSo package
https://espressomd.org
GNU General Public License v3.0
221 stars 183 forks source link

Particle image box is wrong with Verlet lists #4940

Open jngrad opened 1 week ago

jngrad commented 1 week ago

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:

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.

RudolfWeeber commented 1 week ago

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.

jngrad commented 1 week ago

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.