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

TL;DR

This paper introduces a novel, efficient, and accurate random batch sum-of-Gaussians algorithm for 3D Yukawa system molecular dynamics, improving computational speed and force calculation accuracy.

Contribution

The paper presents a new sum-of-Gaussians decomposition and an adaptive importance sampling strategy, enabling efficient, high-accuracy, energy-stable large-scale Yukawa system simulations without FFT.

Findings

Achieves O(N) complexity and high parallel scalability.

Demonstrates high accuracy and stability in large-scale simulations.

Validates method with simulations of plasma and fusion-relevant systems.

Abstract

Yukawa systems have drawn widespread interest across various applications. In this paper, we introduce a novel random batch sum-of-Gaussians (RBSOG) algorithm for molecular dynamics simulations of 3D 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(N) 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. Rigorous theoretical analyses are presented. We validate the performance of RBSOG method through simulations of one-component plasma under weak and strong coupling conditions, using up to 10^6 particles and 1024 CPU cores. As a practical application in fusion ignition, we simulate high-temperature, high-density deuterium-\alpha mixtures to study the energy exchange between deuterium and high-energy \alpha particles. The RBSOG method can be readily extended to other dielectric response functions, offering a promising approach for large-scale simulations.