Resistive random access memories (ReRAM) are expected to play a prominent role in modern computer architectures due to their low cost, simple structure, and unique functionality. The long-range atomic movements inside these devices, which occur over extended timescales under applied fields, can be accurately described by Dynamical Monte Carlo (DMC) simulations. In DMC, the continuum movements of atoms are discretized into ‘events’ on an atomistic graph, which is time-stepped under the influence of external fields (potential, Joule heating). Parallelization can only occur within each step, rendering such simulations highly sensitive to data movement. Here, we present a scalable DMC code that simultaneously optimizes the different computational kernels found in the field solvers (systems of linear equations, matrix-vector multiplication) and event selection (prefix sums). Our implementation leverages preconditioned sparse iterative solvers, graph-based domain decomposition to divide work between nodes, and hybrid CPU-GPU computations to optimize node usage and data transfer in a distributed environment. The acceleration ultimately enables the first investigation of ReRAM crossbar arrays with an atomistic resolution, providing deeper insights into the operating mechanisms of these devices and paving the way for their mainstream adaptation in future memory technologies.