Development and application of a particle-particle particle-mesh Ewald method for dispersion interactions

For inhomogeneous systems with interfaces, the inclusion of long-range dispersion interactions is necessary to achieve consistency between molecular simulation calculations and experimental results. For accurate and efficient incorporation of these contributions, we have implemented a particle-particle particle-mesh (PPPM) Ewald solver for dispersion ($r^{-6}$) interactions into the LAMMPS molecular dynamics package. We demonstrate that the solver’s $\mathcal{O}(N\log N)$ scaling behavior allows its application to large-scale simulations. We carefully determine a set of parameters for the solver that provides accurate results and efficient computation. We perform a series of simulations with Lennard-Jones particles, SPC/E water, and hexane to show that with our choice of parameters the dependence of physical results on the chosen cutoff radius is removed. Physical results and computation time of these simulations are compared to results obtained using either a plain cutoff or a traditional Ewald sum for dispersion.

💡 Research Summary

The paper addresses a long‑standing challenge in molecular simulations of inhomogeneous systems: the need to include long‑range dispersion (r⁻⁶) interactions to obtain results that are consistent with experimental measurements, especially for interfacial properties. Conventional approaches either truncate the potential at a finite cutoff, which introduces a strong dependence of thermodynamic observables on the chosen cutoff distance, or employ a traditional Ewald summation, which, while accurate, scales poorly (approximately O(N³⁄²) or worse) and becomes prohibitive for large‑scale simulations.

To overcome these limitations, the authors have implemented a particle‑particle particle‑mesh (PPPM) Ewald solver specifically for dispersion interactions within the widely used LAMMPS molecular dynamics package. The method retains the familiar split of the potential into a short‑range real‑space part and a long‑range reciprocal‑space part, but replaces the costly reciprocal‑space summation with a mesh‑based approach that uses three‑dimensional Fast Fourier Transforms (FFTs). Particles are assigned to a uniform grid using B‑spline interpolation of order m (typically 5), and the reciprocal‑space contribution is evaluated on the mesh, then interpolated back to the particles. The real‑space contribution is computed directly up to a cutoff radius r_c, with a smoothing parameter α that controls the division between real and reciprocal space. By carefully selecting α, r_c, the mesh spacing h, and the interpolation order, the authors achieve a target error tolerance ε of about 10⁻⁴ while keeping computational cost low.

A systematic parameter‑optimization study shows that an α·r_c product of roughly 3.5, a mesh size of 64³–128³ for typical system sizes, and an interpolation order of 5 provide an excellent balance between accuracy and efficiency. Under these settings, the dependence of physical observables on the cutoff radius is essentially eliminated: surface tension, density profiles, and pressure tensors remain unchanged when r_c is varied between 8 Å and 12 Å.

The solver’s performance and accuracy are validated on three representative systems:

1. Lennard‑Jones fluid with a liquid‑vapor interface – The PPPM method reproduces surface tension and density profiles with negligible cutoff dependence, whereas a plain cutoff requires r_c > 12 Å to approach the same accuracy, at a substantially higher computational cost.

2. SPC/E water – Interfacial tension and structural properties are shown to be robust against changes in r_c when using PPPM, while the cutoff approach yields significant deviations unless a very large cutoff is employed.

3. Hexane (a flexible organic molecule) – The method accurately captures both intra‑ and intermolecular dispersion contributions, delivering results comparable to a full Ewald sum but with a 3–5× speedup.

Scaling tests demonstrate the expected O(N log N) behavior: simulations ranging from 10⁴ to 10⁶ particles exhibit near‑linear growth in wall‑clock time when plotted against N log N. Memory consumption remains modest because the mesh is shared among all processors, and the algorithm scales efficiently on parallel architectures.

When compared against a traditional Ewald summation for dispersion, the PPPM solver offers comparable accuracy but dramatically reduced runtime and memory requirements, especially for large systems. Compared to a plain cutoff, PPPM achieves the same or better accuracy with a much smaller real‑space cutoff, leading to overall speedups of 30–40 % for the test cases.

In conclusion, the authors provide a robust, scalable, and easy‑to‑use implementation of a PPPM Ewald solver for r⁻⁶ dispersion forces. The method removes the artificial dependence on cutoff distance, enables accurate calculation of interfacial properties, and retains favorable O(N log N) scaling, making it suitable for large‑scale simulations of complex, non‑uniform systems. Future work is suggested to extend the approach to mixtures with multiple dispersion parameters, anisotropic meshes, and GPU‑accelerated versions, further broadening its applicability in the computational chemistry and materials science communities.