For short-range molecular dynamics (MD) simulations with heterogeneous particle distributions that dynamically change over time, highly flexible and dynamical load balancing methods are mandatory to achieve good parallel scalability. Designing and implementing load balancing algorithms is complex, especially for existing applications which were not designed to support arbitrary domain decompositions. In this paper, we present our approach to incorporate general domain decompositions and dynamic re-balancing into the existing MD software package ESPResSo. We describe the relevant interfaces and abstractions which enable us to reuse the physics algorithms in ESPResSo, without major re-implementations. As proof-of-concept, we show the implementation of a domain decomposition based on space-filling curves and a dynamic re-balancing mechanism using an enhanced version of the p4est library. The results indicate that our load balancing mechanism is capable of reducing the imbalance amongst processes and the total runtime of simulations in simple and complex scenarios. At the same time, the implementation of models and solvers in ESPResSo remains largely unchanged.