Fast Algorithm for Quasi-2D Coulomb Systems

TL;DR

Quasi-2D Coulomb systems with partial periodicity pose significant computational challenges due to long-range $1/r$ interactions under slab geometry. The paper introduces RBSE2D, a framework that pairs a sum-of-exponentials (SOE) approximation in the non-periodic $z$-direction with random batch importance sampling in the periodic $xy$-directions, achieving deterministic $O(N^{7/5})$ complexity and, with stochastic batching, an effective $O(N)$ per-step cost. Rigorous error analyses under Debye–Hückel theory, along with comprehensive numerical tests, show substantial accuracy and practical speedups of $2$–$3$ orders of magnitude over conventional Ewald2D methods, enabling MD with up to $10^6$ particles on a single core. The RBSE2D approach is kernel-independent (via kernel-independent SOE, e.g., VPMR) and mesh-free, offering a versatile tool for large-scale simulations of confined Coulomb systems and a blueprint for extending to other partially-periodic lattice kernels. Future directions include handling dielectric mismatch, extending to other ensembles, and leveraging CPU/GPU parallelism to push to even larger scales.

Abstract

Quasi-2D Coulomb systems are of fundamental importance and have attracted much attention in many areas nowadays. Their reduced symmetry gives rise to interesting collective behaviors, but also brings great challenges for particle-based simulations. Here, we propose a novel algorithm framework to address the $O(N^2)$ simulation complexity associated with the long-range nature of Coulomb interactions. First, we introduce an efficient Sum-of-Exponentials (SOE) approximation for the long-range kernel associated with Ewald splitting, achieving uniform convergence in terms of inter-particle distance, which reduces the complexity to $O(N^{7/5})$. We then introduce a random batch sampling method in the periodic dimensions, the stochastic approximation is proven to be both unbiased and with reduced variance via a tailored importance sampling strategy, further reducing the computational cost to $O(N)$. The performance of our algorithm is demonstrated via various numerical examples. Notably, it achieves a speedup of $2\sim 3$ orders of magnitude comparing with Ewald2D method, enabling molecular dynamics (MD) simulations with up to $10^6$ particles on a single core. The present approach is therefore well-suited for large-scale particle-based simulations of Coulomb systems under confinement, making it possible to investigate the role of Coulomb interaction in many practical situations.

Figures (7)

• Figure 1: The absolute error of the SOE expansion for (a) $\xi^{\pm}(k,z)$ and (b) $\mathrm{erf}(\alpha z)$ is plotted as a function of $z$, while fixing $k=\alpha=1$; absolute error of the SOE expansion of (c) $\xi^{+}(k,z)$ and (d) $\xi^{-}(k,z)$ as a function of $k^2$, while fixing $z=1$. Data are presented for SOEs with varying numbers of exponentials, $M=4,$ 8 and 16.
• Figure 2: Accuracy in the electrostatic energy by the SOEwald2D method. (a): absolute error as a function of $s$; (b): relative error as a function of total number of ions $N$ with fixed ion density $\rho_{s}$. Results with different number of exponentials $M$ are considered.
• Figure 3: The absolute error in electrostatic energy is evaluated for the SOEwald2D method using three sets of parameters, as well as for the Ewald2D method with $s=5$, as a function of the system's thickness $L_z$.
• Figure 4: The absolute error in the pointwise electrostatic forces calculated using the SOEwald2D versus particles' $z$-coordinates. Two different scenarios are considered: (a) uniformly distributed 50 anions and 50 cations and (b) uniformly distributed 100 cations with surface charge densities $\sigma_{\mathrm{top}}= \sigma_{\mathrm{bot}} = -0.005$.
• Figure 5: (a) The concentration of cations along $z$, with subplot indicating the convergence in the relative error of the average electrostatic energy as a function of batch size $P$; (b) and (c) the MSD profiles in $xy$ and $z$ against time for a $1:1$ electrolyte confined by neutral slabs. Results by using different batch sizes $P=20, 30, 60, 120$ are shown.

Theorems & Definitions (31)

• Theorem 2.1
• proof
• Proposition 2.2
• proof
• Definition 2.3
• Lemma 2.4
• Theorem 2.5
• Proposition 2.6
• Remark 2.7
• Theorem 2.8