Symmetry-Preserving Random Batch Ewald (SRBE)

• SRBE is a randomized long-range electrostatics algorithm that splits Fourier modes into a deterministic low-frequency sector and a stochastic high-frequency sector.
• The method enforces a symmetry-preserving mean-field condition to reduce charge and force fluctuations in anisotropic slab geometries by careful parameter tuning (M and P).
• SRBE achieves O(N) scaling without global FFTs, significantly enhancing performance in simulations of electrochemical systems, such as supercapacitors and membranes.

Searching arXiv for the target paper and closely related context. Symmetry-Preserving Random Batch Ewald (SRBE) is a randomized long-range electrostatics algorithm for constant-potential molecular dynamics in electrochemical systems, introduced for quasi-2D slab geometries with periodicity in $x,y$ and a free boundary in $z$ (Gao et al., 29 Sep 2025). Its defining feature is a split treatment of reciprocal-space modes: low-frequency modes along the direction perpendicular to electrodes are evaluated deterministically, while the remaining Fourier components are approximated by random batch sampling. In the reported formulation, this construction is designed to reduce charge and force fluctuations, satisfy a symmetry-preserving mean field condition in anisotropic systems with large aspect ratios, and remove the FFT-centered scalability bottleneck that limits conventional constant-potential simulations (Gao et al., 29 Sep 2025).

1. Electrostatic formulation in slab geometry

The method is formulated for a slab system in which the total electrostatic energy is written as

$$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$

As in classical Ewald summation, the Coulomb kernel is decomposed through

$\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$

This yields a real-space short-range contribution and a smooth Fourier-space long-range contribution. After zero-padding to a fully periodic box of volume $V=L_xL_yL_z$, the long-range term takes the form

$$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$

Here $S_1(\mathbf k)$ and $S_2(\mathbf k)$ are the structure factors associated with ions and electrode charges. The short-range contribution is described as rapidly decaying, while the reciprocal-space part carries the smooth long-range electrostatics.

The SRBE construction acts entirely within this Fourier representation. Rather than treating all reciprocal modes uniformly, it isolates a subset of low-frequency modes that are geometrically distinguished by the slab anisotropy. This suggests that the method is tailored not merely to accelerate Ewald summation, but to preserve the physically relevant one-dimensional long-range behavior along the nonperiodic direction.

2. Deterministic low-frequency sector and stochastic high-frequency sector

The central decomposition in SRBE is applied to the Fourier sum itself. A low-frequency subset is defined by

$$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$

where $M$ is chosen so that $z$0 captures all “dangerous” low-lying $z$1-modes, typically

$z$2

All remaining $z$3 are treated stochastically (Gao et al., 29 Sep 2025).

The deterministic part requires no FFT. Each mode in $z$4 is summed directly at $z$5 cost, and because $z$6 for fixed aspect ratio, the total cost of this part is $z$7. The complementary high-frequency remainder is

$z$8

This remainder is approximated by an importance-sampling mini-batch of size $z$9, with samples

$$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$0

In practice, a Metropolis–Hastings sampler is used: one proposes $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$1 from a discrete Gaussian neighborhood, forms $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$2, and accepts with the probability given in Eq. (25) of the paper. The resulting estimator is

$$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$3

A common misconception is that random-batch treatment of reciprocal space can ignore geometric anisotropy without qualitative consequences. In the reported slab setting, SRBE is built on the opposite premise: the low-frequency $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$4-axis modes are exceptional and must be protected deterministically.

3. Symmetry-preserving mean field condition

The method’s “symmetry-preserving” designation refers to a specific mean-field requirement in anisotropic slab systems. For the pairwise Fourier force, exact slab-averaging in $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$5 demands

$$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$6

A random-batch scheme that neglects the $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$7 modes violates this condition and yields large variance in $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$8 (Gao et al., 29 Sep 2025). SRBE enforces the condition exactly by deterministically summing the subset $$U_c \;=\;\frac12\int_{\Omega}\rho(\mathbf r')\,d\mathbf r'\;\int_{\mathbb R^2\times[0,H]}\frac{\tilde\rho(\mathbf r'')}{|\mathbf r'-\mathbf r''|}\,d\mathbf r''.$$9 of all $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$0 modes and random-batching only on $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$1. The paper characterizes this as “screening” of low-frequency $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$2-modes, which eliminates the dominant source of variance in both charges $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$3 and forces $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$4.

The variance analysis is stated under a Debye–Hückel assumption. The reported bounds are

$\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$5

and

$\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$6

with

$\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$7

Once $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$8, these variances decay exponentially in $\frac1r \;=\;\frac{\erf(\alpha r)}r\;+\;\frac{\erfc(\alpha r)}r, \qquad \alpha>0.$9 and are $V=L_xL_yL_z$0 in the batch size. This suggests that the parameter $V=L_xL_yL_z$1 is not merely an implementation cutoff; it is the control parameter that restores the correct slab mean field.

4. Algorithmic workflow and implementation

The reported per-step data structures are the positions and velocities of ions $V=L_xL_yL_z$2, electrode Gaussian charges $V=L_xL_yL_z$3, and the structure factors $V=L_xL_yL_z$4 for the current configuration. The workflow given as Algorithm 1 is:

1. Given $V=L_xL_yL_z$5, choose batch size $V=L_xL_yL_z$6 and low-frequency cutoff $V=L_xL_yL_z$7.
2. Sample $V=L_xL_yL_z$8 on $V=L_xL_yL_z$9.
3. Compute the long-range charge-vector entries

$$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$0

1. Update electrode charges $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$1.
2. For each ion $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$2, compute the short-range force by neighbor list; compute the long-range force by splitting $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$3 deterministically and $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$4 stochastically in one batch of size $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$5.
3. Integrate positions and velocities with a Langevin or Verlet scheme.
4. Loop to the next timestep.

Parallelization is described as avoiding a global FFT altogether (Gao et al., 29 Sep 2025). Structure-factor sums over $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$6 are local reductions on each MPI rank. Communication is limited to all-reduce of $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$7 scalars and the small $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$8 contributions to charges and forces. Work per rank remains $$U_{\mathcal F} = \frac1{2V}\sum_{\mathbf k\ne\mathbf0} \frac{4\pi}{k^2} \bigl|S_1(\mathbf k)+S_2(\mathbf k)\bigr|^2 e^{-k^2/(4\alpha^2)}.$$9 with only $S_1(\mathbf k)$0 extra communication. In that sense, SRBE is an algorithmic response to a communication pattern as much as to a floating-point cost model.

5. Complexity and measured performance

The computational comparison drawn in the paper places SRBE against conventional PPPM/FFT methods. Those methods are reported to cost $S_1(\mathbf k)$1 per step and to incur $S_1(\mathbf k)$2 or worse communication in 3D pencil- or slab-decompositions. By contrast, the SRBE cost decomposition is:

• near-field neighbor list: $S_1(\mathbf k)$3;
• low-frequency modes: $S_1(\mathbf k)$4 with $S_1(\mathbf k)$5;
• high-frequency batch: $S_1(\mathbf k)$6 with $S_1(\mathbf k)$7 for fixed density and target variance.

The resulting overall complexity is $S_1(\mathbf k)$8 (Gao et al., 29 Sep 2025).

The performance measurements are summarized below.

Test	SRBE result	PPPM result
Weak scaling, 3776 particles per core, up to 343 cores	long-range part maintains $S_1(\mathbf k)$9 efficiency	degrades to $S_2(\mathbf k)$0
Strong scaling, fixed $S_2(\mathbf k)$1 particles, up to 512 cores	two orders of magnitude faster on the long-range portion	drops below 1% efficiency by 512 cores
CPU time vs. $S_2(\mathbf k)$2 on 64 cores	$S_2(\mathbf k)$3	superlinear cost due to FFT overhead

These measurements are specifically reported for the long-range part in the scaling studies. The article’s broader significance lies in showing that the removal of the global FFT is accompanied by a communication profile that remains favorable at high core counts.

6. Numerical validation, limitations, and application domain

The validation system is a coarse-grained ionic liquid $S_2(\mathbf k)$4 comprising 320 molecules between graphene electrodes containing 3136 C atoms, in a box roughly $S_2(\mathbf k)$5Å (Gao et al., 29 Sep 2025). The reported numerical experiments address both equilibrium and dynamical behavior.

For accuracy, Table 1 compares equilibrium negative-electrode charge mean and standard deviation for truncation parameters $S_2(\mathbf k)$6 and batch sizes $S_2(\mathbf k)$7. Only $S_2(\mathbf k)$8 reproduces the PPPM reference, reported as $S_2(\mathbf k)$9. The case $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$0, identified as classical RBE, fails entirely. Time-series data for atomic charges show that SRBE with $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$1 tracks PPPM, whereas smaller $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$2 produces wild fluctuations. Equilibrium energy per atom versus $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$3 shows that $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$4 suffices for $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$5 accuracy. Ion-density profiles under $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$6, $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$7, and $$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$8 match PPPM, including electric double-layer peaks.

These results delimit both the method’s strengths and its tuning requirements. One must choose the truncation parameter

$$\mathcal I \;=\;\Bigl\{(0,0,k_z)\;:\;k_z=2\pi m_z/L_z,\;1\le|m_z|\le M\Bigr\},$$9

to suppress low-$M$0 variance. Batch size $M$1 trades variance against cost; in practice, the paper states that $M$2 is enough for ionic densities $M$3. A common misconception would be to attribute SRBE’s performance solely to stochastic approximation. The reported failures at small $M$4, especially $M$5, indicate that symmetry preservation is the decisive condition in slab geometries.

The stated application range includes large-scale constant-potential molecular dynamics of supercapacitors, batteries, and electrocatalysts, as well as a hierarchy of anisotropic Coulomb systems including nanopores, membranes, and 2D materials. The paper also identifies adaptation to polarizable models and multiscale grand-canonical charge algorithms as potential applications. A plausible implication is that SRBE is best understood as a geometry-aware reciprocal-space sampling strategy for constant-potential electrochemical simulation, rather than as a generic randomized replacement for Ewald summation.