In this note, we address some issues concerning the accurate pressure calculation of Coulomb syst... more In this note, we address some issues concerning the accurate pressure calculation of Coulomb systems with periodic boundary conditions. First, we prove that the formulas for the excess part of the pressure with Ewald summation also reduce to the ensemble average of one-third of the ratio between the potential energy and the volume so that the comments on our previous work in a recent paper by [Onegin et al., J. Phys. A: Math. Theor. 57 (2024) 205002] are incorrect. Second, we demonstrate that in charge non-neutral systems, the pressure expression must be corrected to include interactions with the neutralizing background. This addresses the issues about pressure computation in LAMMPS raised in the paper by Onegin et al.. Numerical experiments are performed to verify that the pressure obtained via Ewald summation with corrected terms agrees with the average pressure using thermodynamics for the non-neutral OCP system, and are independent of the splitting parameter in the Ewald summation.
Quasi-2D Coulomb systems are of fundamental importance and have attracted much attention in many ... more 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 varies numerical examples. Notably, it achieves a speedup of 2 ∼ 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.
The computational bottleneck of molecular dynamics is the pairwise additive long-range interactio... more The computational bottleneck of molecular dynamics is the pairwise additive long-range interactions between particles. The random batch Ewald (RBE) method provides a highly efficient and superscalable solver for long-range interactions, but the stochastic nature of this algorithm leads to unphysical self-heating effect during the simulation. We propose an energy stable scheme (ESS) for particle systems by employing a Berendsen-type energy bath. The scheme removes the notorious energy drift which exists due to the force error even when a symplectic integrator is employed. Combining the RBE with the ESS, the new method provides a perfect solution of the computational bottleneck of molecular dynamics at the microcanonical ensemble. Numerical results for primitive electrolyte and all-atom pure water systems demonstrate the attractive performance of the algorithm including its dramatically high accuracy, linear complexity and overcoming the energy drift for long-time simulations.
The motile microorganisms such as E. coli, sperm, or some seaweed are usually modelled by self-pr... more The motile microorganisms such as E. coli, sperm, or some seaweed are usually modelled by self-propelled particles that move with the run-and-tumble process. Individual-based stochastic models are usually employed to model the aggregation phenomenon at the boundary, which is an active research field that has attracted a lot of biologists and biophysicists. Self-propelled particles at the microscale have complex behaviors, while characteristics at the population level are more important for practical applications but rely on individual behaviors. Kinetic PDE models that describe the time evolution of the probability density distribution of the motile microorganisms are widely used. However, how to impose the appropriate boundary conditions that take into account the boundary aggregation phenomena is rarely studied. In this paper, we propose the boundary conditions for a 2D confined run-and-tumble model (CRTM) for self-propelled particle populations moving between two parallel plates with a run-and-tumble process. The proposed model satisfies the relative entropy inequality and thus long-time convergence. We establish the relation between CRTM and the confined Fokker-Planck model (CFPM) studied in [22]. We prove theoretically that when the tumble is highly forward peaked and frequent enough, CRTM converges asymptotically to the CFPM. A numerical comparison of the CRTM with aggregation and CFPM is given. The time evolution of both the deterministic PDE model and individual-based stochastic simulations are displayed, which match each other well.
Mathematical Models and Methods in Applied Sciences
The motile micro-organisms such as Escherichia coli, sperm, or some seaweed are usually modeled b... more The motile micro-organisms such as Escherichia coli, sperm, or some seaweed are usually modeled by self-propelled particles that move with the run-and-tumble process. Individual-based stochastic models are usually employed to model the aggregation phenomenon at the boundary, which is an active research field that has attracted a lot of biologists and biophysicists. Self-propelled particles at the microscale have complex behaviors, while characteristics at the population level are more important for practical applications but rely on individual behaviors. Kinetic PDE models that describe the time evolution of the probability density distribution of the motile micro-organisms are widely used. However, how to impose the appropriate boundary conditions that take into account the boundary aggregation phenomena is rarely studied. In this paper, we propose the boundary conditions for a 2D confined run-and-tumble model (CRTM) for self-propelled particle populations moving between two parall...
Coulomb interaction, following an inverse-square force-law, quantifies the amount of force betwee... more Coulomb interaction, following an inverse-square force-law, quantifies the amount of force between two stationary and 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 are developed to calculate 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 develop an accurate, highly efficient and scalable random batch Ewald (RBE) method to conduct ... more We develop an accurate, highly efficient and scalable random batch Ewald (RBE) method to conduct simulations in the isothermal-isobaric ensemble (the NPT ensemble) for charged particles in a periodic box. After discretizing the Langevin equations of motion derived using suitable Lagrangians, the RBE method builds the mini-batch strategy into the Fourier space in the Ewald summation for the pressure and forces so the computational cost is reduced from O(N 2) to O(N) per time step. We implement the method in the LAMMPS package and report accurate simulation results for both dynamical quantities and statistics for equilibrium for typical systems including all-atom bulk water and a semi-isotropic membrane system. Numerical simulations on massive supercomputing cluster are also performed to show promising CPU efficiency of RBE.
We have developed an accurate and efficient method for molecular dynamics simulations of charged ... more We have developed an accurate and efficient method for molecular dynamics simulations of charged particles confined by planar dielectric interfaces. The algorithm combines the image-charge method for near field with the harmonic surface mapping which converts the contribution of infinite far-field charges into a finite number of charges on an auxiliary spherical surface. We approximate the electrostatic potential of far-field charges via spherical harmonic expansion and determine the coefficients by fitting the Dirichlet-to-Neumann boundary condition, which only requires the potential within the simulation cell. Instead of performing the direct evaluation of spherical harmonic series expansion, we use Green's second identity to transform the series expansion into a spherical integral which can be accurately represented by discrete charges on the sphere. Therefore, the fast multipole method can be readily employed to sum over all charges within and on the sphere, achieving truly linear O(N) complexity. Our algorithm can be applied to a broad range of charged complex fluids under dielectric confinement.
We propose a harmonic surface mapping algorithm (HSMA) for electrostatic pairwise sums of an infi... more 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 HSMA 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.
We propose a fast method for the calculation of short-range interactions in molecular dynamics si... more 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 neighborlist method to avoid the construction of a full Verlet 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.
Coulomb interaction, following an inverse-square force-law, quantifies the amount of force betwee... more Coulomb interaction, following an inverse-square force-law, quantifies the amount of force between two stationary and 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 are developed to calculate 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 an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We implement two recently developed fast Coulomb solvers, HSMA3D [J. Chem. Phys. 149 (8) (2018) 0... more 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 develop an accurate, highly efficient, and scalable random batch Ewald (RBE) method to conduct... more We develop an accurate, highly efficient, and scalable random batch Ewald (RBE) method to conduct molecular dynamics simulations in the isothermal–isobaric ensemble (the NPT ensemble) for charged particles in a periodic box. After discretizing the Langevin equations of motion derived using suitable Lagrangians, the RBE method builds the mini-batch strategy into the Fourier space in the Ewald summation for the pressure and forces such that the computational cost is reduced to O(N) per time step. We implement the method in the Large-scale Atomic/Molecular Massively Parallel Simulator package and report accurate simulation results for both dynamical quantities and statistics for equilibrium for typical systems including all-atom bulk water and a semi-isotropic membrane system. Numerical simulations on massive supercomputing cluster are also performed to show promising central processing unit efficiency of the RBE.
We develop an accurate, highly efficient and scalable random batch sum-of-Gaussians (RBSOG) metho... more We develop an accurate, highly efficient and scalable random batch sum-of-Gaussians (RBSOG) method for molecular dynamics simulations of systems with long-range interactions. The idea of the RBSOG method is based on a sum-of-Gaussians decomposition of the Coulomb kernel, and then a random batch importance sampling on the Fourier space is employed for approximating the summation of the Fourier expansion of the Gaussians with large bandwidths (the long-range components). The importance sampling significantly reduces the computational cost, resulting in a scalable algorithm by avoiding the use of communication-intensive fast Fourier transform. Theoretical analysis is present to demonstrate the unbiasedness of the approximate force, the controllability of variance and the weak convergence of the algorithm. The resulting method has O(N) complexity with low communication latency. Accurate simulation results on both dynamical and equilibrium properties of benchmark problems are reported to illustrate the attractive performance of the method. Simulations on parallel computing are also performed to show the high parallel efficiency. The RBSOG method can be straightforwardly extended to more general interactions with long ranged kernels, and thus is promising to construct fast algorithms of a series of molecular dynamics methods for various interacting kernels.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We implement two recently developed fast Coulomb solvers, HSMA3D [J. Chem. Phys. 149 (8) (2018) 0... more 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.
In this note, we address some issues concerning the accurate pressure calculation of Coulomb syst... more In this note, we address some issues concerning the accurate pressure calculation of Coulomb systems with periodic boundary conditions. First, we prove that the formulas for the excess part of the pressure with Ewald summation also reduce to the ensemble average of one-third of the ratio between the potential energy and the volume so that the comments on our previous work in a recent paper by [Onegin et al., J. Phys. A: Math. Theor. 57 (2024) 205002] are incorrect. Second, we demonstrate that in charge non-neutral systems, the pressure expression must be corrected to include interactions with the neutralizing background. This addresses the issues about pressure computation in LAMMPS raised in the paper by Onegin et al.. Numerical experiments are performed to verify that the pressure obtained via Ewald summation with corrected terms agrees with the average pressure using thermodynamics for the non-neutral OCP system, and are independent of the splitting parameter in the Ewald summation.
Quasi-2D Coulomb systems are of fundamental importance and have attracted much attention in many ... more 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 varies numerical examples. Notably, it achieves a speedup of 2 ∼ 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.
The computational bottleneck of molecular dynamics is the pairwise additive long-range interactio... more The computational bottleneck of molecular dynamics is the pairwise additive long-range interactions between particles. The random batch Ewald (RBE) method provides a highly efficient and superscalable solver for long-range interactions, but the stochastic nature of this algorithm leads to unphysical self-heating effect during the simulation. We propose an energy stable scheme (ESS) for particle systems by employing a Berendsen-type energy bath. The scheme removes the notorious energy drift which exists due to the force error even when a symplectic integrator is employed. Combining the RBE with the ESS, the new method provides a perfect solution of the computational bottleneck of molecular dynamics at the microcanonical ensemble. Numerical results for primitive electrolyte and all-atom pure water systems demonstrate the attractive performance of the algorithm including its dramatically high accuracy, linear complexity and overcoming the energy drift for long-time simulations.
The motile microorganisms such as E. coli, sperm, or some seaweed are usually modelled by self-pr... more The motile microorganisms such as E. coli, sperm, or some seaweed are usually modelled by self-propelled particles that move with the run-and-tumble process. Individual-based stochastic models are usually employed to model the aggregation phenomenon at the boundary, which is an active research field that has attracted a lot of biologists and biophysicists. Self-propelled particles at the microscale have complex behaviors, while characteristics at the population level are more important for practical applications but rely on individual behaviors. Kinetic PDE models that describe the time evolution of the probability density distribution of the motile microorganisms are widely used. However, how to impose the appropriate boundary conditions that take into account the boundary aggregation phenomena is rarely studied. In this paper, we propose the boundary conditions for a 2D confined run-and-tumble model (CRTM) for self-propelled particle populations moving between two parallel plates with a run-and-tumble process. The proposed model satisfies the relative entropy inequality and thus long-time convergence. We establish the relation between CRTM and the confined Fokker-Planck model (CFPM) studied in [22]. We prove theoretically that when the tumble is highly forward peaked and frequent enough, CRTM converges asymptotically to the CFPM. A numerical comparison of the CRTM with aggregation and CFPM is given. The time evolution of both the deterministic PDE model and individual-based stochastic simulations are displayed, which match each other well.
Mathematical Models and Methods in Applied Sciences
The motile micro-organisms such as Escherichia coli, sperm, or some seaweed are usually modeled b... more The motile micro-organisms such as Escherichia coli, sperm, or some seaweed are usually modeled by self-propelled particles that move with the run-and-tumble process. Individual-based stochastic models are usually employed to model the aggregation phenomenon at the boundary, which is an active research field that has attracted a lot of biologists and biophysicists. Self-propelled particles at the microscale have complex behaviors, while characteristics at the population level are more important for practical applications but rely on individual behaviors. Kinetic PDE models that describe the time evolution of the probability density distribution of the motile micro-organisms are widely used. However, how to impose the appropriate boundary conditions that take into account the boundary aggregation phenomena is rarely studied. In this paper, we propose the boundary conditions for a 2D confined run-and-tumble model (CRTM) for self-propelled particle populations moving between two parall...
Coulomb interaction, following an inverse-square force-law, quantifies the amount of force betwee... more Coulomb interaction, following an inverse-square force-law, quantifies the amount of force between two stationary and 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 are developed to calculate 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 develop an accurate, highly efficient and scalable random batch Ewald (RBE) method to conduct ... more We develop an accurate, highly efficient and scalable random batch Ewald (RBE) method to conduct simulations in the isothermal-isobaric ensemble (the NPT ensemble) for charged particles in a periodic box. After discretizing the Langevin equations of motion derived using suitable Lagrangians, the RBE method builds the mini-batch strategy into the Fourier space in the Ewald summation for the pressure and forces so the computational cost is reduced from O(N 2) to O(N) per time step. We implement the method in the LAMMPS package and report accurate simulation results for both dynamical quantities and statistics for equilibrium for typical systems including all-atom bulk water and a semi-isotropic membrane system. Numerical simulations on massive supercomputing cluster are also performed to show promising CPU efficiency of RBE.
We have developed an accurate and efficient method for molecular dynamics simulations of charged ... more We have developed an accurate and efficient method for molecular dynamics simulations of charged particles confined by planar dielectric interfaces. The algorithm combines the image-charge method for near field with the harmonic surface mapping which converts the contribution of infinite far-field charges into a finite number of charges on an auxiliary spherical surface. We approximate the electrostatic potential of far-field charges via spherical harmonic expansion and determine the coefficients by fitting the Dirichlet-to-Neumann boundary condition, which only requires the potential within the simulation cell. Instead of performing the direct evaluation of spherical harmonic series expansion, we use Green's second identity to transform the series expansion into a spherical integral which can be accurately represented by discrete charges on the sphere. Therefore, the fast multipole method can be readily employed to sum over all charges within and on the sphere, achieving truly linear O(N) complexity. Our algorithm can be applied to a broad range of charged complex fluids under dielectric confinement.
We propose a harmonic surface mapping algorithm (HSMA) for electrostatic pairwise sums of an infi... more 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 HSMA 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.
We propose a fast method for the calculation of short-range interactions in molecular dynamics si... more 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 neighborlist method to avoid the construction of a full Verlet 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.
Coulomb interaction, following an inverse-square force-law, quantifies the amount of force betwee... more Coulomb interaction, following an inverse-square force-law, quantifies the amount of force between two stationary and 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 are developed to calculate 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 an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We implement two recently developed fast Coulomb solvers, HSMA3D [J. Chem. Phys. 149 (8) (2018) 0... more 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 develop an accurate, highly efficient, and scalable random batch Ewald (RBE) method to conduct... more We develop an accurate, highly efficient, and scalable random batch Ewald (RBE) method to conduct molecular dynamics simulations in the isothermal–isobaric ensemble (the NPT ensemble) for charged particles in a periodic box. After discretizing the Langevin equations of motion derived using suitable Lagrangians, the RBE method builds the mini-batch strategy into the Fourier space in the Ewald summation for the pressure and forces such that the computational cost is reduced to O(N) per time step. We implement the method in the Large-scale Atomic/Molecular Massively Parallel Simulator package and report accurate simulation results for both dynamical quantities and statistics for equilibrium for typical systems including all-atom bulk water and a semi-isotropic membrane system. Numerical simulations on massive supercomputing cluster are also performed to show promising central processing unit efficiency of the RBE.
We develop an accurate, highly efficient and scalable random batch sum-of-Gaussians (RBSOG) metho... more We develop an accurate, highly efficient and scalable random batch sum-of-Gaussians (RBSOG) method for molecular dynamics simulations of systems with long-range interactions. The idea of the RBSOG method is based on a sum-of-Gaussians decomposition of the Coulomb kernel, and then a random batch importance sampling on the Fourier space is employed for approximating the summation of the Fourier expansion of the Gaussians with large bandwidths (the long-range components). The importance sampling significantly reduces the computational cost, resulting in a scalable algorithm by avoiding the use of communication-intensive fast Fourier transform. Theoretical analysis is present to demonstrate the unbiasedness of the approximate force, the controllability of variance and the weak convergence of the algorithm. The resulting method has O(N) complexity with low communication latency. Accurate simulation results on both dynamical and equilibrium properties of benchmark problems are reported to illustrate the attractive performance of the method. Simulations on parallel computing are also performed to show the high parallel efficiency. The RBSOG method can be straightforwardly extended to more general interactions with long ranged kernels, and thus is promising to construct fast algorithms of a series of molecular dynamics methods for various interacting kernels.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We propose an accurate algorithm for a novel sum-of-exponentials (SOE) approximation of kernel fu... more 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 approximating a continuous temporal convolution integral. The SOE method is constructed by a combination of the de la Vallée-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 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 both the new SOE method and the SOEbased 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.
We implement two recently developed fast Coulomb solvers, HSMA3D [J. Chem. Phys. 149 (8) (2018) 0... more 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.
Uploads
Papers by Jiuyang Liang