Error in Newtonian Event Chain algorithm

[edwinb-ai]

Hello,

I have been testing the Newtonian Event Chain (NEC) implementation in hoomd-blue, however, I get an error that says
**ERROR**: Two particles with the same distance.

After looking at the configuration and the radial distribution function in OVITO, I don't see any overlaps whatsoever. Here a screenshot from OVITO.

I really don't know if this is a bug, or if I am doing something wrong, so I thought I should ask first before opening an issue.

Here is the code, and the output of the code.

from itertools import product

import gsd.hoomd
import hoomd
import numpy as np

diameter = 1.0

m = 10
total_particles = 4 * m**3
n_particles = 2048
spacing = 1.4
K = np.ceil(np.sqrt(total_particles)).astype(np.int64)
L = K * spacing
x = np.linspace(-L / 2, L / 2, K, endpoint=False)
# Repeat three times once again for a 3D system
position = np.array(list(product(x, repeat=3)))
position = position[:n_particles]
# position = np.column_stack((position, np.zeros((n_particles,))))
print(position.shape)
orientation = [(1, 0, 0, 0)] * n_particles

frame = gsd.hoomd.Frame()
frame.particles.N = n_particles
frame.particles.position = position
frame.particles.orientation = orientation
frame.particles.typeid = np.zeros((n_particles,))
frame.particles.types = ["A"]
# Set the box dimensions to be 2D
frame.configuration.box = [L, L, L, 0, 0, 0]

with gsd.hoomd.open("initial.gsd", "w") as f:
    f.append(frame)

rng = np.random.default_rng()
device = hoomd.device.CPU()
sim = hoomd.Simulation(device=device, seed=rng.integers(2**16))
sim.create_state_from_snapshot(snapshot=frame)

mc = hoomd.hpmc.integrate.Sphere(default_d=0.4)
mc.shape["A"] = dict(diameter=diameter)
sim.operations.integrator = mc
# We use a `QuickCompress` updater to compress the box to the target
# packing fraction
V_particle = 4.0 * np.pi * (diameter / 2.0) ** 3 / 3.0
initial_box = sim.state.box
final_box = hoomd.Box.from_box(initial_box)
final_volume_fraction = 0.49
final_box.volume = sim.state.N_particles * V_particle / final_volume_fraction
compress = hoomd.hpmc.update.QuickCompress(
    trigger=hoomd.trigger.Periodic(10), target_box=final_box
)
sim.operations.updaters.append(compress)
# Adjust the displacement during then compression steps
compression_steps = 1e6
tune = hoomd.hpmc.tune.MoveSize.scale_solver(
    moves=["d"],
    target=0.4,
    trigger=hoomd.trigger.Periodic(10),
    max_translation_move=0.2,
)
sim.operations.tuners.append(tune)

while not compress.complete and sim.timestep < compression_steps:
    sim.run(1000)

device.notice("Compression done!")

translate_moves = mc.translate_moves
translation_acceptance = translate_moves[0] / sum(translate_moves)
device.notice(f"Translation acceptance: {translation_acceptance}")

sim.operations.updaters.clear()
sim.operations.tuners.clear()

sim.run(1e3)

translate_moves = mc.translate_moves
translation_acceptance = translate_moves[0] / sum(translate_moves)
device.notice(f"Translation acceptance: {translation_acceptance}")

hoomd.write.GSD.write(state=sim.state, filename="compressed.gsd")

# Assign velocities to the particles
sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT=1)

nec = hoomd.hpmc.nec.integrate.Sphere(
    default_d=0.01,
    update_fraction=0.05,
    chain_time=0.05,
    nselect=2,
)
nec.shape["A"] = dict(diameter=diameter)
sim.operations.integrator = nec

# Tuners
tuner_trigger = hoomd.trigger.And(
    [hoomd.trigger.Periodic(50), hoomd.trigger.Before(sim.timestep + 1000)]
)
tune_nec_d = hoomd.hpmc.tune.MoveSize.scale_solver(
    tuner_trigger, moves=["d"], target=0.1, max_translation_move=0.2
)
sim.operations.tuners.append(tune_nec_d)
chain_tuner = hoomd.hpmc.nec.tune.ChainTime.scale_solver(
    trigger=tuner_trigger,
    target=20,
    tol=1,
    gamma=20,
)
sim.operations.tuners.append(chain_tuner)

sim.run(1e4)

sim.run(1e4)

translate_moves = nec.translate_moves
print(
    "Translation collision hits:",
    translate_moves[0] / sum(translate_moves),
    " d=",
    nec.d["A"],
)
print("Particles per chain:", nec.particles_per_chain, " chain_time=", nec.chain_time)

sim.operations.tuners.clear()

sim.run(1e4)

translate_moves = nec.translate_moves
print(
    "Translation collision hits:",
    translate_moves[0] / sum(translate_moves),
    " d=",
    nec.d["A"],
)
print("Particles per chain:", nec.particles_per_chain, " chain_time=", nec.chain_time)

hoomd.write.GSD.write(state=sim.state, filename="nec_snap.gsd")

(2048, 3)
Compression done!
Translation acceptance: 0.5355965576171875
Translation acceptance: 0.4109993896484375
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two pTranslation collision hits: 0.10732273000825619  d= 0.12900415712911925
Particles per chain: 19.4593786407767  chain_time= 0.7204762479147064
Translation collision hits: 0.10730868714359812  d= 0.12900415712911925
Particles per chain: 19.464718932038835  chain_time= 0.7204762479147064
articles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance
**ERROR**: Two particles with the same distance

[joaander]

Thanks for the report and the fully working example. I am able to reproduce the error message you are seeing. I also checked the logs in the validation tests I run on each HOOMD-blue release and find the same messages there. At the same time, these strict validation tests pass within differences less than two standard deviations of the mean (images below).

I am not the original developer of the NEC code, so can only guess at the intent of this message. What is happening is that the particle being moved collides with two other particles at exactly the same time. I don't know how this is happening (the chances of a 3-way collision should be vanishingly small), but it is. In event-driven MD simulations (EDMD), my understanding is that a 3-way collision is an error as the velocities of all three particles cannot be correctly updated reflecting energy and momentum conservation. NEC borrows the idea of energy and momentum conserving collisions from EDMD, but applies them in a Monte Carlo context.

In that context, NEC will produce valid equilibrium distributions no matter what velocities are chosen after the collision, whether they are Newtonian, random, or even fixed (see the original event chain Monte Carlo algorithm as an example). Read the paper on NEC and you will see that the Newtonian collisions significantly decrease the time needed to sample phase space.

In other words, the error message is not actually an error. It is only noting that what should be an extremely rare event occurred. I will remove the message in a future release. When a 3-way collisions occurs, NEC solves the elastic collision between two of the particles and does not modify the velocity of the third.

Hard disks (hoomd-blue 5.1.1):

Hard spheres (hoomd-blue 5.1.1):

[edwinb-ai]

Thanks a lot Joshua, that is great to know! Thank you also for the detailed explanation. For now I will continue to use it, since I am satisfied it is not doing anything out of the ordinary.