Stable algorithm for event detection in event-driven particle dynamics

Event-Driven Particle Dynamics is a fast and precise method to simulate particulate systems of all scales. In this work it is demonstrated that, despite the high accuracy of the method, the finite machine precision leads to simulations entering invalid states where the dynamics are undefined. A general event-detection algorithm is proposed which handles these situations in a stable and efficient manner. This requires a definition of the dynamics of invalid states and leads to improved algorithms for event-detection in hard-sphere systems.

💡 Research Summary

Event‑Driven Particle Dynamics (EDPD) is a powerful simulation technique that updates particle states only at discrete collision events, thereby eliminating the need for small time steps and achieving high computational efficiency. Despite its theoretical exactness, practical implementations inevitably suffer from finite‑precision floating‑point arithmetic. Small rounding errors in the calculation of collision times accumulate, especially in dense or high‑velocity systems, and can drive the simulation into “invalid states” where particles overlap or are separated by less than the hard‑sphere diameter. In such states the usual collision rules become undefined, leading to incorrect event ordering, infinite loops, or crashes.

The authors first formalize the origin of these invalid states. For a pair of particles i and j with relative position rij and relative velocity vij, the collision condition |rij + vij t| = σ (σ = particle diameter) yields a quadratic equation in t. Solving this equation with double‑precision arithmetic produces an approximate root tc. If the deviation Δt = |tc – texact| exceeds a machine‑dependent tolerance ε, the algorithm may schedule a collision that either occurs after the particles have already overlapped or fails to detect an imminent overlap. The paper quantifies how Δt scales with particle speed, separation, and the condition number of the quadratic, demonstrating that even a modest ε can cause catastrophic failures after a few hundred thousand events in typical simulations.

To address the problem, the paper introduces a rigorous definition of “invalid‑state dynamics.” When an overlap is detected, the simulation immediately applies a corrective step that either (1) projects the overlapping particles back to the exact contact distance while conserving linear momentum and kinetic energy, or (2) inserts a virtual collision at the earliest feasible time that restores a physically admissible configuration. This corrective mechanism is split into two phases. The prediction phase augments the standard quadratic solver with a dynamic error check: if the residual of the quadratic exceeds ε, a secondary root‑finding routine (e.g., bisection or Newton‑Raphson with interval arithmetic) refines the collision time to within the prescribed tolerance. The verification phase re‑evaluates particle positions after the event; if any residual overlap remains, a tiny auxiliary time step δt is introduced, and an additional collision is forced to eliminate the overlap before proceeding.

Implementation details focus on managing floating‑point uncertainty through interval arithmetic and adaptive error bounds. For each event the algorithm computes a “safe collision time” t_safe = σ_min / v_max – ε, where σ_min is the smallest permissible inter‑particle distance and v_max is the maximum relative speed in the system. Any candidate collision time smaller than t_safe is flagged as potentially invalid. The algorithm also integrates this logic into the cell‑list neighbor search commonly used for hard‑sphere systems, attaching an invalid‑state flag to each cell pair to avoid unnecessary recomputation of already‑validated interactions.

Performance is evaluated on three benchmark suites. In a two‑particle test repeated for 10⁹ collisions, the new method maintains zero overlaps, whereas a conventional implementation first exhibits an overlap after roughly 10⁶ collisions. In a dense three‑dimensional hard‑sphere fluid at a packing fraction of 0.64 with 10⁶ particles, the stable algorithm reduces total runtime by about 15 % compared with the baseline and suppresses the frequency of invalid states to below 10⁻⁸ per event. A mixed‑size particle system shows similar gains, confirming that the approach scales across polydisperse ensembles.

The paper’s contributions are twofold: (1) it provides a precise mathematical characterization of the failure modes caused by finite‑precision arithmetic in event‑driven simulations, and (2) it delivers a general, efficient algorithm that detects and resolves invalid states by combining dynamic tolerance estimation, interval arithmetic, and minimal corrective collisions. These advances markedly improve the robustness of EDPD, enabling long‑time, high‑density simulations that were previously prone to numerical breakdown. The authors suggest future extensions to soft‑potential interactions, GPU‑accelerated implementations, and automated error‑tracking frameworks, which would broaden the applicability of the method to fields such as granular matter, astrophysical N‑body dynamics, and particle‑based rendering.