Accelerating Fast Ewald Summation with Prolates for Molecular Dynamics Simulations

TL;DR

This work tackles the bottleneck of global FFT communication in long-range electrostatics for molecular dynamics by introducing ESP, which replaces the Gaussian-based Ewald split with a split derived from the first prolate spheroidal wave function (PSWF) and uses PSWFs for the spreading/interpolation kernels. The PSWF-based split dramatically reduces the required FFT grid size, leading to substantial improvements in strong scaling and a 2–3x reduction in total execution time on large-core runs, while maintaining accuracy at standard biological tolerances ($\Delta=10^{-4}$ to $10^{-5}$). ESP demonstrates robust performance gains in large-scale bulk water, lysozyme in solution, and a transmembrane bc1 complex, with up to an ~8x reduction in Fourier modes relative to Gaussian splits. The method is readily integrated into existing MD codes (LAMMPS and GROMACS) and extends to various boundary conditions and potentially to other interaction kernels, offering practical impact for extensive MD simulations and cross-scale modeling.

Abstract

Fast Ewald summation is the most widely used approach for computing long-range Coulomb interactions in molecular dynamics (MD) simulations. While the asymptotic scaling is nearly optimal, its performance on parallel architectures is dominated by the global communication required for the underlying fast Fourier transform (FFT). Here, we develop a novel method, ESP - Ewald summation with prolate spheroidal wave functions (PSWFs) - that, for a fixed precision, sharply reduces the size of this transform by performing the Ewald split via a PSWF. In addition, PSWFs minimize the cost of spreading and interpolation steps that move information between the particles and the underlying uniform grid. We have integrated the ESP method into two widely-used open-source MD packages: LAMMPS and GROMACS. Detailed benchmarks show that this reduces the cost of computing far-field electrostatic interactions by an order of magnitude, leading to better strong scaling with respect to number of cores. The total execution time is reduced by a factor of 2 to 3 when using more than one thousand cores, even after optimally tuning the existing internal parameters in the native codes. We validate the accelerated codes in realistic long-time biological simulations.

Figures (10)

• Figure 1: Performance comparison of ESP with PME and PPPM for bulk water systems. The total simulation time per step, averaged over five runs of two thousand steps each, is shown for four different simulation sizes: a, 3,597,693 atoms; b, 4,214,784 atoms; c, 106,238,712 atoms; and d, 192,000,000 atoms, as a function of the number of CPU cores. Panels a and c were generated using ESP implemented within LAMMPS, compared against the native PPPM option, with an error tolerance of $\Delta = 10^{-4}$. Panels b and d were generated using ESP implemented within GROMACS, compared against its native PME option, with an error tolerance of $\Delta = 10^{-5}$. Blue circles show PPPM/PME with default parameters; brown '×' marks show optimally tuned PPPM/PME; green squares show ESP. Lightly shaded markers indicate confidence intervals assessed by repeated runs.
• Figure 1: Relative error across 100 bulk water configurations randomly sampled from the equilibration state is presented. The results are obtained using implementations in a, LAMMPS with an error tolerance of $\Delta=10^{-4}$; and b, GROMACS with an error tolerance of $\Delta=10^{-5}$. Comparisons are made between the proposed ESP variants of the codes based on PSWF, and the native LAMMPS/GROMACS codes under two settings, the default parameter setup and a near-optimal, fine-tuned setup.
• Figure 2: CPU time per simulation step, averaged over five independent runs with $2000$ steps each, broken down into four components: FFT and IFFT operations (red); spreading, diagonal scaling, and interpolation (yellow); short-range pairwise interactions (light blue); and all other simulation tasks (dark blue). Panels a– d correspond to the large water systems shown in Fig. \ref{['fig:time']}a–d, respectively, each representing a different system size. Data are shown for LAMMPS/PPPM and GROMACS/PME, using both default and optimally tuned parameters, as well as for LAMMPS and GROMACS with the proposed ESP method. Speedups in the long-range interaction are annotated directly within each panel.
• Figure 2: Comparison of additional simulation results obtained using GROMACS with PME- and ESP-based MD on the same lysozyme protein system analyzed in Fig. \ref{['fig:protein']}. a, Root mean square displacement (RMSD) of the protein backbone. b, Radius of gyration of the proteins. Data are shown for both the native PME-based GROMACS implementation and the ESP method over 100 ns of simulation. Light-colored shaded areas and colored markers in a, b indicate confidence intervals based on five independent runs. All results show good agreement between the two methods.
• Figure 3: Comparison of simulation results obtained using GROMACS/PME and GROMACS/ESP on a lysozyme protein system comprising 1,036,152 atoms, including $27$ duplicated lysozyme proteins. a, Simulation snapshot of the local environment around one protein, with coloring based on secondary structure. b, Simulation performance using different methods, including GROMACS with ESP, and the default and optimally tuned setups of the native GROMACS/PME code. c, Root mean square fluctuation (RMSF) across residues; the inset shows two lysozyme proteins from the PME- and ESP-based simulations, colored by residue-specific temperature ($\beta$)-factors. d, Solvent-accessible surface area averaged over all proteins during the simulation. Light-colored shaded areas indicate the standard deviation averaged over five independent runs. e, Distribution of characteristic inter-domain distances between residue pairs C54–C97 and C54–C129. f, The total number of residues belonging to each type of secondary structure.