  To minimise systematic errors in Monte Carlo simulations of charged
particles, long range electrostatic interactions have to be calculated
accurately and efficiently. Standard approaches, such as Ewald summation or the
naive application of the classical Fast Multipole Method, result in a cost per
Metropolis-Hastings step which grows in proportion to some positive power of
the number of particles $N$ in the system. This prohibitively large cost
prevents accurate simulations of systems with a sizeable number of particles.
Currently, large systems are often simulated by truncating the Coulomb
potential which introduces uncontrollable systematic errors. In this paper we
present a new multilevel method which reduces the computational complexity to
$\mathcal{O}(\log(N))$ per Metropolis-Hastings step, while maintaining errors
which are comparable to direct Ewald summation. We show that compared to
related previous work, our approach reduces the overall cost by better
balancing time spent in the proposal- and acceptance- stages of each
Metropolis-Hastings step. By simulating large systems with up to $N=10^5$
particles we demonstrate that our implementation is competitive with
state-of-the-art MC packages and allows the simulation of very large systems of
charged particles with accurate electrostatics.
