Coulomb interaction, following an inverse-square force-law, quantifies the amount of force between two stationary, electrically charged particles. The long-range nature of Coulomb interactions poses a major challenge to molecular dynamics simulations
which are major tools for problems at the nano-/micro- scale. Various algorithms aim to speed up the pairwise Coulomb interactions to a linear scaling but the poor scalability limits the size of simulated systems. Here, we conduct an efficient molecular dynamics algorithm with the random batch Ewald method on all-atom systems where the complete Fourier components in the Coulomb interaction are replaced by randomly selected mini batches. By simulating the N-body systems up to 100 million particles using 10 thousand CPU cores, we show that this algorithm furnishes O(N) complexity, almost perfect scalability and an order of magnitude faster computational speed when compared to the existing state-of-the-art algorithms. Further examinations of our algorithm on distinct systems, including pure water, micro-phase-separated electrolyte and protein solution demonstrate that the spatiotemporal information on all time and length scales investigated and thermodynamic quantities derived from our algorithm are in perfect agreement with those obtained from the existing algorithms. Therefore, our algorithm provides a breakthrough solution on scalability of computing the Coulomb interaction. It is particularly useful and cost-effective to simulate ultra-large systems, which was either impossible or very costing to conduct using existing algorithms, thus would benefit the broad community of sciences.
We propose a fast method for the calculation of short-range interactions in molecular dynamics simulations. The so-called random-batch list method is a stochastic version of the classical neighbor-list method to avoid the construction of a full Verle
t list, which introduces two-level neighbor lists for each particle such that the neighboring particles are located in core and shell regions, respectively. Direct interactions are performed in the core region. For the shell zone, we employ a random batch of interacting particles to reduce the number of interaction pairs. The error estimate of the algorithm is provided. We investigate the Lennard-Jones fluid by molecular dynamics simulations, and show that this novel method can significantly accelerate the simulations with a factor of several fold without loss of the accuracy. This method is simple to implement, can be well combined with any linked cell methods to further speed up and scale up the simulation systems, and can be straightforwardly extended to other interactions such as Ewald short-range part, and thus it is promising for large-scale molecular dynamics simulations.
We implement two recently developed fast Coulomb solvers, HSMA3D [J. Chem. Phys. 149 (8) (2018) 084111] and HSMA2D [J. Chem. Phys. 152 (13) (2020) 134109], into a new user package HSMA for molecular dynamics simulation engine LAMMPS. The HSMA package
is designed for efficient and accurate modeling of electrostatic interactions in 3D and 2D periodic systems with dielectric effects at the O(N) cost. The implementation is hybrid MPI and OpenMP parallelized and compatible with existing LAMMPS functionalities. The vectorization technique following AVX512 instructions is adopted for acceleration. To establish the validity of our implementation, we have presented extensive comparisons to the widely used particle-particle particle-mesh (PPPM) algorithm in LAMMPS and other dielectric solvers. With the proper choice of algorithm parameters and parallelization setup, the package enables calculations of electrostatic interactions that outperform the standard PPPM in speed for a wide range of particle numbers.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel functions, and develop a fast algorithm for convolution quadrature based on the SOE, which allows an order $N$ calculation for $N$ time steps of approximat
ing a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallee-Poussin sums for a semi-analytical exponential expansion of a general kernel, and a model reduction technique for the minimization of the number of exponentials under given error tolerance. We employ the SOE expansion for the finite part of the splitting convolution kernel such that the convolution integral can be solved as a system of ordinary differential equations due to the exponential kernels. The remaining part is explicitly approximated by employing the generalized Taylor expansion. The significant features of our algorithm are that the SOE method is efficient and accurate, and works for general kernels with controllable upperbound of positive exponents. We provide numerical analysis for the SOE-based convolution quadrature. Numerical results on different kernels, the convolution integral and integral equations demonstrate attractive performance of both accuracy and efficiency of the proposed method.
Approximation of interacting kernels by sum of Gaussians (SOG) is frequently required in many applications of scientific and engineering computing in order to construct efficient algorithms for kernel summation or convolution problems. In this paper,
we propose a kernel-independent SOG method by introducing the de la Vallee-Poussin sum and Chebyshev polynomials. The SOG works for general interacting kernels and the lower bound of Gaussian bandwidths is tunable and thus the Gaussians can be easily summed by fast Gaussian algorithms. The number of Gaussians can be further reduced via the model reduction based on the balanced truncation based on the square root method. Numerical results on the accuracy and model reduction efficiency show attractive performance of the proposed method.
We propose a harmonic surface mapping algorithm (HSMA) for electrostatic pairwise sums of an infinite number of image charges. The images are induced by point sources within a box due to a specific boundary condition which can be non-periodic. The HS
MA first introduces an auxiliary surface such that the contribution of images outside the surface can be approximated by the least-squares method using spherical harmonics as basis functions. The so-called harmonic surface mapping is the procedure to transform the approximate solution into a surface charge and a surface dipole over the auxiliary surface, which becomes point images by using numerical integration. The mapping procedure is independent of the number of the sources and is considered to have a low complexity. The electrostatic interactions are then among those charges within the surface and at the integration points, which are all the form of Coulomb potential and can be accelerated straightforwardly by the fast multipole method to achieve linear scaling. Numerical calculations of the Madelung constant of a crystalline lattice, electrostatic energy of ions in a metallic cavity, and the time performance for large-scale systems show that the HSMA is accurate and fast, and thus is attractive for many applications.
