Random batch sum-of-Gaussians method for molecular dynamics simulation of particle systems in the NPT ensemble

Abstract

In this work, we develop a random batch sum-of-Gaussians (RBSOG) method for molecular dynamics simulations of charged systems in the isothermal-isobaric (NPT) ensemble. We introduce an SOG splitting of the pressure-related $1/r^3$ kernel, yielding a smooth short-/long-range decomposition for instantaneous pressure evaluation. The long-range part is treated in Fourier space by random-batch importance sampling. Because the radial and non-radial pressure components favor different proposals, direct sampling either increases structure-factor evaluations and communication or leads to substantial variance inflation. To address this tradeoff, we introduce a measure-recalibration strategy that reuses Fourier modes drawn from the radial proposal and corrects them for the non-radial target, producing an unbiased pressure estimator with significantly reduced variance and negligible extra cost. The resulting method mitigates pressure artifacts caused by cutoff discontinuities in traditional Ewald-based treatments while preserving near-optimal $O(N)$ complexity. We provide theoretical evidence on pressure decomposition error, consistency of stochastic approximation, and convergence of RBSOG-based MD. Numerical experiments on bulk water, LiTFSI ionic liquids, and DPPC membranes show that RBSOG accurately reproduces key structural and dynamical observables with small batch sizes ($P\sim 100$). In large-scale benchmarks up to $10^7$ atoms on $2048$ CPU cores, RBSOG achieves about an order-of-magnitude speedup over particle-particle particle-mesh in electrostatic calculations for NPT simulations, together with a consistent $4\times$ variance reduction relative to random batch Ewald and excellent weak/strong scalability. Overall, RBSOG provides a practical and scalable route to reduce time-to-solution and communication cost in large-scale NPT simulations.

Figures (8)

• Figure 1: Relative errors of the SOG decomposition introduced in Section \ref{['sec::Pressure']} for an SPC/E bulk water system. Panels (A) and (B) show the relative errors of the diagonal and off-diagonal components of the instantaneous pressure tensor, respectively, as functions of the number of Gaussian terms $\widetilde{M}$. The dashed lines in (A)-(B) indicate a fitted decay rate of $\mathcal{O}(\widetilde{b}^{-3\widetilde{M}})$.
• Figure 2: Radial distribution functions of (A) O-O and (B) O-H pairs in the bulk water system. PPPM results are shown as the benchmark (solid line). Results from RBSOG with the batch size $P = 128, \ 256$ and RBE with the batch size $P = 512, \ 1024$ are shown with different markers.
• Figure 3: Diffusion coefficient and shear viscosity of bulk water in the NPT ensemble. Results from RBSOG ($P=128,256$) are compared with those from RBE ($P=512,1024$) and PPPM ($\Delta=10^{-5}$). Each quantity is averaged over multiple independent simulations; solid circles indicate different data points, and error bars denote one standard deviation.
• Figure 4: Radial distribution functions (RDFs) of Li-OW (water oxygen) and Li-N pairs in LiTFSI solutions at different concentrations. Panels (A) and (B) show the PPPM reference results computed with $\Delta=10^{-5}$. Panels (C)-(J) compare the results of random batch (RB)-type methods relative to PPPM: Li-OW RDFs in (C)-(F) and Li-N RDFs in (G)-(J), with RBSOG using $P=128,256$ and RBE using $P=512,1024$. Insets show magnified views of the peak regions.
• Figure 5: Time evolution of (A) the interbilayer thickness and (B) the area of the simulation box projected onto the $X$--$Y$ plane in DPPC membrane simulations. Data are shown for RBSOG with batch size $P=256$, compared with RBE at $P=256,1024$ and the reference PPPM results.

Theorems & Definitions (12)

• Theorem 1
• proof
• Theorem 2
• Proposition 3
• Lemma 4
• Theorem 5
• proof
• Remark 6
• Proposition 7
• Theorem 8