For the integration of the equations of motion we use a fifth order predictor-corrector MD-scheme, see Ref. [18, 19]. Since we are interested in a static situation in 2D, with almost monodisperse particles, no particle has more than six nearest neighbors. For the simulation we keep the neighbors in memory in order to reduce the computational effort.
There are two forces acting on particle i when it overlaps with particle j, i.e. when the distance . We use an elastic force
with the spring constant k, acting on particle i in normal direction . The second force in normal direction is dissipative
accounting for the inelasticity of the contacts. In Eq. 2 the constant is a phenomenological dissipation coefficient, and is the relative velocity of the particles i and j. As mentioned above, we neglect tangential forces. The contact of a particle with a wall or an immobile particle is mimicked by setting the mass of the immobile contact partner to infinity. Finally, the influence of gravity g is readily included into the equations of motion.