Calculation of long-range electrostatic interactions has been a major challenge in improving the speed of molecular dynamics (MD) simulation, especially for periodic systems . While established algorithms such as the particle-mesh Ewald summation (PME) method are widely adopted in simulation softwares, developing new algorithms is needed as the size of simulation systems continues to grow, and also new computing hardware architectures emerge. The popular PME algorithm scales with the number of atoms N as O(N ln N) , and it utilizes the fast Fourier transform method that can pose significant communication bottleneck on multi-processor or multi-GPU systems . On the other hand, while the Fast Multipole Method (FMM) scales approximately as O(N) , it comes with a high cost of calculating prefactors for the multipole expansion.