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 develop a random batch Ewald (RBE) method for molecular dynamics simulations of particle systems with long-range Coulomb interactions, which achieves an $O(N)$ complexity in each step of simulating the $N$-body systems. The RBE method is based on
the Ewald splitting for the Coulomb kernel with a random mini-batch type technique introduced to speed up the summation of the Fourier series for the long-range part of the splitting. Importance sampling is employed to reduce the induced force variance by taking advantage of the fast decay property of the Fourier coefficients. The stochastic approximation is unbiased with controlled variance. Analysis for bounded force fields gives some theoretic support of the method. Simulations of two typical problems of charged systems are presented to illustrate the accuracy and efficiency of the RBE method in comparison to the results from the Debye-Huckel theory and the classical Ewald summation, demonstrating that the proposed method has the attractiveness of being easy to implement with the linear scaling and is promising for many practical applications.
We propose a fast potential splitting Markov Chain Monte Carlo method which costs $O(1)$ time each step for sampling from equilibrium distributions (Gibbs measures) corresponding to particle systems with singular interacting kernels. We decompose the
interacting potential into two parts, one is of long range but is smooth, and the other one is of short range but may be singular. To displace a particle, we first evolve a selected particle using the stochastic differential equation (SDE) under the smooth part with the idea of random batches, as commonly used in stochastic gradient Langevin dynamics. Then, we use the short range part to do a Metropolis rejection. Different from the classical Langevin dynamics, we only run the SDE dynamics with random batch for a short duration of time so that the cost in the first step is $O(p)$, where $p$ is the batch size. The cost of the rejection step is $O(1)$ since the interaction used is of short range. We justify the proposed random-batch Monte Carlo method, which combines the random batch and splitting strategies, both in theory and with numerical experiments. While giving comparable results for typical examples of the Dyson Brownian motion and Lennard-Jones fluids, our method can save more time when compared to the classical Metropolis-Hastings algorithm.
We develop a modified Poisson-Nernst-Planck model which includes both the long-range Coulomb and short-range hard-sphere correlations in its free energy functional such that the model can accurately describe the ion transport in complex environment a
nd under a nanoscale confinement. The Coulomb correlation including the dielectric polarization is treated by solving a generalized Debye-Huckel equation which is a Greens function equation with the correlation energy of a test ion described by the self Greens function. The hard-sphere correlation is modeled through the modified fundamental measure theory. The resulting model is available for problems beyond the mean-field theory such as problems with variable dielectric media, multivalent ions, and strong surface charge density. We solve the generalized Debye-Huckel equation by the Wentzel-Kramers-Brillouin approximation, and study the electrolytes between two parallel dielectric surfaces. In comparison to other modified models, the new model is shown more accurate in agreement with particle-based simulations and capturing the physical properties of ionic structures near interfaces.
In numerical simulations of many charged systems at the micro/nano scale, a common theme is the repeated solution of the Poisson-Boltzmann equation. This task proves challenging, if not entirely infeasible, largely due to the nonlinearity of the equa
tion and the high dimensionality of the physical and parametric domains with the latter emulating the system configuration. In this paper, we for the first time adapt a mathematically rigorous and computationally efficient model order reduction paradigm, the so-called reduced basis method (RBM), to mitigate this challenge. We adopt a finite difference method as the mandatory underlying scheme to produce the {em truth approximations} of the RBM upon which the fast algorithm is built and its performance is measured against. Numerical tests presented in this paper demonstrate the high efficiency and accuracy of the fast algorithm, the reliability of its error estimation, as well as its capability in effectively capturing the boundary layer.
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.
The interaction force between likely charged particles/surfaces is usually repulsive due to the Coulomb interaction. However, the counterintuitive like-charge attraction in electrolytes has been frequently observed in experiments, which has been theo
retically debated for a long time. It is widely known that the mean field Poisson-Boltzmann theory cannot explain or predict this anomalous feature since it ignores many-body properties. In this paper, we develop efficient algorithm and perform the force calculation between two interfaces using a set of self-consistent equations which properly takes into account the electrostatic correlation and the dielectric-boundary effects. By solving the equations and calculating the pressure with the Debye-charging process, we show that the self-consistent equations could be used to study the attraction between like-charge surfaces from weak-coupling to mediate-coupling regime, and that the attraction is due to the electrostatics-driven entropic force which is significantly enhanced by the dielectric depletion of mobile ions. A systematic investigation shows that the interaction forces can be tuned by material permittivity, ionic size and valence, and salt concentration, and that the like-charge attraction exists only for specific regime of these parameters.
