Fast Ewald Summation with Prolates for Charged Systems in the NPT Ensemble

TL;DR

The paper addresses the challenge of efficiently and accurately computing long-range electrostatics and the associated pressure tensor in NPT molecular dynamics simulations of charged, periodic systems. It introduces ESP with prolate spheroidal wave functions (PSWFs) for both splitting and spreading, enabling a spectrally accurate, scalable mesh-Ewald-like method that reduces Fourier-space overhead. A unified, NPT-consistent pressure formulation is derived for isotropic to fully flexible cells, including corrections for non-neutral systems, and a four-step FFT-based algorithm is proposed that requires only a single forward FFT for the long-range pressure with diagonal scaling. The method is implemented in LAMMPS and GROMACS, validated on bulk water, LiTFSI ionic liquids, and a transmembrane bc1 complex, and demonstrates strong scaling and substantial reductions in Fourier-grid requirements, yielding 2–3x speedups in large-scale NPT simulations while preserving thermodynamic accuracy. Overall, ESP provides a practical, high-accuracy, and scalable solution for pressure-sensitive simulations of charged systems in the NPT ensemble.

Abstract

We present an NPT extension of Ewald summation with prolates (ESP), a spectrally accurate and scalable particle-mesh method for molecular dynamics simulations of periodic, charged systems. Building on the recently introduced ESP framework, this work focuses on rigorous and thermodynamically consistent pressure/stress evaluation in the isothermal--isobaric ensemble. ESP employs prolate spheroidal wave functions as both splitting and spreading kernels, reducing the Fourier grid size needed to reach a prescribed pressure accuracy compared with current widely used mesh-Ewald methods based on Gaussian splitting and B-spline spreading. We derive a unified pressure-tensor formulation applicable to isotropic, semi-isotropic, anisotropic, and fully flexible cells, and show that the long-range pressure can be evaluated with a single forward FFT followed by diagonal scaling, whereas force evaluation requires both forward and inverse transforms. We provide production implementations in LAMMPS and GROMACS and validate pressure and force accuracy on bulk water, LiTFSI ionic liquids, and a transmembrane system. Benchmarks on up to $3\times 10^3$ CPU cores demonstrate strong scaling and reduced communication cost at matched accuracy, particularly for NPT pressure evaluation.

Figures (9)

• Figure 1: Relative $L_2$ error of the diagonal components of $\bm{P}_{\mathrm{ins}}$, computed over 100 equilibrium configurations, as a function of the per-dimension grid size $I_d$ (total grid points $I_{\mathcal{F}}=I_d^3$). Error tolerances are (a)$\Delta=10^{-4}$ and (b)$\Delta=10^{-6}$. Results are shown for several spreading orders $P$. The dashed line shows the fitted convergence rate, $O(e^{-I_d/2.5})$.
• Figure 2: Relative $L_2$ error of the off-diagonal components of $\bm{P}_{\mathrm{ins}}$, computed over 100 equilibrium configurations, as a function of the per-dimension grid size $I_d$. Error tolerances are (a)$\Delta=10^{-4}$ and (b)$\Delta=10^{-6}$. Results are shown for several spreading orders $P$. The dashed line shows the fitted convergence rate, $O(e^{-I_d/2.5})$.
• Figure 3: Relative $L_2$ error of the (a) diagonal and (b) off-diagonal components of $\bm{P}_{\mathrm{ins}}$, computed over 100 equilibrium configurations, as a function of the spreading order $P$. Results are shown for different tolerances $\Delta$. The dashed line shows the fitted convergence rate, $O(\mathrm{erfc}(2.05\sqrt{P}))$.
• Figure 4: Heatmaps of the relative error in the (a) diagonal and (b) off-diagonal components of the pressure tensor. The $x$- and $y$-axes denote the splitting tolerance $\Delta_{\mathrm{split}}$ and spreading tolerance $\Delta_{\mathrm{spread}}$, respectively. White circles mark the "optimal" pairs, defined as the largest $\Delta_{\mathrm{spread}}$ that achieves an error within $1.01\times$ the minimum error attained at that $\Delta_{\mathrm{split}}$.
• Figure 5: Relative $L_2$ error of the instantaneous pressure tensor versus the smallest inverse mesh volume $1/(h_xh_yh_z)$ required to reach the specified error level. For ESP, circles and squares denote the errors of the diagonal and off-diagonal components, respectively. For PPPM, crosses and plus signs denote the corresponding diagonal and off-diagonal errors.