Table of Contents
Fetching ...

Kroneckerised Particle Mesh Ewald

Igor Chollet

TL;DR

This work introduces a Kronecker-based reformulation of the PME reciprocal part, replacing FFT-dependent structure with a Sum of Kronecker Products (SKP) to enable scalable parallel evaluation on distributed-memory architectures. It develops two SKP construction routes—a near-optimal SVD-based nearest Kronecker approximation and a quadrature-based expansion for the Erf-based truncation—along with multivariate product-grid interpolation and compressed real Fourier matrices to enable a tensorized, grid-aware representation. A parallel divide-and-conquer strategy (SPKMV) for fast Kronecker products over grids is integrated into a full Kroneckerised PME (KPME) algorithm, with a detailed complexity analysis and numerical validation showing convergence behavior and strong scalability across CPU-based HPC configurations. While sequential complexity may exceed PME in some regimes, the framework offers substantial parallel potential and is extensible with NUFFT/GPU enhancements to close the performance gap and broaden applicability to large-scale molecular dynamics simulations.

Abstract

Particle Mesh Ewald (PME) methods accelerated through Fast Fourier Transforms (FFTs) for their reciprocal part are widely used to solve N -body problems over periodic structures with Laplace-like kernels. The FFT dependence of classical PME may mitigate its performance on parallel distributed-memory architectures. We here introduce a new variant of the reciprocal part based on Sum of Kronecker Products (SKP) instead of FFT. Moreover, our implementation of this new method is not linearithmic (as opposed to classical PME) but has an important parallel potential. We present the different approximation levels exploited in our new scheme and demonstrate to what extent it could be used on parallel distributed-memory architectures. Numerical examples supplement presented assertions.

Kroneckerised Particle Mesh Ewald

TL;DR

This work introduces a Kronecker-based reformulation of the PME reciprocal part, replacing FFT-dependent structure with a Sum of Kronecker Products (SKP) to enable scalable parallel evaluation on distributed-memory architectures. It develops two SKP construction routes—a near-optimal SVD-based nearest Kronecker approximation and a quadrature-based expansion for the Erf-based truncation—along with multivariate product-grid interpolation and compressed real Fourier matrices to enable a tensorized, grid-aware representation. A parallel divide-and-conquer strategy (SPKMV) for fast Kronecker products over grids is integrated into a full Kroneckerised PME (KPME) algorithm, with a detailed complexity analysis and numerical validation showing convergence behavior and strong scalability across CPU-based HPC configurations. While sequential complexity may exceed PME in some regimes, the framework offers substantial parallel potential and is extensible with NUFFT/GPU enhancements to close the performance gap and broaden applicability to large-scale molecular dynamics simulations.

Abstract

Particle Mesh Ewald (PME) methods accelerated through Fast Fourier Transforms (FFTs) for their reciprocal part are widely used to solve N -body problems over periodic structures with Laplace-like kernels. The FFT dependence of classical PME may mitigate its performance on parallel distributed-memory architectures. We here introduce a new variant of the reciprocal part based on Sum of Kronecker Products (SKP) instead of FFT. Moreover, our implementation of this new method is not linearithmic (as opposed to classical PME) but has an important parallel potential. We present the different approximation levels exploited in our new scheme and demonstrate to what extent it could be used on parallel distributed-memory architectures. Numerical examples supplement presented assertions.
Paper Structure (26 sections, 3 theorems, 84 equations, 7 figures, 4 algorithms)

This paper contains 26 sections, 3 theorems, 84 equations, 7 figures, 4 algorithms.

Key Result

Theorem 1

Let $M\in \mathbb{N}^*$, let $\mathcal{H}_M$ be the reciprocal kernel defined in Eq. eq::sumhexpr and let $\epsilon > 0$. Then there exist $\Xi_\mathbb{Y}\subset \mathbb{R}^{3}$ (resp. $\Xi_\mathbb{X}$) a product of finite grids, $\Xi_\mathbb{Y} = \Xi_{\mathbb{Y},0}\otimes \Xi_{\mathbb{Y},1} \otim where the cardinal of $\Lambda$ denoted by $|\Lambda|$ verifies $|\Lambda| \leq (2M+1)^2$, $\mathbf

Figures (7)

  • Figure 1: Relative ranks (i.e. $\frac{r}{(2M+1)^2}$ where $r$ is the number of terms in the SKP approximation) of SKP for $\bm{\alpha}_M$ w.r.t. various values of Ewald parameter$\xi$ for $\tau(p) = erf(\xi R)$ and $M=\xi$ with SVD tolerance $10^{-8}$. In top and parentheses are indicated the maximal rank over all tested $M$ for each corresponding $\xi$ and on the right are given the final ranks for $\xi=100$ for each $M$.
  • Figure 2: For various $|\cdot |_\infty$ precisions in the cardinal sine quadrature (indicated with different colors) resulting in $N_{quad}$ points, compression rates are computed as $\frac{N_{quad}}{(2M+1)^2}$.
  • Figure 3: (Left) Schematic representation of a matrix $\mathcal{A}^{(p)}$ with $K=2$. Matrices with the same color are the same in this representation. (Right) Matrix and vector reshaping for matrix-matrix product.
  • Figure 4: Schematic view of the parallel algorithm at a given quadrature node $\lambda$ on a $3\times 3\times 3$ MPI grid. Elements of an MPI communicator are depicted in the same color (except that the red current process $\mathbf{l}$ which is in the orange communicator). (a)Each process independently applies the Fourier transform $\mathcal{V}[\lambda]_{c_{l_2}}$ to its data reshape by $\mathcal{I}^{\bullet}_2$ on the $Z$-axis of the $\lambda$'s term in SKP decomposition.(b) The information is reduced through the $Z$-axis communicators (in the Fourier domain). Then the Fourier information is locally transformed back into the real one along the $Z$-axis using $\mathcal{V}[\lambda]_{c_{l_2}}^*$. The information (after a reshape) along the $Y$-axis is transformed into the Fourier domain using $\mathcal{V}[\lambda]_{c_{l_1}}$.(c) The information is reduced through the $Y$-axis communicators in the Fourier domain. Then the Fourier information is locally transformed back into the real one along the $X$-axis using $\mathcal{V}[\lambda]_{c_{l_1}}^*$. The information (after a reshape) along the $X$-axis is transformed into the Fourier domain using $\mathcal{V}[\lambda]_{c_{l_0}}$.(d) The information is reduced through the $X$-axis communicators in the Fourier domain. Then the Fourier information is locally transformed back into the real one along the $Z$-axis using $\mathcal{V}[\lambda]_{c_{l_0}}^*$.
  • Figure 5: Convergence plot for various $M$ and cell size ratio $\nu$ (see proof of Thm. \ref{['thm::convratemultcomp']}) using the following code. Red: $M=2$; Grey: $M=4$; Blue: $M=8$; Yellow: $M=16$; Green: $M=32$; Crosses: $\nu=2^{-1}$; Squares: $\nu=4^{-1}$; Circles: $\nu=8^{-1}$(Left) Relative error on a cell w.r.t. one-dimensional interpolation order for various values of $M$(Right) Relative error on final algorithm (interpolation + quadrature using Yarvin-Rokhlin rule for $M=12$ + split compressed real Fourier matrices)
  • ...and 2 more figures

Theorems & Definitions (7)

  • Theorem 1: Existence of KPME
  • Remark 1
  • Remark 2
  • Remark 3
  • Theorem 2
  • Theorem 3
  • proof