The simulation of classical particle systems by means of molecular dynamics techniques requires the evaluation of potentials and forces between particles to explore the phase or configuration space of the system. The interactions can be classified into short-range and long-range contributions. While short-range potentials are evaluated very efficiently by means of neighbor list techniques, which reduce the computational complexity to O(N), the long range interactions must be evaluated over all particle pair contributions in the system. This increases the complexity to O(N2), limiting very often the tractable system size to a few thousand particles. In charged or polarized systems, the evaluation of Coulomb interactions between particles is the dominant part of the computation. To tame the quadratic computational complexity of the problem, a number of different methods with O(N) and O(N log(N)) have been developed over the past. Prominent examples are the linearly scaling fast multipole method and FFT-based methods, like P3M, which exhibit an N log(N) scaling but also tree methods, originally designed for strongly inhomogeneous gravitational problems as well as variants of the multigrid method. To integrate an efficient, yet accurate method into an existing simulation code is often a time consuming task, because of the level of implementation complexity. This is even more pronounced when considering scalability on parallel computers. Furthermore, different types of boundary conditions might be necessary to consider, which are not implemented in a standard way even in community codes. To tackle large particle systems, solvers of low computational complexity and good scalability have to be considered. In order to facilitate the integration of such high level implementations of Coulomb solvers into existing programs, the scalable library ScaFaCoS was developed [1], which contains various methods for treating long range interactions in complex particle systems under various boundary conditions, e.g., open boundaries or periodic boundaries in various cartesian directions [2]. Implemented algorithms comprise the fast multipole method, Barnes–Hut tree method, P3M, multigrid methods, and a newly developed Maxwell solver MEMD. Several methods allow full error control and the optimal adjustment of method parameters to reduce run time. An overview will be given about the methods, their accuracy and stability, as well as results for performance and scalability on parallel clusters. REFERENCES 1. http://www.scafacos.de 2. Arnold, A., Fahrenberger, F. Holm, Chr., Lenz, O., Bolten, M., Dachsel, H., Halver,. R., Kabadshow, I., Gähler, F., Heber, F., Iseringhausen, J., Hofmann, M., Pippig, M., Potts, D., Sutmann, G. Comparison of scalable fast methods for long-range interactions. Phys. Rev. E. 2013, 88, 063308