Monte Carlo simulations of interacting particles with fast and accurate electrostatics

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.