Random batch sum-of-Gaussians algorithm for molecular dynamics simulations of Yukawa systems in three dimensions

Abstract

Yukawa systems have drawn widespread interest across various applications, including plasma physics, colloidal science, and astrophysics, due to their critical role in modeling electrostatic interactions. In this paper, we introduce a novel random batch sum-of-Gaussians (RBSOG) algorithm for molecular dynamics simulations of three-dimensional Yukawa systems with periodic boundary conditions. We develop a sum-of-Gaussians (SOG) decomposition of the Yukawa kernel, dividing the interactions into near-field and far-field components. The near-field component, singular but compactly supported in a local domain, is calculated directly. The far-field component, represented as a sum of smooth Gaussians, is treated using the random batch approximation in Fourier space with an adaptive importance sampling strategy to reduce the variance of force calculations. Unlike the traditional Ewald decomposition, which introduces discontinuities and significant truncation error at the cutoff, the SOG decomposition achieves high-order smoothness and accuracy near the cutoff, allowing for efficient and energy-stable simulations. Additionally, by avoiding the use of the fast Fourier transform, our method achieves optimal $O\left(N\right)$ complexity while maintaining high parallel scalability. Finally, unlike previous random batch approaches, the proposed adaptive importance sampling strategy achieves nearly optimal variance reduction across the regime of the coupling parameters, which is essential for handling varying coupling strengths across weak and strong regimes of electrostatic interactions. Rigorous theoretical analyses are presented, including SOG decomposition construction, variance estimation, and simulation convergence. We validate the performance of RBSOG method through numerical simulations of one-component plasma under weak and strong coupling conditions, using up to 10⁶ particles and 1024 CPU cores. As a practical application in fusion ignition, we simulate high-temperature, high-density deuterium-α mixtures to study the energy exchange between deuterium and high-energy α particles. Due to the flexibility of the Gaussian approximation, the RBSOG method can be readily extended to other dielectric response functions, offering a promising approach for large-scale simulations.